Report an errorErrata
Part I · Module 02

The DFT and FFT, from correlation to butterflies

Every FFT bin is the correlation of your signal with one rotating arrow. This module builds the DFT from that single idea, computes it by hand for N = 4 and N = 8, states exactly what the real and imaginary part of each bin mean, and then derives the radix-2 FFT one butterfly at a time.

You will be able to

  • Compute \(X[k]\) by hand and read \(\mathrm{Re}\,X\), \(\mathrm{Im}\,X\), \(\lvert X\rvert\) and \(\angle X\) as cosine content, sine content, amplitude and phase
  • Predict bin values for a tone (\(A N/2\)), a constant (\(N c\)) and the Nyquist bin, and explain why the upper half mirrors
  • Explain leakage, scalloping (−3.92 dB) and zero padding with the drifting-probe picture
  • Derive the radix-2 FFT, draw the N = 8 butterfly with bit-reversed input, and count \(N\log_2 N\) vs \(N^2\)
  • Build correct frequency axes in MATLAB with fft, fftshift and odd or even N

Assumes

  • Module 01 · Complex numbers and phasors: phasors, \(e^{j\theta} = \cos\theta + j\sin\theta\), conjugate, multiplication = rotate and scale
  • Finite sums and the geometric series
  • Next: Module 03 uses the DFT for the convolution theorem, and Module 04 uses it as a measurement of sampled records
Exact formulaclosed form or exact finite sum, plotted as isSimulationnumerical run in this page (JS) or in the Octave scripts, named in a noteApprox. modelsimplified physics; assumptions statedSchematicdrawing for intuition, not to scale

02.1The inner product: how much of this shape is in my signal?

Take a signal \(x[n]\) and a probe shape \(p[n]\), both N samples long. Multiply sample by sample and add. The result is the inner product:

\[ \begin{aligned} \langle x, p\rangle &= \sum_{n=0}^{N-1} x[n]\,p[n] && (\text{real } p) \\ \langle x, p\rangle &= \sum_{n=0}^{N-1} x[n]\,p^{*}[n] && (\text{complex } p) \\ c_p &= \frac{\langle x, p\rangle}{\langle p, p\rangle} \end{aligned} \]
(02.1)
\(x[n]\)
signal samples, n = 0 … N − 1
\(p[n]\)
probe shape (template), same length
\(p^{*}\)
complex conjugate of the probe
\(c_p\)
projection coefficient: how many copies of p are in x
\(\langle x,p\rangle\)
inner product: multiply sample by sample and sum
\(N\)
record length, samples
\(n\)
sample index, 0 … N−1
\(\sum\)
sum over the index shown
Large and positive: x looks like p. Large and negative: x looks like −p. Near zero: p is not in x. Dividing by \(\langle p, p\rangle\) gives "how many copies of p", the projection coefficient.

You already use this every day. DAS beamforming is an inner product of the channel data with a delay-and-apodization pattern; a matched filter is an inner product with the pulse at each lag. A DFT bin is the same operation with a sampled sinusoid as the probe.

By hand: four samples, one cosine probe

Let \(x = [1, 2, 3, 4]\) and the probe be one cycle of cosine over four samples, \(c[n] = \cos(2\pi n/4) = [1, 0, -1, 0]\). Products: \([1\cdot 1,\ 2\cdot 0,\ 3\cdot(-1),\ 4\cdot 0]\) \(= [1, 0, -3, 0]\). Sum: \(\langle x, c\rangle = -2\). Since \(\langle c, c\rangle = 1 + 0 + 1 + 0 = 2\), x contains \(-2/2 =\) −1 times that cosine. With the sine probe \(s[n] = \sin(2\pi n/4) = [0, 1, 0, -1]\): products \([0, 2, 0, -4]\), sum \(\langle x, s\rangle = -2\), so x also contains −1 times the sine. The figure draws both.

Figure 02.1

Multiply sample by sample, then add

Exact formula
x[n] = [1, 2, 3, 4]probeproduct x[n]·p[n] and its sum
Read: the bars in the right column are the products; their sum is the inner product. Both probes give −2. Section 02.4 shows these two numbers are exactly \(\mathrm{Re}\,X[1] = -2\) and \(-\mathrm{Im}\,X[1] = -2\).
Key idea

An inner product answers one question about one shape. The DFT asks that question N times, once per probe frequency, with probes chosen so that the answers do not interfere (02.3).

02.2One probe is not enough: correlate with a cosine and a sine at once

Suppose \(x[n] = A\cos(2\pi k_0 n/N + \varphi)\) with unknown amplitude A and phase \(\varphi\), and \(0 < k_0 < N/2\) an integer. Correlate it with a cosine at the same frequency. Use \(\cos a\cos b\) \(= \tfrac12[\cos(a - b) + \cos(a + b)]\) with \(\theta = 2\pi k_0 n/N\):

\[ \begin{aligned} \sum_{n} A\cos(\theta + \varphi)\cos\theta &= \frac{A}{2}\sum_{n}\big[\cos\varphi + \cos(2\theta + \varphi)\big] \\ &= \frac{AN}{2}\cos\varphi + 0 \end{aligned} \]
(02.2)
\(A,\ \varphi\)
amplitude and phase of the tone
\(k_0\)
tone frequency in bins (cycles per record)
\(\theta\)
\(2\pi k_0 n/N\), the probe angle at sample n, rad
\(N\)
record length, samples
\(n\)
sample index
\(\sum_n\)
sum over \(n=0\ldots N-1\)
The second term is a sinusoid at frequency \(2k_0\) summed over whole periods, which is zero (02.3). So the cosine-only correlation is \((AN/2)\cos\varphi\): it depends on the unknown phase and is zero when \(\varphi = 90°\).

A cosine probe alone cannot tell "small amplitude" from "wrong phase". Add a sine probe. With \(\cos a\sin b\) \(= \tfrac12[\sin(a + b) - \sin(a - b)]\):

\[ \begin{aligned} \sum_{n} A\cos(\theta + \varphi)\sin\theta &= \frac{A}{2}\sum_{n}\big[\sin(2\theta + \varphi) - \sin\varphi\big] \\ &= -\frac{AN}{2}\sin\varphi \end{aligned} \]
(02.3)
\(A,\ \varphi,\ \theta\)
as in (02.2)
\(N\)
record length, samples
\(n\)
sample index
\(\sum_n\)
sum over \(n=0\ldots N-1\)

Now combine the two numbers as one complex number, cosine correlation minus j times sine correlation. Because \(e^{-j\theta} = \cos\theta - j\sin\theta\), that is one correlation with one complex probe:

\[ \begin{aligned} X[k] &= \sum_{n=0}^{N-1} x[n]\, e^{-j2\pi kn/N} \\ &= \sum_{n=0}^{N-1} x[n]\cos\frac{2\pi kn}{N} \;-\; j\sum_{n=0}^{N-1} x[n]\sin\frac{2\pi kn}{N} \end{aligned} \]
(02.4)
\(X[k]\)
DFT bin k, k = 0 … N − 1 (complex, same units as x times samples)
\(x[n]\)
signal samples, n = 0 … N − 1
\(N\)
record length in samples (FFT length)
\(k\)
bin index: k cycles per record, frequency \(k f_s/N\) for \(k < N/2\)
\(j\)
imaginary unit, \(j^2 = -1\)
For the tone at \(k = k_0\): \(X[k_0] = (AN/2)(\cos\varphi + j\sin\varphi)\) \(= (AN/2)\,e^{j\varphi}\). The magnitude \(AN/2\) no longer depends on \(\varphi\), and the angle is \(\varphi\) itself.
Part of X[k]EqualsPhysical meaning for real x
\(\mathrm{Re}\,X[k]\)\(\sum x[n]\cos(2\pi kn/N)\)How much cosine at frequency k: the in-phase (I) content, \((N/2)\) × cosine amplitude
\(\mathrm{Im}\,X[k]\)\(-\sum x[n]\sin(2\pi kn/N)\)Minus how much sine: the quadrature (Q) content, \(-(N/2)\) × sine amplitude
\(\lvert X[k]\rvert\)\(\sqrt{\mathrm{Re}^2 + \mathrm{Im}^2}\)\((N/2)\) × amplitude A of the component, whatever its phase
\(\angle X[k]\)\(\mathrm{atan2}(\mathrm{Im}, \mathrm{Re})\)Phase \(\varphi\) of the component written as \(A\cos(2\pi kn/N + \varphi)\)

The minus sign in the exponent is a convention chosen so that \(\angle X[k] = +\varphi\). With \(e^{+j}\) in the forward transform you would read \(-\varphi\). MATLAB, numpy and almost every text use \(e^{-j}\) forward.

Demo 02.1

The correlation wheel: signal × rotating probe, running sum

Exact formula
x[n]probe cos (solid), −sin (dashed)products x[n]e−j2πkn/N tip to tailX[k]
Try: with \(k = k_0 = 3\), sweep \(\varphi\): the cosine-only sum (\(\mathrm{Re}\,X\)) swings between ±8 and crosses 0 at \(\varphi = \pm 90°\), while \(\lvert X\rvert\) stays at \(AN/2 = 8\). Then set k = 4 with \(k_0 = 3\): the arrows curl into a closed polygon and X = 0. Set \(k_0 = 3.5\): the curl is only half closed (leakage, 02.7).
Spin the signal backwards

Multiplying by \(e^{-j2\pi kn/N}\) rotates sample n back by k turns per N samples. A component at frequency k stops rotating, so its arrows all point the same way and add up to \((AN/2)e^{j\varphi}\). Every other component keeps spinning and its arrows cancel. This is the I/Q demodulator of Module 08 (mix with \(e^{-j\omega_0 t}\), then low-pass), with the low-pass replaced by a plain sum over N samples.

02.3Why the probes do not interfere: orthogonality of the roots of unity

The DFT uses the N probes \(e^{j2\pi kn/N}\), k = 0 … N − 1. Correlating probe k with probe l gives a sum of the arrows \(e^{j2\pi(k-l)n/N}\). Let \(d = k - l\) and \(r = e^{j2\pi d/N}\), so the terms are \(r^0, r^1, \dots, r^{N-1}\). The geometric series gives

\[ \begin{gathered} \sum_{n=0}^{N-1} r^{n} = \frac{1 - r^{N}}{1 - r} = 0 \quad (r^N = e^{j2\pi d} = 1,\ r \ne 1) \\ \Longrightarrow\quad \sum_{n=0}^{N-1} e^{\,j2\pi(k-l)n/N} = N\,\delta[k-l] \end{gathered} \]
(02.5)
\(d = k - l\)
difference of the two probe indices, integer
\(r\)
\(e^{j2\pi d/N}\), one step of d/N turn on the unit circle
\(\delta[\cdot]\)
1 when the argument is a multiple of N, otherwise 0
\(N\)
record length, samples
\(n\)
sample index, 0 … N−1
\(k,\ l\)
two DFT bin indices (integers)
\(j\)
imaginary unit, \(j^2=-1\)
If d is a multiple of N, r = 1, every term is 1 and the sum is N. Otherwise \(r \ne 1\) and the sum is zero. So the N probes are mutually orthogonal, each with energy N.

Geometrically the terms are the N-th roots of unity, visited in steps of d. Placed tip to tail they walk around a closed polygon and return to the start, so the sum is zero. The cards draw it for N = 8.

Two consequences follow immediately. First, a component at bin l contributes nothing to any other bin \(k \ne l\), which is why a bin-centred tone lands in its own bin and nowhere else. Second, the transform can be undone. Multiply \(X[k]\) by \(e^{+j2\pi km/N}\), sum over k, and swap the order of summation:

\[ \begin{aligned} \sum_{k} X[k]\,e^{j2\pi km/N} &= \sum_{n} x[n] \sum_{k} e^{j2\pi k(m-n)/N} \\ &= \sum_{n} x[n]\cdot N\,\delta[m-n] = N\,x[m] \end{aligned} \]
(02.6)
\(X[k]\)
DFT of \(x\), as in (02.4)
\(x[n]\)
signal samples
\(m,\ n\)
output and input sample indices, 0 … N−1
\(k\)
bin index
\(N\)
record length, samples
\(\delta[\cdot]\)
Kronecker delta, as in (02.5)
\(j\)
imaginary unit, \(j^2=-1\)
So \(x[n] = \frac{1}{N}\sum_k X[k]\,e^{+j2\pi kn/N}\): the inverse DFT (02.13). The 1/N comes from the probe energy N.
Pitfall

Orthogonality needs whole periods inside the N-sample record. A tone at 3.5 bins is not one of the N probes, so it is not orthogonal to its neighbours and spreads into all bins (02.7). Sine and cosine at the same bin are also orthogonal, \(\sum\cos\theta\sin\theta = \tfrac12\sum\sin 2\theta = 0\), except at k = 0 and k = N/2, where the sine probe is identically zero.

