Exercises: Filter Design

Practice problems for Chapter 6

These exercises accompany Chapter 6: Filter Design.

Key formulas for this chapter

Formula	Description
FIR group delay: \((N-1)/2\) samples	Linear-phase FIR delay
Bilinear pre-warping: \(\omega_a = \frac{2}{T}\tan(\omega_d T/2)\)	Analog-to-digital frequency mapping
FIR order (Harris rule): \(N \approx A_s \cdot f_s / (22 \cdot \Delta f)\)	Order estimation from specs
Always use SOS form for IIR filters of order \(\geq 4\)	Numerical stability rule
butter(N, fc, fs=fs, output='sos')	Scipy handles pre-warping internally

Basic

Exercise 1: Filter specification vocabulary

A lowpass filter specification states: \(f_p = 500\) Hz, \(f_s = 700\) Hz, \(\delta_p = 0.5\) dB, \(A_s = 50\) dB, at \(f_\text{samp} = 4000\) Hz.

1. What is the transition bandwidth in Hz?

2. What is the normalised transition width \(\Delta f / f_\text{nyq}\)?

3. If you increase the stopband attenuation to 80 dB while keeping all else fixed, what happens to the required filter order?

Solution

1. Transition bandwidth = \(f_s - f_p = 700 - 500 = 200\) Hz.

2. \(f_\text{nyq} = 2000\) Hz. Normalised width = \(200/2000 = 0.1\).

3. The required order increases. The design triangle tells us: with transition width and ripple fixed, demanding more attenuation forces a higher order.

Exercise 2: Window method basics

You design a 31-tap FIR lowpass filter using firwin with a Hamming window at \(f_s = 2000\) Hz and cutoff 300 Hz.

1. What is the group delay of this filter in samples? In milliseconds?

2. Is the impulse response symmetric? What does that imply about the phase?

3. If you switch to a Blackman window (same \(N\)), what changes: the cutoff frequency, the transition width, or the stopband attenuation?

Solution

1. Group delay = \((N-1)/2 = 15\) samples = \(15/2000 = 7.5\) ms.

2. Yes, firwin produces symmetric (linear-phase) impulse responses by default. This means all frequencies experience the same delay: no phase distortion.

3. The cutoff stays at 300 Hz. The Blackman window has a wider main lobe than Hamming, so the transition width increases. But its side lobes are lower (−58 dB vs −43 dB), so stopband attenuation improves.

Exercise 3: IIR filter families

Match each IIR family to its characteristic:

Family	Characteristic
Butterworth	???
Chebyshev I	???
Chebyshev II	???
Elliptic	???

Options: (A) Equiripple in both bands, (B) Maximally flat everywhere, (C) Equiripple passband + monotonic stopband, (D) Flat passband + equiripple stopband.

Solution

Family	Characteristic
Butterworth	(B) Maximally flat everywhere
Chebyshev I	(C) Equiripple passband + monotonic stopband
Chebyshev II	(D) Flat passband + equiripple stopband
Elliptic	(A) Equiripple in both bands

Elliptic achieves the sharpest transition for a given order by allowing ripple in both bands.

Exercise 4: SOS vs direct form

A 10th-order Butterworth filter is designed with butter(10, 200, fs=1000).

1. How many biquad sections does the SOS form have?

2. How many delay elements does each biquad need (DF II)?

3. Why might the direct-form (b, a) implementation become unstable even though all poles are inside the unit circle?

Solution

1. \(\lceil 10/2 \rceil = 5\) sections.

2. 2 delay elements per biquad (DF II is the canonical form with minimum delays).

3. The 11-element a polynomial has coefficients that must represent 10 roots with high precision. Floating-point rounding in the 64-bit coefficients can effectively shift poles outside the unit circle. With SOS, each section only has 2 poles, so coefficient errors stay small and localised.

Intermediate

Exercise 5: Design and compare FIR methods

Design FIR lowpass filters for the same specification (\(f_p = 400\) Hz, \(f_s = 600\) Hz, \(f_\text{samp} = 4000\) Hz) using:

1. Window method: firwin(N, cutoff, fs=fs) with a Hamming window.

2. Parks–McClellan: remez(N, bands, desired, fs=fs).

Use the same order \(N = 61\) for both. Plot both magnitude responses and compare the stopband attenuation.

Solution

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import firwin, remez, freqz

fs = 4000
N = 61

h_win = firwin(N, 500, fs=fs, window='hamming')  # cutoff at midpoint of transition band
h_pm = remez(N, [0, 400, 600, fs/2], [1, 0], fs=fs)

fig, ax = plt.subplots(figsize=(10, 4))
for h, label in [(h_win, 'Window (Hamming)'), (h_pm, 'Parks–McClellan')]:
    w, H = freqz(h, [1], worN=2048, fs=fs)
    ax.plot(w, 20*np.log10(np.maximum(np.abs(H), 1e-10)), linewidth=1.5, label=label)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_ylim(-80, 5)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_title(f'FIR comparison (N={N})')
fig.tight_layout()
plt.show()

The Parks–McClellan filter has equiripple behaviour: the stopband ripple is evenly distributed, so the worst-case attenuation is better than the window method for the same order. The window method has deeper nulls in some places but higher peaks in others.

Exercise 6: Bilinear transform by hand

A first-order analog lowpass filter has transfer function \(H_a(s) = 1 / (1 + s/\omega_c)\) with \(\omega_c = 2\pi \times 500\) rad/s.

1. Apply the bilinear transform \(s = \frac{2}{T}\frac{z-1}{z+1}\) with \(f_s = 4000\) Hz to derive \(H_d(z)\).

2. What are the digital filter coefficients \((b_0, b_1, a_0, a_1)\)?

3. At what digital frequency does the analog cutoff \(\omega_c\) actually appear? (Hint: apply the warping formula.)

Solution

1. With \(T = 1/4000\) and \(\omega_c = 2\pi \times 500 = 3141.6\) rad/s:

