北京邮电大学《数字信号处理》课件-第7章 有限脉冲响应数字滤波器的设计

第7章

有限脉冲响应数字滤波器的设计北京邮电大学《数字信号处理》2本章主要内容利用窗函数法设计FIR滤波器

利用频域采样法设计FIR滤波器

线性相位FIR滤波器的条件和特点IIR与FIR数字滤波器的比较

FIR数字滤波器的Matlab仿真实现

3假设N为FIR数字滤波器的单位脉冲响应h(n)的长度，那么其系统函数H(z)为：IIR数字滤波器设计过程中只考虑了幅频特性，没有考虑相位特性，所设计的滤波器相位特性一般是非线性的。为了得到线性相位特性，则要采用全通网络进行相位校正，从而使得设计复杂，成本变高，也难以得到严格的线性相位。

引言4引言FIR数字滤波器很容易得到严格的线性相位。FIR数字滤波器的单位脉冲响应是有限长的，因此总是稳定的。FIR滤波器的设计方法：窗函数法频率取样法等纹波逼近法7.1线性相位FIR数字滤波器的条件及特点1.线性相位FIR数字滤波器单位脉冲响应h(n)的长度为N，则其频率响应函数为其中，为幅度特性，为相位特性。而线性相位FIR滤波器是指是w的线性函数。其中：第一类线性相位的相位特性函数是ω的严格线性函数:第二类线性相位FIRDF的相位特性函数如下式:

5

67.1线性相位FIR数字滤波器2.线性相位FIR滤波器的时域约束条件1）第一类线性相位对h(n)的约束条件第一类线性相位FIRDF的相位函数

，由(7.1.1)和（7.1.2）得

由式（7.1.5）得到

7.1线性相位FIR数字滤波器将式（7.1.6）中两式相除得：

即移项并用三角公式化简得到要求和h(n)满足如下条件:

77.1线性相位FIR数字滤波器2)第二类线性相位对h(n)的约束条件第二类线性相位FIRDF的相位函数由式（7.1.1）和式（7.1.2），有：经过同样的推导过程可得到：要求和h(n)满足如下条件

8

7.1线性相位FIR数字滤波器3.线性相位FIR滤波器幅度特性Hg(ω)的特点对幅度特性Hg(ω)的约束情况1:h(n)=h(N-n-1),N为奇数时Hg(ω)的特点：Hg(ω)关于0，π，2π三点偶对称，可以实现各种（低通、高通、带通、带阻）滤波器。情况2:h(n)=h(N-n-1),N为偶数时Hg(ω)的特点:Hg(ω)关于0，2π偶对称,而关于π奇对称，因此不能实现高通和带阻滤波器。97.1线性相位FIR数字滤波器情况3:h(n)=-h(N-n-1),N为奇数时Hg(ω)的特点：Hg(ω)关于0，π，2π三点奇对称，只能实现带通滤波器。情况4:h(n)=-h(N-n-1),N为偶数时Hg(ω)的特点：Hg(ω)关于0，2π奇对称,而关于π偶对称，不能实现低通和带阻滤波器。10线性相位FIRDF的时域和频域特性总结11表7.1.1(a)线性相位FIRDF的时域和频域特性一览表线性相位FIRDF的时域和频域特性总结12表7.1.1(b)线性相位FIRDF的时域和频域特性一览表线性相位FIRDF的零点分布特点4.线性相位FIRDF的零点分布特点将代入上式，得到：由(7.1.14)式可以看出：如果是的零点，其倒数也必然是其零点；又因为h(n)是实序列，H(z)的零点必定共轭成对，因此和也是其零点。这样线性相位FIR滤波器的零点确定其中一个，另外三个也确定了，当然也有特殊情况（共轭为本身or倒数为本身的情况，如图所示）1314设计思想寻找一个FIR滤波器，使其频率响应H(ejω)逼近理想FIR滤波器的频率响应Hd(ejω).其单位脉冲响应用hd(n)表示。一般情况下，Hd(ejω)在边界频率处有不连续点，因此hd(n)是无限长的，且是非因果的序列。7.2利用窗函数法设计FIR滤波器157.2利用窗函数法设计FIR滤波器理想低通滤波器的频率响应

