JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

Cooley-Tukey: From O(N^2) to O(N log N)

Guide 1 built the discrete Fourier transform honestly — and showed it costs O(N^2). Now we watch one of the great algorithms of the century halve the problem, then halve it again, until that quadratic wall collapses into O(N log N). The trick is older than the computer, hiding in a single algebraic split.

Why O(N^2) is a wall, not a nuisance

The previous guide built the discrete Fourier transform from scratch and laid its cost bare: to get all N output frequencies you compute N sums, and each sum touches all N inputs, so you do about N^2 complex multiply-adds. For a toy signal of N = 1000 that is a million operations — instant. But Fourier transforms live at scale: a one-megapixel image edge has N near 10^6, and N^2 is then 10^12, a trillion operations for a single transform. At a billion operations a second that is a quarter of an hour, for one image. Audio, radar, MRI, and weather models call the transform millions of times. At O(N^2) the whole enterprise is simply impossible.

So the question is not 'can we shave a constant factor?' — that would not help. We need to change the growth rate itself, the exponent on N. This is the heart of computational complexity: an O(N log N) algorithm is not a bit faster than O(N^2), it is a different universe. At N = 10^6, N^2 is 10^12 but N log N (base 2) is about 2 times 10^7 — a speed-up of roughly fifty thousand. That trillion-operation transform becomes twenty million operations, milliseconds instead of minutes. Finding such an algorithm is what makes the FFT one of the most consequential pieces of mathematics ever written into code.

The one algebraic split that does everything

Here is the whole secret, and it fits in one idea. Assume N is even. Split your N inputs into the ones at EVEN positions (x_0, x_2, x_4, ...) and the ones at ODD positions (x_1, x_3, x_5, ...). Each half has N/2 points. The astonishing fact is that the size-N DFT of the full signal can be rebuilt, for free, from the size-(N/2) DFT of the evens and the size-(N/2) DFT of the odds. The DFT sum, when you separate even and odd index terms and factor out a common phase, reassembles exactly into two smaller transforms glued together. We have turned one problem of size N into two problems of size N/2 — the signature move of divide and conquer.

Let w = exp(-2*pi*i / N)              # the basic twiddle (Nth root of unity)
For each output index k = 0 .. N/2 - 1 :
    E_k = DFT of the even-indexed samples, at frequency k   # size N/2
    O_k = DFT of the odd-indexed samples,  at frequency k   # size N/2
    X_k         = E_k +  w^k * O_k        # one butterfly...
    X_{k + N/2} = E_k -  w^k * O_k        # ...gives TWO outputs
The Cooley-Tukey recombination: each pair of outputs k and k+N/2 comes from the SAME E_k and O_k, differing only in the sign of the twiddle factor w^k. One multiply, two outputs.

Look closely at why two outputs come from one multiplication. The output X_k uses E_k + w^k O_k, and the output half a turn away, X_{k+N/2}, uses E_k - w^k O_k — the SAME two numbers, the SAME twiddle, just a flipped sign. This little plus-and-minus pair is the butterfly, named for the crossing-wings shape when you draw the data flow. Because each butterfly produces two outputs from one complex multiply, the recombination of the two halves costs only about N/2 multiplies, not N^2. The expensive part has been outsourced to the two smaller transforms.

Recurse all the way down: where N log N is born

We have not finished — we have only halved. But the same split applies to each N/2 transform: cut it into its own evens and odds to make two size-N/4 transforms, and so on. If N is a power of two you can keep halving until you reach transforms of size 1, and the DFT of a single number is just that number itself — nothing to compute. So the recursion bottoms out cleanly. Count the work by levels: each level of halving costs about N multiply-adds to recombine, and there are log_2(N) levels because that is how many times you can halve N down to 1. Total work: about N times log_2(N), the celebrated O(N log N).

  1. Base case: if N = 1, the transform of one sample IS that sample. Return it.
  2. Split the N samples into the even-indexed half and the odd-indexed half, each of size N/2.
  3. Recursively FFT each half — this is the algorithm calling itself on a smaller problem.
  4. Combine: for k = 0 .. N/2 - 1, form the butterfly X_k = E_k + w^k O_k and X_{k+N/2} = E_k - w^k O_k.
  5. Return the N combined outputs. The recursion did the heavy lifting; this level only glues.

This is the Cooley-Tukey algorithm (Cooley-Tukey FFT), published in 1965 — though Gauss had quietly used the same halving idea on hand calculations around 1805, before Fourier's own work was even published. The version above is the radix-2 form, the cleanest to understand. It is exactly the same DFT as guide 1 — same outputs, same mathematics — just evaluated in a smarter order. The fast Fourier transform is not an approximation of the DFT; it is the identical answer, reorganized so that shared sub-computations are done once instead of N times.

What actually makes it fast on a real machine

The flop count is the headline, but it is not the whole story, and an honest account matters. The reordering into evens and odds, repeated at every level, ends up permuting the inputs into bit-reversed order — index 6 (binary 110) swaps with index 3 (binary 011), and so on. Production FFTs do this permutation in place, then sweep through log_2(N) passes of butterflies without recursion at all. The recursion is the clean way to understand it; the loop is the fast way to run it. The twiddle factors w^k are usually precomputed into a small table, because recomputing a sine and cosine inside the inner loop would dominate the cost.

And the inverse? The inverse DFT is the same shape of sum with the sign of the exponent flipped and a 1/N scaling, so the very same butterflies, with conjugated twiddles, run it at the same O(N log N) cost. One algorithm, run forwards and backwards, gives you both directions of the transform — which is exactly what the next two guides will exploit to make convolution fast and to differentiate functions with astonishing accuracy.

Honest caveats: it is fast, not magic, and not exact

Two honest qualifications keep you from over-trusting the FFT. First, the clean radix-2 story needs N to be a power of two. Real signals are not so polite — you might have N = 1000, or a prime. Libraries handle this with mixed-radix splits (any factorization N = p times q works the same way) and, for awkward N, Bluestein's trick that re-expresses the transform as a convolution. The upshot: a good FFT runs in O(N log N) for essentially every N, but a prime length is slower than the nearest power of two, which is why people often pad to a convenient size when they can.

Second — and this is the lesson the whole ladder keeps teaching — the FFT computes in floating-point arithmetic, so its output is APPROXIMATE, not the exact DFT. Here, though, is the genuinely good news: the FFT is more accurate than the naive O(N^2) sum, not less. The direct sum accumulates round-off across N additions, so its error grows like O(N) times the unit roundoff. The FFT only ever adds in a balanced tree of depth log_2(N), so its round-off grows like O(log N) — the same divide-and-conquer that saves time also limits error accumulation. Fewer, better-organized operations means both faster AND more accurate. That is a rare and beautiful win.