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

1、數(shù)字信號處理課程研究性學(xué)習(xí)報告IIR和FIR濾波器設(shè)計(jì)專題研討【目的】(1) 掌握根據(jù)濾波器指標(biāo)設(shè)計(jì)IIR和FIR數(shù)字濾波器的原理和方法。(2) 熟悉通過IIR和FIR數(shù)字濾波器進(jìn)行實(shí)際系統(tǒng)設(shè)計(jì)的方法。(3) 培養(yǎng)學(xué)生自主學(xué)習(xí)能力，以及發(fā)現(xiàn)問題、分析問題和解決問題的能力?！狙杏戭}目】 1 設(shè)計(jì)一個數(shù)字濾波器，在頻率低于的范圍內(nèi)，低通幅度特性近似為常數(shù)，并且不低于0.75dB。在頻率和之間，阻帶衰減至少為20dB。(1)試求滿足這些條件的最低階Butterworth濾波器。(2)試求滿足這些條件的最低階Chebyshev I濾波器。(3)自主選擇一段帶限信號，通過所設(shè)計(jì)的（1）、（2）兩種濾波器

2、，比較各自的輸入和輸出信號。討論兩種濾波器在結(jié)構(gòu)和性能上的差異?！绢}目分析】本題討論模擬濾波器和數(shù)字濾波器的設(shè)計(jì)。結(jié)合課本分析各類濾波器的特性?！綢IR模擬濾波器設(shè)計(jì)的基本方法】【仿真結(jié)果】1)N=8N=2N=12)N=4N=2N=1【結(jié)果分析】BW的N CB1的N可以根據(jù)分式的關(guān)系來確定N的最小值BW的N和Ap成正比CB1的N和Wp成反比在對其取值的時候可以根據(jù)這樣的規(guī)律實(shí)驗(yàn)，模擬的結(jié)果和理論是相同的。從圖像上可以看出N的值就是通頻帶的起伏的次數(shù)，N=1就起伏一次，N=2就起伏2次，N=8就接近連成一條線?！咀灾鲗W(xué)習(xí)內(nèi)容】【閱讀文獻(xiàn)】【仿真程序】1)Wp=0.2613*pi; Ws=0.40

