實(shí)驗(yàn)三-用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器

實(shí)驗(yàn)三用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器實(shí)驗(yàn)工程名稱：用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器實(shí)驗(yàn)工程性質(zhì)：驗(yàn)證性實(shí)驗(yàn)所屬課程名稱：數(shù)字信號(hào)處理實(shí)驗(yàn)方案學(xué)時(shí)：2實(shí)驗(yàn)?zāi)康氖煜び秒p線性變換法設(shè)計(jì)IIR數(shù)字濾波器的原理與方法。掌握數(shù)字濾波器的計(jì)算機(jī)仿真方法。通過(guò)觀察對(duì)實(shí)際心電圖信號(hào)的濾波作用，獲得數(shù)字濾波的感性知識(shí)。實(shí)驗(yàn)內(nèi)容和要求用雙線性變換法設(shè)計(jì)一個(gè)巴特沃斯低通IIR數(shù)字濾波器。設(shè)計(jì)指標(biāo)參數(shù)為：在通帶內(nèi)頻率低于0.2π時(shí)，最大衰減小于1dB；在阻帶內(nèi)[0.3π，π]頻率區(qū)間上，最小衰減大與15dB。以0.02π為采樣間隔，打印出數(shù)字濾波器在頻率區(qū)間[0，π/2]的幅頻響應(yīng)特性曲線。用所設(shè)計(jì)的濾波器對(duì)實(shí)際心電圖信號(hào)采樣序列〔在本實(shí)驗(yàn)后面給出〕進(jìn)行仿真濾波處理，并打印出濾波前后的心電圖信號(hào)波形圖，觀察總結(jié)濾波作用與效果。實(shí)驗(yàn)主要儀器設(shè)備和材料計(jì)算機(jī)，MATLAB6.5或以上版本實(shí)驗(yàn)方法、步驟及結(jié)果測(cè)試復(fù)習(xí)有關(guān)巴特沃斯模擬濾波器設(shè)計(jì)和用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的內(nèi)容，用雙線性變換法設(shè)計(jì)數(shù)字濾波器系統(tǒng)函數(shù)。其中滿足本實(shí)驗(yàn)要求的數(shù)字濾波器系統(tǒng)函數(shù)為： 〔3.1〕式中： 〔3.2〕根據(jù)設(shè)計(jì)指標(biāo)，調(diào)用MATLAB信號(hào)處理工具箱buttord和butter，也可以得到。圖3-1濾波器的組成由公式〔3.1〕和〔3.2〕可見(jiàn)，濾波器由三個(gè)二階濾波器、和級(jí)聯(lián)而成，如圖3-1所示。圖3-1濾波器的組成編寫(xiě)濾波器仿真程序，計(jì)算對(duì)心電圖采樣序列x(n)的響應(yīng)序列y(n)。設(shè)yk(n)為第k級(jí)二階濾波器Hk(z)的輸出序列，yk-1(n)為輸入序列，如圖3-1所示。由〔3.2〕式可得到差分方程為：〔3.3〕當(dāng)k=1時(shí)，yk-1(n)=x(n)。所以H(z)對(duì)x(n)的總響應(yīng)序列y(n)可以用順序迭代算法得到。即依次對(duì)k=1，2，3，求解差分方程〔3.3〕，最后得到y(tǒng)3(n)=y(n)。仿真程序就是實(shí)現(xiàn)上述求解差分方程和順序迭代算法的通用程序。也可以直接調(diào)用MATLABfilter函數(shù)實(shí)現(xiàn)仿真。在通用計(jì)算機(jī)上運(yùn)行仿真濾波程序，并調(diào)用通用繪圖子程序，完成實(shí)驗(yàn)內(nèi)容〔2〕和〔3〕。本實(shí)驗(yàn)中用到的心電圖信號(hào)采用序列x(n)人體心電圖信號(hào)在測(cè)量過(guò)程往往受到工業(yè)高頻干擾，所以必須經(jīng)過(guò)低通濾波處理后，才能作為判斷心臟功能的有用信息。下面給出一實(shí)際心電圖信號(hào)采樣序列樣本x(n)，其中存在高頻干擾。在實(shí)驗(yàn)中，以x(n)作為輸入序列，濾除其中的干擾成分。{x(n)}={-4，-2，0，-4，-6，-4，-2，-4，-6，-6， -4，-4，-6，-6，-2，6，12，8，0，-16， -38，-60，-84，-90，-66，-32，-4，-2，-4，8， 12，12，10，6，6，6，4，0，0，0， 0，0，-2，-4，0，0，0，-2，-2，0，0，-2，-2，-2，-2，0}T=1;Fs=1/T;wpz=0.2;wsz=0.3;wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=15;[N,wc]=buttord(wp,ws,rp,rs,'s');[B,A]=butter(N,wc,'s');fk=0:1/512:1;wk=2*pi*fk;Hk=freqs(B,A,wk);subplot(2,2,1);plot(fk,20*log10(abs(Hk)));gridon;xlabel('\omega/\pi');ylabel('幅度(dB)');axis([0,1,-100,5]);title('(b)');[N,wc]=buttord(wpz,wsz,rp,rs);[Bz,Az]=butter(N,wc);wk=0:pi/512:pi;Hz=freqz(Bz,Az,wk);subplot(2,2,4);plot(wk/pi,20*log10(abs(Hz)));gridon;xlabel('\omega/\pi');ylabel('幅度(dB)');axis([0,1,-100,5]);title('(b)');二；x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];subplot(2,2,1);n=0:55;stem(n,x,'.');xlabel('n');ylabel('x(n)');title('x(n)的脈沖響應(yīng)');A=0.09036;b1=[A,2*A,A];a1=[1,-1.2686,0.7051];h1=filter(b1,a1,x);[H1,w]=freqz(b1,a1,100);b2=[A,2*A,A];a2=[1,-1.0106,0.3583];h2=filter(b2,a2,h1);[H2,w]=freqz(b2,a2,100);b3=[A,2*A,A];a3=[1,-0.9044,0.2155];h3=filter(b3,a3,h2);[H3,w]=freqz(b3,a3,100);subplot(2,2,2);stem(n,h3,'.');xlabel('n');ylabel('y(n)');title('通過(guò)濾波器H1(z),H2(z),H3(z)后的y3(n)函數(shù)');subplot(2,2,3);H4=H1.*(H2);H=H4.*(H3);mag=abs(H);db=20*log10((mag+eps)/max(mag));plot(w/pi,db);xlabel('ω/π');ylabel('20log[Ha3(ejw)]');title('通過(guò)濾波器H1(z),H2(z),H3(z)后的對(duì)數(shù)頻率響應(yīng)20log[Ha3(ejw)]函數(shù)');grid;figure(2);N=1024;n=0:N/2-1;Xk=fft(x,N);AXk=abs(Xk(1:N/2));f=(0:N/2-1)*Fs/N;f=f/Fs;subplot(211);plot(f,AXk);title('x(n)的頻譜');xlabel('f');ylabel('|X(k)|');axis([0,0.5,0,400]);Yk=fft(y,N);AYk=abs(Yk(1:N/2));subplot(212);plot(f,AYk);title('y(n)的頻譜');xlabel('f');ylabel('|Y(k)|');axis([0,0.5,0,400])%x(n)的心電脈沖函數(shù)x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]subplot(2,2,1);n=0:55;stem(n,x,'.');xlabel('n');ylabel('x(n)');title('x(n)的心電脈沖函數(shù)');%通過(guò)濾波器H1(z)后的y1(n)函數(shù)A=0.09036;b1=[A,2*A,A];a1=[1,-1.2686,0.7051];h1=filter(b1,a1,x);[H1,w]=freqz(b1,a1,100);%通過(guò)濾波器H1(z),H2(Z)后的y2(n)函數(shù)b2=[A,2*A,A];a2=[1,-1.0106,0.3583];h2=filter(b2,a2,h1);[H2,w]=freqz(b2,a2,100);%通過(guò)濾波器H1(z),H2(Z),H3(Z)后的y3(n)函數(shù)b3=[A,2*A,A];a3=[1,-0.9044,0.2155];h3=filter(b3,a3,h2);[H3,w]=freqz(b3,a3,100);subplot(2,2,2);stem(n,h3,'.');xlabel('n');ylabel('y(n)');title('通過(guò)濾波器H1(Z),H2(Z),H3(Z)后的y3(n)函數(shù)');subplot(2,2,3);H4=H1.*(H2);H=H4.*(H3);mag=abs(H);db=20*log10((mag+eps)/max(mag));plot(w/pi,db);xlabel('w/pi');ylabel('20log[Ha3(ejw)]');title('通過(guò)濾波器H1(z)、H2(z)、H3(z)后的對(duì)數(shù)頻率響應(yīng)20log[Ha3(ejw)]函數(shù)');grid;MATLAB程序清單：%實(shí)驗(yàn)三，用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];k=1;closeall;figure(1);subplot(2,2,1);n=0:55;%更正stem(n,x,'.');axis([0,56,-100,50]);%更正holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('心電圖信號(hào)采集序列x(n)');B=[0.09036,2*0.09036,0.09036];A=[1.2686,-0.7051];A1=[1.0106,-0.3583];A2=[0.9044,-0.2155];while(k<=3)y=filter(B,A,x);%Thefunctionistofiltethesingalxx=y;ifk==2A=A1;endifk==3A=A2;endk=k+1;endsubplot(2,2,3)n=0:55;%更正stem(n,y,'.');axis([0,56,-15,5]);holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('y(n)');title('三級(jí)濾波后的心電圖信號(hào)');%求數(shù)字濾波器的幅頻特性A=[0.09036,0.1872,0.09036];B1=[1,-1.2686,0.7051];B2=[1,-1.0106,0.3583];B3=[1,-0.9044,0.2155];[H1,w]=freqz(A,B1,100);[H2,w]=freqz(A,B2,100);[H3,w]=freqz(A,B3,100);H4=H1.*(H2);H=H4.*(H3);mag=abs(H);db=20*log10((mag+eps)/max(mag));subplot(2,2,2)plot(w/pi,db);axis