Fast Fourier Transform (FFT) theory

The Discrete Fourier Transform converts a discrete complex sequence x_n of length N into its orthogonal frequency components X_k via the summation X_k = Σ_{n=0}^{N-1} x_n W_N^{kn}, where W_N = exp(-2πi/N) is the primitive N-th root of unity. In matrix notation, this maps to multiplying an N-dimensional vector by a dense N×N Vandermonde matrix F_N whose entries are (F_N)_{k,n} = W_N^{kn}. Because every entry in F_N is non-zero, computing all N output coefficients directly requires N independent inner products, each performing N complex multiplications and N-1 additions, totaling exactly N² complex multiplications and N(N-1) complex additions, establishing a strict Θ(N²) quadratic computational complexity.

Cooley-Tukey Radix-2 Decimation-in-Time (DIT) partitions the N-point summation into two (N/2)-point summations by separating even-indexed samples (n = 2m) from odd-indexed samples (n = 2m + 1): X_k = Σ_{m=0}^{N/2-1} x_{2m} W_N^{2mk} + W_N^k Σ_{m=0}^{N/2-1} x_{2m+1} W_N^{2mk}. Applying the reduction property W_N^{2mk} = W_{N/2}^{mk}, both inner summations become independent (N/2)-point DFTs, denoted E_k (even) and O_k (odd). Furthermore, using the half-period anti-symmetry property W_N^{k + N/2} = -W_N^k, the algorithm computes X_k = E_k + W_N^k O_k and X_{k + N/2} = E_k - W_N^k O_k for 0 ≤ k < N/2. This recursive formulation obeys the recurrence T(N) = 2T(N/2) + Θ(N), which solves via the Master Theorem to Θ(N log₂ N).

A Radix-2 Butterfly Operation is the fundamental arithmetic building block of the Cooley-Tukey algorithm, named after the visual resemblance of its intersecting signal flow lines. It takes two complex sub-transform inputs E_k and O_k and a complex twiddle factor W_N^k, producing two updated outputs: an upper output X_k = E_k + W_N^k O_k and a lower output X_{k + N/2} = E_k - W_N^k O_k. The product T = W_N^k O_k is evaluated once via a single complex multiplication (4 real multiplications and 2 real additions). Then, T is added to E_k to yield the upper branch and subtracted from E_k to yield the lower branch (2 complex additions, or 4 real additions). Thus, a single butterfly requires 1 complex multiplication and 2 complex additions.

In the decimation-in-time decomposition, input elements are recursively separated into even and odd sets at each level of the tree. An element's routing decisions down the recursion tree depend on its binary index bits from least significant to most significant: even indices (bit 0 = 0) route left, odd indices (bit 0 = 1) route right; then bit 1 routes within sub-trees, and so on. At the leaf level of the recursion tree, the element originally at index i resides at an array position whose binary representation is the exact bit-reversal of i (e.g. index 001₂ moves to 100₂ = 4). By applying this bit-reversal permutation prior to butterfly execution, the algorithm can proceed entirely iteratively from stage 1 to log₂(N) using contiguous in-place array updates with zero call-stack overhead and O(1) auxiliary memory.

Multiplying two polynomials A(x) = Σ a_j x^j and B(x) = Σ b_k x^k of degree bound n yields a product polynomial C(x) = A(x)B(x) whose coefficients c_p = Σ a_j b_{p-j} represent the discrete linear convolution of coefficient vectors a and b. Direct polynomial multiplication takes O(n²) operations. According to the Convolution Theorem, circular convolution in the time domain is mathematically isomorphic to element-wise pointwise multiplication in the Fourier frequency domain: F{a * b} = F{a} ⊙ F{b}. By zero-padding the coefficient vectors to length N ≥ 2n - 1 (the next power of two) to prevent circular aliasing, the algorithm performs: (1) forward FFT on a and b in 2 · O(N log N) time; (2) pointwise multiplication C_freq[k] = A_freq[k] · B_freq[k] in O(N) time; and (3) inverse FFT on C_freq in O(N log N) time. The entire multiplication completes in O(N log N) operations.

Decimation-in-Time (DIT) and Decimation-in-Frequency (DIF) are the two dual formulations of the Cooley-Tukey algorithm. DIT decomposes the input sequence in the time domain by splitting even and odd indices recursively, requiring a bit-reversed input array and producing normally ordered frequency outputs. Its butterfly multiplies the odd input by the twiddle factor before addition and subtraction. In contrast, DIF decomposes the output spectrum in the frequency domain, taking normally ordered inputs and producing bit-reversed frequency outputs. Its butterfly performs addition and subtraction first, and then multiplies the difference by the twiddle factor. Both architectures have identical computational complexity Θ(N log N) and can be coupled together to perform convolutions without any bit-reversal permutations (DIF forward FFT followed by DIT inverse FFT).

While Radix-2 decomposes sequences into pairs, Radix-4 decomposes sequences into 4 sub-sequences of length N/4. A Radix-4 butterfly processes 4 inputs simultaneously and requires only 3 complex multiplications instead of 4 in Radix-2, reducing total complex multiplications by 25% across the transform. Split-Radix FFT blends Radix-2 and Radix-4 decompositions, achieving the absolute theoretical minimum arithmetic operation count for power-of-two lengths. For hardware performance, libraries like FFTW (Frigo & Johnson) implement cache-oblivious recursive Cooley-Tukey, decomposing large transforms into small kernels that fit entirely within L1/L2 caches, eliminating TLB thrashing and maximizing SIMD instruction pipeline utilization.