Cooley–Tukey FFT Algorithm

Cooley–Tukey FFT Algorithm is the foundational divide-and-conquer algorithm in computation and algorithms for computing the Discrete Fourier Transform (DFT) of a complex sequence of composite length N = N₁ N₂ and its inverse in O(N log N) arithmetic operations. Published by American mathematicians James W. Cooley (of IBM Research) and John W. Tukey (of Princeton University and Bell Labs) in Mathematics of Computation in April 1965—and later discovered by historical scholarship (Heideman, Johnson, & Burrus, 1984) to have been anticipated in an unpublished 1805 astronomical manuscript by Carl Friedrich Gauss—the Cooley–Tukey algorithm is the dominant realization of the Fast Fourier Transform across scientific computing, digital signal processing, spectral partial differential equation solvers, and large-integer polynomial multiplication.

Within computational complexity theory, direct evaluation of the N-point Discrete Fourier Transform from its matrix-vector definition requires N² complex multiplications and N(N − 1) complex additions. For a signal of N = 2²⁰ ≈ 1.05 × 10⁶ samples, direct O(N²) evaluation demands over 10¹² operations, whereas the Radix-2 Cooley–Tukey decomposition requires only (N/2) log₂ N ≈ 1.05 × 10⁷ complex multiplications—a 100,000-fold reduction that made real-time spectral analysis, OFDM wireless communications (Wi-Fi, 5G LTE), medical MRI reconstruction, and fast integer multiplication for the Miller–Rabin primality test computationally practical.

The Discrete Fourier Transform and Primitive Roots of Unity

Given a sequence of N complex numbers x₀, x₁, ..., x_N−1 ∈ ℂ, the N-point Discrete Fourier Transform (DFT) maps the time-domain (or spatial-domain) vector x to the frequency-domain spectrum X₀, X₁, ..., X_N−1 ∈ ℂ defined by:

where ω_N = exp(−2πi / N) = cos(2π/N) − i sin(2π/N) is the principal primitive N-th root of unity in ℂ, satisfying the periodicity identity ω_N^m+N = ω_N^m, the half-period negation identity ω_N^N/2 = −1 (for even N), and the subgroup decimation identity ω_N^d·m = ω_N/d^m whenever d divides N.

General Mixed-Radix Factorization (N = N₁ · N₂)

Suppose the transform length N is composite: N = N₁ · N₂ with N₁, N₂ > 1. The Cooley–Tukey algorithm re-indexes the one-dimensional time index n ∈ {0, ..., N − 1} and the one-dimensional frequency index k ∈ {0, ..., N − 1} as two-dimensional mixed-radix coordinates (n₁, n₂) and (k₁, k₂):

Substituting n · k = (N₂ n₁ + n₂)(k₁ + N₁ k₂) = N₂ n₁ k₁ + n₂ k₁ + N₁ n₂ k₂ + N n₁ k₂ into ω_N^nk and applying ω_N^{N n₁ k₂} = 1^n₁k₂ = 1 along with ω_N^N₂ = ω_N1 and ω_N^N₁ = ω_N2 separates the double sum into nested inner and outer DFTs:

Thus, an N-point DFT reduces to three structured matrix steps:

1. Compute N₂ inner DFTs, each of length N₁ (cost: N₂ · T(N₁)).
2. Multiply each intermediate element (n₂, k₁) by the complex phase rotation factor ω_N^n₂k₁, called by Gentleman and Sande (1966) a twiddle factor (cost: N multiplications).
3. Compute N₁ outer DFTs, each of length N₂ (cost: N₁ · T(N₂)).

Recursively applying this factorization yields the recurrence T(N₁ N₂) = N₂ T(N₁) + N₁ T(N₂) + O(N), which solves to T(N) = O(N ∑ p_i) = O(N log N) whenever N factors into small prime radices.

The Danielson–Lanczos Lemma and Radix-2 DIT Butterfly

In the most common special case where N = 2^m is a power of two, choosing N₂ = 2 and N₁ = N/2 yields the Radix-2 Decimation-in-Time (DIT) algorithm (independently derived for X-ray crystallography by G. C. Danielson and Cornelius Lanczos in 1942). Splitting the input x_n into its even-indexed samples x_2r and its odd-indexed samples x_2r+1 (r = 0 .. N/2 − 1) gives:

where E_k is the (N/2)-point DFT of the even samples and O_k is the (N/2)-point DFT of the odd samples. Because both E_k and O_k are periodic with period N/2 (E_k+N/2 = E_k and O_k+N/2 = O_k) while the twiddle factor satisfies ω_N^k+N/2 = −ω_N^k, a single complex multiplication t = ω_N^k O_k simultaneously yields both frequency outputs X_k and X_k+N/2 via the Radix-2 DIT Butterfly:

Unrolling this recursion across m = log₂ N stages performs (N/2) log₂ N butterflies—totaling (N/2) log₂ N complex multiplications and N log₂ N complex additions.

Gentleman–Sande Decimation-in-Frequency (DIF) and Bit-Reversal Permutation

