Prev
Unit 7 · Interpolation and Approximation
Trigonometric Interpolation and the FFT
65m100 XP
Next
📈Interpolation and Approximation

Trigonometric Interpolation and the FFT

The discrete Fourier transform, the fast Fourier transform, spectral analysis, leakage and aliasing

1963–1965: listening for bombs
The partial test-ban talks between the United States, the United Kingdom and the Soviet Union stumble on one technical question: atmospheric tests can be detected from afar, but can an underground explosion be told apart from an earthquake, using seismometers outside the Soviet Union? The answer lies in the frequency content of the recorded waves — and computing the spectrum of a long record with the obvious formula takes N2N^2 operations. John Tukey, a Princeton statistician serving on President Kennedy's Science Advisory Committee, sketches a way to compute the transform in far fewer steps. The physicist Richard Garwin, who wants it for exactly this purpose, takes the idea to James Cooley, a programmer at IBM's research centre, and asks him to make it work. In 1965 Cooley and Tukey publish An Algorithm for the Machine Calculation of Complex Fourier Series — five pages. For N=106N = 10^6 samples the fast algorithm is about 100,000 times faster. It changed signal processing, communications, imaging and computing itself. Years later, historians found that Gauss had derived the same idea in 1805 while interpolating asteroid orbits — in a manuscript published only after his death.
For a periodic signal, the natural interpolating functions are sines and cosines. Sample a signal at NN equally spaced times x0,…,xN−1x_0, \dots, x_{N-1} over one period. The discrete Fourier transform (DFT) gives the coefficients of the trigonometric polynomial through those samples:
Xk=∑n=0N−1xn e−2πikn/N,xn=1N∑k=0N−1Xk e2πikn/NX_k = \sum_{n=0}^{N-1} x_n\,e^{-2\pi i kn/N}, \qquad x_n = \frac{1}{N}\sum_{k=0}^{N-1}X_k\,e^{2\pi i kn/N}

Xₖ measures how much of frequency k cycles per record is in the signal. With sampling rate fs, bin k corresponds to f = k fs/N; the resolution is fs/N = 1/(record length). Computed directly: N² complex multiplications.

