Spot the Bug

Ten broken DSP snippets. Can you find the bugs?

Ten DSP code snippets. Each has exactly one bug. Can you find them all?

Difficulty increases as you go. Solutions are hidden — try before you peek.

Rating	Meaning
★	Warm-up
★★	Intermediate
★★★	Tricky
★★★★	Devious
★★★★★	Evil

Keep a tally as you go. Scoring at the bottom.

Bug 1 (★) — dB conversion

You measure signal power = 0.5 W and noise power = 0.005 W. You compute SNR in decibels:

import numpy as np

signal_power = 0.5
noise_power = 0.005
snr_db = 20 * np.log10(signal_power / noise_power)
print(f"SNR = {snr_db:.1f} dB")  # prints 40.0 dB

That SNR looks suspiciously high. What went wrong?

Reveal the bug

The bug: 20 * np.log10 is for amplitude (voltage) ratios. For power ratios, use 10 * np.log10.

The relationship is: power is proportional to amplitude squared, so $10 \log_{10}(P) = 10 \log_{10}(A^2) = 20 \log_{10}(A)$. When your quantities are already powers, the factor is 10, not 20.

Fixed code:

snr_db = 10 * np.log10(signal_power / noise_power)
print(f"SNR = {snr_db:.1f} dB")  # prints 20.0 dB  ✓

Rule of thumb: Ask yourself — is this a power ratio or an amplitude ratio? Power $\to$ 10 * log10. Amplitude $\to$ 20 * log10.

Bug 2 (★) — Aliasing check

You want to check whether a 5000 Hz tone will alias at your sample rate:

fs = 8000  # sample rate
max_freq = fs  # Nyquist frequency
tone = 5000

if tone <= max_freq:
    print("Safe, no aliasing")
else:
    print("Warning: will alias!")

The check says “Safe, no aliasing.” But when you listen to the sampled signal, the tone sounds completely wrong. It has the wrong pitch. Why?

Reveal the bug

The bug: The Nyquist frequency is $f_s / 2$, not $f_s$.

The sampling theorem says you can represent frequencies up to half the sample rate. At $f_s = 8000$ Hz, the Nyquist frequency is 4000 Hz. A 3500 Hz tone is fine — but the check would also pass a 5000 Hz tone, which will alias.

Fixed code:

max_freq = fs / 2  # Nyquist frequency

This is one of the most common beginner mistakes in DSP. The sample rate and the Nyquist frequency are not the same thing.

Bug 3 (★★) — IIR filter coefficients

You want to implement a simple first-order recursive filter: $y[n] = x[n] + 0.9\,y[n-1]$.

from scipy.signal import lfilter

# y[n] = x[n] + 0.9*y[n-1]
x = np.random.randn(1000)
y = lfilter([1], [1, 0.9], x)

The output decays where it should grow, and grows where it should decay. What happened?

Reveal the bug

The bug: lfilter uses the convention $a_0 y[n] + a_1 y[n-1] + \ldots = b_0 x[n] + b_1 x[n-1] + \ldots$, which means it negates the feedback coefficients internally.

Your difference equation $y[n] = x[n] + 0.9\,y[n-1]$ rearranges to $y[n] - 0.9\,y[n-1] = x[n]$, so the a vector should be [1, -0.9].

Fixed code:

y = lfilter([1], [1, -0.9], x)

This is a classic SciPy trap. The sign in the a coefficients is the opposite of what you’d write in the difference equation. See the gotcha log for more on coefficient conventions.

Bug 4 (★★) — FFT frequency axis

You compute an FFT and build a frequency axis for plotting:

fs = 1000
x = np.sin(2 * np.pi * 50 * np.arange(1024) / fs)

X = np.fft.rfft(x)
freqs = np.arange(len(X)) * fs / len(X)

import matplotlib.pyplot as plt
plt.plot(freqs, np.abs(X))
plt.xlabel("Frequency (Hz)")
plt.show()

The peak shows up at the wrong frequency. Where’s the bug?

Reveal the bug

