8.5 The fast Fourier transform

The DFT is remarkably simple to implement. The definition is a sum, and a fully vectorized version is essentially a single matrix multiplication in NumPy:

def dft(x: np.ndarray) -> np.ndarray:
    N = len(x)
    k, n = np.arange(N).reshape(-1, 1), np.arange(N).reshape(1, -1)
    phasors = np.exp(-2j * np.pi * k * n / N)  # (N, N): row k, column n
    return (phasors * x).sum(axis=1)           # weight each row by x, sum over n

Writing the same computation as an explicit double loop makes its cost visible. For each of the \(N\) output bins, we sum over all \(N\) input samples:

def dft_unrolled(x: np.ndarray) -> np.ndarray:
    N = len(x)
    out = np.zeros(N, dtype=np.complex128)
    for k in range(N):
        for n in range(N):
            out[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return out

Those nested loops reveal that the DFT is an \(O(N^2)\) computation. For a short window this is fine, but audio windows are often thousands of samples long, and we may compute the DFT thousands of times per second of audio. Quadratic cost quickly becomes prohibitive. Can we do better?

The same idea was, it was later discovered, known to Gauss around 1805, a century and a half before Cooley and Tukey.

It turns out we can do asymptotically better. The key insight, popularized by James Cooley and John Tukey in 1965, is that an \(N\)-point DFT can be expressed in terms of two \(N/2\)-point DFTs: one over the even-indexed samples and one over the odd-indexed samples. Computer science students will recognize this as divide and conquer, the same recursive strategy behind algorithms like merge sort. We split the problem in half, solve each half recursively, and combine the results.

Fig. 44 A radix-2 FFT drawn as a butterfly diagram for \(N = 8\). The even- and odd-indexed samples are transformed by two \(N/2\)-point DFTs, producing \(E[k]\) and \(O[k]\). Each pair then feeds two outputs, \(X[k] = E[k] + W_N^k\, O[k]\) and \(X[k+4] = E[k] - W_N^k\, O[k]\), where \(W_N^k = e^{-2\pi j k / N}\) is a “twiddle factor”. The crossing lines that combine each pair give the diagram its butterfly shape. Applying this split recursively is the FFT.

The combining step, known as the butterfly, merges the two half-size results in \(O(N)\) time. Recursing all the way down gives \(\log_2 N\) levels, each costing \(O(N)\), for a total of \(O(N \log N)\). That is an enormous improvement over \(O(N^2)\): for a 4096-sample window, it is the difference between roughly 16 million operations and about 50 thousand. In code, the recursion handles the \(\log N\) levels while a vectorized butterfly handles the \(O(N)\) combining at each level:

def fft(x: np.ndarray) -> np.ndarray:
    x = np.asarray(x, dtype=np.complex128)
    N = len(x)
    if N == 1:
        return x
    even, odd = fft(x[::2]), fft(x[1::2])          # divide and recurse
    twiddle = np.exp(-2j * np.pi * np.arange(N // 2) / N) * odd
    return np.concatenate([even + twiddle, even - twiddle])  # butterfly

The full runnable code, including a check that all three implementations agree with NumPy’s optimized FFT, is in code/dft.py.

The FFT is probably the most consequential algorithm in all of digital signal processing, underpinning not just audio analysis but multimedia compression, wireless communication, and much more. Understanding the high-level behavior of the algorithm (divide-and-conquer) and its asymptotic \(O(N \log N)\) performance is far more important than actually implementing the algorithm or understanding the “butterfly” details. In practice you will call a highly-tuned library routine such as np.fft.fft (or np.fft.rfft for real signals), which combines these high-level ideas with additional low-level optimizations.