數(shù)字信號處理課程設(shè)計窗函數(shù)法線性相位FIR數(shù)字濾波器設(shè)計

1、汕 頭 大 學(xué) 工 學(xué) 院三級項目報告課 程 名 稱 ： 數(shù)字信號處理 三級項目題目：窗函數(shù)法線性相位fir數(shù)字濾波器設(shè)計指 導(dǎo) 教 師 ： * 系 別： 電子工程系 專 業(yè)： 姓名學(xué)號完成時間： 2010年 12月12 日 成績： 評閱人：* 1 內(nèi)容與要求1、線性相位fir數(shù)字濾波器應(yīng)滿足的條件。2、頻率選擇性線性相位fir數(shù)字濾波器的設(shè)計算法。2 報告正文一、線性相移fir 數(shù)字濾波器的條件：第一類線性相位：由偶對稱h(n)=h(n-1-n),可得n=(n1)/2為h(n)的偶對稱中心幅度函數(shù)：1）h(n)偶對稱，n為奇數(shù)2）h(n)偶對稱，n為偶數(shù)第二類線性相位由奇對稱h(n)=-h(

2、n-1-n),可得n =(n1)/2為h(n)的奇對稱中心幅度函數(shù)：3）h(n)奇對稱，n為奇數(shù)4）h(n)奇對稱，n為偶數(shù)總結(jié)四種線性相位fir的特點：當h(n)為實數(shù)且偶對稱時，fir為恒相時延，相位曲線是一條過原點、以-(n-1)/2為斜率的直線。信號通過這類濾波器后，各種頻率分量的時延都是(n-1)/2。當n為奇數(shù)時，時延(n-1)/2是整數(shù)，是采樣間隔的整數(shù)倍，采樣點時延后仍是采樣點。但當n為偶數(shù)時，時延(n-1)/2不是整數(shù)，采樣點時延后不在采樣點位置上。同時，n為偶數(shù)時，點是幅度的零點，不能做高通、帶阻濾波器。一般情況下，第一類fir特別適合做各種濾波器。當h(n)為實數(shù)且奇對稱

