1、IIR和FIR濾波器設計專題研討第八組：彭然 09214040 海驍飛09214032 季凱09214035 劉炘09214036 1 設計一個數字濾波器，在頻率低于 的范圍內，低通幅度特性近似為常數，并且不低于0.75dB。在頻率 和 之間，阻帶衰減至少為20dB。(1)試求滿足這些條件的最低階Butterworth濾波器。(2)試求滿足這些條件的最低階Chebyshev I濾波器。(3)自主選擇一段帶限信號，通過所設計的（1）、（2）兩種濾波器，比較各自的輸入和輸出信號。討論兩種濾波器在結構和性能上的差異。0.2613 0.4018 【題目分析題目分析】本題討論模擬濾波器和數字濾波器的設計

2、。結合課本分析各類濾波器的特性。分析：Butterworth濾波器幅度響應在通帶和阻帶都是隨頻率單調變化，在通帶和阻帶都存在裕量，并且階數較高；Chebyshev I濾波器在通帶具有等波紋特性，在阻帶單調下降，因此通帶裕量較Butterworth濾波器要小，但是阻帶較Butterworth濾波器要大，階數要低。【IIR模擬濾波器設計的基本方法模擬濾波器設計的基本方法】模擬濾波器的設計就是求出滿足技術指標的連續系統函數H（s）,模擬濾波器的技術指標有通帶截頻wp、通帶最大衰減Ap、阻帶截頻ws和阻帶最小衰減As。在設計模擬濾波器時，先將待設計的模擬濾波器技術指標轉換為模擬低通濾波器技術指標，然后

3、設計模擬低通濾波器，再通過頻率變換將模擬低通濾波器轉換為所需的濾波器。其中模擬低通濾波器的設計有一些模型，比如Butterworth、Chebyshev等等模擬濾波器模型。（1）Butterworth濾波器:%ButterworthWp=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);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(

4、numd,dend,w);norm=max(abs(h);plot(w/pi,20*log10(abs(h)/norm);xlabel(Normalized frequency);ylabel(Gain,dB);grid on;w=Wp Ws;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2); N=8 Ap= 0.4192 As= 19.9999(2) Chebyshev I濾波器:%Chebyshev IWp=0.261

5、3*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);numa,dena=cheby1(N,Ap,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);plot(w/pi,20*log10(abs(h)/norm);xlabel(Normalized frequency);ylabel(Gain,dB);grid on;w=Wp W

6、s;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);N=4 Ap= 0.7523 As= 21.3368（3）將帶限信號x=cos(0.5*pi*t)通過Butterworth濾波器觀察濾波效果：%ButterworthWp=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);numa,d