\[H_d(z) = \frac{1}{1 + \frac{1}{\omega_c}\frac{2}{T}\frac{z-1}{z+1}} = \frac{1 + z^{-1}}{(1 + \frac{2}{\omega_c T}) + (1 - \frac{2}{\omega_c T})z^{-1}}\]

Compute \(\frac{2}{\omega_c T} = \frac{2 \times 4000}{3141.6} = 2.5465\):

\[H_d(z) = \frac{1 + z^{-1}}{3.5465 + (-1.5465)z^{-1}} = \frac{0.2820(1 + z^{-1})}{1 - 0.4361\,z^{-1}}\]

2. \(b_0 = b_1 = 0.2820\), \(a_0 = 1\), \(a_1 = -0.4361\).

3. The warping formula gives \(\omega_d = \frac{2}{T}\arctan\frac{\omega_c T}{2} = 8000 \times \arctan(0.393) = 8000 \times 0.3745 = 2996\) rad/s, or \(f_d = 476.8\) Hz. The cutoff is slightly lower than the analog 500 Hz due to frequency warping.

import numpy as np
from scipy.signal import freqz

fs = 4000
wc = 2 * np.pi * 500
T = 1/fs
alpha = 2 / (wc * T)

b0 = 1 / (1 + alpha)
a1 = (1 - alpha) / (1 + alpha)
b = [b0, b0]
a = [1, a1]

print(f"b = [{b[0]:.4f}, {b[1]:.4f}]")
print(f"a = [1, {a[1]:.4f}]")

# Verify -3 dB point
w, H = freqz(b, a, worN=2048, fs=fs)
mag_dB = 20 * np.log10(np.abs(H))
idx_3dB = np.argmin(np.abs(mag_dB - (-3)))
print(f"−3 dB frequency: {w[idx_3dB]:.1f} Hz (expected ~476 Hz)")

b = [0.2820, 0.2820]
a = [1, -0.4361]
−3 dB frequency: 475.6 Hz (expected ~476 Hz)

Exercise 7: Minimum order design

You need a lowpass filter with: \(f_p = 1\) kHz, \(f_s = 1.2\) kHz, \(\delta_p = 1\) dB, \(A_s = 60\) dB, \(f_\text{samp} = 8\) kHz.

1. Use buttord, cheb1ord, and ellipord to find the minimum order for each family.

2. Design all three filters (SOS form) and plot their magnitude responses on the same axes.

3. Which filter would you choose for a battery-powered embedded device? Why?

Solution

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import buttord, cheb1ord, ellipord, butter, cheby1, ellip, sosfreqz

fs = 8000
fp, fst = 1000, 1200
rp, As = 1, 60

Nb, Wb = buttord(fp, fst, rp, As, fs=fs)
Nc, Wc = cheb1ord(fp, fst, rp, As, fs=fs)
Ne, We = ellipord(fp, fst, rp, As, fs=fs)

print(f"Butterworth:  N = {Nb}")
print(f"Chebyshev I:  N = {Nc}")
print(f"Elliptic:     N = {Ne}")

sos_b = butter(Nb, Wb, fs=fs, output='sos')
sos_c = cheby1(Nc, rp, Wc, fs=fs, output='sos')
sos_e = ellip(Ne, rp, As, We, fs=fs, output='sos')

fig, ax = plt.subplots(figsize=(10, 4))
for sos, name in [(sos_b, f'Butterworth (N={Nb})'),
                   (sos_c, f'Chebyshev I (N={Nc})'),
                   (sos_e, f'Elliptic (N={Ne})')]:
    w, H = sosfreqz(sos, worN=4096, fs=fs)
    ax.plot(w, 20*np.log10(np.maximum(np.abs(H), 1e-10)), linewidth=1.5, label=name)

ax.axhline(-rp, color='gray', linewidth=0.5, linestyle='--')
ax.axhline(-As, color='gray', linewidth=0.5, linestyle='--')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_xlim(0, 3000)
ax.set_ylim(-80, 5)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

Butterworth:  N = 37
Chebyshev I:  N = 13
Elliptic:     N = 7

3. Elliptic: it needs the fewest biquad sections, meaning fewer multiply-accumulate operations per sample and lower power consumption. The cost is ripple in both bands, which is usually acceptable for sensor data.

Challenge

Exercise 8: Notch filter for power-line interference

A biomedical signal sampled at \(f_s = 500\) Hz is contaminated by 50 Hz power-line interference.

1. Design a second-order IIR notch filter by placing zeros on the unit circle at \(\theta = 2\pi \times 50/500\) and poles at the same angle with radius \(R = 0.95\). Write out the transfer function \(H(z)\).

2. Implement the filter with lfilter and verify the frequency response.

3. Generate a test signal (a 10 Hz ECG-like pulse train + 50 Hz interference) and demonstrate the filter’s effect.

Solution

1. \(\theta = 2\pi \times 50/500 = \pi/5\) rad.

Zeros at \(z = e^{\pm j\pi/5}\): numerator \((1 - e^{j\theta}z^{-1})(1 - e^{-j\theta}z^{-1}) = 1 - 2\cos\theta\,z^{-1} + z^{-2}\).

Poles at \(z = 0.95\,e^{\pm j\pi/5}\): denominator \(1 - 2R\cos\theta\,z^{-1} + R^2 z^{-2}\).

\[H(z) = \frac{1 - 2\cos(\pi/5)\,z^{-1} + z^{-2}}{1 - 1.9\cos(\pi/5)\,z^{-1} + 0.9025\,z^{-2}}\]

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz, lfilter

fs = 500
f_notch = 50
theta = 2 * np.pi * f_notch / fs
R = 0.95

b = [1, -2*np.cos(theta), 1]
a = [1, -2*R*np.cos(theta), R**2]

# b) Frequency response
w, H = freqz(b, a, worN=4096, fs=fs)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(w, 20*np.log10(np.maximum(np.abs(H), 1e-10)), 'C0-', linewidth=1.5)
axes[0].axvline(f_notch, color='C3', linewidth=0.5, linestyle='--')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_ylim(-50, 5)
axes[0].set_title('Notch filter response')
axes[0].grid(True, alpha=0.3)

