《數(shù)字信號處理》課程研究性學(xué)習(xí)報告IIR和FIR濾波器設(shè)計專題研討

1、數(shù)字信號處理課程研究性學(xué)習(xí)報告姓名 學(xué)號 同組成員 指導(dǎo)教師 時間 2011/5/10 iir和fir濾波器設(shè)計專題研討【目的】(1) 掌握根據(jù)濾波器指標(biāo)設(shè)計iir和fir數(shù)字濾波器的原理和方法。(2) 熟悉通過iir和fir數(shù)字濾波器進行實際系統(tǒng)設(shè)計的方法。(3) 培養(yǎng)學(xué)生自主學(xué)習(xí)能力，以及發(fā)現(xiàn)問題、分析問題和解決問題的能力?！狙杏戭}目】 1 設(shè)計一個數(shù)字濾波器，在頻率低于的范圍內(nèi)，低通幅度特性為常數(shù)，并且不低于0.75db。在頻率和之間，阻帶衰減至少為20db。(1)試求滿足這些條件的最低階butterworth濾波器。(2)試求滿足這些條件的最低階chebyshev i濾波器。(3)自主

2、選擇一段帶限信號，通過所設(shè)計的（1）、（2）兩種濾波器，比較各自的輸入和輸出信號。討論兩種濾波器在結(jié)構(gòu)和性能上的差異。【題目分析】本題討論模擬濾波器和數(shù)字濾波器的設(shè)計。結(jié)合課本分析各類濾波器的特性?！緄ir模擬濾波器設(shè)計的基本方法】在設(shè)計iir模擬濾波器的時候，1） 先將待設(shè)計的模擬濾波器技術(shù)指標(biāo)轉(zhuǎn)換為原型模擬低通濾波器的指標(biāo)；2） 然后設(shè)計原型模擬低通濾波器（有butterworth、chebyshev和橢圓等設(shè)計方法）；3） 再通過頻率變換將模擬低通濾波器轉(zhuǎn)換為所需的濾波器?！痉抡娼Y(jié)果】【結(jié)果分析】butterworth濾波器:從圖上可以看出，所設(shè)計的濾波器滿足指標(biāo)，通帶內(nèi)為恒定常數(shù)0db

3、，在過渡帶內(nèi)單調(diào)遞減，在pi附近急劇下降，阻帶衰減有裕量。chebyshev i濾波器：從圖上可以看出，所設(shè)計的濾波器滿足指標(biāo)，通帶波動，過渡帶內(nèi)單調(diào)遞減，相對于butterworth濾波器的阻帶衰減比較均勻?！咀灾鲗W(xué)習(xí)內(nèi)容】1） 用butterworth、chebyshev和橢圓等設(shè)計方法設(shè)計模擬低通濾波器；2） 如何用matlab實現(xiàn)上述模擬濾波器的實現(xiàn)；3） 通過頻率變換將模擬低通濾波器轉(zhuǎn)換為所需的濾波器；4） 用單位脈沖響應(yīng)法和雙線性變換法實現(xiàn)模擬濾波器向數(shù)字濾波器的轉(zhuǎn)換?！鹃喿x文獻】數(shù)字信號處理（第二版） 陳后金主編 高等教育出版社【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題

