数字滤波器

前言

  本文主要介绍数字滤波器的设计原理

---

一、数字滤波器

  过去用模拟滤波器（如用R、L、C元件）完成的事，现在可以使用数字方法实现。利用计算机完成数字相加、乘以常数和延时，将输入序列按既定的要求转换为输出序列。设线性时不变离散系统的输入信号为$f(n)$，其零状态响应为$y(n)$，则响应与输入的Z变换之比为系统函数：
$$H(z)=\frac{Y(z)}{F(z)}$$
  设N阶LTI系统的差分方程为
$$\sum _{k=0}^N{a}_{k}y(n-k)=\sum _{r=0}^M{b}_{r}f(n-r)$$
  输入为因果信号，在零状态下，对其取Z变换，从而有系统函数
$$H(z)=\frac{Y(z)}{F(z)}=\frac{\sum _{r=0}^M{b}_r{z}^{-r}}{\sum _{k=0}^N{a}_k{z}^{-k}}$$
  则数字滤波器就是求出一组系数$a_k$和$b_r$，使得滤波器具有所需特性，按单位响应样式可分为两类：IIR型和FIR型

1.1 IIR滤波器

  若系统函数满足
$$H(z)=\frac{\sum _{r=0}^M{b}_r{z}^{-r}}{1+\sum _{k=1}^N{a}_k{z}^{-k}}$$
  则称系统的脉冲响应$h(n)$是无限长的，该系统的差分方程为
$$\begin{gathered} y(n)=b_{0}f(n)+b_{1}f(n-1)+\cdots +b_{M}f(n-M) \\ -a_{1}y(n-1)-\cdots -a_{N}y(n-N) \end{gathered}$$
  也就是说IIR系统的输出$y(n)$不仅取决于输入值，还取决于输出值$y(n-1)$，$y(n-2)$,…这样的滤波器又称递归型滤波器

1.2 FIR滤波器

  若系统函数满足
$$H(z)=1+b_1{z}^{-1}+b_2{z}^{-2}+\cdots +b_M{z}^{-M}=1+\sum _{r=1}^M{b}_r{z}^{-r}$$
  则系统的脉冲响应$h(n)$为有限长，该系统的差分方程为
$$y(n)=f(n)+\sum _{r=1}^M{b}_{r}f(n-r)$$
  上式表明，FIR滤波器的输出$y(n)$只取决于输入值$f(n)$，$f(n-1)$……，与其他移位的输出无关。且他仅在$z=0$处有极点，故FIR滤波器是稳定的

二、无限脉冲响应数字滤波器

  按照不同的分类方法，数字滤波器有许多种类，但总起来可以分成两大类： 经典滤波器和现代滤波器。经典滤波器的特点是其输入信号中有用的频率成分和希望滤除的频率成分各占有不同的频带，通过一个合适的选频滤波器滤除干扰，得到纯净信号，达到滤波的目的。

  经典数字滤波器从滤波特性上分类，可以分成低通、高通、带通和带阻等滤波器。它们的理想幅频特性如图所示

  这种理想滤波器是不可能实现的，因为它们的单位脉冲响应均是非因果且无限长的，我们只能按照某些准则设计滤波器，使之在误差容限内逼近理想滤波器，理想滤波器可作为逼近的标准。

  无限长单位脉冲响应(IIR)滤波器和有限长单位脉冲响应(FIR)滤波器设计完全不同。IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是： 先设计过渡模拟滤波器得到系统函数$H_a(s)$，然后将$H_a(s)$按某种方法转换成数字滤波器的系统函数$H(z)$。

  模拟滤波器的理论和设计方法已发展得相当成熟，且有多种典型的模拟滤波器供我们选择，这些典型的滤波器各有特点：

• 巴特沃斯(Butterworth)滤波器具有单调下降的幅频特性
• 切比雪夫(Chebyshev)滤波器的幅频特性在通带或者阻带有等波纹特性，可以提高选择性
• 贝塞尔(Bessel)滤波器通带内有较好的线性相位特性
• 椭圆(Ellipse)滤波器的选择性相对前三种是最好的，但通带和阻带内均呈现等波纹幅频特性，相位特性的非线性也稍严重。

  这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用，而且所设计的系统函数都满足电路实现条件。由于模拟滤波器通常更适用硬件，所以现实IIR滤波器完全可以用硬件来实现，这里就不多介绍，重点是FIR滤波器的原理和实现

三、有限脉冲响应数字滤波器

  对于线性相位滤波器，经常采用FIR滤波器。FIR滤波器的设计方法和IIR滤波器的设计方法有很大差别。FIR滤波器设计任务是选择有限长度的h(n)，使频率响应函数H(ejω)满足技术指标要求。

3.1 线性相位FIR数字滤波器

对于长度为N的h(n)，频率响应函数为