# c) Test signal
rng = np.random.default_rng(42)
t = np.arange(5 * fs) / fs
ecg_pulse = np.zeros_like(t)
for beat_time in np.arange(0, 5, 0.1):  # 10 Hz pulse train
    idx = int(beat_time * fs)
    if idx < len(ecg_pulse) - 10:
        ecg_pulse[idx:idx+10] += np.hanning(10)

interference = 0.5 * np.sin(2 * np.pi * 50 * t)
x = ecg_pulse + interference + 0.05 * rng.standard_normal(len(t))
y = lfilter(b, a, x)

t_plot = slice(0, 500)  # first second
axes[1].plot(t[t_plot], x[t_plot], 'b-', linewidth=0.5, alpha=0.5, label='Noisy')
axes[1].plot(t[t_plot], y[t_plot], 'C0-', linewidth=1, label='Filtered')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('50 Hz removal')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)

fig.tight_layout()
plt.show()

Exercise 9: Audio equalizer with biquads

Design a 3-band audio equalizer using cascaded biquad sections. The sampling rate is \(f_s = 44100\) Hz.

• Bass boost: +6 dB peaking filter centred at 100 Hz, \(Q = 1\).
• Mid cut: −4 dB peaking filter centred at 1000 Hz, \(Q = 2\).
• Treble shelf: +3 dB high-shelf above 8000 Hz.

1. Use the Audio EQ Cookbook peaking EQ formula to design the biquad coefficients for the bass and mid bands.

2. Cascade the three sections and plot the combined frequency response.

3. Generate a white noise signal, filter it, and plot the before/after PSD to verify the equalizer curve.

Solution

The peaking EQ biquad from Robert Bristow-Johnson’s Audio EQ Cookbook:

\[A = 10^{G_\text{dB}/40}, \quad \alpha = \frac{\sin\omega_0}{2Q}\]

\[b_0 = 1 + \alpha A, \quad b_1 = -2\cos\omega_0, \quad b_2 = 1 - \alpha A\] \[a_0 = 1 + \alpha/A, \quad a_1 = -2\cos\omega_0, \quad a_2 = 1 - \alpha/A\]

This produces unity gain at DC and Nyquist, with a boost or cut of \(G_\text{dB}\) at \(\omega_0\).

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz, lfilter, welch

fs = 44100

def peaking_eq(f0, gain_dB, Q, fs):
    """Audio EQ Cookbook peaking EQ biquad."""
    A = 10**(gain_dB / 40)
    w0 = 2 * np.pi * f0 / fs
    alpha = np.sin(w0) / (2 * Q)
    b = np.array([1 + alpha * A, -2 * np.cos(w0), 1 - alpha * A])
    a = np.array([1 + alpha / A, -2 * np.cos(w0), 1 - alpha / A])
    return b / a[0], a / a[0]  # normalise so a[0] = 1

def high_shelf(f0, gain_dB, Q, fs):
    """Audio EQ Cookbook high shelf biquad."""
    A = 10**(gain_dB / 40)
    w0 = 2 * np.pi * f0 / fs
    alpha = np.sin(w0) / (2 * Q)
    cos_w0 = np.cos(w0)
    sqrtA = np.sqrt(A)
    b = np.array([
        A * ((A + 1) + (A - 1)*cos_w0 + 2*sqrtA*alpha),
        -2*A * ((A - 1) + (A + 1)*cos_w0),
        A * ((A + 1) + (A - 1)*cos_w0 - 2*sqrtA*alpha)
    ])
    a = np.array([
        (A + 1) - (A - 1)*cos_w0 + 2*sqrtA*alpha,
        2 * ((A - 1) - (A + 1)*cos_w0),
        (A + 1) - (A - 1)*cos_w0 - 2*sqrtA*alpha
    ])
    return b / a[0], a / a[0]

# a) Design each band
b_bass, a_bass = peaking_eq(100, +6, 1.0, fs)
b_mid, a_mid = peaking_eq(1000, -4, 2.0, fs)
b_treble, a_treble = high_shelf(8000, +3, 0.707, fs)

# Combined frequency response
w_freq = np.linspace(0, np.pi, 4096)
H_total = np.ones(len(w_freq), dtype=complex)
for b, a in [(b_bass, a_bass), (b_mid, a_mid), (b_treble, a_treble)]:
    _, H = freqz(b, a, worN=w_freq)
    H_total *= H

f_hz = w_freq * fs / (2 * np.pi)

fig, axes = plt.subplots(2, 1, figsize=(10, 7))

# b) Combined response
axes[0].semilogx(f_hz[1:], 20*np.log10(np.maximum(np.abs(H_total[1:]), 1e-10)),
                  'C0-', linewidth=2)
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Gain [dB]')
axes[0].set_xlim(20, fs/2)
axes[0].set_ylim(-10, 10)
axes[0].set_title('3-band equalizer response')
axes[0].grid(True, alpha=0.3, which='both')
axes[0].axhline(0, color='k', linewidth=0.5)

# c) Filter white noise and compare PSD
rng = np.random.default_rng(42)
N = 5 * fs
x = rng.standard_normal(N)

y = x.copy()
for b, a in [(b_bass, a_bass), (b_mid, a_mid), (b_treble, a_treble)]:
    y = lfilter(b, a, y)

f_x, psd_x = welch(x, fs, nperseg=4096)
f_y, psd_y = welch(y, fs, nperseg=4096)

axes[1].semilogx(f_x[1:], 10*np.log10(np.maximum(psd_y[1:], 1e-20)
                                        / np.maximum(psd_x[1:], 1e-20)),
                  'C0-', linewidth=1.5, label='Measured EQ curve')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Gain [dB]')
axes[1].set_xlim(20, fs/2)
axes[1].set_ylim(-10, 10)
axes[1].set_title('Measured equalizer effect (output PSD / input PSD)')
axes[1].grid(True, alpha=0.3, which='both')
axes[1].axhline(0, color='k', linewidth=0.5)

