Exercises — NumPy FFT — np.fft module

• np.fft.fft(x) → the forward DFT: a complex array of length $N$; each entry $X_k$ is the complex amplitude (DFT coefficient) of frequency-bin $k$ (including negative-frequency bins).
• np.fft.rfft(x) → forward DFT of a real input, returning only the non-redundant bins (for even $N$ that is $N/2+1$ bins, 0 up to Nyquist).
• np.fft.fftfreq(N, d) → the frequency in Hz of each fft bin, including negative frequencies, with sample spacing $d=\Delta t$.
• np.fft.rfftfreq(N, d) → the frequencies for rfft (only $0$ through Nyquist, all $\ge 0$).
• np.fft.ifft(X) → the inverse DFT, taking a spectrum back to the time domain.

$X_k={\sum }_{n=0}^{N-1}{x}_n{e}^{-2\pi ikn/N},k=0,\dots ,N-1$ The exponent carries a minus sign (it is the conjugate of the probe wave we correlate against), and $k$ runs from $0$ to $N-1$ — exactly $N$ output bins for $N$ input samples.

$T=N/f_s=800/400===2\text{s}==$. $\Delta f=f_s/N=400/800===0.5\text{Hz}==$ (equivalently $1/T=1/2=0.5$ Hz). The two formulas agree, as they must.

First get the rate: $f_s=1/\Delta t=1/0.002=500$ Hz. Then $f_k=k{f}_s/N=30\times 500/500===30\text{Hz}==$. (Here $N/f_s=1$ s, so bin index numerically equals frequency — a lucky coincidence of $T=1$ s, not a general rule.)

Where does the $2/N$ come from — the full "why". Write the real sine using Euler's formula: $A\sin (2\pi f_0{t}_n)=\frac{A}{2i}\left(e^{+2\pi i{f}_0{t}_n}-e^{-2\pi i{f}_0{t}_n}\right).$ So a real sine is literally two complex probe-waves: one at $+f_0$ with weight $A/(2i)$ and one at $-f_0$ with weight $-A/(2i)$ — each carries half the amplitude, $A/2$. Now recall the orthogonality fact (proved in L5.1): when you correlate against the matching probe, every mismatched term cancels and the matching term is summed $N$ times. That multiplies whichever weight sits at that frequency by $N$. Hence the bin at $+f_0$ reads $|X_{\text{peak}}|=\frac{A}{2}\times N=\frac{AN}{2}.$ To invert this and get $A$ back, multiply by $2/N$: $A=|X_{\text{peak}}|\times \frac{2}{N}=5000\times \frac{2}{2000}===5==.✓$ The DC bin is different: a constant offset is one term, not a $\pm f$ pair, so it is scaled by $N$ (not $N/2$) and you divide by just $1/N$.

$X_0={\sum }_n{x}_n=1+1+1+1===4==$. For $k\ne 0$: $X_k={\sum }_{n=0}^{3}1\cdot e^{-2\pi ikn/4}={\sum }_{n=0}^3{r}^n,r=e^{-2\pi ik/4}.$ This is a geometric series with $r\ne 1$ and $r^4=e^{-2\pi ik}=1$, so $\sum =\frac{1-r^4}{1-r}=\frac{0}{1-r}=0$. Result: $[4,0,0,0]$. A constant signal has only zero-frequency content — no oscillation exists to detect.

Figure s01 (below): the left panel plots the flat time signal $[1,1,1,1]$; the right panel plots its magnitude spectrum $|X|=[4,0,0,0]$. The point of the picture is that all the energy collapses into the single DC bin ($k=0$) and every oscillating probe returns exactly zero — a flat line "rings" only for the zero-frequency probe.

The Nyquist frequency (see Sampling and Nyquist Theorem) is half the sampling rate: $f_{\text{Nyq}}=f_s/2=8000/2===4000\text{Hz}==.$ rfft returns bins $0,\dots ,N/2$ (for even $N$), so the Nyquist frequency sits in the last bin, index $N/2===4000==$. Any real tone above 4000 Hz would alias (fold back) and appear at the wrong frequency.

