有限脉冲响应数字滤波器设计课件_第1页
有限脉冲响应数字滤波器设计课件_第2页
有限脉冲响应数字滤波器设计课件_第3页
有限脉冲响应数字滤波器设计课件_第4页
有限脉冲响应数字滤波器设计课件_第5页
已阅读5页,还剩256页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

学习目标1.掌握线性相位FIR数字滤波器的特点2.掌握窗函数设计法3.了解频率抽样设计法4.了解FIR滤波器与IIR滤波器的比较第1页/共261页IIR滤波器的优缺点IIR数字滤波器的优点:

可以利用模拟滤波器设计的结果,而模拟滤波器的设计有大量图表可查,方便简单。IIR数字滤波器的缺点:

相位的非线性,将引起频率的色散,若须线性相位,则要采用全通网络进行相位校正,使滤波器设计变得复杂,成本也高。第2页/共261页FIR滤波器优点

*H(z)永远稳定:设FIR滤波器单位冲激响应h(n)长度为N,其系统函数H(z)为:

H(z)是z-1的N-1次多项式,它在z平面上有N-1个零点,原点z=0是N-1阶重极点。

因此,H(z)永远稳定。稳定和线性相位特性是FIR滤波器突出的优点*可以做到严格线性相位*可以具有任意的幅度特性第3页/共261页为何要设计FIR滤波器(1)语音处理,图象处理以及数据传输要求线性相位,任意幅度。(即要求信道具有线性相位特性)而FIR数字滤波器具有严格的线性相位,而且同时可以具有任意的幅度特性。(2)另外FIR数字滤波器的单位抽样响应是有限长的,因而滤波器一定是稳定的,只要经过一定的延时,任何非因果有限长序列都变成因果的有限序列。(3)FIR可以用FFT算法来实现过滤信号。第4页/共261页FIRDF设计任务、设计方法、设计思路设计思路:

根据设计指标,求出所选运算结构要求的h(n)或H(z):

线性卷积和快速卷积型结构,求FIRDF的h(n).

级联和频率采样型结构,求FIRDF的H(z).设计方法:*窗口设计法*频率采样设计法*切比雪夫等波纹逼近法X*计算机辅助设计法

X设计任务:选择有限长度的h(n)

,使传输特性H(ej)满足技术指标要求第5页/共261页7.1线性相位FIR数字滤波器的条件和特点

内容:

FIR滤波器具有线性相位的条件及幅度特性以及零点、网络结构的特点。第6页/共261页Hg(ω)称为幅度特性;()称为相位特性。1.线性相位滤波器

------线性相位概念(7.1.1)(7.1.2)系统的单位脉冲响应h(n)长度为N,系统传输函数为:第7页/共261页(1)幅度特性Hg(ω)

注意:Hg(ω)不同于|H(ejω)|,区别:

*Hg(ω)为ω的实函数,是标量函数,可以包括正值、负值和零;*|H(ejω)|总是正值,而且是ω的偶对称函数和周期函数。

两者在某些ω值上相位相差π第8页/共261页(A)H(ejω)线性相位是指()是ω的线性函数,即

()=-,

:群时延,为常数

(7.1.3)(B)如果θ(ω)满足下式:

()=-+0,0是起始相位

(7.1.4)

严格地说,此时θ(ω)不具有线性相位,但以上两种情况都满足群时延是一个常数,也称这种情况为线性相位。(2)相位特性()()=-

:第一类线性相位;A类()=-+0

:第二类线性相位;B类(groupdelay)为群延迟函数θ0=-π/2是第二类线性相位特性常用的情况,所以本章仅介绍这种情况第9页/共261页

线性相位FIR滤波器的时域约束条件是指满足线性相位时,对h(n)的约束条件。一般要求是h(n)实序列。2.线性相位FIR的时域约束条件第10页/共261页由式(7.1.5)得到:第一类线性相位FIR数字滤波器的相位函数

θ(ω)=-ωτ(7.1.6)(7.1.5)1)第一类线性相位对h(n)的约束条件第11页/共261页移项并用三角公式化简得到:

将(7.1.6)式中两式相除得到:

函数h(n)sinω(n-τ)关于求和区间的中心(N-1)/2奇对称,是满足(7.1.7)式的一组解。因为sinω(n-τ)关于n=τ奇对称,如果取τ=(N-1)/2,则要求h(n)关于(N-1)/2偶对称。即(7.1.7)第12页/共261页(7.1.8)

如果要求单位脉冲响应为h(n)、长度为N的FIR数字滤波器具有第一类线性相位特性(严格线性相位特性),则h(n)应当关于n=(N-1)/2点偶对称。

当N确定时,FIR数字滤波器的相位特性是一个确知的线性函数,即θ(ω)=-ω(N-1)/2。

N为奇数和偶数时,h(n)的对称情况分别如表7.1.1中的情况1和情况2所示。第13页/共261页表7.1.1线性相位FIR数字滤波器的时域和频域特性一览第一类线性相位第14页/共261页第二类线性相位FIR数字滤波器的相位函数

θ(ω)=-π/2-ωτ函数h(n)cos[ω(n-τ)]关于求和区间的中心(N-1)/2奇对称,是满足式(7.1.9)的一组解,因为cos[ω(n-τ)]关于n=τ偶对称,则要求h(n)关于(N-1)/2奇对称。(7.1.9)2)第二类线性相位对h(n)的约束条件第15页/共261页(7.1.10)

如果要求单位脉冲响应为h(n)、长度为N的FIR数字滤波器具有第二类线性相位特性,则h(n)应当关于n=(N-1)/2点奇对称。

N为奇数和偶数时h(n)的对称情况分别如表7.1.1中情况3和情况4所示。第16页/共261页表7.1.1线性相位FIR数字滤波器的时域和频域特性一览第二类线性相位第17页/共261页

