實驗5 脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器

實驗5脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器一、實驗?zāi)康?.掌握利用脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的原理及具體方法。2.加深理解數(shù)字濾波器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化。3.掌握脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用范圍。二、實驗上機環(huán)境計算機、MATLAB軟件環(huán)境三、實驗原理1、基本原理從時域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h(n)模仿模擬濾波器的單位沖激響應(yīng)ha(t),h(n)等于2、變換方法思路：（1）將Ha(s)（2）對H（3）對ha(t)（4）對h(n)進行z變換3.設(shè)計步驟(1)確定數(shù)字濾波器的性能指標(biāo)。(2)將數(shù)字濾波器頻率指標(biāo)轉(zhuǎn)換成響應(yīng)的模擬濾波器頻率指標(biāo)(3)根據(jù)指標(biāo)，，和設(shè)計模擬濾波器。(4)將展成部分分式形式。(5)把模擬極點轉(zhuǎn)換成數(shù)字極點，得到數(shù)字濾波器可見至H(z)間的變換關(guān)系為在MATLAB中有兩種方法可以實現(xiàn)實現(xiàn)上述變換。方法1：利用residue函數(shù)和residuez函數(shù)實現(xiàn)脈沖響應(yīng)不變變換法，實用方法如下：實現(xiàn)多項式形式和部分分式形式之間的轉(zhuǎn)換。實現(xiàn)多項式形式和部分分式形式之間的轉(zhuǎn)換。方法2：matlab中提供了impinvar函數(shù)采用脈沖響應(yīng)不變法實現(xiàn)模擬濾波器到數(shù)字濾波器的變換，其使用如下：[bz,az]=impinvar(b,a,fs)采用脈沖響應(yīng)不變法將模擬濾波器系統(tǒng)函數(shù)的系數(shù)向量b和a變換成為數(shù)字濾波器系統(tǒng)函數(shù)的系數(shù)向量bz和az，fs為采樣頻率（默認為1）。[bz,az]=impinvar(b,a)采樣頻率默認為1的情況下，采用脈沖響應(yīng)不變法將模擬濾波器變換為數(shù)字濾波器。四、實驗內(nèi)容1.設(shè)采樣頻率為fs=4kHz，采用脈沖響應(yīng)不變法設(shè)計一個三階巴特沃斯數(shù)字低通濾波器，其3dB截止頻率為fc=1kHz。2.設(shè)采樣頻率為fs=10kHz，設(shè)計數(shù)字低通濾波器，滿足如下指標(biāo)通帶截止頻率：fp=1kHz，通帶波動：Rp=1dB阻帶截止頻率：fst=1.5kHz，阻帶衰減：As=15dB要求分別采用巴特沃斯、切比雪夫I型、切比雪夫II型和橢圓模擬原型濾波器及脈沖響應(yīng)不變法進行設(shè)計。結(jié)合實驗結(jié)果，分別討論采用上述設(shè)計的數(shù)字濾波器是否都能滿足給定指標(biāo)要求，分析脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用范圍。五、實驗結(jié)果及分析1.實驗代碼：>>fs=4000;>>fc=1000;>>Wc=fc*2*pi;>>b=Wc^3;>>a=[12*Wc2*Wc^2Wc^3];>>[bz,az]=impinvar(b,a,fs)bz=0.00000.58130.2114az=1.0000-0.39840.2475-0.0432>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(221);plot(w/pi,abs(H));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(222);plot(w/pi,20*log10((abs(H))/max(abs(H))));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(223);plot(w/pi,angle(H)/pi);>>gridon;xlabel('\omega(\pi)'),ylabel('PhaseofH(e^j^\omega)');>>grd=grpdelay(bz,az,w);>>subplot(224);plot(w/pi,grd);>>gridon;xlabel('\omega(\pi)'),ylabel('Groupdelay');>>G=tf(bz,az,0.1,'variable','z^-1')G=0.5813z^-1+0.2114z^-2-----------------------------------------------------------1-0.3984z^-1+0.2475z^-2-0.04321z^-3波形圖如下：2.（1）巴特沃斯濾波器實驗代碼：>>fs=10000;>>fp=1000;>>fst=1500;>>Rp=1;>>As=15;>>Ap=Rp;>>Wp=fp*2*pi;>>Wst=fst*2*pi;>>N=ceil(log10((10^(Ap/10)-1)/(10^(As/10)-1))/(2*log10(Wp/Wst)));>>Wc=Wp/(10^(Ap/10)-1)^(1/2/N);>>[b,a]=butter(N,Wc,'s');>>[bz,az]=impinvar(b,a,fs)bz=-0.00000.00060.01010.01610.00410.00010az=1.0000-3.36355.0684-4.27592.1066-0.57060.0661>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(2,2,1);plot(w/pi,abs(H));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,2);plot(w/pi,20*log10((abs(H))/max(abs(H))));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,3);plot(w/pi,angle(H)/pi);>>gridon;xlabel('\omega(\pi)'),ylabel('PhaseofH(e^j^\omega)');>>grd=grpdelay(bz,az,w);>>subplot(2,2,4);plot(w/pi,grd);>>gridon;xlabel('\omega(\pi)'),ylabel('Groupdelay');>>G=tf(bz,az,0.1,'variable','z^-1')G=-5.457e-15+0.000631z^-1+0.0101z^-2+0.01614z^-3+0.004101z^-4+0.0001033z^-5----------------------------------------------------------------------------------------1-3.364z^-1+5.068z^-2-4.276z^-3+2.107z^-4-0.5706z^-5+0.06607z^-6Sampletime:0.1secondsDiscrete-timetransferfunction.波形圖：由圖線可知，結(jié)果滿足設(shè)計指標(biāo)要求。（2）切比雪夫I型實驗代碼：>>fs=10000;>>fp=1000;>>fst=1500;>>Rp=1;>>As=15;>>Ap=Rp;>>epc=sqrt(10^(Ap/10)-1);>>Wp=fp*2*pi;>>Wc=Wp;>>Wst=fst*2*pi;>>N=ceil(acosh(sqrt(10^(0.1*As)-1)/epc)/acosh(Wst/Wc));>>[b,a]=cheby1(N,Rp,Wp,'s');>>[bz,az]=impinvar(b,a,fs)bz=0.00000.00540.01810.00400az=1.0000-3.05913.8323-2.29190.5495>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(2,2,1);plot(w/pi,abs(H));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,2);plot(w/pi,20*log10((abs(H))/max(abs(H))));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,3);plot(w/pi,angle(H)/pi);>>gridon;xlabel('\omega(\pi)'),ylabel('PhaseofH(e^j^\omega)');>>grd=grpdelay(bz,az,w);>>subplot(2,2,4);plot(w/pi,grd);>>gridon;xlabel('\omega(\pi)'),ylabel('Groupdelay');>>G=tf(bz,az,0.1,'variable','z^-1')G=3.411e-17+0.005373z^-1+0.0181z^-2+0.003985z^-3-------------------------------------------------------1-3.059z^-1+3.832z^-2-2.292z^-3+0.5495z^-4Sampletime:0.1secondsDiscrete-timetransferfunction.波形圖：由圖線可知，滿足設(shè)計指標(biāo)要求。（3）切比雪夫II型實驗代碼：>>fs=10000;>>fp=1000;>>fst=1500;>>Rp=1;>>As=15;>>epc=sqrt(10^(Rp/10)-1);>>Wp=fp*2*pi;>>Wst=fst*2*pi;>>Wc=Wst;>>N=ceil(acosh(sqrt(10^(0.1*As)-1)/epc)/acosh(Wc/Wp));>>[b,a]=cheby2(N,As,Wst,'s');>>[bz,az]=impinvar(b,a,fs)bz=-0.42140.9724-0.81190.42050.0000az=1.0000-1.66581.4289-0.51930.0935>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(2,2,1);plot(w/pi,abs(H));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,2);plot(w/pi,20*log10((abs(H))/max(abs(H))));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,3);plot(w/pi,angle(H)/pi);>>gridon;xlabel('\omega(\pi)'),ylabel('PhaseofH(e^j^\omega)');>>grd=grpdelay(bz,az,w);>>subplot(2,2,4);plot(w/pi,grd);>>gridon;xlabel('\omega(\pi)'),ylabel('Groupdelay');>>G=tf(bz,az,0.1,'variable','z^-1')G=-0.4214+0.9724z^-1-0.8119z^-2+0.4205z^-3+1.663e-06z^-4------------------------------------------------------------------1-1.666z^-1+1.429z^-2-0.5193z^-3+0.09352z^-4Sampletime:0.1secondsDiscrete-timetransferfunction.波形圖:由圖線可知，阻帶衰減小于15dB，不滿足設(shè)計要求。（4）橢圓模擬原型濾波器實驗代碼：>>fs=10000;>>fp=1000;>>fst=1500;>>Rp=1;>>As=15;>>epc=sqrt(10^(Rp/10)-1);>>Wp=fp*2*pi;>>Wst=fst*2*pi;>>Wc=Wp;>>A=10^(As/20);>>k1=epc/(sqrt(A^2-1));>>k=Wp/Wst;>>N=ceil(ellipke(k)*ellipke(sqrt(1-k1^2))/ellipke(k1)/ellipke(sqrt(1-k^2)));>>[bz,az]=ellip(N,Rp,As,fp/fs*pi)bz=0.18550.04960.04960.1855az=1.0000-1.41071.2357-0.3549>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(2,2,1);plot(w/pi,abs(H));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,2);plot(w/pi,20*log10((abs(H))/max(abs(H))));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,3);plot(w/pi,angle(H)/pi);>>gridon;xlabel('\omega(\pi)'),ylabel('PhaseofH(e^j^\omega)');>>grd=grpdelay(bz,az,w);>>subplot(2,2,4);plot(w/pi,grd);>>gridon;xlabel('\omega(\pi)'),ylabel('Groupdelay');>>G=tf(bz,az,0.1,'variable','z^-1')G=0.1855+0.04956z^-1+0.04956z^-2+0.1855z^-3--------------------------------------------------1-1.411z^-1+1.236z^-2-0.3549z^-3Sampletime:0.1secondsDiscrete-timetransferfunction.波形圖：由圖線可知，滿足設(shè)計指標(biāo)要求。（5）脈沖相應(yīng)不變法：實驗代碼：>>fs=10000;>>fp=1000;>>fst=1500;>>Rp=1;>>As=15;>>Ap=Rp;>>Wp=fp*2*pi;>>Wst=fst*2*pi;>>N=ceil(log10((10^(Ap/10)-1)/(10^(As/10)-1))/(2*log10(Wp/Wst)));>>Wc=Wp/(10^(Ap/10)-1)^(1/2/N);>>b=Wc^3;>>a=[12*Wc2*Wc^2Wc^3];>>[bz,az]=impinvar(b,a,fs)bz=00.10570.0663az=1.0000-1.64911.0663-0.2450>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(2,2,1);plot(w/pi,abs(H));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,2);plot(w/pi,20*log10((abs(H))/max(abs(H))));>>gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');>>subplot(2,2,3);plot(w/pi,angle(H)/pi);>>gridon;xlabel('\omega(\pi)'),ylabel('PhaseofH(e^j^\omega)')