The Discrete Fourier Transform and FFT

Chapter 16 defined the discrete-time Fourier transform, and Chapters 23 to 25 built the z-transform on top of it. Both are complete answers to the question "what frequencies is this sequence made of". Neither can be evaluated by a computer, and the reason is worth stating precisely, because the fix is forced on us by the difficulty rather than chosen for convenience.

Two things stand in the way. The sum runs over infinitely many samples, and no machine has an infinite record. And \(\Omega\) is a continuous variable, so \(X(e^{j\Omega})\) is a function on an interval — uncountably many numbers, and no machine has memory for that either. A computable transform therefore has to be finite twice over: finitely many samples going in, finitely many numbers coming out.

Both restrictions are natural in practice. Any real measurement is a finite record: you switch on the instrument, you collect \(N\) samples, you switch it off. And any real display of a spectrum is a finite list: an analyser draws a few hundred or a few thousand points, not a continuum. So the question is not whether to make these two restrictions but what they cost, and the honest answer — worked out in Section 26-4 and returned to in Section 26-9 — is that they cost rather more than a beginner expects.

The transform that results is the discrete Fourier transform. It is not a fifth transform to be memorised alongside the Fourier series, the Fourier transform, the DTFT and the z-transform. It is the DTFT of a finite-length sequence, read off at finitely many frequencies, and every one of its peculiarities can be traced back to that single sentence.

Take a sequence that is zero outside \(0 \le n \le N-1\). Its DTFT sum then has only \(N\) non-zero terms, and the first difficulty disappears by itself.

This is a finite sum, but still a function of a continuous \(\Omega\). So sample it. The DTFT is periodic in \(\Omega\) with period \(2\pi\), so one period carries everything; divide that period into \(N\) equal parts and evaluate there.

Why exactly \(N\) points? Because the sequence is specified by \(N\) numbers, so \(N\) numbers of spectrum is precisely enough to specify it back. Fewer would lose information irrecoverably; more would be redundant, since the extra values could be reconstructed from the first \(N\). The count is not a compromise — it is the exact bookkeeping that makes the transform invertible.

The exponential appears so often that it is given a name. The twiddle factor is

With it the pair reads \(X[k] = \sum_n x[n]W_N^{kn}\) and \(x[n] = \tfrac1N\sum_k X[k]W_N^{-kn}\). Geometrically, \(W_N^{kn}\) walks around the unit circle in steps of \(2\pi/N\); the three identities say that after \(N\) steps you are home, that the walk repeats, and that half a lap is a sign change. The third identity is small, obvious and the entire reason the FFT of Section 26-7 exists.

Two coefficients can be written down without any arithmetic. Setting \(k=0\) kills the exponential, so \(X[0] = \sum_n x[n]\) is the sum of the samples — the DC value, and a free check on every DFT you compute. If \(N\) is even, setting \(k = N/2\) gives \(e^{-j\pi n} = (-1)^n\), so \(X[N/2] = \sum_n (-1)^n x[n]\) is the alternating sum, the coefficient at \(\Omega = \pi\), the highest frequency a sampled signal can carry.

The figure is worth carrying around. The grey curve is the true spectrum of the finite record; the gold stems are the \(N\) numbers the DFT returns. Nothing in the gold stems says how the curve behaves between them, and a great many practical mistakes come from forgetting that the DFT reports a sampled spectrum, not a continuous one.

The inversion formula was asserted above; it deserves a proof, and the proof is short. Everything rests on one orthogonality statement about the twiddle factors, which is itself nothing more than the sum of a geometric series.

To see it, put \(r = e^{j2\pi(n-m)/N}\) and read the sum as \(\sum_{k=0}^{N-1} r^{k}\). If \(n-m\) is a multiple of \(N\), then \(r = 1\) and every term is 1, so the sum is \(N\) and the normalised value is 1. Otherwise \(r \neq 1\) and the finite geometric series applies:

The numerator vanishes because \(n-m\) is an integer, and the denominator does not vanish because \(r \ne 1\). Now substitute the forward transform into the claimed inverse and exchange the two finite sums, which is always legal:

The bracket is the orthogonality sum, which is 1 when \(m = n\) and 0 for every other \(m\) in the range, so the outer sum collapses to the single term \(x[n]\). That is the whole proof, and it explains where the \(1/N\) comes from: the exponentials are orthogonal but not normalised, and each has squared length \(N\).

So far the derivation assumed the sequence really was zero outside \(0 \le n \le N-1\). Suppose it is not — suppose \(x[n]\) runs on forever and we sample its DTFT at \(N\) points anyway. What sequence does the inverse DFT hand back?

The calculation is the same one as in the previous section, but now the inner sum over \(m\) runs over all integers, and the orthogonality relation says the bracket is 1 whenever \(m\) and \(n\) differ by a multiple of \(N\) — not only when they are equal.

Chapter 20 established the dual of this statement: sampling a signal in time makes its spectrum periodic, and if the signal is not band-limited the spectral copies overlap and alias. Here the roles are exchanged. Sampling the spectrum in frequency makes the signal periodic, and if the signal is not time-limited to \(N\) samples the time copies overlap. The phenomenon is called time-domain aliasing, and the two statements are the same theorem read in two directions.

