數(shù)字濾波器設(shè)計(jì)和應(yīng)用

1、實(shí)驗(yàn)四 數(shù)字濾波器設(shè)計(jì)及應(yīng)用201030021717 電信轉(zhuǎn)專業(yè)班02號(hào) 王澤軍一、實(shí)驗(yàn)內(nèi)容（一）1，用窗函數(shù)法設(shè)計(jì)如下三種類型的FIR濾波器，濾波器的性能指標(biāo)： 低通濾波器：通帶截止頻率1000Hz，阻帶截止頻率1200Hz，通帶波紋1dB，最小阻帶衰減60dB； 高通濾波器：通帶截止頻率5000Hz，阻帶截止頻率4800Hz，通帶波紋1dB，最小阻帶衰減60dB； 帶通濾波器：通帶截止頻率fp1=1200Hz, fp2=3000Hz，阻帶截止頻率fs1=1000Hz， fs2=3200Hz, 通帶波紋1dB，最小阻帶衰減60dB。 (1)設(shè)計(jì)原理令希望設(shè)計(jì)的濾波器的傳輸函數(shù)是，是與其對(duì)應(yīng)的

2、單位脈沖響應(yīng)。一般情況下，由求出，然后由Z變換求出濾波器的系統(tǒng)函數(shù)。但是通常在邊界頻率處有不連續(xù)點(diǎn)，這使得是無限長的非因果序列，所以實(shí)際是不能實(shí)現(xiàn)的。為了構(gòu)造一個(gè)長度為N的線性相位濾波器，可以將截取一段來近似，并且根據(jù)線性相位的特點(diǎn)，需要保證截取后的序列關(guān)于對(duì)稱。設(shè)截取的一段為，則其中稱為矩形窗函數(shù)。當(dāng)?shù)膶?duì)稱中心點(diǎn)取值為時(shí)，就可以保證所設(shè)計(jì)的濾波器具有線性相位。常用的窗函數(shù)除了矩形窗還有三角形窗、漢寧窗、哈明窗和布萊克曼窗。在MATLAB中提供了相應(yīng)的子程序來實(shí)現(xiàn)這些窗函數(shù)。 w=boxcar(N)，在數(shù)組w中產(chǎn)生N點(diǎn)的矩形窗函數(shù)。 w=bartlett(N)，在數(shù)組w中產(chǎn)生N點(diǎn)的三角形窗函數(shù)

3、。 w=hanning(N)，在數(shù)組w中產(chǎn)生N點(diǎn)的漢寧窗函數(shù)。 w=hamming(N)，在數(shù)組w中產(chǎn)生N點(diǎn)的哈明窗函數(shù)。 w=blackman(N)，在數(shù)組w中產(chǎn)生N點(diǎn)的布萊克曼窗函數(shù)。這5種窗函數(shù)的具體性能見表1。表1 種窗函數(shù)性能比較窗函數(shù)旁瓣峰值dB主瓣寬度(近似過渡帶寬)精確過渡帶寬阻帶最小衰減dB矩形窗-1341.8-21三角形窗-2586.1-25漢寧窗-3186.2-44哈明窗-4186.6-53布萊克曼窗-571211-74(2) FIR濾波器的設(shè)計(jì) 設(shè)計(jì)低通FIR濾波器，設(shè)計(jì)條件：通帶截止頻率1000Hz，阻帶截止頻率1200Hz，通帶波紋1dB，最小阻帶衰減60dB；最小

4、阻帶衰減為60dB，對(duì)照表1，只能選擇布萊克曼窗函數(shù)進(jìn)行設(shè)計(jì)，假設(shè)采樣率，則歸一化的通帶截止角頻率。歸一化的阻帶截止角頻率。由此確定FIR低通濾波器的截止頻率為，設(shè)計(jì)的，編寫如下MATLAB程序：% 數(shù)字低通FIR濾波器的設(shè)計(jì)clearFs=16000;fp=1000;fs=1200;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;wc=(wp+ws)/2; % 截止頻率的確定dw=ws-wp;M=ceil(5.56*pi/dw); % 布萊克曼窗的參數(shù)Mb=blackman(2*M+1);hlp=fir1(2*M,wc/pi,b);wvtool(hlp)程序運(yùn)行后可以得到低通濾波器的

