用PythonNumPy实战解析LTI系统的零极点与频率响应数字信号处理的理论常常让初学者感到抽象难懂尤其是当教科书堆满数学公式时。但如果我们换一种方式——用代码和可视化来探索这些概念一切突然变得清晰起来。本文将带你用Python和NumPy库通过动手实践理解线性时不变(LTI)系统的核心特性。1. 从代码构建LTI系统理解LTI系统的最佳方式是从构建一个开始。在Python中我们可以用NumPy数组来表示系统的差分方程系数。假设我们有一个简单的二阶系统import numpy as np # 定义系统差分方程系数 b [0.1, 0.2, 0.1] # 分子系数 (x[n], x[n-1], x[n-2]) a [1, -1.2, 0.8] # 分母系数 (y[n], y[n-1], y[n-2])这个系统对应的差分方程为 y[n] 0.1x[n] 0.2x[n-1] 0.1x[n-2] 1.2y[n-1] - 0.8y[n-2]注意分母系数a[0]通常为1其他系数表示输出的延迟项权重我们可以用scipy.signal中的tf2zpk函数将这个传递函数转换为零极点形式from scipy import signal zeros, poles, gain signal.tf2zpk(b, a) print(零点位置:, zeros) print(极点位置:, poles) print(系统增益:, gain)2. 可视化零极点图零极点图是理解系统特性的强大工具。极点位置决定了系统的稳定性而零点影响频率响应特性。让我们绘制这个系统的零极点图import matplotlib.pyplot as plt def plot_pole_zero(zeros, poles): plt.figure(figsize(8,6)) # 绘制单位圆 theta np.linspace(0, 2*np.pi, 100) plt.plot(np.cos(theta), np.sin(theta), k--, alpha0.5) # 绘制零极点 plt.scatter(np.real(zeros), np.imag(zeros), markero, facecolorsnone, edgecolorsr, s100, label零点) plt.scatter(np.real(poles), np.imag(poles), markerx, colorb, s100, label极点) plt.axhline(0, colork, alpha0.2) plt.axvline(0, colork, alpha0.2) plt.grid(True) plt.legend() plt.xlabel(实部) plt.ylabel(虚部) plt.title(系统零极点图) plt.axis(equal) plt.show() plot_pole_zero(zeros, poles)从图中可以直观看出所有极点都在单位圆内 → 系统稳定零点位置影响特定频率的抑制极点靠近单位圆会导致共振峰3. 计算并绘制频率响应频率响应展示了系统对不同频率信号的增益和相位变化。我们可以计算并可视化它def plot_frequency_response(b, a, fs1.0): w, h signal.freqz(b, a) plt.figure(figsize(12, 6)) # 幅频响应 plt.subplot(2, 1, 1) plt.plot(w/np.pi, 20 * np.log10(abs(h)), b) plt.ylabel(幅度 (dB)) plt.title(频率响应) plt.grid(True) # 相频响应 plt.subplot(2, 1, 2) angles np.unwrap(np.angle(h)) plt.plot(w/np.pi, angles, g) plt.ylabel(相位 (弧度)) plt.xlabel(归一化频率 (×π rad/sample)) plt.grid(True) plt.tight_layout() plt.show() plot_frequency_response(b, a)关键观察点幅频响应的峰值对应极点位置幅频响应的谷值对应零点位置相位变化率与群延迟相关4. 零极点位置对系统的影响通过修改零极点位置我们可以直观看到系统特性的变化。让我们创建几个不同的系统并比较它们的响应4.1 极点位置的影响# 创建三个不同极点位置的系统 systems [ {name: 稳定系统, poles: [0.8*np.exp(1j*np.pi/4), 0.8*np.exp(-1j*np.pi/4)]}, {name: 临界稳定, poles: [0.99*np.exp(1j*np.pi/4), 0.99*np.exp(-1j*np.pi/4)]}, {name: 不稳定系统, poles: [1.1*np.exp(1j*np.pi/4), 1.1*np.exp(-1j*np.pi/4)]} ] plt.figure(figsize(15, 10)) for i, sys in enumerate(systems): # 从极点创建系统 b [1] a np.poly(sys[poles]) # 绘制冲激响应 t, imp signal.impulse((b, a), N50) plt.subplot(3, 2, 2*i1) plt.stem(t, imp) plt.title(f{sys[name]} - 冲激响应) # 绘制零极点图 plt.subplot(3, 2, 2*i2) plot_pole_zero([], sys[poles]) # 只绘制极点 plt.title(f{sys[name]} - 零极点图) plt.tight_layout() plt.show()从结果可以看出稳定系统的冲激响应逐渐衰减临界稳定系统的冲激响应持续振荡不稳定系统的冲激响应发散4.2 零点位置的影响# 创建三个不同零点位置的系统 systems [ {name: 低通特性, zeros: [-1]}, {name: 带阻特性, zeros: [np.exp(1j*np.pi/4), np.exp(-1j*np.pi/4)]}, {name: 高通特性, zeros: [1]} ] plt.figure(figsize(15, 8)) for i, sys in enumerate(systems): # 从零点创建系统固定极点 b np.poly(sys[zeros]) a [1, -0.8, 0.64] # 固定极点 # 绘制频率响应 w, h signal.freqz(b, a) plt.subplot(3, 2, 2*i1) plt.plot(w/np.pi, 20 * np.log10(abs(h))) plt.title(f{sys[name]} - 幅频响应) plt.grid(True) # 绘制零极点图 plt.subplot(3, 2, 2*i2) plot_pole_zero(sys[zeros], [0.8*np.exp(1j*np.pi/6), 0.8*np.exp(-1j*np.pi/6)]) plt.title(f{sys[name]} - 零极点图) plt.tight_layout() plt.show()观察结果零点在z1(ω0)时抑制低频 → 高通零点在z-1(ωπ)时抑制高频 → 低通零点在单位圆上特定角度抑制对应频率5. 实际应用设计简单滤波器理解了零极点与频率响应的关系后我们可以尝试设计简单的数字滤波器。例如设计一个带通滤波器# 设计带通滤波器 center_freq np.pi/3 # 中心频率 bandwidth np.pi/6 # 带宽 # 极点在中心频率附近单位圆内 r 0.9 # 极点半径 poles [r*np.exp(1j*center_freq), r*np.exp(-1j*center_freq)] # 零点在低频和高频处 zeros [1, -1] b np.poly(zeros) a np.poly(poles) # 绘制响应 plt.figure(figsize(12, 8)) plot_pole_zero(zeros, poles) plot_frequency_response(b, a) # 测试滤波器效果 t np.linspace(0, 10, 500) x np.sin(0.5*t) np.sin(center_freq*t) np.sin(2*t) # 混合频率信号 # 滤波 y signal.lfilter(b, a, x) plt.figure(figsize(12, 4)) plt.plot(t, x, b, alpha0.5, label输入信号) plt.plot(t, y, r, label滤波后信号) plt.legend() plt.title(带通滤波器效果) plt.grid(True) plt.show()这个简单的带通滤波器保留了中心频率附近的信号衰减了低频和高频成分极点半径r控制带宽(r越大带宽越窄)6. 高级话题全通系统与最小相位系统6.1 全通系统实现全通系统的幅频响应平坦只改变相位。我们可以用零极点对来构建# 创建全通系统 pole 0.8*np.exp(1j*np.pi/4) zero 1/np.conj(pole) # 零点与极点共轭倒数 b np.poly([zero]) a np.poly([pole]) # 验证频率响应 w, h signal.freqz(b, a) plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(w/np.pi, 20*np.log10(abs(h))) plt.title(幅频响应 (应为平坦)) plt.grid(True) plt.subplot(1, 2, 2) plt.plot(w/np.pi, np.unwrap(np.angle(h))) plt.title(相频响应) plt.grid(True) plt.show()6.2 最小相位系统特性最小相位系统具有最小的群延迟所有零点都在单位圆内# 比较最小相位和非最小相位系统 zeros_mp [0.50.5j, 0.5-0.5j] # 单位圆内 zeros_nmp [22j, 2-2j] # 单位圆外 b_mp np.poly(zeros_mp) b_nmp np.poly(zeros_nmp) a [1, -0.2, 0.9] # 固定极点 # 计算群延迟 w, gd_mp signal.group_delay((b_mp, a)) _, gd_nmp signal.group_delay((b_nmp, a)) plt.figure(figsize(12, 4)) plt.plot(w/np.pi, gd_mp, label最小相位系统) plt.plot(w/np.pi, gd_nmp, label非最小相位系统) plt.legend() plt.title(群延迟比较) plt.grid(True) plt.show()实际项目中我经常需要将非最小相位系统转换为最小相位形式这可以通过将单位圆外的零点反射到圆内来实现同时保持幅频响应不变。