RF-301 Week 1-2 — Advanced DSP: Filter Design (FIR / IIR / Adaptive / FFT-Based)

"The Kaiser window method of FIR filter design offers a near-optimal trade-off between stopband attenuation and transition-band width. The parametric Kaiser window allows you to adjust the transition width and the stopband attenuation independently." — Understanding Digital Signal Processing, Richard Lyons (3rd ed., Ch 7)

Lecture (90 min × 2)

1.1 Why Filter Design Matters at RF-301 Scale

In RF-201, you used filters as black boxes: SciPy's butter() to remove noise from a LoRa chirp, the GNU Radio Low Pass Filter block with an empirically-chosen cutoff. At RF-201 register, "pick a filter that works" is sufficient.

At RF-301 register, that approach fails. Three scenarios make the distinction clear:

Scenario A: A cognitive-radio spectrum sensor. The sensor must continuously scan a 200 MHz band looking for occupancy. Its filter bank needs to resolve 200 kHz sub-channels across the band in real time at the SNR where occupancy decisions must be reliable. Filter transition bandwidth, passband ripple, and group delay all constrain what the cognitive radio can detect and how fast it can respond. "Empirically chosen" doesn't compose into a provable false-alarm rate.

Scenario B: LTE OFDMA receiver. An LTE base station's OFDMA receiver must separate subcarriers spaced 15 kHz apart at 50-75 dB of adjacent-subcarrier rejection to extract a single user's resource blocks. The filter spec is derived from the 3GPP standard, not from experimentation.

Scenario C: SIGINT capture. A wideband SIGINT receiver captures at 40 MSPS. Channelizing that wideband capture into 100 kHz channels requires a polyphase filter bank with thousands of taps. The design is an exercise in computational efficiency -- how many multiplications per sample -- not just frequency response.

These three scenarios define the three filter-design problems RF-301 opens: frequency-selective FIR design, IIR design from analog prototypes, and adaptive filtering. Week 1 opens FIR (windowed and Parks-McClellan); Week 2 opens IIR and adaptive.

Academy Flowgraph note: Use the in-browser Academy Flowgraph tool (available at portal.virtuscyberacademy.org/workbench/static/academy-flowgraph.html) to sketch filter chains visually before implementing in GNU Radio. The Filter block lets you set cutoff and observe the block-graph structure; the FFT block lets you visualize the frequency response. These are pedagogical sketches, not DSP-accurate simulations -- move to GNU Radio for quantitative verification.

1.2 FIR Filter Design: The Design-by-Specification Problem

A linear-phase FIR filter of order N is:

H(z) = h[0] + h[1]z^{-1} + h[2]z^{-2} + ... + h[N]z^{-N}

The coefficients h[n] are the filter's impulse response (because the output is the convolution of input with h[n]). Linear phase is guaranteed when the coefficients are symmetric: h[n] = h[N-n]. That symmetry gives you constant group delay -- all frequency components are delayed equally, preserving waveshape. For demodulation, constant group delay is essential: a filter with frequency-dependent delay distorts the recovered symbol.

The design problem: Given a specification (passband edge, stopband edge, ripple tolerances), find the coefficient sequence h[n] that meets the spec with minimum filter order (minimum N = minimum computation).

Four FIR design methods:

Parks-McClellan is the method of choice for tight specs because it minimizes filter order for given ripple constraints. For any combination of passband ripple δ_p and stopband attenuation δ_s, Parks-McClellan gives you the minimum N filter -- the mathematically optimal equiripple solution. Lyons Ch 7 walks the Remez exchange algorithm; the key insight is that the optimal filter's approximation error equiripples (alternates) between the passband and stopband bands at exactly N+2 extrema (Chebyshev equiripple property).

SciPy implementation:

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

# Specification
fs = 1e6          # sample rate 1 MHz
f_pass = 100e3    # passband edge 100 kHz
f_stop = 150e3    # stopband edge 150 kHz
atten_db = 60     # stopband attenuation 60 dB
ripple_db = 0.1   # passband ripple 0.1 dB

# Convert to normalized (0..0.5) frequency
nyq = fs / 2