设h(n)为实序列,将时域约束条件h(n)=±h(N-n-1)代入式(7.1.1),即可推导出线性相位条件对FIR数字滤波器的幅度特性Hg(ω)的约束条件。当N取奇数和偶数时对Hg(ω)的约束不同,因此,对于两类线性相位特性,分四种情况讨论其幅度特性的特点。这些特点对正确设计线性相位FIR数字滤波器具有重要的指导作用3.线性相位FIR滤波器幅度特性Hg(ω)的特性幅度特性的特点就是线性相位FIR滤波器的频域约束条件线性相位时域约束条件:h(n)=±h(N-n-1)第18页/共261页式中,[(N-1)/2]

表示取不大于(N-1)/2的最大整数。仅当N为奇数时,M=τ=(N-1)/2。为了推导方便,引入两个参数符号:

对于两类线性相位、h(n)长度N的奇偶分为四类:1)h(n)=h(N-n-1)--偶对称,N为奇数2)h(n)=h(N-n-1)--偶对称,N为偶数3)h(n)=-h(N-n-1)--奇对称,N为奇数4)h(n)=-h(N-n-1)--奇对称,N为偶数第19页/共261页线性相位滤波器时域特性图例1、N为奇数的偶对称:例如N=11,对称中心为n0123456789102、N为偶数的偶对称:例如N=10,对称中心为n01234567893、N为奇数的奇对称:例如N=11,对称中心为n0123456789104、N为偶数的奇对称:例如N=10,对称中心为n0123456789第20页/共261页情况1:h(n)=h(N-n-1),N为奇数h(n)=h(N-n-1)和θ(ω)=-ωτ代入式(7.1.1)和(7.1.2),得到:Hg()第21页/共261页

因为cos[ω(n-τ)]关于ω=0,π,2π三点偶对称,所以Hg(ω)关于ω=0,π,2π三点偶对称。

因此情况1可以实现各种(低通、高通、带通、带阻)滤波器。对于N=13的低通情况,Hg(ω)的一种例图如表7.1.1中情况1所示。(7.1.11)幅度特性:第22页/共261页仿照情况1的推导方法得到:(7.1.12)式中,情况2:h(n)=h(N-n-1),N为偶数幅度特性:第23页/共261页

而且cos[ω(n-τ)]关于过零点奇对称,关于ω=0和2π偶对称。所以Hg(ω)关于ω=π奇对称,关于ω=0和2π偶对称。情况2不能实现高通和带阻滤波器。

对N=12的低通情况,Hg(ω)如表7.1.1中情况2所示因为N是偶数,所以当=时第24页/共261页将时域约束条件h(n)=-h(N-n-1)和θ(ω)=-

π/2-

代入式(7.1.1)和(7.1.2),并考虑h[(N-1/2)]=0,得到:情况3:h(n)=-h(N-n-1),N为奇数第25页/共261页式中,N是奇数,τ=(N-1)/2是整数。

所以,当ω=0,π,2π时,sin[ω(n-τ)]=0,而且sin[ω(n-τ)]关于过零点奇对称。因此Hg(ω)关于ω=0,π,2π三点奇对称。由此可见,情况3只能实现带通滤波器。对N=13的带通滤波器举例,Hg(ω)如表7.1.1中情况3所示。幅度特性:第26页/共261页

式中,N是偶数,τ=(N-1)/2=N/2-1/2。所以,当ω=0,2π时,sin[ω(n-τ)]=0;当ω=π时,sin[ω(n-τ)]=(-1)n-N/2,为峰值点。而且sin[ω(n-τ)]关于过零点ω=0和2π两点奇对称,关于峰值点ω=π偶对称。因此Hg(ω)关于ω=0和2π两点奇对称,关于ω=π偶对称。

情况4不能实现低通和带阻滤波器。对N=12的高通滤波器举例,Hg(ω)如表7.1.1中情况4所示。情况4:h(n)=-h(N-n-1),N为偶数用情况3的推导过程可以得到:(7.1.13)幅度特性:第27页/共261页为了便于比较,将上面四种情况的h(n)及其幅度特性需要满足的条件列于表7.1.1中。应当注意,对每一种情况仅画出满足幅度特性要求的一种例图。例如,情况1仅以低通的幅度特性曲线为例。当然也可以画出满足情况1的幅度约束条件(Hg(ω)关于ω=0,π,2π三点偶对称)的高通、带通和带阻滤波器的幅度特性曲线。所以,仅从表7.1.1就认为情况1只能设计低通滤波器是错误的。第28页/共261页频域特性表Typeh(n)NI情况1N为奇数II情况2N为偶数III情况3N为奇数IV情况4N为偶数第29页/共261页相位特性辅助序列适用幅度特性及特点N为奇数N为偶数N为奇数N为偶数低通高通带通带阻低通带通带通高通带通四种线性相位FIR滤波器第30页/共261页四种线性相位滤波器第31页/共261页四种线性相位滤波器奇对称单位冲激响应h(n)=-h(N-n-1)第32页/共261页例1

如果系统的单位脉冲响应为0≤n≤4其他n求FIR数字滤波器传输函数解:这是第一种情况的线性相位FIR数字滤波器该系统的频率响应为因为h(n)的长度N=5,群延迟是整数,τ(ω)=(N-1)/2=2。振幅相位群延迟第33页/共261页MATLABw=-0*pi:0.1:2*pi;n=0:4;x=1.^n;X=(exp(j*2*w)).*(sin(2.5*w))./(sin(0.5*w));subplot(221)stem(n,x,'.');title('Sequencex(n)');xlabel('Timeindexn');subplot(222)plot(w/pi,abs(X));title('MagnitudeofDTFTofx(n)');axis([0,2,0,max(abs(X))]);subplot(224)plot(w/pi,angle(X)),;title('PhaseofDTFTofx(n)');xlabel('Frequency*piinratians');axis([0,2,min(angle(X)),max(angle(X))]);第34页/共261页例2

系统的单位脉冲响应为0≤n≤5其他n

求FIR数字滤波器传输函数振幅相位群延迟

解:h(n)为偶对称且长度N=6,是第二种情况的线性相位FIR数字滤波器该系统的频率响应为h(n)的长度N=6,群延迟:τ(ω)=(N-1)/2=2.5第35页/共261页例3

系统的单位脉冲响应为h(n)=δ(n)-δ(n-2),求FIR数字滤波器传输函数。振幅相位群延迟解:

h(n)为奇对称且长度N=3,这是第三种情况的线性相位FIR数字滤波器该系统的频率响应为h(n)的长度N=3,群延迟:τ(ω)=(N-1)/2=1第36页/共261页例4系统的单位脉冲响应为h(n)=δ(n)-δ(n-1)求FIR数字滤波器传输函数。振幅相位群延迟解:h(n)为奇对称且长度N=2,这是第四种情况的线性相位FIR数字滤波器该系统的频率响应为h(n)的长度N=2,群延迟:τ(ω)=(N-1)/2=0.5第37页/共261页将h(n)=±h(N-1-n)代入上式,得到:(7.1.14)4.线性相位FIR数字滤波器的零点分布特点

M=N-1-n第38页/共261页

如z=zi是H(z)的零点,其倒数1/zi也必然是其零点;又因为h(n)是实序列,H(z)的零点必定共轭成对,因此zi*和(zi-1)*也是其零点。

线性相位FIR滤波器零点必定是互为倒数的共轭对,确定其中一个,另外三个零点也就确定了,如图7.1.1中

。当然,也有一些特殊情况,如图7.1.1中z1、z2和z4情况。(7.1.14)第39页/共261页图7.1.1线性相位FIR数字滤波器的零点分布第40页/共261页

(1)zi既不在实轴上,也不在单位圆上,则零点是互为倒数的两组共轭对,如图(a)所示。(4个)(2)zi不在实轴上,但是在单位圆上,则共轭对的倒数是它们本身,故此时零点是一组共轭对,如图(b)所示。(2个)(3)zi在实轴上但不在单位圆上,只有倒数部分,无复共轭部分。故零点对如图(c)所示。(2个)(4)zi既在实轴上又在单位圆上,此时只有一个零点,有两种可能,或位于z=1,或位于z=-1,如图(d)、(e)所示。

(1个)互为倒数的共轭对的四种可能性:第41页/共261页线性相位FIR滤波器的零点位置图第42页/共261页线性相位滤波器的零点位置分布图第43页/共261页由幅度响应的讨论可知:*第二种情况的线性相位滤波器由于Hg(π)=0,因此必然有单根z=-1。*第四种情况的线性相位滤波器由于Hg(0)=0,因此必然有单根z=1。*而第三种情况的线性相位滤波器由于Hg(0)=Hg(π)=0,因此这两种单根z=±1都必须有。

了解了线性相位FIR滤波器的特点,便可根据实际需要选择合适类型的FIR滤波器,同时设计时需遵循有关的约束条件。讨论线性相位FIR滤波器的设计方法时,都要用到这些特点。

第44页/共261页

例一个FIR线性相位滤波器的单位脉冲响应是实数的,且n<0和n>6时h(n)=0。如果h(0)=1且系统函数在z=0.5ejπ/3和z=3各有一个零点,H(z)的表达式是什么?H1(z)=(1-0.5ejπ/3

z-1)(1-0.5e-jπ/3z-1)=1-0.5z-1+0.25z-2

解:因为n<0和n>6时h(n)=0,且h(n)是实值,所以当H(z)在z=0.5ejπ/3

有一个复零点时,则在它的共轭位置z=0.5e-jπ/3

处一定有另一个零点。这个零点共轭对产生如下的二阶因子:第45页/共261页

线性相位的约束条件需要在这两个零点的倒数位置上有零点,所以H(z)同样必须包括如下的有关因子:

系统函数还包含一个z=3的零点,同样线性相位的约束条件需要在z=1/3也有一个零点。于是,H(z)还具有如下因子:由此,我们有

最后,多项式中零阶项的系数为A,为使h(0)=1,必定有:A=1。第46页/共261页FIR第一类线性相位网络结构x(n-(N/2-1))x(n-(N-1)/2)第47页/共261页FIR第二类线性相位网络结构第48页/共261页

设计FIR数字滤波器最简单的方法是窗函数法。一般是先给定所要求的理想滤波器的频率响应Hd(ejω),要求设计一个FIR滤波器频率响应 ,去逼近理想的频率响应Hd(ejω)。7.2利用窗函数法设计FIR滤波器

7.2.1窗函数法设计原理hd

(n)一般情况下是非因果无限长序列,需对其进行截短和因果化处理。第49页/共261页改变窗形状改变窗长度误差?窗函数法设计思路:时域方法无限长有限长第50页/共261页

设希望设计的滤波器传输函数为Hd(ejω),hd(n)是与其对应的单位脉冲响应,因此一、窗函数法设计原理

1.构造希望逼近的频率响应函数Hd(ejω)第51页/共261页2.求出hd(n)(7.2.2)相应的单位取样响应hd(n)为:(7.2.1)设计的理想低通滤波器传输函数以设计低通滤波器为例:确定Hd(ejω)

a?a-延迟常数非因果无限长序列第52页/共261页w(n)——为窗函数,w(n)长度为N。

如果要求设计线性相位FIRDF,则要求h(n)关于(N-1)/2点偶对称,而hd(n)关于n=a点偶对称,所以要求a=(N-1)/2,且w(n)关于(N-1)/2点偶对称。各种常用窗函数都满足这种要求。

为了构造一个长度为N的线性相位滤波器,只有将hd(n)截取一段,并保证截取的一段对(N-1)/2对称。设截取的一段用h(n)表示,即

h(n)=hd(n)·w

(n)

(7.2.3)3.加窗得到FIRDF单位脉冲响应h(n)第53页/共261页

实际实现的滤波器的单位取样响应为h(n),长度为N:

系统函数为H(z):

传输函数:H(ejω):图7.2.1理想低通的单位脉冲响应及矩形窗w(n)=RN(n)aa第54页/共261页H(ejω)

与Hd(ejω)的误差:截断误差

傅立叶级数系数h(n)是H(ejω)对应的单位取样响应:设计FIR滤波器是根据要求找到有限项个傅立叶级数系数,以有限项傅立叶级数近似代替无限项傅立叶级数。

-----窗函数法也称为傅立叶级数法周期函数第55页/共261页

