ch7-FIR滤波器设计课件

数字信号处理

第七章FIR滤波器设计授课教师：胡双红联系电话：1FIR滤波器设计数字滤波器预备知识线性相位窗口设计法加窗对频率响应的影响频率采样设计法2第一次课数字滤波器预备知识线性相位3数字频率w的概念定义：

其中：Ω=2πf为模拟角频率T:抽样时间间隔，fs:抽样频率

所以数字滤波器设计必须给出抽样频率数字频率的2π等价于模拟抽样频率Ωs=2πfs由于滤波器的频率特性具有2π的周期性和对称性，其频率特性只限于|w|<ws/2=π的范围56数字滤波器的技术指标幅度要求：绝对指标要求：对幅度响应|H(ejw)|给出要求相对指标要求：以分贝dB形式给出7绝对指标(2)频带[0,wp]称为通带passband,δp是在理想通带响应上可以接受的容度（或波纹）频带[ws,pi]称为阻带stopband,

δs

是相应的阻带容度（或波纹）频带[wp,ws]称为过渡带transitionband,在这个频带内幅度响应不作要求9二、相对指标要求(1)010相对指标(2)Rp:以dB计的通带波纹As:以dB计的阻带衰减两种指标之间的关系：

Rp

和As的计算见P214ex7.1&ex7.211为什么只讨论低通滤波器(LPF)上述指标都是针对低通滤波器的其他类型的频率选择性滤波器(如高通或带通)也能给出类似要求滤波器设计最重要的参数是频带容限和频带边缘频率13FIR滤波器的优点相位响应可以真正线性系统绝对稳定，设计相对容易高效实现可用DFT实现实际应用时，我们感兴趣的是线性相位的FIR滤波器14线性相位响应的优点设计问题中仅有实数运算时延固定，没有时延失真对长为M的滤波器，运算次数只有M/2量级15线性相位FIR滤波器的h(n)设h(n)是长为M的脉冲响应，0≤n≤M-1，则在原点z=0处有(M-1)阶极点，在z平面其它处有M-1个零点，频率响应函数可写为17线性相位的脉冲响应形状(1)因为频率响应函数具有线性相位这里是恒定相位延迟(

constantphasedelay)，由第6章知，h(n)是对称脉冲响应因此，h(n)关于对称，根据M的奇偶有两种对称类型18线性相位的脉冲响应形状(1)19线性相位的脉冲响应形状(2)21对应频率响应特性H(ejw)将M为奇和偶数结合对称和反对称的情况，得到四种类型的线性FIR滤波器对应每种类型其频率响应特性都有独特性质，令其中，Hr(w)是连续的振幅响应函数，可正可负的实函数相位响应是一个不连续函数22例：设脉冲响应为h(n)={1，1，1，1}，求出并画出频率响应解：频率响应函数为由方程可得：23I种线性相位：对称脉冲响应，M为奇数这种情况下，beta=0，alpha=(M-1)/2是整数h(n)=h(M-1-n)，0≤n≤M-1将两式比较可得：25II种线性相位:对称脉冲响应，M为偶数这种情况下，beta=0，alpha=(M-1)/2不是整数h(n)=h(M-1-n)，0≤n≤M-1注意：Hr(pi)=0,因此不能采用这种类型设高通or带阻滤波器26III种线性相位：反对称脉冲响应，M为奇数这种情况下，beta=pi/2，alpha=(M-1)/2是整数h(n)=-h(M-1-n)，0≤n≤M-1Hr(0)=Hr(pi)=0,因此不适合设计低通或高通滤波器exp(jpi/2)=j，非常适合设计希尔伯特变换器和微分器27MATLAB实现Hr_type1:求I种线性相位的Hr(w)调用格式：[Hr,w,a,L]=Hr_type1(h)Hr_type2:求II种线性相位的Hr(w)调用格式：[Hr,w,b,L]=Hr_type2(h)Hr_type3:求III种线性相位的Hr(w)调用格式：[Hr,w,c,L]=Hr_type3(h)Hr_type4:求IV种线性相位的Hr(w)调用格式：[Hr,w,d,L]=Hr_type4(h)29小结了解了线性相位FIR滤波器的各种特性，便可根据实际需要选择合适的FIR滤波器，同时设计时要遵循有关约束条件。如：第3、4种情况，对于任何频率都有固定的π/2相移，一般微分器及90°相移器采用这两种情况，而选频性滤波器则用第1、2种情况。30(1)设计线性相位的低通DigtalFilter从幅度特性考虑，只能选择第1种或第2种第一种：第二种31(2)设计线性相位的高通DF从幅度特性看，可用第一种或第四种第一种第四种32(3)设计线性相位的带阻DF从幅度特性考虑，只能选择第一种33(4)设计线性相位的带通DF从幅度特性考虑，可以选择任一种34线性相位滤波器的零点位置对实序列而言，零点是共轭出现的;对对称序列而言，零点是镜像出现的;令q=z–1,f(q)的系数与f(z)刚好倒序.由于h(n)的系数是对成的，倒序并不会改变系数.如果zk是多项式的根，则pk=zk-1也是.35对称系数多项式的镜像零点如果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)=0361/z11/conj(z1)z1conj(z1)37特殊的如果零点为实数，则只有两个零点：z2,1/z2如果零点在单位圆上且为虚数，则只有两个零点z3,z3*如果零点在单位圆上且为实数，则只有一个零点z438第二次课窗口设计法加窗对频率特性的影响常用窗函数窗口设计法例题397.2窗口设计法设计思想设计步骤理想低通滤波器如何截断窗口设计法40设计思想时域逼近使所设计滤波器的h(n)逼近理想滤波器的hd(n)41设计步骤根据滤波器类型得理想滤波器Hd(ejw)求IDTFT得hd(n)截断得实际滤波器h(n)42理想滤波器的频率响应Hd(ejw)设希望设计的滤波器传输函数为Hd(ejw)，对应脉冲响应为hd(n)，则它们满足关系：一般情况下，Hd(ejw)逐段恒定，在边界频率处有不连续点，因而hd(n)是无限时宽的，且是非因果序列。43例：理想低通滤波器的传输函数Hd(ejw)无失真的理想低通的传输函数为相应的单位取样响应hd(n)由上式可知，hd(n)无限长，且为非因果序列44理想低通滤波器的Hd(ejw)和h(n)波形45设实际实现的低通滤波器单位取样响应为h(n)，长为N，其系统函数设计过程相当于找到一个有限长序列h(n)，去逼近理想低通的hd(n)，这必然会引入误差——频域的吉布斯(Gibbs)效应(截断效应)后果：引起通带和阻带内的波动效应，尤其是使阻带衰减减小设计实现一个FIR滤波器H(ejw)46问题：如何截？砍头去尾。要设计的FIR滤波器h(n)必须满足：因果性：t<0时，h(n)=0-->砍头线性相位：要求h(n)中心对称或反对称，由于砍头，所以必须去尾，让它们中心对称。即用有限长的h(n)去逼近无限长的hd(n).47例：设计截止频率wc=/3时延为6的具有线性相位的FIR低通滤波器为了构造一个长为N的线性相位滤波器，只有将hd(n)截取一段，并保证对(N-1)/2对称设截取的段用h(n)表示，则其中W(n)：长为N的窗函数(这里取矩形序列)当τ=(N-1)/2时，截取的h(n)对(N-1)/2对称，保证设计的滤波器具有线性相位48这里，hd(n)是以n=6为中心偶对称的无限长序列现用一个有限长N=13的因果序列h(n)逼近它最简单的方法：给hd(n)加矩形窗RN(n)，即令W(n)=RN(n)，则49低通滤波器脉冲响应波形截断处理示意图截断处理后，由于h(n)满足对称脉冲响应，所以一定满足第一类线性相位50小结先由Hd(ejw)求付里叶反变换hd(n).截断得h(n)后果：吉布斯效应性能分析方法时域乘积对应频域卷积h(n)=W(n)×hd(n)可见窗函数序列的形状及长度的选择是设计关键51窗口法主要任务寻找最有效的方法截断hd(n)即用一个有限长度的窗口函数序列W(n)来截取hd(n)，使H(ejw)最逼近Hd(ejw)下面分析窗口形状和长度对频响的影响527.3加窗对系统频率响应的影响根据频域卷积定理，加窗后，滤波器的频率响应现在我们以低通滤波器为例来讨论：加窗后，频率响应发生了什么变化加什么样的窗，可以使变化减至最小537.3.1矩形窗矩形窗口的频率特性为用幅度响应和相位响应的乘积表示为54矩形窗（2）当w很小时，当w很大时，WR(w)为周期函数主瓣55矩形窗处理后的频率响应根据频域卷积定理可得56H(wc)=0.5H(0)H(wc-2π/N)H(wc+2π/N)Hd(ejθ)W(ej（w-θ)H(0)卷积就是乘积下的面积57加窗后的低通滤波器频谱58几个特殊频率点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(θ)的通带之外，通带内的旁瓣负的面积大于正的面积，因此卷积结果有最负值，频率响应出现负肩峰59几个特殊频率点（2）当w>wc+2π/N后，Wr(w-θ)的左边旁瓣的起伏部分扫过通带，卷积值围绕零值而波动当w<wc+2π/N时，Wr(w-θ)的右旁瓣进入通带，卷积值围绕H(0)而波动60加矩形窗对理想低通滤波器的影响使理想频率特性不连续点处边沿加宽，形成过渡带，过渡带的宽度等于窗的频谱主瓣宽度4π/N在截止频率wc的两边处，H(w)出现肩峰，肩峰的两侧形成起伏振荡，振荡幅度取决于旁瓣相对幅度，振荡多少，取决于旁瓣的多少增加截取长度，则主瓣附近的窗的频率响应