7、ena=butter(N,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);H=freqz(numd,dend,w);norm=max(abs(H);numd=numd/norm;H=freqz(numd,dend,w);t=0:0.01:2*pi;x=cos(0.5*pi*t);x1=filter(numd,dend,x);subplot(2,1,1);plot(t,x);hold on;plot(t,x1,r);title(Butterworth);grid on;subplot(2,1,2);plot(t,x1-x)

8、;title();grid on;將帶限信號x=cos(0.5*pi*t)通過Chebyshev I濾波器觀察濾波效果：%Chebyshev IWp=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);numa,dena=cheby1(N,Ap,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);H=freqz(numd,dend,w);norm=max(abs(H);numd=numd/n

9、orm;H=freqz(numd,dend,w);t=0:0.01:2*pi;x=cos(0.5*pi*t);x1=filter(numd,dend,x);subplot(2,1,1);plot(t,x);hold on;plot(t,x1,r);title(Chebyshev I);grid on;subplot(2,1,2);plot(t,x1-x);title();grid on;【結果分析結果分析】(1) Butterworth濾波器幅度響應在通帶和阻帶都是隨頻率單調變化，在通帶和阻帶都存在裕量，設計要求是Ap=0.75，As=20，而結果是Ap= 0.4192 As= 19.9999

10、，并且階數較高，N=4；（2）Chebyshev I濾波器在通帶具有等波紋特性，在阻帶單調下降，因此通帶裕量較Butterworth濾波器要小，但是阻帶較Butterworth濾波器要大，但是結果是Ap= 0.7495 As= 25.5334，通帶衰減不滿足題目要求，阻帶存在較大的裕量。（3）通過分析濾波器后時域誤差可知：Butterworth濾波器產生的誤差大概是0.08498，要小于Chebyshev I濾波器的0.096。同時可知，時域波形在0處誤差幾乎是1，我覺得這個是由于余弦信號時域直接加了矩形窗，在頻域為二者卷積，可看成是兩個sa函數形式的波，其幅值與矩形窗長度有關，其N越大，幅值

11、越大，那么即使在濾波時幅值下降了，那么這個可以補償，其時域也幾乎一樣。當t=0:0.01:2*pi;x=cos(0.5*pi*t);誤差上圖可知，當t=0:0.001:2*pi;x=cos(0.5*pi*t);2 分別用Hamming窗，Blackman窗和Kaiser窗設計，滿足下列指標的FIR低通濾波器：0.4 ,0.6 ,0.5,45pspsAdBAdB (1)畫出所設計濾波器的幅頻響應。(2)比較這三種窗的設計結果。【FIR模擬濾波器設計的基本方法模擬濾波器設計的基本方法】FIR數字濾波器的設計通常是根據理想濾波器的頻率響應 ，采用窗函數法、頻率取樣法或優化方法使設計的濾波器的頻率響應

12、 逼近 。jH ejH ejdHe%hammingFIRWp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;N=ceil(7*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;fprintf(N=%.0fn,N);%generate windoww=hamming(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb)

13、;xlabel(hammingFIR);grid on; （1）Hamming窗：（1）Hamming窗：N=35 As=55（2）Blackman窗：%blackmanFIRWp=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=%.0fn,N);w=blackman(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,ome

14、ga);magdb=20*log10(abs(mag);plot(omega/pi,magdb);xlabel(blackmanFIR);grid on;（2）Blackman窗：N=27 As=45%Program:KaiserFIRAp=0.5;As=45;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*As);f=0.4,0.6;a=1,0;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f,a,dev);M=mod(M,2)+M;fprintf(N=%.0fn,M+1);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);xlabel(KaiserFIR);grid on;（3）Kaiser窗：（3）Kaiser窗：N=59 As=78.11【結果分析結果分析】Hamming窗，Blackman窗和Kaiser窗三種窗中Blackman窗的階數最小，衰減也較快，Kaiser窗衰減最大3附件給出了一段含有噪聲的音頻信號。（1）分析該信號的頻譜特點。（2）通過脈沖響應不變法設計一個濾波器對其進行處理，得到有用信息，自主確定各項指標。（3）使用雙線性變換法重做（

16、2）。（4）試定量比較上述兩種濾波器的各項性能，畫出能說明性能差異的相關圖形，對比并解釋。（5）通過不同窗函數法設計FI R濾波器對其進行處理，比較不同窗函數的濾波效果，給出分析和結論。（選做）(6) 請嘗試采用其它的音頻信號，混入不同的噪聲，利用所學的濾波方法進行分析，會得到什么樣的效果？y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);%subplot(2,1,2);plot(abs(Y);xlabel();subplot(2,1,1);plot(y);xlabel();（1）信號的時域波形和頻譜通過對頻譜分析確定設計低通濾波器，其中通帶截止

17、頻率Fp=1400，阻帶截止頻率是Fs=1600Ap=1; As=20，用脈沖響應不變法：y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;Wp=2*pi*Fp/Ft; Ws=2*pi*Fs/Ft; Ap=1; As=20;fs=1; wp=Wp*fs; ws=Ws*fs;N,wc=buttord(wp,ws,Ap,As,s);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,fs);w=linspace(0,pi,1024);h=f

