Lab 1: Filter Design Comparative -- Parks-McClellan FIR vs IIR vs Adaptive LMS

Week: 1-2 -- Advanced DSP
Points: 25
Time estimate: 90 min lab + 3 hr independent
Deliverable: lab-1-report.md + Python scripts

Objectives

1. Design a lowpass filter meeting a given specification using three methods: Parks-McClellan FIR, IIR-from-Butterworth, and adaptive LMS.
2. Measure and compare the frequency response, group delay, and order (coefficient count) for each filter.
3. Apply all three filters to the same test signal (tone in noise); compare SNR improvement.
4. Build the Parks-McClellan filter in GNU Radio and verify the response against the SciPy design.

Specification

Design a lowpass filter with these requirements:

• Sample rate: 1 MHz
• Passband: 0 - 100 kHz (ripple ≤ 0.5 dB)
• Transition band: 100 kHz - 150 kHz
• Stopband: 150 kHz and above (attenuation ≥ 60 dB)
• Linear phase required for FIR; not required for IIR

Part A: Parks-McClellan FIR and Kaiser Window (45 min)

#!/usr/bin/env python3
"""Lab 1 Part A: Parks-McClellan FIR and Kaiser Window FIR comparison."""
from scipy.signal import remez, freqz, kaiserord
import numpy as np
import matplotlib.pyplot as plt

fs = 1e6
f_pass = 100e3
f_stop = 150e3
atten_db = 60
ripple_db = 0.5
nyq = fs / 2

# ---- Parks-McClellan (Remez exchange) ----
# Use 'numtaps' as a starting point; adjust until spec is met
# Trial: start with N=80 and verify
for N in [60, 80, 100, 120]:
    h_pm = remez(
        numtaps=N,
        bands=[0, f_pass/nyq, f_stop/nyq, 1.0],
        desired=[1, 0],
        weight=[1, 10**(atten_db/20)]
    )
    w, H = freqz(h_pm, worN=8192)
    f_axis = w * fs / (2 * np.pi)
    
    # Check passband ripple
    passband_mask = f_axis <= f_pass
    passband_ripple = 20*np.log10(np.max(np.abs(H[passband_mask]))) - \
                      20*np.log10(np.min(np.abs(H[passband_mask])))
    
    # Check stopband attenuation
    stopband_mask = f_axis >= f_stop
    stopband_atten = -20*np.log10(np.max(np.abs(H[stopband_mask])) + 1e-15)
    
    print(f"PM N={N:3d}: passband ripple={passband_ripple:.2f} dB, "
          f"stopband attenuation={stopband_atten:.1f} dB")
    
    if passband_ripple <= ripple_db and stopband_atten >= atten_db:
        print(f"  → SPEC MET at N={N}")
        h_pm_final = h_pm
        N_pm = N
        break
else:
    h_pm_final = remez(numtaps=120, bands=[0, f_pass/nyq, f_stop/nyq, 1.0],
                       desired=[1, 0], weight=[1, 10**(atten_db/20)])
    N_pm = 120

# ---- Kaiser Window ----
N_kaiser, beta = kaiserord(atten_db, (f_stop - f_pass) / nyq)
h_kaiser = np.sinc(2 * f_pass/fs * (np.arange(N_kaiser) - (N_kaiser-1)/2))
h_kaiser *= np.kaiser(N_kaiser, beta)
h_kaiser /= np.sum(h_kaiser)

print(f"\nFilter orders: Parks-McClellan N={N_pm}, Kaiser N={N_kaiser}")
print(f"PM uses {N_pm/(N_kaiser):.2f}× fewer taps than Kaiser")

# ---- Group delay comparison ----
from scipy.signal import group_delay

w_gd, gd_pm = group_delay((h_pm_final, 1), w=8192)
w_gd, gd_kaiser = group_delay((h_kaiser, 1), w=8192)
f_gd = w_gd * fs / (2 * np.pi)

# Report: constant group delay for linear-phase FIR
passband_gd = f_gd <= f_pass
print(f"\nGroup delay in passband:")
print(f"  PM:     {np.mean(gd_pm[passband_gd]):.1f} ± {np.std(gd_pm[passband_gd]):.3f} samples")
print(f"  Kaiser: {np.mean(gd_kaiser[passband_gd]):.1f} ± {np.std(gd_kaiser[passband_gd]):.3f} samples")

# Save filter coefficients for Part D (GNU Radio)
np.save('lab1/h_pm.npy', h_pm_final)

Record:

1. Minimum filter order (N) for each method that meets the spec
2. Passband ripple achieved (dB)
3. Stopband attenuation achieved (dB)
4. Group delay variance (should be ~0 for linear-phase FIR)

Part B: IIR-from-Butterworth (20 min)

#!/usr/bin/env python3
"""Lab 1 Part B: IIR filter design comparison."""
from scipy.signal import butter, cheby1, ellip, sosfilt, sosfreqz, group_delay
import numpy as np
import matplotlib.pyplot as plt

fs = 1e6
f_cutoff = 100e3
nyq = fs / 2

