经典滤波器的设计原理

前言

想写一篇文章，包含经典滤波器的万象。

1、滤波器基础知识

总的来说，滤波器可分为经典滤波器和现代滤波器两大类。经典滤波器是假定输入信号x(n)中的有用成分和希望去除的成分各自占有不同的频带。这样，当x(n)通过一个线性系统（滤波器）之后可以将欲去除的成分有效去除。如果信号和噪声的频谱相互重叠，那么经典滤波器将无能为力。

现代滤波器的主要内容是从含有噪声的数据记录中估计出信号的某些特征或信号本身。一旦信号被估计出，那么利用它们的统计特征，如自相关函数、功率谱等，导出一套最佳的估值算法，所得到的估计信号的信噪比会比原信号更高。典型的现代滤波器是维纳滤波器。

对于常用的数字滤波器，一般分为IIR滤波器和FIR滤波器。
IIR滤波器的转移函数为
$$H(z)=\frac{{\sum }_{r=0}^M{b}_r{z}^{-r}}{1+{\sum }_{k=1}^N{a}_k{z}^{-k}}$$
FIR滤波器的转移函数为
$$H(z)=\sum _{n=0}^{N-1}h(n)z^{-n}$$
这两者滤波器无论是在性能还是设计方法上都有很大区别。FIR滤波器可以对给定频率特性直接进行设计，而IIR滤波器目前最通用的方法是利用已经很成熟的模拟滤波器的设计方法来进行设计。而模拟滤波器的设计方法有ButterWorth滤波器，Chebyshev滤波器。本文讨论的主要是经典滤波器的设计。

2、FIR滤波器和IIR滤波器的选择思路

这两类的滤波器特点如下所示：
1）FIR单位冲击响应为有限长，而IIR单位冲激响应为无限长。
2）从滤波器结构看，IIR滤波器为递归结构，即系统存在反馈环节，受有限字长效应影响较大，可能产生震荡。而FIR滤波器采用非递归结构，无反馈环节，系统只有零值极点，系统永远稳定，有限字长效应很小，运算误差较小。
3）从成本来看，IIR滤波器的极点可位于单位圆内任意一点，因此可用较低的阶数获得较高的选择性。由于阶数低，所需要的的储存单元较少，比较经济。而在相同指标条件下，FIR滤波器的阶数远高于IIR滤波器。其原因是FIR滤波器只有原点处的极点，智能用高阶数实现高频率选择性。通常实现同样指标，FIR滤波器需要的阶数是IIR滤波器的5-10倍。
4）从性能上看，IIR滤波器的相位特性一般是非线性，而FIR滤波器的相位特性则往往是线性。
5）从设计上看，IIR滤波器可借助模拟滤波器的计算结果，而FIR滤波器没有现成的设计公式。
6）IIR滤波器受模拟滤波器类型的约束，只适合设计具有片段性特性的滤波器，如低通、高通、带通、带阻等。而FIR滤波器除这些外，还可设计具有某些特殊用途的滤波器，如微分器、正交变换器等。
7） 从使用要求看，IIR适用于对相位要求不敏感的情况，如语音、通信等，而FIR可用于对线性相位特性要求高的地方，如图像信号处理等。
综上所述，没有哪一种类型的滤波器在所有情况下都是最好的，必须根据实际情况选择合适的滤波器。有两点比较有用的经验是：
1) 当过度带宽这个指标非常重要且要求苛刻，也就是滤波器要求锐截止时，首选IIR滤波器，因为此时FIR滤波器所需的系数非常多，效率很低。
2）当相位的线性要求非常高时，首选FIR滤波器，因为即便增加相位补偿处理，IIR滤波器在边沿仍有比较明显的非线性。

3、滤波器的设计思路

不论是FIR滤波器还是IIR滤波器的设计都包括以下3个步骤：
1、给出所需要的滤波器的技术指标；
2、设计一个H(z)使其逼近所需要的技术指标；
3、实现所设计的H(z)。

