1#!/usr/bin/env python3
2"""
3Digital Filter Fundamentals: FIR vs IIR
4========================================
5
6This script compares Finite Impulse Response (FIR) and Infinite Impulse
7Response (IIR) filters, covering:
8
9- Magnitude response, phase response, and group delay
10- Simple FIR filter: moving average (lowpass)
11- Simple IIR filter: 1st-order recursive (exponential smoother)
12- Impulse response and step response comparison
13- Applying both filters to a noisy signal
14
15Key Concepts:
16 FIR: y[n] = sum_{k=0}^{M} b_k * x[n-k] (always stable)
17 IIR: y[n] = sum b_k x[n-k] - sum a_k y[n-k] (can be unstable)
18
19 Group delay: τ(ω) = -dφ(ω)/dω
20 FIR with symmetric coefficients → linear phase → constant group delay
21
22Author: Educational example for Signal Processing
23License: MIT
24"""
25
26import numpy as np
27from scipy import signal
28import matplotlib.pyplot as plt
29
30# ============================================================================
31# FILTER DEFINITIONS
32# ============================================================================
33
34def moving_average_fir(M):
35 """
36 Moving average FIR filter of length M.
37
38 H(z) = (1/M) * (1 + z^{-1} + ... + z^{-(M-1)})
39 Frequency response: H(e^{jω}) = (1/M) * sin(Mω/2) / sin(ω/2) * e^{-j(M-1)ω/2}
40
41 This is a symmetric (linear-phase) FIR lowpass filter.
42 The cutoff frequency is approximately 0.9/M * π rad/sample.
43
44 Args:
45 M (int): Filter length (number of taps)
46
47 Returns:
48 ndarray: Filter coefficients b (length M), a = [1]
49 """
50 b = np.ones(M) / M
51 a = np.array([1.0])
52 return b, a
53
54
55def first_order_iir(alpha):
56 """
57 1st-order IIR recursive (exponential smoothing) filter.
58
59 y[n] = α * x[n] + (1 - α) * y[n-1]
60
61 Transfer function: H(z) = α / (1 - (1-α) z^{-1})
62 Pole at z = (1 - α) — stable for 0 < α < 2.
63
64 α close to 0 → heavy smoothing (narrow bandwidth)
65 α close to 1 → minimal smoothing (wide bandwidth)
66
67 Args:
68 alpha (float): Smoothing factor, 0 < alpha < 1
69
70 Returns:
71 tuple: (b, a) coefficient arrays
72 """
73 b = np.array([alpha])
74 a = np.array([1.0, -(1.0 - alpha)])
75 return b, a
76
77
78# ============================================================================
79# ANALYSIS FUNCTIONS
80# ============================================================================
81
82def compute_group_delay(b, a, worN=512):
83 """
84 Compute group delay τ(ω) = -dφ(ω)/dω using scipy.signal.group_delay.
85
86 Args:
87 b, a: Filter coefficients
88 worN (int): Number of frequency points
89
90 Returns:
91 tuple: (frequencies in rad/sample, group delay in samples)
92 """
93 w, gd = signal.group_delay((b, a), w=worN)
94 return w, gd
95
96
97def impulse_response(b, a, n_samples=60):
98 """
99 Compute the filter's impulse response h[n] by filtering δ[n].
100
101 Args:
102 b, a: Filter coefficients
103 n_samples (int): Number of output samples
104
105 Returns:
106 ndarray: h[n] for n = 0, 1, ..., n_samples-1
107 """
108 delta = np.zeros(n_samples)
109 delta[0] = 1.0
110 return signal.lfilter(b, a, delta)
111
112
113def step_response(b, a, n_samples=60):
114 """
115 Compute the filter's step response s[n] by filtering u[n].
116
117 Args:
118 b, a: Filter coefficients
119 n_samples (int): Number of output samples
120
121 Returns:
122 ndarray: s[n] for n = 0, 1, ..., n_samples-1
123 """
124 step = np.ones(n_samples)
125 return signal.lfilter(b, a, step)
126
127
128# ============================================================================
129# MAIN DEMONSTRATIONS
130# ============================================================================
131
132def demo_frequency_responses(M=15, alpha=0.2):
133 """Compare FIR and IIR frequency responses."""
134 print("=" * 65)
135 print("1. FREQUENCY RESPONSE COMPARISON: FIR vs IIR")
136 print("=" * 65)
137
138 b_fir, a_fir = moving_average_fir(M)
139 b_iir, a_iir = first_order_iir(alpha)
140
141 w, H_fir = signal.freqz(b_fir, a_fir, worN=512)
142 _, H_iir = signal.freqz(b_iir, a_iir, worN=512)
143
144 _, gd_fir = compute_group_delay(b_fir, a_fir)
145 _, gd_iir = compute_group_delay(b_iir, a_iir)
146
147 print(f"\nFIR Moving Average (M={M} taps)")
148 print(f" DC gain (ω=0): {np.abs(H_fir[0]):.4f} (should be 1.0)")
149 print(f" Half-power freq: ~{0.9 / M:.3f} × π rad/sample")
150 print(f" Group delay: constant = {np.mean(gd_fir[:50]):.1f} samples (linear phase)")
151
152 print(f"\nIIR Exponential Smoother (α={alpha})")
153 print(f" DC gain (ω=0): {np.abs(H_iir[0]):.4f} (should be 1.0)")
154 print(f" Pole location: z = {1 - alpha:.2f} (inside unit circle → stable)")
155 print(f" Group delay at ω=0: {gd_iir[0]:.2f} samples (non-linear phase)")
156
157 return (w, H_fir, H_iir, gd_fir, gd_iir), (b_fir, a_fir), (b_iir, a_iir)
158
159
160def demo_impulse_step_response(b_fir, a_fir, b_iir, a_iir):
161 """Show impulse and step responses for both filters."""
162 print("\n" + "=" * 65)
163 print("2. IMPULSE AND STEP RESPONSES")
164 print("=" * 65)
165
166 h_fir = impulse_response(b_fir, a_fir, n_samples=50)
167 h_iir = impulse_response(b_iir, a_iir, n_samples=50)
168 s_fir = step_response(b_fir, a_fir, n_samples=50)
169 s_iir = step_response(b_iir, a_iir, n_samples=50)
170
171 print(f"\nFIR impulse response: finite support, {np.sum(h_fir != 0)} nonzero samples")
172 print(f" Final value h[49] = {h_fir[49]:.6f} (exactly zero after M taps)")
173 print(f" Steady-state step response: {s_fir[49]:.4f} (DC gain = 1)")
174
175 alpha = b_iir[0]
176 print(f"\nIIR impulse response: infinite support (exponentially decaying)")
177 print(f" h[0] = {h_iir[0]:.4f} (= α = {alpha})")
178 print(f" h[10] = {h_iir[10]:.6f} (= α * (1-α)^10 = {alpha * (1-alpha)**10:.6f})")
179 print(f" h[49] = {h_iir[49]:.2e} (never exactly zero)")
180 print(f" Steady-state step response: {s_iir[49]:.4f} (DC gain = 1)")
181
182 return h_fir, h_iir, s_fir, s_iir
183
184
185def demo_noise_filtering(b_fir, a_fir, b_iir, a_iir):
186 """Apply both filters to a noisy signal and compare."""
187 print("\n" + "=" * 65)
188 print("3. FILTERING A NOISY SIGNAL")
189 print("=" * 65)
190
191 # Generate test signal: slow 5 Hz sine + high-frequency noise
192 fs = 1000 # sampling frequency (Hz)
193 t = np.linspace(0, 0.5, int(fs * 0.5), endpoint=False)
194 clean = np.sin(2 * np.pi * 5 * t) # 5 Hz signal
195 np.random.seed(0)
196 noise = 0.5 * np.random.randn(len(t)) # broadband noise
197 noisy = clean + noise
198
199 y_fir = signal.lfilter(b_fir, a_fir, noisy)
200 y_iir = signal.lfilter(b_iir, a_iir, noisy)
201
202 # SNR calculation (over second half to avoid transient)
203 half = len(t) // 2
204 snr_orig = 10 * np.log10(np.var(clean[half:]) / np.var(noise[half:]))
205 snr_fir = 10 * np.log10(np.var(clean[half:]) / np.var(y_fir[half:] - clean[half:]))
206 snr_iir = 10 * np.log10(np.var(clean[half:]) / np.var(y_iir[half:] - clean[half:]))
207
208 print(f"\nTest signal: 5 Hz sine + Gaussian noise (SNR = {snr_orig:.1f} dB)")
209 print(f" FIR output SNR : {snr_fir:.1f} dB (improvement: {snr_fir - snr_orig:.1f} dB)")
210 print(f" IIR output SNR : {snr_iir:.1f} dB (improvement: {snr_iir - snr_orig:.1f} dB)")
211
212 M_fir = len(b_fir)
213 print(f"\nFIR introduces a constant delay of {(M_fir - 1) // 2} samples (linear phase)")
214 print(f"IIR introduces a signal-dependent (non-linear) phase distortion")
215
216 return t, noisy, clean, y_fir, y_iir
217
218
219# ============================================================================
220# ENTRY POINT
221# ============================================================================
222
223if __name__ == "__main__":
224 print("DIGITAL FILTER FUNDAMENTALS: FIR vs IIR")
225 print("=" * 65)
226
227 M = 15 # FIR moving average taps
228 alpha = 0.2 # IIR smoothing factor (smaller = more smoothing)
229
230 freq_data, (b_fir, a_fir), (b_iir, a_iir) = demo_frequency_responses(M, alpha)
231 w, H_fir, H_iir, gd_fir, gd_iir = freq_data
232
233 h_fir, h_iir, s_fir, s_iir = demo_impulse_step_response(b_fir, a_fir, b_iir, a_iir)
234 t, noisy, clean, y_fir, y_iir = demo_noise_filtering(b_fir, a_fir, b_iir, a_iir)
235
236 # -----------------------------------------------------------------------
237 # VISUALIZATION — 3×2 grid
238 # -----------------------------------------------------------------------
239 fig, axes = plt.subplots(3, 2, figsize=(13, 12))
240 fig.suptitle('Digital Filter Fundamentals: FIR vs IIR', fontsize=14, fontweight='bold')
241
242 norm_freq = w / np.pi # normalized frequency 0..1
243
244 # (0,0) Magnitude response
245 axes[0, 0].plot(norm_freq, 20 * np.log10(np.maximum(np.abs(H_fir), 1e-12)),
246 'b-', linewidth=2, label=f'FIR (M={M})')
247 axes[0, 0].plot(norm_freq, 20 * np.log10(np.maximum(np.abs(H_iir), 1e-12)),
248 'r--', linewidth=2, label=f'IIR (α={alpha})')
249 axes[0, 0].set_xlabel('Normalized Frequency (×π rad/sample)')
250 axes[0, 0].set_ylabel('Magnitude (dB)')
251 axes[0, 0].set_title('Magnitude Response')
252 axes[0, 0].legend()
253 axes[0, 0].grid(True, alpha=0.3)
254 axes[0, 0].set_xlim(0, 1)
255 axes[0, 0].set_ylim(-50, 5)
256
257 # (0,1) Phase response
258 phase_fir = np.unwrap(np.angle(H_fir))
259 phase_iir = np.unwrap(np.angle(H_iir))
260 axes[0, 1].plot(norm_freq, phase_fir, 'b-', linewidth=2, label=f'FIR (linear phase)')
261 axes[0, 1].plot(norm_freq, phase_iir, 'r--', linewidth=2, label='IIR (non-linear phase)')
262 axes[0, 1].set_xlabel('Normalized Frequency (×π rad/sample)')
263 axes[0, 1].set_ylabel('Phase (radians)')
264 axes[0, 1].set_title('Phase Response')
265 axes[0, 1].legend()
266 axes[0, 1].grid(True, alpha=0.3)
267 axes[0, 1].set_xlim(0, 1)
268
269 # (1,0) Group delay
270 axes[1, 0].plot(norm_freq, gd_fir, 'b-', linewidth=2, label='FIR (constant)')
271 axes[1, 0].plot(norm_freq, np.clip(gd_iir, -2, 30), 'r--', linewidth=2,
272 label='IIR (frequency-dependent)')
273 axes[1, 0].set_xlabel('Normalized Frequency (×π rad/sample)')
274 axes[1, 0].set_ylabel('Group Delay (samples)')
275 axes[1, 0].set_title('Group Delay')
276 axes[1, 0].legend()
277 axes[1, 0].grid(True, alpha=0.3)
278 axes[1, 0].set_xlim(0, 1)
279
280 # (1,1) Impulse response
281 n = np.arange(len(h_fir))
282 axes[1, 1].stem(n, h_fir, basefmt='k-', linefmt='C0-', markerfmt='C0o',
283 label=f'FIR h[n]')
284 axes[1, 1].stem(n, h_iir, basefmt='k-', linefmt='C1--', markerfmt='C1^',
285 label='IIR h[n]')
286 axes[1, 1].set_xlabel('Sample n')
287 axes[1, 1].set_ylabel('Amplitude')
288 axes[1, 1].set_title('Impulse Response h[n]\n(FIR: finite; IIR: infinite)')
289 axes[1, 1].legend()
290 axes[1, 1].grid(True, alpha=0.3)
291 axes[1, 1].set_xlim(-1, 50)
292
293 # (2,0) Step response
294 axes[2, 0].plot(n, s_fir, 'b-o', markersize=3, linewidth=1.5, label='FIR s[n]')
295 axes[2, 0].plot(n, s_iir, 'r--^', markersize=3, linewidth=1.5, label='IIR s[n]')
296 axes[2, 0].axhline(1.0, color='k', linestyle=':', linewidth=1, label='Steady state = 1')
297 axes[2, 0].set_xlabel('Sample n')
298 axes[2, 0].set_ylabel('Amplitude')
299 axes[2, 0].set_title('Step Response s[n]')
300 axes[2, 0].legend()
301 axes[2, 0].grid(True, alpha=0.3)
302 axes[2, 0].set_xlim(-1, 50)
303
304 # (2,1) Noisy signal filtering
305 axes[2, 1].plot(t, noisy, 'gray', linewidth=0.8, alpha=0.6, label='Noisy signal')
306 axes[2, 1].plot(t, clean, 'k-', linewidth=2, alpha=0.5, label='Clean 5 Hz signal')
307 axes[2, 1].plot(t, y_fir, 'b-', linewidth=1.8, label=f'FIR filtered')
308 axes[2, 1].plot(t, y_iir, 'r--', linewidth=1.8, label=f'IIR filtered')
309 axes[2, 1].set_xlabel('Time (s)')
310 axes[2, 1].set_ylabel('Amplitude')
311 axes[2, 1].set_title('Filtering a Noisy Signal')
312 axes[2, 1].legend(fontsize=8)
313 axes[2, 1].grid(True, alpha=0.3)
314
315 plt.tight_layout()
316 out_path = '/opt/projects/01_Personal/03_Study/examples/Signal_Processing/08_filter_response.png'
317 plt.savefig(out_path, dpi=150)
318 print(f"\nSaved visualization: {out_path}")
319 plt.close()
320
321 print("\n" + "=" * 65)
322 print("Key Takeaways:")
323 print(" - FIR: finite impulse response, always stable, linear phase")
324 print(" - IIR: infinite impulse response, computationally efficient")
325 print(" - Linear phase (FIR) → constant group delay → no waveform distortion")
326 print(" - Group delay variation (IIR) distorts signal shape")
327 print(" - Both achieve unity DC gain but differ in roll-off shape")
328 print("=" * 65)