Prev
Unit 6 · Lesson 6-2
The DFT and the FFT
1h 55m1300 XP
Next

The DFT and the FFT

The DFT as samples of the DTFT, frequency bins, the DFT matrix, circular convolution, the radix-2 FFT and N² against N log N

1963 — seismometers and a test ban
A meeting of President Kennedy's Science Advisory Committee turns to a hard question: could seismometers placed around the Soviet Union, outside its borders, reveal its nuclear tests? Analysing their records would take an enormous number of Fourier transforms. John Tukey, a statistician at the meeting, comes up with a faster way to compute them (Fast Fourier transform). In discussion with Tukey, Richard Garwin of IBM sees that the method matters far beyond national security — he needs it himself, to find periodicities in a crystal of helium-3 — and hands the idea to James Cooley at IBM's Watson laboratory. Cooley and Tukey publish it within six months, in 1965. Because Tukey did not work for IBM, the method's patentability was doubtful, and it went into the public domain, where the computing revolution made it one of the indispensable algorithms of signal processing (Fast Fourier transform). This lesson is the transform a computer actually computes — the DFT — and the fast Fourier transform that computes it.
The DTFT of an NN-sample record is a continuous function of Ω\Omega. The discrete Fourier transform keeps NN samples of it, equally spaced around the circle: Ωk=2πk/N\Omega_k = 2\pi k/N, k=0,1,…,N−1k = 0, 1, \dots, N - 1. That is just enough: NN numbers in, NN numbers out, and the inverse transform gives the record back exactly.
X[k]=∑n=0N−1x[n] e−j2πkn/N,x[n]=1N∑k=0N−1X[k] ej2πkn/NX[k] = \sum_{n=0}^{N-1} x[n]\,e^{-j2\pi kn/N}, \qquad x[n] = \frac1N\sum_{k=0}^{N-1} X[k]\,e^{j2\pi kn/N}

The DFT and its inverse: N samples of the DTFT, and the N-term sum that rebuilds the record.

💡Frequency bins
Bin kk sits at Ωk=2πk/N\Omega_k = 2\pi k/N, which for samples taken at fsf_s is the frequency fk=k fsN.f_k = \frac{k\,f_s}{N} . The bins are fs/Nf_s/N apart — the reciprocal of the record's duration. Bins above N/2N/2 are the negative frequencies: bin N−kN - k is frequency −fk-f_k. For a real record X[N−k]=X∗[k]X[N - k] = X^*[k], so the second half adds nothing, and numpy.fft.rfft computes only the first. A one-second record sampled at 8 kHz has 8000 bins, 1 Hz apart; a 1 kHz tone lands in bin 1000, and its mirror image in bin 7000.
A matrix. Written out, the DFT is a matrix product X=Fx\mathbf{X} = \mathbf{F}\mathbf{x}, with Fkn=WknF_{kn} = W^{kn} and W=e−j2π/NW = e^{-j2\pi/N}. Its columns are orthogonal — FHF=NI\mathbf{F}^H\mathbf{F} = N\mathbf{I} — so the inverse is FH/N\mathbf{F}^H/N. Computed this way, the DFT costs N2N^2 complex multiplications: about 101210^{12} for a million samples.
Derivation
The radix-2 FFT
Let NN be a power of two and split the sum into even and odd samples, n=2mn = 2m and n=2m+1n = 2m + 1:
X[k]=∑m=0N/2−1x[2m] W2mk+Wk∑m=0N/2−1x[2m+1] W2mk.X[k] = \sum_{m=0}^{N/2-1} x[2m]\,W^{2mk} + W^k\sum_{m=0}^{N/2-1} x[2m+1]\,W^{2mk} .
Since W2=e−j2π/(N/2)W^2 = e^{-j2\pi/(N/2)}, each sum is an N/2N/2-point DFT: X[k]=E[k]+WkO[k]X[k] = E[k] + W^k O[k], where EE and OO are the transforms of the even and odd samples. Both repeat with period N/2N/2, and Wk+N/2=−WkW^{k + N/2} = -W^k, so
X[k]=E[k]+WkO[k],X[k+N/2]=E[k]−WkO[k],k=0,…,N/2−1.X[k] = E[k] + W^k O[k], \qquad X[k + N/2] = E[k] - W^k O[k], \qquad k = 0, \dots, N/2 - 1 .
One multiplication by WkW^k serves two outputs — a butterfly. Combining the halves costs N/2N/2 multiplications; splitting again and again gives log⁡2N\log_2 N stages, and N2log⁡2N\tfrac N2\log_2 N multiplications in all.
📐How much faster
  • N=1024N = 1024: the direct sum needs 10242≈1.051024^2 \approx 1.05 million multiplications; the FFT 512×10=5120512 \times 10 = 5120 — 205 times fewer.
  • N=220N = 2^{20}, about a million samples: 1.1×10121.1 \times 10^{12} against 1.05×1071.05 \times 10^7 — a hundred thousand times fewer. A transform that would take a day takes about a second.
The saving grows with NN, which is why the FFT made spectral analysis of long records routine.
💡Circular convolution
The DFT treats a record as one period of a periodic sequence. Multiplying two DFTs therefore computes a circular convolution, in which the tail of the result wraps round onto its start: y[n]=∑m=0N−1x[m] h[(n−m) mod N]    ↔    Y[k]=X[k] H[k].y[n] = \sum_{m=0}^{N-1} x[m]\,h[(n - m) \bmod N] \;\;\leftrightarrow\;\; Y[k] = X[k]\,H[k] . Zero-pad both records to at least L+M−1L + M - 1 samples and nothing wraps: the result is the ordinary linear convolution (the lab of Lesson 4-2). Done this way, a long convolution costs on the order of Nlog⁡NN\log N operations instead of N2N^2 — fast convolution.
the most important numerical algorithm of our lifetime
– Gilbert Strang on the FFT (1994), quoted in Fast Fourier transform
🕰Gauss got there first
In 1805, interpolating the orbits of the asteroids Pallas and Juno from observations, Carl Friedrich Gauss used a method very similar to Cooley and Tukey's. He never published it — it appeared only in his collected works — and its connection to the FFT was recognised in 1984, when Heideman, Johnson and Burrus showed that Cooley and Tukey had independently rediscovered an algorithm Gauss knew around 1805 (Fast Fourier transform).
⚠The DFT thinks your record repeats
Every DFT assumes the record is one period of a periodic signal. If its two ends do not join smoothly — and for a real recording they almost never do — the jump at the seam is part of what gets transformed, and it spreads energy across the spectrum. That spreading is leakage, the subject of Lesson 6-3; zero padding does not remove it.
💡Where you will use this
Spectrum analysers, audio codecs, the OFDM modulation of Wi-Fi and mobile networks, and MRI reconstruction all run on FFTs. Lesson 6-3 reads DFT spectra honestly, Lesson 6-4 strings them together into spectrograms, and the correlations of Lesson 10-3 are computed as fast convolutions.
📚Further reading
  • James W. Cooley and John W. Tukey, "An Algorithm for the Machine Calculation of Complex Fourier Series", Mathematics of Computation 19 (1965), 297–301.
  • Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time Signal Processing (3rd ed., Pearson, 2010), chapters 8 and 9.
  • Fast Fourier transform on Wikipedia.
↓But there is a problem
Put a tone exactly on a bin and the DFT shows one clean line. Move it half a bin and it lights up every bin in the spectrum — and a faint tone nearby can vanish in the glow. In 1978 a San Diego engineer published a long paper on how to stop it, and it became the most famous thing he ever wrote.