3、18*pi; Ap=0.75; As=20;Fs=1; wp=Wp*Fs; ws=Ws*Fs;N=buttord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);numa,dena=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;plot(w/pi,20*log10(abs(h/norm);xlabel( butter );ylabel(Ap=0.75,dB);d

4、isp(N=);disp(N);2)Wp=0.2613*pi; Ws=0.4018*pi; Ap=0.75; As=20;Fs=1; wp=Wp*Fs; ws=Ws*Fs;N=cheb1ord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);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/norm;plot(w/pi,20*log10(abs(h/

5、norm);xlabel(cheby1);ylabel(Ap=0.75,dB);disp(N=);disp(N);2 分別用Hamming窗，Blackman窗和Kaiser窗設(shè)計(jì)，滿足下列指標(biāo)的FIR低通濾波器：(1)畫出所設(shè)計(jì)濾波器的幅頻響應(yīng)。(2)比較這三種窗的設(shè)計(jì)結(jié)果?！绢}目分析】本題討論窗函數(shù)法設(shè)計(jì)數(shù)字FIR濾波器。結(jié)合課本分析不同窗函數(shù)法的設(shè)計(jì)結(jié)果?！綟IR模擬濾波器設(shè)計(jì)的基本方法】【仿真結(jié)果】【結(jié)果分析】【仿真程序】Wp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;%hanningfigure(1)N=floor(6.2*pi/(Ws-Wp);N=mod(N+1,

6、2)+N;M=N-1;w=hanning(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);plot(omega/pi,20*log10(abs(mag);grid onxlabel(Normalized frequency)ylabel(Gain in dB)title(hanning) %hammingfigure(2)N=floor(7*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1w=hamming(N

7、);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);plot(omega/pi,20*log10(abs(mag);grid onxlabel(Normalized frequency)ylabel(Gain in dB)title(hamming) %blackmanfigure(3)N=floor(11.4*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;w=blackman(N);Wc=(Wp+Ws)/2;k

8、=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);plot(omega/pi,20*log10(abs(mag);grid onxlabel(Normalized frequency)ylabel(Gain in dB)title(blackman) %kaiserfigure(4)f=Wp/pi Ws/pi;a=1,0;dev=10(-As/20)*ones(1,length(a);M,Wc,beta,ftype=kaiserord(f,a,dev);M=m

9、od(M,2)+M;h=fir1(M,Wc,ftype,kaiser(M+1,beta);omega=linspace(0,pi,512);mag=freqz(h,1,omega);plot(omega/pi,20*log10(abs(mag);grid onzoom onxlabel(Normalized frequency)ylabel(Gain in dB)title(kaiser)3附件給出了一段含有噪聲的音頻信號。（1）分析該信號的頻譜特點(diǎn)。（2）通過脈沖響應(yīng)不變法設(shè)計(jì)一個濾波器對其進(jìn)行處理，得到有用信息，自主確定各項(xiàng)指標(biāo)。（3）使用雙線性變換法重做（2）。（4）試定量比較上述兩種濾

10、波器的各項(xiàng)性能，畫出能說明性能差異的相關(guān)圖形，對比并解釋。（5）通過不同窗函數(shù)法設(shè)計(jì)FIR濾波器對其進(jìn)行處理，比較不同窗函數(shù)的濾波效果，給出分析和結(jié)論。（選做）(6) 請嘗試采用其它的音頻信號，混入不同的噪聲，利用所學(xué)的濾波方法進(jìn)行分析，會得到什么樣的效果？【題目分析】本題討論用IIR和FIR數(shù)字濾波器進(jìn)行實(shí)際系統(tǒng)設(shè)計(jì)的方法?！痉抡娼Y(jié)果】z1= wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);plot(z1);z1= wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);n=0:81

11、91;plot(n,Y1);（1）窗函數(shù)設(shè)計(jì)低通濾波器 程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;Fs=8000;wc=2*pi*fc/Fs; wp=2*pi*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7);N=ceil(As-8)/2.285/wdel);wn= kaiser(N+1,beta); ws=(wp+wc)/2/pi;b=fir1(N,ws,wn

12、);figure(1);freqz(b,1);x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(濾波后信號頻譜);subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot(2,2,4);plot(x);title(濾波后信號波形);sound(x,fs,bits);wavplay(x,44100);wav

13、write(x,44100,1b.wav);fname=sprintf(E:1b.wav); wavwrite(x,44100,fname);圖形分析如下：（2）窗函數(shù)設(shè)計(jì)高通濾波器程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fp=2800;fc=3000;As=100;Ap=1;Fs=8000;wc=2*pi*fc/Fs; wp=2*pi*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7);N=ceil(As-8)/2

14、.285/wdel);wn= kaiser(N,beta); ws=(wp+wc)/2/pi;b=fir1(N-1,ws,high,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(濾波后信號頻譜);subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot

15、(2,2,4);plot(x);title(濾波前信號波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,2b.wav);fname=sprintf(E:2b.wav); wavwrite(x,44100,fname);（3）窗函數(shù)設(shè)計(jì)帶通濾波器程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fp1=1200 ;fp2=3000 ;fc1=1000 ;fc2=3200 ;As=100 ;A

16、p=1 ;Fs=8000 ;wp1=2*pi*fp1/Fs; wc1=2*pi*fc1/Fs; wp2=2*pi*fp2/Fs; wc2=2*pi*fc2/Fs;wdel=wp1-wc1;beta=0.112*(As-8.7);N=ceil(As-8)/2.285/wdel);ws =(wp1+wc1)/2/pi,(wp2+wc2)/2/pi;wn= kaiser(N+1,beta); b=fir1(N,ws,wn);figure(1);freqz(b,1)x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);a

17、xis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,2000,0,0.0003);title(濾波后信號頻譜)subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot(2,2,4);plot(x);title(濾波前信號波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,3b.wav);fname=sprintf(E:3b.wav); wavwrite(x,44100,fname);圖形分析如下：（4）雙線性變換法設(shè)計(jì)低

18、通濾波器選用butter程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;Fs=8000;wc=2*fc/Fs; wp=2*fp/Fs; N,ws=buttord(wc,wp,Ap,As);b,a=butter(N,ws);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);pl

19、ot(abs(Y1);axis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(濾波后信號頻譜);subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot(2,2,4);plot(x);title(濾波前信號波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,4b.wav);fname=sprintf(E:4b.wav); wavwrite(x,44100,fname);圖形分析選用che

20、by1程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1; ;Fs=8000;wc=2*fc/Fs;wb=2*fp/Fs;n,wp=cheb1ord(wc,wb,Ap,As);b,a=cheby1(n,Ap,wp);figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);a

21、xis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,4000,0,0.03);title(濾波后信號頻譜);subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot(2,2,4);plot(x);title(濾波前信號波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,5b.wav);fname=sprintf(E:5b.wav); wavwrite(x,44100,fname);圖形分析如下：5），雙線性變換法設(shè)計(jì)高通

22、濾波器選用butter程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;As=100;Ap=1;Fs=8000;wc=2*fc/Fs; wp=2*fp/Fs; N,ws=buttord(wc,wp,Ap,As);b,a=butter(N,ws,high);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2

23、,1);plot(abs(Y1);axis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(濾波后信號頻譜);subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot(2,2,4);plot(x);title(濾波前信號波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,6b.wav);fname=sprintf(E:6b.wav); wavwrite(x,44100,fname);圖形分

24、析如下：選用cheby1程序設(shè)計(jì)如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;As=100;Ap=1; Fs=8000;wc=2*fc/Fs;wb=2*fp/Fs;n,wp=cheb1ord(wc,wb,Ap,As);b,a=cheby1(n,Ap,wp,high);figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1)

25、;plot(abs(Y1);axis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,4000,0,0.03);title(濾波后信號頻譜);subplot(2,2,3);plot(z1);title(濾波前信號波形);subplot(2,2,4);plot(x);title(濾波前信號波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,7b.wav);fname=sprintf(E:7b.wav); wavwrite(x,44100,fname);圖形分析如

26、下：（6），雙線性變換法設(shè)計(jì)帶通濾波器選用butter程序設(shè)計(jì)如下clear;close allz1,fs,bits=wavread(C:UserszhangDesktop課件DSP你好加噪聲信號);y1=z1(1:8192);Y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;As=100;Ap=1; Fs=8000; wc=2*fc1/Fs,2* fc2/Fs;wp=2*fp1/Fs,2*fp2/Fs; N,ws=buttord(wc,wp,Ap,As);b,a=butter(N,ws,stop);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,1.0);title(濾波前信號頻譜);subplot(2,2,2);plot(abs(X);axis(0,100