别再死记硬背了！用Python+零极点法，5分钟搞定一个IIR低通滤波器

用Python玩转零极点设计：5分钟实现IIR低通滤波器的实战指南

数字信号处理课程里那些抽象的零极点图是否让你头疼不已？当课本上的理论公式遇上Python代码，一切都变得直观起来。本文将带你用SciPy和NumPy，通过调整极点在Z平面上的位置，亲手打造一个可听、可看、可调的IIR低通滤波器。不同于传统教材的数学推导，我们将用交互式可视化让你真正理解极点位置如何影响滤波器的带宽和陡峭度。

1. 零极点法的核心原理可视化

零极点设计法的精髓在于：极点是频率的放大器，零点是频率的消除器。在Z平面上，极点靠近单位圆的位置决定了滤波器要增强的频率范围。对于低通滤波器，我们需要在低频区域（靠近z=1）放置极点。

让我们用代码生成一个极点在0.9的一阶低通滤波器：

import numpy as np import matplotlib.pyplot as plt from scipy import signal # 创建极点在0.9的一阶系统 b = [1] # 分子系数（无零点） a = [1, -0.9] # 分母系数（极点在0.9） sys = signal.TransferFunction(b, a, dt=1) # 绘制零极点图 plt.figure(figsize=(10,4)) plt.subplot(121) plt.scatter(np.real(sys.poles), np.imag(sys.poles), marker='x', color='r') plt.scatter(np.real(sys.zeros), np.imag(sys.zeros), marker='o', color='b') unit_circle = plt.Circle((0,0), 1, fill=False, linestyle='--') plt.gca().add_patch(unit_circle) plt.axis('equal') plt.grid() # 绘制频率响应 plt.subplot(122) w, h = signal.freqz(b, a) plt.plot(w, 20*np.log10(np.abs(h))) plt.title('频率响应') plt.ylabel('幅度 (dB)') plt.xlabel('频率 [rad/sample]') plt.tight_layout() plt.show()

运行这段代码，你会看到两个关键信息：

1. 左图显示极点在0.9+0j的位置（红色×），没有零点（蓝色○）
2. 右图显示该滤波器确实让低频通过，高频衰减

极点位置与带宽的关系可以用以下经验公式估算：

带宽 ≈ 1 - |极点半径|

所以当极点在0.9时，带宽约为0.1 rad/sample。将极点改为0.95再运行，会发现通带变窄。

2. 二阶Butterworth滤波器的实战设计

一阶滤波器过渡带太缓？二阶滤波器能提供更陡峭的滚降。Butterworth滤波器以其最大平坦的通带响应著称，让我们用零极点法实现它：

def design_butterworth_lp(cutoff, dt=1): """设计二阶Butterworth低通滤波器""" theta = 2 * np.pi * cutoff * dt pole_radius = 1 - theta/2 # 经验公式 # 共轭极点对 pole1 = pole_radius * np.exp(1j*theta/2) pole2 = np.conj(pole1) # 系统函数 b = [1, 0, 0] # 无有限零点 a = np.poly([pole1, pole2]) return signal.TransferFunction(b, a, dt=dt) # 设计截止频率0.2的二阶Butterworth cutoff = 0.2 sys = design_butterworth_lp(cutoff) # 绘制响应曲线 w, h = signal.freqz(sys.num, sys.den) plt.plot(w/np.pi, 20*np.log10(np.abs(h)), label=f'截止={cutoff}') plt.axvline(cutoff, color='r', linestyle='--') plt.legend(); plt.grid()

关键设计参数对照表：

提示：Butterworth滤波器的极点均匀分布在s平面的左半圆上，经过双线性变换映射到z平面

3. 滤波器性能的量化评估

设计完滤波器后，我们需要量化评估三个关键指标：

1. 通带波纹：通带内最大波动（理想应<1dB）
2. 阻带衰减：阻带最小衰减（理想>40dB）
3. 过渡带宽：通带到阻带的转变速度

改进前面的Butterworth设计，添加性能评估代码：

def evaluate_filter(b, a, cutoff): w, h = signal.freqz(b, a) freq = w / np.pi # 通带指标(0到0.8*cutoff) passband_mask = freq <= 0.8*cutoff passband_ripple = np.max(np.abs(h[passband_mask])) - np.min(np.abs(h[passband_mask])) # 阻带指标(1.2*cutoff到Nyquist) stopband_mask = freq >= 1.2*cutoff stopband_atten = -20*np.log10(np.max(np.abs(h[stopband_mask]))) # 过渡带(0.8*cutoff到1.2*cutoff) transition = 0.4*cutoff print(f"通带波纹: {20*np.log10(1+passband_ripple):.2f} dB") print(f"阻带衰减: {stopband_atten:.2f} dB") print(f"过渡带宽: {transition:.3f}π rad/sample") # 评估Butterworth设计 evaluate_filter(sys.num, sys.den, cutoff)

典型输出结果：

通带波纹: 0.02 dB 阻带衰减: 24.15 dB 过渡带宽: 0.080π rad/sample

4. 高阶滤波器的级联实现

当需要更严格的滤波特性时，可以采用多个二阶节（Biquad）级联的方式实现高阶滤波器。这种方法数值稳定性更好：

def design_iir_lp(order, cutoff): """设计高阶IIR低通滤波器""" # 确定极点位置 poles = [] for k in range(1, order+1, 2): angle = np.pi/2 * (2*k + order - 1)/(2*order) radius = 1 - cutoff/2 # 经验公式 pole = radius * np.exp(1j*angle) poles.extend([pole, np.conj(pole)]) # 转换为二阶节 sos = [] for i in range(0, len(poles), 2): b = [1, 0, 0] # 无有限零点 a = np.poly(poles[i:i+2]) sos.append([b, a]) return sos # 设计8阶IIR滤波器 sos = design_iir_lp(8, 0.2) # 绘制每个二阶节的响应 plt.figure(figsize=(10,6)) for i, (b, a) in enumerate(sos): w, h = signal.freqz(b, a) plt.plot(w/np.pi, 20*np.log10(np.abs(h)), label=f'Biquad {i+1}') plt.title('各二阶节频率响应'); plt.grid(); plt.legend()

高阶滤波器设计要点：

• 每增加一阶，过渡带斜率增加约20dB/decade
• 极点应均匀分布在Z平面低频区域
• 实际实现时应采用二阶节级联，避免直接高阶多项式运算

5. 实时音频滤波演示

最后让我们用PyAudio实现一个实时低通滤波器，直观感受极点位置对声音的影响：

import pyaudio def audio_filter_effect(pole_radius=0.9): p = pyaudio.PyAudio() stream = p.open(format=pyaudio.paFloat32, channels=1, rate=44100, input=True, output=True, frames_per_buffer=1024, stream_callback=make_callback(pole_radius)) print("实时滤波运行中... 按Ctrl+C停止") try: while True: pass except KeyboardInterrupt: stream.stop_stream() stream.close() p.terminate() def make_callback(pole_radius): # 初始化滤波器状态 x_prev = 0 y_prev = 0 def callback(in_data, frame_count, time_info, status): nonlocal x_prev, y_prev x = np.frombuffer(in_data, dtype=np.float32) y = np.zeros_like(x) # 应用一阶IIR滤波: y[n] = x[n] + pole_radius*y[n-1] for n in range(len(x)): y[n] = x[n] + pole_radius * y_prev y_prev = y[n] return (y.tobytes(), pyaudio.paContinue) return callback # 启动实时滤波（极点在0.95） audio_filter_effect(pole_radius=0.95)

操作建议：

1. 对着麦克风说话或播放音乐
2. 尝试不同pole_radius值（0.8-0.99）
3. 观察极点越接近1时，低频保留越完整，高频衰减越明显

这个实时演示生动展示了：极点的位置直接决定了哪些频率成分能被保留。当你想快速验证一个滤波器设计时，这种听觉反馈比任何图表都更直观。