#! https://zhuanlan.zhihu.com/p/673537350

FIR滤波器设计–窗函数法

引言

不像IIR滤波器利用经典模拟滤波器的设计，FIR滤波器设计是直接逼近离散时间系统理想的频率响应或单位抽样响应。

最简单的FIR滤波器设计方法是窗函数法。窗函数法的基本思想是：先设计一个理想线性相位滤波器，然后用一个窗函数对理想线性相位滤波器的单位抽样响应进行截断，从而得到一个有限长的单位抽样响应。

FIR窗函数设计法过程

我们以数字低通线性相位的FIR滤波器设计为例，介绍窗函数法的设计过程。

理想线性相位数字低通滤波器的频率响应如下：

\[\begin{split} H_{d}(e^{j\omega})= \begin{cases} e^{-j\omega \alpha}, \qquad |\omega| \leq \omega_c,\\ 0, \qquad \omega_c<|\omega|\leq \pi\\ \end{cases} \end{split}\]

我们画出其频率响应图如下：

理想线性相位数字低通滤波器的频率响应

其单位抽样响应\(h_d[n]\)为：

\(\displaystyle h_d[n] = \frac{1}{2 \pi}\int_{-\omega_c}^{\omega_c}e^{-j\omega \alpha}e^{j\omega n} d\omega = \frac{\omega_c}{\pi}\frac{\sin(\omega_c (n-\alpha)}{\omega_c (n-\alpha)}\)

同样地，我们画出了单位抽样响应的时域图如下：

理想线性相位数字低通滤波器的单位抽样响应

我们可以看到，\(h_d[n]\)是中心点在\(\alpha\)的偶对称无限长的非因果序列。要获得有限长的，因果的，线性相位的单位抽样响应，我们需要对\(h_d[n]\)进行截断。截断的方法是让\(h_d[n]\)乘上一个窗函数\(w[n]\)，从而得到一个有限长的序列\(h[n]\)，这样一个过程也叫加窗。

\(\displaystyle h[n] = h_d[n]w[n], \quad n=0,1,\cdots, N-1\)

\(\displaystyle \alpha = \frac{N-1}{2}\)

其中，\(N\)为截断的长度。窗函数\(w[n]\)的长度为\(N\)。

我们首先假定窗函数\(w[n]\)是一个长度为\(N\)的矩形窗函数，即：

\(\displaystyle w[n] = \begin{cases} 1, \qquad n=0,1,\cdots, N-1\\ 0, \qquad \text{其他} \end{cases}\)

则加窗后，\(h[n]\)的时域图如下：

时域的截断在频域上等同于两个信号频率的周期卷积。因此，截断后系统的频率响应为：

\(\displaystyle H(e^{j\omega}) = \frac{1}{2 \pi}\int_{-\pi}^{\pi}H_d(e^{j\theta})W(e^{j(\omega-\theta)})d\theta\)

我们可以看到，窗函数法设计的FIR滤波器的频率响应等于理想线性相位滤波器的频率响应与窗函数的频率响应的周期卷积。因此，窗函数的选择对FIR滤波器的频率响应有很大的影响。

矩形窗的频率响应为：

\(\displaystyle W_R(e^{j\omega}) = \sum_{n=0}^{N-1}e^{-j\omega n} = e^{-j\omega(\frac{N-1}{2})}\frac{\sin(\frac{\omega N}{2})}{\sin{\frac{\omega}{2}}}\)

矩形窗函的相位也是线性的，其幅度函数可以表示如下：

\(\displaystyle W_R(\omega) = \frac{\sin(\frac{\omega N}{2})}{\sin{\frac{\omega}{2}}}\)

幅度函数的图形如下：

矩形窗的幅度函数

窗函数的主瓣是指窗函数在频域中的主要能量集中区域，以\(\omega=0\)为中心，左右第1过零点之间的频谱称为主瓣，如上图所示，矩形窗的主瓣宽度为\(\frac{4\pi}{N}\)。主瓣之外的频谱称为旁瓣，第1过零点与第2过零点的频谱称为第一旁瓣，第2过零点与第3过零点的频谱称为第二旁瓣，旁瓣的能量越来越小。

将理想线性相位低通滤波器频率响应写成：

\(\displaystyle H_d(e^{j\omega})=H_d(\omega)e^{-j(\frac{N-1}{2})\omega}\)

则其幅度函数为

\[\begin{split} H_{d}(e^{j\omega})= \begin{cases} 1, \qquad |\omega| \leq \omega_c,\\ 0, \qquad \omega_c<|\omega|\leq \pi\\ \end{cases} \end{split}\]

因此窗函数法实际设计出来滤波器的幅度函数表示为：

\(\displaystyle H(\omega)= H_d(\omega) \ast W_R(\omega)\)

圆周卷积过程用下面的示意图表示，我们给出几个关键点的卷积的结果，这些结果相当于图形阴影部分的面积。

最后我们得出滤波器如下图：

  • 从上图可以看到滤波器波峰的位置在\(\omega_c-\frac{2\pi}{N}\)处，波谷的位置在\(\omega_c+\frac{2\pi}{N}\)，因此波峰和波谷的间距刚好是矩形窗的主瓣宽度\(\frac{4\pi}{N}\)，过渡带宽在波峰和波谷之间。如果要减小过渡带宽，可以增加\(N\)。

  • 波峰和波谷振幅的相对值不会随着\(N\)的增大而改变，也就是说滤波器阻带的衰减固定了，如果要改变阻带的衰减，只能更换衰减更快的窗函数。

FIR窗函数设计法与Gibbs现象的联系

FIR窗函数法设计的核心思想是用窗函数截断理想线性相位滤波器的单位抽样响应，其实际的滤波器的频率响应表示如下：

\(\displaystyle H(e^{j\omega}) = \sum_{n=0}^{N-1}\frac{\sin(\omega_c (n-\alpha))}{\pi (n-\alpha)}w[n]e^{-jn\omega}\)

理想的数字低通滤波器用复指数序列表示如下：

\(\displaystyle H_M(e^{j\omega}) = \sum_{n=-M}^{M}\frac{\sin\omega_c n}{\pi n}e^{-j\omega n}\)

假设窗函数为矩形窗，上面两个表达式除了在时域上时延不同外，其他都一样。因此，对理想线性相位滤波器的单位抽样响应的截断过程，等同于傅里叶级数的收敛问题。傅里叶级数收敛过程中会产生吉布斯现象（Gibbs Phenomenon）， 吉布斯现象的示意图如下：

动态图表示的吉布斯现象如下：

吉布斯现象表明了在跳跃的不连续点上出现振荡，随着\(M\)的增大，振荡向跳跃点靠近，但是振荡的幅度不会减小。这些分析和我们前面滤波器设计的分析结果一致，随着滤波器阶次的增大，振荡幅度不会变化，振荡向理想滤波器截止频率处靠近，这意味着过渡带宽减少。

窗函数

从上面的分析中可以得知，为了加大阻带衰减，只能使用衰减更快的窗函数。本节给出了几种常用的窗函数。

三角形窗

\(\displaystyle w[n] = \begin{cases} \frac{2n}{N-1}, \qquad 0 \leq n \leq \frac{N-1}{2}\\ 2-\frac{2n}{N-1}, \qquad \frac{N-1}{2} \leq n \leq N-1 \end{cases}\)

汉宁窗（Hanning）（又称升余弦窗）

\(\displaystyle w[n] = \frac{1}{2}[1-\cos(\frac{2\pi n}{N-1})]R_N[n]\)

海明窗（Hamming）（又称改进的升余弦窗）

\(\displaystyle w[n] = [0.54-0.46\cos(\frac{2\pi n}{N-1})]R_N[n]\)

布拉克曼（Blackman）窗（又称二阶升余弦窗）

\(\displaystyle w[n] = [0.42-0.5\cos(\frac{2\pi n}{N-1})+0.08\cos(\frac{4\pi n}{N-1})]R_N[n]\)

各种窗函数的时域和频域图表示如下：

各种窗函数基本参数的比较如下：

窗函数

旁瓣峰值（dB）

主瓣宽度

阻带最小衰减（dB）

过渡带宽

矩形窗

-13

2

-21

0.9

三角形窗

-25

4

-25

2.1

汉宁窗

-31

4

-44

3.1

海明窗

-25

4

-53

3.3

布拉克曼窗

-57

6

-74

5.5

  • 主瓣宽度和过渡带宽的单位为\(\frac{2\pi}{N}\)，其中\(N\)为窗函数的长度。

设计技巧：如果要保证衰减的最小值，可以选择衰减更快的窗函数，但从上面的图形和表格可以看到，衰减更快的窗函数，在相同滤波器长度下，其过渡带宽更宽，有利的是过渡带宽也可以通过增加滤波器的阶次来调整。因此对于窗函数法设计来说，第一步是要根据衰减的要求选择窗函数，因为窗函数的形状唯一确定了最小的衰减值，第二步是根据过渡带宽的要求选择滤波器的长度。

例：用窗函数法设计一个线性相位的数字低通滤波器，其通带截止频率为\(\omega_p = 0.3\pi\)，阻带截至频率为\(\omega_{st} = 0.5\pi\)，阻带衰减为\(-40dB\)。

解：

  • 第一步：首先构造一个理想线性相位数字低通滤波器，其单位抽样响应为： \(\displaystyle h_d[n] = \begin{cases} \frac{\sin(\omega_c (n-\alpha))}{\pi(n-\alpha)},\quad n \neq \alpha \\ \frac{\omega_c}{\pi}, \quad n = \alpha \end{cases}\)

    其中，\(\alpha = \frac{N-1}{2}\)。 我们假定理想数字低通滤波器的截止频率在所设计的滤波器的通带截至频率和阻带截至频率之间，即\(\omega_c = \frac{\omega_p+\omega_{st}}{2}=0.4\pi\)。

  • 第二步：选择窗函数。根据阻带衰减的要求，我们选择接近于\(-40dB\)衰减的汉宁窗，通过查表我们可以知道汉宁窗的过渡带宽为\(\bigtriangleup \omega= \frac{6.2\pi}{N}\)，我们题目给出的过渡带宽要求为\(\bigtriangleup \omega = |\omega_{st}-\omega_p|= 0.2\pi\)。我们可以通过下面的公式计算出滤波器的阶次和时延：

    \(\displaystyle N = \frac{6.2}{0.2} = 31\)

    \(\alpha = \frac{N-1}{2} = 15\)

    于是我得出汉宁窗\(w[n]\)的表达式为：

    \(\displaystyle w[n] = \frac{1}{2}[1-\cos(\frac{2\pi n}{N-1})]R_N[n]= \frac{1}{2}[1-\cos(\frac{\pi n}{15})]R_{31}[n]\)

  • 第三步：计算所设计滤波器的单位抽样响应。根据上面的公式，我们可以计算出滤波器的单位抽样响应为：

    \(\displaystyle h[n] = h_d[n]w[n] = \begin{cases}\frac{\sin(0.4\pi(n-15)}{\pi(n-15)}\frac{1}{2}[1-\cos(\frac{\pi n}{15})]R_{31}[n],\quad n \neq 15 \\ 0.4,\quad n = 15 \end{cases}\)

  • 第四步：计算滤波器的频率响应\(H(e^{j\omega})=\sum_{n=0}^{30} h[n]e^{-j\omega n}\)。根据频率响应，检验各项指标是否满足要求。比如阻带衰减是否满足要求，过渡带宽是否满足要求等。我们画出了滤波器的频率响应图如下：

  • 第五步：如果不满足要求，可以通过增加滤波器的阶次来调整，即增大\(N\)。或者换衰减更快的窗函数。然后重复第四步，直到满足要求为止。 在第四步的图中，我们发现所设计出来的滤波器不满足给定的奇数指标要求，所以我们通过增加\(N\)来改善滤波器，这里我们取\(N=32\)，我们画出了滤波器的频率响应图如下：

    从图中可以看出设计出来的滤波器是满足给定的技术指标要求的。

例：用窗函数法设计一个线性相位的数字高通滤波器，其阻带截至频率为\(\omega_{st} = 0.3\pi\)，过渡带宽为\(\bigtriangleup \omega = 0.2\pi\), 阻带衰减为\(-50dB\)。

解：

  • 首先构造一个理想的线性相位的数字高通滤波器，其频率响应如下：

    \[\begin{split}H_{d}(e^{j\omega})= \begin{cases} e^{-j\omega \alpha}, \qquad \omega_c \leq |\omega| \leq \pi,\\ 0, \qquad |\omega|\leq \omega_c\\ \end{cases} \end{split}\]

    通过离散时间傅里叶反变换，求出其单位抽样响应为：

    \[\begin{split}\begin{align*} h_d[n] & = \frac{1}{2 \pi}(\int_{-\pi}^{-\omega_c}e^{-j\omega \alpha}e^{j\omega n} d\omega + \int_{\omega_c}^{\pi}e^{-j\omega \alpha}e^{j\omega n} d\omega) \\ & = \begin{cases} -\frac{sin((n-\alpha)\omega_c)}{\pi(n-\alpha)}, \quad n \neq \alpha\\ 1-\frac{\omega_c}{\pi}, \quad n = \alpha \end{cases} \end{align*}\end{split}\]

    其中\(\alpha = \frac{N-1}{2}\)。我们假定理想数字高通滤波器的截止频率在所设计的滤波器的通带截至频率和阻带截至频率之间，即\(\omega_c = \omega_p+\frac{\bigtriangleup \omega}{2}=0.4\pi\)。

  • 选择窗函数。根据阻带衰减的要求，我们选择接近于\(-50dB\)衰减的海明窗，通过查表我们可以知道海明窗的过渡带宽为\(\bigtriangleup \omega= \frac{6.6\pi}{N}\)，我们题目给出的过渡带宽要求为\(\bigtriangleup \omega= 0.2\pi\)。我们可以通过下面的公式计算出滤波器的长度和时延：

    \(\displaystyle N = \frac{6.6}{0.2} = 33\)

    \(\alpha = \frac{N-1}{2} = 16\)

    海明窗\(w[n]\)的时域表达式为：

    \(\displaystyle w[n] = [0.54-0.46\cos(\frac{2\pi n}{N-1})]R_N[n]= [0.54-0.46\cos(\frac{2\pi n}{32})]R_{33}[n]\)

  • 计算所设计的线性相位FIR高通滤波器的单位抽样响应。根据上面的公式，我们可以计算出滤波器的单位抽样响应为：

    \(\displaystyle h[n] = h_d[n]w[n] = \begin{cases}-\frac{\sin(0.4\pi(n-16)}{\pi(n-16)}[0.54-0.46\cos(\frac{2\pi n}{32})]R_{33}[n], \quad n \neq 16 \\ 0.6, \quad n = 16\end{cases}\)

  • 计算滤波器的频率响应。根据频率响应，检验各项指标是否满足要求。比如阻带衰减是否满足要求，过渡带宽是否满足要求等。这里我们借助于python中的scipy包画出了滤波器的频率响应图如下：

    从图中可以看出设计出来的滤波器的阻带衰减只有\(-45.45dB\)，不满足所给定的技术指标要求。我们需要通过调整滤波器的长度来达到技术指标要求。这里我们增加滤波器长度到\(N=35\)，注意对于高通滤波器来说，滤波器的长度不能为偶数，这是由之前我们的推导的幅度函数的性质决定的，如果\(N\)为偶数，那么幅度函数在\(\omega=\pi\)处的值为0，这是不符合高通滤波器的要求的。我们画出了增加长度后的滤波器的频率响应图如下：

    从图中可以看出，增加长度后的滤波器的阻带衰减为\(-52.9dB\)，满足所给定的技术指标要求。

例：用窗函数法设计一个线性相位FIR带通滤波器，其下阻带截至频率为\(f_{st1}=20kHz\)，上阻带截至频率为\(f_{st2}=60kHz\)，通带下截止频率为\(f_{p1}=30kHz\)，通带上截止频率为\(f_{p2}=50kHz\)，阻带衰减为\(-60dB\), 抽样频率\(f_s = 200kHz\)。

解：

  • 首先将模拟参数转化成数字参数。根据数字滤波器的频率变换公式，我们可以得到数字滤波器的通带截止频率为：

    \(\displaystyle \omega_{p1} = \frac{2\pi f_{p1}}{f_s} = 0.3\pi\)

    \(\displaystyle \omega_{p2} = \frac{2\pi f_{p2}}{f_s} = 0.5\pi\)

    数字滤波器的阻带截至频率为：

    \(\displaystyle \omega_{st1} = \frac{2\pi f_{st1}}{f_s} = 0.2\pi\)

    \(\displaystyle \omega_{st2} = \frac{2\pi f_{st2}}{f_s} = 0.6\pi\)

  • 构造一个理想线性相位数字带通滤波器，其频率响应如下：

    \[\begin{split} H_{d}(e^{j\omega})= \begin{cases} e^{-j\omega \alpha}, \qquad \omega_{c1} \leq |\omega| \leq \omega_{c2},\\ 0, \qquad 其它 \\ \end{cases} \end{split}\]

    通过离散时间傅里叶反变换，求出其单位抽样响应为：

    \[\begin{split}\begin{align*} h_d[n] & = \frac{1}{2 \pi}(\int_{-\omega_{c2}}^{-\omega_{c1}}e^{-j\omega \alpha}e^{j\omega n} d\omega + \int_{\omega_{c1}}^{\omega_{c2}}e^{-j\omega \alpha}e^{j\omega n} d\omega) \\ & = \begin{cases} \frac{1}{\pi(n-\alpha)}\{sin((n-\alpha)\omega_{c2})-sin((n-\alpha)\omega_{c1})\}, \quad n \neq \alpha\\ \frac{\omega_{c2}-\omega_{c1}}{\pi}, \quad n = \alpha\ (\alpha \text{为整数}) \end{cases} \end{align*}\end{split}\]

    其中\(\alpha = \frac{N-1}{2}\)。我们假定理想数字带通滤波器的截止频率在所设计的滤波器的通带截至频率和阻带截至频率之间，即\(\omega_{c1} = \frac{\omega_{p1}+\omega_{st1}}{2}=0.25\pi\)，\(\omega_{c2} = \frac{\omega_{p2}+\omega_{st2}}{2}=0.55\pi\)。

  • 选择窗函数。根据阻带衰减的要求，我们选择接近于\(-55dB\)衰减的布拉克曼窗，通过查表我们可以知道布拉克曼窗的过渡带宽为\(\bigtriangleup \omega= \frac{11\pi}{N}\)。我们题目给出的过渡带宽为\(0.1\pi\)，我们可以通过下面的公式计算出滤波器的长度和时延：

    \(\displaystyle N = \frac{11}{0.1} = 110\)

    \(\alpha = \frac{N-1}{2} = 54.5\)

    布拉克曼窗\(w[n]\)的时域表达式为：

    \(\displaystyle w[n] = [0.42-0.5\cos(\frac{2\pi n}{N-1})+0.08\cos(\frac{4\pi n}{N-1})]R_N[n]= [0.42-0.5\cos(\frac{\pi n}{54.5})+0.08\cos(\frac{\pi n}{27.25})]R_{110}[n]\)

  • 计算所设计的线性相位FIR带通滤波器的单位抽样响应。根据上面的公式，我们可以计算出滤波器的单位抽样响应为：

    \(\displaystyle h[n] = h_d[n]w[n] = \frac{(\sin(0.55\pi(n-54.5)-\sin(0.25\pi(n-54.5))}{\pi(n-54.5)}[0.42-0.5\cos(\frac{\pi n}{54.5})+0.08\cos(\frac{\pi n}{27.25})]R_{110}[n]\)

  • 计算滤波器的频率响应。根据频率响应，检验各项指标是否满足要求。比如阻带衰减是否满足要求，过渡带宽是否满足要求等。这里我们画出了滤波器的频率响应图如下：

从图中可以看出，设计出来的滤波器满足所给定的技术指标要求。

例：用窗函数法设计一个线性相位FIR带阻滤波器，其通带截止频率为\(\omega_{p1} = 0.2\pi\)，\(\omega_{p2} = 0.6\pi\)，阻带截止频率为\(\omega_{st1} = 0.3\pi\)，\(\omega_{st2} = 0.45\pi\)，阻带衰减为\(-50dB\)。

解：

  • 首先构造一个理想线性相位数字带阻滤波器，其频率响应如下：

    \[\begin{split} H_{d}(e^{j\omega})= \begin{cases} 0, \qquad \omega_{c1} \leq |\omega| \leq \omega_{c2},\\ e^{-j\omega \alpha}, \qquad 其它 \\ \end{cases} \end{split}\]

    通过离散时间傅里叶反变换，求出其单位抽样响应为：

    \[\begin{split}\begin{align*} h_d[n] & = \frac{1}{2 \pi}(\int_{-\pi}^{-\omega_{c2}}e^{-j\omega \alpha}e^{j\omega n} d\omega + \int_{-\omega_{c1}}^{\omega_{c1}}e^{-j\omega \alpha}e^{j\omega n} d\omega+\int_{\omega_{c2}}^{\pi}e^{-j\omega \alpha}e^{j\omega n} d\omega) \\ & = \begin{cases} \frac{1}{\pi(n-\alpha)}\{sin((n-\alpha)\pi)-sin((n-\alpha)\omega_{c2})+sin((n-\alpha)\omega_{c1})\}, \quad n \neq \alpha\\ 1-\frac{\omega_{c2}-\omega_{c1}}{\pi}, \quad n = \alpha \end{cases} \end{align*}\end{split}\]

    其中\(\alpha = \frac{N-1}{2}\)。我们假定理想数字带阻滤波器的截止频率在所设计的滤波器的通带截至频率和阻带截至频率之间，即\(\omega_{c1} = \frac{\omega_{p1}+\omega_{st1}}{2}=0.25\pi\)，\(\omega_{c2} = \frac{\omega_{p2}+\omega_{st2}}{2}=0.525\pi\)。

  • 选择窗函数。根据阻带衰减的要求，我们选择接近于\(-50dB\)衰减的海明窗，通过查表我们可以知道海明窗的过渡带宽为\(\bigtriangleup \omega= \frac{6.6\pi}{N}\)，我们题目给出的过渡带宽求为\(\bigtriangleup \omega_1=\omega_{st1}-\omega_{p1}= 0.1\pi\)和\(\bigtriangleup \omega_2=\omega_{p2}-\omega_{st2}= 0.15\pi\)。这里我们选择过渡带宽小的\(\bigtriangleup \omega_1\)，因为此时所要求的\(N\)会比较大。我们可以通过下面的公式计算出滤波器的长度：

    \(\displaystyle N = \frac{6.6}{0.1} = 66\)

    这里要特别注意的是，对于带阻滤波器，我们所要求滤波器的长度\(N\)必须为奇数，这是由之前我们的推导的幅度函数的性质决定的，如果\(N\)为偶数，那么幅度函数在\(\omega=\pi\)处的值为0，这是不符合带阻滤波器的要求的。因此我们取\(N=67\)

    \(\alpha = \frac{N-1}{2} = 33\)

    海明窗\(w[n]\)的时域表达式为：

    \(\displaystyle w[n] = [0.54-0.46\cos(\frac{2\pi n}{N-1})]R_N[n]= [0.54-0.46\cos(\frac{\pi n}{33})]R_{67}[n]\)

  • 计算所设计的线性相位FIR带通滤波器的单位抽样响应。根据上面的公式，我们可以计算出滤波器的单位抽样响应为：

    \(\displaystyle h[n] = h_d[n]w[n] = \begin{cases}\frac{1}{\pi(n-33)}\{sin((n-33)\pi)-sin((n-33)0.525\pi)+sin((n-\alpha)0.25\pi)\}[0.54-0.46\cos(\frac{\pi n}{33})]R_{67}[n], \quad n \neq 33 \\ 0.725, \quad n = 33\end{cases}\)

  • 计算滤波器的频率响应。根据频率响应，检验各项指标是否满足要求。比如阻带衰减是否满足要求，过渡带宽是否满足要求等。这里我们画出了滤波器的频率响应图如下：

从图中可以看出，设计出来的滤波器满足所给定的技术指标要求。

小结

从窗函数法的设计过程中，我们可以发现窗函数法不能精准的控制通带和阻带的衰减，在设计过程中，有时候需要不断利用软件累试找到满足要求的滤波器的长度。这样设计出来的滤波器可能会产生冗余，也就是说有些频率点的指标会远远超过给定的要求。在最优化的设计方法中，我们可以找到满足指标要求的阶次更低的滤波器