Choosing the dual factorization N₁ = 2 and N₂ = N/2 yields the Radix-2 Decimation-in-Frequency (DIF) algorithm (Gentleman & Sande, 1966), which folds the first half x_n and second half x_n+N/2 of the input array first, and applies the twiddle factor ω_Nⁿ to the difference before recursing on the even- and odd-indexed frequencies:

In an in-place iterative implementation of Radix-2 DIT (or DIF), repeatedly separating even and odd indices permutes the array indices by reversing the binary representation (b_m−1 b_m−2 ··· b₁ b₀)₂ ↦ (b₀ b₁ ··· b_m−2 b_m−1)₂, known as the bit-reversal permutation. Notably, connecting a DIF forward transform (natural input → bit-reversed frequency output) to a DIT inverse transform (bit-reversed frequency input → natural output) for fast convolution via the Convolution Theorem cancels both bit-reversal passes!

Split-Radix FFT and Cache-Oblivious / FFTW Architecture

In 1984, Pierre Duhamel and Henk Hollmann introduced the Split-Radix FFT, which applies a Radix-2 decomposition to the even-indexed frequencies X_2k and a Radix-4 decomposition to the odd-indexed frequencies X_4k+1 and X_4k+3:

Split-Radix reduces the total real arithmetic count (real additions + real multiplications) from 5 N log₂ N to 4 N log₂ N − 6N + 8 (further refined to (34/9) N log₂ N by Johnson and Frigo in 2007 via conjugate-pair factor scaling).

On modern multi-core CPUs and GPUs, memory hierarchy bandwidth rather than raw floating-point count dominates runtime. David H. Bailey's (1990) Four-Step / Six-Step FFT chooses N₁ ≈ N₂ ≈ √N so that each √N-length sub-transform fits fully inside L1/L2 cache, while FFTW (Matteo Frigo and Steven G. Johnson, 2005) uses a dynamic-programming planner and codelet compiler (genfft) to select optimal mixed-radix Cooley–Tukey factorizations tailored to the target processor's SIMD vector registers (AVX-512, ARM NEON) alongside divide-and-conquer algorithms such as merge sort.

Carl Friedrich Gauss (1805) and the 1965 Nuclear Test Ban Context

Although the algorithm is named after Cooley and Tukey's 1965 publication—which arose when Richard Garwin on President Kennedy's Science Advisory Committee asked Tukey for a fast method to analyze seismic sensor arrays monitoring subterranean Soviet nuclear tests, and Cooley programmed and generalized it at IBM—Carl Friedrich Gauss had derived the precise same mixed-radix factorization N = N₁ N₂ in 1805 ("Theoria interpolationis methodo nova tractata", published posthumously in Latin in 1866) to interpolate the periodic orbital ephemeris of the asteroid 2 Pallas from 12 sample observations (N = 12 = 3 × 4).

How does the Cooley–Tukey FFT compute the Inverse Discrete Fourier Transform (IDFT) using the precise same code?

The Inverse Discrete Fourier Transform x_n = (1/N) ∑_k=0..N−1 X_k ω_N^−nk differs from the forward DFT only in the sign of the exponent (+2πi/N instead of −2πi/N) and the final 1/N normalization factor. Taking the complex conjugate of X_k, running the unmodified forward Cooley–Tukey FFT, taking the complex conjugate of the result, and dividing by N (or simply swapping the real and imaginary parts of input and output) computes the precise IDFT in O(N log N) time.

What happens if the input length N is a prime number where N = N₁ N₂ cannot be factored?

If N is prime, the Cooley–Tukey index factorization N = N₁ N₂ does not apply directly. Instead, Bluestein's Chirp-Z algorithm (1970) rewrites the product 2*nk* = k² + n² − (k − n)² to express the prime-length DFT as a cyclic convolution, pads that convolution to the next power of two M ≥ 2N − 1, and evaluates it in O(N log N) time using three power-of-two Cooley–Tukey FFTs!

1. Cooley, J. W., & Tukey, J. W. (1965). "An algorithm for the machine calculation of complex Fourier series." Mathematics of Computation, 19(90), 297–301. https://doi.org/10.1090/S0025-5718-1965-0178586-1
2. Gentleman, W. M., & Sande, G. (1966). "Fast Fourier transforms: For fun and profit." Proceedings of the 1966 Fall Joint Computer Conference (AFIPS '66), Vol. 29, 563–578. https://doi.org/10.1145/1464291.1464352
3. Duhamel, P., & Hollmann, H. (1984). "Split radix FFT algorithm." Electronics Letters, 20(1), 14–16. https://doi.org/10.1049/el:19840012
4. Heideman, M. T., Johnson, D. H., & Burrus, C. S. (1984). "Gauss and the history of the fast Fourier transform." IEEE ASSP Magazine, 1(4), 14–21. https://doi.org/10.1109/MASSP.1984.1162257
5. Frigo, M., & Johnson, S. G. (2005). "The design and implementation of FFTW3." Proceedings of the IEEE, 93(2), 216–231. https://doi.org/10.1109/JPROC.2004.840301