Appearance
5.1 — The DFT and the FFT
The Fourier transform of Chapter 4.4 is an integral over infinite time. A computer has neither infinite time nor a continuous signal — it has N numbers. So the transform has to be rebuilt for that situation, and the result is the discrete Fourier transform.
Then comes the part that changed the world. Computing the DFT directly takes N^2 operations, which for a million points is 10^{12} multiplications — hours of work. In 1965 Cooley and Tukey published an algorithm that does it in N\log_2 N, which for a million points is 2\times10^7 — fifty thousand times faster. Almost every piece of digital signal processing that exists is downstream of that single result.
1. The DFT
\boxed{X[k]=\sum_{n=0}^{N-1}x[n]\,e^{-j2\pi kn/N}}, \qquad k=0,1,\ldots,N-1
x[n]=\frac1N\sum_{k=0}^{N-1}X[k]\,e^{j2\pi kn/N}
Read the forward transform aloud: "X of k equals the sum from n equals zero to N minus one of x-of-n times e to the minus j two pi k n over N."
The meaning is the same as always. e^{-j2\pi kn/N} is a probe — a complex sinusoid that completes exactly k cycles over the N samples. Multiplying by it and summing measures how much of that frequency the signal contains. It is a correlation with each candidate frequency in turn, exactly as in Chapter 4.4.
The quantity W_N = e^{-j2\pi/N} is called the twiddle factor, and it is a point on the unit circle. W_N^{kn} walks around that circle, and the whole DFT is a sum of the signal's samples each rotated by a different amount.
What the index k means in hertz
f_k = k\cdot\frac{f_s}{N}
The spacing between adjacent bins is f_s/N, called the bin width or frequency resolution.
Worked example. 1024 samples at 8 kHz:
\Delta f=\frac{8000}{1024}=7.81\ \text{Hz}
So bin 0 is DC, bin 1 is 7.81 Hz, bin 100 is 781 Hz, and so on.
The record length is what sets resolution:
\Delta f = \frac{f_s}{N}=\frac{1}{T_{record}}
To resolve two tones 1 Hz apart you need at least one second of signal. No amount of processing beats this — it is the time-frequency trade of Chapter 4.4 arriving as a hard number. Sampling faster does not help; you must record longer.
The symmetry, and where the useful bins stop
For a real input, the output has conjugate symmetry:
X[N-k]=X^*[k]
So the second half of the output is a mirror image of the first, carrying no new information. Bins 0 to N/2 are the useful ones; bin N/2 is Nyquist.
That halves the storage and the computation for real signals, and it is why an FFT of 1024 real samples is usually plotted as 513 points.
2. The three assumptions the DFT makes
The DFT is not a sampled version of the continuous Fourier transform. It is its own thing, and three properties of it cause every practical difficulty in this Part.
1. The signal is assumed periodic with period N. The DFT does not know your signal stops — it treats the N samples as one period of something that repeats forever. Section 3 is the consequence.
2. The frequency axis is also periodic. X[k+N]=X[k], which is the unit-circle wrapping of Chapter 4.6 restated.
3. Frequencies are quantised to the bin grid. A signal component that falls exactly on a bin is captured perfectly. One that falls between bins is not, and section 3 explains what happens instead.
3. Spectral leakage, and windowing
Take 1024 samples at 8 kHz of a 781.25 Hz sine. That is exactly bin 100 (100\times7.8125 = 781.25), so the sine completes exactly 100 whole cycles in the record. The assumed periodic extension joins up seamlessly. The result is one clean spike at bin 100 and essentially zero everywhere else.
Now make it 800 Hz. That is bin 102.4 — not an integer. The sine completes 102.4 cycles, so the periodic extension has a discontinuity where the end meets the beginning.
A discontinuity has broadband content (Chapter 4.4), so energy appears in every bin. Instead of one spike you get a smeared peak with sidelobes falling away slowly on both sides. This is spectral leakage, and it is not a numerical error — it is the correct transform of the signal the DFT was actually given, which was the discontinuous periodic extension.
The consequence that matters: a small signal near a large one gets buried under the large one's leakage. You cannot see a -60 dB tone 100 Hz away from a full-scale tone if the leakage at that offset is -30 dB.
The fix: windows
Multiply the record by a shape that tapers smoothly to zero at both ends. Then whatever the frequency, the periodic extension joins at zero and there is no discontinuity.
x_w[n]=x[n]\cdot w[n]
The cost, and it is unavoidable: tapering makes the main peak wider. You trade resolution for dynamic range, and every window is a different point on that trade.
| Window | Sidelobes | Main lobe width | Use for |
|---|---|---|---|
| Rectangular (none) | -13 dB | 1 bin | signals exactly on a bin |
| Hann | -31 dB | 2 bins | the general default |
| Hamming | -43 dB | 2 bins | slightly better first sidelobe |
| Blackman | -58 dB | 3 bins | finding small tones near big ones |
| Flat-top | -93 dB | 5 bins | measuring amplitude accurately |
| Kaiser | adjustable | adjustable | when you want to choose the trade |
The Hann window is w[n]=0.5\left(1-\cos\frac{2\pi n}{N-1}\right) — a raised cosine, zero at both ends. It is the sensible default and covers most work.
The flat-top window exists for a specific reason worth knowing. With other windows, a tone falling between bins reads low in amplitude — up to 1.4 dB low with a Hann window — because the peak lands in the trough between two bin responses. The flat-top window has a deliberately flattened main lobe so the amplitude reads correctly regardless of where the tone falls, at the cost of a very wide peak. Use it when you care about how big something is; use Blackman when you care about whether something small is there.
4. Zero-padding — what it does and does not do
Append zeros to make the record longer before transforming.
What it does: produces more output bins, so the peak is drawn with more points. The spectrum looks smoother, and the peak's location can be read more precisely by eye.
What it does not do: improve resolution. Two tones that were unresolved remain unresolved. No new information was added — the zeros contain nothing.
The distinction stated plainly: resolution is the ability to separate two nearby components, and it is set entirely by the record length in seconds. Zero-padding is interpolation of the existing spectrum, drawing the same curve with more points.
Where it is genuinely useful:
- Getting to a power of two so a radix-2 FFT can be used.
- Reading a single peak's frequency more precisely, when you already know there is only one peak.
- Avoiding the circular convolution wrap of Chapter 4.3.
5. The FFT
Direct computation of the DFT needs N multiply-accumulates per output bin, and there are N bins: N^2 operations.
The Cooley-Tukey idea: split the sum into even-indexed and odd-indexed samples.
X[k]=\sum_{n \text{ even}}x[n]W_N^{kn}+\sum_{n\text{ odd}}x[n]W_N^{kn}
Substitute n=2r in the first and n=2r+1 in the second:
=\sum_{r=0}^{N/2-1}x[2r]W_N^{2kr}+W_N^{k}\sum_{r=0}^{N/2-1}x[2r+1]W_N^{2kr}
Now the crucial observation: W_N^2 = e^{-j4\pi/N}=e^{-j2\pi/(N/2)}=W_{N/2}.
\boxed{X[k]=E[k]+W_N^kO[k]}
Each half-sum is itself a DFT of length N/2 — one of the even samples, one of the odd. So an N-point DFT is two N/2-point DFTs plus N multiplications to combine them.
And there is a second economy. Because W_N^{k+N/2}=-W_N^k:
X[k+N/2]=E[k]-W_N^kO[k]
So one multiplication produces two outputs. That pair of operations — one multiply, one add, one subtract — is the butterfly, named for the shape of its data-flow diagram.
The cost
Recursing all the way down, there are \log_2N stages, each doing N/2 butterflies:
\text{operations}=\frac N2\log_2N \quad\text{complex multiplies}
| N | Direct N^2 | FFT \frac N2\log_2N | Speedup |
|---|---|---|---|
| 64 | 4,096 | 192 | 21× |
| 1,024 | 1,048,576 | 5,120 | 205× |
| 65,536 | 4.3\times10^9 | 524,288 | 8,192× |
| 1,048,576 | 1.1\times10^{12} | 1.05\times10^7 | 104,858× |
The advantage grows without limit. A million-point transform that would take an hour takes 30 milliseconds.
Bit reversal
The recursive splitting into even and odd reorders the input. After \log_2 N splits, sample n ends up at the position whose binary representation is n's reversed.
For N=8: index 1 (001) goes to position 4 (100); index 3 (011) goes to 6 (110); indices 0 and 7 stay put.
This is why FFT implementations have a bit-reversal step, and why some accept the input in normal order and produce output in bit-reversed order, or the reverse. For applications that only multiply spectra together — fast convolution — you can skip the reordering entirely at both ends, since the two reversals cancel.
6. The FFT's history, which is worth telling
Cooley and Tukey published in 1965, at IBM and Princeton. Tukey had described the idea at a meeting of President Kennedy's Science Advisory Committee, where the problem under discussion was detecting Soviet nuclear tests by analysing seismometer records — a computation that was impractical at the time.
But Gauss had it in 1805. He used the method to interpolate the orbits of the asteroids Pallas and Juno, wrote it up in Latin, and never published; it appeared in his collected works in 1866 and went unnoticed for a century. Gauss's version predates Fourier's own 1807 paper.
Between them, at least a dozen people rediscovered pieces of it, including Danielson and Lanczos in 1942 for X-ray crystallography.
Why 1965 was the moment it mattered: computers had just become common enough that a fast algorithm changed what was possible rather than merely what was convenient. The FFT is routinely called one of the most important algorithms of the twentieth century, and the claim is defensible — digital audio, mobile phones, medical imaging, radio astronomy and JPEG all depend on it.
7. Fast convolution
The convolution theorem plus the FFT gives a genuinely important practical result.
Direct convolution of a length-N signal with a length-M filter costs NM operations.
Via the FFT: transform both, multiply, transform back. Cost is 3\times\frac L2\log_2L + L where L=N+M-1 rounded up to a power of two.
Worked comparison. N=4096 samples, M=512 taps.
Direct: 4096\times512 = 2.1 million multiplies.
FFT: L=8192. Three transforms at \frac{8192}{2}\log_2 8192 = 4096\times13 = 53{,}248 each, so 159,744, plus 8192 multiplies. Total about 168,000 — twelve times fewer.
The break-even is around 30 to 60 taps on typical hardware. Below that, direct convolution wins because the FFT's overhead dominates.
Remember the circular convolution trap from Chapter 4.3: pad to at least N+M-1 before transforming, or the result's beginning is contaminated by its end.
For continuous streams, overlap-add is the standard method: break the input into blocks, convolve each block with the filter using the FFT, and add the overlapping tails where consecutive blocks meet. Each block's output is M-1 samples longer than its input, and those extra samples overlap into the next block's region.
8. The variants worth knowing
Radix-4 and split-radix. Splitting into four parts rather than two, or mixing the two strategies, reduces the operation count by a further 20 to 30%. Split-radix has the lowest known operation count for power-of-two sizes.
Mixed radix and Bluestein's algorithm. For sizes that are not powers of two. Any N can be transformed in O(N\log N), though power-of-two sizes remain fastest. Modern libraries handle any size, but choosing a power of two, or a product of small primes, is still worth a factor of two or three.
Real-input FFT. Exploiting the conjugate symmetry of section 1 halves the work for real signals. Always use it when your data is real, which it usually is.
The DCT. The discrete cosine transform uses only cosines by treating the signal as if it were mirrored at the boundary rather than wrapped. That mirroring removes the discontinuity that causes leakage, so the DCT concentrates energy into fewer coefficients than the DFT for typical real-world signals.
That energy compaction is why JPEG, MP3, AAC and MPEG video all use the DCT rather than the DFT. Chapter 5.4 and Chapter 5.5 both build on it.
The Goertzel algorithm. When you want only a few specific frequencies rather than the whole spectrum, a small recursive filter computes each one in N operations. For K frequencies the cost is KN against the FFT's \frac N2\log_2N, so Goertzel wins whenever K\lt\frac12\log_2N.
The classic application is telephone touch-tone decoding. Only eight frequencies matter, so eight Goertzel filters are far cheaper than an FFT, and this is what every telephone exchange did for decades.
The FFT tells you what frequencies are present. Chapter 5.2 is about changing them — designing filters that keep what you want and remove what you do not.
Every formula above, built from scratch
None of the results in this chapter are worth memorising, because each one can be rebuilt in under a minute from something simpler. What follows is that rebuilding, one result at a time, so the formula and the reason for it sit on the same page as the explanation that needed them.
The DFT
X[k]=\sum_{n=0}^{N-1}x[n]W_N^{kn}, \qquad W_N=e^{-j2\pi/N}
x[n]=\frac1N\sum_{k=0}^{N-1}X[k]W_N^{-kn}
Bin frequency and resolution:
f_k=k\frac{f_s}{N}, \qquad \Delta f=\frac{f_s}{N}=\frac{1}{T_{record}}
Read the second form: resolution is the reciprocal of how long you recorded. Sampling faster does not improve it; recording longer is the only way.
Conjugate symmetry for real input:
X[N-k]=X^*[k]
so only bins 0 to N/2 carry information.
Amplitude scaling. A sinusoid of amplitude A exactly on bin k gives |X[k]|=AN/2, so
A = \frac{2|X[k]|}{N\cdot\overline{w}}
where \overline w is the window's mean value. Both divisions are needed and forgetting either is the standard error.
The FFT
The decimation-in-time split:
X[k]=E[k]+W_N^kO[k], \qquad X[k+N/2]=E[k]-W_N^kO[k]
Where the second line comes from: W_N^{k+N/2}=W_N^kW_N^{N/2}=W_N^ke^{-j\pi}=-W_N^k. One multiplication produces two outputs, and that is the butterfly.
Cost:
\frac N2\log_2N \text{ complex multiplies}, \qquad N\log_2N \text{ complex adds}
against N^2 for the direct computation.
Fast convolution cost, for length L = next power of two above N+M-1:
3\cdot\frac L2\log_2L + L
Break-even against direct convolution is around 30 to 60 taps.
Processing gain:
G=10\log_{10}\frac N2 \ \text{dB}
A tone's power stays in one bin while noise spreads across all of them, so a longer FFT lowers the apparent noise floor.
Welch averaging:
\text{variance reduction}=\sqrt{K} \text{ for } K \text{ segments}
\text{resolution}=\frac{Kf_s}{N} \text{ with } K \text{ segments of a fixed record}
Peak interpolation:
\delta=\frac12\cdot\frac{\alpha-\gamma}{\alpha-2\beta+\gamma}, \qquad f=(k+\delta)\frac{f_s}{N}
with \alpha,\beta,\gamma the levels in dB of the bins either side of and at the peak.
Goertzel is cheaper than the FFT when:
K \lt \tfrac12\log_2N
for K wanted frequencies.
What the next chapter fixes
The transform tells you what frequencies are present. It does not remove any of them. Chapter 5.2 does: the two families of digital filter, why one can be built with perfectly linear phase and the other cannot, and how to pick between them when both would meet the specification.