5、時(shí)域和頻域特性如下圖1：圖1：低通FIR濾波器的時(shí)頻特性 設(shè)計(jì)高通FIR濾波器，設(shè)計(jì)條件：通帶截止頻率5000Hz，阻帶截止頻率4800Hz，通帶波紋1dB，最小阻帶衰減60dB；最小阻帶衰減為60dB，所以只能選擇布萊克曼窗函數(shù)，假設(shè)采樣率，編寫MATLAB程序如下：% 數(shù)字高通FIR濾波器的設(shè)計(jì)clearFs=16000;fp=5000;fs=4800;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;wc=(wp+ws)/2; % 截止頻率的確定dw=wp-ws;M=ceil(5.56*pi/dw); % 布萊克曼窗的參數(shù)Mb=blackman(2*M+1);hlp=fir1(2*

6、M,wc/pi,'high',b);wvtool(hlp)程序運(yùn)行后可以得到高通濾波器的時(shí)域和頻域特性如下圖2：圖2：高通FIR濾波器的時(shí)頻特性 設(shè)計(jì)FIR帶通濾波器，設(shè)計(jì)條件：通帶截止頻率fp1=1200Hz, fp2=3000Hz，阻帶截止頻率fs1=1000Hz，fs2=3200Hz,通帶波紋1dB，最小阻帶衰減60dB。最小阻帶衰減為60dB，所以只能選擇布萊克曼窗函數(shù)，假設(shè)采樣率，編寫MATLAB程序如下：% 數(shù)字帶通FIR濾波器的設(shè)計(jì)clearFs=16000;fp1=1200;fp2=3000;fs1=1000;fs2=3200;wp1=2*pi*fp1/Fs;w