Hd(ejω)是一个以2π为周期的函数,可以展为傅氏级数,即根据复卷积定理,得到:(7.2.4)Hd(ejω)和RN(ejω)分别是hd(n)和RN(n)的傅里叶变换二、窗函数法的性能分析H(ejω)=FT[h(n)]=FT[hd(n)RN(n)]w(n)=RN(n)第56页/共261页(7.2.5)矩形窗的幅度函数RN(n)的傅里叶变换:RN(ejω)——条件矩形窗的相位函数第57页/共261页理想低通滤波器的幅度特性Hd(ω)为Hd(ejω)写成下式:hd(n)的傅里叶变换:Hd(ejω)幅度特性第58页/共261页

将H(ejω)写成下式:

(7.2.6)滤波器幅度特性H

(ω)等于理想低通滤波器Hd(ω)与窗函数的幅度特性函数RN(ω)的卷积滤波器幅度特性第59页/共261页1000幅度特性加窗后的影响图示窗函数的主瓣旁瓣旁瓣第60页/共261页

图7.2.2矩形窗对理想低通幅度特性的影响

对hd(n)加矩形窗处理后,H(ω)和原理想低通Hd(ω)差别有以下两点:

(1)在理想特性不连续点ω=ωc附近形成过渡带。过渡带的宽度,近似等于RN(ω)主瓣宽度,即4π/N。

(2)通带内增加了波动,最大的峰值在ωc-2π/N处。阻带内产生了余振,最大的负峰在ωc+2π/N处。在主瓣附近,RN(ω)可近似为第61页/共261页理想低通加窗后对幅度特性的影响:(1)理想幅度特性的陡直边沿被加宽,形成过渡带;过渡带:是由窗函数的主瓣宽度引起的。

N越大-窗函数的主瓣宽度越小-过渡带带宽越小(2)过渡带两侧附近产生起伏的肩峰和纹波;肩峰和纹波:是由窗函数的旁瓣引起的。

旁瓣相对值越大-纹波起伏就越大旁瓣相对值:完全取决与窗函数的形状,与N无关。(吉布斯效应)选用旁瓣幅度较小的窗函数结论:为了提高滤波器的性能,尽可能要求窗函数:i.主瓣宽度尽量窄-获得较陡的过渡带ii.旁瓣相对值尽量小-改善通带的平稳度和增大阻带的衰减

第62页/共261页吉布斯效应

所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象,且最大波纹大约为幅度的9%,这个现象称为Gibbs现象。

第63页/共261页7.2.2常用的窗函数

设h(n)=hd(n)w(n)式中w(n)表示窗函数。

本节主要介绍几种常用窗函数的时域表达式、时域波形、幅度特性函数(衰减用dB计量)曲线,以及用各种窗函数设计的FIR数字滤波器的单位脉冲响应和损耗函数曲线。为了叙述简单,我们把这组波形图简称为“四种波形”。下面均以低通为例,Hd(ejω)取理想低通,ωc=π/2,窗函数长度N=31。为了描述方便,定义窗函数的几个参数:旁瓣峰值n—窗函数的幅频函数|Wg(ω)|的最大旁瓣的最大值相对主瓣最大值的衰减值(dB);过渡带宽度Bg—用该窗函数设计的FIR数字滤波器(FIRDF)的过渡带宽度;阻带最小衰减s—用该窗函数设计的FIRDF的阻带最小衰减。第64页/共261页1.矩形窗(RectangleWindow)wR(n)=RN(n)频率响应为主瓣宽度为4π/N,第一旁瓣比主瓣低13dB幅度函数(7.2.7)第65页/共261页2.三角形窗(BartlettWindow)

(7.2.8)其频率响应为

(7.2.9)主瓣宽度为8π/N,比矩形窗主瓣宽度增加一倍,但旁瓣却小很多幅度函数为(7.2.10)三角窗的四种波形参数为:

n=-25dB;Bg=8π/N;s=-25dB第66页/共261页图7.2.5三角窗的四种波形三角窗的四种波形参数为:

n=-25dB;Bg=8π/N;s=-25dB第67页/共261页3.汉宁(Hanning)窗——升余弦窗当N1时,N-1≈N,这三部分之和,使旁瓣互相抵消,能量更集中在主瓣,它的最大旁瓣值比主瓣值约低31dB。但是代价是主瓣宽度比矩形窗的主瓣宽度增加一倍,即为8π/N。(7.2.11)第68页/共261页图7.2.6汉宁窗的四种波形参数为:

n=-31dB;Bg=8π/N;

s=-44dB。第69页/共261页汉宁窗的幅度特性第70页/共261页

4.哈明(Hamming)窗——改进的升余弦窗其频域函数WHm(ejω)为其幅度函数WHm(ω)为当N>>1时,可近似表示为

与汉宁窗相比,主瓣宽度相同,为8π/N,但旁瓣又被进一步压低,结果可将99.963%的能量集中在窗谱的主瓣内,它的最大旁瓣值比主瓣值约低41dB。(7.2.12)第71页/共261页图7.2.7哈明窗的四种波形改进的升余弦窗的能量更加集中在主瓣中,主瓣的能量约占99.96%,瓣峰值幅度为40dB,但其主瓣宽度和汉宁窗的相同,仍为8π/N。可见哈明窗是一种高效窗函数,所以MATLAB窗函数设计函数的默认窗函数就是哈明窗。波形参数:n=-41dB;Bg=8π/N;

s=-53dB第72页/共261页5.布莱克曼(Blackman)窗(7.2.13)频谱函数为幅度函数为(7.2.14)由5部分组成,使旁瓣互相抵消,能量更集中在主瓣,但代价是主瓣宽度比矩形窗的主瓣宽度增加3倍,即为12π/N。第73页/共261页图7.2.8布莱克曼窗的四种波形波形参数:

n=-57dB;Bg=12π/N;

s=-74dB。幅度函数由五部分组成,它们都是移位不同,且幅度也不同的WRg(ω)函数,使旁瓣再进一步抵消。旁瓣峰值幅度进一步增加,其幅度谱主瓣宽度是矩形窗的3倍。第74页/共261页6.凯塞—贝塞尔窗(Kaiser-BaselWindow)

