The Key to the FFT

Christian Zuniga, PhD

The Fast Fourier Transform (FFT) is the most important algorithm of the 20th century [1]. It was popularized by Cooley and Tukey [2] but may have originally been discovered by Gauss [3]. The FFT is critical in areas such as communications, compression, signal processing, medical imaging, and numerical simulation. Without the FFT, many applications would not be practical. Figure 1 shows the difference in computation time between computing the DFT with the FFT and a direct computation.

Figure 1 Difference in computation time between straightforward computation of DFT (blue) and the FFT (red) (Image by author)

As the name implies, the FFT implements a fast version of the Discrete Fourier Transform (DFT).

The most straightforward algorithm computes the DFT using 2 for-loops, an inner loop over n and an outer over k as shown in the function DFT(x,N) in Figure 2. Computation time scales quadratically with N, the length of the input x, or the complexity is O(N^2).

Figure 2 DFT computed with 2-for loops (Image by author)

The Cooley-Tukey FFT speeds up computation using divide and conquer [2][4]. The DFT is broken up into DFTs of smaller size, but of the same form as the original DFT. The idea is to relate DFTs of size N/2 to a DFT of size N. To simplify the problem, N is assumed to be a power of 2, N = 2^r. As a first step, we would like to break up the summation into 2 sums, each of length N/2. The DFT sum would reach N/2–1 and the complex exponential has an N/2 instead of N. The key insight in obtaining the FFT is the halving lemma.

By a simple rule of division, 1/(N/2) = 2/N, the factor of 2 multiplies the index n instead of dividing N. Writing out the summation explicitly allows to see how this identity can be useful. We are looking for a factor of the form 2n in the exponentials.

After examining the summation we will find the factor 2n appears in the even elements of the summation or n = 0,2, 4,..

The subsequence of even elements calculates a DFT of length N/2. Placing all the even elements first and using the halving lemma show this more clearly.

What about the odd elements? A flash of insight (or a lot of thought) will reveal that factoring out e^(-i2πk/N) will also give a DFT of length N/2.

The original N point DFT has been broken up into computing 2 DFTs of length N/2 and then combining their result or in terms of the function, DFT(x, N) = DFT(xe,N/2)+ e^(-i2πk/N) DFT(xo,N/2).

There is a slight problem with this approach, however. In the smaller DFTs, the index k would only reach N/2–1, not N-1. How to compute the DFT for the remaining k=N/2…N-1? The properties of the complex exponentials resolve this issue. Calling the DFT of the even sequence Xe, and that of the odd Xo.

Assuming the DFTs Xe and Xo are available, the DFT of x can be computed with an amount of work linearly proportional to N as shown in Figure 3 where M=N/2. Technically it’s N/2 but this can be ignored for the general speedup analysis.

Figure 3 Computing N point DFT from two N/2 point DFTs. (Image by author)

The same principle can be used to obtain Xe and Xo. The respective sequences are again split into 2 additional even and odd sequences and do 4, N/4 length DFTs. The process is recursive and reaches a base case of length 1 after r function calls as shown in Figure 4.

Figure 4 FFT implementation with recursion (Image by author)

Figure 5 shows an example recursive tree of function calls for N=8. The tree helps to explain the speedup obtained by the FFT. The total time depends on (the number of levels) x (work per level). There are r+1 levels or logN + 1. Not counting the recursive calls, the amount of work at each level j depends on the number of DFTs to compute, which is 2^j, times the work per DFT which is N/2^j , the sequence length at that level. This gives a constant amount of work N per level. The end result is that total time scales as N logN giving the speedup observed in Figure 1. In practice an iterative version of the FFT is slightly faster and also leads to a parallel version of the algorithm.

Figure 5 FFT recursion tree for N=8, r= 3. The number of elements at each level j is 2j where the j=0 is the topmost level (Image by author).

References

[1] Strang “Computational Science and Engineering” Wellesley-Cambridge Press 2007

[2] Cooley, James W.; Tukey, John W. (1965). “An algorithm for the machine calculation of complex Fourier series”. Math. Comput. 19 (90): 297–301. doi:10.2307/2003354. JSTOR 2003354.

[3] https://en.wikipedia.org/wiki/Fast_Fourier_transform#Cooley%E2%80%93Tukey_algorithm

[4] Cormen, et. al “Introduction to Algorithms” Third Edition, MIT Press 2009