使用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?\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的可視化設計，還是有明顯的不足，需要不斷完善。

---

研究學習不易，點贊易。
工作生活不易，收藏易，點收藏不迷茫 ：）

---