7、p2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;wc1=(wp1+ws1)/2;wc2=(wp2+ws2)/2; % 截止頻率的確定dw=max(abs(wp1-ws1),abs(wp2-ws2);M=ceil(5.56*pi/dw); % 布萊克曼窗的參數(shù)Mb=blackman(2*M+1);fc=wc1,wc2/pi;hlp=fir1(2*M,fc,'high',b);wvtool(hlp)程序運(yùn)行后可以得到帶通濾波器的時(shí)域和頻域特性如下圖3：圖3：帶通FIR濾波器的時(shí)頻特性2，再用雙線性變換法設(shè)計(jì)上面的三種IIR濾波器（可以選

8、擇butterworth or chebyshey），并用Freqz函數(shù)畫出各個(gè)濾波器的頻率響應(yīng)。答：在實(shí)驗(yàn)三中已經(jīng)詳細(xì)給出了如何利用雙線性變換法設(shè)計(jì)IIR濾波器，這里不再贅述，只給出相應(yīng)的程序和結(jié)果。(1)分別用butterworth、chebyshey1型和chebyshey2型設(shè)計(jì)低通IIR濾波器，設(shè)計(jì)條件：通帶截止頻率1000Hz，阻帶截止頻率1200Hz，通帶波紋1dB，最小阻帶衰減60dB；取采樣率。編寫MATLAB程序如下：% 用butteworth設(shè)計(jì)低通IIR濾波器cleara=1;b=60;Fs=8000;fp=1000;fs=1200;wp=2*pi*fp/Fs;ws=2

9、*pi*fs/Fs;Wp=tan(wp/2); % 模擬頻帶邊界頻率的確定Ws=tan(ws/2);N,Wn=buttord(Wp,Ws,a,b,'s'); % 確定巴特沃茲低通濾波器的截止頻率num,den=butter(N,Wn,'s'); p2,p1= bilinear(num,den,0.5); % 從s域到z域進(jìn)行雙線性變換Hlp,w=freqz(p2,p1);figure(1)plot(w/2/pi*Fs,abs(Hlp)grid onxlabel('f/Hz'),ylabel('幅度')title('用but

10、teworth設(shè)計(jì)低通IIR濾波器')% 用chebyshey1型設(shè)計(jì)低通IIR濾波器N1,Wc1=cheb1ord(Wp,Ws,a,b,'s'); % 確定切比雪夫1型濾波器的階數(shù)和截止頻率num1,den1=cheby1(N1,a,Wc1,'s'); p4,p3=bilinear(num1,den1,0.5); % 從s域到z域進(jìn)行雙線性變換Hlp1,w1=freqz(p4,p3);figure(2)plot(w1/2/pi*Fs,abs(Hlp1)grid onxlabel('f/Hz'),ylabel('幅度')t

11、itle('用chebyshey1型設(shè)計(jì)低通IIR濾波器')% 用chebyshey2型設(shè)計(jì)低通IIR濾波器N2,Wc2=cheb2ord(Wp,Ws,a,b,'s'); % 確定切比雪夫2型濾波器的階數(shù)和截止頻率num2,den2=cheby2(N2,b,Wc2,'s');p6,p5=bilinear(num2,den2,0.5); % 從s域到z域進(jìn)行雙線性變換Hlp2,w2=freqz(p6,p5);figure(3)plot(w2/2/pi*Fs,abs(Hlp2)grid onxlabel('f/Hz'),ylabel(

12、'幅度')title('用chebyshey2型設(shè)計(jì)低通IIR濾波器')程序運(yùn)行后可以分別得到用butterworth、chebyshey1型和chebyshey2型設(shè)計(jì)的低通IIR濾波器頻域特性如下圖4：圖4-1：用butterworth設(shè)計(jì)低通IIR濾波器的頻域特性圖4-2：用chebyshey1型設(shè)計(jì)低通IIR濾波器的頻域特性圖4-3：用chebyshey2型設(shè)計(jì)低通IIR濾波器的頻域特性(2) 分別用butterworth、chebyshey1型和chebyshey2型設(shè)計(jì)高通IIR濾波器，設(shè)計(jì)條件：通帶截止頻率5000Hz，阻帶截止頻率4800Hz，通

13、帶波紋1dB，最小阻帶衰減60dB；取采樣率。編寫MATLAB程序如下：% 用butteworth設(shè)計(jì)高通IIR濾波器cleara=1;b=60;Fs=16000;fp=5000;fs=4800;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp=tan(wp/2); % 模擬頻帶邊界頻率的確定Ws=tan(ws/2);N,Wn=buttord(Wp,Ws,a,b,'s'); % 確定巴特沃茲高通濾波器的截止頻率num,den=butter(N,Wn,'high','s'); p2,p1= bilinear(num,den,0.5);

14、% 從s域到z域進(jìn)行雙線性變換Hhp,w=freqz(p2,p1);figure(1)plot(w/2/pi*Fs,20*log10(abs(Hhp)grid onxlabel('f/Hz'),ylabel('幅度')title('用butteworth設(shè)計(jì)高通IIR濾波器')% 用chebyshey1型設(shè)計(jì)高通IIR濾波器N1,Wc1=cheb1ord(Wp,Ws,a,b,'s'); % 確定切比雪夫1型濾波器的階數(shù)和截止頻率num1,den1=cheby1(N1,a,Wc1,'high','s'

15、); p4,p3= bilinear(num1,den1,0.5); % 從s域到z域進(jìn)行雙線性變換Hhp1,w1=freqz(p4,p3);figure(2)plot(w1/2/pi*Fs,abs(Hhp1)grid onxlabel('f/Hz'),ylabel('幅度')title('用chebyshey1型設(shè)計(jì)高通IIR濾波器')% 用chebyshey2型設(shè)計(jì)高通IIR濾波器N2,Wc2=cheb2ord(Wp,Ws,a,b,'s'); % 確定切比雪夫1型濾波器的階數(shù)和截止頻率num2,den2=cheby2(N2,b

16、,Wc2,'high','s'); p6,p5=bilinear(num2,den2,0.5); % 從s域到z域進(jìn)行雙線性變換Hhp2,w2=freqz(p6,p5);figure(3)plot(w2/2/pi*Fs,abs(Hhp2)grid onxlabel('f/Hz'),ylabel('幅度')title('用chebyshey2型設(shè)計(jì)高通IIR濾波器')程序運(yùn)行后可以分別得到用butterworth、chebyshey1型和chebyshey2型設(shè)計(jì)的低通IIR濾波器頻域特性如下圖5：圖5-1：用but

17、terworth設(shè)計(jì)高通IIR濾波器的頻域特性圖5-2：用chebyshey1型設(shè)計(jì)高通IIR濾波器的頻域特性圖5-3：用chebyshey2型設(shè)計(jì)高通IIR濾波器的頻域特性(3) 用butterworth設(shè)計(jì)帶通IIR濾波器，設(shè)計(jì)條件：通帶截止頻率fp1=1200Hz, fp2=3000Hz，阻帶截止頻率fs1=1000Hz，fs2=3200Hz,通帶波紋1dB，最小阻帶衰減60dB。設(shè)計(jì)思路是：先用butterworth設(shè)計(jì)低通的IIR濾波器，然后使用lp2bp函數(shù)完成從低通到帶通的轉(zhuǎn)換。取采樣率，編寫MATLAB程序如下：clearFs=8000;a=1;b=60;fp1=1200;fp

18、2=3000;fs1=1000;fs2=3200;wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;Wp1=tan(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);Wp=Wp1*Wp2;Ws=Ws1*Ws2;Bw=Wp2-Wp1;if Wp>Ws Ws1=Wp/Ws2;end;WWp=1;WWs=(Wp-Ws12)/(Bw*Ws1);N,Wn=buttord(WWp,WWs,a,b,'s');B A=butter(N,Wn,'s&#

19、39;);num,den=lp2bp(B,A,sqrt(Wp),Bw);p2,p1=bilinear(num,den,0.5);Hbp,w=freqz(p2,p1);figure(1)plot(w/2/pi*Fs,abs(Hbp)grid onxlabel('f/Hz'),ylabel('幅度')title('用butteworth設(shè)計(jì)帶通IIR濾波器')程序運(yùn)行后可以得到用butteworth設(shè)計(jì)帶通IIR濾波器的頻域特性如圖6所示：圖6：用butteworth設(shè)計(jì)帶通IIR濾波器的頻域特性二、實(shí)驗(yàn)內(nèi)容（二）選擇3個(gè)不同頻段的信號(hào)對(duì)其進(jìn)行頻譜分

20、析，根據(jù)信號(hào)的頻譜特征設(shè)計(jì)3個(gè)不同的數(shù)字濾波器。將三路信號(hào)合成一路信號(hào)，分析合成信號(hào)的時(shí)域和頻域特點(diǎn)，然后將合成信號(hào)分別通過設(shè)計(jì)好的3個(gè)數(shù)字濾波器，分離出原來的三路信號(hào)，分析得到的三路信號(hào)的時(shí)域波形和頻譜，與原始的三路信號(hào)進(jìn)行比較說明頻分復(fù)用的特點(diǎn)。頻分復(fù)用結(jié)構(gòu)如下圖所示：答：(1) 首先產(chǎn)生3個(gè)頻率分別為50Hz、150Hz、250Hz的正弦信號(hào)，實(shí)現(xiàn)的MATLAB程序如下：% 產(chǎn)生三個(gè)不同頻段的信號(hào)clccleart=0:.001:0.2;x1=cos(100*pi*t);x2=cos(300*pi*t);x3=cos(500*pi*t);figure(1)subplot(3,1,1)pl

21、ot(x1),axis(0,200,-1,1)xlabel('t'),ylabel('x1')title('頻率為50Hz的正弦信號(hào)')subplot(3,1,2)plot(x2),axis(0,200,-1,1)xlabel('t'),ylabel('x2')title('頻率為150Hz的正弦信號(hào)')subplot(3,1,3)plot(x3),axis(0,200,-1,1)xlabel('t'),ylabel('x3')title('頻率為250Hz的

22、正弦信號(hào)')程序運(yùn)行后可以得到3個(gè)不同頻率的正弦信號(hào)如圖7：圖7：3個(gè)不同頻段信號(hào)的時(shí)域波形 (2) 對(duì)3個(gè)不同頻段的正弦信號(hào)進(jìn)行頻域分析，實(shí)現(xiàn)的MATLAB程序如下：% 對(duì)三個(gè)不同頻率信號(hào)進(jìn)行頻譜分析N=512;X1=abs(fft(x1,N);X2=abs(fft(x2,N);X3=abs(fft(x3,N);f=1000*(0:255)/512;figure(2)subplot(3,1,1)plot(f,X1(1:256),axis(0,500,0,100)xlabel('f/Hz'),ylabel('X1')title('50Hz正弦信號(hào)

23、的幅頻特性')subplot(3,1,2)plot(f,X2(1:256),axis(0,500,0,120)xlabel('f/Hz'),ylabel('X2')title('150Hz正弦信號(hào)的幅頻特性')subplot(3,1,3)plot(f,X3(1:256),axis(0,500,0,120)xlabel('f/Hz'),ylabel('X3')title('250Hz正弦信號(hào)的幅頻特性')程序運(yùn)行后可以得到3個(gè)不同頻率的正弦信號(hào)的幅頻特性如圖8：圖8：3個(gè)不同頻段信號(hào)的幅頻特性

24、 (3) 將三路不同頻段正弦信號(hào)合成一路信號(hào)并進(jìn)行時(shí)頻分析，實(shí)現(xiàn)的MATLAB程序如下：% 將三路信號(hào)合成一路信號(hào)并進(jìn)行時(shí)頻分析x=x1+x2+x3;figure(3)subplot(2,1,1)plot(t,x)xlabel('t'),ylabel('x')title('合成信號(hào)的時(shí)域波形')X=abs(fft(x,N);subplot(2,1,2)plot(f,X(1:256),axis(0,500,0,120)xlabel('f/Hz'),ylabel('X')title('合成信號(hào)的幅頻特性'

25、;)程序運(yùn)行后可以得到合成信號(hào)的時(shí)域和頻域特性如圖9：圖9：合成信號(hào)的時(shí)域和頻域特性 (4) 設(shè)計(jì)三個(gè)數(shù)字IIR濾波器，它們分別是低通、帶通和高通。其中50Hz在低通IIR濾波器的通帶內(nèi)，150Hz、250Hz在其阻帶內(nèi)；150Hz在帶通濾波器的通帶內(nèi)，50Hz、250Hz在其阻帶內(nèi)；250Hz在高通濾波器的通帶內(nèi)，50Hz、150Hz在其阻帶內(nèi)。實(shí)現(xiàn)的MATLAB程序如下：% 設(shè)計(jì)三個(gè)數(shù)字濾波器(IIR或FIR)% 設(shè)計(jì)IIR數(shù)字低通濾波器% 設(shè)計(jì)條件Fs=600;fp=70;fs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp=tan(wp/2);Ws=tan(ws/

26、2);a=1;b=40;% 用butterworth設(shè)計(jì)低通IIR濾波器N1,Wc1=buttord(Wp,Ws,a,b,'s');num1,den1=butter(N1,Wc1,'s');p2,p1=bilinear(num1,den1,0.5);H1,w1=freqz(p2,p1);figure(4)plot(w1/2/pi*Fs,abs(H1)grid onxlabel('f/Hz'),ylabel('幅度')title('用butterworth設(shè)計(jì)低通IIR濾波器')% 設(shè)計(jì)IIR數(shù)字帶通濾波器% 設(shè)計(jì)條

27、件Fs=800;fp1=100;fp2=180;fs1=80;fs2=200;wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;a=1;b=30;Wp1=tan(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);Wp=Wp1*Wp2;Ws=Ws1*Ws2;% 用butterworth設(shè)計(jì)帶通IIR濾波器if Wp>Ws Ws1=Wp/Ws2;end;Bw=Wp2-Wp1;WWp=1;WWs=(Wp-Ws12)/(Bw*Ws1);N2,Wc2=buttord(

28、WWp,WWs,a,b,'s');B A=butter(N2,Wc2,'s');num,den=lp2bp(B,A,sqrt(Wp),Bw);p4,p3=bilinear(num,den,0.5);H2,w2=freqz(p4,p3);figure(5)plot(w2/2/pi*Fs,abs(H2)grid onxlabel('f/Hz'),ylabel('幅度')title('用butterworth設(shè)計(jì)帶通IIR濾波器')axis(0,300,0,1.4)% 設(shè)計(jì)IIR數(shù)字高通濾波器% 設(shè)計(jì)條件Fs=1000;

29、fs=200;fp=220;ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;a=1;b=60;Wp=tan(wp/2);Ws=tan(ws/2);% 用chebyshey2型設(shè)計(jì)高通IIR濾波器N3,Wc3=cheb2ord(Wp,Ws,a,b,'s');num,den=cheby2(N3,b,Wc3,'high','s');p6,p5=bilinear(num,den,0.5);H3,w3=freqz(p6,p5);figure(6)plot(w3/2/pi*Fs,abs(H3),axis(0,500,0,1.4)grid onxlab

30、el('f/Hz'),ylabel('幅度')title('用chebyshey2型設(shè)計(jì)高通IIR濾波器')程序運(yùn)行后可以得到3個(gè)IIR濾波器的頻域特性如圖10所示：圖10-1：用butterworth設(shè)計(jì)低通IIR濾波器的幅頻特性 圖10-2：用butterworth設(shè)計(jì)帶通IIR濾波器的幅頻特性 圖10-3：用chebyshey2型設(shè)計(jì)高通IIR濾波器的幅頻特性 (5) 將合成信號(hào)依次通過低通、帶通、高通IIR濾波器，濾波后得到了3路信號(hào)。實(shí)現(xiàn)的MATLAB程序如下：% 將合成信號(hào)通過三個(gè)IIR濾波器% 用低通IIR濾波器對(duì)信號(hào)進(jìn)行濾波y1=filter(p2,p1,x);Y1=abs(f