式中，$H_g(\omega )$称为幅度特性; θ(ω)称为相位特性。注意，这里$H_g(\omega )$不同于$|H(e_{j\omega })|$，$H_g(\omega )$为ω的实函数，可能取负值，而$|H(e_{j\omega })|$总是正值。线性相位FIR滤波器是指θ(ω)是ω的线性函数，即

如果θ(ω)满足下式：

严格地说，此时θ(ω)不具有线性相位特性，但以上两种情况都满足群时延是一个常数，即

也称这种情况为线性相位。一般称满足（7.1.3）式是第一类线性相位；满足（7.1.4）式为第二类线性相位。${\theta }_0=-\pi /2$是第二类线性相位特性常用的情况，所以本章仅介绍这种情况。

3.1.1 线性相位FIR的时域约束条件

线性相位FIR滤波器的时域约束条件是指满足线性相位时，对h(n)的约束条件。

1. 第一类线性相位对h(n)的约束条件

第一类线性相位FIR数字滤波器的相位函数θ(ω)=－ωτ，由式（7.1.1）和（7.1.2）得到:

由式（7.1.5）得到:

将（7.1.6）式中两式相除得到：

即

移项并用三角公式化简得到:

函数h(n)sinω(n－τ)关于求和区间的中心(N－1)/2奇对称，是满足（7.1.7）式的一组解。因为sinω(n－τ)关于n=τ奇对称，如果取τ=(N－1)/2，则要求h(n)关于(N－1)/2偶对称，所以要求τ和h(n)满足如下条件:

由以上推导结论可知，如果要求单位脉冲响应为h(n)、长度为N的FIR数字滤波器具有第一类线性相位特性（严格线性相位特性），则h(n)应当关于n=(N－1)/2点偶对称。当N确定时，FIR数字滤波器的相位特性是一个确知的线性函数，即θ(ω)=－ω(N－1)/2。N为奇数和偶数时, h(n)的对称情况分别如表7.1.1中的情况1和情况2所示。

2. 第二类线性相位对h(n)的约束条件

第二类线性相位FIR数字滤波器的相位函数θ(ω)=－π/2－ωτ，由式（7.1.1）和（7.1.2），有：

经过同样的推导过程可得到

函数h(n)cos［ω(n－τ)］关于求和区间的中心(N－1)/2奇对称，是满足式（7.1.9）的一组解，因为cos［ω(n－τ)］关于n=τ偶对称，所以要求τ和h(n)满足如下条件：

由以上推导结论可知，如果要求单位脉冲响应为h(n)、长度为N的FIR数字滤波器具有第二类线性相位特性，则h(n)应当关于n=(N－1)/2点奇对称。N为奇数和偶数时h(n)的对称情况分别如表7.1.1中情况3和情况4所示。

3.1.2 线性相位FIR滤波器幅度特性Hg(ω)的特点

实质上，幅度特性的特点就是线性相位FIR滤波器的频域约束条件。将时域约束条件h(n)=±h(N－n－1)代入式（7.1.1），设h(n)为实序列，即可推导出线性相位条件对FIR数字滤波器的幅度特性$H_g(\omega )$的约束条件。当N取奇数和偶数时对$H_g(\omega )$的约束不同，因此，对于两类线性相位特性，下面分四种情况讨论其幅度特性的特点。这些特点对正确设计线性相位FIR数字滤波器具有重要的指导作用。为了推导方便，引入两个参数符号：

式中，[(N-1)/2]表示取不大于(N－1)/2的最大整数。显然，仅当N为奇数时，M=τ＝(N－1)/2

• 情况1： h(n)=h(N－n－1), N为奇数。

将时域约束条件h(n)=h(N－n－1)和θ(ω)=－ωτ代入式（7.1.1）和（7.1.2），得到:

所以

因为cos［ω(n-τ)］关于ω=0, π, 2π三点偶对称，所以由式（7.1.11）可以看出，$H_g(\omega )$关于ω=0, π, 2π三点偶对称。因此情况1可以实现各种（低通、高通、带通、带阻）滤波器。对于N=13的低通情况，$H_g(\omega )$的一种例图如表7.1.1中情况1所示。

• 情况2： h(n)=h(N－n－1), N为偶数

仿照情况1的推导方法得到:

式中，$\tau =(N-1)/2=N/2-1/2$。因为Ｎ是偶数，所以当ω=π时

而且cos［ω(n－τ)］关于过零点奇对称，关于ω=0和2π偶对称。所以$H_g(\omega )=0$，$H_g(\omega )$关于ω=π奇对称，关于ω=0和2π偶对称。因此，情况2不能实现高通和带阻滤波器。对N=12 的低通情况，$H_g(\omega )$如表7.1.1中情况2所示。

• 情况3： h(n)=－h(N－n－1)，N为奇数。
  　　将时域约束条件h(n)=－h(N－n－1)和θ(ω)=－π/2－ωτ代入式（7.1.1）和（7.1.2）， 并考虑
  　　

得到:

式中，N是奇数，τ=(N－1)/2是整数。所以，当ω=0，π, 2π时，sin［ω(n－τ)］=0, 而且sin［ω(n－τ)］关于过零点奇对称。因此$H_g(\omega )$关于ω=0, π, 2π三点奇对称。由此可见，情况3只能实现带通滤波器。对N=13的带通滤波器举例，Hg(ω)如表7.1.1中情况3所示。

• 情况4： h(n)=－h(N－n－1), N为偶数。
  　　用情况3的推导过程可以得到:
  
  式中，N是偶数，τ=(N－1)/2=N/2－1/2。所以，当ω=0, 2π时，sin［ω(n－τ)］=0；当ω=π时，sin［ω(n－τ)］=(－1)n－N/2, 为峰值点。而且sin［ω(n－τ)］关于过零点ω=0和2π两点奇对称，关于峰值点ω=π偶对称。因此$H_g(\omega )$关于ω=0和2π两点奇对称，关于ω=π偶对称。由此可见，情况4不能实现低通和带阻滤波器。对N=12的高通滤波器举例，Hg(ω)如表7.1.1中情况4所示。

3.1.3 线性相位FIR数字滤波器的零点分布特点

将h(n)=±h(N－1－n)代入上式, 得到:

由（7.1.14）式可以看出，如z=zi是H(z)的零点，其倒数$Z_i^{-1}$也必然是其零点；又因为h(n)是实序列，H(z)的零点必定共轭成对，因此$z_i^{\ast }$和$(z_i^{-1}{)}^{\ast }$也是其零点。这样，线性相位FIR滤波器零点必定是互为倒数的共轭对，确定其中一个，另外三个零点也就确定了, 如图7.1.1中$z_3、z_3^{-1}、z_3^{\ast }和(z_3^{\ast }{)}^{-1}$。当然，也有一些特殊情况，如图7.1.1中$z_1$、$z_2$和$z_4$情况。

3.2 窗函数法

3.2.1 窗函数法设计原理

设希望逼近的滤波器频率响应函数为$H_d(e^{j\omega })$，其单位脉冲响应是$h_d(n)$。

如果能够由已知的$H_d(e^{j\omega })$求出$h_d(n)$，经过Z变换可得到滤波器的系统函数。但通常以理想滤波器作为$H_d(e^{j\omega })$，其幅度特性逐段恒定，在边界频率处有不连续点，因而$h_d(n)$是无限时宽的，且是非因果序列。例如，线性相位理想低通滤波器的频率响应函数$H_d(e^{j\omega })$为

其单位脉冲响应$h_d(n)$为

由上式看到，理想低通滤波器的单位脉冲响应$h_d(n)$是无限长，且是非因果序列。$h_d(n)$的波形如图7.2.1（a）所示。为了构造一个长度为N的第一类线性相位FIR滤波器，只有将$h_d(n)$截取一段，并保证截取的一段关于n=（N－1）/2偶对称。设截取的一段用h(n)表示，即

式中, $R_N(n)$是一个矩形序列，长度为N，波形如图7.2.1（b）所示。由该图可知，当取值为（N－1）/2时，截取的一段h(n) 关于n=（N－1）/2偶对称，保证所设计的滤波器具有线性相位。　我们实际设计的滤波器的单位脉冲响应为h(n)，长度为N，其系统函数为H(z)，

这样用一个有限长的序列h(n)去代替$h_d(n)$，肯定会引起误差，表现在频域就是通常所说的吉布斯（Gibbs）效应。该效应引起过渡带加宽以及通带和阻带内的波动，尤其使阻带的衰减小，从而满足不了技术上的要求，如图7.2.2所示。这种吉布斯效应是由于将$h_d(n)$直接截断引起的，因此，也称为截断效应。下面讨论这种截断效应的产生，以及如何构造窗函数w(n), 用来减少截断效应，设计一个能满足技术要求的FIR线性相位滤波器。窗函数设计法的时域波形（矩形窗，N=30）:

吉普斯效应

以上就是用窗函数法设计FIR滤波器的思想。另外，我们知道$H_d(e^{j\omega })$是一个以2π为周期的函数，可以展为傅里叶级数，即

傅里叶级数的系数为$h_d(n)$，当然就是$H_d(e^{j\omega })$对应的单位脉冲响应。设计FIR滤波器就是根据要求找到N个傅里叶级数系数h(n)，n=1, 2, … , N－1，以N项傅氏级数去近似代替无限项傅氏级数，这样在一些频率不连续点附近会引起较大误差，这种误差就是前面说的截断效应，如图7.2.2所示。

因此，从这一角度来说，窗函数法也称为傅氏级数法。显然，选取傅氏级数的项数愈多，引起的误差就愈小，但项数增多即h(n)长度增加，也使成本和滤波计算量加大，应在满足技术要求的条件下，尽量减小h(n)的长度。