fig.tight_layout()
plt.show()

The PSD ratio confirms the designed EQ curve: +6 dB bass boost around 100 Hz, −4 dB dip near 1 kHz, and +3 dB treble lift above 8 kHz.

Exercise 10: Debug: wrong frequency normalisation (Basic)

A student wants a 4th-order Butterworth lowpass filter with a −3 dB cutoff at 200 Hz, sampled at 1000 Hz. Here is their code:

from scipy.signal import butter, sosfilt
import numpy as np

fs = 1000
fc = 200  # Hz

sos = butter(4, fc, btype='low', output='sos')

1. Run the code. What cutoff frequency does this actually produce? Why?

2. Fix the single error so the filter has its −3 dB point at 200 Hz.

3. How would the answer change if fs = 8000 Hz?

Solution

1. butter with no fs argument expects the cutoff as a normalised frequency in the range \([0, 1]\), where 1.0 maps to the Nyquist frequency. Passing fc = 200 is interpreted as \(200 \times f_\text{nyq}\), which is far outside the valid range. SciPy will raise a ValueError: Digital filter critical frequencies must be 0 < Wn < 1.

2. Supply the sampling rate so SciPy handles the normalisation:

from scipy.signal import butter, sosfilt, sosfreqz
import numpy as np
import matplotlib.pyplot as plt

fs = 1000
fc = 200

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

w, H = sosfreqz(sos, worN=2048, fs=fs)
mag_dB = 20 * np.log10(np.maximum(np.abs(H), 1e-10))
idx = np.argmin(np.abs(mag_dB - (-3)))
print(f"−3 dB point: {w[idx]:.1f} Hz  (target: {fc} Hz)")

−3 dB point: 200.0 Hz  (target: 200 Hz)

3. With fs = 8000 Hz the correct call becomes butter(4, 200, btype='low', fs=8000, output='sos'). The normalised frequency changes to \(200/4000 = 0.05\), so the same physical cutoff now sits much closer to DC on the normalised scale; the filter is “sharper” relative to Nyquist but still correct in Hz.

Exercise 11: Choosing a filter family (Basic)

For each scenario below, name the most appropriate IIR filter family (Butterworth, Chebyshev I, Chebyshev II, or Elliptic) and give a one-sentence justification.

1. A medical-grade ECG monitor where no passband distortion is allowed but you can tolerate slight stopband ripple.

2. A consumer audio application where you want the sharpest possible roll-off for a given filter order and can accept ripple in both bands.

3. A control loop anti-aliasing filter where the passband must be as flat as possible and the stopband rejection only needs to meet a minimum spec.

4. A software-defined radio front-end where strict stopband rejection is needed but any passband ripple is unacceptable.

Solution

1. Chebyshev II: the passband is maximally flat (monotonic), and the equiripple behaviour is confined entirely to the stopband.

2. Elliptic: equiripple in both passband and stopband gives the steepest possible transition for a given order, so the cutoff between bands is sharpest.

3. Butterworth: maximally flat magnitude response in both bands; no ripple anywhere, which simplifies characterisation of the control-loop gain margin.

4. Chebyshev II: the stopband meets a hard equiripple floor while the passband remains ripple-free. (Elliptic achieves a sharper roll-off but sacrifices passband flatness, which is disallowed here.)

Exercise 12: Debug: lfilter vs sosfilt for high-order filters (Basic)

A student designs a 12th-order Chebyshev I lowpass filter and notices the output looks wrong: saturated or all zeros. Here is the code:

from scipy.signal import cheby1, lfilter
import numpy as np

fs = 1000
rng = np.random.default_rng(0)
x = rng.standard_normal(5000)

b, a = cheby1(12, 0.5, 200, fs=fs)
y = lfilter(b, a, x)
print(y[:5])

1. What is likely causing the incorrect output?

2. Fix the code using the numerically stable approach.

3. What is the rule of thumb for when to switch from lfilter(b, a, ...) to sosfilt?

Solution

1. High-order filters in direct-form (b, a) representation suffer from coefficient quantisation errors. A 12th-order polynomial has 13 coefficients that must collectively represent 12 poles. Small floating-point rounding errors shift pole locations, often placing poles outside the unit circle. The filter becomes unstable and the output diverges (or underflows to zero after overflow).

2. Design and filter using second-order sections:

from scipy.signal import cheby1, sosfilt, lfilter
import numpy as np

fs = 1000
rng = np.random.default_rng(0)
x = rng.standard_normal(5000)

# Broken: direct form
b, a = cheby1(12, 0.5, 200, fs=fs)
y_bad = lfilter(b, a, x)

# Fixed: SOS form
sos = cheby1(12, 0.5, 200, fs=fs, output='sos')
y_good = sosfilt(sos, x)

print(f"Direct form output — first 5 values: {y_bad[:5]}")
print(f"SOS output        — first 5 values:  {y_good[:5]}")
print(f"Direct form has NaNs or Infs: {not np.all(np.isfinite(y_bad))}")
print(f"SOS output is finite:          {np.all(np.isfinite(y_good))}")

Direct form output — first 5 values: [4.88494588e-07 8.44580017e-06 7.17938764e-05 4.03723456e-04
 1.70045105e-03]
SOS output        — first 5 values:  [4.88494588e-07 8.44580017e-06 7.17938764e-05 4.03723456e-04
 1.70045105e-03]
Direct form has NaNs or Infs: False
SOS output is finite:          True

3. The rule of thumb: use SOS whenever the filter order exceeds about 4–6. Each biquad section has only two poles, so coefficient errors remain small and localised rather than accumulating across a high-degree polynomial.

Intermediate (continued)

Exercise 13: Data-driven: clean a noisy vibration signal (Intermediate)

A vibration sensor on a rotating machine runs at \(f_s = 2000\) Hz. The useful content is below 200 Hz; everything above 500 Hz is broadband noise. Design and apply an appropriate filter.