# Parks-McClellan (Remez exchange)
# bands: [0, f_pass, f_stop, nyq] normalized to [0, 1]
desired = [1, 0]     # passband gain=1, stopband gain=0
h_pm = remez(
    numtaps=64,      # initial guess; verify with freqz
    bands=[0, f_pass/nyq, f_stop/nyq, 1],
    desired=desired,
    weight=[1, 10**((atten_db - ripple_db)/20)]  # weight stopband heavier
)

# Frequency response
w, H = freqz(h_pm, worN=2048)
f = w * fs / (2 * np.pi)

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(f/1e3, 20*np.log10(np.abs(H) + 1e-10))
plt.xlabel('Frequency (kHz)')
plt.ylabel('Magnitude (dB)')
plt.title('Parks-McClellan FIR')
plt.axvline(f_pass/1e3, color='g', linestyle='--', label='Passband edge')
plt.axvline(f_stop/1e3, color='r', linestyle='--', label='Stopband edge')
plt.legend()
plt.xlim([0, 400])
plt.ylim([-80, 5])
plt.grid(True)

# Kaiser window method (comparison)
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)

w2, H2 = freqz(h_kaiser, worN=2048)
f2 = w2 * fs / (2 * np.pi)

plt.subplot(1, 2, 2)
plt.plot(f2/1e3, 20*np.log10(np.abs(H2) + 1e-10))
plt.xlabel('Frequency (kHz)')
plt.ylabel('Magnitude (dB)')
plt.title(f'Kaiser Window FIR (N={N_kaiser}, β={beta:.1f})')
plt.axvline(f_pass/1e3, color='g', linestyle='--', label='Passband edge')
plt.axvline(f_stop/1e3, color='r', linestyle='--', label='Stopband edge')
plt.legend()
plt.xlim([0, 400])
plt.ylim([-80, 5])
plt.grid(True)

plt.tight_layout()
plt.savefig('filter_comparison.png', dpi=150)
print(f'Parks-McClellan N=64, Kaiser N={N_kaiser}')

The key comparison: for the same spec, Parks-McClellan achieves the spec with fewer taps than Kaiser window in most cases. Fewer taps = fewer multiplications per sample = less computation.

1.3 IIR Filter Design: From Analog Prototype to Digital

The IIR design approach: start with a classical analog filter prototype (Butterworth, Chebyshev Type I/II, Elliptic), then transform it to discrete-time using bilinear transform or impulse invariance.

The Butterworth prototype: Maximally flat passband; no ripple. The Nth-order Butterworth filter has magnitude:

|H(jω)|² = 1 / (1 + (ω/ω_c)^(2N))

At the cutoff frequency ω_c, the response is exactly -3 dB. The slope in the stopband is -20N dB/decade. A 5th-order Butterworth achieves -100 dB/decade.

Chebyshev Type I: Ripple in the passband; maximally steep rolloff for given passband ripple and filter order. Better stopband attenuation than Butterworth for same order.

Chebyshev Type II: Ripple in the stopband (not passband); sometimes preferable when passband flatness is critical.

Elliptic (Cauer): Ripple in both passband and stopband. Steepest rolloff of all four for given order. Minimum order for given spec -- but the group delay is highly non-linear near the band edge.

Bilinear transform: Maps analog s-plane poles to z-plane with frequency warping. The transformation:

s = (2/T) · (z-1)/(z+1)

maps the left s-half-plane to the inside of the unit circle (stability preserved). It warps the frequency axis -- the relationship between analog frequency Ω and digital frequency ω is Ω = (2/T)·tan(ω/2). Pre-warp the critical frequencies before design to compensate.

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

fs = 1e6
f_cutoff = 100e3   # -3 dB cutoff for Butterworth

# Design 5th-order filters using second-order sections (sos) for numerical stability
sos_butter = butter(5, f_cutoff / (fs/2), btype='low', output='sos')
sos_cheby1 = cheby1(5, 0.5, f_cutoff / (fs/2), btype='low', output='sos')  # 0.5 dB ripple
sos_ellip  = ellip(5, 0.5, 60, f_cutoff / (fs/2), btype='low', output='sos')  # 0.5 dB/60 dB

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, sos, label in zip(axes,
                          [sos_butter, sos_cheby1, sos_ellip],
                          ['Butterworth N=5', 'Chebyshev-I N=5\n(0.5 dB ripple)', 'Elliptic N=5\n(0.5/60 dB)']):
    w, H = sosfreqz(sos, worN=4096)
    f = w * fs / (2 * np.pi)
    ax.plot(f/1e3, 20*np.log10(np.abs(H) + 1e-10))
    ax.axvline(f_cutoff/1e3, color='r', linestyle='--', label=f'f_c={f_cutoff/1e3:.0f} kHz')
    ax.set_xlabel('Frequency (kHz)')
    ax.set_ylabel('Magnitude (dB)')
    ax.set_title(label)
    ax.set_xlim([0, 400])
    ax.set_ylim([-100, 5])
    ax.grid(True)

