最优FIR滤波器设计

   1. 四种FIR滤波器对称结构
   设FIR滤波器的输入和输出序列分别为$x[n]$和$y[n]$，滤波器系数为${\omega }_m,0\le m<M$，则滤波过程可用如下差分方程表示
$$y[n]=\sum _{m=0}^{M-1}{\omega }_{m}x[n-m]\text{\hspace{1em}(1)}$$
  此外，可将输出序列写成系统冲激响应$h[n]$与输入序列的卷积，如下
$$y[n]=\sum _{m=0}^{M-1}h[m]x[n-m]\text{\hspace{1em}(2)}$$
  从上面的表达式可以看出，当以样点为单位进行延时时，FIR滤波器的抽头系数和其冲激响应本质上是一致的。
   FIR滤波器最吸引人的性质是可以实现线性相位(保证信号过滤波器后不失真)。具有线性相位的FIR滤波器的冲激响应满足如下对称条件：
$$h[n]=\pm h[M-1-n]\text{\hspace{1em}(4)}$$
  上面的表达式包含4种情况，分别称为1/2/3/4型FIR滤波器，对应关系如下：
$$\{\begin{array}{c}1型：h[n]=h[M-1-n],M为奇数\\ 2型：h[n]=h[M-1-n],M为偶数\\ 3型：h[n]=-h[M-1-n],M为奇数\\ 3型：h[n]=-h[M-1-n],M为偶数\end{array}$$
  下面以1型FIR滤波器为例，推导其频率响应的表达式，其他形式的FIR滤波器可按类似方式进行推导。
$$H(e^{j\omega })=\sum _{m=0}^{N-1}h[n]e^{-jn\omega }=(\sum _{n=0}^{(N-3)/2}h[n]e^{-jn\omega })+h[\frac{N-1}{2}]e^{-j(N-1)\omega /2}+(\sum _{n=(N+1)/2}^{N-1}h[n]e^{-jn\omega })$$
  由于$N$为奇数，所以上式以中间点为对称重新，将求和操作分成前面部分+中间部分+后面部分。对后面部分，令$m=N-1-n$，则${\sum }_{n=(N+1)/2}^{N-1}h[n]e^{-jn\omega }={\sum }_{m=0}^{(N-3)/2}h[N-1-m]e^{-j(N-1-m)\omega }={\sum }_{m=0}^{(N-3)/2}h[m]e^{-j(N-1-m)\omega }$，将其带入上式得到：
$$\begin{gathered} H(e^{j\omega })=(\sum _{n=0}^{(N-3)/2}h[n]e^{-jn\omega })+h[\frac{N-1}{2}]e^{-j(N-1)\omega /2}+(\sum _{n=0}^{(N-3)/2}h[n]e^{-j(N-1-n)\omega }) \\ =(\sum _{n=0}^{(N-3)/2}h[n](e^{-jn\omega }+e^{-j(N-1-n)\omega }))+h[\frac{N-1}{2}]e^{-j(N-1)\omega /2} \\ =e^{-j(N-1)\omega /2}[\sum _{n=0}^{(N-3)/2}h[n](e^{-j\omega [n-(N-1)/2]}+e^{j\omega [n-(N-1)/2]})+h[\frac{N-1}{2}]] \\ =e^{-j(N-1)\omega /2}[2\sum _{n=0}^{(N-3)/2}h[n]cos(\frac{N-1}{2}-n)\omega +h[\frac{N-1}{2}]] \end{gathered}$$
  上面的表达式中，令$m=(\frac{N-1}{2}-n$，则可进一步得到：
$$H(e^{j\omega })=e^{-j(N-1)\omega /2}[2\sum _{n=1}^{(N-1)/2}h[\frac{N-1}{2}-n]cos(n\omega )+h[\frac{N-1}{2}]]$$
  令
$$a[n]=\{\begin{array}{c}h[\frac{N-1}{2}],n=0\\ 2h[\frac{N-1}{2}-n],0<n\le \frac{N-1}{2}\end{array}$$
  则 $H(e^{j\omega })=e^{-j(N-1)\omega /2}{\sum }_{n=0}^{(N-1)/2}a[n]cos(n\omega )$。显然，此时$H(e^{j\omega })$具有线性相位。令
$$\begin{gathered} \varphi (\omega )=arg(H(e^{j\omega }))=-(N-1)\omega /2 \\ {H}_r(e^{j\omega })=\sum _{n=0}^{(N-1)/2}a[n]cos(n\omega ) \end{gathered}$$
  则$H(e^{j\omega })=e^{\varphi (\omega )}{H}_r(e^{j\omega })$。另外3种FIR滤波器的频响可以按类似的方式进行推导，直接将推导结果总结在下面：

表1. 4种FIR滤波器的相位和频响实部总结

滤波器类型	$\varphi (\omega )$	$H_r(e^{j\omega })$
1型	$-\frac{N-1}{2}\omega$	$2{\sum }_{n=1}^{(N-1)/2}h[\frac{N-1}{2}-n]cos(n\omega )+h[\frac{N-1}{2}]={\sum }_{n=0}^{(N-1)/2}a[n]cos(n\omega )$
2型	$-\frac{N-1}{2}\omega$	$2{\sum }_{n=0}^{(N-1)/2}h[n]cos\omega (\frac{N-1}{2}-n)={\sum }_{n=1}^{N/2}b[n]cos\omega (n-1/2)$ $b[n]=2h[\frac{N}{2}-n]$
3型	$\frac{\pi }{2}-\frac{N-1}{2}\omega$	$2{\sum }_{n=0}^{(N-3)/2}h[n]sin\omega (\frac{N-1}{2}-n)={\sum }_{n=1}^{(N-1)/2}c[n]sin\omega n$ $c[n]=2h[\frac{N-1}{2}-n]$
4型	$\frac{\pi }{2}-\frac{N-1}{2}\omega$	$2{\sum }_{n=0}^{N/2-1}h[n]sin\omega (\frac{N-1}{2}-n)={\sum }_{n=1}^{N/2}d[n]sin\omega (n-1/2)$ $d[n]=2h[\frac{N}{2}-n]$

  从上表中可见，3/4型滤波器相比1/2型滤波器需要增加额外的$\pi /2$的相移器，此外，对3/4型滤波器来说，$H_r(0)=0$，并不适用于低通滤波器的设计。一般情况下，我们设计时更多采用I型滤波器。
  基于上表中的表达式，我们期望将$H_r(e^{j\omega })$用统一的形式，即$H_r(\omega )=Q(\omega )P(\omega )$，下面直接给出推导结果

表2. 4种FIR滤波器频响实部分别因子总结

滤波器类型	$Q(\omega )$	$P(\omega )$
1型	$1$	${\sum }_{n=0}^{(N-1)/2}a[n]cos(\omega n)$
2型	$cos(\frac{\omega }{2})$	${\sum }_{n=0}^{N/2-1}\overline{b}[n]cos(\omega n)$ $\overline{b}[0]=\frac{1}{2}b[1]$ $\overline{b}[n]=2b[n]-b[n-1],n=1,...N/2$ $\overline{b}[\frac{N}{2}-1]=2b[\frac{N}{2}]$
3型	$sin(\omega )$	${\sum }_{n=0}^{(N-3)/2}\overline{c}[n]cos(\omega n)$ $\overline{c}[\frac{M-3}{2}]=c[\frac{M-1}{2}]$ $\overline{c}[\frac{M-5}{2}]=2c[\frac{M-3}{2}]$ $\overline{c}[n-1]-\overline{c}[n+1]=2c[n],2\le n\le \frac{M-5}{2}$ $\overline{c}[0]-\frac{1}{2}\overline{c}[2]=c[1]$
4型	$sin(\frac{\omega }{2})$	${\sum }_{n=0}^{N/2-1}\overline{d}[n]cos(\omega n)$ $\overline{d}[\frac{M}{2}-1]=2d[\frac{M}{2}]$ $\overline{d}[n-1]-\overline{d}[n]=2d[n],2\le n\le \frac{M}{2}-1$ $\overline{d}[0]-\frac{1}{2}\overline{d}[1]=d[1]$

   2. 最优滤波器设计
  在上一节中，我们知道四类FIR滤波器的频响实部可统一表示成如下形式：
$$H_r(\omega )=Q(\omega )P(\omega )=Q(\omega )(\sum _{n=0}^L\alpha (n)cos(\omega n))$$
  基于上述通用表达式，我们进行最优FIR滤波器设计。最优FIR滤波器设计的基本思路是在某种优化目标下，让设计的滤波器频响和理论频响的误差最小。由于上面4中类型的FIR滤波器的相频特性与滤波器系数无关，故优化过程无需考虑相频的作用。设$H_{dr}(\omega )$为期望的频响的实部，通常地，期望频响在通带内的增益为1，在阻带内的增益为0。此外，定义一个加权误差函数$W(\omega )$用于控制通带和阻带的波动。具体地
$$W(\omega )=\{\begin{array}{c}{\delta }_2/{\delta }_1,通带内\\ 1,阻带内\end{array}$$
  其中,${\delta }_1$,${\delta }_2$的定义如下图所示

  定义误差函数：
$$\begin{gathered} E(\omega )=W(\omega )[H_{dr}(\omega )-H_r(\omega )]=W(\omega )[H_{dr}(\omega )-Q(\omega )P(\omega )] \\ =W(\omega )Q(\omega )[\frac{H_{dr}(\omega )}{Q(\omega )}-P(\omega )]=\hat{W}(\omega )[\widehat{H_{dr}}(\omega )-P(\omega )] \end{gathered}$$
  我们期望找到一组系数$\alpha (n)$(根据具体滤波器类型，可对应到表中的$\overline{a}[n]、\overline{b}[n]、\overline{c}[n]、\overline{d}[n]$，并进一步对应到滤波器的冲激响应$h[n]$)，使得在要求的频带内$E(\omega )$的最大值可以最小化，所以优化问题可描述为
$$min[\underset{\omega \in S}{\max }|E(\omega )|]\text{\hspace{1em}(5)}$$
  其中，$S$表示最优化的频带集合，通常由通带和阻带组成。下面考虑对$P(\omega )$进行形式转换。对于三角函数有$cos(2\omega )=2cos\omega -1、cos(3\omega )=4co{s}^3\omega -3cos\omega ....$所以$P(\omega )={\sum }_{n=0}^L\alpha (n)cos(\omega n)={\sum }_{n=0}^L\beta (n)co{s}^n\omega$，即$P(\omega )$为$cos\omega$的$L$阶多项式，由于$cos\omega$在$0<\omega <\pi$区间上是单调函数，所以可以将$P(\omega )$按常规的$L$阶多项式进行处理。所以，优化问题转化为找到一组多项式系数使其满足（5）式的目标。这个问题即为切比雪夫最佳逼近问题。下面以一阶多项式逼近为例，直观说明该优化问题的基本思路。
  当多项式为1阶时，优化问题可写成
$$min[\underset{a,b\in [M,N]}{\max }|f(x)-ax-b|]\text{\hspace{1em}(6)}$$
  上面的优化问题可以用下图进行形象说明，我们期望期望找到一条直线，使与图中曲线的最大偏差最小。首先画两条平行线将曲线夹在其中，则期望的最佳直线即为这两条平行线的中线，最小误差即为这两条平行线与y轴截距的中点。

  首先我们定义偏差点的概念。偏差值$E$定义为满足(6)式的最大偏差的最小值。若$x_0\in [M,N]$使得$|f(x_0)-a{x}_0-b|=E$，则称$x_0$为正偏差点，若$x_0\in [M,N]$使得$|f(x_0)-a{x}_0-b|=-E$，则称$x_0$为负偏差点。由上图我们可以总结存在上述最佳逼近直线的充要条件：

01、$f(x)$为连续函数；
02、至少存在3个偏差点，且正负偏差点依次出现。

  不加证明地将上述1次多项式的情况推广到本文的$L$次多项式上，可知若存在$P(\omega )$为$\widehat{H_{dr}}(\omega )$的唯一逼近，则$E(\omega )$在$S$上至少存在$L+2$个偏差点，且这些偏差点正负依次出现。该结论被称为交错定理。
  对于函数$P(\omega )={\sum }_{n=0}^L\alpha (n)cos(\omega n)={\sum }_{n=0}^L\beta (n)co{s}^n\omega$，若$cos(\omega )$进行求导，显然在$0<\omega <\pi$区间内至多有$L-1$个导数值为0，对应$L-1$个偏差点。另外，$\omega =0$、$\omega =\pi$、$\omega ={\omega }_p$和$\omega ={\omega }_s$通常也可为偏差点，故至多存在$L+3$个偏差点。设这些偏差点对应的频率为${\omega }_n$，则根据交错定理，可得
$$\hat{W}({\omega }_n)[\widehat{H_{dr}}({\omega }_n)-P({\omega }_n)]=(-1{)}^n\delta ,n=0,1,...,L+1$$
  其中，$\delta$表示偏差值。将上面的表达式调整为
$$P({\omega }_n)+\frac{(-1{)}^n\delta }{\hat{W}({\omega }_n)}=\widehat{H_{dr}}({\omega }_n)\to \sum _{n=0}^L\alpha (n)cos({\omega }_{n}n)+\frac{(-1{)}^n\delta }{\hat{W}({\omega }_n)}=\widehat{H_{dr}}({\omega }_n)$$
  可将上式写成矩阵的形式:
$$\left[\begin{array}{ccccc}1& cos({\omega }_0)& ...& cos(L{\omega }_0)& \frac{1}{\hat{W}({\omega }_0)}\\ 1& cos({\omega }_1)& ...& cos(L{\omega }_1)& \frac{1}{\hat{W}({\omega }_1)}\\ ⋮& ...& ⋮& ⋮& \\ 1& cos({\omega }_{L+1})& ...& cos(L{\omega }_{L+1})& \frac{(-1{)}^{L+1}}{\hat{W}({\omega }_{L+1})}\end{array}\right]\left[\begin{array}{c}\alpha (0)\\ \alpha (1)\\ ⋮\\ \alpha (L)\\ \delta \end{array}\right]=\left[\begin{array}{c}\widehat{H_{dr}}({\omega }_0)\\ \widehat{H_{dr}}({\omega }_1)\\ ⋮\\ \widehat{H_{dr}}({\omega }_L)\\ \widehat{H_{dr}}({\omega }_{L+1})\end{array}\right]\text{\hspace{1em}(7)}$$
  在上面的矩阵中，$w_n$、$\delta$、$\alpha (n)$均为未知数，无法进行直接求解。一个直观的思路是先猜一组$w_n$，然后确定$P(\omega )$和$\delta$，最后计算$E(\omega )$，如结果不满足要求就重新选择新的偏差点频率，不断迭代求解，但上述求解过程中涉及矩阵求逆，计算量较大，故考虑用Remez算法进行求解。
  首先根据矩阵(7)可以求得$\delta$的表达式如下：
$$\delta =\frac{{\gamma }_0\widehat{H_{dr}}({\omega }_0)+{\gamma }_1\widehat{H_{dr}}({\omega }_1)+...+{\gamma }_{L+1}\widehat{H_{dr}}({\omega }_{L+1})}{\frac{{\gamma }_0}{\hat{W}({\omega }_0)}-\frac{{\gamma }_1}{\hat{W}({\omega }_1)}+...+\frac{(-1{)}^{L+1}{\gamma }_{L+1}}{\hat{W}({\omega }_{L+1})}}\text{\hspace{1em}(8)}$$
其中
$${\gamma }_k=\prod _{n=0,n\ne k}^{L+1}\frac{1}{cos{\omega }_k-cos{\omega }_n}$$
  显然，若已知所有极值点频率，则$\delta$是可求的。到这里，我们再来重新梳理一下我们想干什么事，又已知什么东西。
  首先，从上面的分析我们知道如果目标函数$min[\underset{\omega \in S}{\max }|E(\omega )|]$有唯一解，那么一定至少存在$L+2$个频点满足$\hat{W}({\omega }_n)[\widehat{H_{dr}}({\omega }_n)-P({\omega }_n)]=(-1{)}^n\delta$，其中$\delta$表示最大误差。但问题是，我们不知道${\omega }_n$具体是多少，所以只能随便猜几个。如果运气爆棚，一次就猜准了，那么根据(8)式算出来的$\delta$就是真实的我们想要的最大误差，此时对于所有$\omega \in S$(注意$\omega$是连续的)都会满足$|E(\omega )|\le \delta$。但是如果极值频率点选择错误，计算出来的${\delta }^{\prime }$将比真实$\delta$小（因为真实$\delta$本身表示的就是最大误差）。所以这个时候误差$|E(\omega )|$曲线上一定存在一些局部峰值点使得$|E(\omega )|>\delta$。基于这个分析，我们就可以不断去迭代极值频率点，直到找到所有正确的极值频率点。

图1. 误差曲线示意图(图中红色w表示正确的极值频率点，黑色w表示随机选择的错误的极值频率点，频率点选择错误计算出来的delta比真实delta小)

  将整体操作步骤总结如下：
（1）设置滤波器类型和阶数$N$，并根据表2计算$L$。设置通带/阻带纹波大小${\delta }_1$/${\delta }_2$，设置理想频响$H_{dr}(\omega )$；
（2）基于$L+2$个极值频率点${\omega }_0,...{\omega }_{L+1}$(第一次随机选择)，并根据(8)计算$\delta$；
（3）对$P({\omega }_n)$进行插值(算法中用的是拉格朗日插值) 得到$P(\omega )$，并计算$E(\omega )=\hat{W}(\omega )[\widehat{H_{dr}}(\omega )-P(\omega )]$，显然插值的目的是用密集采样点来逼近连续的情况，一般插值到$16N$点就足够；
（4）找出$E(\omega )$中满足$|E(\omega )|\ge \delta$的全部局部极值点(特别注意这边找的是极值点，不是一般的点，因为对于所有预设频率均能满足$|E({\omega }_n)|=\delta$，但这些预设频率点不一定是极值点)；
（5）如果找到满足步骤(4)要求的极值点，且找到的极值频率点与预设的极值频率点不一致，则将预设极值频率点替换成找到的极值频率点，若找到的极值频率点数大于$L+2$个，则取前$L+2$个极值频率点替换预设极值频率点。重新跳转到步骤(2)进行迭代处理。
（6）当找到的极值频率点数为$L+2$且与预设极值频率点一致时，跳出迭代过程。
（7）当找到正确的极值频率点后，步骤(3)插值得到的$P(\omega )$即为正确的，据此可得所设计的滤波器频响的实部为$H_r(\omega )=Q(\omega )P(\omega )$，结合其线性相位的虚部，即可得到完整的滤波器频响。

---

【附录】

1. 拉格朗日插值