其对应的单位脉冲响应可以看到hd(n)是一个无限长，非因果的序列。167.2利用窗函数法设计FIR滤波器将hd(n)截取长度为N的一段，构成h(n)为了保证设计的滤波器具有第一类线性相位，必须满足截取的一段关于n=(N-1)/2偶对称可以将h(n)看作是hd(n)与矩形窗相乘

17窗函数设计法的频域解释时域加窗频域卷积18窗函数时域表示频域表示幅度相位19FIR滤波器的幅频特性FIR滤波器的幅频特性20理想低通与矩形窗频谱函数卷积过程

21加窗对Hd(ω)的影响(1)在理想特性不连续点ωc附近形成过渡带。过滤带的宽度近似等于WR(θ)主瓣宽度，Δω=4π/N

。(2)通带内增加了波动，最大的峰值在ωc-2π/N

处。阻带内产生了余振，最大的负峰在ωc+2π/N处。通带与阻带中波动的情况与窗函数的幅度谱有关。WR(θ)波动愈快（加大时），通带与阻带内波动愈快，WR(θ)旁瓣的大小直接影响波动的大小。这些影响是对hd(n)加矩形窗引起的，称之为吉布斯效应。

22减小吉布斯效应的方法

增加矩形窗口的宽度N不能减少吉布斯效应的影响。

N的改变只能改变ω坐标的比例和的绝对大小，不能改变主瓣和旁瓣幅度相对值。加大N并不是减少吉布斯效应的有效方法。寻找合适的窗函数形状，使其谱函数的主瓣包含更多的能量，相应旁瓣幅度就变小了；旁瓣的减少可使通带与阻带波动减少，从而加大阻带的衰减。但这样总是以加宽过渡带为代价的。23几种常用的窗函数

矩形窗三角（Bartlett）窗汉宁（Hanning）窗哈明（Hamming）窗布莱克曼（Blackman）窗凯塞-贝塞尔（Kaiser-Basel-）窗24几种常用的窗函数

矩形窗图7.2.3矩形窗的四种波形25三角（Bartlett）窗三角（Bartlett）窗其频谱函数为：其幅度函数为：26汉宁（Hanning）窗汉宁（Hanning）窗27哈明（Hamming）窗

哈明（Hamming）窗28布莱克曼（Blackman）窗布莱克曼（Blackman）窗29凯塞-贝塞尔（Kaiser-Basel）窗凯塞-贝塞尔（Kaiser-Basel）窗30常用窗函数的波形31常用窗函数的频谱

6种窗函数的基本参数窗函数类型旁瓣峰值

n（dB）过渡带宽度

B阻带最小衰减

s（dB）近似值精确值矩形窗-134

/N1.8

/N-21三角窗-258

/N6.1

/N-25汉宁窗-318

/N6.2

/N-44哈明窗-418

/N6.6

/N-53布莱克曼窗-5712

/N11

/N-74凯塞窗(

=7.865)-57

10

/N-8032表中过渡带宽和阻带最小衰减是用对应的窗函数设计的FIR数字滤波器的频率相响应指标。除了以上6种窗函数外，比较有名的还有Chebyshev窗、Gaussian窗。337.2.3用窗函数法设计FIR滤波器的步骤

1.根据允许的过渡带宽及阻带衰减，选定窗函数和N值。2.给出希望设计的滤波器的频率响应函数Hd(ejω)。3.计算hd(n)

如果Hd(ejω)不能用简单函数表示，可以用求和代替积分。34用窗函数法设计FIR滤波器的步骤

4.将hd(n)与窗函数相乘得FIR数字滤波器的冲激响应h(n)5.计算FIR数字滤波器的频率响应，并验证是否达到所要求的指标35例窗函数法设计FIR滤波器例7.2.1用窗函数法设计线性相位高通FIRDF，要求通带截止频率ωp=

/2rad，阻带截止频率ωs=

/4rad，通带最大衰减

p=1dB，阻带最小衰减

s=40dB。解：（1）选择窗函数w(n)，计算窗函数长度N。已知阻带最小衰减

s=40dB，由表(7.2.2)可知汉宁窗和哈明窗均满足要求，这里选择汉宁窗。本例中过渡带宽度

