Worked examples — scipy.signal — filtering, convolution, FFT-based analysis

The scenario matrix

Before any example, let's enumerate every kind of situation this topic can throw at you. Each row is a cell. Every cell gets covered by at least one worked example below.

Notation we will use everywhere

Example 1 — Plain convolution (cell C1)

Step 1 — Output length. Why this step? Getting the length wrong means the whole array is misaligned. Full convolution length is $\text{len}(x)+\text{len}(h)-1=3+3-1=5$.

Step 2 — Flip $h$. Why? Convolution is defined with $h[n-k]$; the $-k$ means we read $h$ backwards. Flipping $[1,1,1]$ gives $[1,1,1]$ — symmetric, so it looks the same, but the concept is "flip then slide".

Step 3 — Slide and sum each overlap position. Why? Each output sample is the sum of products of samples that currently line up.

$y[0]=2\cdot 1=2$ $y[1]=2\cdot 1+0\cdot 1=2$ $y[2]=2\cdot 1+0\cdot 1+(-1)\cdot 1=1$ $y[3]=0\cdot 1+(-1)\cdot 1=-1$ $y[4]=(-1)\cdot 1=-1$

So $y=[2,2,1,-1,-1]$.

The figure shows the blue circles as the input $x[n]$ and the orange squares as the output $y=x\ast h$ shifted slightly to the right so you can compare them. Notice the output has 5 stems where the input had 3 — that visual "spreading" is exactly the $+\text{len}(h)-1$ extra length, and each orange value is labelled so you can trace it back to a step above.

Verify: A moving sum of a signal must have total (sum of all outputs) equal to $(\sum x)\cdot (\sum h)$. Here $\sum x=2+0-1=1$ and $\sum h=3$, so the product is $3$. Indeed $2+2+1-1-1=3$. ✓ In code: convolve([2,0,-1],[1,1,1]).

Example 2 — Degenerate kernel of length 1 (cell C2)

Step 1 — Output length. $3+1-1=3$. Why? The kernel occupies zero width beyond a single sample, so nothing spreads out.

Step 2 — There is only one overlap term at every position. Why? With $h=[3]$, the sum ${\sum }_{k}x[k]h[n-k]$ has exactly one non-zero term, $x[n]\cdot 3$.

$y=[4\cdot 3,7\cdot 3,-2\cdot 3]=[12,21,-6]$

Verify: convolve([4,7,-2],[3]) == [12,21,-6]. Also $\sum y=3\cdot \sum x=3\cdot 9=27$ and $12+21-6=27$. ✓

Example 3 — Kernel LONGER than the signal (cell C3)

Step 1 — Output length, from counting overlaps. Why this step? Instead of memorising the formula, count how many slide positions produce at least one overlapping pair. Start with $h$ just touching $x$ on the left (1 pair overlaps) and end when it just leaves on the right (1 pair overlaps). As you slide by one each time, the number of distinct positions is $\text{len}(x)+\text{len}(h)-1=2+4-1=5$. The "$-1$" is because the two end positions where a single sample overlaps are the boundaries you count once each. The formula is symmetric, so it doesn't matter which array is longer.

Step 2 — Use commutativity. Why this step? $x\ast h=h\ast x$ — convolution does not care which one you flip. So slide the shorter one ($[5,1]$) across the longer one; less bookkeeping.

$y[0]=5\cdot 1=5$ $y[1]=5\cdot 2+1\cdot 1=11$ $y[2]=5\cdot 3+1\cdot 2=17$ $y[3]=5\cdot 4+1\cdot 3=23$ $y[4]=1\cdot 4=4$

So $y=[5,11,17,23,4]$.

Verify: $\sum y=(\sum x)(\sum h)=6\cdot 10=60$, and $5+11+17+23+4=60$. ✓

Example 4 — Correlation, positive AND negative lag (cells C4, C5)

Step 1 — No flip. Why this step? Convolution flips the kernel (it models a system). Correlation asks "how similar is the template to each window?" — similarity needs the template kept in its true orientation, exactly as the definition above shows.

Step 2 — Compute the sliding dot products where $h$ fully overlaps $x$. The bump in $x$ sits at indices 2,3,4 as $[1,2,1]$ — an exact copy of $h$.

$\text{align }h\text{ at }x[2:5]:1\cdot 1+2\cdot 2+1\cdot 1=6$

That value $6$ is the maximum possible ($=\sum h^2=1+4+1$), a perfect match.

In the figure, the blue stems are the signal $x$ (a bump centred on index 3) and the orange curve is the full correlation. Watch where the orange curve reaches its peak (marked red): that peak sits over index 3, telling you the template best matches when centred there.