以上五种窗函数都称为参数固定窗函数,每种窗函数的旁瓣幅度都是固定的。凯塞—贝塞尔窗是一种参数可调的窗函数,是一种最优窗函数。式中

I0(x)是零阶第一类修正贝塞尔函数,可用下面级数计算:(7.2.15)第75页/共261页*一般I0(x)取15~25项,便可以满足精度要求。*α参数可以控制窗的形状。一般α加大,主瓣加宽,旁瓣幅度减小,典型数据为4<α<9。当α=5.44时,窗函数接近哈明窗。α=7.865时,窗函数接近布莱克曼窗。在设计指标给定时,可以调整值,使滤波器阶数最低,所以其性能最优。凯塞(Kaiser)给出的估算β和滤波器阶数M(h(n)的长度N=M+1)的公式如下:(7.2.17)(7.2.16)第76页/共261页式中,Bt=|ωs-ωp|,是数字滤波器过渡带宽度。应当注意,因为式(7.2.17)为阶数估算,所以必须对设计结果进行检验。另外,凯塞窗函数没有独立控制通带波纹幅度,实际中通带波纹幅度近似等于阻带波纹幅度。凯塞窗的幅度函数为(7.2.18)对的8种典型值,将凯塞窗函数的性能列于表7.2.1中,供设计者参考。由表可见,当

=5.568时,各项指标都好于哈明窗。6种典型窗函数基本参数归纳在表7.2.2中,可供设计时参考。第77页/共261页表7.2.1凯塞窗参数对滤波器的性能影

第78页/共261页表7.2.26种窗函数的基本参数Bt=A/NA=4A=8表中过渡带宽和阻带最小衰减是用对应的窗函数设计的FIR数字滤波器的频率响应指标。随着数字信号处理的不断发展,学者们提出的窗函数已多达几十种,除了上述6种窗函数外,比较有名的还有Chebyshev窗、Gaussian窗[5,6]。第79页/共261页图7.2.4常用的窗函数第80页/共261页

图7.2.5常用窗函数的幅度特性(a)矩形窗;(b)巴特利特窗(三角形窗);(c)汉宁窗;(d)哈明窗;(e)布莱克曼窗第81页/共261页

图7.2.6理想低通加窗后的幅度特性(N=51,ωc=0.5π)(a)矩形窗;(b)巴特利特窗(三角形窗);(c)汉宁窗;

(d)哈明窗;(e)布莱克曼窗第82页/共261页wn=boxcar(N)

%列向量wn中返回长度为N的矩形窗函数w(n)wn=bartlett(N)

%列向量wn中返回长度为N的三角窗函数w(n)wn=hanning(N)

%列向量wn中返回长度为N的汉宁窗函数w(n)wn=hamming(N) %列向量wn中返回长度为N的哈明窗函数w(n)wn=blackman(N)

%列向量wn中返回长度为N的布莱克曼窗函数w(n)wn=kaiser(N,beta)

%列向量wn中返回长度为N的凯塞—贝塞尔窗函数w(n)

MATLAB信号处理工具箱提供了14种窗函数的产生函数,下面列出上述6种窗函数的产生函数及其调用格式:MATLAB第83页/共261页用窗函数法设计FIR滤波器的步骤如下:

(1)根据对过渡带及阻带衰减的指标要求,选择窗函数的类型,并估计窗口长度N。

先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下,尽量选择主瓣窄的窗函数。然后根据过渡带宽度估计窗口长度N。待求滤波器的过渡带宽度Bt近似等于窗函数主瓣宽度,且近似与窗口长度N成反比,N≈A/Bt,A取决于窗口类型例如:矩形窗的A=4π,哈明窗的A=8π等,参数A的近似和精确取值参考表7.2.2。7.2.3用窗函数法设计FIR滤波器的步骤第84页/共261页(2)构造希望逼近的频率响应函数Hd(ejω),即

对所谓的“标准窗函数法”,就是选择Hd(ejω)为线性相位理想滤波器(理想低通、理想高通、理想带通、理想带阻)。以低通滤波器为例,幅度特性Hdg(ω)应满足:(7.2.19)

理想滤波器的截止频率ωc近似位于最终设计的FIRDF的过渡带的中心频率点,幅度函数衰减一半(约-6dB)。如果设计指标给定通带边界频率和阻带边界频率ωp和ωs,3dB截止频率一般取:(7.2.20)第85页/共261页(3)计算hd(n)。

I

如果给出待求滤波器的频响函数为Hd(ejω),那么单位脉冲响应用下式求出:(7.2.21)II

如果Hd(ejω)较复杂,或者不能用封闭公式表示,则不能用上式求出hd(n)。对Hd(ejω)从ω=0到ω=2π采样M点,采样值为

,k=0,1,2,…,M-1,进行M点IDFT(IFFT),得到:(7.2.22)根据频域采样理论,hdM(n)与hd(n)应满足如下关系:如果M选得较大,可以保证在窗口内hdM(n)有效逼近hd(n)线性相位理想低通滤波器=(N-1)/2第86页/共261页(4)加窗得到设计结果:h(n)=hd(n)w(n)。第87页/共261页

【例7.2.1】用窗函数法设计线性相位高通FIRDF,要求通带截止频率ωp=π/2rad,阻带截止频率ωs=π/4rad,通带最大衰减

p=1dB,阻带最小衰减

s=40dB。第88页/共261页(2)构造Hd(ejω):式中解

(1)选择窗函数w(n),计算窗函数长度N。已知阻带最小衰减

s=40dB,由表(7.2.2)可知汉宁窗和哈明窗均满足要求,选择汉宁窗。本例中过渡带宽度Bt≤ωp-ωs=π/4,汉宁窗的精确过渡带宽度Bt=6.2π/N,所以要求Bt=6.2π/N≤π/4,解之得N≥24.8。对高通滤波器N必须取奇数,取N=25。由式(7.2.11),有第89页/共261页(3)求出hd(n):将τ=12代入得第一项δ(n-12)对应全通滤波器,第二项是截止频率为3π/8的理想低通滤波器的单位脉冲响应,二者之差就是理想高通滤波器的单位脉冲响应。这就是求理想高通滤波器的单位脉冲响应的另一个公式。第90页/共261页(4)加窗:第91页/共261页