前面已经讲过，IIR滤波器设计最通用的方法是借助于模拟滤波器的设计方法。设计步骤如下：
1）按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的设计指标，注意，是模拟低通滤波器，低通。
2) 根据转换后的指标设计模拟低通滤波器G(s)
3）按一定的规则将G(s)转换为H(z)；
若设计的数字滤波器本身就是低通的，则前面三步即可。如果所设计的数字滤波器是高通、带通、带阻滤波器，则还有步骤4.
4）将高通、带通、带阻数字滤波器的技术指标先转换为低通模拟滤波器的技术指标，然后再重复2-4步骤。

4、由低通滤波器转换到高通、带通、带阻

用$H_L(z)$ 表示低通滤波器的传递函数，通道边沿频率为$w_{lp}$，阻带边沿频率为$w_{ls}$，则高通滤波器的传递函数$H_H(z)$可写为
$$H_H(z)=1-H_L(z)$$
，这样，所得到的高通滤波器的通带边沿频率及阻带边沿频率为
$w_{hp}=w_{ls}$,$w_{hs}=w_{lp}$。
其实思路非常好理解，相当于把原始信号减去进行完低通滤波后的信号，不就得到了高通滤波的信号了。
有了这个思路之后，用两个低通滤波器就可以得到带通滤波器，和带阻滤波器了，思路一模一样。

5、IIR滤波器的设计

5.1 IIR滤波器的主要特征

IIR的单位冲激响应的$h(n)$用数学公式表示如下：
$$h(n)\ne 0,0\le n<\infty$$
输入输出关系用差分方程描述为
$$\begin{gathered} y(n)=\sum _0^{\infty }h(m)x(n-m) \\ =-\sum _{m=1}^{M-1}{a}_{m}y(n-m)+\sum _{m=0}^{N-1}{b}_{m}x(n-m) \end{gathered}$$
$b_m$和$a_m$分别称为滤波器的系数，利用Z变换，可得到IIR滤波器的传递函数H(z)
$$H(z)=\frac{Y(z)}{X(z)}=\frac{{\sum }_{m=0}^{N-1}{b}_m{z}^{-m}}{1+{\sum }_{m=1}^{M-1}{a}_m{z}^{-m}}$$

5.2模拟低通滤波器的设计

5.2.1概述

给定模拟低通滤波器的技术指标${\alpha }_p,{\Omega }_p,{\alpha }_s,{\Omega }_s$，分别为通带允许的最大衰减、通带上限角频率，阻带应达到的最小衰减，阻带下限角频率。

插一下，有关于模拟频率、数字频率、角频率的区别看这里。数字频率与模拟频率的区别

现在，希望设计一个低通滤波器$G(s)$，
$$G(s)=\frac{d_0+d_{1}s+...+d_{N-1}{s}^{N-1}+d_N{s}^N}{c_0+c_{1}s+...+c_{N-1}{s}^{N-1}+c_N{s}^N}$$
使得其对数幅频响应$10lg|G(j\Omega ){|}^2$在${\Omega }_p$和${\Omega }_s$处分别达到${\alpha }_p$和${\alpha }_s$的要求。

我们知道，${\alpha }_p$和${\alpha }_s$都是$\Omega$的函数，它们的大小取决于$|G(j\Omega )|$的形状。因此，我们定义一个衰减函数$\alpha (\Omega )$，即
$$\alpha (\Omega )=10lg|\frac{X(j\Omega )}{Y(j\Omega )}{|}^2=10lg\frac{1}{|G(j\Omega ){|}^2}$$

