Understanding the FFT Algorithm (2013)

This post summarizes the core idea behind FFT, which computes the Discrete Fourier Transform (DFT) quickly. A naive computation of the DFT is O(N^2), but FFT reduces this to O(N log N).

First, the DFT is expressed as a matrix-vector product, and a slow Python implementation DFT_slow is built and compared against numpy.fft.fft. This implementation is correct but is much slower for 1024 inputs, differing by roughly more than 1000x.

The key is the periodic symmetry of the DFT. The author shows the property X_{N+k} = X_k, and based on this, recursively decomposes the input into two smaller DFTs split by even indices and odd indices. This process is the essence of the Cooley-Tukey FFT, and by continuing to split the problem in half until it becomes small enough, the overall complexity becomes O(N log N).

The implementation is presented in two stages.

• FFT: a pure Python + NumPy version that performs recursive divide-and-conquer
• FFT_vectorized: a vectorized version that reduces recursive calls and computes multiple subproblems at once

The vectorized version is one step faster again than the recursive version, and shows fairly close performance to numpy.fft.fft for an input of size 1024×16. However, it still falls short of FFTPACK-based implementations, which is explained as being due to the costs of Python recursion, temporary array creation, and memory copying.

Finally, the reasons FFTPACK is faster are summarized.

• It maximizes the reuse of intermediate computations
• It controls memory usage more precisely using a low-level language like Fortran
• It utilizes splitting schemes other than just radix-2
• Alternative FFT algorithm families such as Bluestein and Rader also exist

Ultimately, the goal of this post is to go beyond using FFT as a black box in practice, and instead give an intuition for why it is fast, what symmetry it exploits, and how it can be implemented directly in Python.