3、時，fir僅是恒群時延。相位曲線是一條截距為2，以-(n-1)/2為斜率的直線。信號通過該濾波器產(chǎn)生的時延也是(n-1)/2個采樣周期，但另外對所有頻率分量均有一個附加的90度的相移。程序代碼：圖一：h=-4 1 -1 -2 5 6 5 -2 -1 1 -4m=length(h);n=0:m-1;hr,w,a,l=hr_type1(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脈沖響應(yīng))subplot(2,2,3);stem(0:l,a);xlabel(n);ylabel(a(n);title(a(n)系數(shù))subplot(2,

4、2,2);plot(w/pi,hr);xlabel(頻率 單位pi);ylabel(hr);title(1型幅度響應(yīng))subplot(2,2,4);pzplotz(h,1);%畫極零圖圖二：h=-4 1 -1 -2 5 6 6 5 -2 -1 1 -4m=length(h);n=0:m-1;hr,w,b,l=hr_type2(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脈沖響應(yīng))subplot(2,2,3);stem(1:l,b);xlabel(n);ylabel(b(n);title(b(n)系數(shù))subplot(2,2,2

5、);plot(w/pi,hr);xlabel(頻率 單位pi);ylabel(h);title(2型幅度響應(yīng))axis(0 2 -20 40);subplot(2,2,4);pzplotz(h,1);%畫極零圖圖三：h=-4 1 -1 -2 5 0 -5 2 1 -1 4m=length(h);n=0:m-1;hr,w,c,l=hr_type3(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脈沖響應(yīng))subplot(2,2,3);stem(0:l,c);xlabel(n);ylabel(c(n);title(c(n)系數(shù))sub

6、plot(2,2,2);plot(w/pi,hr);xlabel(頻率 單位pi);ylabel(h);title(3型幅度響應(yīng))subplot(2,2,4);pzplotz(h,1);%畫極零圖圖四：h=-4 1 -1 -2 5 6 -6 -5 2 1 -1 4m=length(h);n=0:m-1;hr,w,d,l=hr_type4(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脈沖響應(yīng))subplot(2,2,3);stem(1:l,d);xlabel(n);ylabel(d(n);title(d(n)系數(shù))subplot

7、(2,2,2);plot(w/pi,hr);xlabel(頻率 單位pi);ylabel(h);title(4型幅度響應(yīng))subplot(2,2,4);pzplotz(h,1);%畫極零圖function hr,w,a,l=hr_type1(h)%計算1型濾波器的振幅響應(yīng)hr（w）%hr=振幅響應(yīng)%w=數(shù)字頻率%a=1型低通濾波器系數(shù)%l=hr的階次%h=1型低通濾波器的脈沖響應(yīng)m=length(h);l=(m-1)/2;a=2*h(l+1:-1:1);n=0:1:l;w=0:1:1000.*pi/500;hr=cos(w*n)*a.;function hr,w,b,l=hr_type2(h)

8、%計算2型濾波器的振幅響應(yīng)hr（w）%hr=振幅響應(yīng)%w=數(shù)字頻率%b=2型低通濾波器系數(shù)%l=hr的階次%h=2型低通濾波器的脈沖響應(yīng)m=length(h);l=m/2;b=2*h(l:-1:1);n=1:1:l;n=n-0.5;w=0:1:1000.*pi/500;hr=cos(w*n)*b.;function hr,w,c,l=hr_type3(h)%計算3型濾波器的振幅響應(yīng)hr（w）%hr=振幅響應(yīng)%w=數(shù)字頻率%c=2型低通濾波器系數(shù)%l=hr的階次%h=3型低通濾波器的脈沖響應(yīng)m=length(h);l=(m-1)/2;c=2*h(l+1:-1:1);n=0:1:l;w=0:1:1

9、000.*pi/500;hr=sin(w*n)*c.;function hr,w,d,l=hr_type4(h)%計算2型濾波器的振幅響應(yīng)hr（w）%hr=振幅響應(yīng)%w=數(shù)字頻率%d=4型低通濾波器系數(shù)%l=hr的階次%h=4型低通濾波器的脈沖響應(yīng)m=length(h);l=m/2;d=2*h(l:-1:1);n=1:1:l;n=n-0.5;w=0:1:1000.*pi/500;hr=sin(w*n)*d.;function pzplotz(b,a)% pzplotz(b,a) 按給定系數(shù)向量b,a在z平面上畫出零極點分布圖% b - 分子多項式系數(shù)向量% a - 分母多項式系數(shù)向量% a,b

10、向量可從z的最高冪降冪排至z0,也可由z0開始,按z-1的升冪排至z的最負冪.n = length(a); m = length(b); pz = ; zz = ; if (n m) zz = zeros(n-m),1);elseif (m n) pz = zeros(m-n),1); endpz = pz;roots(a); zz = zz;roots(b);pzr = real(pz); pzi = imag(pz);zzr = real(zz); zzi = imag(zz);rzmin = min(pzr,zzr,-1)-0.5; rzmax = max(pzr,zzr,1)+0.5;

11、izmin = min(pzi,zzi,-1)-0.5; izmax = max(pzi,zzi,1)+0.5;zmin = min(rzmin,izmin); zmax = max(rzmax,izmax); zmm = max(abs(zmin,zmax);%uc=exp(j*2*pi*0:1:500/500); % 單位圓plot(real(uc),imag(uc),b,-zmm,zmm,0,0,b,0,0,-zmm,zmm,b);axis(-zmm,zmm,-zmm,zmm);axis(square);holdplot(zzr,zzi,bo,pzr,pzi,rx);holdtext(z

12、mm*1.1,zmm*0.95,z-平面)xlabel(實軸);ylabel(虛軸)title(零極點圖)二、頻率選擇性線性相位fir數(shù)字濾波器的設(shè)計算法及設(shè)計步驟fir數(shù)字濾波器設(shè)計的基本思想(以窗函數(shù)法為例)如果希望得到的濾波器的理想特性為 hd(ejw),那么fir 濾波器的設(shè)計就是尋求一個系統(tǒng)函數(shù)h(z),用其頻響去逼近 hd(ejw)。理想特性的單位沖激響應(yīng)可由傅里葉反變換求取, 即，為了使hd(n)變成有限長, 可將hd ( n) 利用有限長度的窗函數(shù)w(n)直接截短成有限項n , 即?；趉aiser窗的fir數(shù)字濾波器的設(shè)計步驟1.有理想低通濾波器的頻率響應(yīng)得到理想的hd(n)

13、:2.把hd(n)進行截斷移位得到有限長、因果的物理可實現(xiàn)的單位抽樣響應(yīng)h(n)即3.根據(jù)技術(shù)指標計算窗函數(shù)的設(shè)計參數(shù)n和4.將h(n)乘上kaiser窗函數(shù)，既可得到基于kaiser窗的fir數(shù)字濾波器5.濾波器的技術(shù)指標頻帶指標依采樣頻率的歸一化頻率起伏指標簡化為對于較小的進一步的工程簡化二、matlab程序代碼clear all clfy=klh(1,20000,4000,5000,0.1,80); plot(y)title(h(n) pause clf h,w=freqz(y,1,512); f=w*pi;plot(f,20*log10(abs(h)xlabel(f/khz)title

14、(frequency response of fir filter)clfy=klh(-1,20000,5000,4000,0.1,80);plot(y)title(h(n)pauseclfh,w=freqz(y,l,512);f=w*pi;plot(f,20*log10(abs(h)xlabel(f/khz)title(frequency response of fir filter)clfy=kbp(20000,5000,8000,4000,9000,0.1,80,1);plot(y)h,w=freqz(y,1,512);f=w*pi;plot(f,20*log10(abs(h) xlab

15、el(f/khz)title(frequency response of fir filter)%kaiser窗%alpha（即）和n分別是kaiser窗的控制參數(shù)function w=kwind(alpha,n)m=(n-1)/2;den=io(alpha);for n=0:n-1, w(n+1)=io(alpha*sqrt(n*(n-1-n)/m)/den;end%io為bessel方程function s=io(x)eps=10(-9);n=1;s=1;d=1;while d(eps*s), t=x/(2*n); n=n+1; d=d*t2; s=s+d;end%計算kaiser窗控制參

16、數(shù)function w=kwind(alpha,n)m=(n-1)/2;den=io(alpha);for n=0:n-1, w(n+1)=io(alpha*sqrt(n*(n-1-n)/m)/den;end%移位并截斷理想低/高通濾波器單位沖激響應(yīng)%wc為截止頻率%n為濾波器長度%s=1為低通，s=-1為高通；function h=dlh(s,wc,n)m=(n-1)/2;for k=-m:m, if k=0, h(k+m+1)=(1-s)/2+s*wc/pi; else h(k+m+1)=s*sin(wc*k)/(pi*k); endend%低/高通基于kaiser窗的fir濾波器單位沖激

17、響應(yīng)%fs為采樣頻率，fpass為通帶截止頻率，fstop為阻帶截止頻率%function h=klh(s,fs,fpass,fstop,apass,astop)fc=(fpass+fstop)/2;wc=2*pi*fc/fs;df=s*(fstop-fpass); df=df/fs;dpass=(10(apass/20)-1)/(10(apass/20)+1);dstop=10(-astop/20);d=min(dpass,dstop);a=-20*log10(d);alpha,n=kparm(df,a);h=dlh(s,wc,n).*kwind(alpha,n);end%帶通濾波器原理為高

18、低通濾波器的級聯(lián)%fs為采樣頻率，fpa為高通的通帶截止頻率%fpb為低通的通帶截止頻率%fsa,fsb分別為高/低阻帶截止頻率function h=kbp(fs,fpa,fpb,fsa,fsb,apass,astop,s)df=min(fpa-fsa,fsb-fpb);df=df/fs;fa=(1+s)*fpa+(1-s)*fsa-s*df)/2;wa=2*pi*fa/fs;fb=(1+s)*fpb+(1-s)*fsb+s*df)/2;wb=2*pi*fb/fs;dpass=(10(apass/20)-1)/(10(apass/20)+1);dstop=10(-astop/20);d=min

19、(dpass,dstop);a=-20*log10(d);alpha,n=kparm(df,a);h=dbp(wa,wb,n).*kwind(alpha,n);%理想帶通濾波器function h=dbp(wa,wb,n)m=(n-1)/2;for k=-m:m, if k=0, h(k+m+1)=(wb-wa)/pi; else h(k+m+1)=sin(wb*k)/(pi*k)-sin(wa*k)/(pi*k); endend%3 結(jié)果與分析各個技術(shù)指標如下fs=20khz fpass=4khz fstop=5khz apass=0.1db astop=80db通過仿真可以得到基于kais

20、er窗的fir濾波器的單位沖激響應(yīng)為如下:系統(tǒng)的對數(shù)頻率響應(yīng)為：此時n的階數(shù)為103由于階數(shù)n與f成反比，改變f可以實現(xiàn)n的控制；當參數(shù)截止頻率改為：fpass=4khz fstop=4.5khz 此時計算得n=203濾波器的對數(shù)頻率響應(yīng)為：當參數(shù)截止頻率改為：fpass=4khz fstop=6khz 此時計算得n=203濾波器的對數(shù)頻率響應(yīng)為：結(jié)果分析從以上基于kaiser窗的低通fir濾波器的結(jié)果可以得知：相對于矩形窗和hamming窗，kaiser窗的階數(shù)及形狀可以實現(xiàn)控制，而前者兩個窗其、apass、astop和d是一定的，實現(xiàn)不了了任意需求的濾波器，而凱瑟窗的各個參數(shù)是是可以變化的，可以實現(xiàn)任意技術(shù)指標需求的濾波器。kaiser窗是性能取決于窗的階數(shù)n（窗口長度）和其形狀控制參數(shù)，而階數(shù)n和形狀控制參數(shù)可以由各個技術(shù)指標計算得出，從而實現(xiàn)了滿足任意需求的濾波器。通過改變通帶截止頻率和阻帶截止頻率，可以實現(xiàn)不改變kaiser窗的形狀控制參數(shù)而改濾波器的階數(shù)n，上面分別做了fpass=4khz