Step 3 — Map output index to lag. Why this step? scipy.signal.correlate(x, h, mode='full') returns an array of length $\text{len}(x)+\text{len}(h)-1=8$. Its index $i$ maps to a lag by $\text{lag}=i-(\text{len}(h)-1)=i-2.$ The peak sits at output index $i=4$, so $\text{lag}=4-2=2$. A positive lag of 2 means the template, whose left edge starts at 0, best fits when its window starts 2 samples in — i.e. the bump is centred at index 3 (cell C4).

Step 4 — Negative-lag case (C5). Now let $x_2=[1,2,1,0,0,0]$ — bump right at the start. Correlate $x_2$ with the same template but treat the template as a centred reference by using correlate(x2, h) and reading lag relative to the signal centre. The best overlap now occurs at output index $i=2$, giving $\text{lag}=2-2=0$ relative to the start, which is earlier than the signal's mid-point (index 2–3). Measured against that centre, the reported shift is negative: the pattern appears earlier than the reference. The sign of the lag literally tells you "template is ahead of (−) / behind (+) the reference".

Verify: peak value of correlate([0,0,1,2,1,0],[1,2,1]) equals $6$, equals $\sum h^2$, and its output index is 4 → lag $=4-2=2$. ✓

Example 5 — FFT of an on-bin pure tone (cell C6)

Step 1 — Bin spacing. Why this step? Each bin sits at $f_k=k{f}_s/N$ (from the DFT definition above). Here $\Delta f=f_s/N=8/8=1$ Hz. A 2 Hz tone lands exactly on bin $k=2$ — this is why it's called "on-bin". No leakage.

Step 2 — Where the $2|X[k]|/N$ rule comes from. Why this step? From the notation box: a real sine of amplitude $A$ deposits $|X[k]|=\frac{A}{2}N$ into bin $k$ (its energy is split with the mirror bin $N-k$). Invert that: $A=2|X[k]|/N$. Since $0<k=2<N/2=4$, the mirror partner exists, so the $\times 2$ is correct here. Predicted amplitude $=3$.

Left panel: the 8 time samples of the sine. Right panel: the one-sided (rfft) magnitude, already scaled by $2/N$. A single orange stem rises to exactly 3 at 2 Hz (marked red), and every other bin is flat at zero — the textbook picture of an on-bin tone.

Verify: the scaled one-sided magnitude at $k=2$ equals $3$ (to floating-point tolerance) and all other bins are ≈0. ✓

Example 6 — Frequency exactly at Nyquist (cell C7)

Step 1 — Locate the bin. $f_{Nyq}=f_s/2=4$ Hz $=$ bin $k=N/2=4$. Why care? The Nyquist bin has no distinct mirror partner — it is its own mirror ($N-k=8-4=4=k$).

Step 2 — Compute $X[4]$ directly from the DFT sum. Why this step? Rather than trust a half-remembered rule, plug into $X[k]={\sum }_{n}x[n]e^{-i2\pi kn/N}$ with $k=4$: $e^{-i2\pi \cdot 4\cdot n/8}=e^{-i\pi n}=(-1{)}^n,X[4]={\sum }_{n=0}^{7}2(-1{)}^n(-1{)}^n={\sum }_{n=0}^{7}2=16.$ So $|X[4]|=16=A\cdot N$ (here $A=2$, $N=8$). Note this is $A\cdot N$, not $A\cdot N/2$: at Nyquist all the energy piles into the single self-mirrored bin instead of being split into two.

Step 3 — Correct scaling at Nyquist. Why? Because there is no separate mirror bin, the right recovery is $A=|X[4]|/N=16/8=2$ ✓. The naive one-sided rule $2|X[4]|/N=2\cdot 16/8=4$ double-counts and reports twice the true amplitude.

Verify: for $x[n]=2(-1{)}^n$, $N=8$: $|X[4]|=16=A\cdot N$; correct amplitude $|X[4]|/N=2$, while the naive $2|X[4]|/N=4$ (wrong). ✓

Example 7 — Off-bin frequency: spectral leakage (cell C8)

Step 1 — Is it on a bin? $f_k=k\cdot 1$ Hz gives bins at 0,1,2,3,4 Hz. 2.5 Hz matches none. Why this matters: the DFT assumes the signal repeats every $N$ samples; a 2.5 Hz tone does not complete a whole number of cycles in 8 samples, so the periodic assumption is violated.

Step 2 — Energy leaks. Why? The tone's energy spreads into neighbouring bins (2 and 3 get the most, others get some). No single bin equals 1. This is spectral leakage.

