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

1、數(shù)字信號(hào)處理課程研究性學(xué)習(xí)報(bào)告姓名 初步草稿 學(xué)號(hào) 同組成員 指導(dǎo)教師 時(shí)間 數(shù)字濾波器設(shè)計(jì)專(zhuān)題研討【目的】(1) 掌握IIR和FIR數(shù)字濾波器的設(shè)計(jì)方法及各自的特點(diǎn)。(2) 掌握各種窗函數(shù)的時(shí)頻特性及對(duì)濾波器設(shè)計(jì)的影響。(3) 培養(yǎng)學(xué)生自主學(xué)習(xí)能力，以及發(fā)現(xiàn)問(wèn)題、分析問(wèn)題和解決問(wèn)題的能力。【研討題目】中等題 Dhexian.wav是對(duì)頻率為293.66, 369.99, 440Hz的D大調(diào)和弦以8000Hz抽樣所得的數(shù)字音樂(lè)信號(hào)，試設(shè)計(jì)一數(shù)字濾波器從和弦中分離出440Hz的音符。要求：（1）設(shè)計(jì)IIR數(shù)字高通濾波器，通過(guò)實(shí)驗(yàn)研究不同、過(guò)渡帶、對(duì)濾波器設(shè)計(jì)的影響，確定本題最合適的濾波器指標(biāo)。（

2、2）能否用IIR數(shù)字帶通濾波器從和弦中分離出440Hz的音符？試參照（1）確定的最合適的高通濾波器指標(biāo)，給出數(shù)字帶通濾波器的指標(biāo)。設(shè)計(jì)IIR數(shù)字帶通濾波器，并將結(jié)果與高通濾波器比較，給出你的結(jié)論。（3）用窗函數(shù)法設(shè)計(jì)FIR數(shù)字高通濾波器，分別利用矩形窗、漢納窗、哈明窗、布萊克曼窗、凱澤窗截?cái)?。討論用窗函?shù)法設(shè)計(jì)FIR數(shù)字高通濾波器時(shí)如何確定濾波器的指標(biāo)，比較相同過(guò)渡帶時(shí)用矩形窗、漢納窗、哈明窗、布萊克曼窗、凱澤窗設(shè)計(jì)濾波器的階數(shù)。（4）采用Parks-McClellan算法，設(shè)計(jì)FIR數(shù)字高通濾波器。試參照（1）確定的最合適的高通濾波器指標(biāo)，給出FIR數(shù)字高通濾波器的指標(biāo)。將設(shè)計(jì)結(jié)果與（1）中

3、的IIR數(shù)字濾波器，從幅度響應(yīng)、相位響應(yīng)、濾波器階數(shù)等方面進(jìn)行比較。【溫磬提示】在IIR數(shù)字濾波器的設(shè)計(jì)中，不管是用雙線(xiàn)性變換法還是沖激響應(yīng)不變法，其中的參數(shù)T的取值對(duì)設(shè)計(jì)結(jié)果沒(méi)有影響。但若所設(shè)計(jì)的數(shù)字濾波器要取代指定的模擬濾波器時(shí)，則抽樣頻率（或抽樣間隔T）將對(duì)設(shè)計(jì)結(jié)果有影響?！驹O(shè)計(jì)步驟】:首先以8000Hz抽樣時(shí)，293.66, 369.99, 440Hz對(duì)應(yīng)的數(shù)字頻率分別為0.07341*pi,0.0924975*pi,0.11*pi。此處在設(shè)計(jì)低通濾波器時(shí)，Wp顯然要小于0.11*pi，而Ws要大于0.0924975*pi。另外此處設(shè)計(jì)的IIR濾波器均為BW型濾波器，采用雙線(xiàn)性變換法。