Derivation
Why the inverse formula undoes the transform
Complex exponentials. Euler's formula eiθ=cos⁡θ+isin⁡θe^{i\theta} = \cos\theta + i\sin\theta packs a cosine and a sine into one number, so e2πikn/Ne^{2\pi ikn/N} is a sampled wave that makes kk complete turns over the NN samples. Working with them instead of separate sines and cosines makes the algebra short.
The key sum. Fix two indices nn and mm and add the waves over all kk. With r=e2πi(n−m)/Nr = e^{2\pi i(n - m)/N} the sum is a geometric series:
∑k=0N−1e2πik(n−m)/N=∑k=0N−1rk={Nn=m1−rN1−r=0n≠m\sum_{k=0}^{N-1} e^{2\pi ik(n - m)/N} = \sum_{k=0}^{N-1} r^k = \begin{cases} N & n = m \\[2pt] \dfrac{1 - r^N}{1 - r} = 0 & n \ne m \end{cases}
since rN=e2πi(n−m)=1r^N = e^{2\pi i(n - m)} = 1. Different discrete frequencies are orthogonal.
Inverting. Substitute the definition of XkX_k into the inverse formula and swap the order of the sums:
1N∑kXke2πikn/N=1N∑mxm∑ke2πik(n−m)/N=1N xn N=xn.\frac1N\sum_k X_ke^{2\pi ikn/N} = \frac1N\sum_m x_m\sum_k e^{2\pi ik(n - m)/N} = \frac1N\,x_n\,N = x_n.
So the inverse formula reproduces every sample: it is a trigonometric polynomial passing through the data, and the XkX_k are its coefficients — trigonometric interpolation.
💡The fast Fourier transform
Split the sum into even and odd samples: Xk=∑mx2me−2πikm/(N/2)⏟Ek+e−2πik/N∑mx2m+1e−2πikm/(N/2)⏟Ok.X_k = \underbrace{\sum_{m} x_{2m}e^{-2\pi i km/(N/2)}}_{E_k} + e^{-2\pi i k/N}\underbrace{\sum_{m}x_{2m+1}e^{-2\pi i km/(N/2)}}_{O_k}. EkE_k and OkO_k are DFTs of length N/2N/2, and since they repeat with period N/2N/2, Xk+N/2=Ek−e−2πik/NOkX_{k+N/2} = E_k - e^{-2\pi ik/N}O_k. One transform of length NN becomes two of length N/2N/2 plus N/2N/2 multiplications by 'twiddle factors'. Recurse down to length 1: about N2log⁡2N\tfrac{N}{2}\log_2 N multiplications instead of N2N^2 — 5,120 instead of a million for N=1024N = 1024, and 10710^7 instead of 101210^{12} for N=220N = 2^{20}. numpy.fft uses refined versions that work for any NN, fastest when NN has only small prime factors.
Derivation
Counting the FFT's multiplications
Why the halves repeat. EkE_k is a DFT of length N/2N/2, so Ek+N/2=EkE_{k+N/2} = E_k (adding N/2N/2 to kk adds a whole number of turns to every exponential); the same for OkO_k. The twiddle factor changes sign, e−2πi(k+N/2)/N=−e−2πik/Ne^{-2\pi i(k + N/2)/N} = -e^{-2\pi ik/N}. That is why one pair Ek,OkE_k, O_k gives two outputs, XkX_k and Xk+N/2X_{k+N/2}, for one multiplication.
The count. Let C(N)C(N) be the multiplications for length N=2mN = 2^m. Two half-length transforms plus N/2N/2 twiddle multiplications:
C(N)=2 C(N/2)+N2⟹C(N)N=C(N/2)N/2+12.C(N) = 2\,C(N/2) + \frac N2 \qquad\Longrightarrow\qquad \frac{C(N)}{N} = \frac{C(N/2)}{N/2} + \frac12.
Each halving adds 12\tfrac12 to C(N)/NC(N)/N, and there are m=log⁡2Nm = \log_2N halvings down to length 1, which costs nothing. So C(N)/N=12log⁡2NC(N)/N = \tfrac12\log_2N:
C(N)=N2log⁡2N.C(N) = \frac N2\log_2N.
💡Aliasing and the Nyquist frequency
Sampled at rate fsf_s, the signals cos⁡(2πft)\cos(2\pi ft) and cos⁡(2π(f−kfs)t)\cos(2\pi(f - kf_s)t) have identical samples. Frequencies above the Nyquist frequency fs/2f_s/2 fold back into [0,fs/2][0, f_s/2]: a 70 Hz vibration sampled at 100 Hz looks exactly like a 30 Hz one, and so do 130 and 170 Hz. No processing after sampling can undo it. The sampling theorem (Nyquist 1928, Shannon 1949): a signal containing no frequencies above fmaxf_{max} is determined completely by samples taken faster than 2fmax2f_{max}. That is why every data-acquisition system has an analogue anti-aliasing filter before the converter.
Derivation
Why two frequencies can share the same samples
Sampling at rate fsf_s means looking only at the times tn=n/fst_n = n/f_s. Shift the frequency by a whole multiple of fsf_s:
cos⁡(2π(f−kfs)tn)=cos⁡(2πfnfs−2πkn)=cos⁡(2πfnfs),\cos\big(2\pi(f - kf_s)t_n\big) = \cos\Big(2\pi f\frac{n}{f_s} - 2\pi kn\Big) = \cos\Big(2\pi f\frac{n}{f_s}\Big),
because knkn is a whole number of turns. The two signals differ between the samples but agree at every one of them. Since cosine is even, ff and fs−ff_s - f also agree: cos⁡(2π(fs−f)n/fs)=cos⁡(2πn−2πfn/fs)=cos⁡(2πfn/fs)\cos(2\pi(f_s - f)n/f_s) = \cos(2\pi n - 2\pi fn/f_s) = \cos(2\pi fn/f_s). So every frequency has a twin in [0,fs/2][0, f_s/2] — the 70 Hz vibration and its 30 Hz alias at fs=100f_s = 100 Hz — and from samples alone the two are indistinguishable.
📐Listening to a gearbox
An accelerometer on a gearbox, sampled at 1000 Hz for 4 s (N=4000N = 4000, resolution 0.25 Hz), records a shaft vibration at 23.6 Hz and a ten-times-weaker gear-mesh tone at 147.3 Hz. Neither falls on a bin. With no window the FFT finds peaks at 23.50 and 147.25 Hz, but with amplitudes of only 0.76 and 0.093 (true: 1 and 0.1), and energy leaks across the spectrum — 1.5×10−31.5 \times 10^{-3} at 60 Hz, where there is nothing. A Hann window, which tapers the record smoothly to zero at both ends, gives amplitudes 0.90 and 0.097 and cuts the leakage at 60 Hz by a factor of 17. Windowing trades a slightly wider peak for a much cleaner spectrum — essential when a weak tone sits near a strong one.
mathematical analysis is as extensive as nature itself; it defines all perceptible relations, measures times, spaces, forces, temperatures
– Joseph Fourier, Théorie analytique de la chaleur (1822)
🕰Gauss, again
In 1805 Gauss needed to interpolate the orbits of the newly discovered asteroids Pallas and Juno from a handful of observations — trigonometric interpolation of a periodic function. His notes describe splitting a transform of size N=N1N2N = N_1N_2 into smaller ones, the essential idea of the FFT. They were published only in 1866, in Latin, in his collected works, and were noticed by Heideman, Johnson and Burrus in 1984. Ceres opened Unit 5; the asteroids found after it close this unit.
💡Where you will use this
Structures and Rotating machinery: modal testing and condition monitoring are spectral analysis of accelerometer records. GNSS Navigation: receivers acquire satellite signals by FFT-based correlation over code phase and Doppler. Signals and Systems, Power Electronics: harmonic analysis of converter waveforms. Numerical PDEs (Unit 10): spectral methods solve periodic problems by FFT, and fast Poisson solvers use it for whole grids.
📚Further reading
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 12, 'Fast Fourier Transform'.
  • James Cooley and John Tukey, 'An algorithm for the machine calculation of complex Fourier series', Mathematics of Computation 19(90), 297–301 (1965).
  • Michael Heideman, Don Johnson and C. Sidney Burrus, 'Gauss and the history of the fast Fourier transform', IEEE ASSP Magazine 1(4), 14–21 (1984).
  • Steven W. Smith, The Scientist and Engineer's Guide to Digital Signal Processing (California Technical Publishing, 1997) — sampling, aliasing, windows and the FFT for engineers; free online.
↓But there is a problem
This unit was about knowing a function between the points where it was measured — by polynomials, splines and sines. With an interpolant in hand, we can do more than evaluate it. A rocket motor's thrust is recorded every millisecond on a test stand. The number the customer pays for is the total impulse — the integral of that curve. A structural test records displacement at a few gauges; the engineer needs strain — a derivative. Both are calculus applied to numbers. Differentiation, it turns out, is treacherous with measured data, and integration is forgiving. The next unit explains why — and how to squeeze astonishing accuracy out of a few function values.