This is the fact that explains everything strange about the DFT. Because the inverse transform can only ever return a periodic sequence, the DFT necessarily treats its input as one period of a periodic signal. It does not know the difference between the record you gave it and an endless repetition of that record. Shifts therefore wrap round, convolutions wrap round, and the index \(n\) is really \(n \bmod N\) everywhere.

Write \((n)_N\) for \(n\) reduced modulo \(N\). A circular shift by \(m\) is the sequence \(x[(n-m)_N]\): samples pushed off the right-hand end of the record reappear at the left. Picturing the \(N\) samples around a ring, as in the figure, makes it obvious — a circular shift is a rotation of the ring, and nothing is ever lost or created.

The transform of a circular shift is easy, and the proof shows where the wrap-around earns its keep. Put \(p = (n-m)_N\); as \(n\) runs once through \(0,\dots,N-1\), so does \(p\). Then \(W_N^{kn} = W_N^{k(p+m)} = W_N^{km}W_N^{kp}\), where the step is legitimate because \(W_N^{kN}=1\) absorbs the modulo. Hence

Circular convolution follows the same pattern. Define, for two length-\(N\) sequences,

Substituting into the DFT and exchanging the sums gives the theorem in two steps. With \(p = (n-m)_N\) again,

The remaining properties are collected below. Linearity needs no comment. The symmetry property is the DFT version of the conjugate symmetry of Chapter 14: for a real signal the second half of the spectrum is the mirror image of the first, so of the \(N\) complex numbers only about \(N/2\) are independent — which is why a real-input FFT can be run at roughly half the cost.

Filtering is convolution, and convolution is expensive: an \(L\)-point signal through an \(M\)-tap filter costs about \(LM\) multiplications done directly. The circular convolution theorem promises to replace that with two forward transforms, \(N\) multiplications and one inverse transform, which — once the FFT of the next section is available — is dramatically cheaper. The obstacle is the word circular.

The fix comes straight from Section 26-4. Circular convolution is linear convolution wrapped modulo \(N\); if the linear result is short enough to fit inside \(N\) samples, there is nothing to wrap and the two coincide. Chapter 7 gave the length of a linear convolution: an \(L\)-point sequence convolved with an \(M\)-point sequence produces \(L+M-1\) points. So the condition writes itself.

That handles one block. A real filtering problem has a signal that is far longer than the filter, or has no end at all — audio arriving from a microphone, samples arriving from an ADC. Transforming the whole thing is impossible, and would in any case introduce an unacceptable delay. The signal is therefore cut into blocks, and the two standard ways of stitching the blocked results back together are worth knowing by name.

Overlap-add cuts the input into non-overlapping blocks of length \(L\). Each block is convolved with the \(M\)-tap filter using an \(N\)-point DFT with \(N \ge L+M-1\), producing a block of output of length \(L+M-1\). Because convolution is linear and the input blocks add up to the input, the output blocks must add up to the output — but consecutive output blocks overlap by \(M-1\) samples, and in that overlap region they are added. The name is a literal description.

Overlap-save arranges the bookkeeping the other way. Input blocks of length \(N\) are taken with an overlap of \(M-1\) samples between consecutive blocks. Each block is circularly convolved with the padded filter, and the first \(M-1\) output samples of each block — precisely the ones corrupted by wrap-around — are thrown away. Nothing is added; the surviving pieces are simply concatenated. It needs no output buffer, which is why it is the usual choice in hardware.

Count the arithmetic in the DFT as written. Each of the \(N\) values \(X[k]\) is a sum of \(N\) products, so the direct computation costs \(N^2\) complex multiplications and \(N(N-1)\) complex additions. The quadratic growth is the problem. Doubling the record length quadruples the work, and a 4096-point transform needs almost seventeen million complex multiplications — perfectly respectable for a modern processor once, and hopeless if it must be done a thousand times a second.

Cooley and Tukey published the escape in 1965, though Gauss had found the same trick in 1805 and left it unpublished in Latin. The idea is divide and conquer, and the lever is the third twiddle identity, \(W_N^{k+N/2} = -W_N^{k}\).

Let \(N\) be even and split the sum by the parity of the sample index — the even-numbered samples in one group, the odd-numbered in the other. Writing \(n = 2r\) for the even terms and \(n = 2r+1\) for the odd:

Now the crucial observation about the exponent. Since \(W_N^{2} = e^{-j4\pi/N} = e^{-j2\pi/(N/2)} = W_{N/2}\), both remaining sums are exactly \(N/2\)-point DFTs — of the even-indexed subsequence and of the odd-indexed subsequence respectively. Call them \(G[k]\) and \(H[k]\).

This looks like it only covers \(k = 0,\dots,N/2-1\), since \(G\) and \(H\) are \(N/2\)-point transforms and are only defined that far. But they are periodic with period \(N/2\), so \(G[k+N/2] = G[k]\) and \(H[k+N/2] = H[k]\); and the twiddle identity gives \(W_N^{k+N/2} = -W_N^{k}\). Substituting \(k+N/2\) therefore costs nothing new:

