A Small Cooley-Tukey FFT

Compute an 8-point discrete Fourier transform by recursively splitting even and odd samples - the divide-and-conquer that made spectral methods practical.

The DFT and its cost

The discrete Fourier transform X[k] = sum over n of x[n] exp(-2 pi i n k / N) turns samples into frequencies. Computed directly it costs order N^2 operations. The Cooley-Tukey FFT cuts this to order N log N by exploiting symmetry.

Divide and conquer

Split the sequence into even-indexed and odd-indexed samples. The DFT of the whole is built from the DFTs of the two halves combined with twiddle factors: X[k] = E[k] + W^k O[k] and X[k+N/2] = E[k] - W^k O[k], where W = exp(-2 pi i / N). Recurse until length 1.

import numpy as np
def fft(x):
    N=len(x)
    if N==1: return x
    E=fft(x[0::2]); O=fft(x[1::2])
    W=np.exp(-2j*np.pi*np.arange(N//2)/N)
    return np.concatenate([E+W*O, E-W*O])
x=np.array([1,2,3,4,4,3,2,1],float)
print(np.round(fft(x),3))
print(np.allclose(fft(x),np.fft.fft(x)))  # matches numpy

Why N log N matters

For N = 8 the saving is modest, but for a million-point transform, N^2 is a trillion operations while N log N is about twenty million - the difference between infeasible and instant. This speedup made digital signal processing, spectral PDE solvers, and the Poisson solves inside plasma codes routine.

Requirements and variants

This radix-2 form needs N to be a power of two; mixed-radix and Bluestein algorithms handle arbitrary lengths. For real input, symmetry halves the work again. The inverse transform is the same algorithm with a conjugated twiddle and a 1/N scaling.