02.4The DFT by hand: N = 4, then N = 8

Write \(W = e^{-j2\pi/N}\) (one step clockwise around the unit circle). Then \(X[k] = \sum x[n]\,W^{kn}\). For N = 4, \(W = e^{-j\pi/2} = -j\), so every probe value is one of \(1, -j, -1, j\). A 4-point DFT needs no multiplications at all, only sign changes and swapping real with imaginary.

\[ \begin{gathered} W^{0} = 1,\quad W^{1} = -j,\quad W^{2} = -1,\quad W^{3} = j \\ W^{4} = 1 = W^{0}\ \dots \qquad (\text{exponent } kn \bmod 4) \end{gathered} \]
(02.7)
\(W\)
\(e^{-j2\pi/N}\), the twiddle base (here N = 4)
\(W^{kn}\)
probe value for bin k at sample n
\(W\)
twiddle factor, \(e^{-j2\pi/N}\) with \(N=4\)
\(k,\ n\)
bin and sample indices
\(\bmod\)
remainder after division
\(j\)
imaginary unit, \(j^2=-1\)

Every multiply, x = [1, 2, 3, 4]

kn = 0: \(x\cdot W^{0}\)n = 1: \(x\cdot W^{k}\)n = 2: \(x\cdot W^{2k}\)n = 3: \(x\cdot W^{3k}\)X[k] = row sum
01·1 = 12·1 = 23·1 = 34·1 = 410
11·1 = 12·(−j) = −2j3·(−1) = −34·(j) = 4j−2 + 2j
21·1 = 12·(−1) = −23·(1) = 34·(−1) = −4−2
31·1 = 12·(j) = 2j3·(−1) = −34·(−j) = −4j−2 − 2j

In matrix form \(\mathbf{X} = \mathbf{F}_4\mathbf{x}\), where row k of \(\mathbf{F}_4\) is the probe for bin k:

\[ \begin{bmatrix} 10 \\ -2+2j \\ -2 \\ -2-2j \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} \]
(02.8)
\((\mathbf{F}_N)_{kn}\)
\(W^{kn}\), row k = probe for bin k, column n = sample n
\(\text{left vector}\)
DFT \(X[k]\), \(k=0\ldots3\)
\(\text{matrix}\)
DFT matrix \(\mathbf{F}_4\) of probe values \(W^{kn}\)
\(\text{right vector}\)
signal \(x[n]=(1,2,3,4)\)
\(j\)
imaginary unit, \(j^2=-1\)
\(\mathbf{F}_4\) is symmetric (\(F_{kn} = W^{kn} = F_{nk}\)) and \(\mathbf{F}_4^{H}\mathbf{F}_4 = 4\mathbf{I}\). Module 04 treats \(\mathbf{F}\) as a change of basis.
Figure 02.4

The 8 × 8 DFT matrix drawn as clock hands

Exact formulaEach hand points at angle −45°·(kn mod 8): the exact entry \(W^{kn}\) of \(\mathbf{F}_8\) on the unit circle.
Read: row k is the probe for bin k; moving right, its hand turns k eighths of a turn clockwise per sample. Row 1 (highlighted) makes one full turn; row 4 alternates ±1 (Nyquist); row 7 turns the other way, the conjugate of row 1. Column n read downward is the same pattern (the matrix is symmetric). The DFT is "multiply x by each row and add".

What each of the four numbers means

BinX[k]Re: cosine partIm: −(sine part)|X|, ∠XComponent of x[n] it encodes
0 (DC)10100 (always)10, 0°mean = \(X[0]/N\) = 2.5
1 (fs/4)−2 + 2j−2 → cos amplitude \(2\cdot(-2)/4 = -1\)+2 → sin amplitude \(-2\cdot 2/4 = -1\)2.83, 135°\(\tfrac{2}{4}\cdot 2.83\cos(\pi n/2 + 135°)\) = \(-\cos(\pi n/2) - \sin(\pi n/2)\)
2 (Nyquist)−2−20 (always)2, 180°\(X[2](-1)^n/N\) = \(-0.5(-1)^n\)
3 (−fs/4)−2 − 2j−2−22.83, −135°\(X[1]^{*}\): the negative-frequency half of bin 1, no new information

Check by adding the pieces: \(x[n] = 2.5 - \cos(\pi n/2)\) \(- \sin(\pi n/2) - 0.5(-1)^n\). At n = 0: 2.5 − 1 − 0 − 0.5 = 1. At n = 1: 2.5 − 0 − 1 + 0.5 = 2. At n = 2: 2.5 + 1 − 0 − 0.5 = 3. At n = 3: 2.5 − 0 + 1 + 0.5 = 4. The cosine and sine amplitudes −1 and −1 are exactly the two inner products of Figure 02.1 divided by \(\langle c, c\rangle = 2\).

N = 8: the probes need √2/2

For N = 8, \(W = e^{-j\pi/4} = \frac{\sqrt2}{2}(1 - j)\) \(\approx 0.7071 - 0.7071j\), and the probe values are the eight points \(1,\ 0.707 - 0.707j,\ -j,\) \(-0.707 - 0.707j,\ -1, \dots\) on the unit circle. Take the 4-sample pulse \(x = [1, 1, 1, 1, 0, 0, 0, 0]\):

\[ \begin{aligned} X[1] &= 1 + (0.707 - 0.707j) + (-j) \\ &\quad + (-0.707 - 0.707j) = 1 - 2.414j \\ \lvert X[1]\rvert &= 2.613, \qquad \angle X[1] = -67.5° \end{aligned} \]
(02.9)
\(X[1]\)
DFT bin 1 of the 8-sample example
\(\lvert X[1]\rvert,\ \angle X[1]\)
magnitude and phase (degrees) of that bin
\(j\)
imaginary unit, \(j^2=-1\)
The phase \(-67.5° = -2\pi\cdot 1\cdot 1.5/8\) says the pulse is centred at n = 1.5 (a delay is a linear phase, 02.9). \(X[2] = 1 - j - 1 + j = 0\): the pulse holds exactly one whole period of bin 2, so it cancels. Step through every bin in the demo.
Demo 02.2

The N = 8 DFT table, one row at a time

Exact formula
probe point e−j2πkn/8running sumRe XIm X
Try: "sin(2πn/8)" at k = 1: the cosine column sums to 0 and the −sin column to −4, so \(X[1] = -4j\) (pure imaginary means pure sine). At k = 7 you get \(+4j\), the conjugate. "impulse at n = 1": every bin has \(\lvert X\rvert = 1\) and the phase falls by 45° per bin, a pure delay.

02.5Reading a spectrum: real, imaginary, magnitude, phase, and the mirrored half

For real x and even N, the inverse DFT regroups into real cosines and sines. Bins k and N − k carry the same information (shown next), so they combine into one real sinusoid:

\[ \begin{aligned} x[n] = {} &\frac{X[0]}{N} + \frac{X[N/2]\,(-1)^{n}}{N} \\ &+ \frac{2}{N}\sum_{k=1}^{N/2-1}\mathrm{Re}\,X[k]\cos\frac{2\pi kn}{N} \\ &- \frac{2}{N}\sum_{k=1}^{N/2-1}\mathrm{Im}\,X[k]\sin\frac{2\pi kn}{N} \end{aligned} \]
(02.10)
\(X[0]/N\)
mean (DC) of x
\(\tfrac{2}{N}\mathrm{Re}\,X[k]\)
cosine amplitude at frequency \(k f_s/N\)
\(-\tfrac{2}{N}\mathrm{Im}\,X[k]\)
sine amplitude at frequency \(k f_s/N\)
\(X[N/2]/N\)
amplitude of the Nyquist alternation \((-1)^n\) (even N only)
\(x[n]\)
reconstructed real signal
\(n\)
sample index
\(N\)
record length (even), samples
\(X[k]\)
DFT of \(x\), as in (02.4)
\(k\)
bin index
\(\mathrm{Re},\ \mathrm{Im}\)
real and imaginary parts
\((-1)^n\)
alternating sign: the Nyquist cosine
Equivalently \(\frac{2}{N}\lvert X[k]\rvert\cos(2\pi kn/N + \angle X[k])\) for each \(1 \le k \le N/2 - 1\). This is the precise statement of what Re and Im mean.

Why the upper half mirrors: \(X[N-k] = X[k]^{*}\)

\[ \begin{aligned} X[N-k] &= \sum x[n]\,e^{-j2\pi(N-k)n/N} \\ &= \sum x[n]\,e^{-j2\pi n}\,e^{+j2\pi kn/N} = \sum x[n]\,e^{+j2\pi kn/N} \\ &= \Big(\sum x[n]\,e^{-j2\pi kn/N}\Big)^{*} = X[k]^{*} \end{aligned} \]
(02.11)
\(X[k]\)
DFT of a real record \(x[n]\)
\(N\)
record length, samples
\(n\)
sample index, 0 … N−1
\(k\)
bin index
\((\cdot)^{*}\)
complex conjugate
\(\sum\)
sum over \(n=0\ldots N-1\)
\(j\)
imaginary unit, \(j^2=-1\)
The last step needs x real (\(x^{*} = x\)). So Re X is even about N/2, Im X is odd, \(\lvert X\rvert\) is even, \(\angle X\) is odd. Bin N − k is frequency −k: \(e^{-j2\pi(N-k)n/N}\) and \(e^{+j2\pi kn/N}\) are the same samples, the aliasing of Module 04 inside the DFT itself.

Physically: a real cosine is the sum of two counter-rotating arrows, \(\cos\theta = (e^{j\theta} + e^{-j\theta})/2\) (Module 01). The DFT reports each arrow in its own bin, k for the positive one and N − k for the negative one, each with half the amplitude. Their imaginary parts are opposite and cancel in x[n].

DC and Nyquist are real

The probe for k = 0 is all ones and the probe for k = N/2 is \((-1)^n\); both are real, so \(X[0]\) and \(X[N/2]\) are real for real x. The sine probe at Nyquist is \(\sin(\pi n) = 0\) at every sample, so the sine part of a Nyquist tone is invisible: \(A\cos(\pi n + \varphi) = A\cos\varphi\cdot(-1)^n\). Its phase cannot be recovered, only \(A\cos\varphi\).

Scaling conventions

ConventionForward factorInverse factor|X[k]| of A cos, bin-centredX[0] of constant cUsed by
Unnormalized1\(1/N\)\(AN/2\)\(Nc\)MATLAB fft/ifft, numpy, FFTW
Unitary\(1/\sqrt N\)\(1/\sqrt N\)\(A\sqrt N/2\)\(\sqrt N\,c\)Linear algebra (\(\mathbf{F}/\sqrt N\) is unitary, Parseval without factors)
Amplitude\(1/N\)1\(A/2\)\(c\)Fourier-series coefficients; one-sided amplitude = \(2\lvert X[k]\rvert/N\)

Noise behaves differently from a tone. White noise and a bin-centred tone scale differently with N (unnormalized DFT):

\[ \begin{gathered} \mathrm{E}\,\lvert X_{\text{noise}}[k]\rvert^{2} = N\sigma^{2}, \qquad \lvert X_{\text{tone}}[k_0]\rvert^{2} = \frac{A^{2}N^{2}}{4} \\ \Longrightarrow\quad \mathrm{SNR}_{\text{bin}} = \frac{N}{2}\cdot\frac{A^{2}}{2\sigma^{2}} \end{gathered} \]
(02.12)
\(\sigma^2\)
per-sample variance of real white noise
\(A^2/(2\sigma^2)\)
per-sample SNR of the tone (signal power over noise power)
\(N/2\)
coherent processing gain of the DFT
\(\mathrm{E}\)
expected value (ensemble average)
\(X_{\text{noise}}[k]\)
DFT bin of white noise
\(X_{\text{tone}}[k_0]\)
DFT bin at the tone
\(k_0\)
tone frequency in bins
\(N\)
record length, samples
\(A\)
tone amplitude
\(\mathrm{SNR}_{\text{bin}}\)
power SNR in the tone bin (ratio)
The tone grows as \(N^2\), the noise as N: a coherent processing gain of N/2, \(10\log_{10}512 = 27.1\) dB for N = 1024. That is why a weak tone can sit far below the per-sample noise and still stand out in an FFT.
Demo 02.3

Spectrum explorer: Re, Im, |X| and ∠X of a tone you build

Exact formula
top: x[n]; Re X or |X|Im X or ∠Xfaded stems: mirror half (negative frequencies)
Try: set \(\varphi = 0\) (pure cosine): only Re is non-zero. Set \(\varphi = -90°\) (pure sine, \(\cos(\theta - 90°) = \sin\theta\)): Re vanishes and \(\mathrm{Im}\,X[k_0] = -AN/2\). Move \(k_0\) to 8 (Nyquist): the bin is real and equals \(NA\cos\varphi\); at \(\varphi = 90°\) the tone disappears entirely.
Pitfall: phase of an empty bin

