3数字信号处理讲义ch5_2fir窗函数法

1、数字信号处理(Digital Signal Processing)信号与系统系列课程组国家电工电子教学基地FIR数字滤波器线性相位FIR数字滤波器的性质窗函数法设计FIR数字滤波器频率取样法设计线性相位FIR数字滤波器线性相位FIR数字滤波器的优化设计利用MATLAB设计FIR数字滤波器窗函数法设计FIR数字滤波器设计方法吉伯斯现象常用窗函数Kaiser(凯泽)窗利用Kaiser窗实现多带滤波器FIR窗函数法吉伯斯(1839-1903 )出生于康涅狄格新港口市，1854-1858年在耶鲁大学学习,获学士学位。1863年(美国提供授予博士学位的第3年),在耶鲁大学的设菲尔德学院获得工程学博士学1

2、871- 1903在耶鲁作数学物理学的教授（18711880年约翰斯霍普金斯雇用他时没有工资）。1879年，被选入国家科学协会。FIR窗函数法FIR窗函数法FIR窗函数法问题提出1WWjkWh k =jH (e)eddd2hd k一般情况下是非因果无限长序列，需对其进行截短和因果化处理。FIR窗函数法问题：理想滤波器的Hd (ejW )已知，设计一物理可实现的滤波器，使其频率响应逼近Hd (ejW )。设计方法可实现性处理方案方案一：只能设计I型或III型线性相位FIR滤波器(1) 将hdk 关于k =0对称截短，即hk= hdk, -Lk L(2) 将其右移L使其成为因果系统hk= hdk

3、- L， 0k 2LFIR窗函数法设计方法可实现性处理方案方案二：可设计四种类型的线性相位FIR滤波器(1) 确定线性相位因子ej(-0.5MW+b)和幅度函数Ad (ejW(2) 计算hdk ，取其在0 k M范围的值hk= hdk， 0k M)FIR窗函数法设计方法设计结果分析积分平方误差由Parseval定理12dW(e jW ) - H (e jW )e 2=Hd-2hk - hk 2e 2最小=k =-d-1Mhk 2 +hk - hk 2 +hk 2=dddk =-k =0k =M +1等于零，即hk= hdk, 0k M结论：上述方法所设计的滤波器是在积分平方误差最小意义下的最佳

