北交大DSP研究性學(xué)習(xí)報告數(shù)字濾波器設(shè)計2014

1、數(shù)字信號處理課程研究性學(xué)習(xí)報告姓名 學(xué)號 同組成員 指導(dǎo)教師 時間 2014年5月18日星期日 數(shù)字濾波器設(shè)計專題研討【目的】(1) 掌握IIR和FIR數(shù)字濾波器的設(shè)計方法及各自的特點。(2) 掌握各種窗函數(shù)的時頻特性及對濾波器設(shè)計的影響。(3) 培養(yǎng)學(xué)生自主學(xué)習(xí)能力，以及發(fā)現(xiàn)問題、分析問題和解決問題的能力。【研討題目】 基本題 分析矩形窗、漢納窗、哈明窗、布萊克曼窗、凱澤窗的頻域特性，并進行比較?！绢}目分析】本題分析不同的窗函數(shù)的頻域特性，預(yù)計可以看出不同的窗有不同的過渡帶大小和不同的旁瓣寬度，可以滿足對不同的設(shè)計設(shè)計要求。【仿真結(jié)果】【結(jié)果分析】各種窗有何特點?計算過程中，不同的窗都采用了

2、相同的長度（均為10），進行fft計算的長度均為512點，從以上結(jié)果可以看出，矩形窗的幅度最大，主瓣最窄，旁瓣幅度大；其余幾個窗函數(shù)的主瓣幅度均小于矩形窗的主瓣幅度，主瓣寬度約為矩形窗主瓣寬度的2倍多；其中漢納窗和哈明窗的旁瓣有較小幅度，而布萊克曼窗和凱澤窗旁瓣基本為0。【閱讀文獻】1數(shù)字信號處理.陳后金【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：matlab提供的函數(shù)產(chǎn)生的窗為N行1列的列向量，不是1行N列的一行數(shù)據(jù)，計算時必須將其進行轉(zhuǎn)置。使用fft計算的結(jié)果是原序列周期化后求出的頻譜的主值區(qū)間，必須使用fftshift將其整理，然后通過向左平移L/2才能求出原序列的頻譜?！締?/p>

3、題探究】在譜分析中如何選擇窗函數(shù)，在濾波器設(shè)計中如何選擇窗函數(shù)？在譜分析中窗函數(shù)的選擇需要考慮泄露現(xiàn)象，當需要分辨的頻率為f時，若使用矩形窗，其寬度必須大于fsamf。當要求分辨較小的頻率而計算量較小可以選擇其他窗函數(shù)來改善。在FIR濾波器的設(shè)計中，使用矩形窗會由于其銳截止而產(chǎn)生Gibbs現(xiàn)象，會有9%的過沖，為了改善效果，可以在滿足過渡帶寬度和增益要求的前提下適當選擇其他種類的窗函數(shù)?！痉抡娉绦颉?·ÖÎö¸÷ÖÖ´°º¯ÊýµÄ&#

4、198;µÓòÌØÐÔclear all;L=512;%L±ØÐè´óÓÚÏßÐÔ¾í»ýµÄ³¤¶È£¬2µÄÕýÊý´ÎÃÝN=10;%ÒÔÏÂ&

5、#184;÷¸ö´°µÄ³¤¶È%¾ØÐδ°figure(1)w1=zeros(1,50);w2=ones(1,N);w3=zeros(1,50);wh1=w1 w2 w3;WH1=fftshift(fft(wh1,L);w=(0:L-1)-L/2;plot(w,abs(WH1);title('¾ØÐδ°µÄ·ù&#