可见改变N，只能改变窗的主瓣宽度，w坐标的比例和Wr(w)的绝对值大小，而不能改变主瓣与旁瓣的相对比例61矩形窗的频谱示意图626364各种窗函数矩形窗截断造成负肩峰为8.95%，阻带最小衰减为21dB，不符合工程需要为了加大阻带衰减，只能改善窗函数形状，使窗谱尽量逼近冲击函数，即绝大部分能量集中在频谱中点一般窗函数满足两项要求：窗谱主瓣尽可能的窄，以得到较陡的过渡带尽量减少最大旁瓣的相对幅度一般而言，上面两项要求不能同时满足65矩形窗截断后滤波器阻带衰减最小值为21dB66三角形(BARTLETT)窗截断后滤波器阻带衰减最小值为25dB67升余弦窗（汉宁Hanning窗）-1左移右移倒余弦68升余弦窗（汉宁Hanning窗）-269升余弦窗（汉宁Hanning窗）-3频谱由三个互有频移、不同幅值的矩形窗谱相加而成，旁瓣大大抵消，能量有效集中在主瓣内。代价:主瓣加宽一倍，过渡带加大优点：减少肩峰，余振，提高阻带衰减。截断后滤波器阻带衰减最小值为44dB70改进的升余弦窗(汉明Hanning窗)-1其频谱函数为其幅度函数为7172改进的升余弦窗（汉明Hamming窗）-2截断后滤波器阻带衰减最小值为53dB73二阶升余弦窗(布拉克曼Blackman窗)74截断后滤波器阻带衰减最小值为71dB75滤波器阶数(长度)M的选择