Left panel: the 8 time samples — note they do not line up into a whole number of cycles. Right panel: instead of one clean stem you see a smeared hump; bins 2 and 3 (either side of 2.5 Hz) are the tallest but neither reaches the true amplitude 1. That smear is leakage.

Step 3 — The peak underestimates amplitude. The largest bin will read less than the true 1, because the energy is split. Windowing (Hann, etc.) reshapes but never fully removes leakage for off-bin tones.

Verify: the max scaled magnitude for the 2.5 Hz tone is strictly less than 1 (leakage), and the energy in bins 2 and 3 exceeds any other neighbouring pair. ✓

Example 8 — Resolving two close peaks (cell C9)

Step 1 — Current bin spacing. $\Delta f=f_s/N=1000/200=5$ Hz. Why this step? Two peaks only appear separate if they fall in different bins, i.e. their gap must exceed $\Delta f$. The gap here is $52-50=2$ Hz $<5$ Hz → they blur into one bump.

Step 2 — Faster sampling doesn't help. Why? Increasing $f_s$ raises the highest visible frequency (Nyquist), but $\Delta f=f_s/N$ — if you double $f_s$ while keeping the time fixed you also double $N$, so $\Delta f$ is unchanged. Resolution is set by record length in seconds $T=N/f_s$, since $\Delta f=1/T$.

Step 3 — Record longer. Why? $\Delta f=1/T$. To get $\Delta f\le 2$ Hz we need $T\ge 1/2=0.5$ s, i.e. $N\ge f_s\cdot T=1000\cdot 0.5=500$ samples. With $N=500$, $\Delta f=1000/500=2$ Hz — just enough; go a bit longer (say $N=1000$, $\Delta f=1$ Hz) for a clean split.

Verify: with $N=200$, $\Delta f=5>2$ (unresolved). With $N=500$, $\Delta f=2\le 2$ (just resolved). And required $T=1/2=0.5$ s → $N=500$. ✓

Example 9 — Real-world: remove 50 Hz mains hum (cell C10)

Step 1 — Choose cutoff. Why this step? We must pass the ECG (<40 Hz) and block 50 Hz. Put the cutoff between: $f_{cut}=45$ Hz is a safe choice.

Step 2 — Normalised cutoff (old API). Why? Butterworth's classic signature takes cutoff as a fraction of Nyquist. $W_n=\frac{f_{cut}}{f_s/2}=\frac{45}{500}=0.09$. (Or pass fs=1000 and 45 directly in modern scipy.)

Step 3 — Choose filtfilt, not lfilter, and write the code. Why filtfilt? A causal filter (lfilter) delays the signal — an ECG peak would move in time, ruining timing measurements. filtfilt runs the filter forward then backward; the two equal-and-opposite delays cancel → zero phase, so peaks stay put. Since we have the full recording offline, causality isn't needed. The full call:

from scipy.signal import butter, filtfilt
b, a = butter(4, 45, btype='low', fs=1000)   # or Wn=0.09 with the old API
clean = filtfilt(b, a, ecg)                   # zero-phase; peaks stay put

This connects to Digital Filter Design (why Butterworth) and Linear Time-Invariant Systems (why a filter is convolution).

Verify: $W_n=45/(1000/2)=0.09$, and $0<W_n<1$ so it is a legal digital cutoff. ✓

Example 10 — Exam twist: the illegal cutoff (cell C11)

Step 1 — What 80 means without fs=. Why this step? In the old signature, the second argument is a fraction of Nyquist, i.e. a number in $(0,1)$. So 80 is read as "80× Nyquist" — nonsense; scipy rejects any $W_n>1$.

Step 2 — Fix A: normalise. $W_n=\frac{80}{1000/2}=\frac{80}{500}=0.16$. Then butter(4, 0.16, btype='low') is valid because $0<0.16<1$.

Step 3 — Fix B: pass real Hz. Supply fs and let scipy normalise internally:

from scipy.signal import butter
b, a = butter(4, 80, btype='low', fs=1000)   # Hz allowed once fs is given
# equivalently: butter(4, 0.16, btype='low')  # old-style normalised cutoff

Verify: $80/(1000/2)=0.16$, which lies in $(0,1)$ — legal. The raw 80 does not. ✓

Connections

• Fourier Transform — leakage and on-bin/off-bin theory behind Ex 5–7.
• Aliasing and Sampling Theorem — the Nyquist edge case of Ex 6.
• Linear Time-Invariant Systems — why convolution (Ex 1–3) is filtering.
• Cross-correlation and Template Matching — the lag reading of Ex 4.
• Digital Filter Design — Butterworth choice in Ex 9–10.
• numpy.fft — the FFT engine under Ex 5–8.