显然，
$$\begin{gathered} {\alpha }_p=\alpha ({\Omega }_p)=-10lg|G(j{\Omega }_p){|}^2, \\ {\alpha }_s=\alpha ({\Omega }_s)=-10lg|G(j{\Omega }_s){|}^2 \end{gathered}$$
这样的话，通过衰减函数 $\alpha (\Omega )$,就把低通模拟滤波器的四个技术指标和滤波器的幅平方特性联系了起来。
这里只介绍一下butterworth滤波器的幅平方特性的表达式。
$$|G(j\Omega {)}^2|=\frac{1}{1+C^2({\Omega }^2{)}^N}$$
其中，C为待定常数，N为待定的滤波器的阶次。
Butterworth滤波器也是最简单的滤波器。
由于每一个滤波器的频率范围将取决于设计者所应用的目的，因此必然是千差万别。为了使设计规范化，我们需要将滤波器的频率参数进行归一化处理。设所给的实际频率为$\Omega$，归一化后的频率为$\lambda$，对低通模拟滤波器，令
$$\lambda =\Omega /{\Omega }_p$$
显然，${\lambda }_p=1$，${\lambda }_s={\Omega }_s/{\Omega }_p$

5.2.2butterworth 模拟低通滤波器的设计

butterworth模拟低通滤波器的设计可按下面三个步骤来进行，

1)将实际频率$\Omega$归一化

得到
$$|G(j\lambda ){|}^2=\frac{1}{1+C^2{\lambda }^{2N}}$$
由此可以看出，在$|G(j\lambda ){|}^2$中只有两个参数C和N，N是滤波器的阶次。

2)求出C和N

将归一化后的频率$\lambda$带入衰减函数，我们得到
$$\alpha (\lambda )=10lg(1+C^2{\lambda }^{2N})$$
进一步，我们得到
$$C^2{\lambda }^{2N}=10^{\alpha (\lambda )/10}-1$$
我们再将${\lambda }_p$和${\lambda }_s$带入，得到
$$\begin{gathered} C^2{\lambda }_p^{2N}=10^{{\alpha }_p/10}-1 \\ {C}^2{\lambda }_s^{2N}=10^{{\alpha }_s/10}-1 \end{gathered}$$
这个时候，我们不要忘记了前面有一条非常好用的东西，就是${\lambda }_p=1$。于是，我们便可以解出C和N
$$\begin{gathered} C^2=10^{{\alpha }_p/10}-1 \\ N=lg\sqrt{\frac{10^{{\alpha }_s/10}-1}{10^{{\alpha }_p/10}-1}}/lg({\lambda }_s) \end{gathered}$$
这样，我们便求出了C和N。
如果令${\alpha }_p=3dB$，这时C=1。这样butterworth滤波器的设计就只剩下一个参数N，这时
$$|G(j\lambda ){|}^2=\frac{1}{1+{\lambda }^{2N}}=\frac{1}{1+(\Omega /{\Omega }_p{)}^{2N}}$$
可以利用这个式子，简单讨论一下butterworth滤波器响应的一些特点。
（1）当$\Omega =0$，$\lambda =0$，$|G(j\lambda ){|}^2=1$，$\alpha (0)=0$，即在$\Omega =0$处无衰减。
（2）当$\Omega ={\Omega }_p$，${\lambda }_p=1$，$|G(j\lambda ){|}^2=0.5$，$|G(j\lambda )|=0.707$，即${\alpha }_p=3dB$
(3) 当$\lambda$由0增加到1时，$|G(j\lambda ){|}^2$单调减小，N越大，$|G(j\lambda ){|}^2$减小的越慢，即在通带内$|G(j\lambda ){|}^2$越平缓。
（4）当$\Omega >{\Omega }_p$，即${\lambda }_p>1$时，$|G(j\lambda ){|}^2$也是单调减小。但是，这个时候由于$\lambda >1$，因此，N越大，衰减速度越大。
如下图所示。

(3)确定G(s)

与归一化频率$\lambda =\Omega /{\Omega }_p$相对应，归一化复数变量为p，$p=j\lambda =j\Omega /{\Omega }_p=s/{\Omega }_p$,代入得到
$$G(p)G(-p)=\frac{1}{1+(p/j{)}^{2N}}=\frac{1}{1+(-1{)}^N{p}^{2N}}$$
要确定G(s)，我们需要确定G(s)的极点所在。令
$$1+(-1{)}^N{p}^{2N}=0$$
解得，
$$p_k=exp(j\frac{2k+N-1}{2N}\pi ),k=1,2,...2N$$

