Advanced Algorithms & Data Structures: Lecturer: Karimzhan Nurlan Berlibekuly

ADVANCED
ALGORITHMS &
DATA STRUCTURES
Lecture-05
Fast Fourier Transform (FFT)
Lecturer: Karimzhan Nurlan Berlibekuly

[email protected]
• The Fourier transform is one of the most fundamental concepts
in the information sciences. It’s a method for representing an
irregular signal — such as the voltage fluctuations in the wire
that connects an MP3 player to a loudspeaker — as a
combination of pure frequencies. It’s universal in signal
processing, but it can also be used to compress image and audio
files, solve differential equations and price stock options, among
other things.
• The reason the Fourier transform is so prevalent is an algorithm
called the fast Fourier transform (FFT), devised in the mid-1960s,
which made it practical to calculate Fourier transforms on the fly.
Ever since the FFT was proposed, however, people have
wondered whether an even faster algorithm could be found.
Fast Fourier Transform
•A polynomial A(x) can be written in the following forms:
A(x) = a0 + a1x + a2x2 + · · · + an-1xn-1

=
= (a0, a1, a2, . . . , an-1 ) (coefficient vector)
The degree of A is n – 1.
Fast Fourier Transformation for polynomial
multiplication
A(x) * B(x) = C(x) ?
O(n2)
Example:
A(x) = 6x3 + 7x2 - 10x + 9
B(x) = -2x3 + 4x – 5
Coefficient representation of A(x) = (9, -10, 7, 6)
Coefficient representation of B(x) = (-5, 4, 0, -2)
Input :
A[] = {9, -10, 7, 6}
B[] = {-5, 4, 0, -2}
Output :
We can do better, if we represent the polynomial in
another form.
Idea is to represent polynomial in point-value form

and then compute the product. A point-value
representation of a polynomial A(x) of degree-
bound n is a set of n point-value pairs is{ (x0, y0),
(x1, y1), …, (xn-1, yn-1)} such that all of the xi are
distinct and yi = A(xi) for i = 0, 1, …, n-1.
Example:
A(x) = x3 - 2x + 1
xi : 0, 1, 2, 3
A(xi) : 1, 0, 5, 22
Point-value representation of above polynomial is { (0, 1), (1, 0), (2, 5), (3, 22) }.
Using Horner’s method, n-point evaluation takes time O(n2). It’s just calculation
of values of A(x) at some x for n different points, so time complexity is O(n 2).
Now that the polynomial is converted into point value, it can be easily calculated
C(x) = A(x)*B(x) again using horner’s method. This takes O(n) time. An
important point here is C(x) has degree bound 2n, then n points will give only n
points of C(x), so for that case we need 2n different values of x to calculate 2n
different values of y. Now that the product is calculated, the answer can to be
converted back into coefficient vector form. To get back to coefficient vector
form we use inverse of this evaluation. The inverse of evaluation is called
interpolation. Interpolation using Lagrange’s formula gives point value-form to
coefficient vector form of the polynomial.
Lagrange’s formula is –
O(n2)
Discrete Fourier transform
(DFT)
This idea still solves the problem in O(n2) time complexity. We can
use any points we want as evaluation points, but by choosing the
evaluation points carefully, we can convert between
representations in only O(n log n) time. If we choose “complex
roots of unity” as the evaluation points, we can produce a point-
value representation by taking the discrete Fourier transform
(DFT) of a coefficient vector. We can perform the inverse
operation, interpolation, by taking the “inverse DFT” of point-value
pairs, yielding a coefficient vector. Fast Fourier Transform (FFT) can
perform DFT and inverse DFT in time O(nlogn).
Discrete Fourier transform
(DFT)
DFT is evaluating values of polynomial at n complex nth roots of
unity
So, for
k = 0, 1, 2, …, n-1, y = (y0, y1, y2, …, yn-1) is Discrete fourier
Transformation (DFT) of given polynomial.
The product of two polynomials of degree-bound n is a polynomial
of degree-bound 2n. Before evaluating the input polynomials A
and B, therefore, we first double their degree-bounds to 2n by
adding n high-order coefficients of 0. Because the vectors have 2n
elements, we use “complex 2nth roots of unity, ” which are
denoted by the W2n (omega 2n). We assume that n is a power of
2; we can always meet this requirement by adding high-order zero
coefficients.
Here is the Divide-and-conquer strategy to solve this problem.
Define two new polynomials of degree-bound n/2, using even-index and odd-index
coefficients of A(x) separately
The problem of evaluating A(x) at
reduces to evaluating the degree-bound n/2 polynomials A0(x) and A1(x) at the
points
Now combining the results by A(x) = A0(x^2) + xA1(x^2)