\(\angle X\) of a bin that should be zero is atan2 of rounding noise (≈ \(10^{-15}\)) and is a random number between −180° and 180°. MATLAB happily plots it. Mask phase wherever \(\lvert X\rvert\) is below a threshold relative to the peak, as the demo does.

Ultrasound link

An RF line of N = 1024 samples at \(f_s\) = 40 MHz has bins \(\Delta f = f_s/N\) = 39.06 kHz apart. A 5 MHz probe with a 3–7 MHz band occupies bins 77–179 and their mirror 845–947; the rest is noise. Complex I/Q data break the mirror symmetry: after demodulation (Module 08) the bins above N/2 are genuine negative frequencies, which is how Doppler separates flow toward and away from the probe.

02.6The inverse DFT: rebuilding x[n] from N rotating arrows

The inverse transform proved in 02.3 reads as a recipe for drawing the signal:

\[ x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k]\, e^{+j2\pi kn/N}, \qquad \mathbf{x} = \frac{1}{N}\mathbf{F}^{H}\mathbf{X} \]
(02.13)
\(X[k]/N\)
complex amplitude of arrow k: length \(\lvert X[k]\rvert/N\), start angle \(\angle X[k]\)
\(e^{+j2\pi kn/N}\)
arrow k turning k revolutions every N samples (counter-clockwise)
\(\mathbf{F}^{H}\)
conjugate transpose of the DFT matrix \((\mathbf{F})_{kn} = W^{kn}\)
\(x[n]\)
reconstructed signal
\(n,\ k\)
sample and bin indices
\(N\)
record length, samples
\(\mathbf{x},\ \mathbf{X}\)
column vectors of \(x[n]\) and \(X[k]\)
\(\mathbf{F}\)
DFT matrix, \(F_{kn}=e^{-j2\pi kn/N}\)
\(j\)
imaginary unit, \(j^2=-1\)
Add the N arrows tip to tail; at sample n the tip of the chain is x[n].

In matrix form, \(\mathbf{X} = \mathbf{F}\mathbf{x}\) with \(F_{kn} = W^{kn}\), and orthogonality says \(\mathbf{F}^{H}\mathbf{F} = N\mathbf{I}\). Therefore \(\mathbf{x} = \frac{1}{N}\mathbf{F}^{H}\mathbf{X}\): the inverse is the conjugate transpose divided by N, and nothing has to be inverted. The energy relation follows from the same identity (Parseval):

\[ \sum_{n=0}^{N-1} \lvert x[n]\rvert^{2} = \frac{1}{N}\sum_{k=0}^{N-1} \lvert X[k]\rvert^{2} \]
(02.14)
\(x[n]\)
signal samples
\(X[k]\)
its DFT, as in (02.4)
\(N\)
record length, samples
\(n,\ k\)
sample and bin indices
\(\lvert\cdot\rvert^2\)
power (squared magnitude)
Example, x = [1, 2, 3, 4]: 1 + 4 + 9 + 16 = 30 and \((\lvert 10\rvert^2 + \lvert -2+2j\rvert^2 + \lvert -2\rvert^2 + \lvert -2-2j\rvert^2)/4\) \(= (100 + 8 + 4 + 8)/4 = 30\).
Demo 02.4

Inverse DFT: add the arrows, trace the signal

Exact formula
arrows X[k]ej2πkt/N/NRe(tip) vs tIm(tip) vs tsamples x[n]
Try: compare the two arrow sets. Both pass through every x[n] at integer t, because \(e^{j2\pi(k-N)n/N} = e^{j2\pi kn/N}\) at integer n. Between samples the "as stored" chain leaves the real axis (all arrows turn the same way), while the signed chain pairs each arrow with its counter-rotating conjugate and stays exactly real: that curve is the trigonometric interpolant of the samples. Reduce "arrows used" to see a partial Fourier sum.
What the imaginary part of the inverse means

If the rebuilt x[n] has an imaginary part larger than rounding error, the spectrum you fed in was not conjugate-symmetric: you edited a bin k without editing N − k. For a real signal every edit must be mirrored (or use ifft(X,'symmetric') in MATLAB, which assumes the mirror).

02.7Leakage: when the tone is not on a bin

Let the tone be \(e^{j2\pi k_0 n/N}\) with \(k_0\) not an integer, and write \(\delta = k_0 - k\) for its distance from bin k. The product with the probe is \(e^{j2\pi\delta n/N}\): the arrows no longer point one way, they turn by \(2\pi\delta/N\) per sample. Summing the geometric series:

\[ \begin{aligned} X[k] &= \sum_{n=0}^{N-1} e^{j2\pi\delta n/N} = e^{j\pi\delta(N-1)/N}\,\frac{\sin(\pi\delta)}{\sin(\pi\delta/N)} \\ &= e^{j\pi\delta(N-1)/N}\,D_N(\delta) \end{aligned} \]
(02.15)
\(\delta\)
\(k_0 - k\), offset of the tone from bin k, in bins
\(D_N(\delta)\)
Dirichlet kernel (periodic sinc), \(D_N(0) = N\)
\(e^{j\pi\delta(N-1)/N}\)
linear phase: the record is centred at n = (N − 1)/2
\(X[k]\)
DFT bin k of a unit-amplitude tone
\(k_0\)
tone frequency, bins
\(k\)
bin index
\(N\)
record length, samples
\(n\)
sample index
\(j\)
imaginary unit, \(j^2=-1\)
\(D_N = 0\) at every non-zero integer \(\delta\) (the bins of an on-bin tone); \(D_N(\delta) \approx N\,\mathrm{sinc}(\delta)\) for \(\lvert\delta\rvert \ll N\). Every bin now receives some of the tone: that is leakage.

The drifting arrows give the worst case without any algebra. At \(\delta = \tfrac12\) the product turns half a revolution over the record, so the N unit arrows bend into a half circle of arc length N. The chord of a half circle is its diameter, \((2/\pi)\) × arc, so \(\lvert X\rvert \approx (2/\pi)N\):

\[ \begin{gathered} \frac{\lvert X[k]\rvert}{N}\bigg|_{\delta = 1/2} = \frac{1}{N\sin\!\big(\pi/(2N)\big)} \;\xrightarrow{\;N \to \infty\;}\; \frac{2}{\pi} \\ 20\log_{10}\frac{2}{\pi} = -3.92\ \mathrm{dB} \end{gathered} \]
(02.16)
\(X[k]\)
DFT bin, as in (02.4)
\(N\)
record length, samples
\(\delta\)
offset of the tone from the bin center, bins (as in 02.15)
\(\log_{10}\)
base-10 logarithm
\(\mathrm{dB}\)
decibel, amplitude ratio \(20\log_{10}\)
The scalloping loss of the rectangular window (N = 16 already gives −3.91 dB). A Hann window \(w[n] = \tfrac12 - \tfrac12\cos(2\pi n/N)\) reduces it to 1.42 dB and pushes the first sidelobe from −13.3 dB to −31.5 dB, at the cost of a main lobe twice as wide (null to null 4 bins instead of 2) and an equivalent noise bandwidth of 1.5 bins. Module 04 covers window choice in depth.

Zero padding samples the same curve more finely

Append \((P - 1)N\) zeros and take a PN-point DFT. The zeros add nothing to the sum, so

\[ \begin{aligned} X_P[m] &= \sum_{n=0}^{N-1} x[n]\, e^{-j2\pi mn/(PN)} = X(e^{j\omega})\Big|_{\omega = 2\pi m/(PN)} \\ X(e^{j\omega}) &= \sum_{n} x[n]\,e^{-j\omega n} \end{aligned} \]
(02.17)
\(P\)
zero-padding factor (transform length PN)
\(X(e^{j\omega})\)
DTFT of the N-sample record, continuous in \(\omega\) (rad/sample)
\(m\)
bin index of the padded DFT, m = 0 … PN − 1
\(X_P[m]\)
DFT of the zero-padded record, \(PN\) points
\(x[n]\)
signal samples, \(n=0\ldots N-1\)
\(N\)
record length, samples
\(n\)
sample index
\(\omega\)
normalized angular frequency, rad/sample
\(j\)
imaginary unit, \(j^2=-1\)
The same DTFT evaluated on a grid P times finer. The width of the Dirichlet main lobe is still set by N, i.e. by the record duration \(T = N/f_s\). Zero padding interpolates; it cannot make two tones closer than about 1/T reliably separable.
Demo 02.5

Leakage explorer: off-bin tones, windows, zero padding

Exact formula
windowed x[n] and its periodic copyDTFT (continuous)N-point DFT binszero-padded DFT
Try: slide \(k_0\) from 6.0 to 7.0. At integer \(k_0\) the N bins sit on the peak and on the zeros of the kernel, so leakage is invisible; at 6.5 the two nearest bins are each about 3.9 dB down (3.92 dB for a complex tone in the large-N limit; the negative-frequency image of this real cosine moves them by a few tenths of a dB) and every other bin is lit. The faint copy in the top panel is the periodic extension the DFT implicitly assumes: an off-bin tone has a jump at the wrap. Switch to Hann and the jump disappears.
Ultrasound link

Spectral QUS estimates (attenuation by spectral shift, backscatter coefficient) use short gates. A 2 µs gate at \(f_s\) = 40 MHz is N = 80 samples, so \(\Delta f = f_s/N\) = 500 kHz. Zero padding to 256 gives a 156 kHz grid, but the resolution is still ≈ 1/T = 500 kHz, and leakage from the strong in-band energy biases the centroid unless the gate is windowed (Module 27).

02.8The FFT: split even and odd, reuse, repeat

The direct DFT costs N complex multiplies per bin, \(N^2\) in total. The radix-2 decimation-in-time FFT (Cooley and Tukey, 1965) gets the same numbers with \((N/2)\log_2 N\) multiplies. It takes four steps, each one identity.

Step 1. Split the sum into even and odd samples (n = 2m and n = 2m + 1, m = 0 … N/2 − 1):

\[ \begin{aligned} X[k] &= \sum_{m} x[2m]\,W_N^{2mk} + \sum_{m} x[2m+1]\,W_N^{(2m+1)k} \\ &= \sum_{m} x[2m]\,W_N^{2mk} + W_N^{k}\sum_{m} x[2m+1]\,W_N^{2mk} \end{aligned} \]
(02.18)
\(W_N\)
\(e^{-j2\pi/N}\), twiddle base for length N
\(m\)
index of the half-length sequences, m = 0 … N/2 − 1
\(X[k]\)
N-point DFT
\(x[n]\)
signal samples
\(N\)
record length (even), samples
\(k\)
bin index
\(m\)
pair index, \(0\ldots N/2-1\) (even sample \(2m\), odd \(2m+1\))
\(W_N^{p}\)
twiddle \(e^{-j2\pi p/N}\)

Step 2. Two steps on an N-point circle are one step on an N/2-point circle:

\[ \begin{gathered} W_N^{2mk} = e^{-j2\pi\cdot 2mk/N} = e^{-j2\pi mk/(N/2)} = W_{N/2}^{mk} \\ \Longrightarrow\quad X[k] = E[k] + W_N^{k}\,O[k] \end{gathered} \]
(02.19)
\(E[k]\)
(N/2)-point DFT of the even samples x[0], x[2], …
\(O[k]\)
(N/2)-point DFT of the odd samples x[1], x[3], …
\(W_N^{k}\)
twiddle factor
\(X[k]\)
N-point DFT
\(N\)
record length (even), samples
\(k\)
bin index
\(m\)
pair index
\(W_{N/2}\)
twiddle of the half-length DFT
\(j\)
imaginary unit, \(j^2=-1\)
The twiddle \(W_N^{k}\) rotates \(O[k]\) by \(-2\pi k/N\): the phase of the one-sample delay (on the original grid) of the odd samples relative to the even ones.

Step 3. Reuse for the upper half. E and O have period N/2, so \(E[k + N/2] = E[k]\) and \(O[k + N/2] = O[k]\). The twiddle factor goes half a turn further, \(W_N^{k+N/2} = W_N^{k}e^{-j\pi} = -W_N^{k}\). So one product \(W^{k}O[k]\) serves two outputs:

