Every FFT tutorial shows you the same trick: split a signal of length n into its even and odd samples, transform each half, combine. It is elegant, it is fast, and it only works when n is a power of two. Real signals are rarely so obliging — you record 1000 samples, or 44,100, or whatever a sensor happened to hand you — and a library that only transforms powers of two isn't a library, it's a demo.

So Equana's fft has two engines behind one function. Power-of-two lengths run the classic iterative radix-2 Cooley–Tukey algorithm. Every other length runs Bluestein's chirp-z transform, which is a small piece of algorithmic sleight-of-hand: it turns an arbitrary-length transform into a power-of-two convolution, and then uses the fast engine to do that convolution. This post is about how those two fit together, with just enough math to see why it works.

The problem: the DFT is O(n²)

The discrete Fourier transform is a single formula. Given n complex samples x₀ … xₙ₋₁, it produces n coefficients:

X_k = Σ_j x_j · exp(-2πi · j·k / n) for k = 0 … n-1

Read it literally and it is a matrix–vector product: each of the n outputs is a sum over all n inputs, so computing the whole transform is n² complex multiply-adds. That is fine at n = 64 and a catastrophe at n = 1,000,000 — a quadratic curve that makes the "obvious" DFT useless for anything real. Every fast algorithm below computes exactly the same numbers; they just refuse to do n² work to get them.

Cooley–Tukey: divide, conquer, and twiddle

The insight is that the transform of a power-of-two-length signal contains two smaller transforms hiding inside it. Split the sum into even-indexed and odd-indexed samples:

X_k = Σ_j x_2j · exp(-2πi·(2j)k/n) + Σ_j x_2j+1 · exp(-2πi·(2j+1)k/n) = E_k + exp(-2πi·k/n) · O_k

E_k is the length-n/2 DFT of the even samples and O_k is the DFT of the odd samples. The factor exp(-2πi·k/n) — the "twiddle factor" — is the only glue you need to stitch two half-size transforms into a full one. Because the half-size transforms are periodic, one pair (E_k, O_k) produces two outputs: X_k = E_k + w·O_k and X_{k+n/2} = E_k − w·O_k, the "butterfly".

Recurse until the transforms are length 1 (which are trivially themselves), and the cost solves to O(n log n): log₂ n stages, each doing O(n) butterflies. That is the whole difference between quadratic and near-linear — the algorithm that made the FFT the most important numerical routine of the 20th century.

Equana's kernel does this iteratively rather than recursively. It first applies the bit-reversal permutation — the recursion's even/odd splitting, unrolled, reorders the input so index j moves to the position given by reversing its bits — and then runs the butterflies in-place, doubling the transform length each stage from 2 up to n:

for len = 2, 4, 8, … , n: # each stage doubles the block size w = exp(-2πi / len) # stage twiddle for each block of size len: combine the two half-blocks with butterflies

No recursion, no scratch allocation, cache-friendly linear passes over one array. The catch is right there in the loop: len only ever takes power-of-two values. Feed this kernel a length of 1000 and the doubling never lands on it.

Bluestein: turn any length into a convolution

Here is where the second algorithm earns its keep. Bluestein's trick starts from a piece of algebra that looks like nothing at first. Use the identity j·k = (j² + k² − (k−j)²) / 2 inside the exponent of the DFT. Substituting and factoring out the parts that depend only on k:

X_k = exp(-πi·k²/n) · Σ_j [ x_j · exp(-πi·j²/n) ] · exp(+πi·(k−j)²/n)

Stare at the sum. The term in brackets depends only on j; call it a_j. The last factor depends only on the difference k − j; call it b_{k−j}. A sum over j of a_j · b_{k−j} is, by definition, a convolution. So the DFT — for any n, prime or not — is:

1. Multiply the input by a "chirp", exp(-πi·j²/n), whose frequency ramps quadratically (that's the chirp-z name).
2. Convolve the result with the conjugate chirp.
3. Multiply by the chirp once more.

Steps 1 and 3 are cheap element-wise passes. Step 2 is the interesting one, because a convolution has its own fast route: by the convolution theorem, a convolution is a pointwise product in the frequency domain. And here is the payoff — that convolution can be any length you like, as long as it's at least 2n−1. So you zero-pad both sequences up to the next power of two, and now the fast radix-2 engine can do the transforms.

That is exactly what Equana's kernel does. Given a length n that isn't a power of two, it picks the smallest power of two m ≥ 2n + 1, builds the chirp c_k = exp(-πi·k²/n) (computing k² mod 2n first, so the angle stays accurate for large k), forms the two padded sequences, and runs three radix-2 FFTs — forward on each sequence, pointwise-multiply, inverse — to get the circular convolution. A final chirp multiply recovers the spectrum:

# The user never sees any of this — one function, both algorithms x = randn(1000) # length 1000 = 2^3 · 5^3, NOT a power of two X = fft(x) # routes to Bluestein automatically y = ifft(X) # round-trips: y ≈ x norm(real(y) - x) # ~1e-13 p = randn(1024) # a power of two P = fft(p) # routes to radix-2 Cooley–Tukey instead

The cost accounting is worth appreciating. Bluestein does three FFTs of length m ≈ 2n, so it is a constant factor slower than a native power-of-two FFT of the same size — but it is still O(n log n), not O(n²). A length-1009 prime transform that would be a million-operation quadratic slog becomes three ~2048-point radix-2 FFTs plus a few linear passes. Arbitrary lengths stop being a special case and become just… slightly more work.

One function, two engines

The dispatch is a single branch on whether the length is a power of two — a two-line check n & (n − 1) === 0 — sending the transform to fftRadix2 or fftBluestein. From the notebook there is no seam: fft([…]) returns the same DFT coefficients X_k = Σ_j x_j exp(-2πi·j·k/n) whether the length is 1024 or 1000, and ifft inverts either one through the conjugation identity ifft(x) = conj(fft(conj(x))) / n. The companion functions fftfreq, fftshift, and ifftshift handle the bookkeeping of turning bin indices into real frequencies.

The pattern here is one I keep coming back to when building numerical code: make the common case fast, then find the algebraic bridge that reduces the general case to the common case rather than writing a second implementation from scratch. Bluestein is that bridge for the FFT — a reduction from "arbitrary length" to "power-of-two convolution" — and it means the library never has to tell a user "sorry, record a power-of-two number of samples next time."

This is the transform behind the applied walkthrough in FFT signal analysis in the browser, where a noisy signal gets decomposed into its clean tones in ten lines. And if convolution-as-a-primitive interests you, the same "spend FLOPs only where they matter" instinct runs through Equana's sparse matrices. Everything runs in a browser tab at equana.dev.