4、滤波器。FIR窗函数法设计方法（采用方案二）最小积分平方误差设计FIR的步骤1. 根据要求确定线性相位FIR滤波器的类型2. 确定理想滤波器的幅度函数Ad(W)和相位fd(W)fd(W)= -0.5MW +b3. 计算IDTFT得hdk1W )ef (W )Wh k =2 WjjkA (edddd4截短hdkhk= hdk， 0k MFIR窗函数法例：设计一个幅度响应能逼近理想带通滤波器的线性相位FIR滤波器。解：(1)确定线性相位FIR滤波器类型：可选用I型或II型确定理想滤波器的幅度函数Ad(W)和相位fd(W)(2)A(W ) = 1W c 1W W c 2 0d其他fd(W)= -0.

5、5MW例：设计一个幅度响应能逼近理想带通滤波器的线性相位FIR滤波器。解：(3)计算IDTFT得hdk1(W )e jfd (W ) e jkW dWhk =Add-211-WWe jW (k -0.5M ) dW +e jW (k -0.5M ) dW=c1c2-Wc2Wc122= W c2Sa(W(k - 0.5M ) - W c1 Sa(W(k - 0.5M )c2c1(4)截断hdkhk = hd kwN k长度为N=M+1的矩形窗例：设计一个幅度响应能逼近理想带通滤波器的线性相位FIR滤波器。Wc1=0.3p，Wc2=0.7p，M=30解：hkA(W)k0Wp00.7p0.3p单位脉

6、冲响应幅度函数例：设计一个幅度响应能逼近理想带通滤波器的线性相位FIR滤波器。讨论：1.理想低通和高通滤波器是带通滤波器的特殊情况hk = W c2 Sa (W(k - 0.5M ) - W c1 Sa (W(k - 0.5M )dc2c1 当Wc1=0，滤波器为理想低通hk = W c2 Sa(W(k - 0.5M )dc2 当Wc2=p，滤波器为理想高通(M必须为偶数)hk = d k - 0.5M - W c1 Sa(W(k - 0.5M )dc1吉伯斯现象矩形窗设计的FIR低通滤波器在不同阶数时的幅度函数(Wc=p/2)0-10-21-30-4000.250.50.751所设计滤波器的

7、幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的9%，这个现象称为Gibbs 现象。M=14 M=30吉伯斯现象产生吉伯斯现象原因分析由于所设计滤波器的单位脉冲响应hk = hd kwN kN=M+1利用DFTF的性质可得所设计FIR滤波器的幅度函数1A(W ) =Ad (q )W (W - q )dq2 -A(W)逼近Ad(W)的好坏，取决于窗函数的W(W)FIR窗函数法吉伯斯现象矩形窗的幅度函数W(W)矩形窗函数的频谱sin( NW / 2)sin( W / 2)W (e jW ) = e-jW ( N -1) / 2W (W )1. W(W)的主瓣宽度4p/NN2. 旁瓣相对

8、衰减为常数3pWN2p4pp0NNFIR窗函数法吉伯斯现象矩形窗设计FIR滤波器的幅度函数W (W - q )0 W Wc- 2 / NA(q )dq- W- WccA( W )1W- - WcA(W)主要由主瓣的面积确定A(W)的波动F由IR窗旁函瓣数法引起吉伯斯现象矩形窗设计FIR滤波器的幅度函数W (W - q )Wc- 2 / N W W c+ 2/ NAd (q )qW- W- A( W )cc10.5W- - WcA(W)完全由旁瓣的面积确定FIR窗函数法吉伯斯现象矩形窗设计FIR滤波器的幅度函数A( W )1W- - WWcc4N矩形窗截断产生的波峰大约是9%，故用矩形窗设计出的

9、滤波器阻带最大衰减-20lg (9%)21dBFIR窗函数法0.5吉伯斯现象结论窗函数的主瓣宽度决定了H(ejW )过渡带的宽度,1.窗函数长度N增大，过渡带减小。2.旁瓣的大小决定了FIR滤波器在阻带的衰减用矩形窗设计出的滤波器阻带最大衰减为20lg (9%)-21dB如何提高阻带衰减?选用旁瓣幅度较小的窗函数FIR窗函数法常用窗函数矩形窗1.091Ap =-20lg(1-dp) 0.82dB,As = -20lg(ds ) 21dB0.5dp=ds=0.090.09WWcFIR窗函数法1.8Nwk = 10 k M0其他常用窗函数Hann(汉纳)窗(w=hanning(M+1)1.0064

10、dp=ds =0.0064wk111-0.00640.0064WkWc6.2 N0102030A A 44dB0.056dB,psFIR窗函数法wk = 0.5 - 0.5 cos(2k / M )0 k M0其他由Hann窗与矩形窗的比较(Wc=p/2 , M=38)0-20-44-60-8000.250.50.751Gain dBSquare Hann常用窗函数Hamming(哈明)窗(w=hamming(M+1)1.0022dp=ds =0.0022wk111-0.00220.0022Wk0102030WcApAs0.019dB,53dBFIR窗函数法7.0Nwk = 0.54 - 0.

11、46 cos(2k / M )0 k M0其它Hamming窗与Hann窗的比较(Wc=p/2 , M=38)0-20-44-52-7000.10.20.30.40.50.60.70.80.91HannHamming常用窗函数Blackman窗(w=blackman(M+1)1.0002dp=ds =0.0002wk11-0.00020.0002Wk0102030WcAp0.0017dB，As74dBFIR窗函数法11.4Nwk = 0.42 - 0.5 cos(2k / M ) + 0.08 cos(4k / M )0 k M0其它Blackman窗与Hamming窗的比较(Wc=p/2 ,

12、 M=38)0-20-52-75-10000.10.20.30.40.50.60.70.80.91 Hamming BlackmanKaiser(凯泽)窗( w=kaiser(M+1,beta) )I( b1 - 1 - 2k / M 2 )定义wk =0 k M 0,I 0 ( b )b 是一可调参数，调节窗函数的形状。I0(x ):the modified zeroth-order Bessel function.I0(x )可用幂级数表示为2 ( x / 2)n I 0 ( x) = 1 + n!n=1一般求20项就能达到所需精度。FIR窗函数法Kaiser(凯泽)窗b取不同值时窗的形状

13、b=010.500246810b=1121416182010.500246810b=3121416182010.500246810FIR窗函数法1214161820Kaiser(凯泽)窗M与b 的确定0.1102( A - 8.7),0.5842( A - 21)0.4A 5021 A 50A 21,2.285 W p - WsFIR窗函数法例：用Kaiser窗设计满足下列指标的I型线性相位FIR低通滤波器。Wp=0.2p, Ws=0.4p，Ap=0.3dB, As=50dB解：b = 0.1102(As-8.7)估计b1.2.估计滤波器的阶数MAs - 7.95M 2.285 W p- Ws

14、3. 设定理想低通的截频Wc = (Wp4.hk=hdk wk+ Ws) / 2wp=0.3*pi;ws=0.4*pi;As=50;M=ceil(As-7.95)/(ws-wp)/2.285) M=M+mod(M,2)beta=0.1102*(As-8.7);w=kaiser(M+1,beta); wc=(wp+ws)/2; alpha=M/2;k=0:M;hd=(wc/pi)*sinc(wc/pi)*(k-alpha); omega=linspace(0,pi,512); mag=freqz(h,1,omega); magdB=20*log10(abs(mag); plot(omega/pi

15、,magdB);axis(0,1,-70,0);grid;h=hd.*w;例：用Kaiser窗设计满足下列指标的I型线性相位FIR低通滤波器。Wp=0.2p, Ws=0.4p，Ap=0.3dB, As=50dB解：设计结果Gain response of lowpass FIR filter0-30-50M=60Ap=0.0105 dBAs=50.7524 dB-8000.10.20.30.40.5Normalized frequencyGain, dB例：理想数字微分器的频率响应为HDIF(ejW)=jW,|W|p。试用窗口法设计一线性相位FIR滤波器，使其幅度响应逼近理想数字微分器。解：1

16、)确定线性相位FIR滤波器类型：由于HDIF(ejW)=jW，含有90相移，故选用III型或IV型2)确定理想滤波器的幅度函数Ad(W)和相位fd(W) = WW (e jWf (W)= -0.5MW+0.5pA0dd其他3)计算IDTFT得hdk1(W )e jfd (W ) e jkW dWk =hADIFd-2例：理想数字微分器的频率响应为HDIF(ejW)=jW,|W|p。试用窗口法设计一线性相位FIR滤波器，使其幅度响应逼近理想数字微分器。解：3)计算IDTFT得hdkjW e jW ( k -0.5 M ) dWk =hDIF-2(-1)( k -0.5 M )M为偶, k 0.5

17、M M为偶, k = 0.5MM为奇,k - 0.5M=0,(-1)( k -0.5 M +0.5 ) p(k - 0.5M ) 2hk = hDIFkwN k4)截断hDIFk例：理想数字微分器的频率响应为HDIF(ejW)=jW,|W|p。试用窗口法设计一线性相位FIR滤波器，使其幅度响应逼近理想数字微分器。A(W)A(W)M=10M=9WW-p-p0pp0结论1.由于III型性相位FIR滤波器在W=p存在固有的零点，所以M=10的微分器的幅度函数在W=p附近误差较大。2. 由于M为偶时hDIFk衰减较慢(按|k-0.5|-1速率衰减其幅度函数有较大的波动。)，所以窗函数法的MATLAB实

18、现若FIR滤波器有4个频带，其分别为pf4Wpf5pf6Wp0Wpf1pf2Wpf3FIR滤波器在4个频带中的幅度值为a1a2a3a4(通带取1，阻带取0)FIR滤波器在4个频带中的波动d1d2d3d4则可以利用MATLAB中的kaiserord和fir1函数实现f= f1 a= a1f2f3 a2f4f5 a3f6;a4;dev= d1d2d3d4;M,Wc,beta,ftype = kaiserord(f,a,dev);h = fir1(M,Wc,ftype,kaiser(M+1,beta)FIR窗函数法例:试用Kaiser窗设计满足下列指标的具有2个通带FIR滤波器Ws1=0.1p,Wp

19、1=0.2p， Wp2=0.4p,Ws2=0.5p,Wp3=0.7p,Wp4=0.8p,Ws4=0.9p,ds=0.01Ws3=0.6p,f=0.10.20.40.50.60.70.80.9;a=0,1,0,1,0;Rs=0.01;dev=Rs*ones(1,length(a);N,Wc,beta,ftype=kaiserord(f,a,dev);h=fir1(N,Wc,ftype,kaiser(N+1,beta);omega=linspace(0,pi,512); mag=freqz(h,1,omega); plot(omega/pi,20*log10(abs(mag);xlabel(Normalizedfrequency);ylabel(Gain,db);grid;5);axis(01-80例:试用Kaiser窗设计满足下列指标的