名称窗函数主瓣宽度旁瓣峰值衰减滤波器过渡带宽最小阻带衰减矩形4π/M13dB1.8π/M21dB巴特利特8π/M26dB6.1π/M25dB汉宁8π/M31dB6.2π/M44dB汉明8π/M41dB6.6π/M51dB布莱克曼12π/M57dB11π/M74dB取Kaiser窗时设定beta,再用kaiserord函数求得M76例：设计一个FIR理想低通滤波器，wc=0.25π,

分别取N=11,21,31,观察加窗后对滤波器幅频特性的影响理想低通滤波器的脉冲响应将hd(n)截短为N=2τ+1,并将截短后的hd(n)移位并将wc=0.25π代入得77（1）加矩形窗N=11h(n)={-0.045,0,0.075,0.1592,0.2251,0.25,0.2251,0.1592,0.075,0,-0.045}N=21h(n)={0.0318,0.025,0,-0.0322,-0.0531,-0.045,0,0.075,0.1592,0.2251,0.25,0.2251,0.1592,0.075,0,-0.045,-0.0531,-0.0322,0,0.025,0.0318}N=31h(n)={-0.0150,-0.0227,-0.0173,0,0.0205,0.0318,0.025,0,-0.0322,-0.0531,-0.045,0,0.075,0.1592,0.2251,0.25,0.2251,0.1592,0.075,0,-0.045,-0.0531,-0.0322,0,0.025,0.0318,0.0205,0,-0.0173,-0.0227,-0.0150}78加汉明窗N=11h(n)={-0.003600.02980.10860.20530.250.20530.10860.02980-0.0036}N=21h(n)={0.00250.00260-0.0087-0.0211-0.024300.06080.14520.220.250.220.14520.06080-0.0243-0.0211-0.008700.00260.0025}N=31h(n)={-0.0012-0.002-0.002100.00480.00990.00990-0.0189-0.0362-0.034700.06840.15280.22280.250.22280.15280.06840-0.0347-0.0362-0.018900.00990.00990.00480-0.0021-0.002-0.0012}7980818283例：设计一个数字FIR低通滤波器，技术指标如下:fp=2kHz，fstop＝3kHz，As＝40dB，fs=10kHz查表选择满足阻带衰减的窗函数As＝40dB,至少应选择汉宁窗根据过渡带宽确定滤波器的阶数△f=1kHz，转换为数字频率过渡带宽△w=2πfc/fs=0.2π84确定截止频率和延时fc=(fp+fstop)/2=2.5kHzwc=0.5π延时τ=(N-1)/2=15所以滤波器的冲激响应为8586第三次课带宽换性能的凯塞窗滤波器设计的MATLAB实现Hilbert变换器和微分器设计频率采样设计法87凯塞窗（Kaiser窗）以上几种窗函数是各以一定主瓣加宽为代价，来换取某种程度的旁瓣抑制，而凯窗则是：全面地反映主瓣与旁瓣衰减之间的交换关系，可以在它们两者之间自由地选择它们的比重。888990Matlab实现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窗Examples91理想低通脉冲响应计算函数functionhd=ideal_lp(wc,M);%理想低通滤波器计算%--------------------------------%[hd]=ideal_lp(wc,M)%hd=0toM-1之间的理想脉冲响应%wc=截止频率(弧度)%M=理想滤波器的长度%alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha+eps;hd=sin(wc*m)./(pi*m);92计算绝对和相对幅度响应、相位响应和群时延响应的函数function[db,mag,pha,grd,w]=freqz_m(b,a);%freqz子程序的改进版本%------------------------------------%[db,mag,pha,grd,w]=freqz_m(b,a);%db=[0到pi弧度]区间内的相对振幅(db)%mag=[0到pi弧度]区间内的绝对振幅%pha=[0到pi弧度]区间内的相位响应%grd=[0到pi弧度]区间内的群迟延%w=[0到pi弧度]区间内的501个频率样本向量%b=Ha(z)的分子多项式系数(对FIRb=h)%a=Ha(z)的分母多项式系数(对FIR:a=[1])93[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);grd=grpdelay(b,a,w);94Hilbert变换器和微分器引导分析设计思路对照书上例题讲解957.4频率采样设计法设计思想数学推导直接设计法最优设计法96设计思想频域逼近原理：系统函数H(z)能够从频率响应H(ejw)的样本H(k)中恢复97设计步骤已知理想低通滤波器Hd(ejw)，选取滤波器长度为M[0，2π]对Hd(ejw)

