實(shí)驗(yàn)二 數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告IIR濾波器設(shè)計(jì)

-PAGE1-數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱：IIR濾波器設(shè)計(jì)實(shí)驗(yàn)日期：姓名：學(xué)號(hào)：哈爾濱工業(yè)大學(xué)（威海）數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告一、實(shí)驗(yàn)?zāi)康?、掌握沖激響應(yīng)法和雙線性變換法設(shè)計(jì)IIR濾波器的原理及具體設(shè)計(jì)方法，熟悉用雙線性設(shè)計(jì)法設(shè)計(jì)低通、帶通和高通IIR數(shù)字濾波器的計(jì)算機(jī)程序；2、熟悉模擬Butterworth濾波器的設(shè)計(jì)，掌握沖激響應(yīng)法和雙線性變換法設(shè)計(jì)數(shù)字IIR濾波器的方法。二、實(shí)驗(yàn)原理IIR數(shù)字濾波器的設(shè)計(jì)方法可以概括為圖8-1，本實(shí)驗(yàn)主要掌握IIR濾波器的第一種方法，即利用模擬濾波器設(shè)計(jì)IIR數(shù)字濾波器，這是IIR數(shù)字濾波器設(shè)計(jì)最常用的方法。利用模擬濾波器設(shè)計(jì)，需要將模擬域的轉(zhuǎn)換為數(shù)字域，最常用的轉(zhuǎn)換方法為脈沖響應(yīng)不變法和雙線性變換法。

圖8-1

IIR數(shù)字濾波器的設(shè)計(jì)方法圖

（1）脈沖響應(yīng)不變法

用數(shù)字濾波器的單位脈沖響應(yīng)序列模仿模擬濾波器的沖激響應(yīng)，讓正好等于的采樣值，即其中T為采樣間隔。如果以及分別表示的拉氏變換及的Z變換，則在MATLAB中，可用函數(shù)impinvar實(shí)現(xiàn)從模擬濾波器到數(shù)字濾波器的脈沖響應(yīng)不變映射。

（2）雙線性變換法s平面與z平面之間滿足下列映射關(guān)系：

或s平面的虛軸單值地映射于z平面的單位圓上，s平面的左半平面完全映射到z平面的單位圓內(nèi)。雙線性變換不存在頻率混疊問(wèn)題。

在MATLAB中，可用函數(shù)bilinear實(shí)現(xiàn)從模擬濾波器到數(shù)字濾波器的雙線性變換映射。

雙線性變換是一種非線性變換，即，這種非線性引起的幅頻特性畸變可通過(guò)預(yù)畸得到校正。