1. Synthesise a test signal: a 50 Hz sinusoid (amplitude 1.0) plus white noise at 0 dB SNR.

2. Design a filter that meets: passband edge 200 Hz (≤ 1 dB ripple), stopband edge 500 Hz (≥ 40 dB attenuation). Use the minimum-order IIR family that achieves this, chosen for flat passband.

3. Apply the filter with zero-phase processing (sosfiltfilt) and plot the input PSD and output PSD on the same axes.

4. Estimate the SNR of the 50 Hz tone before and after filtering.

Solution

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import cheb2ord, cheby2, sosfilt, sosfiltfilt, sosfreqz, welch

rng = np.random.default_rng(42)
fs = 2000
N = 4 * fs  # 4 seconds

t = np.arange(N) / fs
signal = np.sin(2 * np.pi * 50 * t)           # 50 Hz tone, power = 0.5
noise  = rng.standard_normal(N) / np.sqrt(2)   # noise power = 0.5 → 0 dB SNR
x = signal + noise

# b) Minimum order — flat passband → Chebyshev II (flat passband, equiripple stopband)
order, Wn = cheb2ord(200, 500, 1, 40, fs=fs)
print(f"Chebyshev II order: {order}")
sos = cheby2(order, 40, Wn, fs=fs, output='sos')

# c) Zero-phase filtering
y = sosfiltfilt(sos, x)

# PSD
f_x, Pxx = welch(x, fs, nperseg=1024)
f_y, Pyy = welch(y, fs, nperseg=1024)

fig, ax = plt.subplots(figsize=(10, 4))
ax.semilogy(f_x, Pxx, 'C7-', linewidth=1, alpha=0.6, label='Input')
ax.semilogy(f_y, Pyy, 'C0-', linewidth=1.5, label='Filtered')
ax.axvline(200, color='C3', linewidth=0.8, linestyle='--', label='Passband edge')
ax.axvline(500, color='C1', linewidth=0.8, linestyle='--', label='Stopband edge')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('PSD [V²/Hz]')
ax.set_title('Vibration signal: before and after filtering')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, fs/2)
fig.tight_layout()
plt.show()

# d) SNR estimate: signal power at 50 Hz bin vs. noise floor
# Approximate SNR: ratio of peak PSD bin to median noise floor (density ratio, not integrated power)
def tone_snr(signal_clean, x_noisy, fs, f_tone, nperseg=1024):
    f, Pxx = welch(x_noisy, fs, nperseg=nperseg)
    idx = np.argmin(np.abs(f - f_tone))
    signal_power = Pxx[idx]
    noise_floor  = np.median(Pxx)
    return 10 * np.log10(signal_power / noise_floor)

snr_before = tone_snr(signal, x,  fs, 50)
snr_after  = tone_snr(signal, y,  fs, 50)
print(f"SNR before filtering: {snr_before:.1f} dB")
print(f"SNR after  filtering: {snr_after:.1f} dB")

Chebyshev II order: 4

SNR before filtering: 24.6 dB
SNR after  filtering: 105.1 dB

The zero-phase filter suppresses the broadband noise above 500 Hz without introducing phase distortion; the 50 Hz tone is recovered with clearly improved SNR. Note that sosfiltfilt doubles the effective order (applies the filter twice), so a lower design order suffices for the same attenuation spec.

Exercise 14: FIR vs IIR: latency trade-off (Intermediate)

You are filtering a real-time audio stream at \(f_s = 48000\) Hz for live performance monitoring. The specification requires a lowpass filter at 4 kHz (−3 dB) with 1 dB passband ripple, stopband ≥ 60 dB above 6 kHz.

1. Estimate the minimum FIR order using the Harris rule of thumb \(N \approx A_s / (22 \Delta f / f_s)\) where \(\Delta f = f_\text{stop} - f_\text{pass}\) and \(A_s\) is stopband attenuation in dB. What is the one-way latency in milliseconds?

2. Find the minimum-order elliptic IIR filter with the same specification. What is the filter order and how many biquad sections does it need?

3. For a live concert system the round-trip latency budget is 10 ms total. Allowing 4 ms for the filter, can you use FIR? Can you use IIR? Justify your choice.

Solution

1. Harris rule of thumb: \(N \approx A_s \cdot f_s / (22 \cdot \Delta f)\).

\[\Delta f = 6000 - 4000 = 2000 \text{ Hz}, \quad A_s = 60 \text{ dB}\]

\[N \approx \frac{60 \times 48000}{22 \times 2000} = \frac{2{,}880{,}000}{44{,}000} \approx 65\]

Group delay = \((N-1)/2 \approx 32\) samples. Latency = \(32/48000 \approx 0.67\) ms one-way.

import numpy as np
from scipy.signal import ellipord, ellip

fs = 48000
fp, fst = 4000, 6000
rp, As = 1.0, 60

# a) Harris approximation
delta_f = fst - fp
N_fir = int(np.ceil(As * fs / (22 * delta_f)))
latency_fir_ms = ((N_fir - 1) / 2) / fs * 1000
print(f"FIR order (Harris): ~{N_fir}, one-way latency: {latency_fir_ms:.2f} ms")

# b) Minimum-order elliptic IIR
order, Wn = ellipord(fp, fst, rp, As, fs=fs)
n_biquads = int(np.ceil(order / 2))
print(f"Elliptic IIR order: {order}, biquad sections: {n_biquads}")
print(f"IIR one-way latency: effectively 0 (causal, real-time)")

FIR order (Harris): ~66, one-way latency: 0.68 ms
Elliptic IIR order: 6, biquad sections: 3
IIR one-way latency: effectively 0 (causal, real-time)

2. The elliptic filter requires far fewer sections and introduces only one sample of computational latency (the processing delay of the algorithm), not a structural group delay.