4、【仿真結(jié)果】：考慮到一般濾波器對(duì)通阻帶截頻的要求，我們可以先讓As=0.5,Ap=15(此兩個(gè)指標(biāo)暫不作考慮）【結(jié)果分析】 1.Wp=0.11*pi,Ws=0.092497*pi,Ap=0.5,As=15 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看評(píng)價(jià)：聽(tīng)聲音感覺(jué)濾波前后有一些變化，嘈雜、渾濁的聲音變得較為清脆，顯然是頻譜成分減少了。仔細(xì)觀察濾波前后音樂(lè)的頻譜，440Hz的信號(hào)保留了，但300400Hz間也有部分信號(hào)沒(méi)有被濾掉，盡管其幅度較小。猜測(cè)是阻帶截頻太小導(dǎo)致此情況，所以增大阻帶截頻。 改變Ws2. Wp= 0.11*pi Ws=0.0944975 Ap=0.5,As=15 音樂(lè)濾波后的頻譜關(guān)鍵位置

5、放大看評(píng)價(jià)：聽(tīng)聲音感覺(jué)和上次一樣，濾波前后有一些變化，嘈雜、渾濁的聲音變得較為清脆，顯然是頻譜成分減少了。仔細(xì)觀察濾波前后音樂(lè)的頻譜，440Hz的信號(hào)保留了，但300400Hz間也有部分信號(hào)沒(méi)有被濾掉，但其幅度較上次有所減小。3. Wp= 0.11*pi Ws=0.0954975 Ap=0.5,As=15音樂(lè)濾波后的頻譜關(guān)鍵位置放大看評(píng)價(jià)：聽(tīng)聲音感覺(jué)濾波前后有一些變化，嘈雜、渾濁的聲音變得較為清脆，顯然是頻譜成分減少了。仔細(xì)觀察濾波前后音樂(lè)的頻譜，440Hz的信號(hào)保留了，但300400Hz間的信號(hào)基本被濾掉，剩下的部分幅度很小，結(jié)果較為滿(mǎn)意。4. Wp= 0.11*pi Ws=0.096497

6、5 Ap=0.5,As=15音樂(lè)濾波后的頻譜關(guān)鍵位置放大看評(píng)價(jià)：聽(tīng)聲音感覺(jué)不如前幾次效果好，無(wú)論從濾波器的幅度響應(yīng)還是濾波后信號(hào)的頻譜看，效果均較差。綜上，確定Ws=0.0954975pi 改變Wp5. Wp= 0.1099*pi Ws=0.0954975 *pi Ap=0.5,As=15 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看評(píng)價(jià)：聽(tīng)聲音感覺(jué)不如前幾次效果好，無(wú)論從濾波器的幅度響應(yīng)還是濾波后信號(hào)的頻譜看，效果均較差。Wp減小一點(diǎn)就會(huì)造成較大誤差。 再次總結(jié)：取Wp=0.11*pi，Ws=0.0954975*pi 探究Ap、As的影響6. Wp= 0.11*pi Ws=0.0954975*pi Ap=

7、0.5,As=13音樂(lè)濾波后的頻譜關(guān)鍵位置放大看 關(guān)鍵位置放大看7. Wp= 0.11*pi Ws=0.0954975*pi Ap=0.5,As=14 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看8. Wp= 0.11*pi Ws=0.0954975*pi Ap=0.5,As=15 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看 9.Wp= 0.11*pi Ws=0.0954975*pi Ap=0.5,As=16 關(guān)鍵位置放大看再次總結(jié)：從聲音的變化、濾波器的幅度響應(yīng)和濾波后的頻譜除440Hz其他成分的幅度三方面來(lái)看，確定As=1510.Wp= 0.11*pi Ws=0.0954975*pi Ap=0.3,As=1511.

8、 Wp= 0.11*pi Ws=0.0954975*pi Ap=1.0,As=15 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看12. Wp= 0.11*pi Ws=0.0954975*pi Ap=1.5,As=15 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看綜上：確定高通濾波器的指標(biāo)為Wp= 0.11*pi（440Hz） Ws=0.0954975*pi （381.99Hz） Ap=0.5,As=15(2) 從理論上分析，只要440Hz在數(shù)字帶通濾波器的通帶范圍內(nèi)，293.66Hz,369.99Hz處在其阻帶內(nèi)，就能從和弦中分離出440Hz的音符。根據(jù)以上設(shè)計(jì)的高通濾波器的指標(biāo)先確定帶通濾波器的以下指標(biāo)，W1.11*p

9、i,Ws1=0.0954975*pi,As=15,Ap=0.5,至于Wp2和Ws2的確定如下：1. Wp2=0.190*pi,Ws2=0.999*pi實(shí)驗(yàn)程序：fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音樂(lè)頻譜');%設(shè)計(jì)數(shù)字濾波器Ws=0.0954975;W

10、p=0.11;Ap=0.5;As=15;N,Wc=buttord(Wp,Ws,Ap,As)num,den=butter(N,Wc,'high');%設(shè)計(jì)數(shù)字濾波器w=linspace(0,pi,1024);h=freqz(num,den,w);%求濾波器的增益響應(yīng)figure(2);plot(w/pi,20*log10(abs(h);gridtitle('高通數(shù)字濾波器的增益響應(yīng)');y=filter(num,den,x);%濾波操作wavplay(y,fs);%聽(tīng)濾波后的音樂(lè)Y=fftshift(fft(y,L);figure(3);f=(0:L-1)*fs/

11、L-0.5*fs;plot(f,abs(Y);title('音樂(lè)濾波后的頻譜'); 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看2. Wp2= 0.19*pi,Ws2=0.8*pi音樂(lè)濾波后的頻譜關(guān)鍵位置放大看3. Wp2=0.19 *pi,Ws2=0.5*pi 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看4. Wp2= 0.4*pi,Ws2=0.5*pi5.Wp2=0.4 *pi,Ws2=0.7*pi6. Wp2=0.89 *pi,Ws2=0.9075025*pi總結(jié)分析：由以上實(shí)驗(yàn)結(jié)果可知不同的Wp2,Ws2指標(biāo)對(duì)實(shí)驗(yàn)結(jié)果的影響很小，究其原因是所給的音樂(lè)的頻譜較為單一，它只在293.66Hz,369.

12、99Hz和440Hz及其附近幅度較大，在大于440Hz處幅度基本為零，所以帶通濾波器的Wp2,Ws2指標(biāo)對(duì)實(shí)驗(yàn)結(jié)果的影響非常小。帶通濾波器確實(shí)能從和弦中分離出440Hz的音符，在此處從濾波器指標(biāo)的對(duì)稱(chēng)性考慮，我們選擇帶通濾波器的指標(biāo)為Ws1=0.0954975*pi(381.99Hz);Wp1=0.11*pi(440Hz);Wp2=0.89*pi(3560Hz);Ws2=0.9045025*pi(3619Hz);As=15;Ap=15。實(shí)驗(yàn)程序：fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian

13、.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音樂(lè)頻譜');%設(shè)計(jì)數(shù)字濾波器Ws1=0.0954975;Wp1=0.11;Wp2=0.89;Ws2=0.9075025;Wp=Wp1,Wp2;Ws=Ws1,Ws2;Ap=0.5;As=15;N,Wc=buttord(Wp,Ws,Ap,As)num,den=butter(N,Wc);%設(shè)計(jì)數(shù)字濾波器w=linspace(0,pi,1024);h=freqz(num,

14、den,w);%求濾波器的增益響應(yīng)figure(2);plot(w/pi,20*log10(abs(h);gridtitle('高通數(shù)字濾波器的增益響應(yīng)');y=filter(num,den,x);%濾波操作wavplay(y,fs);%聽(tīng)濾波后的音樂(lè)Y=fftshift(fft(y,L);figure(3);f=(0:L-1)*fs/L-0.5*fs;plot(f,abs(Y);title('音樂(lè)濾波后的頻譜');（3）1. 用矩形窗（階數(shù)M=124) 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看%用矩形窗fs=8000;nbits=16;L=8192;x,fs,nbits

15、=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音樂(lè)頻譜');N=ceil(1.8*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1w=linspace(1,N,1);Wc=(Wp+Ws)/2;k=0:M;hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)

16、+1;h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);figure(2);plot(omega/pi,magdb);title('用矩形窗設(shè)計(jì)的高通濾波器的幅度響應(yīng)');y=filter(h,1,x);wavplay(y,fs);%聽(tīng)濾波后音樂(lè)的頻譜f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);title('濾波后音樂(lè)的頻譜');2. 用Hanning窗(階數(shù)M=428）

17、音樂(lè)濾波后的頻譜關(guān)鍵位置放大看漢寧窗法的程序：%用hanning窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音樂(lè)頻譜');Ws=0.0954972*pi;Wp=0.11*pi;Wc=(Wp+Ws)/(2*pi);%使濾波器為I型N=ceil(6.2*pi/(W

18、p-Ws);N=mod(N+1,2)+N;M=N-1h = fir1(M,Wc,'high',hanning(N);omega=linspace(0,pi,512);mag=freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用hanning窗設(shè)計(jì)的高通濾波器的幅度響應(yīng)');y=filter(h,1,x);wavplay(y,fs);%聽(tīng)濾波后音樂(lè)的頻譜f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);

19、title('濾波后音樂(lè)的頻譜');3. 用哈明窗（階數(shù)M=482) 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看哈明窗法的程序：%用Hamming窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音樂(lè)頻譜');Ws=0.0954972*pi;Wp=0.11*pi

20、;Wc=(Wp+Ws)/(2*pi);%使濾波器為I型N=ceil(7*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1h = fir1(M,Wc,'high',hamming(N);omega=linspace(0,pi,512);mag=freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用Hamming窗設(shè)計(jì)的高通濾波器的幅度響應(yīng)');y=filter(h,1,x);wavplay(y,fs);%聽(tīng)濾波后音樂(lè)的頻譜f=(0:L-1)*fs/L-0.5*fs;Y=f

21、ftshift(fft(y,L);figure(3);plot(f,abs(Y);title('濾波后音樂(lè)的頻譜');4. 布萊克曼窗(階數(shù)M=786） 音樂(lè)濾波后的頻譜關(guān)鍵位置放大看Blackman窗法的程序：%用blackman窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);t

22、itle('原始音樂(lè)頻譜');Ws=0.0954972*pi;Wp=0.11*pi;Wc=(Wp+Ws)/(2*pi);%使濾波器為I型N=ceil(11.4*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1h = fir1(M,Wc,'high',blackman(N);omega=linspace(0,pi,512);mag=freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用blackman窗設(shè)計(jì)的高通濾波器的幅度響應(yīng)');y=filter(h,

23、1,x);wavplay(y,fs);%聽(tīng)濾波后音樂(lè)的頻譜f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);title('濾波后音樂(lè)的頻譜');5. 凱澤窗（階數(shù)M=166）音樂(lè)濾波后的頻譜關(guān)鍵位置放大看采用凱澤窗的程序：%用kaiser窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%聽(tīng)濾波前的音樂(lè)f=(0:L-1)*fs/L-0.5*fs;X

24、=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音樂(lè)頻譜');%Program:利用Kaiser窗設(shè)計(jì)FIR高通濾波器Ap=0.5;As=15;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*As);f=0.0954975,0.11;a=0,1;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f,a,dev);%使濾波器為I型M=mod(M,2)+Mh = fir1(M,Wc,ftype,kaiser(M+1,beta);omega=linspace(0,pi,512);mag=

25、freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用kaiser窗設(shè)計(jì)的高通濾波器的幅度響應(yīng)');y=filter(h,1,x);wavplay(y,fs);%聽(tīng)濾波后音樂(lè)的頻譜f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);title('濾波后音樂(lè)的頻譜');【結(jié)果分析】：一、用窗函數(shù)法設(shè)計(jì)FIR濾波器的方法及步驟：1.根據(jù)Wp,Ws，Ap,As確定逼近的理想濾波器及窗函數(shù)類(lèi)型及窗的長(zhǎng)度。由Wp,

26、Ws確定理想濾波器Wc及窗長(zhǎng)度N,由As確定窗類(lèi)型，2. 根據(jù)要求確定線(xiàn)性相位FIR濾波器的類(lèi)型3. 確定理想濾波器的幅度函數(shù)Ad()和相位d()，d()= -0.5M+ 4.計(jì)算IDTFT得hdk 5.截短hdk, hk= hdkwk， 0<=k<= M二、用MATLAB實(shí)現(xiàn)的方法： h = fir1(M,Wc,'ftype',window)M表示濾波器的階數(shù)；Wc表示理想FIR濾波器的B個(gè)頻帶； ftype表示濾波器的類(lèi)型，缺省值為空；window是一長(zhǎng)度為N+1的向量, Hamming窗3、 本題中相同過(guò)渡帶時(shí)，矩形窗、漢納窗、哈明窗、布萊克曼窗、凱澤窗的階數(shù)

27、分別是 矩形窗（124）<凱澤窗（166）<漢納窗（428）<哈明窗（482）<布萊克曼窗（786）(4) 采用Parks-McClellan算法 濾波后音樂(lè)頻譜關(guān)鍵位置放大看【結(jié)果分析】：比較采用Parks-McClellan算法和IIR數(shù)字濾波器。從幅度響應(yīng)來(lái)看，前者在通帶和阻帶都會(huì)波動(dòng)，而后者只在阻帶波動(dòng)，通帶不波動(dòng)從相位響應(yīng)來(lái)看，前者在0到pi間波動(dòng)，且較均勻，而后者也在0到pi間波動(dòng)，但是在低頻附近波動(dòng)更厲害，隨著頻率升高，波動(dòng)減弱。從階數(shù)來(lái)看，前者為102階，后者為20階，IIR數(shù)字濾波器的階數(shù)明顯低于采用Parks-McClellan算法設(shè)計(jì)的濾波器的階數(shù)。%采用Parks-McClellan算法fs=8000;nbits=16;L=8192;x,fs,nbits=wavread(&#