2. IIR梳状滤波器实现

任务目标

实现一个反馈结构的IIR梳状滤波器，推导其传递函数，并通过零极点图与频率响应分析其稳定性及滤波特性。滤波器差分方程形式为：

y[n] = x[n] + a * y[n - M]

其中M为延迟样本数，a为反馈系数。本案例使用Python实现，通过数值计算验证理论分析。

核心实现

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

def iir_comb_filter(x, M, a):
    """IIR梳状滤波器（反馈结构）
    y[n] = x[n] + a * y[n-M]
    """
    y = np.zeros_like(x, dtype=np.float64)
    for n in range(len(x)):
        if n - M >= 0:
            y[n] = x[n] + a * y[n - M]
        else:
            y[n] = x[n]  # 初始M个样本无反馈
    return y

# 参数设置
M = 8          # 延迟样本数
a = 0.8        # 反馈系数（|a|<1 保证稳定）
fs = 8000      # 采样率(Hz)

# 1. 频率响应分析（理论）
w, h = signal.freqz([1], [1, 0, 0, 0, 0, 0, 0, 0, -a], worN=2048, fs=fs)
# 注意：分母系数长度为M+1，索引0对应z^0，索引M对应z^{-M}

# 2. 零极点图
z = np.roots([1, 0, 0, 0, 0, 0, 0, 0, -a])  # 分母多项式求根
p = np.zeros(M)  # M重零点（分子为1，即z^M）

# 3. 验证：脉冲响应
impulse = np.zeros(64)
impulse[0] = 1.0
y_imp = iir_comb_filter(impulse, M, a)

# 4. 验证：正弦信号滤波
t = np.arange(256) / fs
x_test = np.sin(2 * np.pi * 500 * t) + 0.5 * np.sin(2 * np.pi * 2000 * t)
y_test = iir_comb_filter(x_test, M, a)

# 绘图
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 幅频响应
axes[0, 0].plot(w, 20 * np.log10(np.abs(h)))
axes[0, 0].set_title('幅度响应 (a=0.8, M=8)')
axes[0, 0].set_xlabel('频率 (Hz)')
axes[0, 0].set_ylabel('增益 (dB)')
axes[0, 0].grid(True)

# 零极点图
axes[0, 1].scatter(np.real(z), np.imag(z), marker='x', s=100, label='极点')
axes[0, 1].scatter(np.real(p), np.imag(p), marker='o', s=50, facecolors='none', edgecolors='b', label='零点(M重)')
axes[0, 1].axhline(0, color='k', lw=0.5)
axes[0, 1].axvline(0, color='k', lw=0.5)
axes[0, 1].add_patch(plt.Circle((0, 0), 1, fill=False, linestyle='--'))
axes[0, 1].set_title('零极点分布')
axes[0, 1].legend()
axes[0, 1].axis('equal')
axes[0, 1].grid(True)

# 脉冲响应
axes[1, 0].stem(range(len(y_imp)), y_imp)
axes[1, 0].set_title(f'脉冲响应 (a={a})')
axes[1, 0].set_xlabel('样本')
axes[1, 0].grid(True)

# 输入输出对比
axes[1, 1].plot(t, x_test, label='输入')
axes[1, 1].plot(t, y_test, label='输出', alpha=0.8)
axes[1, 1].set_title('滤波效果 (500Hz + 2000Hz 混合)')
axes[1, 1].set_xlabel('时间 (s)')
axes[1, 1].legend()
axes[1, 1].grid(True)

plt.tight_layout()
plt.show()

# 稳定性判定
print(f"极点模长最大值: {np.max(np.abs(z)):.4f}")
print(f"稳定性: {'稳定' if np.max(np.abs(z)) < 1 else '不稳定'}")

运行与验证

环境依赖：需要numpy、scipy、matplotlib。运行前请确认已安装：pip install numpy scipy matplotlib

运行方式：将上述代码保存为iir_comb.py，在终端执行python iir_comb.py。

验证方法：

1. 稳定性验证：观察控制台输出的极点模长最大值。当|a| < 1时，极点位于单位圆内，输出应显示“稳定”。可修改a = 1.2重新运行，此时应显示“不稳定”，且脉冲响应发散。
2. 频率响应验证：幅度响应图中，在频率f = k * fs / M（k为整数）处应出现峰值（如fs/M = 1000Hz的整数倍处），这是梳状滤波器的特征——在延迟频率的整数倍处产生谐振峰。
3. 滤波效果验证：输入信号包含500Hz和2000Hz两个分量。当M=8, fs=8000时，谐振频率为1000Hz的整数倍，因此2000Hz分量被增强，500Hz分量被抑制。观察输出波形中2000Hz成分的幅度变化。

稳定性与频率响应分析

传递函数为H(z) = 1 / (1 - a·z⁻ᴹ)，极点为z = a^(1/M) · e^(j·2πk/M)（k=0,1,...,M-1），即M个极点均匀分布在半径为|a|^(1/M)的圆上。当|a| < 1时所有极点位于单位圆内，系统稳定；当|a| ≥ 1时系统不稳定。反馈系数a越接近1，谐振峰越尖锐，频率选择性越强，但稳定性裕度越小。

关键结论：反馈系数a的值必须小于1才能保证系统稳定。极点均匀分布在复平面圆周上，其角度间隔为2π/M，对应频域中周期性出现的谐振峰。延迟M决定谐振峰之间的频率间隔（fs/M），a决定峰的尖锐程度。本实现中a=0.8, M=8时系统稳定，谐振峰位于1000Hz的整数倍处。