3. At ~0.67 ms one-way the FIR fits within the 4 ms budget, and importantly it is linear-phase, so all frequencies in the passband arrive simultaneously. The IIR introduces frequency-dependent group delay (phase distortion) but essentially zero structural latency. For live performance monitoring the IIR is preferred because it minimises perceptible latency, and small phase differences across 0–4 kHz are generally inaudible in a monitoring context. If the task were audio recording where phase accuracy matters, the FIR (or zero-phase IIR offline) would be better.

Exercise 15: Debug: forgotten frequency pre-warping (Intermediate)

A student designs a digital Butterworth lowpass filter with a cutoff at exactly 3 kHz (\(f_s = 8000\) Hz) by directly converting an analog prototype. The code below produces the wrong digital cutoff. Identify the bug and fix it.

from scipy.signal import butter, freqz
import numpy as np

fs = 8000
fc_digital = 3000  # Hz — desired digital cutoff

# Student converts Hz to rad/s for analog design
wc_analog = 2 * np.pi * fc_digital   # ← student's pre-warp step

# Design analog prototype and bilinear transform manually via butter
sos = butter(4, wc_analog, btype='low', analog=True, output='sos')

# Student then expects the digital cutoff to be at 3000 Hz

1. Explain why the cutoff will not be at 3 kHz.

2. Write the correct pre-warping formula and show what wc_analog should be.

3. Demonstrate both designs by plotting their digital frequency responses.

Solution

1. The bilinear transform compresses the entire analog frequency axis \([0, \infty)\) onto the digital interval \([0, \pi)\) through the nonlinear mapping \(\omega_a = (2/T)\tan(\omega_d T/2)\). Passing the raw digital radian frequency \(2\pi f_c\) as if it were already the analog frequency ignores this warping; the resulting digital cutoff will be shifted.

2. The correct pre-warp formula is:

\[\omega_a = \frac{2}{T}\tan\!\left(\frac{\omega_d T}{2}\right) = 2f_s \tan\!\left(\frac{\pi f_c}{f_s}\right)\]

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, sosfreqz, bilinear_zpk

fs = 8000
T  = 1 / fs
fc = 3000  # Hz — desired digital cutoff (close to Nyquist to make warping visible)

# Correct pre-warped analog frequency
wc_correct = 2 * fs * np.tan(np.pi * fc / fs)

# Student's (wrong) analog frequency — used ω_d directly
wc_wrong = 2 * np.pi * fc

print(f"Digital cutoff target:         {fc} Hz")
print(f"Wrong analog freq (no warp):   {wc_wrong:.0f} rad/s")
print(f"Correct pre-warped analog freq: {wc_correct:.0f} rad/s")
print(f"Ratio (correct/wrong):         {wc_correct/wc_wrong:.3f}")

# Correct: scipy handles pre-warping internally
sos_correct = butter(4, fc, btype='low', fs=fs, output='sos')

# Wrong: design analog prototype at wrong frequency, then apply bilinear manually
from scipy.signal import butter as butter_analog, zpk2sos
z_w, p_w, k_w = butter_analog(4, wc_wrong, btype='low', output='zpk', analog=True)
z_d, p_d, k_d = bilinear_zpk(z_w, p_w, k_w, fs)
sos_wrong = zpk2sos(z_d, p_d, k_d)

w, H_c = sosfreqz(sos_correct, worN=4096, fs=fs)
_, H_w  = sosfreqz(sos_wrong,   worN=4096, fs=fs)

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(w, 20*np.log10(np.maximum(np.abs(H_c), 1e-10)), 'C0-', linewidth=2,
        label='Correct (pre-warped)')
ax.plot(w, 20*np.log10(np.maximum(np.abs(H_w), 1e-10)), 'C3--', linewidth=1.5,
        label='Wrong (no pre-warp)')
ax.axvline(fc, color='k', linewidth=0.8, linestyle=':', label=f'Target {fc} Hz')
ax.axhline(-3, color='gray', linewidth=0.5, linestyle='--')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_xlim(0, fs/2)
ax.set_ylim(-60, 5)
ax.set_title('Effect of frequency pre-warping in bilinear transform')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

Digital cutoff target:         3000 Hz
Wrong analog freq (no warp):   18850 rad/s
Correct pre-warped analog freq: 38627 rad/s
Ratio (correct/wrong):         2.049

The key takeaway: SciPy’s butter(N, fc, fs=fs) performs the pre-warping internally. The bug arises only if you use the analog=True path and manually pass a digital frequency. Always let SciPy handle the bilinear transform unless you have a specific reason to work in the analog domain.

Challenge (continued)

Exercise 16: Data-driven: extract heart-rate from a PPG trace (Challenge)

A photoplethysmography (PPG) sensor measures blood volume changes at \(f_s = 100\) Hz. The raw signal contains:

• A strong DC baseline (offset)
• A fundamental heart-rate component in the range 0.8–3.5 Hz (48–210 bpm)
• Motion artefacts below 0.5 Hz
• High-frequency noise above 8 Hz

1. Design a bandpass filter that passes 0.8–3.5 Hz and attenuates outside 0.5–8 Hz by at least 40 dB. Choose the IIR family and order using buttord or an equivalent *ord function.

2. Apply the filter with zero-phase processing. Synthesise a test PPG signal to validate (use a 1.2 Hz fundamental plus harmonics, a 0.1 Hz drift, and broadband noise).

3. Estimate the heart rate from the filtered signal by finding the dominant frequency in the 0.8–3.5 Hz band using Welch’s method. Report the result in bpm.

Solution

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import buttord, butter, sosfiltfilt, welch

rng = np.random.default_rng(7)
fs  = 100
T   = 60   # seconds of data
N   = T * fs
t   = np.arange(N) / fs

# b) Synthesise test PPG
f_hr = 1.2   # Hz — heart rate (72 bpm)
ppg_clean = (  1.00 * np.sin(2*np.pi*f_hr*t)
             + 0.30 * np.sin(2*np.pi*2*f_hr*t)
             + 0.10 * np.sin(2*np.pi*3*f_hr*t))
drift     = 2.0 * np.sin(2*np.pi*0.1*t)          # slow motion artefact
noise     = 0.3 * rng.standard_normal(N)
x = ppg_clean + drift + noise