The bug: The frequency spacing is fs / len(x), not fs / len(X). The rfft output has len(x)//2 + 1 bins, which is a different length than x. Using len(X) in the denominator stretches the frequency axis.

For len(x) = 1024, len(X) = 513. The frequency resolution is $f_s / N = 1000/1024 \approx 0.977$ Hz, but the buggy code uses $1000/513 \approx 1.95$ Hz — almost double. Every frequency reads nearly twice its true value.

Fixed code:

freqs = np.fft.rfftfreq(len(x), d=1/fs)

Or manually: np.arange(len(X)) * fs / len(x). The key is that the DFT length is determined by the input, not the output.

Bug 5 (★★★) — Windowed FFT amplitude

You apply a Hanning window and try to recover the true amplitude of a sine wave:

fs = 1000
N = 1024
t = np.arange(N) / fs
x = 3.0 * np.sin(2 * np.pi * 50 * t)  # amplitude = 3.0

window = np.hanning(N)
X = np.fft.rfft(x * window)
amplitude = 2 * np.abs(X).max() / N  # factor 2 for one-sided
print(f"Recovered amplitude: {amplitude:.2f}")  # should be 3.0

The recovered amplitude is about 1.5 instead of 3.0. What went wrong?

Reveal the bug

The bug: When you apply a window, the normalization factor is the sum of the window, not $N$. The Hanning window has a mean of 0.5, so np.sum(window) $\approx N/2$. Dividing by $N$ instead of np.sum(window) cuts the amplitude in half.

Fixed code:

amplitude = 2 * np.abs(X).max() / np.sum(window)
print(f"Recovered amplitude: {amplitude:.2f}")  # ≈ 3.0  ✓

Why? The window attenuates the signal. The total attenuation equals the sum of the window values. To undo it, divide by that sum. For a rectangular window (no windowing), np.sum(window) = N, so the two formulas agree — which is why you only notice this bug when you start windowing.

Bug 6 (★★★) — One-sided PSD

You estimate the power spectral density of white noise and compare to the theoretical value:

fs = 1000
N = 4096
x = np.random.randn(N)  # unit variance white noise

X = np.fft.rfft(x)
psd = np.abs(X)**2 / (N * fs)

# Theoretical one-sided PSD of unit-variance white noise: 2/fs
import matplotlib.pyplot as plt
plt.semilogy(np.fft.rfftfreq(N, 1/fs), psd)
plt.axhline(2/fs, color='r', ls='--', label='theoretical')
plt.legend()
plt.show()

The estimated PSD sits about 3 dB below the theoretical line. Why?

Reveal the bug

The bug: When converting from a two-sided to a one-sided spectrum, you must double the power in the interior bins (all bins except DC and Nyquist). The negative-frequency energy folds into the positive side.

Fixed code:

psd = np.abs(X)**2 / (N * fs)
psd[1:-1] *= 2  # double interior bins for one-sided spectrum

DC (index 0) and Nyquist (index -1) have no mirror image, so they stay as-is. Every other bin represents two symmetric frequencies, so their power doubles. Missing this step loses exactly half the power — which is 3 dB.

Bug 7 (★★★★) — Unstable notch filter

You design a notch filter to remove 50 Hz mains hum from a signal sampled at 1000 Hz:

from scipy.signal import lfilter

fs = 1000
f0 = 50
w0 = 2 * np.pi * f0 / fs

# Notch filter: zeros on unit circle, poles just inside
r = 1.01  # "close to 1 for a sharp notch"
b = [1, -2*np.cos(w0), 1]
a = [1, -2*r*np.cos(w0), r**2]

# Apply to a test signal
t = np.arange(2000) / fs
x = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*200*t)
y = lfilter(b, a, x)

The output explodes. Can you spot why before running it?

Reveal the bug

The bug: r = 1.01 puts the poles outside the unit circle. The pole radius is $r = 1.01 > 1$ (the poles sit at $r\,e^{\pm j\omega_0}$, so $r$ itself is the radius; $r^2 = 1.0201$ is the $a_2$ coefficient), so the filter is unstable. The output grows exponentially.