5.2.3 例题

试设计一个模拟低通butterworth滤波器，要求截止频率$f_p=5000Hz$，通带最大衰减为${\alpha }_p=3dB$，阻带起始频率为$f_s=10000Hz$，阻带最小衰减${\alpha }_s=30dB$。
解答：
首先，我们将频率归一化，${\Omega }_p=2\pi 5000Hz,{\lambda }_p=1,{\lambda }_s=2$。

5.3 模拟滤波器到数字滤波器

将模拟滤波器映射到数字滤波器的方法有2种，分别是冲激不变法和双线性变换法。

5.3.1冲激不变法

冲激不变法本质上是一种将模拟滤波器转换成数字滤波器的方法，它的基本思想是对给定的模拟滤波器$H_a(s)$进行LapLace反变换，得到滤波器的冲激响应$h_a(t)$，然后对$h_a(t)$进行采样得到h(n)，再对得到的h(n)进行Z变换，得到数字滤波器的传递函数H(z)，从而得到滤波器的系数。基本实现框图如图所示：

假设模拟滤波器只有单阶极点，那么可以用部分分式分解的方法，将$H_a(s)$分解为：
$$H_a(s)=\sum _{i=1}^N\frac{R_i}{s-s_i}$$
$R_i$表示第i个极点的系数。得到反变换后的单位冲激响应$h_a(t)$为：
$$h_a(t)=\sum _{i=1}^N{R}_i{e}^{s_{i}t}u(t)$$
这步的推导过程用到了$e^{at}$的LapLace变换为$\frac{1}{s-a}$.
对$h_a(t)$进行采样，可得到数字滤波器的单位冲激响应h(n)为
$$h(n)=h_a(n{T}_s)=\sum _{i=1}^N{R}_i{e}^{s_{i}n{T}_s}u(n{T}_s)$$
其中$T_s$为采样周期。对h(n)再进行Z变换，得到最后的数字滤波器的传递函数H(z)为
$$H(z)=\sum _{R_i}^{1-e^{s_i{T}_s}}$$

6、FIR滤波器的设计

FIR滤波器是指单位冲激响应长度有限的滤波器。主要特征包括冲激响应长度有限、传递函数不包含极点、常用卷积来实现以及可保证严格的线性相位4个主要特征。

6.1FIR滤波器的基本特征

单位冲激响应$h(n)$定义为
$$h(n)\ne 0,0\le n\le N-1$$
长度为N的FIR滤波器零点个数为N-1。N-1也为FIR滤波器的阶数。
h(n)的Z变换为
$$\begin{gathered} H(z)=\sum _{n=0}^{N-1}{z}^{-n} \\ =h(0)+h(1)z^{-1}+h(2)z^{-2}+...+h(N-1)z^{-(N-1)} \end{gathered}$$
我们可以知道，对FIR滤波器来说，单位冲激响应h(n),n=0,1,2,…N-1也即是滤波器的系数。因此，将输入及输出的差分方程可描述为：
$$\begin{gathered} y(n)=h(0)x(n)+h(1)x(n-1)+...+h(N-1)x(n-(N-1)) \\ y(n)=\sum _{m=0}^{N-1}h(m)x(n-m)=h(n)\ast x(n) \end{gathered}$$
我们可以看出，FIR滤波器当前时刻的输出仅与当前及之前时刻的输入有关。

6.2 FIR滤波器系数的计算方法

最常用的方法为窗函数法、频率采样法以及最优化方法。其中以最优化方法应用最为广泛。

参考资料

1、《数字信号处理 理论、算法与实现》 胡广书 著
2、《学以致用：深入浅出数字信号处理》 江志宏著