plt.tight_layout()
plt.savefig('iir_comparison.png', dpi=150)

When to use IIR vs FIR: IIR filters achieve steeper rolloff at lower order (fewer coefficients, less computation). The cost: non-linear phase response (group delay varies with frequency -- problematic for demodulation). FIR with linear phase is preferred for demodulation; IIR is preferred for anti-aliasing, noise rejection where phase linearity is not critical, and real-time processing where coefficient count is constrained.

1.4 Adaptive Filters: LMS and RLS

A static FIR or IIR filter has fixed coefficients -- the response is the same regardless of the input signal statistics. An adaptive filter updates its coefficients in real time to minimize an error signal.

Use cases in RF:

• Noise cancellation: A secondary antenna captures the noise only; the adaptive filter subtracts it from the primary antenna's received signal
• Channel equalization: The channel distorts the transmitted signal; the adaptive filter estimates and inverts the channel response
• Interference suppression: A narrowband interferer appears in a wideband channel; the adaptive filter notches it without a priori frequency knowledge

Least Mean Squares (LMS) algorithm:

The LMS algorithm updates the filter weights w[n] after each sample:

y[n]   = w[n]^T · x[n]          # filter output
e[n]   = d[n] - y[n]            # error (desired - output)
w[n+1] = w[n] + μ · e[n] · x[n] # weight update

where μ is the step size (learning rate). The LMS rule is the gradient descent of the mean squared error E[e²[n]]. Smaller μ = slower convergence but more stable; larger μ = faster convergence but can diverge if μ > 2/λ_max where λ_max is the maximum eigenvalue of the input autocorrelation matrix.

import numpy as np
import matplotlib.pyplot as plt

def lms_filter(x, d, mu, N):
    """LMS adaptive filter.
    x:  input signal
    d:  desired signal
    mu: step size
    N:  filter length (taps)
    Returns: y (output), e (error), W (final weights)
    """
    n_samples = len(x)
    w = np.zeros(N)      # initial weights
    y = np.zeros(n_samples)
    e = np.zeros(n_samples)
    
    for n in range(N, n_samples):
        x_vec = x[n:n-N:-1]  # input buffer (most recent first)
        y[n] = np.dot(w, x_vec)
        e[n] = d[n] - y[n]
        w = w + mu * e[n] * x_vec  # LMS weight update
    
    return y, e, w

# Demonstration: noise cancellation
np.random.seed(42)
n = np.arange(2000)
signal = np.sin(2 * np.pi * 0.05 * n)           # desired 50 Hz sine
noise = 0.5 * np.random.randn(len(n))            # additive noise
x = signal + noise                                # noisy input
d = signal                                        # desired output (clean signal)

y, e, w = lms_filter(x + 0.3*np.random.randn(len(n)), d, mu=0.01, N=32)

plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(n[:200], x[:200], alpha=0.5, label='Noisy input')
plt.plot(n[:200], d[:200], label='Desired (clean)')
plt.legend()
plt.title('LMS Adaptive Filter: Noise Cancellation')
plt.subplot(2, 1, 2)
plt.plot(n[:200], e[:200], label='Error (before convergence)')
plt.plot(n[1800:], e[1800:], label='Error (after convergence)', color='r')
plt.legend()
plt.tight_layout()
plt.savefig('lms_convergence.png', dpi=150)

RLS (Recursive Least Squares): Converges faster than LMS at the cost of more computation per sample (O(N²) vs O(N)). Uses the exact least-squares solution updated recursively via the matrix inversion lemma. Preferred when convergence speed matters more than computational efficiency.

1.5 FFT-Based Processing: Spectral Analysis and Overlap-Save

For block processing -- when you have a batch of samples, not a sample-by-sample stream -- the FFT offers O(N log N) convolution instead of O(N²) direct convolution.