Start from the definition and substitute $N-k$ for $k$: $X_{N-k}={\sum }_{n=0}^{N-1}{x}_n{e}^{-2\pi i(N-k)n/N}={\sum }_n{x}_n{e}^{-2\pi in}{e}^{+2\pi ikn/N}.$ Since $e^{-2\pi in}=1$ for every integer $n$, this is $X_{N-k}={\sum }_n{x}_n{e}^{+2\pi ikn/N}.$ Now $x_n$ is real so $x_n=\overline{x_n}$, and the whole sum is the conjugate of $X_k$ (taking the conjugate flips the sign of the exponent, $\overline{e^{-i\theta }}=e^{+i\theta }$): ${\sum }_n{x}_n{e}^{+2\pi ikn/N}=\overline{{\sum }_n{x}_n{e}^{-2\pi ikn/N}}=\overline{X_k}.■$ Consequence: the top half of the spectrum is the mirror-conjugate of the bottom half, so rfft keeps only the non-redundant half with no loss.

Resolution is $\Delta f=f_s/N=1/T$.

• Alice: $\Delta f=1000/1000=1\text{Hz}$ (window $T=1$ s).
• Bob: $\Delta f=1000/4000=0.25\text{Hz}$ (window $T=4$ s). Bob's resolution is finer by a factor of $1/0.25===4==$. Longer recording ⇒ tighter bins. Note the sampling rate is identical, so both see the same maximum frequency (Nyquist 500 Hz) — a longer record buys resolution, not range.

Figure s02 (below): both spectra are plotted near 40 Hz. The red dashed line marks the true 40 Hz peak — notice both curves peak at the same place (same range), but Bob's blue curve (longer window $T=4$ s) has bins spaced 4× closer together than Alice's yellow curve, so it can separate tones that Alice would merge. This is the visual meaning of "resolution = $1/T$."

The "mystery" peak is the negative-frequency twin of your tone. fft lays the bins out as $[0,1,\dots ,N/2,\underset{\text{negative freqs}}{\underbrace{-N/2+1,\dots ,-1}}],$ so the top of the array holds negative frequencies masquerading as high indices. Fix: call np.fft.fftshift(X) (and fftshift(fftfreq(...))) to move zero frequency to the centre, giving a spectrum symmetric about 0. Then the two peaks sit at $+f$ and $-f$, mirror images, exactly as Euler's formula predicts for a real cosine.

import numpy as np
fs = 500.0
N  = 1000                 # 2 s * 500 Hz
t  = np.arange(N) / fs
x  = 2*np.sin(2*np.pi*60*t) + 1
 
X   = np.fft.rfft(x)
f   = np.fft.rfftfreq(N, d=1/fs)
amp = np.abs(X) * 2 / N            # sinusoid scaling
 
k60 = np.argmin(np.abs(f - 60))    # bin nearest 60 Hz
print(f[k60], amp[k60])            # ~60.0 Hz, ~2.0
print(np.abs(X[0]) / N)            # DC = 1.0

Predictions:

• The tone bin recovers amplitude $===2.0==$ (the $2/N$ factor undoes the $N/2$ scaling of a real sine).
• DC bin: $|X_0|/N=(\sum x_n)/N=\text{mean}(x)$. The sine averages to 0 over whole cycles, so the mean is the offset $===1.0==$. Note DC uses $1/N$, not $2/N$.

The gap most people leave is implementing the special case. Here is the pattern:

import numpy as np
X   = np.fft.rfft(x)              # real input, length N
N   = x.size
amp = np.abs(X) * 2 / N           # start: assume every bin is a sine pair
amp[0] = np.abs(X[0]) / N         # DC bin has NO mirror twin -> use 1/N
if N % 2 == 0:                    # even N: last rfft bin is exactly Nyquist
    amp[-1] = np.abs(X[-1]) / N   # Nyquist bin also has no twin -> use 1/N

Why: the $2/N$ factor assumes the energy was split into a $+f$ and $-f$ mirror pair. Bin 0 (DC) is a single constant term, and — only when $N$ is even — the last bin lands exactly on Nyquist and is likewise unpaired. Both must be divided by $N$, not $2N$. For odd $N$ the rfft output does not include an exact Nyquist bin, so only the DC correction is needed (the if N % 2 == 0 guard handles this automatically). Run it on L4.1's signal and the DC entry now reads exactly $1.0$ while the 60 Hz entry reads $2.0$.

