1、Digtlal Signal Processing Using MATLAB第七章 FIR滤波器设计数字频率w的概念定义:其中:=2f为模拟角频率 T:抽样时间间隔,fs:抽样频率 所以数字滤波器设计必须给出抽样频率数字频率的2等价于模拟抽样频率s=2fs按照Nyquist抽样定理,基带信号的频率特性只能限于|w|ws/2=的范围数字滤波器幅度响应(1)数字低通滤波器1数字高通滤波器数字滤波器幅度响应(2)11数字带通数字带阻7.1 概论滤波器设计:给定技术要求设计系统设计步骤:n确定技术要求:由具体应用条件决定n提供一个逼近要求的滤波器的表述n根据表述实现滤波器下面讨论时我们均假设技术要求已
2、知7.1.1 技术要求的给定幅度要求:n绝对指标要求:对幅度响应|H(ejw)|给出要求n相对指标要求:以分贝dB形式给出相位要求:线性相位一、绝对指标要求(1)Passband ripple通带波纹Stopband ripple阻带波纹Transition band过渡带绝对指标(2)频带 0,wp 称为 通带passband,1 是在理想通带响应上可以接受的容度(或波纹)频带ws,pi 称为 阻带stopband,2 是相应的阻带容度(或波纹)频带wp,ws 称为 过渡带transition band,在这个频带内幅度响应不作要求二、相对指标要求(1)0相对指标(2)Rp:以dB计的通带波
3、纹As:以dB计的阻带衰减两种指标之间的关系:Rp和As的计算见P214ex7.1&ex7.2三、为什么只讨论低通滤波器(LPF)上述指标都是针对低通滤波器的其他类型的频率选择性滤波器(如高通或带通)也能给出类似要求滤波器设计最重要的参数是频带容限和频带边缘频率四、技术指标举例设计一个低通滤波器,它具有一个通带 0,wp,通带内频带容限为1(或Rp,单位 dB),一个阻带ws,pi,阻带内容度为2(或As,单位dB)最后求得结果是得出滤波器的系统函数H(z)或差分方程五、FIR滤波器的优点相位响应可以真正线性系统绝对稳定,设计相对容易高效实现可用DFT实现实际应用时,我们感兴趣的是线性相位的F
4、IR滤波器六、线性相位响应的优点设计问题中仅有实数运算时延固定,没有时延失真对长为M的滤波器,运算次数只有M/2量级7.2 线性相位FIR滤波器性质包括脉冲和频率响应的形状,系统函数零点的位置设h(n)是长为M的脉冲响应,0nM-1,则在原点z=0处有(M-1)阶零点,在z平面其它处有 M-1个零点,频率响应函数可写为线性相位的脉冲响应形状(1)因为频率响应函数具有线性相位这里是恒定相位延迟(constant phase delay),由第6章知,h(n)是对称脉冲响应因此,h(n)关于对称,根据M的奇偶有两种对称类型线性相位的脉冲响应形状(1)线性相位的脉冲响应形状(2)第二类线性相位满足条
5、件相位响应不通过原点,但斜率恒为常数,此时称群时延(group delay),可知h(n)是反对称脉冲响应h(n)仍然关于对称,根据M的奇偶有两种对称类型线性相位的脉冲响应形状(2)对应频率响应特性H(ejw)将M为奇和偶数结合对称和反对称的情况,得到四种类型的线性FIR滤波器对应每种类型其频率响应特性都有独特性质,令其中,Hr(w)是连续的振幅响应函数,可正可负的实函数 相位响应是一个不连续函数例:设脉冲响应为h(n)=1,1,1,1,求出并画出频率响应解:频率响应函数为由方程可得:I类线性相位:对称脉冲响应,M为奇数这种情况下,beta=0,alpha=(M-1)/2是整数 h(n)=h(
6、M-1-n),0nM-1将两式比较可得:II类线性相位:对称脉冲响应,M为偶数这种情况下,beta=0,alpha=(M-1)/2不是整数 h(n)=h(M-1-n),0nM-1注意:Hr(pi)=0,因此不能采用这种类型设高通or带阻滤波器III类线性相位:反对称脉冲响应,M为奇数这种情况下,beta=pi/2,alpha=(M-1)/2是整数 h(n)=-h(M-1-n),0nM-1Hr(0)=Hr(pi)=0,因此这种滤波器不适合设计低通或高通滤波器exp(jpi/2)=j,这种特性非常适合设计希尔伯特变换器和微分器IV类线性相位:反对称脉冲响应,M为偶数这种情况和II类似,有Hr(0)
7、=0 and exp(jpi/2)=j.因此这种类型适合用于设计数字希尔伯特变换器和微分器MATLAB实现Hr_type1:求I类线性相位的Hr(w)n调用格式:Hr,w,a,L=Hr_type1(h)Hr_type2:求II类线性相位的Hr(w)n调用格式:Hr,w,b,L=Hr_type2(h)Hr_type3:求III类线性相位的Hr(w)n调用格式:Hr,w,c,L=Hr_type3(h)Hr_type4:求IV类线性相位的Hr(w)n调用格式:Hr,w,d,L=Hr_type4(h)小结了解了线性相位FIR滤波器的各种特性,便可根据实际需要选择合适的FIR滤波器,同时设计时要遵循有关
8、约束条件。如:第3、4种情况,对于任何频率都有固定的/2相移,一般微分器及90相移器采用这两种情况,而选频性滤波器则用第1、2种情况。(1)设计线性相位的低通Digtal Filter从幅度特性考虑,只能选择第1种或第2种第一种:第二种(2)设计线性相位的高通DF从幅度特性看,可用第一种或第四种第一种第四种(3)设计线性相位的带阻DF从幅度特性考虑,只能选择第一种(4)设计线性相位的带通DF从幅度特性考虑,可以选择任一种线性相位滤波器的零点位置对实序列而言,零点是共轭出现的;对对称序列而言,零点是镜像出现的;n令q=z 1,f(q)的系数与f(z)刚好倒序.n由于h(n)的系数是对成的,倒序并
9、不会改变系数.n如果zk是多项式的根,则pk=zk-1也是.对称系数多项式的镜像零点如果 zk 满足多项式:h0+h1zk-1+h2zk-2+.+hM-2zk-M+2+hM-1zk-M+1=0 此时 hM-1=h0,hM-2=h1,那么 rk=zk 1 同样会满足方程 h0+h1rk+h2rk2+h1rkM-2+h0rkM-1=h0zkM-1+h1zkM-2+h2zk2+h1zk+h0=zkM-1(h0+h1zk-1+h1zk-M+2+h0zk M+1)=0 1/z11/conj(z1)z1conj(z1)特殊的如果零点为实数,则只有两个零点:z2,1/z2如果零点在单位圆上且为虚数,则只有两
10、个零点z3,z3*如果零点在单位圆上且为实数,则只有一个零点z47.3 窗口设计法设计步骤n给定要求设计的理想滤波器的频率响应Hd(ejw)n设计一个FIR滤波器频率响应H(ejw)n由于设计是在时域中进行,使所设计滤波器的h(n)去逼近理想单位取样响应序列hd(n)理想滤波器的频率响应Hd(ejw)设希望设计的滤波器传输函数为Hd(ejw),对应脉冲响应为hd(n),则它们满足关系:若已知Hd(ejw),即可求出hd(n),再经过z变换,就可以求出系统函数H(z),从而设计出系统一般情况下,Hd(ejw)逐段恒定逐段恒定,在边界频率处有不连续点边界频率处有不连续点,因而hd(n)是无限时宽无
11、限时宽的,且是非因果序列。例:理想低通滤波器的传输函数Hd(ejw)无失真的理想低通的传输函数为相应的单位取样响应hd(n)由上式可知,hd(n)无限长,且为非因果序列理想低通滤波器的Hd(ejw)和h(n)波形设实际实现的低通滤波器单位取样响应为h(n),长为N,其系统函数设计过程相当于找到一个有限长序列h(n),去逼近理想低通的hd(n),这必然会引入误差频域的吉布斯(Gibbs)效应(截断效应)后果:引起通带和阻带内的波动效应,尤其是使阻带衰减减小设计实现一个FIR滤波器H(ejw)例:设计截止频率wc=/3时延为6的具有线性相位的FIR低通滤波器为了构造一个长为N的线性相位滤波器,只有
12、将hd(n)截取一段,并保证对(N-1)/2对称设截取的段用h(n)表示,则其中W(n):长为N的窗函数(这里取矩形序列)当=(N-1)/2时,截取的h(n)对(N-1)/2对称,保证设计的滤波器具有线性相位这里,hd(n)是以n=6为中心偶对称的无限长序列现用一个有限长N=13的因果序列h(n)逼近它最简单的方法:给hd(n)加矩形窗RN(n),即令W(n)=RN(n),则低通滤波器脉冲响应波形截断处理示意图截断处理后,由于h(n)满足对称脉冲响应,所以一定满足第一类线性相位设计步骤先由Hd(ejw)求付里叶反变换hd(n).砍头去尾。n因为我们要设计FIR滤波器h(n)必须满足:n因果性:
13、t砍头n线性相位:要求h(n)中心对称或反对称,由于砍头,所以必须去尾,让它们中心对称中心对称。即用有限长的h(n)去逼近无限长的hd(n).利用卷积过程。即h(n)=W(n)hd(n)可见窗函数序列的形状及长度的选择是设计关键。窗口法设计数字滤波器主要任务:寻找最有效的方法截断hd(n),即用一个有限长度的窗口函数序列W(n)来截取hd(n),使H(ejw)最逼近Hd(ejw)通过加窗可得到不同类型的数字滤波器n数字低通n数字高通n数字带通n数字带阻数字低通设h(n)是长为N,以=(N-1)/2为中心偶对称的函数h(n)的设计根据前面讨论可知,低通滤波器只能选择对称脉冲响应当N为奇数时,设计
14、第一种情况的线性相位低通DF当N为偶数时,设计第二种情况的线性相位低通DF设选用矩形窗,即设计举例:用矩形窗设计截止频率wc=/3的具有线性相位的FIR低通滤波器若取N=13,为奇数,则对称脉冲响应若取N=12,为偶数,则数字高通理想的线性相位高通DF的频率特性为:其幅度特性:1冲激响应理想高通滤波器冲激响应加窗处理后的数字滤波器冲激响应分析因为 为偶函数,W(n)为常数当N为奇数时,对应第一种线性相位,h(n)=h(N-1-n)为对称脉冲响应当N为偶数时,对应第四种线性相位,h(n)=-h(N-1-n)为反对称脉冲响应取矩形窗时,W(n)=RN(n)取N=13,为奇数,则取N=12,为偶数,
15、则理想数字带通滤波器理想的线性相位带通DF的频率特性为:1其冲激响应hd(n)加矩形窗处理后,得到分析:若选择相位有相移的理想带通DF频率特性为:此时的脉冲响应加矩形窗处理后分析此时的h(n)一定为反对称序列当N为奇数时,对应第三种线性相位,当N为偶数时,对应第四种线性相位,7.4 加窗对系统频率响应的影响根据频域卷积定理,加窗后,滤波器的频率响应现在我们以低通滤波器为例来讨论:n加窗后,频率响应发生了什么变化n加什么样的窗,可以使变化减至最小7.4.1 矩形窗矩形窗口的频率特性为用幅度响应和相位响应的乘积表示为矩形窗(2)当w很小时,当w很大时,WR(w)为周期函数主瓣矩形窗处理后的频率响应
16、根据频域卷积定理可得Wr(w-wc)Wr(w-wc+2/N)Wr(w-wc-2/N)加窗后的低通滤波器频谱几个特殊频率点w=0处,响应值 为窗函数频谱Wr(w-)和理想低通滤波器频率特性Hd()的乘积的积分,可近似看作Wr()在-到的全部积分面积w=wc处,Hd()刚好与Wr(w-)的一半重叠,因此H(wc)=0.5H(0)w=wc-2/N处,Wr(w-)的全部主瓣在Hd()的通带之内,因此卷积结果有最大值,频率响应出现正肩峰w=wc+2/N处,Wr(w-)的全部主瓣在Hd()的通带之外,通带内的旁瓣负的面积大于正的面积,因此卷积结果有最负值,频率响应出现负肩峰几个特殊频率点(2)当wwc+2
17、/N后,Wr(w-)的左边旁瓣的起伏部分扫过通带,卷积值围绕零值而波动当wwc+2/N时,Wr(w-)的右旁瓣进入通带,卷积值围绕H(0)而波动加矩形窗对理想低通滤波器的影响使理想频率特性不连续点处边沿加宽,形成过渡带,过渡带的宽度等于窗的频谱主瓣宽度4/N在截止频率wc的两边 处,H(w)出现肩峰,肩峰的两侧形成起伏振荡,振荡幅度取决于旁瓣相对幅度,振荡多少,取决于旁瓣的多少增加截取长度,则主瓣附近的窗的频率响应 可见改变N,只能改变窗的主瓣宽度,w坐标的比例和Wr(w)的绝对值大小,而不能改变主瓣与旁瓣的相对比例矩形窗的频谱示意图各种窗函数矩形窗阶段造成府肩峰为8.95%,阻带最小衰减为2
18、1dB,不符合工程需要为了加大阻带衰减,只能改善窗函数形状,使窗谱尽量逼近冲击函数,即绝大部分能量集中在频谱中点一般窗函数满足两项要求:n窗谱主瓣尽可能的窄,以得到较陡的过渡带n尽量减少最大旁瓣的相对幅度n一般而言,上面两项要求不能同时满足矩形窗三角形(BARTLETT)窗升余弦窗(汉宁Hanning窗)-1左移右移倒余弦升余弦窗(汉宁Hanning窗)-2升余弦窗(汉宁Hanning窗)-3由于频谱是由三个互有频移的不同幅值的矩形窗函数相加而成,这样使旁瓣大大抵消,从而能量相当有效地集中在主瓣内。其代价:主瓣加宽一倍,可达到减少肩峰,余振,提高阻带衰减。缺点:过滤带加大改进的升余弦窗(汉明H
19、anning窗)-1其频谱函数为其幅度函数为改进的升余弦窗(汉明Hamming窗)-2二阶升余弦窗(布拉克曼Blackman窗)凯塞窗(Kaiser窗)以上几种窗函数是各以一定主瓣加宽为代价,来换取某种程度的旁瓣抑制,而凯窗则是:全面地反映主瓣与旁瓣衰减之间的交换关系,可以在它们两者之间自由地选择它们的比重。滤波器阶数(长度)M的选择名称近似过渡带宽精确过渡带宽最小阻带衰减矩形4/M1.8/M21dB巴特利特8/M6.1/M25dB汉宁8/M6.2/M44dB哈明8/M6.6/M51dB布莱克曼12/M11/M74dB取Kaiser窗时设定beta,再用kaiserord函数求得MMatlab
20、 实现W=boxcar(M):产生M点的矩形窗W=triang(M):产生M点的Bartlett窗W=hanning(M)产生M点的Hanning窗W=hamming(M)产生M点的Hamming窗W=blackman(M)产生M点的Blackman窗W=kaiser(M,beta)产生beta值的M点Kaiser窗Examples例:设计一个数字FIR低通滤波器,技术指标如下:wp=0.2,Rp0.25dB,ws0.3,As50dB首先查表选择满足阻带衰减的窗函数从中选择最合适的(过渡带较小的)调用MATLAB函数进行设计验证设计滤波器是否满足通带波纹如果不能满足,则换用过渡带较大的,重复以
21、上步骤如果满足,则设计成功频率采样设计法(1)设计原理:系统函数H(z)能够从频率响应H(ejw)的样本H(k)中恢复基本思想:已知理想低通滤波器Hd(ejw),选取滤波器长度为M,在0,2区间以M等分频率对Hd(ejw)采样得H(k),再由其离散傅里叶反变换h(n)得系统函数H(z)h(n)=IDFTH(k)可用函数ifft 计算频率采样设计法(2)特点:n采样频率点上近似误差为0n其它频率点上的近似误差取决于理想响应的形状,理想响应愈陡峭,近似误差越大n靠近通带边缘的误差较大,通带内误差较小分类:n直接设计法:直接利用基本思想,在近似误差上不给出任何条件n最优设计法:通过改变过渡带内的样本
22、值将阻带内误差减至最小Phase for Type 1&2Phase for Type 3&4直接设计法(Naive design methods)设计思想:n令H(k)=Hd(e j2k/M),k=0,1,M-1,n用h(n)=IDFTH(k)求得脉冲响应h(n)例:用频率采样法设计一个数字FIR低通滤波器,技术指标如下:wp=0.2,Rp0.25dB,ws0.3,As50dB分析:取M=20,使在wp处有一个样本,即k=2 wp=0.2=(2/20)2 下一个样本在ws,即在k=3 wp=0.3=(2/20)3 这样通带内0wwp内有3个样本,在阻带wsw内有7个样本Hr(k)=1,1,1
23、,0,0,1,1共15个零由于M=20,=(M-1)/2=9.5,为II类线性相位滤波器再由IDFT可得h(n)MATLAB编程解得M=20;alpha=(M-1)/2;l=0:M-1;w1=(2*pi/M)*l;Hrs=1,1,1,zeros(1,15),1,1;Hdr=1,1,0,0;wdl=0,0.25,0.25,1;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);db,mag,pha,gr
24、d,w=freqz_m(h,1);Hr,ww,a,L=Hr_Type2(h);最优设计法(Optimum design method)设计方法:增大取样点数M,并让过渡样本作为自由样本,改变他们的值以得到在给定M的条件下的最大衰减及过渡带宽例:利用最优设计法设计一个比上例更好的低通滤波器增加取样点数M=40,以使过渡带内(0.2w0.3)有一个样本,在k=5和k=35处,用T表示这两个样本值,其中0T1,则以采样的振幅响应Hr=1,1,1,1,1,T,0,0,T,1,1,1,1共29个零由于alpha=(M-1)/2=19.5,相位响应的样本是现在我们考虑如何选取T值,以得到更好的最小阻带衰减
25、首先我们选取通带和阻带幅度的中值0.5用MATLAB编程解得MATLAB程序M=40;alpha=(M-1)/2;l=0:M-1;w1=(2*pi/M)*l;Hrs=1,1,1,1,1,0.5,zeros(1,29),0.5,1,1,1,1;Hdr=1,1,0,0;wdl=0,0.25,0.25,1;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);db,mag,pha,grd,w=freqz_m(
26、h,1);Hr,ww,a,L=Hr_Type2(h);结论:通过改变一个样本值,我们得到一种更好的设计实际系统的过渡带往往很小,只有一到两个样本,只需优化较少的样本就可以获得最大的最小阻带衰减这等小于使最大旁瓣幅度最小化,因此这类优化问题也称最大最小化问题(minimax problem)最优过渡值表见文献19的附录B最优等波纹设计法窗口设计法和频率采样设计法的缺陷n设计过程不能将边缘频率wp和wc精确给定n不能够同时标定波纹因子1和2,n近似误差在频带区间上不是均匀分布的,靠近频带边缘误差愈大,远离频带边缘误差愈小上述缺陷的克服对线性相位FIR滤波器而言,可以导出一组条件,使最大近似误差最小
27、化的意义下设计的解达到最优(最大值最小误差或Chebyshev误差)满足这种条件的滤波器称为等波纹滤波器,其近似误差在通带和阻带均匀分布,且实现相同性能滤波器时阶数更低最大最小问题的建立线性相位FIR滤波器4种情况的频率响应都能写为如下形式其中 beta 和 Hr(w)的表达式在表 7.2(P.265)中给出利用三角函数恒等式,可将上面每个Hr(w)表达式写成一个w的函数Q(w)和一个余弦和的函数P(w).其中四种情况下的Q(w),L和P(w)由表7.3(P.279)给出Chebyshev 近似问题分析的目的是为了对4种情况有Hr(w)的共同形式.这将使问题的描述更容易为了将问题归结为Cheb
28、yshev近似问题,必须定义期望的振幅响应Hdr(w)和通带与阻带内定义的加权函数W(w)(用以独立控制)加权误差定义:Chebyshev 近似问题加权误差响应E(w)与加权函数W(w)关系详见P266定义Chebyshev 近似问题的陈述给定准确的wp,ws,1和2,确定一组系数a(n)、b(n)、c(n)或d(n)以使在通带内和阻带内E(w)的最大绝对值最小,即极值数目的的确定问题的提出:对某一给定的M点滤波器而言,在误差函数E(w)内存在多少个局部最大值和最小值结论在P(w)表达式中交错点(Alternation)定理设S为闭区间0,pi内任意闭合子集,为使P(w)是在S上对Hdr(w)
29、的唯一最大值最小近似,其充要条件是E(w)在S内至少出现(L+2)个交错点(alternations)或极值频率,即在S内一定存在(L+2)个频率wi 使之满足最优等波纹滤波器在S内的误差函数要么有(L+2)个,要么有(L+3)个Parks-McClellan 算法交错点定理确保最大最小近似问题的解存在且唯一,Parks-McClellan 算法完成求解工作它由Remez交换算法提供迭代解.1.滤波器的阶数M由(7.48)计算2.估猜极值频率wi(i=1:L+2)3.根据这些极值频率点拟合一个L阶多项式4.在一个很细的密度上确定局部最大误差,并在新的极值上调整得到新的极值频率wi5.重复第3步
30、6.一直到最优一组频率和全局最大误差找到为止,最后得出多项式P(w),再确定系数(n),计算出a(n)和脉冲响应h(n)等波纹设计函数remez调用格式:h=remez(N,f,n,weights,ftype)当 weight=1,同时ftype不是Hilbert滤波器或微分器时h=remez(N,f,m)nh 是滤波器脉冲响应的系数,长为 M=N+1nN 定义了滤波器的阶数nf 一个数组,定义了以 为单位的频带边缘频率.nm 一个数组,定义了在f 给定频率上的期望幅度响应用Remez 函数作等波纹设计实例例题7.23:低通滤波器设计n同窗函数设计法比较(例7.8)n同频率采样设计法比较(例7.14,7.15,7.16)例题7.24:带通滤波器设计n同窗函数设计法比较(例7.10)n同频率采样设计法比较(例7.17)例题7.25:高通滤波器设计例题7.26:阶梯形滤波器设计其它等波纹设计实例例题7.27:使用remez函数设计数字微分器例题7.28:使用remez函数设计数字Hilbert滤波器作业(P282):p7.14 p7.19