在（7.2.3）式中，$R_N(n)$（矩形序列）就是起对无限长序列的截断作用，可以形象地把$R_N(n)$看做一个窗口，h(n)则是从窗口看到的一段$h_d(n)$序列，所以称$h(n)=h_d(n)R_N(n)$为用矩形窗对$h_d(n)$进行加窗处理。下面分析用矩形窗截断的影响和改进的措施。为了叙述方便，用w(n)表示窗函数，用下标表示窗函数类型，矩形窗记为$w_R(n)$。用N表示窗函数长度。

根据傅里叶变换的时域卷积定理，得到（7.2.3）式的傅里叶变换：

式中，$H_d(e^{j\omega })$和$W_R(e^{j\omega })$分别是$h_d(n)$和$R_N(n)$的傅里叶变换，即

式中

$W_{Rg}(\omega )$称为矩形窗的幅度函数，如图7.2.3（b）所示，将图中［－2π/N, 2π/N］区间上的一段波形称为WRg(ω)的主瓣，其余较小的波动称为旁瓣。将$H_d(ej\omega )$写成$H_d(e^{j\omega })=H_{dg}(\omega )e^{-j\omega }$， 则按照（7.2.1）式，理想低通滤波器的幅度特性函数（如图7.2.3（a）所示）为

矩形窗加窗效应

将$H_d(e^{j\omega })$和$W_R(e^{j\omega })$代入（7.2.4）式，得到：

将$H(e^{j\omega })$写成$H(e^{j\omega })=H_g(\omega )e^{-j\omega }$ ，则

式中$H_g(\omega )$是$H(e^{j\omega })$的幅度特性。该式说明加窗后的滤波器的幅度特性等于理想低通滤波器的幅度特性$H_{dg}(\omega )$与矩形窗幅度特性$W_{Rg}(\omega )$的卷积。图7.2.3（f）表示$H_{dg}(\omega )$与$W_{Rg}(\omega )$卷积形成的$H_g(\omega )$波形。当ω=0时，$H_g(0)$等于图7.2.3（a）与（b）两波形乘积的积分，相当于对$W_{Rg}(\omega )$在$\pm {\omega }_c$之间一段波形的积分，当${\omega }_c>>2\pi /N$时，近似为±π之间波形的积分。 将H(0)值归一化到1。当$\omega ={\omega }_c$时，情况如图7.2.3（c）所示，当${\omega }_c>>2\pi /N$时，积分近似为$W_{Rg}(\theta )$一半波形的积分，对$H_g(0)$归一化后的值近似为1/2。

当$\omega ={\omega }_c－2\pi /N$时，情况如图7.2.2（d）所示，$W_R(\omega )$主瓣完全在区间$[-{\omega }_c,{\omega }_c]$之内，而最大的一个负旁瓣移到区间$[-{\omega }_c,{\omega }_c]$之外，因此$H_g(\omega c－2\pi /N)$有一个最大的正峰。当$\omega ={\omega }_c+2\pi /N$时，情况如图7.2.2（e）所示，$W_{Rg}(\omega )$主瓣完全移到积分区间外边，由于最大的一个负旁瓣完全在区间$[-{\omega }_c,{\omega }_c]$内，因此$H_g({\omega }_c+2\pi /N)$形成最大的负峰。图7.2.2表明，$H_g(\omega )$最大的正峰与最大的负峰对应的频率相距4π/N。通过以上分析可知，对hd(n)加矩形窗处理后，$H_g(\omega )$与原理想低通$H_{dg}(\omega )$的差别有以下两点：

（1） 在理想特性不连续点ω=ωc附近形成过渡带。过渡带的宽度近似等于WRg(ω)主瓣宽度4π/N。
（2） 通带内产生了波纹，最大的峰值在ωc－2π/N处。阻带内产生了余振，最大的负峰在ωc+2π/N处。通带与阻带中波纹的情况与窗函数的幅度谱有关， WRg(ω)旁瓣幅度的大小直接影响Hg(ω)波纹幅度的大小。

以上两点就是对hd(n)用矩形窗截断后，在频域的反映，称为吉布斯效应。这种效应直接影响滤波器的性能。通带内的波纹影响滤波器通带的平稳性，阻带内的波纹影响阻带内的衰减，可能使最小衰减不满足技术指标要求。当然，一般滤波器都要求过渡带愈窄愈好。下面研究如何减少吉布斯效应的影响，设计一个满足要求的FIR滤波器。

直观上，好像增加矩形窗的长度，即加大N，就可以减少吉布斯效应的影响。只要分析一下N加大时WRg(ω)的变化，就可以看到这一结论不是完全正确。我们讨论在主瓣附近的情况。在主瓣附近，按照式（7.2.5），WRg(ω)可近似为

该函数的性质是随x加大（N加大），主瓣幅度加高，同时旁瓣也加高，保持主瓣和旁瓣幅度相对值不变；另一方面，N加大时，WRg(ω)的主瓣和旁瓣宽度变窄，波动的频率加快。三种不同长度的矩形窗函数的幅度特性WRg(ω)曲线如图(a)、(b)、©所示。