实际设计时一般用MATLAB工具箱函数。可调用工具箱函数fir1实现窗函数法设计步骤(2)~(4)的解题过程。

(1)fir1是用窗函数法设计线性相位FIR数字滤波器的工具箱函数,以实现线性相位FIR数字滤波器的标准窗函数法设计。这里的所谓“标准”,是指在设计低通、高通、带通和带阻FIR滤波器时,Hd(ejω)分别表示相应的线性相位理想低通、高通、带通和带阻滤波器的频率响应函数。因而将所设计的滤波器的频率响应称为标准频率响应。7.2.4窗函数法的MATLAB设计函数简介第92页/共261页

Fir1的调用格式及功能:

·hn=fir1(M,wc),返回6dB截止频率为wc的M阶(单位脉冲响应h(n)长度N=M+1)FIR低通(wc为标量)滤波器系数向量hn,默认选用哈明窗。滤波器单位脉冲响应h(n)与向量hn的关系为

h(n)=hn(n+1)n=0,1,2,…,M而且满足线性相位条件:h(n)=h(N-1-n)。其中wc为对π归一化的数字频率,0≤wc≤1。当wc=[wcl,wcu]时,得到的是带通滤波器,其-6dB通带为wcl≤ω≤wcu。第93页/共261页

hn=fir1(M,wc,′ftype′),可设计高通和带阻FIR滤波器。当ftype=high时,设计高通FIR滤波器;当ftype=stop,且wc=[wcl,wcu]时,设计带阻FIR滤波器。应当注意,在设计高通和带阻FIR滤波器时,阶数M只能取偶数(h(n)长度N=M+1为奇数)。不过,当用户将M设置为奇数时,fir1会自动对M加1。

·hn=fir1(M,wc,window),可以指定窗函数向量window。如果缺省window参数,则fir1默认为哈明窗。例如:hn=fir1(M,wc,bartlett(M+1)),使用Bartlett窗设计;

hn=fir1(M,wc,blackman(M+1)),使用blackman窗设计;

hn=fir1(M,wc,'ftype',window),通过选择wc、ftype和window参数(含义同上),可以设计各种加窗滤波器。第94页/共261页(2)fir2为任意形状幅度特性的窗函数法设计函数,用fir2设计时,可以指定任意形状的Hd(ejω),它实质是一种频率采样法与窗函数法的综合设计函数。主要用于设计幅度特性形状特殊的滤波器(如数字微分器和多带滤波器等)。用help命令查阅其调用格式及调用参数的含义。第95页/共261页

例7.2.1

的设计程序ep721.m如下:%ep721.m:例7.2.1用窗函数法设计线性相位高通FIR数字滤波器

wp=pi/2;ws=pi/4;

Bt=wp-ws; %计算过渡带宽度

N0=ceil(6.2*pi/Bt);%根据表7.2.2汉宁窗计算所需h(n)长度N0,ceil(x)取大于等于x的最小整数

N=N0+mod(N0+1,2);%确保h(n)长度N是奇数

wc=(wp+ws)/2/pi;%计算理想高通滤波器通带截止频率(关于π归一化)

hn=fir1(N-1,wc,'high',hanning(N)); %调用fir1计算高通FIR数字滤波器的h(n)

第96页/共261页

%绘图部分M=1024;hk=fft(hn,M);n=0:N;subplot(2,2,1);stem(n,hn,'.');line([0,30],[0,0])xlabel('n');ylabel('h(n)');k=1:M/2;w=2*(0:M/2-1)/M;subplot(2,2,2);plot(w,20*log10(abs(hk(k))));axis([0,1,-80,5]);xlabel('w/pi');ylabel('20lg|Hg(w)|');gridon第97页/共261页运行程序得到h(n)的25个值:

h(n)=[-0.0004 -0.00060.0028

0.0071

-0.0000 -0.0185-0.0210

0.0165

0.06240.0355

0.1061

-0.2898 0.6249-0.2898

-0.1061 0.03550.0624

0.0165

-0.0210

0.0185-0.0000

0.0071 0.0028-0.0006-0.0004]图7.2.9高通FIR数字滤波器的h(n)波形及损耗函数曲线第98页/共261页

【例7.2.2】对模拟信号进行低通滤波处理,要求通带0≤f≤1.5kHz内衰减小于1dB,阻带2.5kHz≤f≤∞上衰减大于40dB。希望对模拟信号采样后用线性相位FIR数字滤波器实现上述滤波,采样频率Fs=10kHz。用窗函数法设计满足要求的FIR数字低通滤波器,求出h(n),并画出损耗函数曲线。为了降低运算量,希望滤波器阶数尽量低。第99页/共261页

阻带截止频率为阻带最小衰减为

s=40dB(2)用窗函数法设计FIR数字低通滤波器,为了降低阶数选择凯塞窗。根据式(7.2.16)计算凯塞窗的控制参数为解

(1)确定相应的数字滤波器指标:通带截止频率为第100页/共261页指标要求过渡带宽度Bt=ωs-ωp=0.2π,根据式(7.2.17)计算滤波器阶数为取满足要求的最小整数M=23。所以h(n)长度为N=M+1=24。但是,如果用汉宁窗,h(n)长度为N=40。理想低通滤波器的通带截止频率ωc=(ωs+ωp)/2=0.4π,所以由式(7.2.2)和式(7.2.3),得到:式中,w(n)是长度为24(

=3.395)的凯塞窗函数。第101页/共261页实现本例设计的MATLAB程序为ep722.m。%例7.2.2用凯塞窗函数设计线性相位低通FIR数字滤波器

fp=1500;fs=2500;rs=40;

wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;

Bt=ws-wp;%计算过渡带宽度

alph=0.5842*(rs-21)^0.4+0.07886*(rs-21);%根据(7.2.16)式计算kaiser窗的控制参数α