Count what has been achieved. Computing two \(N/2\)-point DFTs directly costs \(2(N/2)^2 = N^2/2\) multiplications, and the recombination costs \(N/2\) more — already roughly half the original \(N^2\). But the two half-length transforms can be split the same way, and their halves after that. If \(N\) is a power of two the splitting continues until the transforms are of length 1, and a 1-point DFT is the identity: \(X[0] = x[0]\), no arithmetic at all.

There are \(\log_2 N\) levels of splitting, each level containing \(N/2\) butterflies, each butterfly costing one complex multiplication and two complex additions. The total is \(\tfrac{N}{2}\log_2 N\) complex multiplications and \(N\log_2 N\) complex additions, and the transform has gone from quadratic to very nearly linear.

The recombination formula is drawn as a two-input, two-output crossing whose shape gives it its name. The upper input passes straight through to both outputs; the lower input is multiplied by the twiddle factor once and then added at the top and subtracted at the bottom.

Stacking \(N/2\) butterflies per stage and \(\log_2 N\) stages produces the flow graph below for \(N = 8\). Read it left to right: the first stage combines samples in pairs, the second combines the pair results, the third combines those. Three stages, four butterflies each, twelve butterflies in all — against the sixty-four multiplications a direct 8-point DFT would need.

The input ordering looks arbitrary and is not. The first split separates even from odd indices; the second splits each half by the next bit of the index; and so on. After \(\log_2 N\) splits, sample \(n\) has migrated to the position whose binary representation is that of \(n\) written backwards. For \(N = 8\): index 1 is \(001\), reversed \(100 = 4\), so \(x[1]\) sits in position 4 — exactly where the figure shows it.

Reading the "position" column downwards gives \(0,4,2,6,1,5,3,7\), the input order in the flow graph. The permutation is its own inverse, so a single in-place pass of pairwise swaps performs it, and the whole transform then proceeds in the same array with no extra storage. That in-place property matters as much as the speed on small hardware.

There is a mirror-image algorithm, decimation in frequency, which splits the output index instead of the input index. It takes the inputs in natural order and delivers the outputs bit-reversed, has the same operation count, and differs only in whether the twiddle multiplication happens before or after the add-subtract. Neither is better; libraries choose whichever fits their memory layout.

The payoff is best appreciated numerically.

Suppose a continuous signal is sampled at \(f_s\) hertz, \(N\) samples are collected, and their DFT is taken. Bin \(k\) sits at \(\Omega_k = 2\pi k/N\) radians per sample, which by Chapter 20's relation \(\Omega = 2\pi f/f_s\) corresponds to the analogue frequency

The bin spacing \(\Delta f\) is the reciprocal of the record duration \(T\), and that single relation governs everything that follows. Watch a signal for a tenth of a second and the spectrum arrives on a 10 Hz grid, whatever the sampling rate.

Now the difficulty. Truncating an infinite sinusoid to \(N\) samples is multiplication by a rectangular window, and Chapter 14 established that multiplication in time is convolution in frequency. The true spectrum of a pure tone is a single line; convolved with the transform of the rectangular window — a Dirichlet kernel, the discrete cousin of the sinc — it becomes a main lobe of finite width surrounded by a fence of side lobes. This smearing is spectral leakage, and it is present in every DFT of every measured signal.

There is one lucky case. If the tone's frequency lands exactly on a bin, then an integer number of its cycles fits inside the record, the periodic extension of Section 26-4 joins up seamlessly, and every side lobe falls precisely on a zero of the kernel. The DFT then shows a single clean spike. Move the tone half a bin off centre and the record no longer joins up: the discontinuity at the seam has broadband content, the side lobes no longer land on zeros, and energy appears across the whole spectrum.

The cure is to taper the record so it joins up smoothly, which means multiplying by a window that falls to zero at both ends instead of cutting abruptly. A gentler taper lowers the side lobes dramatically, at the cost of a wider main lobe — you trade the ability to see a small tone next to a large one against the ability to separate two tones of similar size. That trade-off is the whole subject of window design.

Two further practical points close the chapter. The first is resolution. Two tones can be told apart only if their main lobes do not merge, so with a rectangular window they must be separated by roughly \(1/T\) hertz, and with a Hann or Hamming window by roughly \(2/T\). Resolution is bought with observation time and with nothing else — not with a faster sampling rate, and not with a longer transform of the same data.

That last remark is the second point. Zero-padding a record to a longer transform length produces more bins, spaced more finely, and a visibly smoother plot. It adds no information whatever: the extra samples are zeros, the record duration \(T\) is unchanged, and the finer grid merely interpolates the same underlying DTFT. Zero-padding is useful — it locates a peak more precisely and removes the picket-fence effect, in which a peak falling between bins is under-read by up to 3.9 dB — but a claim that it improves resolution is simply false.

Problems 1 to 3 exercise the definition and its properties; 4 and 5 concern circular convolution; 6 and 7 the FFT structure; 8 the practicalities of measurement. Check every numerical transform with \(X[0]\) and with Parseval before moving on.