python数字信号处理之观测任意数字滤波器的频率响应：从原理解析到代码实现

数字滤波器分为FIR、IIR两种，其均可以表示为抽头系数a，b的形式，数字滤波器的传递函数一般表示为：

$$H(z)=\frac{B(z)}{A(z)}=\frac{b_0+b_1{z}^{-1}+\cdots +b_M{z}^{-M}}{a_0+a_1{z}^{-1}+\cdots +a_N{z}^{-N}}$$

其中，$B(z)$ 和 $A(z)$ 分别是分子和分母多项式，系数由 b 和 a 参数提供。

在拿到a,b之后，可以使用signal.freqz就可以观察任意数字滤波器的频率响应。下面是原理分析。

1. 底层原理：滤波器的频率响应等于传递函数的分子和分母的频域响应之比

freqz 通过计算数字滤波器的传递函数在单位圆（即$z=e^{-jw}$）上的 Z 变换实现频率响应分析。对于分子系数为 b、分母系数为 a 的滤波器，其传递函数为：

$$H(e^{jw})=\frac{B(e^{jw})}{A(e^{jw})}=\frac{b_0+b_1{e}^{-jw}+\cdots +b_M{e}^{-jwM}}{a_0+a_1{e}^{-jw}+\cdots +a_N{e}^{-jwN}}=\frac{{\sum }_{k=0}^M{b}_k{e}^{-j\omega k}}{{\sum }_{k=0}^N{a}_k{e}^{-j\omega k}}=\frac{DTFT(b_k)}{DTFT(a_k)}$$

其中 $w=2\pi f/fs$，$f$ 为实际频率（单位：Hz），$fs$ 为采样频率。此时 Z 变换退化为离散时间傅里叶变换（DTFT）。 可见，频率响应即分子和分母多项式在单位圆上的复频域响应之比。

**2. **freqz函数的实现原理

1. 单位圆上的均匀采样：

• 在单位圆上选取 worN 个均匀分布的点，对应频率范围 $0\le \omega \le 2\pi$（若 whole=True）或 $0\le \omega \le \pi$（默认）。
• 当设置 fs（采样频率）时，频率单位从弧度转换为实际 Hz：$\omega =2\pi f/fs$。

2. 分子/分母多项式系数的FFT计算：

• 对分子系数 b 和分母系数 a 分别进行零填充至长度 worN，然后计算其 FFT，得到 $B(e^{j\omega })$ 和 $A(e^{j\omega })$ 的频域响应。

3. 频率响应为分子分母FFT结果相除：

$$H(e^{j\omega })=\frac{B(e^{j\omega })}{A(e^{j\omega })}$$

4. 归一化与频率映射：

• 若未指定 fs，输出频率单位为归一化弧度（$0\sim \pi$ 对应 $0\sim f_s/2$）；若指定 fs，频率直接映射为实际 Hz 值。

5. 输出结果

返回复数频率响应 $h$ 和对应频率数组 $w$，可通过 np.abs(h) 和 np.angle(h) 分别获取幅频响应（dB或线性单位）和相频响应（弧度或角度）。

---

参数解析

参数	说明
`b`	分子系数数组（对应传递函数分子多项式系数）。
`a`	分母系数数组（默认值为1，即FIR滤波器）。
`worN`	频率点数：若为整数，默认在有效频段内均匀采样；若为数组，直接指定计算频率。
`whole`	若为 `True`，计算全频段（0到fs）；否则仅计算0到fs/2。
`fs`	采样频率（单位：Hz），用于将频率单位从弧度转换为实际Hz（默认值为 $2\pi$，即弧度单位）。
`include_nyquist`	当 `whole=False` 且 `worN` 为整数时，是否包含奈奎斯特频率点（默认不包含）。

---

典型用法示例

1. 基本调用

   import scipy.signal as signal
   b = [0.8, 0.5, 0.6]  # 分子系数（FIR滤波器）
   a = [1, 0.2, 0.4]    # 分母系数
   w, h = signal.freqz(b, a, worN=512, fs=2000)
   

2. 绘制幅频和相频响应

   import matplotlib.pyplot as plt
   plt.figure()
   plt.subplot(2, 1, 1)
   plt.plot(w, 20 * np.log10(np.abs(h)))  # 幅频响应（dB单位）
   plt.ylabel('Amplitude (dB)')
   plt.subplot(2, 1, 2)
   plt.plot(w, np.angle(h) * 180 / np.pi)  # 相频响应（角度单位）
   plt.xlabel('Frequency (Hz)')
   plt.ylabel('Phase (degrees)')
   

3. 自定义频率点

   custom_freqs = np.array([50, 100, 150])  # 指定计算50Hz、100Hz、150Hz处的响应
   w, h = signal.freqz(b, a, worN=custom_freqs, fs=2000)
   
   

---

注意事项

• 系数维度兼容性：b 和 a 的维度需满足广播规则，多维数组需确保形状兼容。

• 频率单位转换：若需以Hz为单位显示频率，必须正确设置 fs 参数，避免混淆弧度与Hz的关系。

• 性能优化：当 worN 为2的幂次时，FFT计算速度更快。