The list
does not contain n distinct values, but n/2 complex n/2th roots of unity.
Polynomials A0 and A1 are recursively evaluated at the n complex nth roots of
unity. Subproblems have exactly the same form as the original problem, but are
half the size. So recurrence formed is T(n) = 2T(n/2) + O(n), i.e complexity
O(nlogn).
Algorithm
1. Add n higher-order zero coefficients to A(x) and B(x)
2. Evaluate A(x) and B(x) using FFT for 2n points
3. Pointwise multiplication of point-value forms
4. Interpolate C(x) using FFT to compute inverse DFT
Pseudo code of recursive FFT
Recursive_FFT(a){ y0 = Recursive_FFT(A0) // local array
n = length(a) // a is the input coefficient y1 = Recursive-FFT(A1) // local array
vector
if n = 1 for k = 0 to n/2 - 1
then return a
// y array stores values of the DFT
// wn is principle complex nth root of unity. // of given polynomial.
wn = e^(2*pi*i/n) do y[k] = y0[k] + w*y1[k]
w=1 y[k+(n/2)] = y0[k] - w*y1[k]
w = w*wn
// even indexed coefficients return y
A0 = (a0, a2, ..., an-2 ) }
// odd indexed coefficients

A1 = (a1, a3, ..., an-1 )
Recursion Tree of Above Execution
Why does this work ?
since,
Thus, the vector y returned by Recursive-FFT is indeed the DFT of the input vector a.
FFT - Source Code (1)
#include <bits/stdc++.h> vector<cd> A0(n / 2), A1(n / 2);
using namespace std; for (int i = 0; i < n / 2; i++) {
// even indexed coefficients

// For storing complex values of nth roots
A0[i] = a[i * 2];
// of unity we use complex<double>
typedef complex<double> cd; // odd indexed coefficients
A1[i] = a[i * 2 + 1];
// Recursive function of FFT }
vector<cd> fft(vector<cd>& a)
{ // Recursive call for even indexed coefficients
int n = a.size(); vector<cd> y0 = fft(A0);
// Recursive call for odd indexed coefficients

// if input contains just one element vector<cd> y1 = fft(A1);
if (n == 1)
return vector<cd>(1, a[0]); // for storing values of y0, y1, y2, ..., yn-1.
vector<cd> y(n);
// For storing n complex nth roots of unity
vector<cd> w(n); for (int k = 0; k < n / 2; k++) {
for (int i = 0; i < n; i++) { y[k] = y0[k] + w[k] * y1[k];
y[k + n / 2] = y0[k] - w[k] * y1[k];
double alpha = 2 * M_PI * i / n;
}
w[i] = cd(cos(alpha), sin(alpha)); return y;
} }
FFT - Source Code (2)
// Driver code Input: 1 2 3 4
int main()
{ Output:
vector<cd> a{1, 2, 3, 4}; (10, 0)
vector<cd> b = fft(a);
for (int i = 0; i < 4; i++) (-2, -2)
cout << b[i] << endl; (-2, 0)
return 0; (-2, 2)
}
Interpolation
Switch the roles of a and y.
Replace wn by wn^-1.
Divide each element of the result by n.
Time Complexity : O(nlogn).
Interpolation
Switch the roles of a and y.
Replace wn by wn^-1.
Divide each element of the result by n.
Time Complexity : O(nlogn).
Applications
Fourier (frequency) space many applications. The polynomial A∗ = FFT(A) is
complex, and the amplitude |a∗ k| represents the amplitude of the frequency-k signal,
while arg(ak ∗) (the angle of the 2D vector) represents the phase shift of that signal.
For example, this perspective is particularly useful for audio processing, as used by
Adobe Audition, Audacity, etc.:
• High-pass filters zero out high frequencies

• Low-pass filters zero out low frequencies
• Pitch shifts shift the frequency vector
• Used in MP3 compression, etc.
Cooley–Tukey FFT algorithm
• The Cooley–Tukey algorithm, named after J. W. Cooley and John Tukey, is the most
common fast Fourier transform (FFT) algorithm. It re-expresses the
discrete Fourier transform (DFT) of an arbitrary composite size N = N1N2 in terms
of N1 smaller DFTs of sizes N2, recursively, to reduce the computation time to
O(N log N) for highly composite N (smooth numbers). Because of the algorithm's
importance, specific variants and implementation styles have become known by their
own names, as described below.
• Because the Cooley–Tukey algorithm breaks the DFT into smaller DFTs, it can be
combined arbitrarily with any other algorithm for the DFT. For example, Rader's or
Bluestein's algorithm can be used to handle large prime factors that cannot be
decomposed by Cooley–Tukey, or the prime-factor algorithm can be exploited for
greater efficiency in separating out relatively prime factors.
• The algorithm, along with its recursive application, was invented by
Carl Friedrich Gauss. Cooley and Tukey independently rediscovered and popularized it
160 years later.
Realisations of code:
https://rosettacode.org/wiki/Fast_Fourier_transform
Cooley–Tukey FFT algorithm
#include <complex> // inverse fft (in-place)

#include <iostream> void ifft(CArray& x)
#include <valarray> {
const double PI = 3.141592653589793238460; // conjugate the complex numbers
typedef std::complex<double> Complex; x = x.apply(std::conj);
typedef std::valarray<Complex> CArray; // forward fft
// Cooley–Tukey FFT (in-place, divide-and-conquer) fft( x );
// Higher memory requirements and redundancy although // conjugate the complex numbers again
more intuitive x = x.apply(std::conj);
void fft(CArray& x) // scale the numbers
{ x /= x.size();
const size_t N = x.size(); }
if (N <= 1) return; int main()
// divide {
CArray even = x[std::slice(0, N/2, 2)]; const Complex test[] = { 1.0, 1.0, 1.0, 1.0, 0.0,
CArray odd = x[std::slice(1, N/2, 2)]; 0.0, 0.0, 0.0 };
// conquer CArray data(test, 8);
fft(even); // forward fft
fft(odd); fft(data);
// combine std::cout << "fft" << std::endl;
for (size_t k = 0; k < N/2; ++k) for (int i = 0; i < 8; ++i)
{ {
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k]; std::cout << data[i] << std::endl;
x[k ] = even[k] + t; }
x[k+N/2] = even[k] - t; // inverse fft
} ifft(data);
} std::cout << std::endl << "ifft" << std::endl;
for (int i = 0; i < 8; ++i)
{
std::cout << data[i] << std::endl;
}
return 0;
}
Additional Resources
https://e-maxx.ru/algo/fft_multiply
https://habr.com/ru/company/otus/blog/449996/

Advanced Algorithms & Data Structures: Lecturer: Karimzhan Nurlan Berlibekuly

Uploaded by

Copyright:

Available Formats

Advanced Algorithms & Data Structures: Lecturer: Karimzhan Nurlan Berlibekuly

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Advanced Algorithms & Data Structures: Lecturer: Karimzhan Nurlan Berlibekuly

Uploaded by

Copyright:

Available Formats

ADVANCED

Lecturer: Karimzhan Nurlan Berlibekuly

A(x) = a0 + a1x + a2x2 + · · · + an-1xn-1

= (a0, a1, a2, . . . , an-1 ) (coefficient vector)

A(x) * B(x) = C(x) ?

A(x) = 6x3 + 7x2 - 10x + 9

Idea is to represent polynomial in point-value form

Now combining the results by A(x) = A0(x^2) + xA1(x^2)

// odd indexed coefficients

// even indexed coefficients

// Recursive call for odd indexed coefficients

• High-pass filters zero out high frequencies

#include <complex> // inverse fft (in-place)

You might also like