The Fast Fourier Transform

15.1 Computing the DFT Directly

15.1.1 The DFT and twiddle factors

Chapter 8 introduced the DFT as the discrete transform of a discrete (finite-length) sequence. Recall its definition, written compactly using the twiddle factor（旋轉因子） $W_N$:

$W_N$ is a unit-magnitude complex number sitting at angle $-2\pi/N$ on the unit circle. Because $W_N^{r}$ only depends on $r \bmod N$, the powers $W_N^r$ take on only $N$ distinct values — the $N$-th roots of unity. Two consequences of this structure drive everything in this chapter:

In matrix form the DFT is a dense $N\times N$ matrix–vector product,

but the matrix is anything but arbitrary: every entry is one of only $N$ values, and the entries are arranged in a highly regular pattern. Figure 15-1 visualizes the real and imaginary parts of the $N=32$ DFT matrix — the structure（規律性結構）is exactly the opportunity for efficiency that the FFT exploits.

15.1.2 Cost of direct evaluation

Evaluate the DFT sum literally. For each output point $X[k]$ we need $N$ complex multiplies and $N-1$ complex adds; there are $N$ output points ($k = 0,\dots,N-1$). Converting to real arithmetic with

$(a+jb)(c+jd) = (ac-bd) + j(ad+bc)$ = 4 real mults + 2 real adds,   1 complex add = 2 real adds,

the total damage is:

15.1.3 Goertzel's algorithm

Before attacking the full DFT, the lecture gives a first taste of exploiting twiddle structure: Goertzel's algorithm（Goertzel 演算法） computes a single DFT bin as the output of a small recursive filter. Start from the DFT sum and pull out a factor $W_N^{kN} = 1$:

The right-hand sum looks like a convolution evaluated at time $N$: it is $\sum_\ell x[\ell]\, h_k[N-\ell]$ with $h_k[n] = W_N^{-kn}u[n]$. Therefore

So one first-order complex recursion, clocked $N+1$ times (feeding $x_e[N]=0$ at the end), delivers $y_k[N] = X[k]$:

Properties claimed on the slide, worth memorizing:

• Separate "filters" for each $X[k]$ — you can compute just a few bins of interest (e.g. DTMF tone detection) without doing a whole DFT.
• No large buffer, no coefficient table — each filter needs one stored coefficient $W_N^{-k}$ and one state variable; samples are consumed as they arrive.
• Same complexity as the direct full DFT if you do compute all $N$ bins: $4N^2$ mults, $4N^2 - 2N$ adds.
• Halving the multiplies: make the denominator real by multiplying numerator and denominator by $(1 - W_N^{k}z^{-1})$:

The recursion now uses the real coefficient $2\cos(2\pi k/N)$ — only 2 real mults per step — and the complex numerator $(1 - W_N^k z^{-1})$ is evaluated only once, at the last step $n=N$.

15.2 Decimation-in-Time FFT

15.2.1 Splitting the DFT into even/odd halves

The decimation-in-time（時間抽取）DIT FFT rearranges the DFT sum into two halves: terms with even time index $n=2m$ and terms with odd index $n = 2m+1$ (assume $N$ even, ultimately $N = 2^M$):

The key step uses $W_N^{2mk} = W_{N/2}^{mk}$ (index halving). $X_0$ is the $N/2$-point DFT of the even-indexed samples, and $X_1$ is the $N/2$-point DFT of the odd-indexed samples. Because an $N/2$-point DFT is periodic in $k$ with period $N/2$, we evaluate both at $\langle k\rangle_{N/2}$ ($k$ modulo $N/2$). The result, in the lecture's compact notation:

Why does this save anything? Count again:

One split halves the computation (the $\varepsilon$ accounts for the $N$ extra twiddle multiplies and adds needed to combine the halves). We can evaluate an $N$-pt DFT as two $N/2$-pt DFTs plus a few mults/adds.

15.2.2 One-stage DIT flowgraph