The poles must be inside the unit circle for a stable causal filter. The pole radius $r$ controls the notch bandwidth: closer to 1 means sharper, but it must be strictly less than 1.

Fixed code:

r = 0.99  # poles INSIDE the unit circle

Intuition: the zeros on the unit circle create the null at 50 Hz. The poles pull the response back up everywhere else. If the poles are outside the circle, they don’t just pull — they blow up.

Bug 8 (★★★★) — Real-time zero-phase filtering

A colleague writes a real-time audio processing function and proudly announces it does “zero-phase filtering”:

from scipy.signal import sosfiltfilt, butter

sos = butter(4, 300, btype='low', fs=16000, output='sos')

def process_block(block, sos, zi):
    """Process one block of real-time audio."""
    filtered = sosfiltfilt(sos, block)  # zero-phase!
    return filtered

The filter works in offline tests but produces glitchy audio when used in a real-time stream. Why?

Reveal the bug

The bug: sosfiltfilt is a non-causal filter. It works by filtering forward, then backward over the entire signal. In a real-time system, you don’t have “the entire signal” — you only have the current block.

Calling sosfiltfilt on each block independently means:

Each block is treated as a separate signal (no state continuity between blocks).
The forward-backward pass creates startup transients at every block boundary.
This is not true zero-phase filtering — it’s a broken approximation.

Fixed code:

from scipy.signal import sosfilt

def process_block(block, sos, zi):
    """Process one block of real-time audio."""
    filtered, zi = sosfilt(sos, block, zi=zi)  # causal, with state
    return filtered, zi

Zero-phase filtering is fundamentally an offline operation. In real-time, you use causal filtering (sosfilt) with state variables (zi) that carry over between blocks. If you need linear phase, use a symmetric FIR filter (which adds latency but is at least causal).

Bug 9 (★★★★★) — Chirp signal generation

You generate a linear chirp sweeping from 1 kHz to 5 kHz over 10 ms:

fs = 44100
f0, f1, T = 1000, 5000, 0.01
t = np.arange(int(T * fs)) / fs

# Instantaneous frequency increases linearly
freq = f0 + (f1 - f0) * t / T
chirp = np.sin(2 * np.pi * freq * t)

A spectrogram reveals the chirp sweeps too fast — it reaches 5 kHz well before $t = T$. The code looks right. What’s the subtle error?

Reveal the bug

The bug: The instantaneous frequency of a sinusoid is the derivative of its phase, not a multiplier of $t$. Writing sin(2*pi*freq*t) where freq varies with t does not produce a signal whose instantaneous frequency equals freq.

The phase of a linear chirp is the integral of the instantaneous frequency:

\[\phi(t) = 2\pi \int_0^t f(\tau)\,d\tau = 2\pi\left(f_0 t + \frac{f_1 - f_0}{2T} t^2\right)\]

The product freq * t gives $f_0 t + (f_1-f_0)t^2/T$ — missing the factor of $1/2$ on the quadratic term. The actual sweep rate is double what you intended.

Fixed code:

phase = 2 * np.pi * (f0 * t + 0.5 * (f1 - f0) * t**2 / T)
chirp = np.sin(phase)

Or use SciPy, which handles this correctly:

from scipy.signal import chirp as scipy_chirp
chirp = scipy_chirp(t, f0, T, f1, method='linear')

This is a classic trap. It catches experienced engineers too — see the gotcha log.

Bug 10 (★★★★★) — Matched filter implementation

You implement a matched filter to detect a known chirp template in a noisy signal:

# Template: a short chirp pulse
from scipy.signal import chirp
template = chirp(t_template, 1000, t_template[-1], 5000)

# Received signal: noisy, with the chirp buried somewhere
received = noise + np.roll(template_padded, delay)

# Matched filter: correlate received with template
output = np.convolve(received, template, mode='full')
peak = np.argmax(np.abs(output))
print(f"Detected at sample {peak}")

The detected delay is consistently wrong — it’s off by the template length. The SNR gain is also lower than expected for an asymmetric template. Why?