用这三种窗函数设计的FIR滤波器的幅度特性Hg(ω)曲线如图(d)、(e)、(f)所示。因此，当N加大时，Hg(ω)的波动幅度没有多大改善，带内最大肩峰比H(0)高8.95%，阻带最大负峰值为H(0)的8.95%，使阻带最小衰减只有21 dB。加大N只能使Hg(ω)过渡带变窄（过渡带近似为主瓣宽度4π/N）。因此加大N, 并不是减小吉布斯效应的有效方法。

以上分析说明，调整窗口长度N只能有效地控制过渡带的宽度，而要减少带内波动以及增大阻带衰减，只能从窗函数的形状上找解决问题的方法。构造新的窗函数形状，使其谱函数的主瓣包含更多的能量，相应旁瓣幅度更小。旁瓣的减小可使通带、阻带波动减小，从而加大阻带衰减。但这样总是以加宽过渡带为代价的。

3.2.2 典型窗函数介绍

本节主要介绍几种常用窗函数的时域表达式、时域波形、幅度特性函数（衰减用dB计量）曲线，以及用各种窗函数设计的FIR数字滤波器的单位脉冲响应和损耗函数曲线。为了叙述简单，我们把这组波形图简称为“四种波形”。下面均以低通为例，Hd(ejω)取理想低通，ωc=π/2，窗函数长度N=31。

1.矩形窗（Rectangle Window）

前面已分析过，按照（7.2.5）式，其幅度函数为

为了描述方便，定义窗函数的几个参数:

旁瓣峰值αn——窗函数的幅频函数|Wg(ω)|的最大旁瓣的最大值相对主瓣最大值的衰减值（dB）；
　　过渡带宽度Bg——用该窗函数设计的FIR数字滤波器（FIRDF）的过渡带宽度；
　　阻带最小衰减αs——用该窗函数设计的FIRDF的阻带最小衰减。

2.三角形窗（Bartlett Window）

其频谱函数为

其幅度函数为

三角窗的四种波形如图7.2.5所示，参数为: αn=－25 dB; Bg=8π/N; αs=－25 dB。

3.汉宁（Hanning）窗——升余弦窗

当N>>1时， N－1≈N

汉宁窗的幅度函数WHng(ω)由三部分相加，旁瓣互相对消，使能量更集中在主瓣中。汉宁窗的四种波形如图7.2.6所示，参数为: α n=－31 dB; Bg=8π/N; α s=－44 dB。

4.哈明（Hamming）窗——改进的升余弦窗

其频谱函数WHm(ejω)为

其幅度函数WHmg(ω)为

当N>>１时，其可近似表示为

这种改进的升余弦窗，能量更加集中在主瓣中，主瓣的能量约占99.96%，瓣峰值幅度为40 dB，但其主瓣宽度和汉宁窗的相同，仍为8π/N。可见哈明窗是一种高效窗函数，所以MATLAB窗函数设计函数的默认窗函数就是哈明窗。哈明窗的四种波形如图7.2.7所示，参数为: α n=－41 dB; Bg=8π/N; α s=－53 dB。

5.布莱克曼（Blackman）窗

其频谱函数为

其幅度函数为

这样其幅度函数由五部分组成，它们都是移位不同，且幅度也不同的WRg(ω)函数，使旁瓣再进一步抵消。旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍。布莱克曼窗的四种波形如图7.2.8所示，参数为: α n=－57 dB; αΔB=12π/N; α s=－74 dB。

6.凯塞—贝塞尔窗（Kaiser-Basel Window）

凯塞—贝塞尔窗是一种参数可调的窗函数，是一种最优窗函数。

式中

I0(β)是零阶第一类修正贝塞尔函数，可用下面级数计算：

一般I0(β)取15~25项，便可以满足精度要求。 α参数可以控制窗的形状。一般加大，主瓣加宽，旁瓣幅度减小，典型数据为4< α <9。当α =5.44时，窗函数接近哈明窗。α =7.865时，窗函数接近布莱克曼窗。在设计指标给定时，可以调整α值，使滤波器阶数最低，所以其性能最优。凯塞（Kaiser）给出的估算β和滤波器阶数N的公式如下:

式中，Bt=|ωs－ωp|, 是数字滤波器过渡带宽度。应当注意，因为式（7.2.17）为阶数估算，所以必须对设计结果进行检验。另外，凯塞窗函数没有独立控制通带波纹幅度，实际中通带波纹幅度近似等于阻带波纹幅度。凯塞窗的幅度函数为

对的8种典型值，将凯塞窗函数的性能列于表7.2.1中，供设计者参考。由表可见, 当α =5.568时, 各项指标都好于哈明窗。6种典型窗函数基本参数归纳在表7.2.2中，可供设计时参考。