4、)：無【問題探究】單位脈沖響應(yīng)法與雙線性變換法的比較：（使用butterworth的設(shè)計方法討論）仿真程序wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=buttord(wp,ws,ap,as,s);num_a,den_a=butter(n,wc,s);num_d1,den_d1=impinvar(num_a,den_a,fs);num_d2,den_d2=bilinear(num_a,den_a,fs);w=linspace(0,pi,1024);h1=freqz(num_d1,den_d1,w

5、);norm1=max(abs(h1);num_d1=num_d1/norm1;h2=freqz(num_d2,den_d2,w);norm2=max(abs(h2);num_d2=num_d2/norm2;plot(w/pi,20*log10(abs(h1/norm1),blue);hold on;plot(w/pi,20*log10(abs(h2/norm2),green);xlabel(normalized frequency);ylabel(gain,db);仿真結(jié)果結(jié)果分析在本次設(shè)計中，通帶內(nèi)，用單位脈沖響應(yīng)法和雙線性變換區(qū)別并不明顯，相比于單位脈沖響應(yīng)法，雙線性變換法，通帶截頻處地

6、衰減略大，在過渡帶和阻帶則出現(xiàn)明顯的差別，即雙線性變換法的衰減明顯大于單位脈沖響應(yīng)法?！痉抡娉绦颉縝utterworth濾波器:wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=buttord(wp,ws,ap,as,s);num_a,den_a=butter(n,wc,s);num_d,den_d=impinvar(num_a,den_a,fs);w=linspace(0,pi,1024);h=freqz(num_d,den_d,w);norm=max(abs(h);num_d=num_d/nor

7、m;plot(w/pi,20*log10(abs(h/norm);xlabel(normalized frequency);ylabel(gain,db);chebyshev i濾波器：wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=cheb1ord(wp,ws,ap,as,s);num_a,den_a=cheby1(n,ap,wc,s);num_d,den_d=impinvar(num_a,den_a,fs);w=linspace(0,pi,1024);h=freqz(num_d,den_d,w

8、);norm=max(abs(h);num_d=num_d/norm;plot(w/pi,20*log10(abs(h/norm);xlabel(normalized frequency);ylabel(gain,db);濾波：1）butterworth濾波器wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=buttord(wp,ws,ap,as,s);num_a,den_a=butter(n,wc,s);num_d,den_d=impinvar(num_a,den_a,fs); w=linspac

9、e(0,pi,1024);h=freqz(num_d,den_d,w);x=w;s=h.*w;plot(w,abs(s);2）chebyshev i濾波器wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=cheb1ord(wp,ws,ap,as,s);num_a,den_a=cheby1(n,ap,wc,s);num_d,den_d=impinvar(num_a,den_a,fs);w=linspace(0,pi,1024);h=freqz(num_d,den_d,w);x=w;s=h.*w;plo

10、t(w,abs(s);2 分別用hamming窗，blackman窗和kaiser窗設(shè)計，滿足下列指標(biāo)的fir低通濾波器：(1)畫出所設(shè)計濾波器的幅頻響應(yīng)。(2)比較這三種窗的設(shè)計結(jié)果?！绢}目分析】本題討論窗函數(shù)法設(shè)計數(shù)字fir濾波器。結(jié)合課本分析不同窗函數(shù)法的設(shè)計結(jié)果?！緁ir模擬濾波器設(shè)計的基本方法】1、 窗函數(shù)法1） 根據(jù)待逼近的理想濾波器的頻率響應(yīng)，確定濾波器的類型；2） 確定理想濾波器的幅度函數(shù)和相位函數(shù)，由idtft求出理想濾波器的單位脈沖響應(yīng)；3） 將上述所求的無限長單位脈沖響應(yīng)加窗截斷得到有限長的序列。2、 頻率取樣法1） 根據(jù)待逼近的理想濾波器的頻率響應(yīng)，確定濾波器的類型；2

11、） 確定理想濾波器的幅度函數(shù)、相位函數(shù)和頻率響應(yīng)上的取樣點，由idft得到所需的濾波器有限長的單位脈沖響應(yīng)。【仿真結(jié)果】【結(jié)果分析】從上面的仿真可以看出，1） 分別用hamming窗、blackman窗、kaiser窗設(shè)計的fir低通濾波器在通帶內(nèi)沒有波動；2） 通帶衰減和阻帶衰減，kaiser窗最小，hamming窗次之，blackman窗最大；3） 三種窗設(shè)計的階數(shù)由kaiser窗、hamming窗、blackman窗依次增加。【自主學(xué)習(xí)內(nèi)容】1) 窗函數(shù)設(shè)計線性相位fir濾波器的思想和步驟；2) 各種窗函數(shù)的具體內(nèi)容、區(qū)別和使用；3) 頻率頻率取樣法設(shè)計線性相位fir濾波器的思想及步驟；

12、4) 用matlab實現(xiàn)上述兩種方法的設(shè)計方法?！鹃喿x文獻】數(shù)字信號處理（第二版） 陳后金主編 高等教育出版社【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：kaiser窗是一種可調(diào)窗，如何調(diào)節(jié)阻帶衰減？【問題探究】kaiser窗rs參數(shù)的調(diào)節(jié)對濾波器設(shè)計的影響（程序與）仿真結(jié)果 結(jié)果分析1） rs影響阻帶衰減，rs越小，阻帶衰減越大；2） rs影響濾波器的階數(shù)，rs越大，階數(shù)越高；3） rs影響濾波器的阻帶波動范圍。【仿真程序】hamming窗:（n=35）wp=0.4*pi;ws=0.6*pi;ap=0.5;as=45;n=ceil(7*pi/(ws-wp);n=mod(n+1,2

13、)+n;m=n-1;fprintf(n=%.0fn,n);w=hamming(n);wc=(wp+ws)/2;k=0:m;hd=(wc/pi)*sinc(wc*(k-0.5*m)/pi);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb);blackman窗:（n=59）wp=0.4*pi;ws=0.6*pi;ap=0.5;as=45;n=ceil(11.4*pi/(ws-wp);n=mod(n+1,2)+n;m=n-1;fprintf(n=%.0f

14、n,n);w=blackman(n);wc=(wp+ws)/2;k=0:m;hd=(wc/pi)*sinc(wc*(k-0.5*m)/pi);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb);kaiser窗:（n=24）rs=0.01;f=0.4 0.6;a=1,0;dev=rs*ones(1,length(a);m,wc,beta,ftype=kaiserord(f,a,dev);h=fir1(m,wc,ftype,kaiser(m+1,beta

15、);omega=linspace(0,pi,512);mag=freqz(h,1,omega);plot(omega/pi,20*log10(abs(mag);三種窗：rs=0.01;f=0.4 0.6;a=1,0;dev=rs*ones(1,length(a);m1,wc1,beta,ftype=kaiserord(f,a,dev);h1=fir1(m1,wc1,ftype,kaiser(m1+1,beta);omega=linspace(0,pi,512);mag1=freqz(h1,1,omega);plot(omega/pi,20*log10(abs(mag1),blue);hold

16、on;wp=0.4*pi;ws=0.6*pi;ap=0.5;as=45;n2=ceil(7*pi/(ws-wp);n2=mod(n2+1,2)+n;m2=n2-1;fprintf(n2=%.0fn,n2);w2=hamming(n2);wc2=(wp+ws)/2;k=0:m2;hd2=(wc2/pi)*sinc(wc2*(k-0.5*m2)/pi);h2=hd2.*w2;mag2=freqz(h2,1,omega);magdb2=20*log10(abs(mag2);plot(omega/pi,magdb2,green);hold on;wp=0.4*pi;ws=0.6*pi;ap=0.5;a

17、s=45;n3=ceil(11.4*pi/(ws-wp);n3=mod(n3+1,2)+n3;m3=n3-1;fprintf(n3=%.0fn,n3);w3=blackman(n3);wc3=(wp+ws)/2;k=0:m3;hd3=(wc3/pi)*sinc(wc3*(k-0.5*m3)/pi);h3=hd3.*w3;mag3=freqz(h3,1,omega);magdb3=20*log10(abs(mag3);plot(omega/pi,magdb3,red);3附件給出了一段含有噪聲的音頻信號。（1）分析該信號的頻譜特點。（2）通過脈沖響應(yīng)不變法設(shè)計一個濾波器對其進行處理，得到有用信息

18、，自主確定各項指標(biāo)。（3）使用雙線性變換法重做（2）。（4）試定量比較上述兩種濾波器的各項性能，畫出能說明性能差異的相關(guān)圖形，對比并解釋。（5）通過不同窗函數(shù)法設(shè)計fir濾波器對其進行處理，比較不同窗函數(shù)的濾波效果，給出分析和結(jié)論。（選做）(6) 請嘗試采用其它的音頻信號，混入不同的噪聲，利用所學(xué)的濾波方法進行分析，會得到什么樣的效果？【題目分析】本題討論用iir和fir數(shù)字濾波器進行實際系統(tǒng)設(shè)計的方法?！痉抡娼Y(jié)果】(1)【結(jié)果分析】對幾種濾波器得到的音頻信號進行對比?！咀灾鲗W(xué)習(xí)內(nèi)容】這條指令主要功能是讀取音頻文件，該音頻文件的格式規(guī)定為microsoft wave格式，也就是以.wav為后綴

19、的音頻文件。有以下幾種使用方法。1、 y = wavread(file);該語句讀取文件名為file的音頻文件，并將音頻數(shù)據(jù)存儲在y中，返回給用戶使用。音頻數(shù)據(jù)的值域為-1,12、y,fs,nbits = wavread(file)同上，file指定音頻文件名，y是音頻數(shù)據(jù)。另外該語句還能夠返回該音頻文件的采樣率fs,和采樣位數(shù)nbits3、. = wavread(file,n)該語句用來讀取文件名為file的音頻文件中的前n個數(shù)據(jù)。如果有多個通道，每個通道的數(shù)據(jù)都會被讀取出來。.表示隨用戶自己定義。如果用戶在這里只寫一個y，該語句就只會返回音頻數(shù)據(jù)，如果用戶寫成y,fs,nbits,該語句就

20、不但會返回音頻數(shù)據(jù)y,還會返回采樣率fs，和采樣位數(shù)nbits4、.=wavread(file,n1,n2)該語句用來讀取文件名為file的音頻文件中的第n1個到第n2個之間的數(shù)據(jù)。5、siz = wavread(file,size)該語句返回文件名為file的音頻文件的數(shù)據(jù)大小 siz = samples channels6、y,fs,nbits,opts=wavread(.) 該語句返回的參數(shù)opts包含了音頻文件的額外信息。【閱讀文獻】數(shù)字信號處理（第二版） 陳后金主編 高等教育出版社【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：當(dāng)有裕量出現(xiàn)時，選擇不同的窗對濾波器選擇有階數(shù)的

21、影響。【問題探究】在iir濾波器設(shè)計過程中，由于利用指標(biāo)數(shù)據(jù)的不同，造成裕量的出現(xiàn)。討論利用不同指標(biāo)出現(xiàn)裕量對濾波器性能的影響，以及如何有效地利用它？在fir濾波器設(shè)計過程中，由于所選用窗函數(shù)的不同，導(dǎo)致對信號濾波的效果不同。本題的語音信號加入的是雙頻噪聲，討論如果對信號加入其它形式的噪聲，采用何種濾波器更合適?！痉抡娉绦颉浚?）含有噪聲音頻分析:x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); plot(abs(x1);title(含有噪聲的音頻信號);仿真結(jié)果結(jié)果分析讀取音頻文件后，得知抽樣頻率為fs=44100hz,共取n=5

22、0688個點。經(jīng)過將圖放大，在8800hz8950hz，1005hz1025hz兩處范圍內(nèi)有噪聲信號，發(fā)現(xiàn)所需要設(shè)計的濾波器參數(shù)應(yīng)滿足wp6000,ws8800,才能更可能將噪聲去除。(2) 利用脈沖響應(yīng)不變法所設(shè)計的低通濾波器: wp=6336;ws=8448;ap=1;as=25;（讓角頻率可以取整數(shù)）x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); w=linspace(0,2*pi,50688);subplot(3,1,1);plot(w/pi,abs(x1);title(含有噪聲的音頻信號頻譜圖);np=6336;ns=8

23、448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;subplot(3,1,2);plot(w/pi,20*log10(abs(h)/norm); title(利用脈沖響應(yīng)不變法

24、所設(shè)計的濾波器幅度分布);x=x1.*(h);subplot(3,1,3);plot(w/pi,abs(x);title(濾除噪聲后的音頻信號頻譜圖);仿真結(jié)果結(jié)果分析要使噪聲消除，則濾波器在噪聲頻率處需要完成衰減，此音頻信號中，需要完成的衰減約為49.54db，(3) 利用雙線性變換法所設(shè)計的低通濾波器: wp=6336;ws=8448;ap=1;as=25;（讓角頻率可以取整數(shù)）x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); w=linspace(0,2*pi,50688);subplot(3,1,1);plot(w/pi,a

25、bs(x1);title(含有噪聲的音頻信號頻譜圖);np=6336;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;subplot(3,1,2);plot(w/pi

26、,20*log10(abs(h)/norm); title(利用雙線性變化法所設(shè)計的濾波器幅度分布);x=x1.*(h);subplot(3,1,3);plot(w/pi,abs(x);title(濾除噪聲后的音頻信號頻譜圖);結(jié)果分析要使噪聲消除，則濾波器在噪聲頻率處需要完成衰減，此音頻信號中，需要完成的衰減約為49.54db，當(dāng)使用雙線性變換法時，在此噪音頻率前已已達到標(biāo)準(zhǔn)，過渡帶衰減速率高。(4) 定量比較上述兩種濾波器的各項性能:從第一題中可以看出，兩方法的最大區(qū)別是過渡帶的差距，則先保持ap=1;as=25不變，調(diào)大wp，使通帶增加， np=6336w=linspace(0,2*pi

27、,50688);np=6336;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) ,blue); xla

28、bel(normalized frequency);ylabel(gain,db); hold on;np=6336;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;

29、plot(w/pi,20*log10(abs(h)/norm) ,red); xlabel(normalized frequency);ylabel(gain,db);legend(脈沖響應(yīng)不變法, 雙線性變化法)title(兩種濾波器np=6336);np=6500w=linspace(0,2*pi,50688);np=6500;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);nu

30、ma,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=6500;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(w

31、p,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脈沖響應(yīng)不變法, 雙線性變化法)title(兩種濾波器np=6500);np=7000w=linspace(0,

32、2*pi,50688);np=7000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue)

33、; xlabel(normalized frequency);ylabel(gain,db); hold on;np=7000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/

34、norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脈沖響應(yīng)不變法, 雙線性變化法)title(兩種濾波器np=7000);np=7500w=linspace(0,2*pi,50688);np=7500;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2

35、/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=7200;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=but

36、tord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脈沖響應(yīng)不變法, 雙線性變化法)title(兩種濾波器np=7200);np=7100w=linsp

37、ace(0,2*pi,50688);np=7100;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) ,

38、 blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=7100;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd

39、=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脈沖響應(yīng)不變法, 雙線性變化法)title(兩種濾波器np=7100);np=8000w=linspace(0,2*pi,50688);np=8000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-

40、1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=8000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs

41、;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脈沖響應(yīng)不變法, 雙線性變化法)title(兩種濾波器np=8000); 結(jié)果分析當(dāng)ws，as,ap不變時，wp增大，兩種方法的過渡帶變小，阻帶變寬，但是用脈沖響應(yīng)法做出的變化劇烈。（5）hamming窗:（n=85）x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); w=linspace(0,2*pi,50688);subplot(3,1,1);plot(w/pi,abs(x1);title(含有噪聲的音頻信號頻譜圖);np=6336;ns=8448;n=