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}\text{\hspace{1em}(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滤波器，对采样率$fs=2000Sa/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 = 2000
    numtaps = 185
    freqs = [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()

图2.窗函数设计法结果示意图

  (2) 最小均方误差法
  该方法的设计目标是使下面的表达式最小化：
$${\int }_0^{\pi }W(\omega )(A(\omega )-D(\omega ){)}^{2}d\omega$$
  在scipy库中，scipy.signal.firls函数实现了上述最优化问题。该函数用三个参数定义了期望滤波器的频响。第一个参数是bands，它是偶数长度的一组两两成对的频点，用于表征期望滤波器频响的重要频带信息。第二个参数是desired，它是跟bands一一对应的表征各个频点内增益的参数。第三个参数是系数向量$W(\omega )$，它是一个可选参数，长度是bands/desired的一半，用于表征各个频带的频响在最优化过程中所占的比重。下面通过一个例子说明该函数的用法：我们期望设计一个43阶的低通滤波器，系统采样率为200Hz。滤波器的通带为[0,15]Hz，阻带是[30,100]Hz，过渡带是[15,30]Hz，过渡带内的增益由1线性递减到0。代码如下。下面代码中bands和desired均为偶数，0~15Hz为滤波器通带，增益为1,；15-30Hz为过渡带，其中15Hz上增益为1，30Hz上增益为0，表示过渡带增益由1递减到0;30-100Hz为阻带，增益为0。weights的长度为bands/desired的一半，weight_1中，3个元素均为1，表示在最优化过程中通带、过渡带和阻带保持相同的权重。weight_2中，第一个元素为100，第二个元素为0.01，第三个元素为1，表示在最优化过程中我们最关注通带的性能，对过渡带的性能要求最低。从图3中也能看出，weight_2对应的结果在通带内的表现优于weight_1，在过渡带内的性能却不如weight_1,与预期是相符的。

from scipy.signal import freqz, firls
import matplotlib.pyplot as plt
if __name__ == "__main__":
    fs = 200
    bands = np.array([0, 15, 15, 30, 30, 100])
    desired = np.array([1, 1, 1, 0, 0, 0])
    numtaps = 43
    weight_1 = np.array([1, 1, 1])
    weight_2 = np.array([100, 0.01, 1])
    taps_1 = firls(numtaps, bands, desired, nyq=fs / 2, weight=weight_1)
    taps_2 = firls(numtaps, bands, desired, nyq=fs / 2, weight=weight_2)
    w_1, h_1 = freqz(taps_1, worN=2000)
    w_2, h_2 = freqz(taps_2, worN=2000)
    plt.figure()
    plt.plot(bands, desired, 'r--')
    plt.plot(w_1 * fs / (2 * np.pi), np.abs(h_1), 'b')
    plt.plot(w_2 * fs / (2 * np.pi), np.abs(h_2), 'k')
    plt.legend(['ideal', 'weights=[1, 1, 1]', 'weights=[100, 0.01, 1]'])
    plt.xlabel('frequency(Hz)')
    plt.ylabel('gain')
    plt.show()

图3.firls函数设计结果

  需要说明的是，基于最小二乘的方法，当weights中所有元素保持一致时，它的效果和boxcar窗函数的设计效果是一致的。基于最小二乘方法的一大好处就是可以根据要求对不同频段设置不同权重。

  (3) Parks-McClellan法
  Parks-McClellan法是基于remze算法实现的。它本质上是一个最小化最大值的优化问题，如下式所示
$$E(\omega )=W(\omega )(A(\omega )-D(\omega ))$$
  Parks-McClellan法的目的是使上式中的$E(\omega )$的最大值最小化，scipy中是通过函数remez来实现的，下面通过一个例子进行说明：我们期望设计一个带通滤波器，系统采样率是2000Hz，阻带是[0,250]Hz和[700,1000]Hz，通带是[350,550]Hz。它和firls函数类似，都需要通过bands参数指定各个频带，且bands均为偶数。不过它的desire参数与firls不一样，它的长度是bands的一半。设计代码如下：

from scipy.signal import freqz, remez
import matplotlib.pyplot as plt
if __name__ == "__main__":
    fs = 2000
    bands = [0, 250, 350, 550, 700, 0.5 * fs]
    desired = [0, 1, 0]
    numtaps = 31
    weight_1 = [1, 1, 1]
    weight_2 = [1, 25, 1]

    taps_1 = remez(numtaps, bands, desired, weight_1, fs=fs)
    taps_2 = remez(numtaps, bands, desired, weight_2, fs=fs)
    w_1, h_1 = freqz(taps_1, worN=2000)
    w_2, h_2 = freqz(taps_2, worN=2000)
    plt.figure()
    plt.plot(bands, [0, 0, 1, 1, 0, 0], 'r--')
    plt.plot(w_1 * fs / (2 * np.pi), np.abs(h_1), 'b')
    plt.plot(w_2 * fs / (2 * np.pi), np.abs(h_2), 'k')
    plt.legend(['ideal', 'weights=[1, 1, 1]', 'weights=[1, 25, 1]'])
    plt.xlabel('frequency(Hz)')
    plt.ylabel('gain')
    plt.show()

图4.remez函数设计结果   利用remez进行滤波器设计时，最好看一下设计出来的滤波器频响，因为它可能出现不收敛的情况，此时可以考虑通过参数maxiter参数增大迭代次数，或通过grid_density增大插值密度。

未完待续…
未完待续…
未完待续…
未完待续…