# Try different orders for Butterworth to meet 60 dB stopband at 150 kHz
# Butterworth -3dB at f_cutoff; check attenuation at 150 kHz
for N in range(3, 12):
    sos = butter(N, f_cutoff/nyq, btype='low', output='sos')
    w, H = sosfreqz(sos, worN=8192)
    f_axis = w * fs / (2 * np.pi)
    
    stopband_mask = f_axis >= 150e3
    if np.any(stopband_mask):
        atten = -20*np.log10(np.max(np.abs(H[stopband_mask])) + 1e-15)
        print(f"Butterworth N={N}: stopband atten at 150 kHz = {atten:.1f} dB")
        if atten >= 60:
            print(f"  → 60 dB SPEC MET at N={N}")
            sos_butter = sos
            N_butter = N
            break

# Elliptic for comparison (minimum order for same spec)
for N in range(3, 12):
    sos_ellip = ellip(N, 0.5, 60, f_cutoff/nyq, btype='low', output='sos')
    w, H = sosfreqz(sos_ellip, worN=8192)
    f_axis = w * fs / (2 * np.pi)
    stopband_mask = f_axis >= 150e3
    if np.any(stopband_mask):
        atten = -20*np.log10(np.max(np.abs(H[stopband_mask])) + 1e-15)
        if atten >= 60:
            print(f"Elliptic N={N}: 60 dB spec met, passband ripple ≤0.5 dB")
            N_ellip = N
            break

print(f"\nIIR order comparison: Butterworth N={N_butter}, Elliptic N={N_ellip}")

# Group delay of IIR (NOT constant -- show the non-linearity)
_, gd_butter = group_delay((sos_butter[0,:3], sos_butter[0,3:]), w=8192)
print(f"\nButterworth group delay in passband: varies from "
      f"{np.min(gd_butter[:100]):.1f} to {np.max(gd_butter[:100]):.1f} samples")
print("(Non-constant group delay is why IIR is not preferred for demodulation)")

Record:

1. Butterworth order N required to meet the spec
2. Elliptic order N required for same spec (fewer taps = more efficient)
3. Group delay range in the passband -- quantify the non-linearity

Part C: Adaptive LMS Filter (25 min)

The adaptive filter task: given a signal with a tonal interferer at an unknown frequency, use LMS to notch it without knowing the interferer frequency in advance.

#!/usr/bin/env python3
"""Lab 1 Part C: Adaptive LMS for interference cancellation."""
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
fs = 1e6
N_samples = 50000
t = np.arange(N_samples) / fs

# Target signal: 50 kHz tone (inside passband)
target = np.sin(2 * np.pi * 50e3 * t)

# Interferer: 200 kHz tone (in stopband) -- unknown frequency
interferer = 2.0 * np.sin(2 * np.pi * 200e3 * t)

# Noise
noise = 0.1 * np.random.randn(N_samples)

# Composite input
x = target + interferer + noise

# LMS adaptive notch filter
# Strategy: use x[n] and x[n-1] as input pair; predict the interferer;
# the prediction error is the interferer-suppressed signal
mu = 0.005  # step size
M = 32      # filter length
w = np.zeros(M)
y = np.zeros(N_samples)
e = np.zeros(N_samples)

for n in range(M, N_samples):
    x_buf = x[n:n-M:-1]   # input buffer
    y[n] = np.dot(w, x_buf)
    e[n] = x[n] - y[n]    # prediction error ≈ target + noise (interferer suppressed)
    w += mu * e[n] * x_buf

# Measure output SNR before and after
def snr_db(signal, noise_estimate):
    sig_power = np.var(signal)
    noise_power = np.var(noise_estimate)
    return 10 * np.log10(sig_power / noise_power) if noise_power > 0 else float('inf')

# SNR of raw x (only include steady-state region)
ss_start = 10000
input_snr = snr_db(target[ss_start:], interferer[ss_start:] + noise[ss_start:])
output_snr = snr_db(target[ss_start:], e[ss_start:] - target[ss_start:])
print(f"Input SNR (target vs. interferer): {input_snr:.1f} dB")
print(f"Output SNR (after LMS adaptation): {output_snr:.1f} dB")
print(f"SNR improvement: {output_snr - input_snr:.1f} dB")

# Check convergence
plt.figure(figsize=(10, 4))
plt.plot(np.arange(N_samples)/fs*1e3, np.abs(e)**2, alpha=0.3, label='Error power')
plt.xlabel('Time (ms)')
plt.ylabel('Error power')
plt.title('LMS Convergence')
plt.yscale('log')
plt.grid(True)
plt.savefig('lab1/lms_convergence.png', dpi=150)

Record:

1. Convergence time (approximately how many samples before error stabilizes)
2. SNR improvement achieved
3. What happens if you increase μ to 0.05? Decrease to 0.001?

Part D: GNU Radio Verification (optional extension, +0 pts)

Load the Parks-McClellan coefficients into a GNU Radio Frequency Xlating FIR Filter block and apply them to a live RTL-SDR capture. Verify the frequency response matches the SciPy design.

# Convert h_pm.npy to GNU Radio compatible format
import numpy as np
h = np.load('lab1/h_pm.npy')
# In GRC: Low Pass Filter → custom taps: paste h as Python list
print("GNU Radio taps (first 10):", h[:10].tolist())

Lab Report

Create lab-1-report.md with:

1. Filter order comparison table:

2. LMS adaptation: convergence plot + SNR improvement table (μ = 0.005, 0.05, 0.001)
3. Analysis paragraph: "For an SDR demodulator (e.g., a QPSK demodulator), which filter would you choose and why? Specifically address group delay and filter order."

Grading