N=ceil((rs-8)/2.285/Bt);%根据(7.2.17)式计算kaiser窗所需阶数Nwc=(wp+ws)/2/pi;%计算理想高通滤波器通带截止频率(关于π归一化)

hn=fir1(N,wc,kaiser(N+1,alph));%调用kaiser计算低通FIRDF的h(n)

第102页/共261页%以下绘图部分M=1024;hk=fft(hn,M);n=0:N;subplot(2,2,1);stem(n,hn,'.');line([0,30],[0,0])xlabel('n');ylabel('h(n)');k=1:M/2;w=2*(0:M/2-1)/M;subplot(2,2,2);plot(w,20*log10(abs(hk(k))));axis([0,1,-80,5]);xlabel('w/pi');ylabel('20lg|Hg(w)|');gridon第103页/共261页运行程序得到h(n)的24个值:h(n)=[0.0039

0.0041-0.0062-0.01470.0000

0.0286

0.0242-0.0332-0.0755

0.0000

0.19660.3724

0.3724

0.1966-0.0000-0.0755-0.0332

0.02420.0286

0.0000-0.0147-0.0062

0.0041

0.0039]图7.2.10低通FIR数字滤波器的h(n)波形及损耗函数曲线第104页/共261页

【例7.2.3】窗函数法设计一个线性相位FIR带阻滤波器。要求通带下截止频率ωlp=0.2π,阻带下截止频率ωls=0.35π,阻通带上截止频率ωus=0.65π,通带上截止频率ωup=0.8π,通带最大衰减

p=1dB,阻带最小衰减

s=60dB。

第105页/共261页解之得N=80。调用参数解

本例直接调用fir1函数设计。因为阻带最小衰减

s=60dB,所以选择布莱克曼窗,再根据过渡带宽度选择滤波器长度N,布莱克曼窗的过渡带宽度Bt=12π/N,所以第106页/共261页设计程序为ep723.m,参数计算也由程序完成。%ep723.m:例7.2.3用窗函数法设计线性相位带阻FIR数字滤波器

wlp=0.2*pi;wls=0.35*pi;wus=0.65*pi;wup=0.8*pi;%设计指标参数赋值

B=wls-wlp;%过渡带宽度

N=ceil(12*pi/B);%计算阶数N,ceil(x)为大于等于x的最小整数wp=[(wls+wlp)/2/pi,(wus+wup)/2/pi];%设置理想带通截止频率

hn=fir1(N,wp,‘stop’,blackman(N+1));%带阻滤波器要求h(n)长度为奇数,所以取N+1%绘图部分第107页/共261页%绘图部分M=1024;hk=fft(hn,M);n=0:N;subplot(2,2,1);stem(n,hn,'.');xlabel('n');ylabel('h(n)');gridonk=1:M/2+1;w=2*(0:M/2)/M;subplot(2,2,2);plot(w,20*log10(abs(hk(k))));axis([0,1,-90,5]);xlabel('w/pi');ylabel('20lg|Hg(w)|');gridon第108页/共261页图7.2.11带阻FIR数字滤波器的h(n)波形及损耗函数曲线程序运行结果:

N=81由于h(n)数据量太大,因而仅给出h(n)的波形及损耗函数曲线,如图7.2.11所示。第109页/共261页

例用矩形窗、汉宁窗和布莱克曼窗设计FIR低通滤波器,设N=11,ωc=0.2πrad。N=11解:用理想低通作为逼近滤波器:用矩形窗设计:第110页/共261页用汉宁窗设计:用布莱克曼窗设计:

低通幅度特性第111页/共261页优点窗口法设计的主要优点是简单,使用方便。窗口函数大多有封闭的公式可循,性能、参数都已有表格、资料可供参考,计算程序简便,所以很实用。缺点是通带和阻带的截止频率不易控制。

窗口法设计的优点缺点第112页/共261页常用的理想滤波器1.理想低通滤波器频率响应:单位脉冲响应:第113页/共261页2.理想高通滤波器频率响应:单位脉冲响应:高通滤波器相当于一个全通滤波器减去一个低通滤波器第114页/共261页频率响应:单位脉冲响应:3.理想带通滤波器带通滤波器相当于两个低通滤波器相减,其中一个低通滤波器截止频率为c2,另一个截止频率为c1第115页/共261页频率响应:单位脉冲响应:4.理想带阻滤波器带阻滤波器相当于一个低通滤波器(截止频率为c1)加一个高通滤波器(截止频率为c2)第116页/共261页解

AFDF技术指标ωp=0.3πrad,ωs=0.46πrad,

as=45dB查表可知,海明窗和布莱克曼窗均可提供大于45dB的阻带衰减,但海明窗具有较小的过渡带从而具有较小的长度N。根据题意,所要设计的滤波器的过渡带为

由表可知,利用海明窗设计的滤波器的过渡带宽

Δω=8π/N例

根据下列技术指标,设计一个FIR低通滤波器。

通带截止频率p=30πrad/s,

阻带截止频率

s=46πrad/s,阻带衰减as=45dB。Ts=0.01s第117页/共261页3dB通带截止频率ωc理想低通滤波器的单位脉冲响应为海明窗为

低通滤波器单位脉冲响应的长度为取取N=51则所设计的滤波器的单位脉冲响应为第118页/共261页解

查表可知,海明窗和布莱克曼窗均可提供大于50dB的衰减。但海明窗具有较小的过渡带从而具有较小的长度N。根据题意,所要设计的滤波器的过渡带为

由表可知,利用海明窗设计的滤波器的过渡带宽

Δω=8π/N,例根据下列技术指标,设计一个FIR低通滤波器。

通带截止频率ωp=0.2π,通带允许波动ap=0.25dB;

阻带截止频率ωs=0.3π,阻带衰减as=50dB。第119页/共261页3dB通带截止频率为

理想低通滤波器的单位脉冲响应为海明窗为低通滤波器单位脉冲响应的长度为或取第120页/共261页则所设计的滤波器的单位脉冲响应为N=80所设计的滤波器的频率响应为