，汉宁窗的精确过渡带宽度，解不等式得到：N≥24.8，而对高通滤波器N要取奇数，所以此处N取25.所以有：36例窗函数法设计FIR滤波器（2）构造式中，（3）求出将τ=12代入，得：公式中第一项对应全通滤波器，第二项是低通，两者之差即为高通。37例窗函数法设计FIR滤波器（4）加窗窗函数的MATLAB设计函数简介实际设计的时候一般用Matlab工具箱函数，可调用工具箱函数fir1实现窗函数法设计步骤的（2）~（4）的过程。1）fir1是用窗函数法设计线性相位FIRDF的工具箱函数，以实现线性相位FIRDF的标准窗函数法设计。2）用fir2函数设计时，可以指定任意形状的，所以称之为任意形状幅度特性窗函数法设计函数，它实质是一种频率采样法与窗函数法的综合设计函数38397.3利用频域采样法设计FIR滤波器1.基本思想频域采样法是在频率域对理想滤波器Hd(ejω)采样，在采样点上设计的滤波器H(ejω)和理想滤波器Hd(ejω)幅度值相等，然后根据频率域的采样值求得实际设计的滤波器的频率特性H(ejω)。对理想滤波器的频率特性Hd(ejω)在[0,2

]范围内等间隔地取样N个点

407.3利用频域采样法设计FIR滤波器根据插值公式此式就是直接利用频率采样值Hd(k)形成滤波器的系统函数，适合频率采样结构412.设计线性相位滤波器时对Hd(k)的约束条件当h(n)是偶对称，h(n)=h(N-n-1)，N为奇数。H(ω)是偶对称的对H(ejω)在[0,2

]范围内等间隔地取样N个点422.设计线性相位滤波器时对Hd(k)的约束条件当h(n)是偶对称，h(n)=h(N-n-1)，N为偶数。H(ω)是偶对称的对H(ejω)在[0,2

]范围内等间隔地取样N个点433.逼近误差及其改进措施采用频域取样法设计的FIR数字滤波器在阻带内的衰减很小，在实际应用中往往达不到要求。产生这种现象的原因是由于在通带边缘采样点的陡然变化而引起的起伏振荡。增加阻带衰减的方法是在通带和阻带的边界处增加一些过渡的采样点，从而减小频带边缘的突变，也就减小了起伏振荡，增大了阻带最小衰减。

3.逼近误差及其改进措施的内插表示形式：式中因此采样点处H(ejwk)与H(k)相等，逼近误差为044453.逼近误差及其改进措施

46图7.3.1频率采样法设计过程中的波形（N=15）图7.3.2采样点数N不同时的逼近误差比较（N=15,75）3.逼近误差及其改进措施m123

s44～54dB65～75dB85～95dB47将过渡带采样点的个数m与滤波器阻带最小衰减

s的经验数据列于表7.3.1中，我们可以根据给定的阻带最小衰减

s选择过渡带采样点的个数m。表7.3.1过渡带采样点的个数m与滤波器阻带最小衰减

s的经验数据484.频率采样法设计步骤(1)根据阻带最小衰减选择过渡带采样点的个数m。(2)确定过渡带宽度，估算频域采样点数（即滤波器长度）N。如果增加m个过渡带采样点，则过渡带宽度近似变成。当N确定时，m越大，过渡带越宽。如果给定过渡带宽度，则要求，即滤波器长度N必须满足以下等式：(3)构造一个希望逼近的频率响应函数：设计标准型片段常数特性的FIR数字滤波器时，一般构造幅度特性函数为相应的理想频响特性，且满足表7.1.1要求的对称性。494.频率采样法设计步骤(4)按照(7.3.1)式进行频域采样：

并加入过渡采样。过渡带采样值可以设置为经验值，或用累试法确定，也可以采用优化方法估算。(5)对进行N点IDFT，得到第一类线性相位FIR数字滤波器的单位脉冲响应：(6)检验设计结果。如果阻带最小衰减未达到指标要求，则要改变过渡带采样值，直到满足指标要求为止。如果滤波器边界频率未达到指标要求，则要微调的边界频率。7.4利用等波纹最佳逼近法设计FIR数字滤波器等波纹最佳逼近法的特点：（1）是一种优化设计法，它克服了窗函数设计法和频率采样法的缺点，使最大误差（即波纹的峰值）最小化，并在整个逼近频段上均匀分布。即幅频响应在通带和阻带都是等波纹的；（2）可以分别控制通带和阻带波纹幅度。（3）在滤波器长度给定的条件下，使加权误差波纹幅度最小化。（4）与窗函数设计法和频率采样法比较，设计的滤波器性能价格比最高。阶数相同时，这种设计法使滤波器的最大逼近误差最小，即通带最大衰减最小，阻带最小衰减最大；指标相同时，这种设计法使滤波器阶数最低。比较例题7.2.3和例题7.2.1，设计指标要求相同的带阻滤波器，等波纹最佳逼近法需要28阶，而窗函数法需要79阶！507.4.1等波纹最佳逼近法的基本思想