Overlap-save convolution: Divide the input stream into overlapping blocks; compute FFT of each block; multiply by the filter's frequency response; IFFT back. Each output block requires one FFT, one element-wise multiply, and one IFFT -- regardless of filter length.

Polyphase filter bank: A channelizer that decomposes a wideband input into M narrowband channels simultaneously, using M filters running at 1/M the input rate. The computational trick: all M filters share the same coefficient set, just polyphase-decomposed. Equivalent to M FIR filters running in parallel, but implemented as one filter at the original rate followed by an FFT.

import numpy as np
from scipy.signal import remez

def polyphase_channelizer(x, M, taps_per_channel):
    """Polyphase FFT channelizer: M channels, taps_per_channel per channel."""
    # Design prototype lowpass filter for one channel
    h = remez(
        numtaps=M * taps_per_channel,
        bands=[0, 0.4/M, 0.6/M, 1.0],
        desired=[1, 0],
        weight=[1, 100]
    )
    
    # Reshape into polyphase branches
    # h[m, k] = h[m + M*k]
    H = h.reshape(-1, M).T  # shape: M × taps_per_channel
    
    # Apply polyphase filter to decimated input
    x_padded = np.pad(x, (0, M - len(x) % M))
    x_blocks = x_padded.reshape(-1, M)  # each row: one M-sample input block
    
    # Filter each polyphase branch
    from scipy.signal import lfilter
    Y = np.zeros_like(x_blocks, dtype=complex)
    for m in range(M):
        Y[:, m] = lfilter(H[m], [1.0], x_blocks[:, m])
    
    # FFT across branches for each block
    channels = np.fft.fft(Y, axis=1)  # shape: blocks × M
    
    return channels  # channels[block, channel_index]

# Example: channelize a 1 MHz signal into 8 × 125 kHz channels
fs = 1e6
t = np.arange(10000) / fs
# Two tones at 150 kHz and 425 kHz
x = np.exp(2j * np.pi * 150e3 * t) + 0.5 * np.exp(2j * np.pi * 425e3 * t)

channels = polyphase_channelizer(x.real, M=8, taps_per_channel=16)
print(f"Output: {channels.shape[0]} blocks × {channels.shape[1]} channels")

Why this matters for SIGINT: A SIGINT receiver running a polyphase channelizer can simultaneously monitor M frequency channels with the computation of roughly one FIR filter. For the carrier-scale scenario (200 MHz band, 200 kHz channels = 1000 channels), this is the only computationally tractable approach.

1.6 Anchor Weave: Lyons Ch 6-7

Richard Lyons's Understanding Digital Signal Processing (3rd ed.) Chapters 6-7 carry the mathematical foundation for everything in this week's lecture. Chapter 6 (Finite Impulse Response Filters) walks the window method and introduces the Parks-McClellan optimality criterion in the context it belongs in: the equiripple approximation problem. Chapter 7 (Infinite Impulse Response Filters) derives the bilinear transform from the analog prototype starting point.

The reading commitment at Belt-5: Read Ch 6-7 before the Week 2 lab session. Lyons is careful about the transition between time-domain (coefficient design) and frequency-domain (specification) thinking; the two chapters together give you both sides. The lab will ask you to defend your filter order choices, which requires the Ch 7 analytical framework.

Cross-reference: Wyglinski et al. Ch 4 (Receiver Architectures) situates filter design in the receiver chain: each stage in the chain (RF bandpass, IF bandpass, baseband lowpass, channelizer) has a filter with a specification derived from the preceding stage's output. The Wyglinski framing turns Lyons's abstract filter math into engineering decisions you make at each receiver stage.

Lab Introduction

Lab 1 (25 pts): Filter design comparative lab -- Parks-McClellan FIR vs IIR-from-Butterworth vs adaptive LMS against the same channel model. See labs/lab-1.md.

Independent Practice

1. Lyons Ch 6: work through Examples 6-1 through 6-3 (windowed FIR design by hand; verify with Python)
2. Lyons Ch 7: trace the bilinear transform derivation for a 3rd-order Butterworth; confirm the digital pole locations
3. Implement the polyphase channelizer above with M=16 channels; verify that a tone at f_c lands in the correct channel bin
4. Explore the Academy Flowgraph Filter block: set the cutoff and tap count and observe the frequency response; sketch a 3-block signal chain (Source → Filter → FFT) in the browser tool before building it in GNU Radio