數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告4

1、 實(shí)驗(yàn)報(bào)告學(xué)院（系）名稱：計(jì)算機(jī)與通信工程學(xué)院姓名學(xué)號(hào)專業(yè)班級(jí)實(shí)驗(yàn)項(xiàng)目實(shí)驗(yàn)四 IIR數(shù)字濾波器的設(shè)計(jì)課程名稱數(shù)字信號(hào)處理課程代碼實(shí)驗(yàn)時(shí)間2013年05月31日實(shí)驗(yàn)地點(diǎn)主校區(qū)計(jì)算機(jī)基礎(chǔ)實(shí)驗(yàn)室批改意見(jiàn)成績(jī)教師簽字：一，實(shí)驗(yàn)?zāi)康?加深理解IIR數(shù)字濾波器的時(shí)域特性和頻率特性，掌握IIR數(shù)字濾波器的設(shè)計(jì)原理與設(shè)計(jì)方法，以及IIR數(shù)字濾波器的應(yīng)用。二，實(shí)驗(yàn)原理IIR數(shù)字濾波器一般為香型移不變的因果離散系統(tǒng)，N階IIR數(shù)字濾波器的系統(tǒng)函數(shù)可以表達(dá)為的有理多項(xiàng)式形式，即 式中：系數(shù)至少有一個(gè)非零。對(duì)于因果IIR數(shù)字濾波器，應(yīng)滿足。IIR數(shù)字濾波器的設(shè)計(jì)主要通過(guò)成熟的模擬濾波器設(shè)計(jì)方法來(lái)實(shí)現(xiàn)。首先在頻域?qū)?shù)字

2、濾波器設(shè)計(jì)指標(biāo)轉(zhuǎn)換為模擬濾波器設(shè)計(jì)指標(biāo)，然后將任意的模擬濾波器轉(zhuǎn)換為原型模擬低通濾波器指標(biāo)，根據(jù)模擬低通濾波器的設(shè)計(jì)指標(biāo)來(lái)設(shè)計(jì)出模擬低通濾波器，然后由經(jīng)過(guò)相應(yīng)的復(fù)頻移轉(zhuǎn)換得到 ，最后由經(jīng)過(guò)脈沖響應(yīng)不變法或雙線性變換法得到所需要的IIR數(shù)字濾波器。由此可見(jiàn)，IIR數(shù)字濾波器設(shè)計(jì)的重要環(huán)節(jié)是模擬低通濾波器的設(shè)計(jì)。設(shè)計(jì)模擬低通濾波器的主要方法有Butterworth、Chebyshev和橢圓等濾波器設(shè)計(jì)方法。MATLAB信號(hào)處理工具箱中提供了IIR濾波器設(shè)計(jì)的函數(shù)。下面分別來(lái)介紹這些常用的設(shè)計(jì)函數(shù)。1、 Butterworth數(shù)字濾波器（1）Butterworth濾波器是通帶阻都單調(diào)衰減的濾波器。調(diào)

3、用buttord函數(shù)可以確定巴特沃斯濾波器的階數(shù)，其格式為：N，Omegac=buttord(Omegap,Omegas,Rp,As,s)。其中，輸入?yún)?shù)Rp,As分別為通帶最大衰減和帶阻最小衰減，以dB為單位；Omegap、Omegas分別為通帶截止頻率和阻帶截止頻率，s說(shuō)明所設(shè)計(jì)的是模擬濾波器。輸出參數(shù)N為濾波器的階數(shù)，Omegac為3dB截止頻率。（2）調(diào)用歸一化巴特沃斯模擬原型濾波器的函數(shù)，其格式為:z0,p0,k0=buttap(N)其中，只要輸入巴特沃斯了不起的階數(shù)，它就可以返回零點(diǎn)和極點(diǎn)數(shù)組z0、p0以及增益k0。當(dāng)需要具有任意Omegac的未歸一化的巴特沃斯濾波器時(shí)，就要用Om