三、實(shí)驗(yàn)內(nèi)容和步驟不同階次模擬巴特沃茲濾波器的頻率響應(yīng)程序：利用Matlab系統(tǒng)函數(shù)計(jì)算不同階次的巴特沃茲濾波器頻率響應(yīng)選取階數(shù)N為3,5,10,12時(shí)clc,clearall,closeall;freq_axis=[0:0.01:2];%歸一化的頻率N_array=[351012];%濾波器的階次figure(1);forI=1:4N=N_array(I);[zpk]=buttap(N);%獲得5階巴特沃茲濾波器的參數(shù)[ba]=zp2tf(z,p,k);%得到傳遞函數(shù)tf:trans-functiondisp(['N='num2str(N)'極點(diǎn):']);pdisp(['N='num2str(N),'零點(diǎn):']);z[H]=freqs(b,a,freq_axis);magh=abs(H);subplot(2,2,I);plot(freq_axis,magh);axis([0201]);%限制坐標(biāo)軸xlabel('w/wc--歸一化頻率');%標(biāo)注ylabel('幅度譜');title(['巴特沃茲濾波器幅度譜N=',num2str(N)]);end圖形：結(jié)論：由上圖對(duì)比可知巴特沃茲濾波器在通帶和阻帶內(nèi)都很平坦，起伏不大。巴特沃茲濾波器的階數(shù)越大，過(guò)渡帶越窄。2、根據(jù)模擬濾波器指標(biāo)，設(shè)計(jì)低通巴特沃茲濾波器要求通帶截止頻率fp=5kHz,通帶最大衰減ɑp=2dB，阻帶截止頻率fs=12kHz,阻帶最小衰減ɑs=30dB,按照上述指標(biāo)設(shè)計(jì)巴特沃茲低通濾波器。程序：J=sqrt(-1);fp=5e3;%e科學(xué)計(jì)數(shù)法，代表10的幾次方，5KHz：通帶截止頻率alfa_p=2;%最大通帶衰減2dBfs=12e3;%12Kz阻帶截止頻率alfa_s=30;%30dB阻帶最小衰減k_sp=sqrt(10^(alfa_p/10)-1)/sqrt(10^(alfa_s/10)-1);lamda_sp=2*pi*fs/(2*pi*fp);N=ceil(-log(k_sp)/log(lamda_sp));%ceil：加一個(gè)正小數(shù)取整p=zeros(1,N);forI=1:Nk=I-1;p(I)=exp(J*pi/2+J*pi*(2*k+1)/2/N);%e指數(shù)endz=[];%表示無(wú)零點(diǎn)[b1a1]=zp2tf(z,p,1);%得到傳遞函數(shù)tf:trans-functionfc_freq=fp*(10^(alfa_p/10)-1)^(-1/2/N);%3dB截至頻率figure(2);freq_axis1=[0:50:15e3];[H1]=freqs(b1,a1,freq_axis1/fc_freq);%對(duì)Wc作頻率歸一化plot(freq_axis1,20*log10(abs(H1)));H1=freqs(b1,a1,freq_axis);figure(3);plot(freq_axis,abs(H1));holdon;plot(freq_axis,abs(H),'r.')圖形：結(jié)論N為5階的巴特沃茲低通濾波器。3、給定模擬濾波器，按照不同采樣頻率設(shè)計(jì)數(shù)字IIR，觀察采樣頻率對(duì)設(shè)計(jì)結(jié)果的影響已知模擬濾波器傳輸函數(shù)H(s)為按照不同采樣頻率Fs1=1Hz,Fs2=10Hz使用脈沖響應(yīng)不變法將H(s)轉(zhuǎn)換為H(z)數(shù)字IIR濾波器。觀察不同采樣頻率設(shè)計(jì)得到的數(shù)字IIR的差異。程序：J=sqrt(-1);b=[000.5012];%高次項(xiàng)系數(shù)在前a=[10.64490.7079];freq_analog=[0:0.5:30];%模擬角頻率!!!H_analog=freqs(b,a,freq_analog)；%計(jì)算在模擬角頻率處的%頻響figure(1);subplot(2,1,1);plot(freq_analog,20*log10(abs(H_analog)));xlabel('角頻率(rad/s)');ylabel('幅頻響應(yīng)(dB)');title('模擬濾波器頻率響應(yīng)');Fs1=1;%1Hz采樣率[bz1,az1]=impinvar(b,a,Fs1);%根據(jù)模擬濾波器系數(shù)b,a，%按照%采樣率fs轉(zhuǎn)換為系數(shù)為bz,az的數(shù)字IIR濾波器Fs2=10;%10Hz采樣頻率[bz2,az2]=impinvar(b,a,Fs2);%根據(jù)模擬濾波器系數(shù)b,a，%按照%采樣率fs轉(zhuǎn)換為系數(shù)為bz,az的數(shù)字IIR濾波器freq_digital=[0:0.1:2*pi];H_dig1=freqz(bz1,az1,freq_digital);H_dig2=freqz(bz2,az2,freq_digital);subplot(2,1,2);plot(freq_digital,20*log10(abs(H_dig1)));xlabel('數(shù)字頻率');ylabel('幅頻響應(yīng)(dB)');title('不同采樣頻率對(duì)于數(shù)字IIR頻響的影響');holdon;plot(freq_digital,20*log10(abs(H_dig2)),'r');legend(['Fs1='num2str(Fs1)'Hz'],['Fs2='num2str(Fs2)'Hz']);圖形：結(jié)論：不同采樣頻率設(shè)計(jì)得到的數(shù)字IIR的差異。實(shí)際系統(tǒng)不可能是嚴(yán)格帶限的，都會(huì)混疊失真，由上圖可知采樣頻率越大，混疊失真越小。4、給定數(shù)字IIR指標(biāo)，用脈沖響應(yīng)不變法設(shè)計(jì)IIR。數(shù)字低通濾波器的指標(biāo)為:,使用巴特沃茲濾波器設(shè)計(jì)數(shù)字IIR。比較采樣頻率Fs=1，2時(shí)設(shè)計(jì)過(guò)程中的模擬濾波器，數(shù)字濾波器的變化。程序：J=sqrt(-1);freq_p=0.2*pi;%通帶截止頻率0.2pialpha_p=1;%通帶最大衰減1dBfreq_s=0.3*pi;alpha_s=15;Fs_array=[12];Color_array=['b','k'];forIII=1:2color_used=Color_array(III);Fs=Fs_array(III);T=1/Fs;analog_freq_p=freq_p/T%數(shù)字頻率w到模擬角頻率W的轉(zhuǎn)換W=w/T;analog_freq_s=freq_s/T%根據(jù)模擬濾波器的指標(biāo)求巴特沃茲濾波器的階次k_sp=sqrt(10^(alpha_p/10)-1)/sqrt(10^(alpha_s/10)-1);lamda_sp=2*pi*analog_freq_s/(2*pi*analog_freq_p);N=ceil(-log(k_sp)/log(lamda_sp));Wc=analog_freq_p*(10^(analog_freq_p/10)-1)^(-1/2/N);%3dB截至頻率%查表，得到歸一化的傳輸函數(shù)或者根據(jù)公式得到極點(diǎn)pp=zeros(1,N);forI=1:Nk=I-1;p(I)=exp(J*pi/2+J*pi*(2*k+1)/2/N);end%根據(jù)得到的零點(diǎn)、極點(diǎn)得到歸一化的傳輸函數(shù)H(p)z=[];%表示無(wú)零點(diǎn)[ba]=zp2tf(z,p,1);%得到傳遞函數(shù)tf:trans-function,第三個(gè)參數(shù)為增益=1%根據(jù)p=s/Wc把歸一化傳輸函數(shù)H(p)的系數(shù)轉(zhuǎn)換為實(shí)際的模擬傳輸函數(shù)H(s)bs=zeros(1,length(b));%H(s)的分子系數(shù)as=zeros(1,length(a));%對(duì)與N階模擬濾波器，設(shè)其傳輸函數(shù)分子分母系數(shù)分別為bs,astemp=[N:-1:0];temp=Wc.^temp;bs=b./temp;as=a./temp;temp=as(1);as=as/temp;%驗(yàn)證與書上結(jié)果是否一致bs=bs/temp;figure(1);freq_analog=linspace(0,2*analog_freq_s);H_analog=freqs(bs,as,freq_analog);plot(freq_analog,20*log10(abs(H_analog)),color_used);gridon;holdon;xlabel('模擬角頻率(rad/s)');ylabel('幅頻響應(yīng)(dB)');title('模擬濾波器的頻率響應(yīng)');%根據(jù)H(s)獲得H（z），使用脈沖響應(yīng)不變法[bz,az]=impinvar(bs,as,Fs);%根據(jù)模擬濾波器系數(shù)b,a，按照采樣率fs轉(zhuǎn)換為系數(shù)為bz,az的數(shù)字IIR濾波器freq_digital=[0:0.1:2*pi-0.1];H=freqz(bz,az,freq_digital);figure(2);plot(freq_digital/pi,20*log10(abs(H)),color_used);gridon;holdon;xlabel('數(shù)字頻率/pi');ylabel('幅頻響應(yīng)(dB)');title('數(shù)字IIR濾波器頻率響應(yīng)');a=1;end%%%%%%%%%%%%%%%%%%%%%%%%%%%當(dāng)采樣頻率變?yōu)?Hz時(shí)，求得出的數(shù)字IIR頻率響應(yīng)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(1),legend(['Fs='num2str(Fs_array(1))'Hz'],['Fs='num2str(Fs_array(2))'Hz']);figure(2),legend(['Fs='num2str(Fs_array(1))'Hz'],['Fs='num2str(Fs_array(2))'Hz']);圖形：結(jié)論：當(dāng)濾波器設(shè)計(jì)指標(biāo)以數(shù)字域頻率Wc給定時(shí)，采樣頻率的變化會(huì)導(dǎo)致模擬頻率改變，由圖可知fs增大會(huì)導(dǎo)致模擬截止頻率增大。所以當(dāng)濾波器設(shè)計(jì)指標(biāo)以數(shù)字域頻率Wc給定時(shí)，不能通過(guò)提高抽樣頻率來(lái)改善混疊現(xiàn)象。5、用雙線性變換法設(shè)計(jì)IIR，并同脈沖響應(yīng)不變法比較設(shè)計(jì)的結(jié)果。數(shù)字濾波器指標(biāo)要求如下程序：clc,clearall,closeall;%給定數(shù)字濾波器指標(biāo)J=sqrt(-1);freq_p=0.2*pi;%通帶截止頻率0.2pialpha_p=1;%通帶最大衰減1dBfreq_s=0.3*pi;alpha_s=15;%模擬低通指標(biāo)Fs=1;T=1/Fs;analog_freq_p=2/T*tan(freq_p/2)analog_freq_s=2/T*tan(freq_s/2)%計(jì)算巴特沃茲濾波器的指標(biāo)k_sp=sqrt(10^(alpha_p/10)-1)/sqrt(10^(alpha_s/10)-1);lamda_sp=2*pi*analog_freq_s/(2*pi*analog_freq_p);N=ceil(-log(k_sp)/log(lamda_sp));Wc=analog_freq_p*(10^(analog_freq_p/10)-1)^(-1/2/N);%3dB截至頻率%查表，得到歸一化的傳輸函數(shù)或者根據(jù)公式得到極點(diǎn)pp=zeros(1,N);forI=1:Nk=I-1;p(I)=exp(J*pi/2+J*pi*(2*k+1)/2/N);end%根據(jù)得到的零點(diǎn)、極點(diǎn)得到歸一化的傳輸函數(shù)H(p)z=[];%表示無(wú)零點(diǎn)[ba]=zp2tf(z,p,1);%得到傳遞函數(shù)tf:trans-function,第三個(gè)參數(shù)為增益=1%根據(jù)p=s/Wc把歸一化傳輸函數(shù)H(p)的系數(shù)轉(zhuǎn)換為實(shí)際的模擬傳輸函數(shù)H(s)bs=zeros(1,length(b));%H(s)的分子系數(shù)as=zeros(1,length(a));%對(duì)與N階模擬濾波器，設(shè)其傳輸函數(shù)分子分母系數(shù)分別為bs,astemp=[N:-1:0];temp=Wc.^temp;bs=b./tem