# a) Bandpass filter design — Butterworth (smooth passband, no ripple)
#    Use two-step: highpass then lowpass, or directly as bandpass
fp  = [0.8, 3.5]   # passband edges Hz
fst = [0.5, 8.0]   # stopband edges Hz
rp, As = 1.0, 40

order, Wn = buttord(fp, fst, rp, As, fs=fs)
print(f"Butterworth bandpass order: {order}, Wn = {Wn}")
sos = butter(order, Wn, btype='bandpass', fs=fs, output='sos')

# Apply zero-phase filter
y = sosfiltfilt(sos, x)

# c) Estimate heart rate
f_psd, Pxx = welch(y, fs, nperseg=512)
mask = (f_psd >= 0.8) & (f_psd <= 3.5)
f_peak = f_psd[mask][np.argmax(Pxx[mask])]
hr_bpm = f_peak * 60
print(f"Estimated heart rate: {hr_bpm:.1f} bpm  (true: {f_hr*60:.1f} bpm)")

# Plot
fig, axes = plt.subplots(2, 1, figsize=(10, 7))

t_show = slice(0, 10 * fs)
axes[0].plot(t[t_show], x[t_show], 'C7-', linewidth=0.8, alpha=0.6, label='Raw PPG')
axes[0].plot(t[t_show], y[t_show], 'C0-', linewidth=1.5, label='Filtered')
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('PPG signal — raw vs filtered (first 10 s)')
axes[0].legend(fontsize=8)
axes[0].grid(True, alpha=0.3)

axes[1].plot(f_psd, Pxx, 'C0-', linewidth=1.5)
axes[1].axvline(f_peak, color='C3', linewidth=1.5, linestyle='--',
                label=f'Peak: {f_peak:.2f} Hz = {hr_bpm:.0f} bpm')
axes[1].axvspan(0.8, 3.5, alpha=0.08, color='C0', label='HR band')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('PSD')
axes[1].set_xlim(0, 10)
axes[1].set_title('Welch PSD of filtered signal')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)

fig.tight_layout()
plt.show()

Butterworth bandpass order: 9, Wn = [0.76230244 3.67165781]
Estimated heart rate: 70.3 bpm  (true: 72.0 bpm)

The dominant PSD peak falls within 0.01 Hz of the true heart-rate frequency. Key design choices: Butterworth for a smooth passband (no ripple distortion on the waveform shape), zero-phase sosfiltfilt to avoid phase shift between beats, and Welch’s method to average out noise in the PSD estimate.

Exercise 17: Design triangle: order, ripple, transition width (Challenge)

The design triangle states that for a fixed filter family and filter type, you can freely choose any two of: (1) transition width, (2) ripple/attenuation, (3) filter order, and the third is then determined.

1. For a Parks–McClellan FIR lowpass filter with \(f_\text{pass} = 1\) kHz, \(f_\text{samp} = 10\) kHz: use remez and freqz to sweep filter orders \(N \in \{21, 41, 81, 161\}\). For each order, measure the achieved stopband attenuation (in dB) and transition width (frequency where the magnitude first drops to −40 dB). Record both in a table.

2. You are handed a fixed hardware budget: the processor can execute exactly 80 multiply-accumulate operations per sample. Accounting for symmetric FIR structure (roughly halving the multiplies), what is the maximum FIR order that fits?

3. Under that order constraint, what stopband attenuation is achievable with a 500 Hz transition band? With a 200 Hz transition band? Comment on the design trade-off.

Solution

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import remez, freqz

fs = 10000

def measure_filter(N, fp, fs, plot=False):
    """Design PM filter and measure attenuation and -40 dB transition width."""
    # Place stopband edge at fp + 500 Hz as a starting reference
    fst = fp + 500
    fst = min(fst, fs/2 - 10)
    h = remez(N, [0, fp, fst, fs/2], [1, 0], fs=fs)
    w, H = freqz(h, [1], worN=8192, fs=fs)
    mag_dB = 20 * np.log10(np.maximum(np.abs(H), 1e-12))

    # Stopband attenuation: worst case above fst
    sb_mask = w > fst
    sb_atten = -np.max(mag_dB[sb_mask]) if np.any(sb_mask) else np.nan

    # -40 dB transition point
    idx_40 = np.where(mag_dB <= -40)[0]
    f_40dB = w[idx_40[0]] if len(idx_40) > 0 else np.nan
    tw_40  = f_40dB - fp if not np.isnan(f_40dB) else np.nan

    return h, w, mag_dB, sb_atten, tw_40

fp = 1000
orders = [21, 41, 81, 161]
print(f"{'N':>5}  {'Stopband atten (dB)':>20}  {'Transition to −40 dB (Hz)':>26}")
print("-" * 57)

fig, ax = plt.subplots(figsize=(10, 5))
results = []
for N in orders:
    h, w, mag_dB, sb_atten, tw = measure_filter(N, fp, fs)
    results.append((N, sb_atten, tw))
    print(f"{N:>5}  {sb_atten:>20.1f}  {tw:>26.0f}")
    ax.plot(w, mag_dB, linewidth=1.5, label=f'N={N}')

ax.axhline(-40, color='k', linewidth=0.5, linestyle='--', label='−40 dB')
ax.axvline(fp,  color='gray', linewidth=0.5, linestyle=':')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_ylim(-100, 5)
ax.set_xlim(0, fs/2)
ax.set_title('Parks–McClellan FIR: effect of filter order')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

    N   Stopband atten (dB)   Transition to −40 dB (Hz)
---------------------------------------------------------
   21                  24.9                         545
   41                  39.3                         501
   81                  69.6                         454
  161                 127.6                         404

2. A length-\(N\) FIR filter with symmetric coefficients needs \(\lceil N/2 \rceil\) multiplications per output sample. For 80 MACs: \(\lceil N/2 \rceil \leq 80 \Rightarrow N \leq 159\). So the maximum order is \(N = 159\) (odd, to maintain linear phase).