4、egac乘以p0或k0來(lái)進(jìn)行歸一化。（3）調(diào)用脈沖響應(yīng)不變法或雙線性變換法來(lái)設(shè)計(jì)巴特沃斯數(shù)字濾波器的函數(shù)，其格式分別如下。 脈沖響應(yīng)不變法：bd,ad=impinvar(b,a,Fs)。其中：b為模擬濾波器分子系數(shù)向量；a為模擬濾波器分母系數(shù)向量；Fs為采樣頻率；bd為數(shù)字濾波器的分子多項(xiàng)式系數(shù)；ad為數(shù)字濾波器的分母多項(xiàng)式系數(shù)。雙線性變換法：bd,ad=bilinear(b,a,Fs),參數(shù)含義與impinvar一致。2、Chebyshev型數(shù)字濾波器設(shè)計(jì)Chebyshev型濾波器為通帶紋波控制器：在通帶呈現(xiàn)紋波特性，在阻帶單調(diào)衰減。其格式如下：N,Omegac=cheb1ord(Omega

5、p,Omegas,Pp,As,s)z0,p0,k0=cheb1ap(N,Rp)參數(shù)含義與bueeord buttap中參數(shù)一致。3、Chebyshev型濾波器設(shè)計(jì)Chebyshev型濾波器為阻帶紋波濾波控制器：在阻帶呈現(xiàn)紋波特性，在通帶單調(diào)衰減。其格式如下：N,Omegac=cheb2ord(Omegap,Omegas,Rp,As,s)z0,p0,k0=cheb2ap(N,As)其參數(shù)含義同上。4、橢圓濾波器設(shè)計(jì) 橢圓濾波器在通帶和阻帶都呈現(xiàn)紋波特性。N,Omegac=ellipord(Omegap,Omegas,Rp,As.s)z0,p0,k0=ellipap(N,Rp,As)其參數(shù)含義同上

6、。三，實(shí)驗(yàn)內(nèi)容與方法【例3-4-1】利用脈沖響應(yīng)不變法，用巴特沃斯原型設(shè)計(jì)一個(gè)數(shù)字低通濾波器，其中濾波器的技術(shù)指標(biāo)如下：、，濾波器采樣頻率為1kHz。MATLAB程序如下：wp=0.4*pi;ws=0.6*pi;Rp=0.5;As=50;Fs=1000;OmegaP=wp*Fs;OmegaS=ws*Fs;ep=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);N,OmegaC=buttord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=buttap(N);p=p0*OmegaC;z=z0*OmegaC;k

7、=k0*OmegaCN;ba=k*real(poly(z);aa=real(poly(p);bd,ad = impinvar(ba,aa, Fs);H,w=freqz(bd,ad,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(bd,ad,w);subplot(2,2,1);plot(w/pi,mag);title(幅度響應(yīng))xlabel();ylabel(|H|);axis(0,1,0,1.1)set(gca,XTickMode

8、,manual,XTick,0,0.4,0.6,1);grid on;set(gca,YTickmode,manual,YTick,0,Attn,Ripple,1);subplot(2,2,3);plot(w/pi,db);title(幅度（單位：dB）);xlabel(頻率（單位：pi）);ylabel(分貝);axis(0,1,-150,50);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-150,-50,0,50);grid on;subplot(2,2,2);plot(w/pi,p

9、ha/pi);title(相位響應(yīng));xlabel();ylabel(單位：pi);axis(0,1,-1,1)set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-1,0,1);grid on;subplot(2,2,4);plot(w/pi,grd);title(群延遲 );xlabel(頻率（單位：pi）);ylabel(樣本);axis(0,1,0,20);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manua

10、l,YTick,0:2:20);grid on;set(gcf,color,w);程序運(yùn)行結(jié)果如圖所示。【例3-4-2】按照例3-4-1的技術(shù)指標(biāo)，即、，采樣頻率為1kHz。利用脈沖響應(yīng)不變法，要求用切比雪夫型來(lái)實(shí)現(xiàn)數(shù)字低通濾波器。 MATLAB程序如下：wp=0.4*pi;ws=0.6*pi;Rp=0.5;As=50;Fs=1000;OmegaP=wp*Fs;OmegaS=ws*Fs;ep=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);N,OmegaC=cheb1ord(OmegaP,OmegaS,Rp,As,s);z

11、0,p0,k0=cheb1ap(N,Rp);p=p0*OmegaC;z=z0*OmegaC;k=k0*OmegaCN;ba=k*real(poly(z);aa=real(poly(p) bd,ad = impinvar(ba,aa, Fs)H,w=freqz(bd,ad,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(bd,ad,w);subplot(2,2,1);plot(w/pi,mag);title(幅度響應(yīng))ylabel(

12、|H|);axis(0,1,0,1.1);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on;set(gca,YTickmode,manual,YTick,0,Attn,Ripple,1);subplot(2,2,3);plot(w/pi,db);title(幅度（單位：dB）);xlabel(頻率（單位：pi）);ylabel(分貝);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-100,-50,0,50);grid on;%set(g

13、ca,YTickLabelMode,manual,YTickLabels,-100;-50;0;50);subplot(2,2,2);plot(w/pi,pha/pi);title(相位響應(yīng));xlabel();ylabel(單位：pi );axis(0,1,-1,1)set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-1,0,1);grid on;subplot(2,2,4);plot(w/pi,grd);title(群延遲);xlabel(頻率（單位：pi）);ylabel(樣本);axis

14、(0,1,0,20);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,0:2:20);grid on;set(gcf,color,w);程序運(yùn)行結(jié)果如圖所示?！纠?-4-3】按照例3-4-1的技術(shù)指標(biāo)，即、，采樣頻率為1kHz。利用雙線性變換法，并要求用巴特沃斯來(lái)實(shí)現(xiàn)數(shù)字低通濾波器。 MATLAB程序如下：wp=0.4*pi;ws=0.6*pi;Rp=0.5;As=50;Fs=1000;T=1/Fs;OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);e