Hd(ω)——表示希望逼近的幅度特性函数。Hg(ω)——表示实际设计的滤波器幅度特性函数。

——加权误差函数51W(ω)——称为误差加权函数，用来控制不同频段（一般指通带和阻带）的逼近精度。等波纹最佳逼近基于切比雪夫逼近，在通带和阻带以的最大值最小化为准则，采用Remez多重交换迭代算法求解滤波器系数h(n)。

W(ω)取值越大的频段逼近精度越高，开始设计时应根据逼近精度要求确定W(ω)，在Remez多重交换迭代过程中W(ω)是确知函数。由于切比雪夫（Chebyshev）和雷米兹（Remez）对解决该问题做出了贡献，所以又称之为切比雪夫逼近法，或雷米兹逼近法。7.4.1等波纹最佳逼近法的基本思想等波纹最佳逼近设计中，把数字频段分为“逼近（或研究）区域”——一般指通带和阻带“无关区域”——一般指过渡带设计过程中只考虑对逼近区域的最佳逼近。注意：无关区宽度不能为零，即Hd(ω)不能是理想滤波特性。否则，设计必将失败。527.4.1等波纹最佳逼近法的基本思想等波纹滤波器的技术指标及其描述参数：等波纹最佳逼近设计法的数学证明复杂，所以本节略去其复杂的数学推导，只介绍其基本思想和实现线性相位FIR数字滤波器的等波纹最佳逼近设计的MATLAB信号处理工具箱函数remez和remezord。53[M,fo,mo,w]=remezord(f,m,rip,Fs)hn=remez(M,fo,mo,w)求滤波器阶数N和误差加权函数W(ω)时，要求给出

1和2工程实际中设计指标通常以p和s给出两种描述参数之间的换算关系：7.4.1等波纹最佳逼近法的基本思想54由图可见：和N由滤波器设计指标（即

p和s，以及过渡带宽度）确定

期望逼近的通带为[0,

/4]，阻带为[5

/16,]

，过渡带已知等波纹最佳逼近法设计FIR数字滤波器的过程:

（1）根据给定的逼近指标估算滤波器阶数N和误差加权函数；（2）采用remez算法得到滤波器单位脉冲响应h(n)。

MATLAB工具箱函数remezord和remez就是完成以上2个设计步骤的有效函数。557.4.2remez和remezord函数介绍1．Remez：采用remez算法实现线性相位FIR数字滤波器的等波纹最佳逼近设计。其调用格式为:hn=remez(M,f,m,w)返回单位脉冲响应向量hn。remez函数的调用参数（M,f,m,w）一般通过调用remezord函数来计算。调用参数含义如下:

M为FIR数字滤波器阶数，hn长度N=M+1。

f和m给出希望逼近的幅度特性：

f为边界频率向量，0≤f≤1，要求f为单调增向量（即f(k)<f(k+1)），而且从0开始，以1结束，1对应数字频率ω=

（模拟频率Fs/2，Fs表示时域采样频率）。

m是与f对应的幅度向量，m与f长度相等，m(k)表示频点f(k)的幅度响应值。如果用命令Plot(f,m)画出幅频响应曲线，则k为奇数时，频段[f(k)，f(k+1)]上的幅频响应就是期望逼近的幅频响应值，频段[f(k+1)，f(k+2)]为无关区。简言之，Plot(f,m)命令画出的幅频响应曲线中，起始频段为第一段，奇数频段为逼近区，偶数频段为无关区。567.4.2remez和remezord函数介绍例如，通带频段，阻带频段，则，f=[0,1/4,5/16,1]；m=[1,1,0,0]plot(f,m)画出的幅度特性曲线如图所示。图中奇数段（第1、3段）的水平幅度为希望逼近的幅度特性，偶数段（第2段）的下降斜线为无关部分，逼近时形成过渡带，并不考虑该频段的幅频响应形状。