\[ \begin{aligned} X[k] &= E[k] + W_N^{k}\,O[k] \\ X[k + N/2] &= E[k] - W_N^{k}\,O[k], \qquad k = 0 \dots N/2 - 1 \end{aligned} \]
(02.20)
\(E[k],\ O[k]\)
half-length DFTs of even and odd samples
\(W_N^{k}\)
\(e^{-j2\pi k/N}\), twiddle factor
\(X[k],\ X[k+N/2]\)
the two outputs that share one complex multiply
\(N\)
record length (even), samples
\(X[k]\)
N-point DFT output
\(N\)
record length (even), samples
\(k\)
bin index, 0 … N/2−1
This pair is the butterfly: one complex multiply, one add, one subtract.

Step 4. Recurse. Each N/2-point DFT splits the same way, down to 1-point DFTs (a single sample is its own DFT). There are \(\log_2 N\) levels with N/2 butterflies each, so the cost obeys

\[ \begin{gathered} C(N) = 2\,C(N/2) + \frac{N}{2}, \qquad C(1) = 0 \\ \Longrightarrow\quad C(N) = \frac{N}{2}\log_2 N \ \text{complex multiplies} \\ N\log_2 N \ \text{complex additions} \end{gathered} \]
(02.21)
\(C(N)\)
number of complex multiplies for an N-point FFT
\(N\)
record length (power of 2), samples
\(\log_2\)
base-2 logarithm
\(\Longrightarrow\)
solving the recursion gives
Counting the trivial \(W^0 = 1\). Compare \(N^2\) multiplies for the direct sum.

By hand: the N = 4 example again, now as an FFT

Even samples [1, 3]: E = [1 + 3, 1 − 3] = [4, −2]. Odd samples [2, 4]: O = [6, −2]. Twiddles \(W_4^0 = 1\), \(W_4^1 = -j\), so \(W^1 O[1] = (-j)(-2) = 2j\). Butterflies: X[0] = 4 + 6 = 10, X[2] = 4 − 6 = −2, X[1] = −2 + 2j, X[3] = −2 − 2j. Same answer as the table in 02.4, with one non-trivial multiply instead of sixteen.

Figure 02.8

N = 4 FFT of x = [1, 2, 3, 4], every value shown

Exact formulaDataflow layout is a drawing; every number on it is computed by the page's radix-2 arithmetic (exact finite sums).
Read: left to right. Each X-shaped crossing is a butterfly: the top output is a + Wb, the bottom a − Wb. The input is in bit-reversed order [x0, x2, x1, x3].

Why the input ends up bit-reversed

The first split sorts samples by the last bit of n (even: 0, odd: 1). The second split sorts each half by the next bit, and so on. After \(\log_2 N\) splits, position p holds the sample whose index is p with its bits reversed. For N = 8: 0 (000), 4 (100), 2 (010), 6 (110), 1 (001), 5 (101), 3 (011), 7 (111). An in-place FFT permutes the input once, then runs the butterflies with outputs in natural order.

Figure 02.8b

Three even/odd splits produce the bit-reversed order

Schematic
Read: top to bottom. Blue indices are even, pink odd. Each level splits every box by one more bit of n, lowest bit first, so the final left-to-right order is n with its three bits reversed: position 1 (001) holds sample 4 (100), position 3 (011) holds sample 6 (110).
Demo 02.6

The N = 8 butterfly with live values

Exact formulaSignal-flow layout is a drawing; every value on it is computed exactly from (02.20).
a + Wb (top output)a − Wb (bottom output)selected butterfly
Try: step "Stages computed" from 0 to 3. Stage 1 is four 2-point DFTs (W = 1), stage 2 two 4-point DFTs (W = 1, −j), stage 3 one 8-point combination (\(W_8^{0\dots 3}\)). Inspect stage 3, butterfly 1: it combines E[1] and O[1] into X[1] and X[5]. Every output matches the direct DFT of Demo 02.2.
Figure 02.8c

The textbook N = 8 decimation-in-time flow graph

Schematic
Read: every node adds its two incoming branches. A twiddle label multiplies the branch leaving the lower node of a butterfly; −1 marks the branch into the lower output. Stage s uses \(W_{2^s}^{j} = W_8^{\,j\cdot 8/2^s}\): only \(W_8^0 = 1\) in stage 1, \(W_8^0, W_8^2 = -j\) in stage 2, \(W_8^{0\dots 3}\) in stage 3. Blue branches feed top outputs (a + Wb), pink branches bottom outputs (a − Wb), as in Demo 02.6.
NDirect DFT: \(N^2\) multipliesRadix-2: \((N/2)\log_2 N\)Ratio \(2N/\log_2 N\)Ultrasound example
864125.3this page
25665 5361 02464short spectral gate, zero padded
1 0241 048 5765 120205one RF line, 25.6 µs at 40 MHz
4 09616 777 21624 576683deep RF line at 62.5 MHz
65 5364.29 × 10⁹524 2888 192same count as a 256 × 256 fft2
Other FFTs, same idea

Decimation in frequency splits the outputs into even and odd bins instead: input in natural order, output bit-reversed, twiddle applied after the subtraction. Radix-4 and split-radix cut the multiply count further. Lengths that are not powers of two use mixed radices or Bluestein's chirp-z method; MATLAB's fft is built on FFTW, which picks a plan for any N, so padding to a power of two is a speed habit, not a requirement.

02.9FFTs in practice: index order, fftshift, frequency axes, two-for-one, 2D

MATLAB indexes from 1, so X(1) is bin 0 (DC) and X(k+1) is bin k. Bin k is frequency \(k f_s/N\) for \(k < N/2\) and \((k - N) f_s/N\) for \(k > N/2\); for even N, bin N/2 is \(\pm f_s/2\) (both at once, by aliasing).

TaskMATLAB (base)Note
Frequency axis, fft orderf = (0:N-1)*fs/N;upper half are negative frequencies
Centred axisXc = fftshift(X); fc = (-floor(N/2):ceil(N/2)-1)*fs/N;N even: −N/2 … N/2−1; N odd: −(N−1)/2 … (N−1)/2
Undo the shiftX = ifftshift(Xc);differs from fftshift for odd N
One-sided amplitude (N even)A = abs(X(1:N/2+1))/N; A(2:end-1) = 2*A(2:end-1);DC and Nyquist are not doubled
Delay by m samplesX.*exp(-2j*pi*(0:N-1)'*m/N)circular shift; for fractional m use the signed bins and keep the Nyquist bin real
2D DFTfft2(x) == fft(fft(x).').'rows, then columns (separable)

Two real FFTs for the price of one complex FFT

Pack two real signals as \(z = x + jy\) and take one FFT, \(Z = X + jY\). Conjugate symmetry of X and Y separates them again (indices mod N):

\[ \begin{gathered} X[k] = \frac{Z[k] + Z^{*}[N-k]}{2} \\ Y[k] = \frac{Z[k] - Z^{*}[N-k]}{2j} \end{gathered} \]
(02.22)
\(x,\ y\)
two real signals of length N (e.g. two RF lines)
\(Z[k]\)
N-point DFT of \(z = x + jy\)
\(X[k],\ Y[k]\)
N-point DFTs of the real signals \(x\) and \(y\)
\(Z[k]\)
DFT of \(z\), as in the legend
\(k\)
bin index
\(N\)
record length, samples
\((\cdot)^{*}\)
complex conjugate
\(j\)
imaginary unit, \(j^2=-1\)
Because \(Z^{*}[N-k] = X^{*}[N-k] - jY^{*}[N-k]\) \(= X[k] - jY[k]\). Two RF lines per FFT halve the transform count; the same idea computes one real N-point FFT with an N/2-point complex FFT.

The 2D DFT is rows, then columns

\[ \begin{aligned} X[k, l] &= \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} x[m, n]\, W_M^{km}\, W_N^{ln} \\ &= \sum_{m} W_M^{km}\Big(\sum_{n} x[m, n]\, W_N^{ln}\Big) \end{aligned} \]
(02.23)
\(x[m, n]\)
M × N array (e.g. depth samples × elements)
\(k,\ l\)
bin indices along the two dimensions
\(X[k,l]\)
2-D DFT
\(M,\ N\)
array size along the two axes (e.g. depth samples, channels)
\(m,\ n\)
sample indices along the two axes
\(W_M,\ W_N\)
twiddles \(e^{-j2\pi/M}\), \(e^{-j2\pi/N}\)
\(\sum\)
sum over the index shown
The kernel factors, so first transform every row over n, then every column over m. The cost scales as \(MN\log_2(MN)\) instead of \((MN)^2\). This is the k-space of Module 04 and the f-k domain of Module 24.
Pitfall

Plotting abs(fft(x)) against 0 … \(f_s\) and reporting the mirror peak at \(f_s - f_0\) as a second frequency. For real data the upper half is the negative-frequency copy. For complex I/Q data it is not a copy, and dropping it discards half the information (flow direction in Doppler).

02.10Visual summary

Each card is a question you should answer in five seconds, with the picture to keep in mind. The slider sets the tone frequency \(k_0\) (in bins, N = 16) for the first three cards: put it between two integers and watch the lines leak.

k₀ = 3.00 bins: on a bin, so all energy lands in bins 3 and 13.

02.11How experts think about it

An expert sees the DFT as N inner products with N orthogonal rotating arrows: each bin asks "how much of this frequency, at what phase", the real part is the cosine answer and the imaginary part the (negated) sine answer. Everything else follows from three facts: the probes are orthogonal only on the N-point grid (hence leakage), real signals produce conjugate pairs (hence the mirror and the \(AN/2\)), and the arrows for N and N/2 points are related by \(W_N^2 = W_{N/2}\) (hence the FFT).

A DFT bin is a correlation with one rotating arrow.
\(X[k] = \sum x[n]e^{-j2\pi kn/N}\)bin k ↔ \(k f_s/N\)

The probe \(e^{-j2\pi kn/N}\) turns k times over the record (02.4). Multiplying by it spins the signal backwards at frequency k; the component at k stops turning and its products add, everything else keeps turning and cancels. For x = [1, 2, 3, 4], \(X[1] = 1 - 2j - 3 + 4j = -2 + 2j\).

Exact formula
Only frequencies with a whole number of cycles in the record are probes; any other frequency spreads over all bins (leakage). Assumes uniform sampling.
The real part is the cosine content; the imaginary part is minus the sine content.
\(a\cos\theta + b\sin\theta \to \tfrac{N}{2}(a - jb)\)\(\sin(2\pi n/8) \to X[1] = -4j\)

Substitute \(x = a\cos\theta + b\sin\theta\) into \(X = \sum x\cos\theta - j\sum x\sin\theta\); the cross terms vanish by orthogonality and \(\sum\cos^2\theta = \sum\sin^2\theta = N/2\) (\(0 < k < N/2\)). So \(\mathrm{Re}\,X = aN/2\) and \(\mathrm{Im}\,X = -bN/2\) (02.10). A pure sine gives a pure negative-imaginary bin; a pure cosine a pure real one.

Exact formula
Holds with the \(e^{-j}\) forward convention and phase measured from sample n = 0. Move the time origin and Re and Im mix (a delay is a phase rotation). At k = 0 and N/2 there is no sine part.
A bin-centred tone of amplitude A shows up as A·N/2, split between bins k and N − k.
\(\lvert X[k_0]\rvert = AN/2\)DC: \(Nc\)Nyquist: \(NA\cos\varphi\)

\(A\cos\theta = \tfrac{A}{2}(e^{j\theta} + e^{-j\theta})\); each half-amplitude arrow correlates with its own probe for N samples, giving \((A/2)\cdot N\). A unit tone in a 1024-point FFT reads 512, i.e. 54.2 dB re 1. One-sided amplitude is therefore \(2\lvert X\rvert/N\).

Exact formula
Off-bin tones read up to 3.92 dB low (rectangular window). A window multiplies the peak by its coherent gain (Hann 0.5, so \(AN/4\)).
The upper half of the spectrum is the negative frequencies, mirrored and conjugated.
\(X[N-k] = X[k]^{*}\)Re even, Im odd

\(e^{-j2\pi(N-k)n/N} = e^{+j2\pi kn/N}\) at integer n, so bin N − k is frequency −k, and for real x its sum is the conjugate of bin k (02.11). For the 4-sample pulse at N = 8, \(X[1] = 1 - 2.414j\) and \(X[7] = 1 + 2.414j\). Only N/2 + 1 bins are independent.

Exact formula
Real input only. Complex I/Q data have independent negative frequencies (flow direction in Doppler); then all N bins carry information.
DC and Nyquist bins of a real signal are real, and a Nyquist tone has no recoverable phase.
probes 1 and \((-1)^n\)\(A\cos(\pi n + \varphi) = A\cos\varphi\,(-1)^n\)

The probe at k = N/2 is \(e^{-j\pi n} = (-1)^n\), purely real, and the sine at Nyquist is \(\sin(\pi n) = 0\) at every sample. Any sinusoid at exactly \(f_s/2\) is sampled at its zero crossings or peaks depending on \(\varphi\): at \(\varphi = 90°\) it vanishes completely.

Exact formula
Even N only; odd N has no Nyquist bin. Energy at exactly \(f_s/2\) is ambiguous, one reason RF is sampled well above \(2f_{\max}\).
The probes are orthogonal because the roots of unity sum to zero.
\(\sum e^{j2\pi dn/N} = N\delta[d]\)polygon closes

Geometric series: \((1 - r^N)/(1 - r)\) with \(r^N = 1\) (02.5). Tip to tail, the arrows walk around a closed polygon. This single identity gives the inverse DFT, Parseval, and the fact that an on-bin tone lights only its own bin.

Exact formula
Exact only for frequencies on the N-point grid. Two arbitrary frequencies over a finite record are not orthogonal; their overlap is the Dirichlet kernel.
The DFT is a rotation scaled by √N: energy is conserved up to a factor N.
\(\sum\lvert x\rvert^2 = \tfrac{1}{N}\sum\lvert X\rvert^2\)30 = 120/4

\(\mathbf{F}^{H}\mathbf{F} = N\mathbf{I}\), so \(\lVert\mathbf{X}\rVert^2 = \mathbf{x}^{H}\mathbf{F}^{H}\mathbf{F}\mathbf{x} = N\lVert\mathbf{x}\rVert^2\) (02.14). For x = [1, 2, 3, 4]: \(\lvert X\rvert^2\) = 100, 8, 4, 8, total 120 = 4 × 30. With the unitary scaling \(\mathbf{F}/\sqrt N\) the factor disappears.

Exact formula
Factor depends on the convention. A power spectral density per hertz needs \(\lvert X\rvert^2/(N f_s)\), and a window needs division by \(\sum w^2\) instead of N.
Coherent gain: a tone grows as N, noise as √N, so the FFT buys an SNR gain of N/2.
\(\mathrm{SNR}_{\text{bin}} = \tfrac{N}{2}\cdot\tfrac{A^2}{2\sigma^2}\)N = 1024: +27.1 dB

The tone's arrows add in phase: \(\lvert X\rvert^2 = A^2N^2/4\). Noise arrows add in random directions: \(\mathrm{E}\lvert X\rvert^2 = N\sigma^2\). The ratio is \(N/2\) times the per-sample SNR \(A^2/(2\sigma^2)\) (02.12). The figure: A = 0.5, σ = 1 (\(A^2/2\sigma^2\) = −9.0 dB per sample), N = 256, so the bin SNR is 128 × 0.125 = 16, i.e. 12.0 dB.

SimulationMonte Carlo, JS: one realization of white Gaussian noise (seeded) plus the tone, FFT in the page.
Needs a stable (coherent) tone over the whole record. Off-bin position, windowing (ENBW 1.5 bins for Hann) and coloured noise all reduce it.
An off-bin tone leaks: the probe and the tone drift apart and the sum curls up.
\(D_N(\delta) = \sin\pi\delta/\sin(\pi\delta/N)\)half bin: −3.92 dB

With \(\delta = k_0 - k\) the products turn \(\delta\) revolutions over the record (02.15). \(\delta = 0\): a straight line of length N. \(\delta = 1\): a closed circle, zero. \(\delta = \tfrac12\): a half circle, whose chord is \((2/\pi)\) × arc, so the bin reads \(20\log_{10}(2/\pi) = -3.92\) dB.

Exact formula
Large-N limit for a complex tone (N = 16 gives −3.91 dB); a real cosine adds its negative-frequency image (a few tenths of a dB at low k). Hann reduces the loss to 1.42 dB.
A window trades leakage for width: Hann buys −31.5 dB sidelobes with a main lobe twice as wide.
rect −13.3 dBHann −31.5 dBENBW 1.0 vs 1.5 bins

The window multiplies the record, so its transform convolves the spectrum. The rectangle's sharp ends give sidelobes that decay as \(1/\delta\); Hann's smooth ends (zero value and slope) give −31.5 dB first sidelobes decaying as \(1/\delta^3\), but its main lobe spans ±2 bins instead of ±1.

Exact formula
A window lowers the coherent gain (Hann 0.5) and raises the noise bandwidth. For a transient shorter than the record, the window just reshapes the pulse.
Zero padding interpolates the spectrum; only a longer record resolves.
resolution ≈ \(1/T = f_s/N\)grid \(f_s/(PN)\)

Padded zeros add nothing to \(\sum x[n]e^{-j\omega n}\); the padded DFT samples the same DTFT on a finer grid (02.17). Two equal, in-phase tones 0.5 bin apart in 16 samples stay one lobe at any padding (left); 2 bins apart they separate (right). Near 1/T the outcome depends on the tones' relative phase (two peaks for 25 % of phases at 0.5 bin, 58 % at 1 bin, all at 1.5 bins), and padding changes neither.

Exact formula
Padding does help locate an isolated peak, display a smooth spectrum, and avoid circular wrap in FFT convolution (pad to L + M − 1).
The FFT is the DFT split into even and odd samples, recursively.
\(X[k] = E[k] + W^kO[k]\)\(X[k+N/2] = E[k] - W^kO[k]\)

\(W_N^2 = W_{N/2}\) turns each half into an N/2-point DFT, and \(W_N^{k+N/2} = -W_N^k\) lets one product serve two outputs (02.20). For x = [1, 2, 3, 4]: E = [4, −2], O = [6, −2], X = [10, −2 + 2j, −2, −2 − 2j].

Schematic
Radix-2 needs \(N = 2^p\). Other lengths use mixed radix or Bluestein; libraries such as FFTW handle any N.
N log N, not N²: a 1024-point FFT needs about 200 times fewer multiplies.
1 048 576 vs 5 120ratio \(2N/\log_2 N\)

\(\log_2 N\) stages × N/2 butterflies × 1 multiply = \((N/2)\log_2 N\) (02.21). At N = 1024 that is 5 120 against 1 048 576, a ratio of 205; at N = 65 536 it is 8 192.

Exact formula
Real run time at large N is set by memory traffic, not multiplies. To get a few bins only, direct correlation (or the Goertzel recursion) costs N per bin and can win.
Decimation in time reads the input in bit-reversed order.
N = 8: 0 4 2 6 1 5 3 7

Each even/odd split sorts by the lowest remaining bit of n, so after \(\log_2 N\) splits position p holds sample bitrev(p): 1 = 001 goes to 100 = 4, 3 = 011 goes to 110 = 6. One permutation, then in-place butterflies with natural-order output.

Schematic
Decimation in frequency moves the reversal to the output; autosort (Stockham) variants avoid it at the cost of a second buffer.
fft returns bins 0 … N−1; above N/2 they are negative frequencies, and fftshift puts DC in the middle.
X(1) = DC\(f_k = k f_s/N,\ k < N/2\)\((k-N) f_s/N,\ k > N/2\)

For N = 8 at \(f_s\) = 40 MHz the bins are 0, 5, 10, 15, ±20, −15, −10, −5 MHz. After fftshift: −20, −15, …, 15 MHz, with axis (-floor(N/2):ceil(N/2)-1)*fs/N.

Exact formula
For odd N there is no Nyquist bin and ifftshift ≠ fftshift; always undo with ifftshift.
A delay is a linear phase, and the DFT's delays are circular.
\(x[n-m] \leftrightarrow W^{km}X[k]\)−360°·km/N

Shift the impulse to n = 1 and every bin becomes \(W^k\): magnitude 1, phase −45° per bin at N = 8. The slope of the phase is the delay. The 4-sample pulse centred at n = 1.5 has \(\angle X[1] = -67.5°\) (02.9).

Exact formula
The shift wraps around the record end. Fractional delays require a band-limited signal and the signed (−N/2 … N/2) bins with a real Nyquist bin.
Two real FFTs fit in one complex FFT, and a 2D FFT is rows, then columns.
\(X = (Z[k] + Z^{*}[N-k])/2\)fft2 = fft over rows, then columns

With \(z = x + jy\), the conjugate-symmetric part of Z is X and the antisymmetric part is jY (02.22). Two RF lines per transform halve the FFT count. The 2D kernel \(W_M^{km}W_N^{ln}\) factors (02.23), so a 256 × 256 fft2 costs 2 × 256 FFTs of length 256, the same \((N/2)\log_2 N\) count as one 65 536-point FFT.

The two-for-one trick costs some extra additions and loses accuracy when one signal is far weaker than the other. Separability needs a rectangular grid.

02.12Cheat sheet

ConceptFormulaUltrasound use
DFT / inverse\(X[k] = \sum x[n]W^{kn}\), \(x[n] = \frac{1}{N}\sum X[k]W^{-kn}\), \(W = e^{-j2\pi/N}\)RF spectra, filtering, k-space
Meaning of Re, Im\(\mathrm{Re}\,X = \sum x\cos\), \(\mathrm{Im}\,X = -\sum x\sin\); \(a\cos\theta + b\sin\theta \to \frac{N}{2}(a - jb)\)I and Q content at each frequency
Amplitude, phase\(A = 2\lvert X[k]\rvert/N\) (\(0 < k < N/2\)), \(\varphi = \angle X[k]\)Spectral amplitude for QUS; phase for delay
Orthogonality\(\sum_n e^{j2\pi(k-l)n/N} = N\,\delta[k-l]\)Why bins separate; inverse DFT
Conjugate symmetry\(X[N-k] = X[k]^{*}\) (real x); \(X[0]\), \(X[N/2]\) realKeep N/2 + 1 bins of real RF; not for I/Q
Parseval\(\sum\lvert x\rvert^2 = \frac{1}{N}\sum\lvert X\rvert^2\)Energy and power spectra
Coherent gain\(\mathrm{SNR}_{\text{bin}} = \frac{N}{2}\) × per-sample SNRDetecting weak harmonics, Doppler tones
Leakage kernel\(D_N(\delta) = \sin(\pi\delta)/\sin(\pi\delta/N)\); half bin −3.92 dB (Hann −1.42 dB)Bias in spectral centroid and attenuation
Resolution vs gridresolution ≈ \(f_s/N = 1/T\); grid \(f_s/(PN)\) with paddingGate length sets spectral resolution
Radix-2 butterfly\(X[k] = E + W^kO\), \(X[k + N/2] = E - W^kO\)Every FFT in the scanner
FFT cost\((N/2)\log_2 N\) multiplies vs \(N^2\)Real-time spectral Doppler, Fourier beamforming
Delay theorem\(x[n-m] \leftrightarrow e^{-j2\pi km/N}X[k]\)Fractional delays, aberration phase screens
Frequency axisf = (0:N-1)*fs/N; centred (-floor(N/2):ceil(N/2)-1)*fs/N with fftshiftCorrect MHz labels on every plot
Two-for-one\(X = \frac{Z[k] + Z^{*}[N-k]}{2}\), \(Y = \frac{Z[k] - Z^{*}[N-k]}{2j}\)Two RF lines per complex FFT

02.13Interview and qualifying-exam questions

I1What does the imaginary part of an FFT bin physically mean? · screen

It is minus the correlation of the signal with a sine at that bin's frequency, i.e. the quadrature content. For real x, x[n] contains \(\frac{2}{N}[\mathrm{Re}\,X[k]\cos(2\pi kn/N)\) \(- \mathrm{Im}\,X[k]\sin(2\pi kn/N)]\), so \(-\frac{2}{N}\mathrm{Im}\,X[k]\) is the sine amplitude (02.10). A pure sine at bin 1 of an 8-point FFT gives \(X[1] = -4j\). It is not "imaginary" in any physical sense: it is the 90°-shifted projection, and it changes if you move the time origin. Follow-up: what is \(\mathrm{Im}\,X[0]\)? (Zero for real x; there is no sine at DC.)

I2A 1 V amplitude sinusoid sits exactly on bin 100 of a 1024-point FFT. What does abs(X(101)) read, and how do you convert to volts? · screen

\(AN/2 = 512\). Divide by N/2 for the one-sided amplitude (1 V). With a Hann window it reads 256 (coherent gain 0.5). If the tone sits half a bin off, the rectangular peak drops to \(512 \times 2/\pi \approx 326\) (−3.92 dB). Follow-up: where is the other half? (Bin 924, the conjugate.)

I3Why do the upper N/2 bins mirror the lower ones for real data, and when do they not? · screen

\(e^{-j2\pi(N-k)n/N} = e^{+j2\pi kn/N}\), so bin N − k is frequency −k, and for real x it equals \(X[k]^{*}\). A real cosine is two counter-rotating phasors, so each appears twice with half amplitude. Complex I/Q data have no such symmetry: after demodulation the upper bins are genuine negative frequencies, which encode flow direction in Doppler.

I4Derive the radix-2 FFT butterfly. · onsite

Split n into even and odd: \(X[k] = \sum x[2m]W_N^{2mk}\) \(+ W_N^k\sum x[2m+1]W_N^{2mk}\). Use \(W_N^2 = W_{N/2}\), so the sums are N/2-point DFTs \(E[k]\), \(O[k]\). They have period N/2 and \(W_N^{k+N/2} = -W_N^k\), so \(X[k] = E + W^kO\) and \(X[k + N/2] = E - W^kO\). Recursion gives \(\log_2 N\) stages of N/2 butterflies: \((N/2)\log_2 N\) multiplies. Follow-up: why is the input bit-reversed? (Each split sorts by the lowest remaining bit.)

I5Does zero padding improve frequency resolution? · screen

No. It samples the same DTFT more finely; the main-lobe width is set by the record length, ≈ \(1/T = f_s/N\). Padding a 2 µs gate (80 samples at 40 MHz) to 256 gives a 156 kHz grid, but two components closer than about 500 kHz cannot be reliably separated. It does improve peak location by interpolation and prevents wrap-around in FFT convolution. Follow-up: what does improve resolution? (A longer gate, or model-based methods such as AR/Burg with their own bias, Module 27.)

I6What is the scalloping loss, and how would you remove it for amplitude measurements? · onsite

The peak-bin amplitude deficit for a tone between bins: \(\lvert D_N(\tfrac12)\rvert/N \approx 2/\pi\), −3.92 dB, rectangular window; −1.42 dB with Hann. Remedies: a flat-top window (wide main lobe, small amplitude error), zero padding and peak interpolation (parabolic on log magnitude), or fitting the known kernel. Follow-up: what does that cost in noise bandwidth? (Wider windows raise ENBW, 1.5 bins for Hann, more for flat-top.)

I7A tone is 20 dB below the per-sample noise. How long an FFT do you need to see it at 10 dB SNR in its bin? · onsite

Bin SNR = \((N/2)\) × per-sample SNR, so \(N/2 = 10^{(10 + 20)/10} = 1000\), N ≈ 2000; take 2048 (gain 30.1 dB). That assumes the tone is coherent over the whole record, sits on a bin, and the noise is white. Follow-up: with Hann? (ENBW 1.5 costs 1.76 dB, so go to 4096 or accept 8.3 dB.)

I8How do you build the frequency axis for fftshift(fft(x)) with N odd? · screen

f = (-floor(N/2):ceil(N/2)-1)*fs/N, which for N = 7 is \((-3 \dots 3)\cdot f_s/7\): no Nyquist bin. For N even it is −N/2 … N/2 − 1, with \(-f_s/2\) first. Undo with ifftshift, not fftshift, which differ by one sample for odd N. Follow-up: where is DC after the shift? (Index floor(N/2) + 1 in MATLAB.)

I9Show that the DFT's inverse is \(\frac{1}{N}\mathbf{F}^{H}\) without inverting anything. · research

\((\mathbf{F}^{H}\mathbf{F})_{lk} = \sum_n e^{j2\pi(l-k)n/N} = N\,\delta[l-k]\) by the geometric series with \(r^N = 1\). So \(\mathbf{F}^{H}\mathbf{F} = N\mathbf{I}\), \(\mathbf{F}\) is \(\sqrt N\) times a unitary matrix, and Parseval follows. Follow-up: what are the eigenvalues of \(\mathbf{F}/\sqrt N\)? (Only ±1 and ±j: \(\mathbf{F}^4 = N^2\mathbf{I}\).)

I10You need the spectrum of 128 RF lines of 2048 samples each, real-time. What FFT tricks do you use? · research

Pack pairs of lines as \(x + jy\) and split with the conjugate-symmetry formulas (02.22) (64 complex FFTs instead of 128); keep only N/2 + 1 bins; batch the transforms (GPU or FFTW plans reused across lines); if only a band is needed, demodulate and decimate first (Module 08), so N shrinks by the decimation factor. Cost per line ≈ \((N/2)\log_2 N\) = 11 264 complex multiplies vs 4.2 million direct. Follow-up: how does the 2D FFT of the whole frame relate? (Rows then columns: k-space, Module 04.)

02.14MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.

Download the scripts: m02_dft_by_hand.m · m02_fft_radix2.m · m02_leakage_zero_padding.m · m02_fft_practice.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m02_dft_by_hand.mN = 4 and N = 8 DFTs by loops, every multiply printed, vs matrix and fft; Re/Im interpretation · change x · runtime ≈ 1 s 
%% m02_dft_by_hand.m  -- Module 02: the DFT by hand for N = 4 and N = 8
% Purpose : compute X[k] = sum_n x[n] exp(-j*2*pi*k*n/N) three ways (explicit loops with every
%           multiply printed, the DFT matrix F, and fft), then read each bin as cosine content
%           (real part), minus sine content (imaginary part), amplitude and phase.
% Change  : x4 (the N = 4 example) and x8 (the N = 8 example).
% Expect  : x4 = [1 2 3 4] -> X = [10, -2+2j, -2, -2-2j]; x[n] = 2.5 - cos(pi n/2) - sin(pi n/2)
%           - 0.5(-1)^n; Parseval 30 = 120/4. Pulse [1 1 1 1 0 0 0 0] -> X[1] = 1 - 2.414j,
%           |X[1]| = 2.613, angle -67.5 deg, X[2] = X[4] = X[6] = 0. Loops, matrix and fft agree
%           to about 1e-15.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);
cfmt = @(z) sprintf('%7.3f %s %6.3fj', real(z), char(43 + 2*(imag(z) < 0)), abs(imag(z)));  % '+' or '-'

%% 1. N = 4: every multiply
x4 = [1 2 3 4];
N = numel(x4);
X_loop = zeros(1, N);
fprintf('--- N = 4, x = [%s] ---\n', num2str(x4));
fprintf('W = exp(-j*2*pi/4) = -j, so every probe value is 1, -j, -1 or j\n');
for k = 0:N-1
    terms = '';
    for n = 0:N-1
        W = exp(-1j*2*pi*mod(k*n, N)/N);                   % probe value W^(kn)
        W = round(W*1e12)/1e12;                            % clean -0 and 1e-16 for printing
        p = x4(n+1)*W;                                     % one multiply
        X_loop(k+1) = X_loop(k+1) + p;
        terms = [terms sprintf('  %d*(%s)', x4(n+1), strtrim(cfmt(W)))]; %#ok<AGROW>
    end
    fprintf('k = %d:%s\n        X[%d] = %s\n', k, terms, k, cfmt(X_loop(k+1)));
end

%% 2. The DFT matrix and fft
[kk, nn] = ndgrid(0:N-1, 0:N-1);
F4 = exp(-1j*2*pi*kk.*nn/N);                               % F(k+1, n+1) = W^(kn)
X_mat = (F4*x4(:)).';
X_fft = fft(x4);
fprintf('max |loop - fft| = %.1e, max |matrix - fft| = %.1e\n', max(abs(X_loop - X_fft)), max(abs(X_mat - X_fft)));
fprintf('F4''*F4 = %g * I (max off-diagonal %.1e)\n', real(trace(F4'*F4))/N, max(max(abs(F4'*F4 - N*eye(N)))));

%% 3. What each bin means
fprintf('\n--- meaning of each bin (N = 4) ---\n');
fprintf('X[0]/N = %.3f   (mean, DC)\n', real(X_fft(1))/N);
fprintf('bin 1: cos amplitude 2*Re/N = %.3f, sin amplitude -2*Im/N = %.3f, |X| = %.3f, angle = %.1f deg\n', ...
        2*real(X_fft(2))/N, -2*imag(X_fft(2))/N, abs(X_fft(2)), angle(X_fft(2))*180/pi);
fprintf('X[2]/N = %.3f   (Nyquist, multiplies (-1)^n; always real for real x)\n', real(X_fft(3))/N);
fprintf('X[3] = %s = conj(X[1]) (negative frequency)\n', cfmt(X_fft(4)));
n = 0:N-1;
x_rebuilt = real(X_fft(1))/N + (2/N)*(real(X_fft(2))*cos(2*pi*n/N) - imag(X_fft(2))*sin(2*pi*n/N)) + real(X_fft(3))/N*(-1).^n;
fprintf('rebuilt from the pieces: [%s]\n', num2str(x_rebuilt, '%6.3f'));
fprintf('Parseval: sum|x|^2 = %g, sum|X|^2/N = %g\n', sum(abs(x4).^2), sum(abs(X_fft).^2)/N);

%% 4. N = 8: the 4-sample pulse, table for bin 1
x8 = [1 1 1 1 0 0 0 0];
N = numel(x8); k = 1;
fprintf('\n--- N = 8, x = [%s], bin k = %d ---\n', num2str(x8), k);
fprintf('  n  x[n]  theta(deg)   cos     -sin    x*cos   -x*sin   sum Re   sum Im\n');
sre = 0; sim = 0;
for n = 0:N-1
    th = 2*pi*mod(k*n, N)/N;
    sre = sre + x8(n+1)*cos(th); sim = sim - x8(n+1)*sin(th);
    fprintf('  %d  %4.1f  %6.0f   %7.3f %7.3f %7.3f %7.3f  %7.3f  %7.3f\n', n, x8(n+1), th*180/pi, ...
            cos(th), -sin(th), x8(n+1)*cos(th), -x8(n+1)*sin(th), sre, sim);
end
X8 = fft(x8);
fprintf('X[1] = %s, |X[1]| = %.3f, angle = %.1f deg (pulse centred at n = 1.5: -2*pi*1.5/8 = %.1f deg)\n', ...
        cfmt(X8(2)), abs(X8(2)), angle(X8(2))*180/pi, -2*pi*1.5/8*180/pi);
fprintf('all 8 bins:\n');
for kk8 = 0:N-1
    fprintf('  X[%d] = %s   |X| = %.4f\n', kk8, cfmt(X8(kk8+1)), abs(X8(kk8+1)));
end
[kk, nn] = ndgrid(0:N-1, 0:N-1); F8 = exp(-1j*2*pi*kk.*nn/N);
fprintf('max |F8*x - fft(x)| = %.1e\n', max(abs((F8*x8(:)).' - X8)));

%% 5. Figures
figure('Position', [100 100 900 600]);
subplot(2, 2, 1); stem(0:7, x8, 'filled'); xlabel('sample n'); ylabel('x[n]'); title('4-sample pulse, N = 8'); xlim([-0.5 7.5]);
subplot(2, 2, 2); stem((0:7) - 0.1, real(X8), 'filled'); hold on; stem((0:7) + 0.1, imag(X8), 'filled'); hold off;
xlabel('bin k'); ylabel('X[k]'); legend('Re X (cosine content)', 'Im X (-sine content)'); title('Re even, Im odd about k = 4'); xlim([-0.5 7.5]);
subplot(2, 2, 3); imagesc(0:7, 0:7, real(F8)); axis image; colorbar; xlabel('n'); ylabel('k'); title('Re F_8 = cos(2\pi kn/8)');
subplot(2, 2, 4); imagesc(0:7, 0:7, imag(F8)); axis image; colorbar; xlabel('n'); ylabel('k'); title('Im F_8 = -sin(2\pi kn/8)');

What to look for: X = [10, −2+2j, −2, −2−2j] for x = [1 2 3 4]; the reconstruction 2.5 − cos − sin − 0.5(−1)n; X[1] = 1 − 2.414j for the 8-sample pulse; loop, matrix and fft agree to ≈ 10⁻¹⁵; Parseval 30 = 120/4.

m02_fft_radix2.mRecursive and iterative radix-2 FFT checked against fft, bit reversal, timing and operation counts · change N range · runtime ≈ 2 s 
%% m02_fft_radix2.m  -- Module 02: radix-2 FFT, recursive and iterative, checked against fft
% Purpose : implement X[k] = E[k] + W^k O[k], X[k+N/2] = E[k] - W^k O[k] two ways:
%           (a) recursive, written with anonymous functions only (no local functions), and
%           (b) iterative in place: bit-reversal permutation, then log2(N) butterfly stages.
%           Check both against fft, show the bit-reversed order and time them against the
%           direct O(N^2) DFT.
% Change  : Nlist (powers of two), the test signal.
% Expect  : bit-reversed order for N = 8: 0 4 2 6 1 5 3 7; N = 4 by hand E = [4 -2], O = [6 -2],
%           X = [10, -2+2j, -2, -2-2j]; errors vs fft ~1e-13 or smaller; direct DFT time grows
%           ~4x per doubling of N, FFTs ~2x; multiplies (N/2)log2(N) = 5120 vs N^2 = 1048576
%           at N = 1024 (ratio 204.8).
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

%% 1. Recursive radix-2 FFT with anonymous functions
% A cell array of two handles plays the role of if/else: only the selected branch is called.
tw     = @(N) exp(-2j*pi*(0:N/2-1).'/N);                      % twiddle factors W_N^k, k = 0..N/2-1
comb   = @(E, O, N) [E + tw(N).*O; E - tw(N).*O];             % one level of butterflies
branch = {@(f, x) x, @(f, x) comb(f(f, x(1:2:end)), f(f, x(2:2:end)), numel(x))};
rfft   = @(f, x) branch{(numel(x) > 1) + 1}(f, x);            % 1-point DFT is the sample itself
fft_rec = @(x) rfft(rfft, x(:));

x4 = [1; 2; 3; 4];
E = fft([1; 3]); O = fft([2; 4]);
fprintf('--- N = 4 by hand ---\n');
fprintf('E = DFT[1 3] = [%g %g], O = DFT[2 4] = [%g %g], W4^1 = -j\n', E, O);
X4 = fft_rec(x4);
fprintf('recursive FFT: X = [%s]\n', num2str(X4.', '%6.2f'));

%% 2. Iterative in-place radix-2 FFT
bitrev = @(N) bin2dec(fliplr(dec2bin(0:N-1, round(log2(N))))) + 1;   % 1-based bit-reversed index
fprintf('\n--- bit reversal ---\nN = 8 input order: %s\n', num2str((bitrev(8) - 1).'));
fprintf('N = 16 input order: %s\n', num2str((bitrev(16) - 1).'));

%% 3. Accuracy against fft
Nlist = 2.^(6:12);
nN = numel(Nlist);
t_dft = nan(1, nN); t_rec = zeros(1, nN); t_it = zeros(1, nN); t_fft = zeros(1, nN);
err_rec = zeros(1, nN); err_it = zeros(1, nN);
fprintf('\n     N   err(recursive)  err(iterative)  t_dft(ms)  t_rec(ms)  t_iter(ms)  t_fft(ms)\n');
for i = 1:nN
    N = Nlist(i); L = round(log2(N));
    x = randn(N, 1) + 1j*randn(N, 1);
    Xref = fft(x);
    % recursive
    tic; Xr = fft_rec(x); t_rec(i) = toc;
    % iterative: permute once, then stages; group the butterflies of each stage as columns
    tic;
    X = x(bitrev(N));
    for s = 1:L
        m = 2^s; h = m/2;
        w = exp(-2j*pi*(0:h-1).'/m);
        Xs = reshape(X, m, N/m);                     % each column is one group of length m
        a = Xs(1:h, :); b = w .* Xs(h+1:m, :);       % implicit expansion of w over the groups
        Xs = [a + b; a - b];
        X = Xs(:);
    end
    t_it(i) = toc;
    % direct DFT by matrix (only up to N = 2048 to keep memory small)
    if N <= 2048
        tic; [kk, nn] = ndgrid(0:N-1, 0:N-1); Xd = exp(-2j*pi*kk.*nn/N)*x; t_dft(i) = toc; %#ok<NASGU>
    end
    tic; for r = 1:20, Xf = fft(x); end; t_fft(i) = toc/20; %#ok<NASGU>
    err_rec(i) = max(abs(Xr - Xref))/max(abs(Xref));
    err_it(i)  = max(abs(X  - Xref))/max(abs(Xref));
    fprintf('%6d   %12.1e   %12.1e   %9.2f  %9.2f  %10.2f  %9.3f\n', N, err_rec(i), err_it(i), ...
            1e3*t_dft(i), 1e3*t_rec(i), 1e3*t_it(i), 1e3*t_fft(i));
end

%% 4. Operation counts
fprintf('\n     N        N^2   (N/2)log2N   ratio 2N/log2N\n');
for N = [8 256 1024 4096 65536]
    fprintf('%6d %10d %12d %12.1f\n', N, N^2, N/2*log2(N), N^2/(N/2*log2(N)));
end

%% 5. Figure: timing and counts
figure('Position', [100 100 900 380]);
subplot(1, 2, 1);
loglog(Nlist, 1e3*t_dft, 'o-', Nlist, 1e3*t_rec, 's-', Nlist, 1e3*t_it, 'd-', Nlist, 1e3*t_fft, '^-');
xlabel('N'); ylabel('time (ms)'); grid on; legend('direct DFT (matrix)', 'recursive radix-2', 'iterative radix-2', 'built-in fft', 'Location', 'northwest');
title('Run time (interpreter overhead dominates the hand-written FFTs)');
subplot(1, 2, 2);
Nc = 2.^(2:16);
loglog(Nc, Nc.^2, 'o-', Nc, Nc/2.*log2(Nc), 's-'); grid on;
xlabel('N'); ylabel('complex multiplies'); legend('N^2', '(N/2) log_2 N', 'Location', 'northwest'); title('Operation count');

What to look for: bit-reversed order 0 4 2 6 1 5 3 7; both implementations match fft to ≈ 10⁻¹³ or better; the direct DFT's time grows ≈ 4× per doubling of N, the hand-written FFTs ≈ 2× or less (interpreter overhead); multiply ratio 205 at N = 1024.

m02_leakage_zero_padding.mOff-bin tone, scalloping loss rect vs Hann, Dirichlet kernel, zero padding vs resolution · change k0, P · runtime ≈ 1 s 
%% m02_leakage_zero_padding.m  -- Module 02: leakage, scalloping, Dirichlet kernel, zero padding
% Purpose : show that a tone between bins spreads over every bin (leakage), that the bins sample
%           the Dirichlet kernel D_N(delta) = sin(pi*delta)/sin(pi*delta/N), measure the scalloping
%           loss for rectangular and Hann windows, and show that zero padding samples the same
%           DTFT more finely without improving resolution.
% Change  : N, k0 (tone frequency in bins), the padding factor P, the tone separations sep.
% Expect  : complex tone, half-bin offset: rect -3.92 dB, Hann -1.42 dB; measured |X[k]|/N equals
%           |D_N(delta)|/N to ~1e-15; padded bins coincide with the DTFT; two tones 0.7 bin apart
%           show one peak at P = 1, 4, 16; two tones 2 bins apart show two peaks. Near 1/T the
%           outcome depends on the relative phase of the tones (section 4 sweeps it).
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

N  = 64;                                   % record length (samples)
n  = (0:N-1).';
hann = 0.5 - 0.5*cos(2*pi*n/N);            % periodic Hann, written out (no toolbox)
dirichlet = @(d) (abs(d) < 1e-12) + (abs(d) >= 1e-12).*sin(pi*d)./(N*sin(pi*d/N) + (abs(d) < 1e-12));  % D_N(d)/N, = 1 at d = 0

%% 1. Scalloping loss: complex tone at an integer bin and half a bin off
fprintf('--- scalloping loss (complex tone, N = %d) ---\n', N);
for k0 = [10 10.25 10.5]
    z = exp(2j*pi*k0*n/N);
    Xr = fft(z);        Xh = fft(z.*hann);
    pr = max(abs(Xr))/N;                   % rect: on-bin peak = N
    ph = max(abs(Xh))/(N*mean(hann));      % Hann: on-bin peak = N * coherent gain (0.5)
    fprintf('k0 = %5.2f: peak bin rect %6.2f dB (theory %6.2f), Hann %6.2f dB\n', k0, 20*log10(pr), ...
            20*log10(abs(dirichlet(k0 - round(k0)))), 20*log10(ph));
end
fprintf('theory at half bin: 20log10(2/pi) = %.2f dB (rect)\n', 20*log10(2/pi));

%% 2. The bins sample the Dirichlet kernel
k0 = 10.3; z = exp(2j*pi*k0*n/N); X = fft(z);
k = (0:N-1).'; d = k0 - k;
fprintf('\n--- Dirichlet check, k0 = %.2f ---\n', k0);
fprintf('max | |X[k]|/N - |D_N(k0-k)|/N | = %.1e\n', max(abs(abs(X)/N - abs(dirichlet(d)))));
fprintf('bins 9..12: |X|/N = %s\n', num2str(abs(X(10:13)).'/N, '%8.4f'));
fprintf('fraction of energy outside bins 9..11: %.1f %%\n', 100*(1 - sum(abs(X(10:12)).^2)/sum(abs(X).^2)));

%% 3. Zero padding: the same DTFT on a finer grid
x = cos(2*pi*k0*n/N);                      % real tone this time
P = 8;
Xp = fft(x, P*N);                          % pads with (P-1)N zeros
fd = (0:P*N-1)/P;                          % frequency in bins of the original record
Xd = zeros(size(fd));                      % DTFT evaluated directly
for i = 1:numel(fd), Xd(i) = sum(x.*exp(-2j*pi*fd(i)*n/N)); end
fprintf('\n--- zero padding ---\nmax |padded FFT - direct DTFT| = %.1e\n', max(abs(Xp(:).' - Xd)));
[~, ip] = max(abs(Xp(1:P*N/2))); X1 = fft(x); [~, i1] = max(abs(X1(1:N/2)));
fprintf('padded peak at %.3f bins (true %.2f); unpadded peak at bin %d\n', fd(ip), k0, i1 - 1);

%% 4. Resolution is set by N, not by padding
fprintf('\n--- two tones: padding vs resolution (N = 16) ---\n');
Ns = 16; ns = (0:Ns-1).';
cnt = @(v) sum(v(2:end-1) > v(1:end-2) & v(2:end-1) > v(3:end) & v(2:end-1) > 0.5*max(v));   % strong local maxima
for sep = [0.5 2.0]
    xs = cos(2*pi*4*ns/Ns) + cos(2*pi*(4 + sep)*ns/Ns);
    msg = sprintf('separation %.1f bin:', sep);
    for Pp = [1 4 16]
        V = abs(fft(xs, Pp*Ns)); V = V(1:Pp*Ns/2 + 1);
        msg = [msg sprintf('  P = %2d -> %d peak(s)', Pp, cnt(V(:)))]; %#ok<AGROW>
    end
    fprintf('%s\n', msg);
end
fprintf('relative phase sweep (P = 16, 24 phases): fraction of phases with two peaks\n');
for sep = [0.5 1.0 1.5 2.0]
    nres = 0;
    for ph = (0:23)*2*pi/24
        xs = cos(2*pi*4*ns/Ns) + cos(2*pi*(4 + sep)*ns/Ns + ph);
        V = abs(fft(xs, 16*Ns)); V = V(1:8*Ns + 1); nres = nres + (cnt(V(:)) >= 2);
    end
    fprintf('  separation %.1f bin: %4.0f %%\n', sep, 100*nres/24);
end

%% 5. Figures
figure('Position', [100 100 900 600]);
subplot(2, 1, 1);
plot(fd, 20*log10(abs(Xd)/(N/2) + 1e-12), 'Color', [0.5 0.5 0.5]); hold on;
Xb = 20*log10(abs(X1(:))/(N/2) + 1e-12);
stem(0:N/2, Xb(1:N/2+1), 'filled', 'BaseValue', -80);
plot(fd(1:P*N/2), 20*log10(abs(Xp(1:P*N/2))/(N/2) + 1e-12), '.'); hold off;
xlim([0 N/2]); ylim([-80 5]); grid on;
xlabel('frequency (bins of the N-sample record)'); ylabel('dB re A N/2');
title(sprintf('Real tone at k_0 = %.1f: DTFT (grey), N-point bins (stems), %dx padded (dots)', k0, P));
subplot(2, 1, 2);
dd = linspace(-4, 4, 801);
plot(dd, 20*log10(abs(dirichlet(dd + 1e-9))), dd, 20*log10(abs(arrayfun(@(q) sum(hann.*exp(2j*pi*q*n/N)), dd))/(N/2)));
ylim([-80 3]); grid on; xlabel('offset \delta from the bin (bins)'); ylabel('dB');
legend('rectangular |D_N(\delta)|/N', 'Hann'); title('Window kernels: first sidelobe -13.3 dB vs -31.5 dB');

What to look for: complex-tone scalloping −3.92 dB (rect) and −1.42 dB (Hann); measured Dirichlet magnitudes match sin(πδ)/sin(πδ/N); padded bins lie on the same curve; two equal, in-phase tones 0.5 bin apart remain one lobe at every padding, 2 bins apart they separate; the phase sweep gives two peaks for 25 %, 58 % and 100 % of relative phases at 0.5, 1.0 and 1.5 bins.

m02_fft_practice.mFrequency axes and fftshift (odd/even N), one-sided amplitude in MHz, two-for-one real FFT, fft2 = rows then columns, coherent gain with noise · runtime ≈ 1 s 
%% m02_fft_practice.m  -- Module 02: FFTs in practice
% Purpose : build frequency axes for fft and fftshift output (even and odd N), read a one-sided
%           amplitude spectrum in MHz, compute two real FFTs with one complex FFT, verify that
%           fft2 is fft over rows then columns, and measure the coherent processing gain N/2.
% Change  : fs, f0, N, the noise level sigma, the tone amplitude A.
% Expect  : 5 MHz pulse peaks at 5.00 MHz on both axes; fftshift puts DC at index floor(N/2)+1;
%           ifftshift(fftshift(X)) == X for odd N while fftshift(fftshift(X)) does not;
%           two-for-one and rows-then-columns errors ~1e-14; bin SNR within a few tenths of a dB
%           of 10log10(N/2) + per-sample SNR (A = 0.5, sigma = 1: -9.03 + 27.09 = 18.06 dB at N = 1024).
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

%% 1. Frequency axes for an RF pulse
fs = 40;                         % sampling rate (MHz)
f0 = 5;                          % centre frequency (MHz)
N  = 1024;                       % FFT length
t  = ((0:N-1) - N/2)/fs;         % time (us), pulse centred in the record
x  = exp(-(t/0.25).^2) .* cos(2*pi*f0*t);   % Gaussian-envelope RF pulse, 1/e half-width 0.25 us
X  = fft(x);
f_fft   = (0:N-1)*fs/N;                          % fft order
f_shift = (-floor(N/2):ceil(N/2)-1)*fs/N;        % after fftshift
Xs = fftshift(X);
[~, i1] = max(abs(X(1:N/2+1)));  [~, i2] = max(abs(Xs) .* (f_shift > 0));
fprintf('--- frequency axes (N = %d, fs = %g MHz, df = %.2f kHz) ---\n', N, fs, 1e3*fs/N);
fprintf('peak on fft axis: %.2f MHz (index %d); on shifted axis: %.2f MHz\n', f_fft(i1), i1, f_shift(i2));
fprintf('DC after fftshift sits at index %d = floor(N/2)+1\n', find(f_shift == 0));
A1 = abs(X(1:N/2+1))/N; A1(2:end-1) = 2*A1(2:end-1);   % one-sided amplitude
fprintf('one-sided amplitude peak = %.4f (pulse spectrum, not a tone)\n', max(A1));

%% 2. Odd N: fftshift vs ifftshift
No = 7; v = 0:No-1;
fprintf('\n--- odd N = %d ---\n', No);
fprintf('fftshift(0:6)            = [%s]\n', num2str(fftshift(v)));
fprintf('centred axis (bins)      = [%s]\n', num2str(-floor(No/2):ceil(No/2)-1));
fprintf('ifftshift(fftshift(v))   = [%s]  (restores)\n', num2str(ifftshift(fftshift(v))));
fprintf('fftshift(fftshift(v))    = [%s]  (does not)\n', num2str(fftshift(fftshift(v))));

%% 3. Two real FFTs with one complex FFT
a = randn(1, N); b = randn(1, N);
Z = fft(a + 1j*b);
Zc = conj(Z([1, N:-1:2]));                       % Z*[N-k], indices mod N
Xa = (Z + Zc)/2; Xb = (Z - Zc)/(2j);
fprintf('\n--- two-for-one ---\nmax |Xa - fft(a)| = %.1e, max |Xb - fft(b)| = %.1e\n', max(abs(Xa - fft(a))), max(abs(Xb - fft(b))));

%% 4. fft2 = rows then columns
img = randn(64, 128);
R = fft(img, [], 2);                             % every row
C = fft(R, [], 1);                               % then every column
fprintf('max |fft2 - rows-then-columns| = %.1e\n', max(max(abs(fft2(img) - C))));

%% 5. Coherent processing gain
A = 0.5; sigma = 1; trials = 400;           % tone at fs/8, i.e. bin N/8 (always on a bin)
snr_sample = 10*log10(A^2/2/sigma^2);
fprintf('\n--- coherent gain (A = %.1f, sigma = %.0f, per-sample SNR %.2f dB) ---\n', A, sigma, snr_sample);
fprintf('     N   theory bin SNR (dB)   simulated (dB)\n');
for Nn = [64 256 1024 4096]
    nn = 0:Nn-1; ks = Nn/8 + 1;                          % MATLAB index of bin N/8
    s = A*cos(2*pi*(Nn/8)*nn/Nn);
    P_noise = 0;
    for r = 1:trials
        Wn = fft(sigma*randn(1, Nn));
        P_noise = P_noise + mean(abs(Wn([2:Nn/2])).^2);   % noise power per bin (positive half, no DC)
    end
    P_noise = P_noise/trials;
    Ssig = abs(fft(s)); P_tone = Ssig(ks)^2;
    fprintf('%6d   %18.2f   %14.2f\n', Nn, snr_sample + 10*log10(Nn/2), 10*log10(P_tone/P_noise));
end

%% 6. Figure
figure('Position', [100 100 900 600]);
subplot(2, 1, 1); plot(f_fft, abs(X)); xlim([0 fs]); grid on;
xlabel('frequency (MHz), fft order'); ylabel('|X|'); title('fft order: positive half, then the mirrored negative frequencies');
subplot(2, 1, 2); plot(f_shift, abs(Xs)); xlim([-fs/2 fs/2]); grid on;
xlabel('frequency (MHz), after fftshift'); ylabel('|X|'); title('fftshift: DC in the middle, peaks at \pm5 MHz');

What to look for: the 5 MHz pulse peaks at 5.0 MHz on the shifted axis; two-for-one and rows-then-columns errors below 10⁻¹³; simulated bin SNR within a few tenths of a dB of 10 log₁₀(N/2) + per-sample SNR.

02.15Self-check

Q1Compute the 4-point DFT of x = [1, 0, −1, 0] by hand. What does each bin mean?

\(X[0] = 0\), \(X[1] = 1 + 0 + (-1)(-1) + 0 = 2\), \(X[2] = 1 - 0 + (-1)(1) - 0 = 0\), \(X[3] = 2\). x is \(\cos(2\pi n/4)\) with A = 1: two real lines of \(AN/2 = 2\) at bins 1 and 3, no DC, no Nyquist, no sine part (Im = 0).

Q2An 8-point DFT of a real signal has \(X[3] = 2 - 6j\). What is \(X[5]\), and what sinusoid does the pair represent?

\(X[5] = X[3]^{*} = 2 + 6j\). The pair is \(\frac{2}{8}[2\cos(2\pi 3n/8) + 6\sin(2\pi 3n/8)]\) \(= 0.5\cos(3\pi n/4) + 1.5\sin(3\pi n/4)\), amplitude \(\frac{2}{8}\lvert 2 - 6j\rvert = 0.25 \times 6.32 = 1.58\), phase \(\angle X[3] = -71.6°\).

Q3Why is \(X[4]\) of a real 8-point signal always real, and what happens to \(\sin(\pi n)\)?

The bin-4 probe is \(e^{-j\pi n} = (-1)^n\), real, so \(X[4] = \sum x[n](-1)^n\) is real. \(\sin(\pi n) = 0\) at every integer n, so a Nyquist sine is invisible; the bin only records \(A\cos\varphi\).

Q4Using the radix-2 butterfly, compute the 2-point DFT of [a, b], then the 4-point DFT of [0, 1, 0, 0].

2-point: [a + b, a − b]. For [0, 1, 0, 0]: E = DFT[0, 0] = [0, 0], O = DFT[1, 0] = [1, 1]. \(X[0] = 0 + 1 = 1\), \(X[1] = 0 + (-j)(1) = -j\), \(X[2] = 0 - 1 = -1\), \(X[3] = 0 - (-j) = j\). That is \(W^k\): an impulse at n = 1 has unit magnitude and phase −90° per bin.

Q5\(f_s\) = 40 MHz, N = 1000. Which MATLAB index holds 5 MHz, and which holds −5 MHz?

\(\Delta f\) = 40 kHz, so 5 MHz is bin 125, index X(126). −5 MHz is bin N − 125 = 875, index X(876). After fftshift, DC sits at index 501 and 5 MHz at 626.

Q6A tone sits at 10.5 bins in a 64-point rectangular-window FFT. Roughly how much lower is the peak than for 10.0, and what is \(\lvert X[12]\rvert\) relative to \(\lvert X[10]\rvert\)?

Peak loss ≈ 3.92 dB (\(2/\pi\)). Using \(D_N(\delta) \approx N\,\mathrm{sinc}(\delta)\): bins 10 and 11 have \(\delta = \pm 0.5\), \(\lvert\mathrm{sinc}\rvert = 0.637\); bin 12 has \(\delta = -1.5\), \(\lvert\mathrm{sinc}\rvert = 0.212\). Ratio 0.212/0.637 = 1/3, i.e. −9.5 dB below the two peak bins.

Q7How many complex multiplies does a radix-2 FFT of N = 4096 need, compared with the direct DFT?

\((N/2)\log_2 N = 2048 \times 12 = 24\,576\), against \(N^2 = 16\,777\,216\): 683 times fewer.

02.16Figures from the literature

Published figures and tables worth opening next to this module. Links only (no copyrighted images are reproduced); each entry says what to look at and which equation of this module it confirms.

Cooley and Tukey (1965): the operation count, and the only measured data in the paper. The paper has no numbered figures. Look at the general count \(T = N(r_1 + r_2 + \dots + r_m)\) (Eq. 9) and its equal-factor form \(T(r) = rN\log_r N\) (Eq. 12), the remark that r = 3 is formally most efficient but only about 6 % better than r = 2 or 4, and the unlabeled table of measured IBM 7094 run times for 3D Fourier sums. Compare with (02.21): the radix-2 count \((N/2)\log_2 N\) multiplies is the modern refinement of \(2N\log_2 N\) operations.

J. W. Cooley and J. W. Tukey, "An algorithm for the machine calculation of complex Fourier series," Mathematics of Computation 19(90), 297–301 (1965). doi:10.1090/S0025-5718-1965-0178586-1

Harris (1978), Fig. 13 and Figs. 16–19: rectangle and cosα windows with their log-magnitude spectra. Fig. 13 is the rectangle (Dirichlet) window and its transform, the kernel \(D_N\) of (02.15) with its −13 dB first sidelobe and 6 dB/octave decay. Figs. 16–19 show the cosα family for α = 1 … 4; α = 2 is the Hann window of Demo 02.5, with −31.5 dB sidelobes and 18 dB/octave decay. Note Harris's "DFT-even" construction (Figs. 1–2): the periodic window used in this module's demos and scripts.

F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," Proceedings of the IEEE 66(1), 51–83 (1978). doi:10.1109/PROC.1978.10837

Harris (1978), Figs. 10–12 and Table I: leakage, resolution, scalloping. Fig. 10 is the transform of a finite-duration tone (leakage, as in Demo 02.5); Fig. 11 the resolution criterion for two equal lines (compare rule "Zero padding interpolates"); Fig. 12 plots highest sidelobe level against worst-case processing loss for every window. Table I lists ENBW, coherent gain and scalloping loss per window: check rectangle 3.92 dB and Hann 1.42 dB against (02.16), and Hann ENBW 1.50 bins.

Same paper as above; scalloping loss is defined in Section IV-D. doi:10.1109/PROC.1978.10837

Frigo and Johnson (2005), Figs. 1 and 2: measured FFT speed, power-of-two vs other sizes. Fig. 1 benchmarks double-precision 1D complex DFTs of power-of-two sizes on a 2.8 GHz Pentium IV; Fig. 2 repeats it for non-power-of-two sizes. Speed is reported in "mflops" normalised by \(5N\log_2 N\) per transform, so a flat curve means the run time scales as \(N\log N\). The fall-off at large N is the memory-traffic limit noted in the rule "N log N, not N²"; Fig. 2 is why padding to a power of two is a habit, not a requirement.

M. Frigo and S. G. Johnson, "The design and implementation of FFTW3," Proceedings of the IEEE 93(2), 216–231 (2005). doi:10.1109/JPROC.2004.840301 · open copy: fftw.org

02.17Go deeper

  • J. W. Cooley and J. W. Tukey, "An algorithm for the machine calculation of complex Fourier series," Mathematics of Computation 19, 297–301 (1965). The radix-2 FFT.
  • M. T. Heideman, D. H. Johnson and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine 1(4), 14–21 (1984).
  • F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," Proceedings of the IEEE 66(1), 51–83 (1978). Scalloping, ENBW and window tables.
  • A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Prentice Hall (2010). DFT, sampling the DTFT, FFT algorithms.
  • R. G. Lyons, Understanding Digital Signal Processing, 3rd ed., Prentice Hall (2010). The most hand-worked DFT and FFT examples in print.
  • M. Frigo and S. G. Johnson, "The design and implementation of FFTW3," Proceedings of the IEEE 93(2), 216–231 (2005). How production FFTs choose plans for any N.
  • Within this course: Module 01 (phasors), Module 04 (sampling, windows, k-space), Module 08 (I/Q), Module 27 (spectral estimation), Module 03 (convolution).