18、reqz(numd,dend,w);norm=max(abs(h);w=Wp Ws;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);y1=filter(numd,dend,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title(

19、);grid;subplot(2,2,1);plot(abs(Y),r);title();grid;（3）雙線性變化法：y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;wp=2*Ft*tan(Wp/2);ws=2*Ft*tan(Ws/2);N,wc=buttord(wp,ws,Ap,As,s); %num,den=butter(N,wc,s); %S numd,dend=bilinear(

20、num,den,Ft); %SZ h,w=freqz(numd,dend); %w=Wp Ws;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);y1=filter(numd,dend,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title();grid;su

21、bplot(2,2,1);plot(abs(Y),r);title();grid;（4）脈沖響應不變法與雙線性變換法性能比較y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;Wp1=2*pi*Fp/Ft; Ws1=2*pi*Fs/Ft; Ap=1; As=20;fs=1; wp1=Wp1*fs; ws1=Ws1*fs;N1,wc1=buttord(wp1,ws1,Ap,As,s);numa1,dena1=butter(N1,wc1,s);numd1,dend1=impinvar(numa1,

22、dena1,fs);w=linspace(0,pi,1024);h1=freqz(numd1,dend1,w);norm=max(abs(h1);numd1=numd1/norm;h1=freqz(numd1,dend1,w);Wp2=2*pi*Fp/Ft;Ws2=2*pi*Fs/Ft;wp2=2*Ft*tan(Wp2/2);ws2=2*Ft*tan(Ws2/2);N2,wc2=buttord(wp2,ws2,Ap,As,s); num2,den2=butter(N2,wc2,s); numd2,dend2=bilinear(num2,den2,Ft); h2,w=freqz(numd2,de

23、nd2); y1=filter(numd1,dend1,y);y2=filter(numd2,dend2,y);plot(abs(y1-y2);xlabel();（5）通過不同窗函數法設計FI R濾波器處理：Hamming窗%HammingFIRy,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;N=ceil(7*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;fprintf(

24、N=%.0fn,N);%generate windoww=hamming(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;y1=fftfilt(h,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title();grid;subplot(2,2,1);plot(abs(Y),r);title();grid;Blackman窗%B

25、lackmanFIRy,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;N=ceil(11.4*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;fprintf(N=%.0fn,N);w=blackman(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;y1=fftfilt(h,y);figure(1)

26、subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title();grid;subplot(2,2,1);plot(abs(Y),r);title();grid;Kaiser窗%KaiserFIRy,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Rp=1-10.(-0.05*Ap);

27、Rs=10.(-0.05*As);f=2800/16000,3200/16000;a=1,0;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f,a,dev);M=mod(M,2)+M;fprintf(N=%.0fn,M+1);h = fir1(M,Wc,ftype,kaiser(M+1,beta);y1=fftfilt(h,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);ti

28、tle();grid;subplot(2,2,1);plot(abs(Y),r);title();grid;性能對比：y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;N1=ceil(7*pi/(Ws-Wp);N1=mod(N1+1,2)+N1;M1=N1-1;fprintf(N1=%.0fn,N1);%generate windoww1=hamming(N1);Wc=(Wp+Ws)/2;k

29、=0:M1;hd1=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M1);h1=hd1.*w1;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*As);f=2800/16000,3200/16000;a=1,0;dev=Rp,Rs;M2,Wc,beta,ftype = kaiserord(f,a,dev);M2=mod(M2,2)+M2;fprintf(N2=%.0fn,M2+1);h2= fir1(M2,Wc,ftype,kaiser(M2+1,beta); N3=ceil(11.4*pi/(Ws-Wp);N3=mod(N3+1,2)+N3;M3=N3-1;fprintf(N3=%.0fn,N3);w3=blackman(N3);Wc=(Wp+Ws)/2;k=0:M3;hd3=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M3);h3=hd3.*w3; y1=fftfilt(h1,y);