Fast Fourier transform

A synth plays an arpeggio, and the picture shows which frequencies sound at each moment: a bright band is a loud frequency. To draw it, the browser cuts six seconds of sound into 513 windows of 1024 samples and splits each into frequencies. By the definition that is 1024 × 1024 multiplications per window, more than half a billion per picture. The fast Fourier transform does the same in 5120 multiplications per window, two hundred times fewer, and the picture is recomputed faster than you can let go of the button. Change the waveform, the filter, the window length.

A sine has one line in the picture: one frequency. A saw has a ladder: the fundamental and all the overtones that are multiples of it; a square wave only the odd ones, so its ladder is sparser. The FM bell sounds with overtones that are not multiples of the fundamental, so its ladder is crooked. A filter that opens and closes smoothly shows up as a wave of brightness along the top. And the window length is a choice between time and frequency: a window of 256 samples sees each note separately but smears the frequencies, a window of 4096 resolves frequencies down to 5 Hz and smears the notes into one another.

The idea

The discrete Fourier transform (DFT) of N samples x₀, …, x_N−1 is the N numbers

Look at this sum as a polynomial: P(z) = x₀ + x₁z + x₂z² + … + x_N−1z^N−1. Then X_k = P(ω^k), the value of the polynomial at the k-th N-th root of unity. The DFT is evaluating a polynomial at N points, and done head-on this takes N² multiplications.

Split the polynomial into even and odd powers: P(z) = E(z²) + z · O(z²), where E and O have half as many coefficients. The squares of the N roots of unity are only N/2 distinct numbers, the roots of degree N/2. So it is enough to evaluate E and O at N/2 points, which is two DFTs of half the length, and then assemble the answer: X_k = E_k + ω^kO_k and X_k+N/2 = E_k − ω^kO_k. One problem of size N has turned into two of size N/2 and N/2 multiplications.

Repeating the split until the length is one gives log₂ N levels of N/2 multiplications: N/2 · log₂ N. The minus sign in the second formula is because ω^N/2 = −1: the points z and −z have the same square, and one evaluation of E and O serves both. The derivation with a proof of invertibility and the story of Gauss, who came up with this technique around 1805, is in the chapter “Fourier series” of “Mathematics, the Queen of the Sciences”. Here we look at how the algorithm works in practice and what is done with it.

Butterflies

The recursion is conveniently unrolled. If at each level the even samples go left and the odd ones right, at the very bottom the samples end up in the order where the index is written in binary digits back to front: 0, 4, 2, 6, 1, 5, 3, 7. Rearrange them this way at once, and from then on the algorithm goes bottom up, level by level, merging neighbouring blocks with an operation that, after the shape of its diagram, is called a butterfly: the pair a, b becomes a + w·b and a − w·b.

Each circle is a complex number, and the arrow shows its direction. The input is four ones and four zeros. At the output X₀ = 4, the sum of all the samples, and X₂, X₄, X₆ are zero: a signal made of a block of ones does not contain these frequencies. Twelve butterflies, twelve multiplications instead of sixty-four, and the whole computation runs in place, in one array.

Code

The iterative version: a bit-reversal permutation, then levels of butterflies of length 2, 4, 8, … The inverse transform is the same with ω replaced by ω⁻¹ and a division by N. The program splits a signal made of four sinusoids and checks itself: on 128 samples it compares the answer with the DFT by the definition, and the inverse transform brings back the original signal.

n = 1024: 5120 multiplications in the FFT, 1048576 by the definition
bin   55 =    440 Hz, amplitude 1.000
bin  110 =    880 Hz, amplitude 0.500
bin  165 =   1320 Hz, amplitude 0.250
bin  250 =   2000 Hz, amplitude 0.300

The program gets the factors ω^k by multiplying again and again by the same complex number: this is shorter, but rounding errors pile up. In long transforms the factors are computed afresh with sine and cosine or taken from a table, as in the next two listings.

Now real data: the mean daily temperature of Almaty over 85 years, 31,046 days from the ERA5 reanalysis. What periods are hidden in it? Pad the series with zeros to 131,072 samples, more than four times as many, so that the frequency grid is fine, and look at the highest peaks of the spectrum.

31046 days, mean 8.76 C, FFT of 131072 samples: 1114112 butterflies instead of 17179869184 products
peak at bin 359: period 365.2 days, amplitude 13.2 C
peak at bin 4: period 33500.9 days, amplitude 1.0 C
peak at bin 718: period 182.6 days, amplitude 0.8 C

The main peak is 365.2 days: the FFT found the length of the year on its own, knowing nothing about the calendar. The amplitude of 13.2 °C means summer is warmer than average by about as much as winter is colder. The second peak is the longest wave that fits into 85 years: this is how slow changes in the mean temperature over decades look in a spectrum. The third is half a year, 182.6 days, with an amplitude of 0.8 °C: the yearly course is not a pure sinusoid, and the second harmonic describes its unevenness. Computing head-on would have taken 17 billion multiplications; the FFT manages with a million butterflies.

Multiplying polynomials and big numbers

Multiply (1 + 2x + 3x²)(4 + 5x + 6x²). The coefficient of x^k in the product is the sum of a_ib_j over all i + j = k:

This operation on sequences of coefficients is called convolution, and done head-on it costs n² multiplications. Substitute x = 10 and you get 4 + 130 + 2800 + 27,000 + 180,000 = 209,934 = 321 · 654: a number is a polynomial in ten, its digits are the coefficients, lowest first. Long multiplication of numbers is exactly a convolution of digits with carries.