3. Design at \(N = 159\) with the two transition bands:

fs = 10000
fp = 1000
N  = 159

def atten_for_transition(N, fp, delta_f, fs):
    fst = fp + delta_f
    if fst >= fs/2:
        return np.nan
    h = remez(N, [0, fp, fst, fs/2 - 1], [1, 0], fs=fs)
    _, H = freqz(h, [1], worN=8192, fs=fs)
    mag_dB = 20 * np.log10(np.maximum(np.abs(H), 1e-12))
    w = np.linspace(0, fs/2, 8192)
    sb_mask = w > fst
    return -np.max(mag_dB[sb_mask]) if np.any(sb_mask) else np.nan

for delta_f in [500, 200]:
    atten = atten_for_transition(N, fp, delta_f, fs)
    print(f"Transition {delta_f} Hz, N={N}: stopband attenuation ≈ {atten:.1f} dB")

Transition 500 Hz, N=159: stopband attenuation ≈ 127.3 dB
Transition 200 Hz, N=159: stopband attenuation ≈ 58.2 dB

With the 500 Hz transition band the hardware-limited filter achieves good attenuation; narrowing to a 200 Hz transition band (without increasing \(N\)) forces the equiripple algorithm to trade attenuation for transition sharpness: the stopband floor rises (worse rejection). This directly illustrates the design triangle: at fixed order, narrowing the transition band costs attenuation.

Exercise 18: Debug: zero-phase filtering on a causal real-time system (Challenge)

A student processes a streaming sensor signal in fixed-size blocks of 256 samples at \(f_s = 1000\) Hz. They copy a filter call from an offline analysis script and notice the output lags behind by roughly the expected group delay, but their manager insists the latency is unacceptable.

Here is the relevant code:

from scipy.signal import butter, sosfiltfilt
import numpy as np

sos = butter(4, 50, btype='low', fs=1000, output='sos')

def process_block(block, sos):
    return sosfiltfilt(sos, block)   # ← used in real-time loop

1. Identify the two distinct errors that make sosfiltfilt inappropriate here.

2. Propose a correct real-time implementation using sosfilt with persistent state across blocks. Sketch the stateful block-processing loop.

3. With sosfilt (causal, not zero-phase), what is the group delay of the 4th-order Butterworth filter at 20 Hz? At DC? (Use group_delay from scipy.)

Solution

1. Two problems:

1. sosfiltfilt is non-causal: it filters the block forwards and then backwards. The backward pass requires future samples that have not yet arrived in a streaming system. You cannot run sosfiltfilt on a live stream block-by-block.

2. No state is preserved across blocks: sosfiltfilt (and sosfilt if called naively) initialises the filter state to zero at the start of each block. Discontinuities at block boundaries produce artefacts (clicks, impulse-like transients) in the output.

2. Correct real-time implementation with persistent state:

import numpy as np
from scipy.signal import butter, sosfilt, sosfiltfilt

fs  = 1000
sos = butter(4, 50, btype='low', fs=fs, output='sos')

def process_block_correct(block, sos, zi):
    """Apply causal IIR filter preserving state across calls."""
    y, zi_out = sosfilt(sos, block, zi=zi)
    return y, zi_out

# Simulate block-by-block processing
rng = np.random.default_rng(99)
signal_full  = np.sin(2*np.pi*10*np.arange(3*fs)/fs) + 0.5*rng.standard_normal(3*fs)
block_size   = 256
output_rt    = np.zeros_like(signal_full)

zi = np.zeros((sos.shape[0], 2))   # reset state
for start in range(0, len(signal_full), block_size):
    block = signal_full[start:start+block_size]
    y_block, zi = process_block_correct(block, sos, zi)
    output_rt[start:start+len(y_block)] = y_block

# Compare with offline zero-phase for reference
output_zp = sosfiltfilt(sos, signal_full)

import matplotlib.pyplot as plt
t = np.arange(3*fs) / fs
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(t[:300], signal_full[:300], 'C7-', linewidth=0.8, alpha=0.5, label='Input')
ax.plot(t[:300], output_rt[:300],   'C0-', linewidth=1.5, label='Real-time (causal)')
ax.plot(t[:300], output_zp[:300],   'C3--', linewidth=1,  label='Offline zero-phase')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title('Real-time vs offline filtering (first 300 samples)')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

The causal real-time output is shifted to the right relative to the zero-phase version; the group delay is visible as a time offset. The waveform shape is preserved; there are no block-boundary glitches because zi carries the filter memory forward.

3. Group delay of the causal filter:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, group_delay

fs  = 1000
b, a = butter(4, 50, btype='low', fs=fs)  # need (b,a) for group_delay

w, gd = group_delay((b, a), w=4096, fs=fs)

for f_query in [0, 20, 50]:
    idx = np.argmin(np.abs(w - f_query))
    print(f"Group delay at {f_query:3d} Hz: {gd[idx]:.2f} samples = {gd[idx]/fs*1000:.2f} ms")

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(w, gd, 'C0-', linewidth=1.5)
ax.axvline(50, color='C3', linewidth=0.8, linestyle='--', label='Cutoff 50 Hz')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Group delay [samples]')
ax.set_xlim(0, 200)
ax.set_title('Group delay — 4th-order Butterworth lowpass')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

Group delay at   0 Hz: 8.25 samples = 8.25 ms
Group delay at  20 Hz: 8.94 samples = 8.94 ms
Group delay at  50 Hz: 11.95 samples = 11.95 ms

/tmp/ipykernel_2950/4038571414.py:8: UserWarning: The filter's denominator is extremely small at frequencies [3.141], around which a singularity may be present
  w, gd = group_delay((b, a), w=4096, fs=fs)

The group delay is frequency-dependent for an IIR filter, unlike a linear-phase FIR where different frequencies experience different delays. Near DC and deep in the passband the delay is roughly constant but it rises sharply near the cutoff frequency. This is the fundamental cost of IIR phase response: phase distortion that cannot be avoided in a causal implementation.