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

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

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

3、題探究】在譜分析中如何選擇窗函數(shù)，在濾波器設(shè)計(jì)中如何選擇窗函數(shù)？在譜分析中窗函數(shù)的選擇需要考慮泄露現(xiàn)象，當(dāng)需要分辨的頻率為f時，若使用矩形窗，其寬度必須大于fsamf。當(dāng)要求分辨較小的頻率而計(jì)算量較小可以選擇其他窗函數(shù)來改善。在FIR濾波器的設(shè)計(jì)中，使用矩形窗會由于其銳截止而產(chǎn)生Gibbs現(xiàn)象，會有9%的過沖，為了改善效果，可以在滿足過渡帶寬度和增益要求的前提下適當(dāng)選擇其他種類的窗函數(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è)計(jì)數(shù)字濾波器從和弦中分離出369.99Hz的音符。要求：（1）設(shè)計(jì)IIR數(shù)字帶通濾波器，通過實(shí)驗(yàn)研究不同、過渡帶、對濾波器設(shè)計(jì)的影響，確定本題最合適的濾波器指標(biāo)。（2）通過將IIR數(shù)字低通濾波器和數(shù)字高通濾波器級聯(lián)也可實(shí)現(xiàn)369.99Hz音符的分離。分別確定 IIR數(shù)字低通濾波器和數(shù)字高通濾波器的指標(biāo)，并進(jìn)行濾

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

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

12、=20(此兩個指標(biāo)暫不作考慮）【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兩個指標(biāo)為變量研究其對低通濾波器的影響一、研究As對濾波器的影響 As³10dB As³20dB As³32dB As³40dB由上面四幅圖可知，一味增

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

14、0.085p, (3)有表格可以看出，選取不同的窗函數(shù)，要通過確定其過渡帶寬度來確定其階數(shù)，因此，在相同過渡帶寬度的前提下，我們通過不同的窗函數(shù)設(shè)計(jì)了不同的帶通濾波器，其幅度相應(yīng)及階數(shù)如下圖所示。矩形窗 N=301漢納窗 N=1034哈明窗 N=1167布萊克曼窗 N=1901凱澤窗 N=237(beta=3.395)(4)(1)中帶通濾波器的指標(biāo)如下：Wp1=0.085p, Wp2=0.099p, Ws1=0.078p, Ws2=0.105p Ap£1dB As³20dB 轉(zhuǎn)換為FIR數(shù)字帶通濾波器指標(biāo)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)的方法也能達(dá)到濾波的效果，同時成本較帶通而言基本一致，但不易將雜波完全濾去，仍會保留一些高頻、低頻分量。從頻譜來看，濾波效果不如帶通濾波器；從聲音效果來看，兩者幾乎一樣。高低通級聯(lián)濾波效果不如用帶通濾波器好，但在物理上更易實(shí)現(xiàn)；在誤差允許的范圍內(nèi)，高低通級聯(lián)也可以實(shí)現(xiàn)帶通的效果。因此，工程上可以考慮使用高低通級聯(lián)代替帶通濾波器。(3)用窗函數(shù)發(fā)設(shè)計(jì)FIR數(shù)字帶通濾波器時，首先確定過渡帶寬度，然后要選擇窗函數(shù)，根據(jù)窗函數(shù)近似過渡帶

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

18、也很難實(shí)現(xiàn)。所以當(dāng)發(fā)現(xiàn)這種情況時，我們應(yīng)該改變設(shè)計(jì)指標(biāo)，使濾波器更易實(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(&#