FIR-IIR滤波器设计1：窗函数法设计FIR数字滤波器

0、基础知识

FIR滤波器的z变换为：
$$H(z)=\sum _{n=0}^{N-1}h(n)z^{-n}$$
它在z平面上有N−1个零点并在原点z=0处有N−1重极点，所以一般来说H(z)永远为稳定系统。

FIR滤波器设计的基本思想就是寻求一个有限长序列作为FIR系统的单位脉冲响应h(n)，使其频率响应函数$H\left(e^{jw}\right)$逼近所期望的频率响应函数$H_{id}\left(e^{jw}\right)$。这种逼近有两种不同的途径，即从时域角度进行逼近或从频域角度逼近。

• 时域角度: 以$H_{id}\left(e^{jw}\right)$所对应的$h_{id}(n)$为目标，构造一个有限长的h(n)作为所设计滤波器的单位脉冲响应，这种设计方法称为窗函数法。

• 频域角度: 使所设计的滤波器频率响应在某些样本点处与期望特性相同，以此达到逼近的目的，这种设计方法就称为频率采样法。

  另外还可以按照某种最佳逼近准则，使所设计的滤波器在此准则下最优，以此来达到满足性能指标的最优滤波器设计。最优等波纹线性相位滤波器设计方法就是这种方法。。

首先我们介绍窗函数法：

窗函数的基本思想就是拿一个窗函数和理想的频响相乘，从而得到时域逼近：
$$h(k)=h_{id}(k)w(k)$$
窗函数应该符合以下标准:

① 主瓣宽度应尽量窄，以期获得较陡峭的过渡带特性。

② 旁瓣应尽量少，且其幅度与主瓣相比应尽可能小。也就是说，窗口频谱的能量应尽量集中于主瓣，借以减少肩峰与余振，从而增大阻带衰减。

常见窗函数比较如下：

1、调用函数

signal库提供的firwin()采用此算法，它的调用形式如下

scipy.signal.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None)

• 输入参数
  重点关注numtaps，cutoff，window，pass_zero四个参数：
  numtaps——滤波器点数
  cutoff——截至频率
  window——窗函数，可取’hamming’,'hanning’等
  pass_zero——{True, False, ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’}就是过不过0，用来确定是低通还是高通,带通还是带阻，默认为True，当然也可以使用具体的‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’来进行确定。

• 输出参数

  b——滤波器系数(a都等于1)

2、具体调用方式

• 低通

signal.firwin(numtaps, f)

• 高通

signal.firwin(numtaps, f, pass_zero=False)

• 带通

signal.firwin(numtaps, [f1, f2], pass_zero=False)

• 带阻

signal.firwin(numtaps, [f1, f2])

• 多带通:

# (带通范围为[0, f1], [f2, f3] 和 [f4, 1])
signal.firwin(numtaps, [f1, f2, f3, f4])
# (带通范围为[f1, f2] 和 [f3, f4])
signal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)

3、频响曲线

调用完成以后，会返回滤波器参数，然后可以使用freqz计算数字滤波器的频率响应，它的调用方式如下：

freqz(b, a=1, worN=None, whole=0, plot=None)

其中b和a是滤波器的系数，worN为所计算的频率点数，whole为0表示计算频率的上限为pi，whole为1表示计算频率的上限为2*pi。

它返回一个组元 (w,h) ，其中w为所有计算了响应的频率数组，其值为正规化的圆频率。h是滤波器系统在每个对应的频率点的响应。复数的幅值表示滤波器的增益特性，相角表示滤波器的相位特性。

4、实际例子

设计一个数字FIR低通滤波器，其技术指标为
$${\omega }_p=0.2\pi ,{\alpha }_p=0.25dB；{\omega }_s=0.4\pi ,{\alpha }_s=50dB$$

解：

阻带最小衰减为50 dB，查前面的表可知只有汉明窗和布莱克曼窗可提供大于50dB的衰减，故选汉明窗为窗函数，它提供较小的过渡带，因此具有较小的阶数。
$$\Delta w=w_s-w_p=0.2\pi =\frac{6.6\pi }{N}$$
可获得滤波器阶数N。具体代码如下：

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


wp = 0.2 * np.pi
ws = 0.4 * np.pi
deltaw = ws - wp
N = int(np.ceil(6.6 * np.pi / deltaw) + 1)
b = signal.firwin(N, wp / np.pi, window="hamming")
w, h = signal.freqz(b)
plt.title("Digital filter frequency response")

plt.plot(w, 20 * np.log10(np.abs(h)))
plt.ylabel("Amplitude Response (dB)")
plt.xlabel("Frequency (rad/sample)")
plt.grid()
plt.show()