[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 / nyq

def 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, order

def design_fir_filter():
    """设计FIR滤波器（凯塞窗方法）"""
    # 计算过渡带宽度
    width = abs(stopband - passband) / nyq
    # 估算凯塞窗参数
    A = stop_atten
    beta = 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, numtaps

def 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