📚 课程目录第1章基础梳状滤波器设计

1. 基础梳状滤波器设计

任务目标

设计一个阶数为 M 的FIR梳状滤波器，其系统函数为 H(z) = 1 - z^-M，该滤波器会在 2πk/M（k = 0, 1, …, M-1）频率处产生周期性陷波。本案例要求实现该滤波器并绘制其幅度响应，验证陷波位置与理论一致。

实现原理

滤波器 H(z) = 1 - z^-M 的冲激响应为 h[0]=1，h[M]=-1，其余为0。其频率响应为 H(e^jω) = 1 - e^-jωM，幅度响应为 |H(e^jω)| = 2|sin(ωM/2)|。当 ωM/2 = kπ 即 ω = 2πk/M 时幅度为零，形成等间隔分布的陷波。

完整代码

import numpy as np
import matplotlib.pyplot as plt

# 滤波器参数
M = 8          # 延迟阶数（陷波数量 = M）
N = 2048       # 频率采样点数

# 冲激响应：h[0]=1, h[M]=-1
h = np.zeros(M + 1)
h[0] = 1.0
h[M] = -1.0

# 计算频率响应
w = np.linspace(0, np.pi, N)          # 0 ~ π 弧度/采样
H = np.exp(-1j * w * 0) - np.exp(-1j * w * M)   # 等效于 1 - e^{-jωM}
mag = np.abs(H)                        # 幅度响应

# 理论陷波位置（0 ~ π 范围内）
notch_freqs = [2 * np.pi * k / M for k in range(M // 2 + 1)
               if 2 * np.pi * k / M <= np.pi]

# 绘制幅度响应
plt.figure(figsize=(10, 4))
plt.plot(w, mag, 'b-', linewidth=1.5)
plt.xlabel('频率 (rad/sample)')
plt.ylabel('幅度 |H(e^{jω})|')
plt.title(f'FIR梳状滤波器幅度响应 (M={M}, H(z)=1-z⁻ᴹ)')
plt.grid(True, alpha=0.3)
plt.xlim(0, np.pi)
plt.ylim(0, 2.2)

# 标记陷波位置
for f in notch_freqs:
    plt.axvline(f, color='r', linestyle='--', alpha=0.7, linewidth=1)
    plt.text(f, 0.1, f'{f:.2f}', fontsize=8, ha='center', color='red')

plt.tight_layout()
plt.show()

# 验证：打印陷波频率处的幅度值（应接近0）
print("陷波频率处的幅度值验证：")
for f in notch_freqs:
    idx = np.argmin(np.abs(w - f))
    print(f"  ω = {f:.4f} rad → |H| = {mag[idx]:.6f}")

运行方式

安装依赖：pip install numpy matplotlib
将代码保存为 comb_filter.py 并运行：python comb_filter.py

验证方法

运行后应弹出幅度响应图，图中蓝色曲线在红色虚线标记的频率处（ω = 0, π/4, π/2, 3π/4, π）幅度降为0。
控制台输出的幅度值应接近0（数值误差小于1e-5），证明陷波位置准确。
相邻两个陷波之间的间隔为 π/M = π/8 rad/sample，与理论一致。

预期输出

陷波频率处的幅度值验证：
  ω = 0.0000 rad → |H| = 0.000000
  ω = 0.7854 rad → |H| = 0.000000
  ω = 1.5708 rad → |H| = 0.000000
  ω = 2.3562 rad → |H| = 0.000000
  ω = 3.1416 rad → |H| = 0.000000

图中红色虚线标记了5个陷波位置（含边界），相邻陷波间隔为 π/8，曲线峰值均为2，与理论公式 |H(e^jω)| = 2|sin(4ω)| 完全吻合。该滤波器可用于消除特定基频及其谐波干扰，例如去除50Hz工频及其整数倍噪声。