The FFT makes convolution fast. A polynomial of degree n is uniquely determined by its values at n + 1 points. Evaluate both factors at the roots of unity with a forward FFT, multiply the values pairwise, which is n multiplications, and recover the coefficients of the product with an inverse FFT. Three transforms of N log N plus N multiplications instead of n².

9543 digits times 10142 digits: FFT of length 8192
the product has 19684 digits: 399553154987076…763468471600001
rounding error below 0.01: yes
equals the built-in multiplication: yes
multiplications: column method 10754961, FFT about 167936

The digits are grouped in threes: each coefficient is a number from 0 to 999. Taking more is risky: the coefficients of the product reach n · 999², and floating-point computation must tell them apart from neighbouring integers. The program checks how far the coefficients it found are from integers: here less than 0.01. For even longer numbers the error grows, and then one moves to arithmetic of remainders: the number-theoretic transform does the same with no rounding.

Below is a benchmark: long multiplication, our FFT and the built-in multiplication, int in Python and math/big in Go.

Long multiplication grows as a square, the FFT a little faster than linearly, and our FFT overtakes long multiplication in Python by a thousand digits, and in Go between a thousand and three thousand. But the built-in multiplication turns out faster than both: in CPython and in Go it is Karatsuba’s algorithm, asymptotically worse than the FFT at d^1.58, but written in C and assembly, and it wins by its constant on hundreds of thousands of digits. Libraries move to the Fourier transform where numbers are longer still: Python’s decimal module uses a number-theoretic transform for coefficients longer than 1024 machine words, and the GIMPS project squares numbers of tens of millions of digits with a special variant of the FFT.

What it costs

The running time satisfies T(N) = 2T(N/2) + O(N): two problems half the size and a linear merge. By the master theorem this is O(N log N), or more precisely N/2 · log₂ N multiplications and N log₂ N additions. Memory is O(N): the iterative version works in place.

Multiplying polynomials of degree n is three transforms of length N ≥ 2n + 1: O(n log n). For big numbers of d digits the carries are added, also linear, and with groups of fixed size this comes to O(d log d), as long as floating-point precision suffices.

Where the FFT trips up

• The length is a power of two. The simplest version needs N = 2^m. The signal is padded with zeros, as with the weather above, and for arbitrary lengths there are splittings by other factors and Bluestein’s algorithm; libraries like numpy.fft accept any length.
• The convolution is circular. The product of values at the roots of unity gives a convolution around a circle: the tail of the product wraps round to the start. To get an ordinary one, both factors are padded with zeros to a length of at least the sum of the lengths minus one.
• Spectral leakage. A window of a whole number of periods gives sharp peaks, as in the first listing. In real life a window cuts a wave anywhere, and the break at the edges splashes energy over neighbouring frequencies. Switch the window in the first scene to rectangular and the lines blur. The Hann window fades the edges smoothly and almost removes the splashes at the price of slightly wider peaks.
• Time against frequency. A window of length N tells frequencies apart in steps of (sampling rate)/N: the more precise the frequency, the longer the window and the less clear when exactly it sounded. This is not a flaw of the algorithm but a property of waves themselves.
• Rounding. In floating point the error grows with the length and with the size of the coefficients; for exact convolutions of integers there is the number-theoretic transform.

Where it works

• Communications. Wi-Fi of the 802.11a/g standards uses 52 subcarriers in a 20 MHz channel, 48 for data and 4 pilots: the transmitter assembles a symbol with an inverse 64-point FFT, and the receiver takes it apart with a forward one.
• Sound and images. Spectrum analysers, tuners, equalisers, convolution reverb, noise reduction. JPEG and MP3 compression rest on related cosine transforms, covered in the chapter “Fourier series”.
• Scientific computing. numpy.fft runs on the pocketfft library; solvers of partial differential equations, crystallography, telescope data processing.
• Long arithmetic. The search for Mersenne primes: the Lucas–Lehmer test squares numbers of tens of millions of digits, and the GIMPS programs do it with a variant of the FFT proposed by Richard Crandall and Barry Fagin in 1994. Python’s decimal module multiplies long coefficients with a number-theoretic transform.
• Olympiads. Convolutions of counts: in how many ways each sum can be made, matches of strings at all shifts; both problems below.

History

Gauss came up with this technique around 1805 while computing the orbits of asteroids, but it stayed unpublished during Gauss’s lifetime; more on this in “Mathematics, the Queen of the Sciences”. The algorithm’s second birth is tied to seismology. According to those involved, in 1963 the physicist Richard Garwin of IBM was listening at a meeting of science advisers to reports on detecting underground nuclear tests and noticed that the statistician John Tukey had a way to compute the Fourier transform much faster. Analysis of seismic signals came down to exactly this transform. Garwin found James Cooley at IBM, who wrote the program; in 1964 it ran on an IBM 7094, and in April 1965 the five-page paper by Cooley and Tukey, “An Algorithm for the Machine Calculation of Complex Fourier Series”, came out. In 1971 Arnold Schönhage and Volker Strassen built on the Fourier transform a multiplication of integers in O(n log n log log n); how this bound was brought down to n log n is told in the chapter “Divide and conquer” of the basic course.

Problems

In the first a convolution of counts answers all the questions about sums at once. In the second two strings have to be compared at all shifts, and for that the matching of letters has to be turned into a product of numbers.