Writing the combination for all $k = 0,\dots,N-1$ and using periodicity $X_0[\langle k\rangle_{N/2}]$, $X_1[\langle k\rangle_{N/2}]$, the structure for $N=8$ is shown in Fig. 15-3. Note carefully — as the slide stresses — the twiddle factors $W_N^k$ always apply to the odd-half output $X_1[\cdot]$, and they are NOT mirror-imaged: output $X[k+N/2]$ uses $W_N^{k+N/2} = -W_N^k$, not $W_N^{k}$ again.

15.2.3 Recursing to multiple stages

If decomposing one $\mathrm{DFT}_N$ into two smaller $\mathrm{DFT}_{N/2}$'s speeds things up — why not further divide into $\mathrm{DFT}_{N/4}$'s? Apply the identical split to $X_0[k]$ (the even-sample DFT): split its inputs into even-of-even and odd-of-even subsets:

where $X_{00}$ is the $N/4$-pt DFT of the even points within the even subset of $x[n]$ (i.e. $x[0],x[4],x[8],\dots$), $X_{01}$ the $N/4$-pt DFT of the odd points of the even subset, and so on. Note that the second-level twiddle is $W_{N/2}^{k}$ — for $N=8$ that means $W_4^k = W_8^{2k}$, so only $W_8^0$ and $W_8^2$ appear in the second stage.

Keep recursing until we hit 2-point DFTs, which need no multiplies at all:

An $N = 2^M$-point DFT thus reduces to $M = \log_2 N$ stages of twiddle factors and summation — the $O(N^2)$ part vanishes entirely.

15.2.4 The butterfly and its simplification

At every stage the elementary operation pairs two values $r$ and $r + N/2$ apart and produces two outputs. This is the butterfly（蝶形運算）. As first derived, it needs two complex multiplies, $W_N^r$ and $W_N^{r+N/2}$. But by half-period symmetry,

so the second multiply is just the negative of the first — i.e. subtract rather than add. Each butterfly then costs one complex multiply and two complex adds:

15.3 The Complete Radix-2 FFT

15.3.1 The 8-point DIT flowgraph

Assembling all $M = \log_2 8 = 3$ stages of simplified butterflies gives the classic flowgraph of Fig. 15-5. Trace any output back: stage 1 computes four 2-pt DFTs, stage 2 merges them into two 4-pt DFTs (twiddles $W_8^0, W_8^2$ — i.e. $W_4^0, W_4^1$), stage 3 merges those into the final 8-pt DFT (twiddles $W_8^0,\dots,W_8^3$).

15.3.2 Bit-reversed input ordering

Each DIT split sends even-indexed samples up and odd-indexed samples down — i.e. it sorts on the least significant bit (LSB) of the time index. Recursing sorts on the next bit, and so on. After $M$ levels, the sample at input row $r$ is $x[\text{bit-reverse}(r)]$: the address with its $M$-bit binary representation reversed（位元反轉）.

15.3.3 In-place computation

Look again at Fig. 15-5: every butterfly takes two values from rows $(p, q)$ and writes its two outputs back to the same two rows $(p, q)$. No other butterfly in that stage touches those rows. Consequently the whole FFT is an in-place algorithm（原位演算法）: one length-$N$ complex array suffices — after the initial bit-reversal shuffle, the $M$ stages overwrite the array sequentially, and the final array contents are $X[0],\dots,X[N-1]$ in natural order. Memory cost: $N$ complex words plus a small twiddle table (or none for Goertzel-style on-the-fly generation).

15.3.4 Complexity: (N/2) log₂ N

An $N = 2^M$-pt DFT reduces to $M = \log_2 N$ stages. Each stage has $N/2$ butterflies, each with at most 1 complex multiply (4 real mults) and 2 complex adds (4 real adds):

The quadratic term is gone. Figure 15-6 shows just how dramatic the difference is — at $N = 2^{16}$ the FFT needs about 3 orders of magnitude fewer multiplies.

15.4 FFT Variations

15.4.1 FFT for other values of N

Having $N = 2^M$ meant we could divide each stage into 2 halves — the radix-2 FFT（基數 2 FFT）. The same divide-and-conquer approach works more generally:

• $N = 3^M$ → radix-3 FFT (split into thirds);
• $N = 4^M$ → radix-4 FFT — a more optimized radix-2 (fewer twiddle multiplies per output);
• etc. for any prime radix;
• composite $N = a\cdot b\cdot c\cdot d$ → mixed radix（混合基數）: a different $N/r$-point split at each stage;
• … or just zero-pad（補零） the sequence up to the next $N = 2^M$ (recall from Ch. 8 that zero-padding only interpolates the spectrum more finely — it changes which frequencies are sampled, not the underlying DTFT).

15.4.2 Inverse FFT via the forward FFT

Recall the IDFT:

The only differences from the forward DFT are the $1/N$ scale and the sign of the exponent. Conjugating twice fixes the sign:

The bracketed sum is exactly a forward DFT of $x'[n] = X^*[k]\big|_{k=n}$ — a "time" sequence made from the spectrum. So the recipe is: conjugate the spectrum, run the same FFT engine, conjugate the result, divide by $N$. One FFT implementation serves both directions; in a pure real flowgraph this amounts to swapping the roles of real/imaginary parts with a couple of sign flips and $1/N$ scalings at the ports.

15.4.3 DFT of real sequences — two for the price of one

If $x[n]$ is pure-real, the DFT wastes multiplies: by conjugate symmetry（共軛對稱）, $X[k] = X^*[-k] = X^*[N-k]$, so half the outputs are redundant. The lecture's trick: pack two real sequences into one complex FFT. Given real $x[n]$ and $w[n]$, form $y[n] = x[n] + j\,w[n]$ and take a single $N$-pt DFT $V[k] = \mathrm{DFT}\{y[n]\} = X[k] + Y[k]$, where $Y[k] = \mathrm{DFT}\{j\,w[n]\} = jW[k]$. Then

(Here $V^*[-k] = V^*[\langle N-k\rangle_N]$, and the cancellations use the conjugate symmetry of the real-input spectra $X$ and $W$.) That is: compute the DFTs of two $N$-pt real sequences with a single $N$-pt complex DFT — Exercise 5 verifies this numerically to machine precision.

15.5 The FFT in MATLAB

MATLAB's fft/ifft call the FFTW library — mixed-radix $O(N\log N)$ for every length, fastest for $N = 2^M$:

N  = 4096;
x  = randn(1, N);
X  = fft(x);              % N-point FFT (FFTW under the hood)
xr = ifft(X);             % inverse; max(abs(xr - x)) ~ 1e-16
X2 = fft(x, 8192);        % zero-pads x to 8192 points first

% --- timing: power-of-two vs prime length ---
N1 = 2^16;  N2 = 65537;                 % 65536 vs the prime 65537
x1 = randn(1, N1);  x2 = randn(1, N2);
t1 = timeit(@() fft(x1));               % ~ a few ms
t2 = timeit(@() fft(x2));               % noticeably slower
fprintf('N = %6d : %8.3f ms\n', N1, 1e3*t1);
fprintf('N = %6d : %8.3f ms\n', N2, 1e3*t2);

Figure 15-7 shows measured run times on this machine (numpy's FFTPACK-style engine; MATLAB/FFTW behaves the same way qualitatively): both power-of-two and prime lengths scale like $N \log N$, but the prime-length constant is ~5× worse, and the direct $O(N^2)$ matrix product blows past both almost immediately.

15.6 Summary

Exercises

N = 64;  n = 0:N-1;
x = cos(2*pi*5*n/N) + 0.5*cos(2*pi*12*n/N);   % real sequence 1
w = sin(2*pi*9*n/N) + 0.3*randn(1,N);         % real sequence 2

y  = x + 1j*w;                                % pack into ONE complex sequence
V  = fft(y);                                  % single N-pt complex FFT
Vc = conj(V(mod(N-(0:N-1), N) + 1));          % V*[_N]  (= V*[-k])

X = (V + Vc)/2;                               % spectrum of x
W = (V - Vc)/(2j);                            % spectrum of w  (-j/2 form)

errX = max(abs(X - fft(x)))                   % ~ 3.6e-15
errW = max(abs(W - fft(w)))                   % ~ 5.3e-15

subplot(2,1,1); stem(n, abs(X)); ylabel('|X[k]|');
title('Recovered from one complex FFT');
subplot(2,1,2); stem(n, abs(W)); ylabel('|W[k]|'); xlabel('k');