15、p=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);N,OmegaC=buttord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=buttap(N);p=p0*OmegaC;z=z0*OmegaC;k=k0*OmegaCN;ba0=real(poly(z0);ba0=k0*ba0;aa0=real(poly(p0);ba=real(poly(z);ba=k*ba;aa=real(poly(p) bd,ad = bilinear(ba,aa, Fs);bd1,ad1=bilinear(ba0,aa0,Fs

16、/OmegaC);figure(1);subplot(1,1,1)H,w=freqz(bd,ad,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(bd,ad,w);subplot(2,2,1);plot(w/pi,mag);title(幅度響應(yīng))ylabel(|H|);axis(0,1,0,1.1);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on;set(gca,YTickm

17、ode,manual,YTick,0,Attn,Ripple,1);subplot(2,2,3);plot(w/pi,db);title(幅度（單位：dB）);xlabel(頻率（單位：pi）);ylabel(分貝);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-100,-50,0,50);grid on;subplot(2,2,2);plot(w/pi,pha/pi);title(相位響應(yīng));xlabel();ylabel(單位：pi );axis(0,1,-1,1)set(gca,XT

18、ickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-1,0,1);grid on;subplot(2,2,4);plot(w/pi,grd);title(群延遲);xlabel(頻率（單位：pi）);ylabel(樣本);axis(0,1,0,20);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,0:2:20);grid on;set(gcf,color,w);程序運(yùn)行結(jié)果如圖所示?！纠?-4-4】按照例3-4-

19、1的技術(shù)指標(biāo)， 、，采樣頻率為1kHz。利用雙線性變化法，并要求用切比雪夫型來(lái)實(shí)現(xiàn)數(shù)字低通濾波器。MATLAB程序如下：wp=0.4*pi;ws=0.6*pi;Rp=0.5;As=50;T=0.001;Fs=1/T;OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);ep=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);N,OmegaC=cheb1ord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=cheb1ap(N,Rp);ba0=k0*real(poly(z0

20、);aa0=real(poly(p0); bd,ad = bilinear(ba0,aa0, Fs/OmegaC);sos,G=tf2sos(bd,ad);figure(1);subplot(1,1,1)H,w=freqz(bd,ad,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(bd,ad,w);subplot(2,2,1);plot(w/pi,mag);title(幅度響應(yīng))ylabel(|H|);axis(0,1,0,1.

21、1);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on;set(gca,YTickmode,manual,YTick,0,Attn,Ripple,1);subplot(2,2,3);plot(w/pi,db);title(幅度（單位：dB）);xlabel(頻率（單位：pi）);ylabel(分貝);axis(0,1,-100,10);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-100,-50,0,50);grid on;subplo

22、t(2,2,2);plot(w/pi,pha/pi);title(相位響應(yīng));xlabel();ylabel(單位：pi );axis(0,1,-1,1)set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-1,0,1);grid on;subplot(2,2,4);plot(w/pi,grd);title(群延遲);xlabel(頻率（單位：pi）);ylabel(樣本);axis(0,1,0,20);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set

23、(gca,YTickmode,manual,YTick,0:2:20);grid on;set(gcf,color,w);程序運(yùn)行結(jié)果如圖所示。四、程序設(shè)計(jì)實(shí)驗(yàn)利用切比雪夫型設(shè)計(jì)一個(gè)數(shù)字低通濾波器，濾波器的技術(shù)指標(biāo)如下：、，濾波器的采樣頻率為1kHz。MATLAB程序如下：wp=0.6*pi;ws=0.8*pi;Rp=0.3;As=55;T=0.001;Fs=1/T;OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);ep=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);N,Omeg

24、aC=cheb2ord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=cheb2ap(N,As);ba0=k0*real(poly(z0);aa0=real(poly(p0); bd,ad = bilinear(ba0,aa0, Fs/OmegaC);sos,G=tf2sos(bd,ad);figure(1);subplot(1,1,1)H,w=freqz(bd,ad,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdela

25、y(bd,ad,w);subplot(2,2,1);plot(w/pi,mag);title(幅度響應(yīng))ylabel(|H|);axis(0,1,0,1.1);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on;set(gca,YTickmode,manual,YTick,0,Attn,Ripple,1);subplot(2,2,3);plot(w/pi,db);title(幅度（單位：dB）);xlabel(頻率（單位：pi）);ylabel(分貝);axis(0,1,-100,10);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);set(gca,YTickmode,manual,YTick,-100,-50,0,50);grid on