凯塞窗参数对滤波器的性能影

3.2.3 用窗函数法设计FIR滤波器的步骤

6种窗函数的基本参数

用窗函数法设计FIR滤波器的步骤如下：
　　(1) 根据对过渡带及阻带衰减的指标要求，选择窗函数的类型，并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣窄的窗函数。然后根据过渡带宽度估计窗口长度N。待求滤波器的过渡带宽度Bt近似等于窗函数主瓣宽度，且近似与窗口长度N成反比，N≈A/Bt，A取决于窗口类型，例如，矩形窗的A=4π，哈明窗的A=8π等，参数A的近似和精确取值参考表。
　　(2) 构造希望逼近的频率响应函数Hd(ejω)，即
　　
对所谓的“标准窗函数法”，就是选择Hd(ejω)为线性相位理想滤波器（理想低通、理想高通、理想带通、理想带阻）。以低通滤波器为例，Hdg(ω)应满足：

由于理想滤波器的截止频率ωc近似位于最终设计的FIRDF的过渡带的中心频率点，幅度函数衰减一半（约－6 dB）。所以如果设计指标给定通带边界频率和阻带边界频率ωp和ωs， 一般取

(3) 计算hd(n)。 如果给出待求滤波器的频响函数为Hd(ejω)，那么单位脉冲响应用下式求出：

如果Hd(ejω)较复杂，或者不能用封闭公式表示，则不能用上式求出hd(n)。我们可以对Hd(ejω)从ω=0到ω=2π采样M点，采样值为$H_{dM}(k)=H_d(e^{j\frac{2\pi }{M}k})$ ，k=0 ,1, 2, …, M－1，进行M点IDFT(IFFT)，得到：

根据频域采样理论，hdM(n)与hd(n)应满足如下关系：

因此，如果M选得较大，可以保证在窗口内hdM(n)有效逼近hd(n)。对（7.2.19）式给出的线性相位理想低通滤波器作为Hd(ejω)，由（7.2.2）式求出单位脉冲响应hd(n)：

为保证线性相位特性，$\alpha =(N－1)/2$
(4) 加窗得到设计结果：h(n)=hd(n)w(n)。

3.3 频率采样法

3.3.1 基本思想

设希望逼近的滤波器的频响函数用$H_d(e^{j\omega })$表示，对它在ω=0到2π之间等间隔采样N点，得到$H_d(k)$：

再对$H_d(k)$进行N点IDFT，得到h(n)：

式中，h(n)作为所设计的FIR滤波器的单位脉冲响应，其系统函数H(z)为

另外根据频率域采样理论，利用频率域采样值恢复原信号Z变换的（3.3.5）和（3.3.6）式，得到H(z)的内插表示形式：

此式就是直接利用频率采样值Hd(k)形成滤波器的系统函数，（7.3.3）式适合FIR直接型网络结构，（7.3.4）式适合频率采样结构。下面讨论两个问题：一个是为了设计线性相位FIR滤波器，频域采样序列$H_d(k)$应满足的条件；另一个是逼近误差问题及其改进措施。

3.3.2 约束条件

FIR滤波器具有线性相位的条件是h(n)为实序列，且满足h(n)=h(N－n－1)，在此基础上我们已推导出其频响函数应满足的条件是：

在ω=0~2π区间上N个等间隔的采样频点为

将$\omega ={\omega }_k$代入式中，并写成k的函数：

（7.3.9）~（7.3.12）式就是对频率采样值的约束条件。（7.3.11）式说明，N等于奇数时Hg(k)关于N/2点偶对称。（7.3.12）式说明，N等于偶数时，Hg(k)关于N/2点奇对称，且Hg(N/2)=0。

设用理想低通作为希望逼近的滤波器Hd(ejω)，截止频率为ωc，采样点数为N，Hg(k)和θ(k)用下列公式计算：

N=奇数时，

N=偶数时，

上面公式中的kc是通带内最后一个采样点的序号，所以kc的值取不大于［ωcN/(2π)］的最大整数。另外，对于高通和带阻滤波器，这里N只能取奇数。

3.3.3 改进措施

如果待逼近的滤波器为Hd(ejω)，对应的单位脉冲响应为hd(n)，则由频率域采样定理知道，在频域0~2π范围等间隔采样N点，利用IDFT得到的h(n)应是hd(n)以N为周期的周期延拓的主值区序列，即

如果Hd(ejω)有间断点，那么相应的单位脉冲响应hd(n)应是无限长的。这样，由于时域混叠及截断，使h(n)与hd(n)有偏差。所以，频域的采样点数N愈大，时域混叠愈小，设计出的滤波器频响特性愈逼近Hd(ejω)。

上面是从时域分析其设计误差的来源，下面从频域分析。频域采样定理表明，频率域等间隔采样H(k)，经过IDFT得到h(n)，由（3.3.7）和（3.3.8）式，得到H(ejω)=FT［h(n)］的内插表示形式：