M等分采样得H(k)离散傅里叶反变换得h(n)H(z)981.2类线性相位3.4类线性相位99PhaseforType1&2PhaseforType3&4100分类：直接设计法：直接利用基本思想，在近似误差上不给出任何条件最优设计法：通过改变过渡带内的样本值将阻带内误差减至最小101直接设计法(Naivedesignmethods)设计思想：令H(k)=Hd(ej2πk/M)，k=0,1,…,M-1，用h(n)=IDFT[H(k)]求得脉冲响应h(n)例：用频率采样法设计一个数字FIR低通滤波器，技术指标如下:wp=0.2π，Rp＝0.25dB，ws＝0.3π，As＝50dB102分析：取M=20，使在wp处有一个样本，即k=2

wp=0.2π=(2π/20)2

下一个样本在ws，即在k=3

wp=0.3π=(2π/20)3

这样通带内[0≤w≤wp]内有3个样本，在阻带[ws≤w≤π]内有7个样本

Hr(k)=[1,1,1,0,……,0,1,1]共15个零由于M=20，α=(M-1)/2=9.5，为II类线性相位滤波器再由IDFT可得h(n)103MATLAB编程解得M=20;alpha=(M-1)/2;l=0:M-1;w1=(2*pi/M)*l;Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];Hrs=[1,1,1,zeros(1,15),1,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(h,1);[Hr,ww,a,L]=Hr_Type2(h);集中定义参数理想滤波器频率特性采样得到的幅度和相位值合成得到实际H(k)IDFT得h(n)验证性能104直接设计法仿真结果105特点：采样频率点上近似误差为0其它频率点上的近似误差取决于理想响应的形状，理想响应愈陡峭，近似误差越大靠近通带边缘的误差较大，通带内误差较小106最优设计法(Optimumdesignmethod)设计方法：增大取样点数M，并让过渡样本作为自由样本，改变他们的值以得到在给定M的条件下的最大衰减及过渡带宽例：利用最优设计法设计一个比上例更好的低通滤波器增加取样点数M=40，以使过渡带内(0.2

π<w<0.3π)有一个样本，在k=5和k=35处，用T表示这两个样本值，其中0<T<1，则以采样的振幅响应Hr=[1,1,1,1,1,T,0,……,0,T,1,1,1,1]共29个零107由于alpha=(M-1)/2=19.5，相位响应的样本是现在我们考虑如何选取T值，以得到更好的最小阻带衰减首先我们选取通带和阻带幅度的中值0.5用MATLAB编程