Problem 29 : An Efficient way to compute Discrete Fourier Transform

Discrete Fourier Transform (DFT) transforms a sequence of values (could be values generated from a time-dependent function) into another sequence of values that represent the values of the Fourier Transform of the time-dependent function.

Mathematically, the DFT of a vector of N complex numbers X is a vector of N complex numbers Y given by,

\[Y[j]=\sum_{k=0}^{N-1} X[k] W_N^{j k}\]

where \(W_N = e^{-2 \pi i / N}\) denotes the \(N^{th}\) complex root of unity. On the complex plane, these points lie on the unit circle in the complex plane.

Since the sum goes over \(N\) terms for each element of \(Y\) and \(Y\) has \(N\) terms, a naive sum would require O(\(N^2\)) computations. But, can we do better? The answer is a resounding yes. When I was first posed with this problem, I was pretty certain that it’d either employ divide-and-conquer or dynamic programming strategies. The idea here is to break the problem in to sub-problems that can then be recursively solved and combined to solve the full problem!

Cooley-Tukey Algorithm

First note that the formula for DFT can be slightly generalized and written as a polynomial function \(P(x)\) as,

\[P(x)=\sum_{k=0}^{N-1} X[k] x^k\]

where we are interested in evaluating the polynomial at special points which correspond to \(x = W_N^j\) for evaluating \(Y[j]\).

The simplest division that we can do should involve dividing the polynomial into two equal (roughly) parts. Let us try and see if dividing it into terms with even powers and odd powers helps in reducing the complexity.

It is easy to see that one can write the polynomial \(P(x)\) as a sum of even and odd powered terms as follows,

\[P(x) = A(x^2) + x*B(x^2)\]

where \(A(x^2) = \sum_{k=0}^{N/2-1} X[2k] x^k\) and \(B(x^2) = \sum_{k=0}^{N/2-1} X[2k+1] x^k\).

Well, so far things don’t look super interesting. But, they become interesting once we observe the points at which we are interested in evaluating such a polynomial correspond to the \(N^{th}\) complex roots of unity since \(Y[j] = P(W_N^j) = A(W_N^{2j}) + W_N^j*B(W_N^{2j})\).

Now, \(A(W_N^{2j}) = \sum_{k=0}^{N/2-1} X[2k] W_N^{2j.k}\) looks like another FFT but this time it has \(N/2\) terms. So, the \(N/2^{th}\) root of unity applies here i.e. equal to \(e^{-2 \pi 2 i / N}\) which is nothing but \(W_N^2\) ! So, \(A(x^2)\) evaluated at the \(W_N^{2j}\) which, in turn, is the FFT of the sequence \(\left(X_0, X_2, \ldots, X_{N-2}\right)\). We are staring at the core of Cooley-Tukey Algorithm which essentially divides the FFT problems into two sub-problems of half the size.

A small trick

A small but important track to further reduce the complexity is to note that, the following holds true,

\[Y[j+N/2] = P(W_N^{j+N/2}) = A(W_N^{2j+N}) + W_N^{j+N/2}*B(W_N^{2j+N}) = A(W_N^{2j}) - W_N^{j}*B(W_N^{2j})\]

where the last equality can be shown using the Euler’s formula, \(W_N = e^{-2 \pi i / N}\).