滤波器长N=80,实际阻带衰减为As=53dB,通带波动为Ap=0.0316dB,均满足设计要求。dB第121页/共261页低通滤波器设计结果理想低通滤波器的单位脉冲响应hd(n)海明窗函数实际低通滤波器的单位脉冲响应h(n)实际低通滤波器的幅频特性|H(ejω)第122页/共261页例用窗函数法设计线性相位高通FIRDF,技术指标如下:通带截止频率:ωp=/2rad阻带截止频率:ωs=/4rad通带最大衰减:αp=1dB阻带最小衰减:

αs=40dB第123页/共261页解(1)选择窗函数w(n)

已知阻带最小衰减αs=40dB,由表可知汉宁窗和海明窗均满足要求——选海明窗。本例中过渡带宽度Δ≤ωp-ωs=π/4,海明窗的过渡带宽度Δ=8π/N,所以要求

解之得N≥32 对高通滤波器,N必须取奇数,所以,N=33。第124页/共261页其中

第125页/共261页(3)求出hd(n):

δ(n-16)对应全通滤波器:τ=16所以第126页/共261页(4)加窗:第127页/共261页-ωc≤ω-ω0≤ωc

0≤ω<ω0-ωc,ω0+ωc<ω≤π

(1)设计N为奇数时的h(n)。(2)设计N为偶数时的h(n)。(3)若改用海明窗设计,求以上两种形式的h(n)表达式。

例用矩形窗设计一个线性相位带通滤波器第128页/共261页解

根据该线性相位带通滤波器的相位

可知该滤波器只能是h(n)=h(N-1-n)即h(n)偶对称的情况,h(n)偶对称时,可为第一类和第二类滤波器,其频响第129页/共261页

(1)当N为奇数时,h(n)=h(N-1-n),可知H(ejω)为第一类线性相位滤波器,H(ω)关于ω=0,π,2π有偶对称结构。题目中仅给出了Hd(ejω)在0~π上的取值,但用傅里叶反变换求hd(n)时,需要Hd(ejω)在一个周期[-π,π]或[0,2π]上的值,因此,Hd(ejω)需根据第一类线性相位滤波器的要求进行扩展,扩展结果为ω0-ωc≤ω≤ω0+ωc,-ω0-ωc≤ω≤-ω0+ωc-ω0+ωc<ω<ω0-ωc,-π≤ω<-ω0-ωc,ω0+ωc<ω≤π第130页/共261页则

h(n)=hd(n)RN(n)第131页/共261页

(2)N

为偶数时,H(ejω)为第二类线性相位滤波器,H(ω)关于ω=0呈偶对称。所以,Hd(ejω)在[-π,π]之间的扩展同上,则hd(n)也同上,即:第132页/共261页(3)若改用海明窗

N为奇数时

N为偶数时

上面两个表达式形式虽然完全一样,但由于N为奇数时,对称中心点α=(N-1)/2为整数,N为偶数时,α为非整数,因此N在奇数和偶数情况下,滤波器的单位脉冲响应的对称中心不同,在0≤n≤N-1上的取值也完全不同。第133页/共261页MATLAB工具箱函数实现设计指标---fir1()fir2()

fir1是用窗函数法设计线性相位FIRDF的工具箱函数调用格式如下:

hn=fir1(N,wc,′ftype′,window) fir1实现线性相位FIR滤波器的标准窗函数法设计。“标准”是指在设计低通、高通、带通和带阻FIR滤波器时,Hd(ejω)分别取相应的理想低通、高通、带通和带阻滤波器,故而设计的滤波器的频率响应称为标准频率响应。第134页/共261页

hn=fir1(N,wc)可得到截止频率为wc的N阶(单位脉冲响应h(n)长度为N+1)FIR低通滤波器,默认(缺省参数windows)选用hammiing窗。其单位脉冲响应h(n)为

h(n)=hn(n+1),n=0,1,2,…,N

而且满足线性相位条件:

h(n)=h(N-1-n)

其中wc为对π归一化的数字频率,0≤wc≤1。 当wc=[wc1,wc2]时,得到的是带通滤波器,其通带为wc1≤ω≤wc2。第135页/共261页hn=fir1(N,wc,′ftype′)可设计高通和带阻滤波器。

·当ftype=high时,设计高通FIR滤波器;

·当ftype=stop时,设计带阻FIR滤波器。

注意:在设计高通和带阻滤波器时,阶数N只能取偶数(h(n)长度N+1为奇数)。不过,当用户将N设置为奇数时,fir1会自动对N加1。

hn=fir1(N,wc,window)可以指定窗函数向量window。如果缺省window参数,则fir1默认为hamming窗。可用的其他窗函数有Boxcar,Hanning,Bartlett,Blackman,Kaiser和Chebwin窗。这些窗函数的使用很简单(可用help命令查到),例如:

hn=fir1(N,wc,bartlett(N+1))使用Bartlett窗设计;

hn=fir1(N,wc,chebwin(N+1,R))使用Chebyshev窗设计。

hn=fir1(N,wc,′ftype′,window)通过选择wc、ftype和window参数(含义同上),可以设计各种加窗滤波器。第136页/共261页

Fir2可以指定任意形状的Hd(ejω),用help命令查阅其调用格式。

Fir2(N,f,m,window)f频率向量,其值在0~1之间,m

是和f相对应的所希望的幅频相应。

window缺省时自动选用Hamming窗。第137页/共261页7.3利用频率采样法设计FIR滤波器

频率取样法设计思想:基本原理:(频域采样理论)一个有限长序列,可以通过其频谱的等间隔采样值准确地恢复出原有序列。设计过程:在频域中,对理想的频率响应等间隔采样;以Hd(k)作为实际滤波器频率特性的采样值H(k),(k=0,1,…,N-1)由H(k)通过IDFT变换可求出有限长序列h(n)、H(z)、H(ejω)第138页/共261页频率采样结构误差?增加过渡点增加采样点频率采样法思路:频域方法直接型结构第139页/共261页设希望逼近的滤波器的频响函数用Hd(ejω)表示,对它在ω=0到2π之间等间

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论