b = [1,0,0,0] is the identity for circular convolution (a single spike at position 0). Convolving anything with it returns that thing unchanged, so we expect $[1,2,3,4]$. Via the theorem:

a = np.array([1,2,3,4], float)
b = np.array([1,0,0,0], float)
c = np.fft.ifft(np.fft.fft(a) * np.fft.fft(b)).real
print(np.round(c, 6))   # [1. 2. 3. 4.]

Why it works: $\text{fft}(b)=[1,1,1,1]$ (all bins 1, since a unit spike contains every frequency equally), so multiplying leaves $\text{fft}(a)$ untouched, and ifft returns $a$. Result: $==[1,2,3,4]==$.

The artefact is spectral leakage — energy smears into neighbouring bins instead of sitting in one clean spike. Cause: the DFT assumes the window repeats periodically; a non-integer number of cycles creates a discontinuity at the wrap-around, and a sharp jump contains many frequencies. Fix: apply a window function (Hann, Hamming, ...) that tapers the signal smoothly to zero at both ends, removing the discontinuity. See Spectral Leakage and Windowing. Higher-level windows and spectral tools live in scipy.signal.

Let $r=e^{2\pi i(k-m)/N}$.

• Case $k=m$: then $r=e^0=1$, so each of the $N$ terms is $1$ and the sum is $==N==$.
• Case $k\ne m$: $r\ne 1$ but $r^N=e^{2\pi i(k-m)}=1$. The geometric series gives ${\sum }_{n=0}^{N-1}{r}^n=\frac{1-r^N}{1-r}=\frac{1-1}{1-r}===0==.$ Why it matters: when you plug $X_k$ into the inverse formula $x_n=\frac{1}{N}{\sum }_k{X}_k{e}^{+2\pi ikn/N}$, the orthogonality makes every "wrong" frequency cancel and leaves exactly the original sample — so the round trip is exact. This is the algebraic heart of Discrete Fourier Transform.

Use $X_k={\sum }_{n=0}^3{x}_n{e}^{-2\pi ikn/4}$ with $x=[0,1,0,-1]$. Only $n=1$ and $n=3$ contribute.

• $X_0=0+1+0+(-1)===0==$ (DC / mean is zero — a pure oscillation).
• $X_1=1\cdot e^{-i\pi /2}+(-1)\cdot e^{-i3\pi /2}=(-i)+(-1)(i)=-i-i===-2i==$.
• $X_2=1\cdot e^{-i\pi }+(-1)\cdot e^{-i3\pi }=(-1)+(-1)(-1)=-1+1===0==$.
• $X_3=1\cdot e^{-i3\pi /2}+(-1)\cdot e^{-i9\pi /2}=(i)+(-1)(-i)=i+i===2i==$.

Result: $[0,-2i,0,2i]$. Interpretation: energy only in bins 1 and 3, which are the $+$/$-$ frequency pair (note $X_3=\overline{X_1}$, i.e. the conjugate — confirming conjugate symmetry for real input). The purely imaginary values are the signature of a sine (odd function); a cosine would give purely real bins.

Left side (time energy): $1^2+2^2+3^2+4^2=1+4+9+16===30==$. Right side (frequency energy): first fft([1,2,3,4]):

• $X_0=1+2+3+4=10$.
• $X_1=1+2(-i)+3(-1)+4(i)=(1-3)+(-2+4)i=-2+2i$, so $|X_1{|}^2=4+4=8$.
• $X_2=1-2+3-4=-2$, so $|X_2{|}^2=4$.
• $X_3=\overline{X_1}=-2-2i$, so $|X_3{|}^2=8$.

Sum $=100+8+4+8=120$. Divide by $N=4$: $120/4===30==$. Both sides equal 30 ✓. Meaning: the FFT just repackages energy from time into frequency — none is created or destroyed. The $1/N$ is bookkeeping from NumPy's convention.