1#!/usr/bin/env python3
2"""
3IIR Filter Design Methods
4==========================
5
6This script demonstrates the four classical IIR lowpass filter designs:
7
81. Butterworth — maximally flat magnitude (no ripple in either band)
92. Chebyshev I — equiripple in passband, monotone in stopband
103. Chebyshev II — monotone in passband, equiripple in stopband
114. Elliptic — equiripple in both bands (minimum order for given spec)
12
13Key Concepts:
14 - Analog prototype → bilinear transform → digital IIR filter
15 - Filter order N controls transition band sharpness
16 - Passband ripple (Rp) and stopband attenuation (Rs) are design knobs
17 - Pole-zero plot reveals where the filter's gain is amplified (poles)
18 and nulled (zeros, for Chebyshev II and Elliptic)
19
20Design Trade-offs:
21 - Butterworth : Flattest passband, widest transition band, all-pole
22 - Chebyshev I : Sharper roll-off than Butterworth (same order),
23 ripple in passband
24 - Chebyshev II: Sharper roll-off, ripple in stopband,
25 finite zeros on unit circle
26 - Elliptic : Sharpest roll-off, ripple in both bands,
27 most sensitive to coefficient quantization
28
29Author: Educational example for Signal Processing
30License: MIT
31"""
32
33import numpy as np
34from scipy import signal
35import matplotlib.pyplot as plt
36
37# ============================================================================
38# DESIGN PARAMETERS
39# ============================================================================
40
41FS = 1000.0 # Sampling frequency (Hz)
42FC = 150.0 # Passband cutoff (Hz) — 3 dB point for Butterworth
43RP = 1.0 # Passband ripple (dB) — used by Cheby I and Elliptic
44RS = 40.0 # Stopband attenuation (dB) — Cheby II, Elliptic
45ORDER = 5 # Filter order (common across all types for fair comparison)
46FC_NORM = FC / (FS / 2) # Normalized cutoff (0–1, where 1 = Nyquist)
47
48
49# ============================================================================
50# DESIGN FUNCTIONS
51# ============================================================================
52
53def design_all_filters(order, fc_norm, rp_db, rs_db):
54 """
55 Design four classical IIR lowpass filters of the same order.
56
57 All filters use the bilinear transform (analog=False) so that the
58 design frequency fc_norm maps directly to the normalized digital
59 frequency axis (0–1 where 1 = Nyquist).
60
61 Args:
62 order (int) : Filter order N
63 fc_norm (float): Normalized cutoff (0–1)
64 rp_db (float): Passband ripple in dB (Cheby I, Elliptic)
65 rs_db (float): Stopband attenuation in dB (Cheby II, Elliptic)
66
67 Returns:
68 dict: {name: (b, a)} second-order sections not used here for clarity
69 """
70 filters = {}
71
72 # Butterworth: maximally flat, -3 dB exactly at fc_norm
73 b, a = signal.butter(order, fc_norm, btype='low', analog=False)
74 filters['Butterworth'] = (b, a)
75
76 # Chebyshev Type I: equiripple in passband, monotone stopband
77 b, a = signal.cheby1(order, rp_db, fc_norm, btype='low', analog=False)
78 filters['Chebyshev I'] = (b, a)
79
80 # Chebyshev Type II: monotone passband, equiripple stopband
81 # fc_norm here marks the beginning of the stopband (Rs attenuation)
82 b, a = signal.cheby2(order, rs_db, fc_norm, btype='low', analog=False)
83 filters['Chebyshev II'] = (b, a)
84
85 # Elliptic (Cauer): equiripple in both bands — minimum order
86 b, a = signal.ellip(order, rp_db, rs_db, fc_norm, btype='low', analog=False)
87 filters['Elliptic'] = (b, a)
88
89 return filters
90
91
92# ============================================================================
93# ANALYSIS FUNCTIONS
94# ============================================================================
95
96def compute_responses(filters, worN=2048):
97 """
98 Compute frequency responses for all filters.
99
100 Returns:
101 dict: {name: (w, H)} where w is normalized freq (0–π) and H is complex
102 """
103 responses = {}
104 for name, (b, a) in filters.items():
105 w, H = signal.freqz(b, a, worN=worN)
106 responses[name] = (w, H)
107 return responses
108
109
110def print_filter_stats(filters, fc_norm, rp_db, rs_db):
111 """Print key numerical metrics for each filter."""
112 print("=" * 70)
113 print(f"IIR LOWPASS FILTER COMPARISON (order={ORDER}, fc={FC:.0f} Hz, fs={FS:.0f} Hz)")
114 print("=" * 70)
115 print(f"{'Filter':<16} {'Taps':>5} {'3 dB point':>13} {'PB ripple':>12} {'SB atten':>10}")
116 print("-" * 70)
117
118 for name, (b, a) in filters.items():
119 w, H = signal.freqz(b, a, worN=4096)
120 mag = np.abs(H)
121 mag_db = 20 * np.log10(np.maximum(mag, 1e-15))
122
123 # -3 dB crossing
124 idx_3dB = np.argmin(np.abs(mag - (1.0 / np.sqrt(2))))
125 freq_3dB = w[idx_3dB] / np.pi # normalized (fraction of Nyquist)
126
127 # Passband ripple: max deviation from 0 dB up to fc_norm
128 idx_pb = np.searchsorted(w / np.pi, fc_norm)
129 pb_ripple = np.max(mag_db[:idx_pb + 1]) - np.min(mag_db[1:idx_pb + 1])
130
131 # Stopband attenuation: max gain past fc_norm + 20 % transition
132 idx_sb = np.searchsorted(w / np.pi, min(fc_norm * 1.4, 0.99))
133 sb_atten = -np.max(mag_db[idx_sb:])
134
135 taps = len(b) + len(a) - 1 # numerator + denominator - shared gain
136 print(f"{name:<16} {taps:>5} {freq_3dB:>8.3f}×Nyq "
137 f"{pb_ripple:>10.2f} dB {sb_atten:>8.1f} dB")
138
139 print()
140 print(f"Design specs: Rp={rp_db} dB passband ripple, Rs={rs_db} dB stopband atten.")
141 print(f"Note: Chebyshev II fc marks start of stopband, not -3 dB point.")
142
143
144def minimum_order_comparison(fc_norm, rp_db, rs_db):
145 """
146 Show what order each type needs to meet the same specification.
147
148 Specification: Rp dB ripple in passband up to fc_norm,
149 Rs dB attenuation in stopband from 1.3×fc_norm upward.
150 """
151 print("\n" + "=" * 70)
152 print("MINIMUM ORDER TO MEET SPEC (Rp={} dB, Rs={} dB, fc={} Hz)".format(
153 rp_db, rs_db, FC))
154 print("Stopband edge = 1.3 × passband edge")
155 print("=" * 70)
156
157 fs_norm = min(fc_norm * 1.3, 0.99) # stopband edge (normalized)
158
159 types = ['butter', 'cheby1', 'cheby2', 'ellip']
160 labels = ['Butterworth', 'Chebyshev I', 'Chebyshev II', 'Elliptic']
161
162 print(f"{'Filter':<16} {'Min order':>10}")
163 print("-" * 30)
164 for ftype, label in zip(types, labels):
165 try:
166 N, _ = signal.buttord(fc_norm, fs_norm, rp_db, rs_db) # not used
167 except Exception:
168 N = '—'
169
170 if ftype == 'butter':
171 N, Wn = signal.buttord(fc_norm, fs_norm, rp_db, rs_db)
172 elif ftype == 'cheby1':
173 N, Wn = signal.cheb1ord(fc_norm, fs_norm, rp_db, rs_db)
174 elif ftype == 'cheby2':
175 N, Wn = signal.cheb2ord(fc_norm, fs_norm, rp_db, rs_db)
176 else:
177 N, Wn = signal.ellipord(fc_norm, fs_norm, rp_db, rs_db)
178
179 print(f"{label:<16} {N:>10}")
180
181 print("\n → Elliptic always needs the fewest taps for any given spec.")
182 print(" → Butterworth needs the most (but has the flattest passband).")
183
184
185# ============================================================================
186# SIGNAL APPLICATION
187# ============================================================================
188
189def apply_filters_to_signal(filters):
190 """
191 Apply all four filters to a noisy test signal.
192
193 Signal: 50 Hz tone (in passband) + 300 Hz tone (in stopband)
194 + broadband Gaussian noise
195 """
196 print("\n" + "=" * 70)
197 print("FILTERING A NOISY TEST SIGNAL")
198 print("=" * 70)
199
200 rng = np.random.default_rng(42)
201 t = np.linspace(0, 0.5, int(FS * 0.5), endpoint=False)
202 f_pass = 50.0 # Hz — inside passband
203 f_stop = 300.0 # Hz — inside stopband
204
205 # Construct noisy signal
206 clean = np.sin(2 * np.pi * f_pass * t)
207 interference = 0.8 * np.sin(2 * np.pi * f_stop * t)
208 noise = 0.3 * rng.standard_normal(len(t))
209 x = clean + interference + noise
210
211 print(f" Signal : {f_pass} Hz (passband) + {f_stop} Hz (stopband) + noise")
212 print(f" Cutoff : {FC} Hz (filters should remove {f_stop} Hz)")
213
214 outputs = {}
215 for name, (b, a) in filters.items():
216 # zero-phase filtering removes group-delay distortion for comparison
217 y = signal.filtfilt(b, a, x)
218 outputs[name] = y
219
220 # Power of residual at stopband frequency
221 N = len(y)
222 Y = np.fft.rfft(y) / N
223 freqs = np.fft.rfftfreq(N, d=1.0 / FS)
224 idx = np.argmin(np.abs(freqs - f_stop))
225 stop_power_db = 20 * np.log10(max(np.abs(Y[idx]), 1e-15))
226 print(f" {name:<16}: residual at {f_stop:.0f} Hz = {stop_power_db:.1f} dBFS")
227
228 return t, x, outputs
229
230
231# ============================================================================
232# ENTRY POINT
233# ============================================================================
234
235if __name__ == "__main__":
236 print("IIR FILTER DESIGN METHODS")
237
238 # Design all four filters with the same order
239 filters = design_all_filters(ORDER, FC_NORM, RP, RS)
240 responses = compute_responses(filters)
241
242 # Print statistics
243 print_filter_stats(filters, FC_NORM, RP, RS)
244 minimum_order_comparison(FC_NORM, RP, RS)
245
246 # Apply to a test signal
247 t, x, outputs = apply_filters_to_signal(filters)
248
249 # -----------------------------------------------------------------------
250 # VISUALIZATION — 3×2 grid
251 # -----------------------------------------------------------------------
252 colors = {'Butterworth': 'tab:blue',
253 'Chebyshev I': 'tab:orange',
254 'Chebyshev II': 'tab:green',
255 'Elliptic': 'tab:red'}
256 styles = {'Butterworth': '-',
257 'Chebyshev I': '--',
258 'Chebyshev II': '-.',
259 'Elliptic': ':'}
260
261 fig, axes = plt.subplots(3, 2, figsize=(13, 12))
262 fig.suptitle(f'IIR Filter Design Comparison (N={ORDER}, fc={FC} Hz, fs={FS} Hz)',
263 fontsize=13, fontweight='bold')
264
265 # (0,0) Magnitude response (linear)
266 ax = axes[0, 0]
267 for name, (w, H) in responses.items():
268 ax.plot(w / np.pi * (FS / 2), np.abs(H),
269 color=colors[name], ls=styles[name], lw=2, label=name)
270 ax.axvline(FC, color='k', ls=':', lw=1, label=f'fc={FC} Hz')
271 ax.axhline(1.0 / np.sqrt(2), color='gray', ls='--', lw=1, alpha=0.6, label='-3 dB')
272 ax.set_xlabel('Frequency (Hz)')
273 ax.set_ylabel('|H(f)|')
274 ax.set_title('Magnitude Response (linear)')
275 ax.legend(fontsize=8)
276 ax.grid(True, alpha=0.3)
277 ax.set_xlim(0, FS / 2)
278
279 # (0,1) Magnitude response (dB)
280 ax = axes[0, 1]
281 for name, (w, H) in responses.items():
282 mag_db = 20 * np.log10(np.maximum(np.abs(H), 1e-15))
283 ax.plot(w / np.pi * (FS / 2), mag_db,
284 color=colors[name], ls=styles[name], lw=2, label=name)
285 ax.axvline(FC, color='k', ls=':', lw=1)
286 ax.axhline(-RP, color='gray', ls='--', lw=1, alpha=0.6, label=f'-{RP} dB (Rp)')
287 ax.axhline(-RS, color='gray', ls='-.', lw=1, alpha=0.6, label=f'-{RS} dB (Rs)')
288 ax.set_xlabel('Frequency (Hz)')
289 ax.set_ylabel('Magnitude (dB)')
290 ax.set_title('Magnitude Response (dB)')
291 ax.legend(fontsize=8)
292 ax.grid(True, alpha=0.3)
293 ax.set_xlim(0, FS / 2)
294 ax.set_ylim(-60, 5)
295
296 # (1,0) Phase response
297 ax = axes[1, 0]
298 for name, (w, H) in responses.items():
299 phase = np.unwrap(np.angle(H))
300 ax.plot(w / np.pi * (FS / 2), np.degrees(phase),
301 color=colors[name], ls=styles[name], lw=2, label=name)
302 ax.axvline(FC, color='k', ls=':', lw=1)
303 ax.set_xlabel('Frequency (Hz)')
304 ax.set_ylabel('Phase (degrees)')
305 ax.set_title('Phase Response (unwrapped)')
306 ax.legend(fontsize=8)
307 ax.grid(True, alpha=0.3)
308 ax.set_xlim(0, FS / 2)
309
310 # (1,1) Pole-zero plots (2×2 sub-grid within this subplot area)
311 # Use inset axes for each filter's pole-zero diagram
312 ax = axes[1, 1]
313 ax.set_visible(False) # hide the parent; draw four insets manually
314
315 filter_list = list(filters.items())
316 inset_positions = [(0.01, 0.51, 0.48, 0.46), # Butterworth (top-left)
317 (0.51, 0.51, 0.48, 0.46), # Cheby I (top-right)
318 (0.01, 0.01, 0.48, 0.46), # Cheby II (bottom-left)
319 (0.51, 0.01, 0.48, 0.46)] # Elliptic (bottom-right)
320
321 for (name, (b, a)), pos in zip(filter_list, inset_positions):
322 ax_in = fig.add_axes([
323 axes[1, 1].get_position().x0 + pos[0] * axes[1, 1].get_position().width,
324 axes[1, 1].get_position().y0 + pos[1] * axes[1, 1].get_position().height,
325 pos[2] * axes[1, 1].get_position().width,
326 pos[3] * axes[1, 1].get_position().height,
327 ])
328 zeros, poles, _ = signal.tf2zpk(b, a)
329 # Unit circle
330 theta = np.linspace(0, 2 * np.pi, 300)
331 ax_in.plot(np.cos(theta), np.sin(theta), 'k-', lw=0.8, alpha=0.4)
332 ax_in.plot(np.real(poles), np.imag(poles), 'rx', ms=7, mew=2, label='Poles')
333 ax_in.plot(np.real(zeros), np.imag(zeros), 'bo', ms=5, mfc='none', mew=1.5,
334 label='Zeros')
335 ax_in.axhline(0, color='gray', lw=0.5)
336 ax_in.axvline(0, color='gray', lw=0.5)
337 ax_in.set_title(name, fontsize=8)
338 ax_in.set_aspect('equal')
339 ax_in.set_xlim(-1.6, 1.6)
340 ax_in.set_ylim(-1.6, 1.6)
341 ax_in.tick_params(labelsize=6)
342 ax_in.grid(True, alpha=0.2)
343 if name == 'Butterworth':
344 ax_in.legend(fontsize=6, loc='lower right')
345
346 # Add a title for the pole-zero panel area
347 fig.text(
348 axes[1, 1].get_position().x0 + 0.5 * axes[1, 1].get_position().width,
349 axes[1, 1].get_position().y1 + 0.005,
350 'Pole-Zero Diagrams (unit circle shown)',
351 ha='center', va='bottom', fontsize=10, fontweight='bold'
352 )
353
354 # (2,0) Filtered signal — time domain (first 100 ms)
355 ax = axes[2, 0]
356 idx_end = int(0.1 * FS)
357 ax.plot(t[:idx_end] * 1e3, x[:idx_end], color='gray', lw=0.8,
358 alpha=0.5, label='Noisy input')
359 for name, y in outputs.items():
360 ax.plot(t[:idx_end] * 1e3, y[:idx_end],
361 color=colors[name], ls=styles[name], lw=1.8, label=name)
362 ax.set_xlabel('Time (ms)')
363 ax.set_ylabel('Amplitude')
364 ax.set_title('Filtered Signal (first 100 ms, zero-phase)')
365 ax.legend(fontsize=8)
366 ax.grid(True, alpha=0.3)
367
368 # (2,1) Spectrum of filtered outputs
369 ax = axes[2, 1]
370 N_sig = len(x)
371 freqs = np.fft.rfftfreq(N_sig, d=1.0 / FS)
372 X_mag = np.abs(np.fft.rfft(x)) / N_sig
373 ax.plot(freqs, 20 * np.log10(np.maximum(X_mag, 1e-15)),
374 color='gray', lw=1.0, alpha=0.5, label='Input')
375 for name, y in outputs.items():
376 Y_mag = np.abs(np.fft.rfft(y)) / N_sig
377 ax.plot(freqs, 20 * np.log10(np.maximum(Y_mag, 1e-15)),
378 color=colors[name], ls=styles[name], lw=1.8, label=name)
379 ax.axvline(FC, color='k', ls=':', lw=1)
380 ax.set_xlabel('Frequency (Hz)')
381 ax.set_ylabel('Magnitude (dBFS)')
382 ax.set_title('Output Spectrum')
383 ax.legend(fontsize=8)
384 ax.grid(True, alpha=0.3)
385 ax.set_xlim(0, FS / 2)
386 ax.set_ylim(-80, 10)
387
388 plt.tight_layout()
389 out_path = '/opt/projects/01_Personal/03_Study/examples/Signal_Processing/10_iir_design.png'
390 plt.savefig(out_path, dpi=150)
391 print(f"\nSaved visualization: {out_path}")
392 plt.close()
393
394 print("\n" + "=" * 70)
395 print("Key Takeaways:")
396 print(" - Butterworth: flattest passband, but needs highest order for sharp roll-off")
397 print(" - Chebyshev I: sharper than Butterworth, at cost of passband ripple")
398 print(" - Chebyshev II: all ripple pushed into stopband, flat passband")
399 print(" - Elliptic: sharpest roll-off for given order, ripple in both bands")
400 print(" - All IIR filters have nonlinear phase — use filtfilt for zero-phase")
401 print(" - Poles inside unit circle → stable; zeros on unit circle → deep nulls")
402 print("=" * 70)