Reveal the bug

The bug: np.convolve computes convolution, not correlation. Convolution flips (time-reverses) one of the inputs internally. For a matched filter, you need cross-correlation, which uses the template as-is.

For symmetric templates (like a rectangular pulse), convolution and correlation give the same result, so you’d never notice. But a chirp is asymmetric — time-reversing it changes the signal. The convolution output is the correlation with a reversed chirp, which doesn’t match the target.

Fixed code:

# Option 1: time-reverse the template so convolution undoes the flip
output = np.convolve(received, template[::-1], mode='full')

# Option 2: use correlate directly
output = np.correlate(received, template, mode='full')

The math: Convolution is $(f * g)[n] = \sum_k f[k]\,g[n-k]$, which flips $g$. Correlation is $(f \star g)[n] = \sum_k f[k]\,g[n+k]$, no flip. Pre-reversing the template cancels the flip in convolution.

How did you do?

Count the bugs you found before revealing the answer.

Score	Verdict
10/10	DSP ninja. You’ve debugged this code before.
7–9	Solid — you catch the common traps.
4–6	Getting there — review the gotcha log!
0–3	Time to re-read the basics chapters. No shame — these are real bugs from real code.

Many of these bugs are documented in the gotcha log, which tracks common DSP pitfalls we’ve encountered. If you enjoyed this, the exercise sets have more debugging challenges embedded in the problem sets.

