Scipy库中FIR滤波器的应用

在上一篇文章《Scipy库中IIR滤波器的应用》中，我们阐述了利用Scipy库进行IIR滤波器设计的一些基本做法。在这篇文章中我们将进一步总结Scipy库在FIR滤波器设计中1的应用。

1. FIR滤波器基本概念
在上篇文章中，我们在给出线性滤波器的差分方程喝系统函数的一般形式时指出FIR滤波器是一个无反馈的全零点型滤波器。设输入序列为 $x_n$ ，系数为 $b_m$ ，则滤波输出为
$y_n=\sum_{m=0}^M b_m x_{n-m} \tag{1}$
上面的表达式本质上是抽头系数序列和输入序列的卷积。其中， $M$ 称为滤波器的阶数，一般 $M$ 阶FIR滤波器有 $M + 1$ 个抽头系数。
在Scipy中可以调用scipy.signal.lfilter或scipy.signal.convolve函数来对输入序列进行FIR滤波。需要主要的是对于一个有限长输入序列 $x_0,x_1,...,x_N$ ，式(1)并没有指明当 $n < M$ 时如何处理。函数 scipy.signal.convolve中有一个入参mode可以用于指定上述问题的处理方法。比如当mode='valid’时，将严格按照(1)进行计算，当mode='same’时会对输入序列进行补零，使得输入和输出长度保持一致。
FIR滤波器的频响可以调用scipy.signal.freqz得到，它的输入为抽头系数及频率点数，下面以滑动平均滤波器为例进行说明。函数freqz返回的是归一化角频率，范围是 $0-\pi$ ，如果已知信号采样率 $f_s$ ，可将其乘以 $f_s/(2\pi)$ 将横坐标转化为以Hz为单位的频率值。

import matplotlib.pyplot as plt
from scipy.signal import freqz
if __name__ == "__main__":for n in [3, 7, 21]:taps = np.full(n, fill_value=1.0 / n)w, h = freqz(taps, worN=2000)plt.plot(w, np.abs(h), label='n=%d' % n)plt.show()

在这里插入图片描述

图1.滑动平均滤波器的频响

2. FIR滤波器设计
Scipy数据库中提供了四种FIR滤波器设计方法，分别是：
a) 窗函数法：该方法首先计算理想FIR滤波器的冲激响应，然后对其进行加窗得到最终的滤波器冲激响应；
b) 最小均方误差法：该方法使得设计的FIR滤波器频响与理想滤波器频响的均衡误差最小。
c) Parks-McClellan法：该方法是一种等纹波设计方法，它使得设计的FIR滤波器和理想FIR滤波器频响的最大误差最小化。
d) 线性规划法：该方法其实是Parks-McClellan法的一种线性规划解法。
在详细描述上述方法之前，首先定义 $\omega$ 为角频率， $A(\omega)$ 为所设计的滤波器频响的实部。 $D(\omega)$ 为理想滤波器的频响， $W(\omega)$ 为权重向量，用于调整误差 $A(\omega)-D(\omega)$ 的权重。

(1) 窗函数法
窗函数法的基本步骤是首先计算理想FIR滤波器的冲激响应，然后利用窗函数对其进行截断得到最终FIR滤波器的冲激响应。scipy.signal中有两个窗函数设计的函数firwin和firwin2。本小节中只展示firwin2的用法，firwin的用法将在后面介绍Kaiser窗函数设计法的时候讲到。
在这边我们打算设计一个185阶的FIR滤波器，对采样率 $f s = 2000 S a / s$ 的信号进行滤波。这个滤波器整体表现为低通特性，它的过渡带是150Hz~175Hz。我们还希望这个滤波器在48Hz-72Hz之间有一个倒三角的凹坑，凹坑中心在60Hz处，且凹坑顶点处的增益为0.1。实现代码如下。firwin2函数的输入参数包括滤波器阶数、频点(各个频率转折点，范围是0-fs)、各个频点对应的增益、奈奎斯特频率fs/2，窗函数类型。下面的案例中对比了矩形窗、hamming窗和凯撒窗的结果。

import matplotlib.pyplot as plt
from scipy.signal import freqz, firwin2
if __name__ == "__main__":fs = 2000numtaps = 185freqs = [0, 48, 60, 72, 150, 175, 1000]gains = [1, 1, 0.1, 1, 1, 0, 0]taps_rect = firwin2(numtaps, freqs, gains, nyq=0.5 * fs, window=None)taps_hamming = firwin2(numtaps, freqs, gains, nyq=0.5 * fs, window='hamming')taps_kaiser = firwin2(numtaps,freqs, gains, nyq=0.5 * fs, window=('kaiser', 2.7))w_rect, h_rect = freqz(taps_rect, worN=2000)w_hamming, h_hamming = freqz(taps_hamming, worN=2000)w_kaiser, h_kaiser = freqz(taps_kaiser, worN=2000)plt.figure()plt.plot(freqs, gains, 'r--')plt.plot(w_rect * fs / (2 * np.pi), np.abs(h_rect), 'b')plt.plot(w_hamming * fs / (2 * np.pi), np.abs(h_hamming), 'y')plt.plot(w_kaiser * fs / (2 * np.pi), np.abs(h_kaiser), 'k')plt.legend(['ideal', 'rect', 'hamming', 'kaiser'])plt.xlabel('frequency(Hz)')plt.ylabel('gain')plt.show()