1#!/usr/bin/env python3
2"""
3FIR Filter Design Methods
4==========================
5
6This script demonstrates three principal FIR filter design techniques:
7
81. Window Method — Hamming window (lowpass)
92. Window Method — Kaiser window (lowpass, tunable stopband attenuation)
103. Window Method — Hamming window (bandpass)
114. Parks-McClellan (Remez) optimal equiripple design (lowpass)
125. Comparison of window method vs optimal design
136. Applying the designed filters to a test signal
14
15Key Concepts:
16 - Ideal filter → multiply by window → practical FIR filter
17 - Window choice controls sidelobe level vs transition width trade-off
18 - Parks-McClellan minimises the maximum error (Chebyshev criterion)
19 - Kaiser window: adjustable β controls sidelobe attenuation
20
21Design relationships (window method):
22 - Filter order N determined by transition width and window type
23 - Hamming: stopband attenuation ≈ 53 dB
24 - Kaiser: As = 2.285 * (N-1) * Δω + 7.95 (Harris formula)
25
26Author: Educational example for Signal Processing
27License: MIT
28"""
29
30import numpy as np
31from scipy import signal
32import matplotlib.pyplot as plt
33
34# ============================================================================
35# DESIGN HELPERS
36# ============================================================================
37
38def kaiser_parameters(As_dB, delta_omega):
39 """
40 Estimate Kaiser window β and filter order from design specifications.
41
42 Using the Kaiser empirical formulas (Harris 1978):
43 β = 0 if As < 21 dB
44 β = 0.5842*(As-21)^0.4 + 0.07886*(As-21) if 21 ≤ As ≤ 50 dB
45 β = 0.1102*(As-8.7) if As > 50 dB
46
47 N ≈ (As - 7.95) / (2.285 * Δω) + 1 (Δω in rad/sample)
48
49 Args:
50 As_dB (float): Desired stopband attenuation in dB (positive value)
51 delta_omega (float): Transition width in rad/sample
52
53 Returns:
54 tuple: (N, beta) — filter order and Kaiser β parameter
55 """
56 if As_dB < 21:
57 beta = 0.0
58 elif As_dB <= 50:
59 beta = 0.5842 * (As_dB - 21) ** 0.4 + 0.07886 * (As_dB - 21)
60 else:
61 beta = 0.1102 * (As_dB - 8.7)
62
63 N = int(np.ceil((As_dB - 7.95) / (2.285 * delta_omega))) + 1
64 if N % 2 == 0:
65 N += 1 # ensure odd length for symmetric filter (type I)
66
67 return N, beta
68
69
70# ============================================================================
71# DESIGN FUNCTIONS
72# ============================================================================
73
74def design_lowpass_hamming(cutoff_norm, N):
75 """
76 Design FIR lowpass filter using the Hamming window method.
77
78 Steps:
79 1. Compute ideal sinc-based impulse response h_ideal[n]
80 2. Multiply by Hamming window w[n]
81 3. Result is a linear-phase FIR filter
82
83 Args:
84 cutoff_norm (float): Normalized cutoff frequency (0 to 1, where 1 = Nyquist)
85 N (int): Filter order (number of taps = N + 1, should be even for odd length)
86
87 Returns:
88 ndarray: FIR filter coefficients (length N+1)
89 """
90 # scipy.signal.firwin uses half-band convention: cutoff in cycles/sample
91 # fs=2 maps cutoff_norm in [0,1] to Hz, where 1 corresponds to Nyquist
92 b = signal.firwin(N + 1, cutoff_norm, window='hamming')
93 return b
94
95
96def design_lowpass_kaiser(cutoff_norm, As_dB, delta_norm):
97 """
98 Design FIR lowpass filter using the Kaiser window method.
99
100 The Kaiser window allows direct specification of stopband attenuation
101 and automatically computes the required filter order.
102
103 Args:
104 cutoff_norm (float): Normalized cutoff (0 to 1, where 1 = Nyquist)
105 As_dB (float): Minimum stopband attenuation in dB
106 delta_norm (float): Transition band width, normalized (0 to 1)
107
108 Returns:
109 tuple: (b, N, beta) — coefficients, actual order, and Kaiser β
110 """
111 delta_omega = delta_norm * np.pi # convert to rad/sample
112 N, beta = kaiser_parameters(As_dB, delta_omega)
113 b = signal.firwin(N, cutoff_norm, window=('kaiser', beta))
114 return b, N, beta
115
116
117def design_bandpass_hamming(low_norm, high_norm, N):
118 """
119 Design FIR bandpass filter using the Hamming window method.
120
121 Args:
122 low_norm (float): Lower normalized cutoff (0 to 1)
123 high_norm (float): Upper normalized cutoff (0 to 1)
124 N (int): Filter order
125
126 Returns:
127 ndarray: FIR filter coefficients (length N+1)
128 """
129 # firwin with two cutoffs and pass_zero=False → bandpass
130 b = signal.firwin(N + 1, [low_norm, high_norm], pass_zero=False, window='hamming')
131 return b
132
133
134def design_optimal_remez(cutoff_norm, trans_norm, As_dB):
135 """
136 Design optimal FIR lowpass filter using Parks-McClellan (Remez) algorithm.
137
138 The Remez algorithm minimises the maximum weighted error between the
139 desired and actual frequency response — equiripple in both bands.
140
141 Estimate filter order using Kaiser's order formula as a starting point.
142
143 Args:
144 cutoff_norm (float): Normalized passband edge (0 to 1)
145 trans_norm (float): Normalized transition band width (0 to 1)
146 As_dB (float): Approximate stopband attenuation in dB
147
148 Returns:
149 tuple: (b, N) — filter coefficients and filter order
150 """
151 delta_omega = trans_norm * np.pi
152 N, _ = kaiser_parameters(As_dB, delta_omega)
153
154 # Define band edges: [0, f_pass, f_stop, 1] in normalized units (0..1)
155 f_pass = cutoff_norm - trans_norm / 2
156 f_stop = cutoff_norm + trans_norm / 2
157
158 # Clamp to valid range
159 f_pass = max(f_pass, 0.01)
160 f_stop = min(f_stop, 0.99)
161
162 bands = [0, f_pass, f_stop, 1.0]
163 desired = [1, 0] # passband gain = 1, stopband gain = 0
164
165 # Convert passband ripple / stopband attenuation to relative weights
166 # δ_p ≈ δ_s for equal-ripple; weight = 1 everywhere for simplicity
167 b = signal.remez(N, bands, desired, fs=2.0)
168 return b, N
169
170
171# ============================================================================
172# MAIN DEMONSTRATIONS
173# ============================================================================
174
175def demo_lowpass_designs():
176 """Design and compare lowpass FIR filters."""
177 print("=" * 65)
178 print("1. LOWPASS FIR FILTER DESIGNS")
179 print("=" * 65)
180
181 cutoff = 0.3 # normalized cutoff (fraction of Nyquist)
182 N_hamming = 50 # fixed order for Hamming
183
184 # Hamming window design
185 b_hamming = design_lowpass_hamming(cutoff, N_hamming)
186 w, H_hamming = signal.freqz(b_hamming, worN=1024)
187
188 print(f"\nHamming Window Lowpass (N={N_hamming}, cutoff={cutoff}×Nyquist)")
189 # Find -3 dB point
190 mag = np.abs(H_hamming)
191 idx_3dB = np.argmin(np.abs(mag - (1 / np.sqrt(2))))
192 print(f" Taps : {len(b_hamming)}")
193 print(f" -3 dB : {w[idx_3dB] / np.pi:.3f} × π rad/sample")
194 # Stopband attenuation (at 0.4 Nyquist and above)
195 idx_stop = np.where(w / np.pi > cutoff + 0.1)[0][0]
196 stop_atten = -20 * np.log10(np.maximum(np.max(mag[idx_stop:]), 1e-12))
197 print(f" Stopband attenuation: ≥ {stop_atten:.1f} dB")
198
199 # Kaiser window design
200 As_dB = 60.0
201 delta_norm = 0.08
202 b_kaiser, N_kaiser, beta = design_lowpass_kaiser(cutoff, As_dB, delta_norm)
203 w, H_kaiser = signal.freqz(b_kaiser, worN=1024)
204
205 print(f"\nKaiser Window Lowpass (As={As_dB} dB, Δf={delta_norm}×Nyquist)")
206 print(f" Computed N : {N_kaiser}")
207 print(f" Kaiser β : {beta:.4f}")
208 print(f" Taps : {len(b_kaiser)}")
209 mag_k = np.abs(H_kaiser)
210 stop_atten_k = -20 * np.log10(np.maximum(np.max(mag_k[idx_stop:]), 1e-12))
211 print(f" Stopband attenuation: ≥ {stop_atten_k:.1f} dB")
212
213 return b_hamming, b_kaiser, cutoff
214
215
216def demo_bandpass_design():
217 """Design a bandpass FIR filter."""
218 print("\n" + "=" * 65)
219 print("2. BANDPASS FIR FILTER (Hamming Window)")
220 print("=" * 65)
221
222 low_norm = 0.2 # lower cutoff (fraction of Nyquist)
223 high_norm = 0.5 # upper cutoff
224 N = 60
225
226 b_bp = design_bandpass_hamming(low_norm, high_norm, N)
227 w, H_bp = signal.freqz(b_bp, worN=1024)
228
229 print(f"\nBandpass: passband {low_norm}–{high_norm} × Nyquist, N={N}")
230 print(f" Taps: {len(b_bp)}")
231
232 # Report gain at passband center and stopband edges
233 center = (low_norm + high_norm) / 2
234 mag_bp = np.abs(H_bp)
235 idx_center = np.argmin(np.abs(w / np.pi - center))
236 print(f" Gain at center ({center:.2f}×π): {20*np.log10(mag_bp[idx_center]):.2f} dB")
237
238 idx_stop_lo = np.argmin(np.abs(w / np.pi - 0.05))
239 idx_stop_hi = np.argmin(np.abs(w / np.pi - 0.85))
240 print(f" Gain at lower stopband (0.05×π): {20*np.log10(max(mag_bp[idx_stop_lo], 1e-12)):.1f} dB")
241 print(f" Gain at upper stopband (0.85×π): {20*np.log10(max(mag_bp[idx_stop_hi], 1e-12)):.1f} dB")
242
243 return b_bp, low_norm, high_norm
244
245
246def demo_remez_vs_window(cutoff, b_hamming):
247 """Compare Parks-McClellan optimal design with Hamming window method."""
248 print("\n" + "=" * 65)
249 print("3. PARKS-McCLELLAN (REMEZ) vs WINDOW METHOD")
250 print("=" * 65)
251
252 As_dB = 53.0 # target matching Hamming's ~53 dB
253 trans_norm = 0.08
254
255 b_remez, N_remez = design_optimal_remez(cutoff, trans_norm, As_dB)
256 w, H_remez = signal.freqz(b_remez, worN=1024)
257 _, H_hamming = signal.freqz(b_hamming, worN=1024)
258
259 mag_r = np.abs(H_remez)
260 mag_h = np.abs(H_hamming)
261
262 # Passband ripple
263 idx_pb_end = np.argmin(np.abs(w / np.pi - (cutoff - trans_norm / 2)))
264 pb_ripple_r = 20 * np.log10(np.max(mag_r[:idx_pb_end]) / np.min(mag_r[1:idx_pb_end]))
265 pb_ripple_h = 20 * np.log10(np.max(mag_h[:idx_pb_end]) / np.min(mag_h[1:idx_pb_end]))
266
267 # Stopband attenuation
268 idx_sb_start = np.argmin(np.abs(w / np.pi - (cutoff + trans_norm / 2)))
269 sb_atten_r = -20 * np.log10(np.maximum(np.max(mag_r[idx_sb_start:]), 1e-12))
270 sb_atten_h = -20 * np.log10(np.maximum(np.max(mag_h[idx_sb_start:]), 1e-12))
271
272 print(f"\nCutoff = {cutoff}×Nyquist, target As = {As_dB} dB")
273 print(f"\n{'Method':<20} {'Taps':>6} {'PB Ripple':>12} {'SB Atten':>12}")
274 print("-" * 54)
275 print(f"{'Hamming Window':<20} {len(b_hamming):>6} {pb_ripple_h:>11.2f} dB {sb_atten_h:>10.1f} dB")
276 print(f"{'Remez (optimal)':<20} {len(b_remez):>6} {pb_ripple_r:>11.2f} dB {sb_atten_r:>10.1f} dB")
277 print(f"\nRemez achieves equiripple in both bands (Chebyshev criterion).")
278 print(f"Window method is simpler but less flexible in band control.")
279
280 return b_remez, w, H_remez, H_hamming
281
282
283def demo_apply_to_signal(b_hamming, b_kaiser, b_bp):
284 """Apply designed filters to a multi-tone test signal."""
285 print("\n" + "=" * 65)
286 print("4. APPLYING FILTERS TO A TEST SIGNAL")
287 print("=" * 65)
288
289 # Test signal: sum of sinusoids at 5, 15, 30, and 45 Hz
290 # Sampling at 100 Hz → Nyquist = 50 Hz
291 fs = 100.0
292 t = np.linspace(0, 2.0, int(fs * 2), endpoint=False)
293 f1, f2, f3, f4 = 5, 15, 30, 45 # Hz
294 # Normalized: f1/Nyq=0.1, f2=0.3, f3=0.6, f4=0.9
295 x = (np.sin(2 * np.pi * f1 * t) +
296 np.sin(2 * np.pi * f2 * t) +
297 np.sin(2 * np.pi * f3 * t) +
298 np.sin(2 * np.pi * f4 * t))
299
300 # Lowpass (Hamming, cutoff=0.3 → 15 Hz): keep f1, f2; reject f3, f4
301 y_lp = signal.lfilter(b_hamming, 1, x)
302
303 # Lowpass (Kaiser, cutoff=0.3, 60 dB): same goal, better attenuation
304 y_kaiser = signal.lfilter(b_kaiser, 1, x)
305
306 # Bandpass (0.2–0.5 Nyquist → 10–25 Hz): keep f2; reject f1, f3, f4
307 y_bp = signal.lfilter(b_bp, 1, x)
308
309 print(f"\nTest signal: {f1}, {f2}, {f3}, {f4} Hz (fs={fs} Hz)")
310 print(f" Normalized: {f1/50:.2f}, {f2/50:.2f}, {f3/50:.2f}, {f4/50:.2f} × Nyquist")
311 print(f"\nLowpass Hamming (cutoff=0.3×Nyq = {0.3*50:.0f} Hz):")
312 print(f" Expected: pass {f1}+{f2} Hz, reject {f3}+{f4} Hz")
313 print(f"\nKaiser (same cutoff, 60 dB stopband):")
314 print(f" Better stopband attenuation than Hamming")
315 print(f"\nBandpass (0.2–0.5 × Nyq = {0.2*50:.0f}–{0.5*50:.0f} Hz):")
316 print(f" Expected: pass {f2} Hz, reject {f1}+{f3}+{f4} Hz")
317
318 return t, x, y_lp, y_kaiser, y_bp
319
320
321# ============================================================================
322# ENTRY POINT
323# ============================================================================
324
325if __name__ == "__main__":
326 print("FIR FILTER DESIGN METHODS")
327 print("=" * 65)
328
329 b_hamming, b_kaiser, cutoff = demo_lowpass_designs()
330 b_bp, low_norm, high_norm = demo_bandpass_design()
331 b_remez, w_cmp, H_remez, H_hamming_cmp = demo_remez_vs_window(cutoff, b_hamming)
332 t, x, y_lp, y_kaiser, y_bp = demo_apply_to_signal(b_hamming, b_kaiser, b_bp)
333
334 # -----------------------------------------------------------------------
335 # VISUALIZATION — 3×2 grid
336 # -----------------------------------------------------------------------
337 fig, axes = plt.subplots(3, 2, figsize=(13, 12))
338 fig.suptitle('FIR Filter Design Methods', fontsize=14, fontweight='bold')
339
340 w_h, H_h = signal.freqz(b_hamming, worN=1024)
341 w_k, H_k = signal.freqz(b_kaiser, worN=1024)
342 w_bp, H_bp = signal.freqz(b_bp, worN=1024)
343
344 # (0,0) Lowpass magnitude responses (Hamming vs Kaiser)
345 axes[0, 0].plot(w_h / np.pi, 20 * np.log10(np.maximum(np.abs(H_h), 1e-12)),
346 'b-', linewidth=2, label=f'Hamming (N={len(b_hamming)-1})')
347 axes[0, 0].plot(w_k / np.pi, 20 * np.log10(np.maximum(np.abs(H_k), 1e-12)),
348 'r--', linewidth=2, label=f'Kaiser (N={len(b_kaiser)-1}, β)')
349 axes[0, 0].axvline(cutoff, color='k', linestyle=':', linewidth=1, label=f'Cutoff {cutoff}×π')
350 axes[0, 0].set_xlabel('Normalized Frequency (×π rad/sample)')
351 axes[0, 0].set_ylabel('Magnitude (dB)')
352 axes[0, 0].set_title('Lowpass: Hamming vs Kaiser Window')
353 axes[0, 0].legend(fontsize=9)
354 axes[0, 0].grid(True, alpha=0.3)
355 axes[0, 0].set_xlim(0, 1)
356 axes[0, 0].set_ylim(-90, 5)
357
358 # (0,1) Lowpass impulse responses
359 axes[0, 1].plot(b_hamming, 'b-o', markersize=3, linewidth=1.5, label='Hamming')
360 axes[0, 1].plot(np.linspace(0, len(b_hamming) - 1, len(b_kaiser)),
361 b_kaiser, 'r--^', markersize=3, linewidth=1.5, label='Kaiser')
362 axes[0, 1].set_xlabel('Tap index n')
363 axes[0, 1].set_ylabel('Coefficient value')
364 axes[0, 1].set_title('Impulse Response (Filter Coefficients)')
365 axes[0, 1].legend()
366 axes[0, 1].grid(True, alpha=0.3)
367
368 # (1,0) Bandpass magnitude response
369 axes[1, 0].plot(w_bp / np.pi, 20 * np.log10(np.maximum(np.abs(H_bp), 1e-12)),
370 'g-', linewidth=2)
371 axes[1, 0].axvspan(low_norm, high_norm, alpha=0.15, color='green', label='Passband')
372 axes[1, 0].set_xlabel('Normalized Frequency (×π rad/sample)')
373 axes[1, 0].set_ylabel('Magnitude (dB)')
374 axes[1, 0].set_title(f'Bandpass: {low_norm}–{high_norm} × Nyquist (Hamming)')
375 axes[1, 0].legend()
376 axes[1, 0].grid(True, alpha=0.3)
377 axes[1, 0].set_xlim(0, 1)
378 axes[1, 0].set_ylim(-90, 5)
379
380 # (1,1) Remez vs Window comparison
381 axes[1, 1].plot(w_cmp / np.pi, 20 * np.log10(np.maximum(np.abs(H_hamming_cmp), 1e-12)),
382 'b-', linewidth=2, label=f'Hamming (N={len(b_hamming)-1})')
383 axes[1, 1].plot(w_cmp / np.pi, 20 * np.log10(np.maximum(np.abs(H_remez), 1e-12)),
384 'm--', linewidth=2, label=f'Remez/PM (N={len(b_remez)-1})')
385 axes[1, 1].axvline(cutoff, color='k', linestyle=':', linewidth=1)
386 axes[1, 1].set_xlabel('Normalized Frequency (×π rad/sample)')
387 axes[1, 1].set_ylabel('Magnitude (dB)')
388 axes[1, 1].set_title('Window Method vs Parks-McClellan (Remez)')
389 axes[1, 1].legend()
390 axes[1, 1].grid(True, alpha=0.3)
391 axes[1, 1].set_xlim(0, 1)
392 axes[1, 1].set_ylim(-90, 5)
393
394 # (2,0) Test signal time domain — first 0.5 s
395 idx = int(0.5 * 100) # first 0.5 s at fs=100 Hz
396 axes[2, 0].plot(t[:idx], x[:idx], 'gray', linewidth=0.8, alpha=0.5, label='Input (4 tones)')
397 axes[2, 0].plot(t[:idx], y_lp[:idx], 'b-', linewidth=2, label='Lowpass (Hamming)')
398 axes[2, 0].plot(t[:idx], y_bp[:idx], 'g--', linewidth=2, label='Bandpass')
399 axes[2, 0].set_xlabel('Time (s)')
400 axes[2, 0].set_ylabel('Amplitude')
401 axes[2, 0].set_title('Filtered Test Signal (time domain)')
402 axes[2, 0].legend(fontsize=9)
403 axes[2, 0].grid(True, alpha=0.3)
404
405 # (2,1) Spectrum of input vs filtered signals
406 N_fft = len(x)
407 freqs = np.fft.rfftfreq(N_fft, d=1.0 / 100.0)
408 X_mag = np.abs(np.fft.rfft(x)) / N_fft
409 YLP_mag = np.abs(np.fft.rfft(y_lp)) / N_fft
410 YBP_mag = np.abs(np.fft.rfft(y_bp)) / N_fft
411
412 axes[2, 1].plot(freqs, X_mag, 'gray', linewidth=1.5, alpha=0.6, label='Input spectrum')
413 axes[2, 1].plot(freqs, YLP_mag, 'b-', linewidth=2, label='Lowpass out')
414 axes[2, 1].plot(freqs, YBP_mag, 'g--', linewidth=2, label='Bandpass out')
415 axes[2, 1].set_xlabel('Frequency (Hz)')
416 axes[2, 1].set_ylabel('Magnitude')
417 axes[2, 1].set_title('Filtered Test Signal (frequency domain)')
418 axes[2, 1].legend(fontsize=9)
419 axes[2, 1].grid(True, alpha=0.3)
420 axes[2, 1].set_xlim(0, 50)
421
422 plt.tight_layout()
423 out_path = '/opt/projects/01_Personal/03_Study/examples/Signal_Processing/09_fir_design.png'
424 plt.savefig(out_path, dpi=150)
425 print(f"\nSaved visualization: {out_path}")
426 plt.close()
427
428 print("\n" + "=" * 65)
429 print("Key Takeaways:")
430 print(" - Window method: simple, fixed sidelobe shape per window type")
431 print(" - Hamming: ~53 dB stopband; Kaiser: adjustable via β parameter")
432 print(" - Parks-McClellan: optimal equiripple, minimum order for given spec")
433 print(" - All FIR methods produce linear-phase (symmetric) filters")
434 print(" - Higher stopband attenuation requires more taps (longer delay)")
435 print("=" * 65)