6、182;ÈÆ×');%ººÄÉ´°figure(2)wh2=hann(N)'WH2=fftshift(fft(wh2,L);w=(0:L-1)-L/2;plot(w,abs(WH2);title('ººÄÉ´°µÄ·ù¶ÈÆ×');%¹þÃ÷´°figure(3)wh3=ham

7、ming(N)'WH3=fftshift(fft(wh3,L);w=(0:L-1)-L/2;plot(w,abs(WH3);title('¹þÃ÷´°µÄ·ù¶ÈÆ×');%²¼À³¿ËÂü´°figure(4)wh4=blackman(N)'WH4=fftshift(fft(wh4,L);w=(0:L-1)-L/2;plot

8、(w,abs(WH4);title('²¼À³¿ËÂü´°µÄ·ù¶ÈÆ×');%¿­Ôó´°figure(5)beta=6;wh5=kaiser(N,beta)'WH5=fftshift(fft(wh5,L);w=(0:L-1)-L/2;plot(w,abs(WH5);title('¿­Ô

9、43;´°µÄ·ù¶ÈÆ×');【研討題目】中等題 Dhexian.wav是對頻率為293.66, 369.99, 440Hz的D大調(diào)和弦以8000Hz抽樣所得的數(shù)字音樂信號，試設(shè)計數(shù)字濾波器從和弦中分離出369.99Hz的音符。要求：（1）設(shè)計IIR數(shù)字帶通濾波器，通過實驗研究不同、過渡帶、對濾波器設(shè)計的影響，確定本題最合適的濾波器指標。（2）通過將IIR數(shù)字低通濾波器和數(shù)字高通濾波器級聯(lián)也可實現(xiàn)369.99Hz音符的分離。分別確定 IIR數(shù)字低通濾波器和數(shù)字高通濾波器的指標，并進行濾

10、波器設(shè)計。將結(jié)果與帶通濾波器比較，給出你的結(jié)論。（3）用窗函數(shù)法設(shè)計FIR數(shù)字帶通濾波器，分別利用矩形窗、漢納窗、哈明窗、布萊克曼窗、凱澤窗截斷。討論用窗函數(shù)法設(shè)計FIR數(shù)字帶通濾波器時如何確定濾波器的指標，比較相同過渡帶時用矩形窗、漢納窗、哈明窗、布萊克曼窗、凱澤窗設(shè)計濾波器的階數(shù)。（4）采用Parks-McClellan算法，設(shè)計FIR數(shù)字帶通濾波器。試參照（1）確定的最合適的帶通濾波器指標，給出FIR數(shù)字帶通濾波器的指標。將設(shè)計結(jié)果與（1）中的IIR數(shù)字濾波器，從幅度響應(yīng)、相位響應(yīng)、濾波器階數(shù)等方面進行比較?！緶仨嗵崾尽吭贗IR數(shù)字濾波器的設(shè)計中，不管是用雙線性變換法還是沖激響應(yīng)不變法，

11、其中的參數(shù)T的取值對設(shè)計結(jié)果沒有影響。但若所設(shè)計的數(shù)字濾波器要取代指定的模擬濾波器時，則抽樣頻率（或抽樣間隔T）將對設(shè)計結(jié)果有影響。【設(shè)計步驟】1.講所給模擬頻率轉(zhuǎn)換為數(shù)字指標：通過公式=2ffsam可以計算出三個音符對應(yīng)的數(shù)字頻率1=0.073 2=0.092 3=0.11 2.采用雙線性變換法，通過公式=2Ttan(2)將數(shù)字頻率指標轉(zhuǎn)換成模擬指標3.根據(jù)模擬指標設(shè)計模擬濾波器4.利用公式Hz=Hs|s=2T1-z-11+z-1將所得模擬濾波器轉(zhuǎn)換為數(shù)字濾波器【仿真結(jié)果】(1) 確定數(shù)字帶通濾波器的參數(shù)：首先討論通阻帶截頻參數(shù)，考慮到一般濾波器對通阻帶截頻的要求，我們可以先讓As=1,Ap

12、=20(此兩個指標暫不作考慮）【1】取p1=0.089 ,p2=0.095 ,s1=0.073, s2=0.11 從濾波后頻譜可以看出，仍留存有高低頻分量，因此應(yīng)調(diào)整阻帶頻率【2】取p1=0.089 ,p2=0.095 ,s1=0.080, s2=0.10 濾波效果較好，經(jīng)多次調(diào)整后，最終確定的濾波器參數(shù)為p1=0.085 ,p2=0.099 ,s1=0.078, s2=0.105。 下面討論Ap, As兩個指標為變量研究其對低通濾波器的影響一、研究As對濾波器的影響 As³10dB As³20dB As³32dB As³40dB由上面四幅圖可知，一味增

13、大阻帶衰減幅度可能是結(jié)果更糟，只需選取適當阻帶衰減幅度即可，故令A(yù)s³20dB 二、研究Ap對濾波器的影響 Ap£0.5dB Ap£1dBAp£2dB由圖可知Ap£0.5dB時阻帶衰減幅度最大，但Ap£1dB時也可以滿足濾波的要求，考慮到系統(tǒng)成本，選取Ap£1dB綜上所述，帶通濾波器的指標確定為：p1=0.085 ,p2=0.099 ,s1=0.078, s2=0.105,Ap£1dB, As³20dB (2) 指標： Wp2=0.095p, Ws2=0.1p 指標：Wp1=0.0875p, , Ws1=

14、0.085p, (3)有表格可以看出，選取不同的窗函數(shù)，要通過確定其過渡帶寬度來確定其階數(shù)，因此，在相同過渡帶寬度的前提下，我們通過不同的窗函數(shù)設(shè)計了不同的帶通濾波器，其幅度相應(yīng)及階數(shù)如下圖所示。矩形窗 N=301漢納窗 N=1034哈明窗 N=1167布萊克曼窗 N=1901凱澤窗 N=237(beta=3.395)(4)(1)中帶通濾波器的指標如下：Wp1=0.085p, Wp2=0.099p, Ws1=0.078p, Ws2=0.105p Ap£1dB As³20dB 轉(zhuǎn)換為FIR數(shù)字帶通濾波器指標Wp1=0.085p, Wp2=0.099p, Ws1=0.078p,

15、 Ws2=0.105p dp = 0.109 ds = 0.1 等紋濾波器幅度譜和相位譜 N=545IIR帶通濾波器 N=6【結(jié)果分析】(1)(2)采用高低通濾波器級聯(lián)的方法也能達到濾波的效果，同時成本較帶通而言基本一致，但不易將雜波完全濾去，仍會保留一些高頻、低頻分量。從頻譜來看，濾波效果不如帶通濾波器；從聲音效果來看，兩者幾乎一樣。高低通級聯(lián)濾波效果不如用帶通濾波器好，但在物理上更易實現(xiàn)；在誤差允許的范圍內(nèi)，高低通級聯(lián)也可以實現(xiàn)帶通的效果。因此，工程上可以考慮使用高低通級聯(lián)代替帶通濾波器。(3)用窗函數(shù)發(fā)設(shè)計FIR數(shù)字帶通濾波器時，首先確定過渡帶寬度，然后要選擇窗函數(shù)，根據(jù)窗函數(shù)近似過渡帶

16、的寬度確定濾波器的階數(shù)N。由實驗結(jié)果可知，按照矩形窗、漢納窗、哈明窗、布萊克曼窗的順序，阻帶衰減逐漸增大，階數(shù)逐漸增高。而選擇凱澤窗時，既能保證較高的阻帶衰減，又有較小的階數(shù)，與其他窗函數(shù)相比，更為靈活。(4)比較采用Parks-McClellan算法和IIR數(shù)字濾波器。  幅度響應(yīng)來看，前者在通帶和阻帶都會波動，而后者波紋明顯較小，但過渡帶較前者寬； 從相位響應(yīng)來看，前者在0到pi間波動，且較均勻，而后者也在0到pi間波動，但是在低頻附近波動更厲害，隨著頻率升高，波動減弱。  從階數(shù)來看，前者為545階，后者為6階，IIR數(shù)字濾波器的階數(shù)明顯低于采用Parks-McCle

17、llan算法設(shè)計的濾波器的階數(shù)?！景l(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：設(shè)計濾波器的過程中使用BW型有時頻譜會有振蕩產(chǎn)生,并且當阻帶增益As越大時，震蕩越明顯，甚至?xí)霈F(xiàn)嚴重的失真?！締栴}探究】過渡帶越窄，通帶最大衰減越小，阻帶最小衰減越大，則濾波器設(shè)計難度越大；反之，過渡帶越寬，通帶最大衰減越大，阻帶最小衰減越小，濾波器更易設(shè)計。我們發(fā)現(xiàn)，當濾波器設(shè)計難度到達一定程度時，matlab已經(jīng)不能設(shè)計出相應(yīng)濾波器或是出現(xiàn)上圖震蕩的情況。經(jīng)查閱資料，發(fā)現(xiàn)應(yīng)該是濾波器指標不合適導(dǎo)致達到計算機計算的極限，結(jié)果會有相當大的誤差。原本理論可實現(xiàn)的系統(tǒng)計算機無法算出。但這也說明這樣的系統(tǒng)物理上

18、也很難實現(xiàn)。所以當發(fā)現(xiàn)這種情況時，我們應(yīng)該改變設(shè)計指標，使濾波器更易實現(xiàn)?！痉抡娉绦颉?1)Fsam=8000;x1 = audioread('Dhexian.wav');y1,f1=ctft1(x1,Fsam,8001);t=0:8000;figure(1)plot(f1,abs(y1)axis(200,500,0,0.18)sound(x1,Fsam);pause(2); Ws1=0.078*pi; Wp1=0.085*pi; Wp2=0.099*pi ; Ws2=0.105*pi; Ap=1;As=20;%Êý×Öµ

19、05;ͨÂ˲¨Æ÷Ö¸±êws1=tan(Ws1/2); wp1=tan(Wp1/2); wp2=tan(Wp2/2); ws2=tan(Ws2/2);%Âö³åÏìÓ¦²»±ä·¨×ª»¯³ÉÄ£Äâ´øÍ

20、;¨Â˲¨Æ÷Ö¸±êB=wp2-wp1; w0=sqrt(wp1*wp2); wp=1; ws=min(abs(ws12-wp1*wp2)/(B*ws1),abs(ws22-wp1*wp2)/(B*ws2);%Ä£Äâ´øͨÂ˲¨Æ÷Ö¸±ê»¯ÎªÄ

21、£ÄâÔ­Ð͵ÍͨÂ˲¨Æ÷Ö¸±êN,Wc=buttord(wp,ws,Ap,As,'s');num,den = butter(N,Wc,'s');numt,dent = lp2bp(num,den,w0,B);numd,dend=bilinear(numt,dent,0.5);w=linspace(0,pi,1024);h=freqz(numd,

22、dend,w);figure(2)plot(w/pi,20*log10(abs(h),'r');grid;xlabel('Normalized frequency');ylabel('Gain,dB');axis(0,0.2,-100,0)%hold on; figure(3)x2=filter(numd,dend,x1);y2,f2=ctft1(x2,Fsam,8001);plot(f2,abs(y2);axis(-500,500,0,0.18)player1=audioplayer(x2,fs);play(player1);axis(200,

23、500,0,0.18) (2)fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('E:¿Î³ÌÊý×ÖÐźŴ¦ÀíDSP2Dhexian.wav');wavplay(x,fs);f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('ԭʼ&#

24、210;ôÀÖƵÆ×');%hpWp1=365*2/fs;Ws1=330*2/fs;Wp2=3560*2/fs;Ws2=3605*2/fsWp=Wp1,Wp2;Ws=Ws1,Ws2;Ap=1;As=12;N,wc=buttord(Wp,Ws,Ap,As);num,den=butter(N,wc);w=linspace(0,pi,1000);h1=freqz(num,den,w);figure(2);plot(w/pi,20*log10(abs(h1);grid on;title('¸ß&

25、#205;¨Êý×ÖÂ˲¨Æ÷µÄÔöÒæÏìÓ¦');axis(0 0.8 -180 20);y1=filter(num,den,x);%Â˲¨²Ù×÷?wavplay(y1,fs);%ÌýÂ˲¨ºóµ&

26、#196;ÒôÀÖ?Y1=fftshift(fft(y1,L);figure(3);f=(0:L-1)*fs/L-0.5*fs;plot(f,abs(Y1);title('ÒôÀÖÂ˲¨ºóµÄƵÆ×');%lp Wp4=370*2/fs;Ws4=400*2/fsAp=1;As=7;N,wc=buttord(Wp4,Ws4,Ap,As);num,den=butter(N,wc)

27、;w=linspace(0,pi,1000);h2=freqz(num,den,w);figure(4);plot(w/pi,20*log10(abs(h2);grid on;title('µÍͨÊý×ÖÂ˲¨Æ÷µÄÔöÒæÏìÓ¦');y2=filter(num,den,y1);%Â˲¨

28、78;Ù×÷?wavplay(y2,fs);%ÌýÂ˲¨ºóµÄÒôÀÖ?Y2=fftshift(fft(y2,L);figure(5);f=(0:L-1)*fs/L-0.5*fs;plot(f,abs(Y2);title('ÒôÀÖÂ˲¨ºóµÄƵÆ×&

29、#39;);h=h1.*h2;%plot(w/pi,20*log10(abs(h);grid on;(3)Ws1=0.078*pi; Wp1=0.085*pi; Wp2=0.099*pi ; Ws2=0.105*pi; Ap=1;As=20;Wc1=(Wp1+Ws1)/2;Wc2=(Wp2+Ws2)/2;%¾ØÐδ°% N1=ceil(1.8*pi/(Wp1-Ws1);% N2=ceil(1.8*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2,boxcar(N+1);%

30、6;ºÄÉ´°% N1=ceil(6.2*pi/(Wp1-Ws1);% N2=ceil(6.2*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2,hann(N+1);%¹þÃ÷´°% N1=ceil(7*pi/(Wp1-Ws1);% N2=ceil(7*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2);%²¼À³¿Ë&

31、#194;ü´°% N1=ceil(11.4*pi/(Wp1-Ws1);% N2=ceil(11.4*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2,blackman(N+1);%¿­Ôó´°f=Ws1 Wp1 Wp2 Ws2;a=0,1,0;Rs=0.01;dev=Rs*ones(1,length(a);N,Wc,beta,ftype = kaiserord(f,a,dev); numd = fir1(N,Wc,ftype,kaiser(N+1,beta);w=linspace(0,pi,1024);h=freqz(numd,1,w);plot(w/(pi*pi),20*log10(abs(h),'r');grid;xlabel('Normalized frequency');ylabel(&#