w为误差加权向量，其长度为f的一半。w(i)表示对m中第i个逼近频段的误差加权值。图7.4.2(a)中w=[1,10]。缺省w时，默认w为全1（即每个逼近频段的误差加权值相同）。57希望逼近的幅度特性曲线remez函数还可以设计两种特殊滤波器：希尔伯特变换器和数字微分器。调用格式：hn=remez(M,f,m,w,'hilbert')hn=remez(M,f,m,w,'defferentiator')7.4.2remez和remezord函数介绍2．Remezord：根据逼近指标估算FIR数字滤波器的最低阶数M、误差加权向量w和归一化边界频率向量f。返回参数作为remez函数的调用参数。调用格式:

[M,fo,mo,w]=remezord(f,m,rip,Fs)58hn=remez(M,fo,mo,w)参数说明：

f与remez中类似，这里f可以是模拟频率（Hz）或归一化数字频率，但必须以0开始，以Fs/2（用归一化频率时对应1）结束，而且省略了0和Fs/2两个频点。Fs为采样频率，缺省时默认Fs=2Hz。注意：这里f的长度（包括省略的0和Fs/2两个频点）是m的两倍，即m中的每个元素表示f给定的一个逼近频段上希望逼近的幅度值。例如，对上面希望逼近的幅度特性曲线图，f=[1/4,5/16];m=[1,0]。7.4.2remez和remezord函数介绍注意:①省略Fs时，f中必须为归一化频率。②有时估算的阶数M略小，使设计结果达不到指标要求，这时要取M+1或M+2（必须注意对滤波器长度N=M+1的奇偶性要求）。所以必须检验设计结果。③如果无关区（过渡带）太窄，或截止频率太接近零频率和Fs/2时，设计结果可能不正确。

rip表示f和m描述的各逼近频段允许的波纹幅度（幅频响应最大逼近误差），f的长度是rip的两倍。例：对上图所示的低通滤波器，

59对高通滤波器，remezord函数的调用参数f,m,rip和Fs

调用remez和remezord函数设计线性相位FIR数字滤波器的关键问题：设计指标remezord函数的调用参数f,m,rip和Fs其中，Fs一般是根据实际信号处理要求（按照采样定理）确定。下面给出由各种滤波器设计指标确定remezord调用参数f,m和rip的公式，编程时直接套用即可。601、低通滤波器设计指标逼近通带:，通带最大衰减:

逼近阻带:，阻带最小衰减:

remezord调用参数:remezord函数的调用参数f,m,rip和Fs2、高通滤波器设计指标 逼近通带:，通带最大衰减:；逼近阻带:，阻带最小衰减:。调用参数:

61，，

3、带通滤波器设计指标 ，通带最大衰减:

逼近阻带:

，阻带最小衰减:

remezord调用参数:逼近通带:

4、带阻滤波器设计指标 ，阻带最大衰减:

逼近通带:

，通带最小衰减:remezord调用参数:，，

逼近阻带:注意：工程实际中常常给出对模拟信号的滤波指标要求，设计数字滤波器，对输入模拟信号采样后进行数字滤波。这时，调用参数f可以用模拟频率表示。但是，调用remezord时一定要加入采样频率参数Fs。

等波纹最佳逼近法设计举例例7.4.1

利用等波纹最佳逼近法重新设计FIR带阻通滤波器。指标与例7.2.3相同，即:

逼近通带:，通带最大衰减:

逼近阻带:，阻带最小衰减:解:调用remezord和remez函数求解。本例设计程序为ep741.m。％ep741.m:例7.4.1用remez函数设计带阻滤波器f=[0.2,0.35,0.65,0.8]; %省略了0和1m=[1,0,1];rp=1;rs=60;%由式（7.4.4）和式（7.4.5）求通带和阻带波纹幅度dat1、dat2和rip：dat1=(10^(rp/20)-1)/(10^(rp/20)+1);dat2=10^(-rs/20);rip=[dat1,dat2,dat1];[M,fo,mo,w]=remezord(f,m,rip);hn=remez(M,fo,mo,w);％以下绘图检验部分略去62等波纹最佳逼近法设计举例程序运行结果:M=28。即h(n)的长度N=29。注意：例7.2.3中N=80。63等波纹最佳逼近法设计举例例7.4.2利用等波纹最佳逼近法设计FIR数字低通滤波器，实现对模拟信号的数字滤波处理。要求与例7.2.2相同。求出h(n)，并画出损耗函数曲线和相频特性曲线。解:将例7.2.2中所给指标重写如下:通带截止频率: