[python] 使用python设计滤波器

使用python设计滤波器

习惯了python后，matlab逐渐成为了牛夫人，尤其是漂亮国对龙国制裁后，我作为有骨气其的码农，虽然有和谐版的matlab可用，但是终究是放不下码农的尊严，不到万不得一，已经很少用matlab了。绝大部分仿真工作，都已经移步到python，但最近要进行滤波器设计，为了不使用matlab的fdatool，便尝试着用python设计滤波器。

闲话少说，代码说话：

---

完整滤波器设计代码（未经完整验证，博主还在不断完善中）

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from matplotlib.ticker import EngFormatter# 设计参数
fs = 10000       # 采样率 10kHz
nyq = fs / 2     # 奈奎斯特频率
passband = 1500  # 通带截止频率 (Hz)
stopband = 2000  # 阻带截止频率 (Hz)
pass_ripple = 1  # 通带波动 (dB)
stop_atten = 40  # 阻带衰减 (dB)# 计算归一化频率
wp = passband / nyq
ws = stopband / nyqdef design_iir_filter():"""设计椭圆IIR滤波器"""# 计算最小阶数和自然频率order, wn = signal.ellipord(wp, ws, pass_ripple, stop_atten)# 设计椭圆滤波器b, a = signal.ellip(order, pass_ripple, stop_atten, wn, btype='low', analog=False)return b, a, orderdef design_fir_filter():"""设计FIR滤波器（凯塞窗方法）"""# 计算过渡带宽度width = abs(stopband - passband) / nyq# 估算凯塞窗参数A = stop_attenbeta = 0.1102*(A - 8.7) if A > 50 else 0.5842*(A - 21)**0.4 + 0.07886*(A - 21)# 计算所需阶数numtaps = int((A - 8) / (2.285 * 2 * np.pi * width)) + 1# 设计FIR滤波器taps = signal.firwin(numtaps, wn=passband/nyq, window=('kaiser', beta), pass_zero='lowpass')return taps, numtapsdef analyze_filter(b, a=None):"""分析滤波器性能"""fig, (ax_mag, ax_phase, ax_grp, ax_zp) = plt.subplots(4, 1, figsize=(10, 12))# 幅频响应w, h = signal.freqz(b, a, worN=8000, fs=fs)ax_mag.plot(w, 20 * np.log10(np.abs(h)), color='blue')ax_mag.set_title('幅频响应')ax_mag.set_ylabel('幅度 (dB)')ax_mag.grid(True, which='both', linestyle='--')ax_mag.axvline(passband, color='g', linestyle='--', alpha=0.7)ax_mag.axvline(stopband, color='r', linestyle='--', alpha=0.7)ax_mag.set_ylim([-stop_atten-20, 5])# 相频响应angles = np.unwrap(np.angle(h))ax_phase.plot(w, angles, color='green')ax_phase.set_title('相频响应')ax_phase.set_ylabel('相位 (弧度)')ax_phase.grid(True)# 群延迟grp_delay = -np.diff(angles) / np.diff(w)ax_grp.plot(w[:-1], grp_delay, color='red')ax_grp.set_title('群延迟')ax_grp.set_ylabel('采样点数')ax_grp.grid(True)# 零极点图if a is not None:  # IIR滤波器z, p, k = signal.tf2zpk(b, a)ax_zp.scatter(np.real(z), np.imag(z), marker='o', facecolors='none', edgecolors='b', s=80)else:  # FIR滤波器z = np.roots(b)p = np.zeros(len(z))ax_zp.scatter(np.real(z), np.imag(z), marker='o', facecolors='none', edgecolors='b', s=80)ax_zp.scatter(np.real(p), np.imag(p), marker='x', color='r', s=80)unit_circle = plt.Circle((0,0), radius=1, fill=False, color='gray', linestyle='--')ax_zp.add_patch(unit_circle)ax_zp.set_title('零极点图')ax_zp.set_xlabel('实部')ax_zp.set_ylabel('虚部')ax_zp.grid(True)ax_zp.axis('equal')ax_zp.axhline(0, color='black', linewidth=0.5)ax_zp.axvline(0, color='black', linewidth=0.5)plt.tight_layout()return fig# 主程序
if __name__ == "__main__":# 设计IIR滤波器b_iir, a_iir, iir_order = design_iir_filter()print(f"IIR滤波器阶数: {iir_order}")# 设计FIR滤波器fir_taps, fir_order = design_fir_filter()print(f"FIR滤波器阶数: {fir_order}")# 分析IIR滤波器fig_iir = analyze_filter(b_iir, a_iir)fig_iir.suptitle('椭圆IIR滤波器分析', fontsize=16)# 分析FIR滤波器fig_fir = analyze_filter(fir_taps)fig_fir.suptitle('FIR滤波器分析', fontsize=16)# 应用滤波器示例t = np.linspace(0, 1, fs)  # 1秒时长sig = np.sin(2*np.pi*1000*t) + 0.5*np.sin(2*np.pi*3000*t)# IIR滤波filtered_iir = signal.lfilter(b_iir, a_iir, sig)# FIR滤波filtered_fir = signal.lfilter(fir_taps, 1.0, sig)# 绘制结果plt.figure(figsize=(10, 6))plt.plot(t[:500], sig[:500], 'b-', label='原始信号')plt.plot(t[:500], filtered_iir[:500], 'r-', label='IIR滤波后')plt.plot(t[:500], filtered_fir[:500], 'g-', label='FIR滤波后')plt.title('时域滤波效果对比 (前500点)')plt.xlabel('时间 (秒)')plt.ylabel('幅度')plt.legend()plt.grid(True)plt.tight_layout()plt.show()

---

关键原理与代码对应说明

1. 滤波器类型选择

• IIR滤波器：椭圆滤波器（Elliptic）提供最陡峭过渡带

  order, wn = signal.ellipord(wp, ws, pass_ripple, stop_atten)
  b, a = signal.ellip(order, pass_ripple, stop_atten, wn, btype='low')
  

  • ellipord计算最小阶数和自然频率
  • 椭圆滤波器在通带/阻带均有等波纹特性

• FIR滤波器：凯塞窗（Kaiser）提供参数化控制

  beta = 0.1102*(A - 8.7)  # 窗形参数计算
  taps = signal.firwin(numtaps, cutoff, window=('kaiser', beta))
  

  • 凯塞窗通过beta控制主瓣宽度和旁瓣衰减平衡

2. 阶数估算原理

• IIR阶数估算：
  $$N=\frac{K({\omega }_s)K(\sqrt{1-{\omega }_p^2})}{K({\omega }_p)K(\sqrt{1-{\omega }_s^2})}$$
  其中K为第一类完全椭圆积分

• FIR阶数估算：
  $$N\approx \frac{A-8}{2.285\cdot \Delta \omega }$$
  $\Delta \omega$为过渡带宽度，A为阻带衰减(dB)

3. 性能分析技术

• 幅频响应：signal.freqz计算复数频率响应

  w, h = signal.freqz(b, a, worN=8000, fs=fs)
  20*np.log10(np.abs(h))  # 转换为dB
  

• 相位特性：

  angles = np.unwrap(np.angle(h))  # 解卷绕相位
  

• 群延迟：相位导数计算

  -np.diff(angles)/np.diff(w)  # 群延迟 = -dφ/dω
  

• 稳定性分析：

  • IIR：极点是否在单位圆内
  • FIR：恒稳定（全零点系统）

4. 设计参数调整指南

1. 过渡带陡度：

   • IIR：增加阶数
   • FIR：增加窗长度

2. 阻带衰减：

   • IIR：增加stop_atten参数
   • FIR：增大凯塞窗的beta值

3. 计算效率：

   • IIR：阶数低但非线性相位
   • FIR：高阶但线性相位

---

此方案应该能替代MATLAB FDATool的核心功能，提供从设计到分析的完整工作流。可根据具体应用调整滤波器类型（低通/高通/带通）和设计参数。但是比起FDATool的可视化设计，还是有明显的不足，需要不断完善。

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---