--- title: "Spot the Bug" subtitle: "Ten broken DSP snippets. Can you find the bugs?" --- Ten DSP code snippets. Each has **exactly one bug**. Can you find them all? Difficulty increases as you go. Solutions are hidden --- try before you peek. | Rating | Meaning | |--------|---------| | ★ | Warm-up | | ★★ | Intermediate | | ★★★ | Tricky | | ★★★★ | Devious | | ★★★★★ | Evil | Keep a tally as you go. Scoring at the bottom. --- ::: {.callout-warning title="Bug 1 (★) --- dB conversion"} You measure signal power = 0.5 W and noise power = 0.005 W. You compute SNR in decibels: ```python import numpy as np signal_power = 0.5 noise_power = 0.005 snr_db = 20 * np.log10(signal_power / noise_power) print(f"SNR = {snr_db:.1f} dB") # prints 40.0 dB ``` That SNR looks suspiciously high. What went wrong? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** `20 * np.log10` is for **amplitude** (voltage) ratios. For **power** ratios, use `10 * np.log10`. The relationship is: power is proportional to amplitude squared, so $10 \log_{10}(P) = 10 \log_{10}(A^2) = 20 \log_{10}(A)$. When your quantities are already powers, the factor is 10, not 20. **Fixed code:** ```python snr_db = 10 * np.log10(signal_power / noise_power) print(f"SNR = {snr_db:.1f} dB") # prints 20.0 dB ✓ ``` **Rule of thumb:** Ask yourself --- is this a power ratio or an amplitude ratio? Power $\to$ `10 * log10`. Amplitude $\to$ `20 * log10`. ::: ::: ::: {.callout-warning title="Bug 2 (★) --- Aliasing check"} You want to check whether a 5000 Hz tone will alias at your sample rate: ```python fs = 8000 # sample rate max_freq = fs # Nyquist frequency tone = 5000 if tone <= max_freq: print("Safe, no aliasing") else: print("Warning: will alias!") ``` The check says "Safe, no aliasing." But when you listen to the sampled signal, the tone sounds completely wrong. It has the wrong pitch. Why? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** The Nyquist frequency is $f_s / 2$, not $f_s$. The sampling theorem says you can represent frequencies up to **half** the sample rate. At $f_s = 8000$ Hz, the Nyquist frequency is 4000 Hz. A 3500 Hz tone is fine --- but the check would also pass a 5000 Hz tone, which **will** alias. **Fixed code:** ```python max_freq = fs / 2 # Nyquist frequency ``` This is one of the most common beginner mistakes in DSP. The sample rate and the Nyquist frequency are not the same thing. ::: ::: ::: {.callout-warning title="Bug 3 (★★) --- IIR filter coefficients"} You want to implement a simple first-order recursive filter: $y[n] = x[n] + 0.9\,y[n-1]$. ```python from scipy.signal import lfilter # y[n] = x[n] + 0.9*y[n-1] x = np.random.randn(1000) y = lfilter([1], [1, 0.9], x) ``` The output decays where it should grow, and grows where it should decay. What happened? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** `lfilter` uses the convention $a_0 y[n] + a_1 y[n-1] + \ldots = b_0 x[n] + b_1 x[n-1] + \ldots$, which means it **negates** the feedback coefficients internally. Your difference equation $y[n] = x[n] + 0.9\,y[n-1]$ rearranges to $y[n] - 0.9\,y[n-1] = x[n]$, so the `a` vector should be `[1, -0.9]`. **Fixed code:** ```python y = lfilter([1], [1, -0.9], x) ``` This is a classic SciPy trap. The sign in the `a` coefficients is the opposite of what you'd write in the difference equation. See the [gotcha log](https://github.com/ducroq/dsp-workshop/blob/main/memory/gotcha-log.md) for more on coefficient conventions. ::: ::: ::: {.callout-warning title="Bug 4 (★★) --- FFT frequency axis"} You compute an FFT and build a frequency axis for plotting: ```python fs = 1000 x = np.sin(2 * np.pi * 50 * np.arange(1024) / fs) X = np.fft.rfft(x) freqs = np.arange(len(X)) * fs / len(X) import matplotlib.pyplot as plt plt.plot(freqs, np.abs(X)) plt.xlabel("Frequency (Hz)") plt.show() ``` The peak shows up at the wrong frequency. Where's the bug? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** The frequency spacing is `fs / len(x)`, not `fs / len(X)`. The `rfft` output has `len(x)//2 + 1` bins, which is a different length than `x`. Using `len(X)` in the denominator stretches the frequency axis. For `len(x) = 1024`, `len(X) = 513`. The frequency resolution is $f_s / N = 1000/1024 \approx 0.977$ Hz, but the buggy code uses $1000/513 \approx 1.95$ Hz --- almost double. Every frequency reads nearly twice its true value. **Fixed code:** ```python freqs = np.fft.rfftfreq(len(x), d=1/fs) ``` Or manually: `np.arange(len(X)) * fs / len(x)`. The key is that the DFT length is determined by the **input**, not the output. ::: ::: ::: {.callout-warning title="Bug 5 (★★★) --- Windowed FFT amplitude"} You apply a Hanning window and try to recover the true amplitude of a sine wave: ```python fs = 1000 N = 1024 t = np.arange(N) / fs x = 3.0 * np.sin(2 * np.pi * 50 * t) # amplitude = 3.0 window = np.hanning(N) X = np.fft.rfft(x * window) amplitude = 2 * np.abs(X).max() / N # factor 2 for one-sided print(f"Recovered amplitude: {amplitude:.2f}") # should be 3.0 ``` The recovered amplitude is about 1.5 instead of 3.0. What went wrong? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** When you apply a window, the normalization factor is the **sum of the window**, not $N$. The Hanning window has a mean of 0.5, so `np.sum(window)` $\approx N/2$. Dividing by $N$ instead of `np.sum(window)` cuts the amplitude in half. **Fixed code:** ```python amplitude = 2 * np.abs(X).max() / np.sum(window) print(f"Recovered amplitude: {amplitude:.2f}") # ≈ 3.0 ✓ ``` **Why?** The window attenuates the signal. The total attenuation equals the sum of the window values. To undo it, divide by that sum. For a rectangular window (no windowing), `np.sum(window) = N`, so the two formulas agree --- which is why you only notice this bug when you start windowing. ::: ::: ::: {.callout-warning title="Bug 6 (★★★) --- One-sided PSD"} You estimate the power spectral density of white noise and compare to the theoretical value: ```python fs = 1000 N = 4096 x = np.random.randn(N) # unit variance white noise X = np.fft.rfft(x) psd = np.abs(X)**2 / (N * fs) # Theoretical one-sided PSD of unit-variance white noise: 2/fs import matplotlib.pyplot as plt plt.semilogy(np.fft.rfftfreq(N, 1/fs), psd) plt.axhline(2/fs, color='r', ls='--', label='theoretical') plt.legend() plt.show() ``` The estimated PSD sits about 3 dB below the theoretical line. Why? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** When converting from a two-sided to a one-sided spectrum, you must **double the power in the interior bins** (all bins except DC and Nyquist). The negative-frequency energy folds into the positive side. **Fixed code:** ```python psd = np.abs(X)**2 / (N * fs) psd[1:-1] *= 2 # double interior bins for one-sided spectrum ``` DC (index 0) and Nyquist (index -1) have no mirror image, so they stay as-is. Every other bin represents two symmetric frequencies, so their power doubles. Missing this step loses exactly half the power --- which is 3 dB. ::: ::: ::: {.callout-warning title="Bug 7 (★★★★) --- Unstable notch filter"} You design a notch filter to remove 50 Hz mains hum from a signal sampled at 1000 Hz: ```python from scipy.signal import lfilter fs = 1000 f0 = 50 w0 = 2 * np.pi * f0 / fs # Notch filter: zeros on unit circle, poles just inside r = 1.01 # "close to 1 for a sharp notch" b = [1, -2*np.cos(w0), 1] a = [1, -2*r*np.cos(w0), r**2] # Apply to a test signal t = np.arange(2000) / fs x = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*200*t) y = lfilter(b, a, x) ``` The output explodes. Can you spot why before running it? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** `r = 1.01` puts the poles **outside** the unit circle. The pole radius is $r = 1.01 > 1$ (the poles sit at $r\,e^{\pm j\omega_0}$, so $r$ itself is the radius; $r^2 = 1.0201$ is the $a_2$ coefficient), so the filter is **unstable**. The output grows exponentially. The poles must be **inside** the unit circle for a stable causal filter. The pole radius $r$ controls the notch bandwidth: closer to 1 means sharper, but it must be strictly less than 1. **Fixed code:** ```python r = 0.99 # poles INSIDE the unit circle ``` **Intuition:** the zeros on the unit circle create the null at 50 Hz. The poles pull the response back up everywhere else. If the poles are outside the circle, they don't just pull --- they blow up. ::: ::: ::: {.callout-warning title="Bug 8 (★★★★) --- Real-time zero-phase filtering"} A colleague writes a real-time audio processing function and proudly announces it does "zero-phase filtering": ```python from scipy.signal import sosfiltfilt, butter sos = butter(4, 300, btype='low', fs=16000, output='sos') def process_block(block, sos, zi): """Process one block of real-time audio.""" filtered = sosfiltfilt(sos, block) # zero-phase! return filtered ``` The filter works in offline tests but produces glitchy audio when used in a real-time stream. Why? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** `sosfiltfilt` is a **non-causal** filter. It works by filtering forward, then backward over the **entire** signal. In a real-time system, you don't have "the entire signal" --- you only have the current block. Calling `sosfiltfilt` on each block independently means: 1. Each block is treated as a separate signal (no state continuity between blocks). 2. The forward-backward pass creates startup transients at every block boundary. 3. This is not true zero-phase filtering --- it's a broken approximation. **Fixed code:** ```python from scipy.signal import sosfilt def process_block(block, sos, zi): """Process one block of real-time audio.""" filtered, zi = sosfilt(sos, block, zi=zi) # causal, with state return filtered, zi ``` Zero-phase filtering is fundamentally an **offline** operation. In real-time, you use causal filtering (`sosfilt`) with state variables (`zi`) that carry over between blocks. If you need linear phase, use a symmetric FIR filter (which adds latency but is at least causal). ::: ::: ::: {.callout-warning title="Bug 9 (★★★★★) --- Chirp signal generation"} You generate a linear chirp sweeping from 1 kHz to 5 kHz over 10 ms: ```python fs = 44100 f0, f1, T = 1000, 5000, 0.01 t = np.arange(int(T * fs)) / fs # Instantaneous frequency increases linearly freq = f0 + (f1 - f0) * t / T chirp = np.sin(2 * np.pi * freq * t) ``` A spectrogram reveals the chirp sweeps too fast --- it reaches 5 kHz well before $t = T$. The code looks right. What's the subtle error? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** The instantaneous frequency of a sinusoid is the **derivative of its phase**, not a multiplier of $t$. Writing `sin(2*pi*freq*t)` where `freq` varies with `t` does not produce a signal whose instantaneous frequency equals `freq`. The phase of a linear chirp is the **integral** of the instantaneous frequency: $$\phi(t) = 2\pi \int_0^t f(\tau)\,d\tau = 2\pi\left(f_0 t + \frac{f_1 - f_0}{2T} t^2\right)$$ The product `freq * t` gives $f_0 t + (f_1-f_0)t^2/T$ --- missing the factor of $1/2$ on the quadratic term. The actual sweep rate is **double** what you intended. **Fixed code:** ```python phase = 2 * np.pi * (f0 * t + 0.5 * (f1 - f0) * t**2 / T) chirp = np.sin(phase) ``` Or use SciPy, which handles this correctly: ```python from scipy.signal import chirp as scipy_chirp chirp = scipy_chirp(t, f0, T, f1, method='linear') ``` This is a classic trap. It catches experienced engineers too --- see the [gotcha log](https://github.com/ducroq/dsp-workshop/blob/main/memory/gotcha-log.md). ::: ::: ::: {.callout-warning title="Bug 10 (★★★★★) --- Matched filter implementation"} You implement a matched filter to detect a known chirp template in a noisy signal: ```python # Template: a short chirp pulse from scipy.signal import chirp template = chirp(t_template, 1000, t_template[-1], 5000) # Received signal: noisy, with the chirp buried somewhere received = noise + np.roll(template_padded, delay) # Matched filter: correlate received with template output = np.convolve(received, template, mode='full') peak = np.argmax(np.abs(output)) print(f"Detected at sample {peak}") ``` The detected delay is consistently wrong --- it's off by the template length. The SNR gain is also lower than expected for an asymmetric template. Why? ::: {.callout-note collapse="true" title="Reveal the bug"} **The bug:** `np.convolve` computes **convolution**, not **correlation**. Convolution flips (time-reverses) one of the inputs internally. For a matched filter, you need cross-correlation, which uses the template **as-is**. For symmetric templates (like a rectangular pulse), convolution and correlation give the same result, so you'd never notice. But a chirp is asymmetric --- time-reversing it changes the signal. The convolution output is the correlation with a **reversed** chirp, which doesn't match the target. **Fixed code:** ```python # Option 1: time-reverse the template so convolution undoes the flip output = np.convolve(received, template[::-1], mode='full') # Option 2: use correlate directly output = np.correlate(received, template, mode='full') ``` **The math:** Convolution is $(f * g)[n] = \sum_k f[k]\,g[n-k]$, which flips $g$. Correlation is $(f \star g)[n] = \sum_k f[k]\,g[n+k]$, no flip. Pre-reversing the template cancels the flip in convolution. ::: ::: --- ## How did you do? Count the bugs you found **before** revealing the answer. | Score | Verdict | |-------|---------| | **10/10** | DSP ninja. You've debugged this code before. | | **7--9** | Solid --- you catch the common traps. | | **4--6** | Getting there --- review the [gotcha log](https://github.com/ducroq/dsp-workshop/blob/main/memory/gotcha-log.md)! | | **0--3** | Time to re-read the [basics chapters](../basics/index.qmd). No shame --- these are real bugs from real code. | --- *Many of these bugs are documented in the [gotcha log](https://github.com/ducroq/dsp-workshop/blob/main/memory/gotcha-log.md), which tracks common DSP pitfalls we've encountered. If you enjoyed this, the [exercise sets](index.qmd) have more debugging challenges embedded in the problem sets.*