式中

上式表明，在采样频点，ωk=2πk/N, k=0，1，2，…，N－1, Φ(ω－2πk/N)=1，因此采样点处　　　　与H(k)相等，逼近误差为0。在采样点之间，$H(e^{j{\omega }_k})$由N项H(k)Φ(ω－2πk/N)之和形成。频域幅度采样序列Hg(k)及其内插波形Hg(ω)如图7.3.1所示。

图7.3.1的上图 7.3.1(a) 中, 实线表示希望逼近的理想幅度函数Hdg(ω), 黑点表示幅度采样序列Hg(k); 下图7.3.1(b)中, 实线Hg(ω）与虚线Ηdg(ω)的误差与Hdg(ω)特性的平滑程度有关，Hdg(ω)特性愈平滑的区域，误差愈小；特性曲线间断点处，误差最大。表现形式为间断点变成倾斜下降的过渡带曲线，过渡带宽度近似为2π/N。通带和阻带内产生震荡波纹，且间断点附近振荡幅度最大，使阻带衰减减小，往往不能满足技术要求。

当然，增加N可以使过渡带变窄，但是通带最大衰减和阻带最小衰减随N的增大并无明显改善。且N太大，会增加滤波器的阶数，即增加了运算量和成本。N=15和N=75两种情况下的幅度内插波形Hg(ω)如图7.3.2所示，图中的空心圆和实心圆点分别表示N=15和N=75时的频域幅度采样。

在窗函数设计法中，通过加大过渡带宽度换取阻带衰减的增加。频率采样法同样满足这一规律。提高阻带衰减的具体方法是在频响间断点附近区间内插一个或几个过渡采样点，使不连续点变成缓慢过渡带，这样，虽然加大了过渡带，但阻带中相邻内插函数的旁瓣正负对消，明显增大了阻带衰减。

过渡带采样点的个数与阻带最小衰减${\alpha }_s$的关系以及使阻带最小衰减${\alpha }_s$最大化的每个过渡带采样值求解都要用优化算法解决。其基本思想是将过渡带采样值设为自由量，用一种优化算法（如线性规划算法）改变它们，最终使阻带最小衰减${\alpha }_s$最大。

3.3.4 频率采样法设计步骤

（1） 根据阻带最小衰减${\alpha }_s$选择过渡带采样点的个数m。
（2） 确定过渡带宽度Bt，估算频域采样点数（即滤波器长度）N。如果增加m个过渡带采样点，则过渡带宽度近似变成(m+1)2π/N。当N确定时，m越大，过渡带越宽。如果给定过渡带宽度Bt，则要求(m+1)2π/N≤Bt ，滤波器长度N必须满足如下估算公式:

（3） 构造一个希望逼近的频率响应函数:

设计标准型片断常数特性的FIR数字滤波器时，一般构造幅度特性函数Hdg(ω)为相应的理想频响特性，且满足表7.1.1要求的对称性。

（4） 按照（7.3.1）式进行频域采样，并加入过渡带采样。过渡带采样值可以设置为经验值，或用累试法确定，也可以采用优化算法估算。
　
（5） 对H(k)进行N点IDFT，得到第一类线性相位FIR数字滤波器的单位脉冲响应：

（6） 检验设计结果。如果阻带最小衰减未达到指标要求，则要改变过渡带采样值，直到满足指标要求为止。如果滤波器边界频率未达到指标要求，则要微调Hdg(ω)的边界频率。

3.4 等波纹最佳逼近法

用$H_d(\omega )$表示希望逼近的幅度特性函数，要求设计线性相位FIR数字滤波器时，$H_d(\omega )$必须满足线性相位约束条件。用$H_g(\omega )$表示实际设计的滤波器幅度特性函数。定义加权误差函数E(ω)为

式中，W(ω)称为误差加权函数，用来控制不同频段（一般指通带和阻带）的逼近精度。等波纹最佳逼近基于切比雪夫逼近，在通带和阻带以|E(ω)|的最大值最小化为准则，采用Remez多重交换迭代算法求解滤波器系数h(n)。所以W(ω)取值越大的频段, 逼近精度越高，开始设计时应根据逼近精度要求确定W(ω)，在Remez多重交换迭代过程中W(ω)是确知函数。

等波纹最佳逼近设计中，把数字频段分为“逼近（或研究）区域”和“无关区域”。逼近区域一般指通带和阻带，而无关区域一般指过渡带。设计过程中只考虑对逼近区域的最佳逼近。应当注意，无关区宽度不能为零，即$H_d(\omega )$不能是理想滤波特性。图7.4.1给出了等波纹滤波器技术指标的两种描述参数。

图7.4.1 (a)用损耗函数描述，即${\omega }_p=\pi /2$，${\alpha }_p=2dB$，${\omega }_s=11\pi /20$， $\alpha +s=20dB$。这是工程实际中常用的指标描述方法。但是，用等波纹最佳逼近设计法求滤波器阶数N和误差加权函数W(ω)时，要求给出滤波器通带和阻带的振荡波纹幅度${\delta }_1$和${\delta }_2$。图7.4.1(b)给出了用通带和阻带的振荡波纹幅度δ1和δ2描述的技术指标。显然，两种描述参数之间可以换算。如果设计指标以${\alpha }_p$和${\alpha }_s$给出，为了调用MATLAB工具箱函数remezord和remez进行设计，就必须由${\alpha }_p$和${\alpha }_s$换算出通带和阻带的振荡波纹幅度${\delta }_1$和${\delta }_2$。对比图7.4.2(a)和(b)得出关系式:

由式（7.4.2）和（7.4.3）得到

按照式（7.4.4）和（7.4.5）计算得到图7.4.1(b)中的参数: ${\delta }_1=0.1146$ ，${\delta }_2=0.1$。实际中, ${\delta }_1$和${\delta }_2$一般很小，这里为了观察等波纹特性及参数${\delta }_1$和${\delta }_2$的含义，特意取较大值。

下面举例说明误差加权函数W(ω)的作用，以及滤波器阶数N和波纹幅度${\delta }_1$和${\delta }_2$的制约关系。设期望逼近的通带和阻带分别为［0, π/4］和［5π/16, π］，对下面四种不同的控制参数，等波纹最佳逼近的损耗函数曲线分别如图7.4.2(a)、(b)、©和(d)所示。图中，W=［w1, w2］表示第一个逼近区［0, π/4］上的误差加权函数$W(\omega )=w_1$，第二个逼近区［5π/16, π］上的误差加权函数$W(\omega )=w_2$。图7.4.2(a)中, 通带频段［０，π/４］上的W(ω)=1，阻带频段［5π/16, π］上的W(ω)=10。

比较图7.4.2(a)、(b)、©和(d)可以得出结论: 当N一定时，误差加权函数W(ω)较大的频带逼近精度较高，W(ω)较小的频带逼近精度较低，如果改变W(ω)使通（阻）带逼近精度提高，则必然使阻（通）带逼近精度降低。滤波器阶数N增大才能使通带和阻带逼近精度同时提高。所以，W(ω)和N由滤波器设计指标（即${\alpha }_p$和${\alpha }_s$以及过渡带宽度）确定。所以用等波纹最佳逼近法设计FIR数字滤波器的过程是:
（1） 根据给定的逼近指标估算滤波器阶数N和误差加权函数W(ω)；
　（2） 采用remez算法得到滤波器单位脉冲响应h(n)。

3.5 IIR和FIR数字滤波器的比较

首先，从性能上来说，IIR滤波器系统函数的极点可位于单位圆内的任何地方，因此零点和极点相结合，可用较低的阶数获得较高的选择性，所用的存储单元少，计算量小，所以经济高效。但是这个高效率是以相位的非线性为代价的。相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器系统函数的极点固定在原点，因而只能用较高的阶数达到高的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数一般比IIR滤波器高5~10倍，使成本较高，信号延时也较大；如果按相同的选择性和相同的线性相位要求来说，则IIR滤波器就必须加全通网络进行相位校正，同样要大大增加滤波器的阶数和复杂性。

从结构上看，IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。另外，在这种结构中，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引起寄生振荡。相反，FIR滤波器主要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，运算误差引起的输出信号噪声功率也较小。此外，FIR滤波器可以采用FFT算法实现，在相同阶数的条件下，运算速度可以大大提高。

从设计工具看，IIR滤波器可以借助成熟模拟滤波器设计成果，因此一般都有封闭形式的设计公式可供准确计算，计算工作量比较小，对计算工具的要求不高。FIR滤波器计算通带和阻带衰减等仍无显式表达式，其边界频率也不易精确控制。一般，FIR滤波器的设计只有计算程序可循，因此对计算工具要求较高。但在计算机普及的今天，很容易实现其设计计算。

另外，也应看到，IIR滤波器虽然设计简单，但主要是用于设计具有片断常数特性的选频型滤波器，如低通、高通、带通及带阻等，往往脱离不了几种典型模拟滤波器的频响特性的约束。而FIR滤波器则要灵活得多，易于适应某些特殊的应用，如构成微分器或积分器，或用于巴特沃斯、切比雪夫等逼近不可能达到预定指标的情况，例如由于某些原因要求三角形振幅响应或一些更复杂的幅频响应形状，因而FIR滤波器有更大的适应性和更广阔的应用场合。

从上面的简单比较可以看到，IIR与FIR滤波器各有所长，所以在实际应用时应该全面考虑加以选择。例如，从使用要求上看，在对相位要求不敏感的场合，如语音通讯等，选用IIR滤波器较为合适，这样可以充分发挥其经济高效的特点；而对于图像信号处理，数据传输等以波形携带信号的系统，则对线性相位要求较高，采用FIR滤波器较好。