DSP實(shí)驗(yàn)5 IIR濾波器的設(shè)計(jì)與信號(hào)濾波

1、實(shí)驗(yàn)五 IIR濾波器的設(shè)計(jì)與信號(hào)濾波一、實(shí)驗(yàn)?zāi)康脑O(shè)計(jì)IIR濾波器，實(shí)現(xiàn)對(duì)存在加性干擾的時(shí)域離散信號(hào)進(jìn)行濾波。二、實(shí)驗(yàn)內(nèi)容已知帶加性干擾的信號(hào)用x(n)表示，x(n)=xs(n)+(n)，式中xs(n)是有用的信號(hào)，是一個(gè)00.2rad的帶限信號(hào)。(n)是一個(gè)干擾信號(hào)，其頻譜分布在0.3rad以上。要求設(shè)計(jì)一個(gè)巴特沃斯IIR數(shù)字濾波器對(duì)信號(hào)x(n)進(jìn)行濾波，將干擾(n)濾除。要求在xs(n) 所在的通帶內(nèi)濾波器幅度平坦，在0.2rad處幅度衰減不大于1dB，在噪聲所在的0.3rad 以上的頻帶內(nèi)濾波器幅度衰減大于等于40 dB。三、實(shí)驗(yàn)步驟1. 根據(jù)題目要求確定要設(shè)計(jì)的數(shù)字濾波器的技術(shù)指標(biāo) (低

2、通濾波器指標(biāo)：wp=0.2rad, ws=0.3rad, p=1dB, s=40dB)；2. 用雙線性變換法頻率轉(zhuǎn)換公式，將上面的技術(shù)指標(biāo)轉(zhuǎn)換為巴特沃斯模擬濾波器的技術(shù)指標(biāo)；3. 調(diào)用MATLAB函數(shù)buttord和butter，設(shè)計(jì)該模擬濾波器；4. 用雙線性變換法將設(shè)計(jì)的模擬濾波器轉(zhuǎn)換為數(shù)字濾波器，求出該濾波器的系統(tǒng)函數(shù)并畫出它的直接型網(wǎng)絡(luò)結(jié)構(gòu)圖；5. 畫出該濾波器的幅度特性曲線；6. 根據(jù)1所確定的技術(shù)指標(biāo)，調(diào)用MATLAB函數(shù)buttord和butter，直接設(shè)計(jì)數(shù)字濾波器，觀察設(shè)計(jì)結(jié)果與上面用雙線性變換法的設(shè)計(jì)結(jié)果是否相同。8. 濾波仿真：調(diào)用MATLAB工具箱函數(shù)filter對(duì)下面

3、給出的帶加性干擾的信號(hào)x(n)進(jìn)行濾波，觀察濾波效果(繪制濾波前后信號(hào)的時(shí)域和頻域波形)。實(shí)驗(yàn)信號(hào)為心電圖信號(hào)采樣序列，人體心電圖信號(hào)在測(cè)量過程中往往受到工業(yè)高頻干擾，所以必須經(jīng)過低通濾波處理后，才能作為判斷心臟功能的有用信息。下面給出一實(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,

4、 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 四、實(shí)驗(yàn)用MATLAB函數(shù)簡(jiǎn)介filter功能：一維數(shù)字濾波器直接型實(shí)現(xiàn)函數(shù)。調(diào)用格式：yn=filter(B, A, xn)：按直接型實(shí)現(xiàn)結(jié)構(gòu)對(duì)輸入信號(hào)向量xn濾波，返回濾波器輸出信號(hào)向量yn，調(diào)用參數(shù)B和A分別為濾波器系統(tǒng)函數(shù)的分子和分母多項(xiàng)式系數(shù)向量。其實(shí)質(zhì)是求解差分方程：A(1)y(n)=B(1)x(n)+B(2)x(n-1)+B(M+1)x(n-M)-A(2)y(n-1)-A(N+1)y(n-N如果A(1)不

5、等于1時(shí)，則對(duì)系數(shù)關(guān)于A(1)歸一化后計(jì)算輸出信號(hào)y(n)。當(dāng)A=1時(shí)，對(duì)應(yīng)IIR-DF的直接型實(shí)現(xiàn)。五、實(shí)驗(yàn)報(bào)告要求1. 寫出所設(shè)計(jì)的數(shù)字濾波器的主要技術(shù)指標(biāo)、系統(tǒng)函數(shù)、網(wǎng)絡(luò)結(jié)構(gòu)；2. 寫出兩種IIR數(shù)字濾波器的設(shè)計(jì)原理、步驟；3. 寫出仿真濾波試驗(yàn)的結(jié)果。附錄1：用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的程序：%輸入信號(hào)及其波形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

6、,-2,-2,0,0,-2,-2,-2,-2,0;n=0:55;subplot(2,2,1);stem(n,x,'.');axis(0 60 -100 50);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('心電圖信號(hào)采樣序列x(n)');X=fft(x,256);subplot(2,2,2)k=0:255;plot(2*pi*k/256,abs(X)xlabel('k');ylabel('X(k)');t

7、itle('心電圖信號(hào)采樣序列的頻譜');%用雙線性變換法濾波器設(shè)計(jì)IIR-DFT=1;Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=40;Omegap=2*tan(Wp/2)/T; %模擬波器通帶截止頻率Omegas=2*tan(Ws/2)/T; %模擬波器通帶截止頻率Fs=1/T; %采樣頻率N,Wc=buttord(Omegap,Omegas,Rp,Rs,'s'); %計(jì)算巴特沃斯模擬濾波器階數(shù)N及%3dB截止頻率WcZ,P,K=buttap(N); %求butterworth歸一化原型H(p)，零極點(diǎn)增益形式b,a=zp2tf(Z,P,K);

8、% 轉(zhuǎn)換成歸化系統(tǒng)函數(shù)形式，H(p)為分子分母多項(xiàng)式形式bt,at=lp2lp(b,a,Wc); %去歸一化轉(zhuǎn)換為實(shí)際低通Ha(s)% b1,a1=butter(N,Wc,'s');%可以代替以上三個(gè)步驟 bb,ab=bilinear(bt,at,Fs); %雙線性變換得H(z)，F(xiàn)s為采樣頻率%求數(shù)字濾波器的幅頻特性H,w=freqz(bb,ab,100); %對(duì)應(yīng)采樣間隔為0.02*pimag=abs(H);db=20*log10(mag+eps)/max(mag);subplot(2,2,3);plot(w/pi,db);axis(0,0.6,-50,10);grid o

9、ntitle('濾波器的幅頻響應(yīng)曲線');%用設(shè)計(jì)濾波器對(duì)輸入信號(hào)進(jìn)行濾波y=filter(bb,ab,x);subplot(2,2,4);n=0:55;stem(n,y,'.');axis(0 60 -100 50);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('進(jìn)行濾波后的心電圖信號(hào)');附錄2：直接設(shè)計(jì)IIR巴特沃斯數(shù)字濾波器的程序：%直接設(shè)計(jì)巴特沃斯數(shù)字濾波器Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs

10、=40;wp=Wp/pi;ws=Ws/pi;%求數(shù)字濾波器的歸一化頻率N,wc=buttord(wp,ws,Rp,Rs); %計(jì)算巴特沃斯數(shù)字濾波器階數(shù)N及3dB截止頻率b1,a1=butter(N,wc);%設(shè)計(jì)巴特沃斯數(shù)字濾波器%求數(shù)字濾波器的幅頻特性H,w=freqz(b1,a1,100); %對(duì)應(yīng)采樣間隔為0.02*pimag=abs(H);db=20*log10(mag+eps)/max(mag);subplot(3,1,1);plot(w/pi,db);axis(0,0.6,-50,10);grid ontitle('濾波器的幅頻響應(yīng)曲線');%輸入信號(hào)及其波形x=

11、-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;n=0:55;subplot(3,1,2);stem(n,x,'.');axis(0 60 -100 50);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)'

12、);title('心電圖信號(hào)采樣序列x(n)');%用設(shè)計(jì)濾波器對(duì)輸入信號(hào)進(jìn)行濾波y=filter(b1,a1,x);subplot(3,1,3);n=0:55;stem(n,y,'.');axis(0 60 -100 50);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('進(jìn)行濾波后的心電圖信號(hào)');比較直接設(shè)計(jì)IIR數(shù)字濾波器和用雙線性變換法設(shè)計(jì)IIR濾波器的結(jié)果是否相同？脈沖響應(yīng)不變法和雙線性變換法在IIR數(shù)字濾波器

13、的設(shè)計(jì)中起到什么作用？比較他們的優(yōu)缺點(diǎn)脈沖響應(yīng)不變法的優(yōu)點(diǎn)： 1，模擬頻率到數(shù)字頻率的轉(zhuǎn)換時(shí)線性的； 2，數(shù)字濾波器單位脈沖響應(yīng)的數(shù)字表示近似原型的模擬濾波器單位脈沖響應(yīng)，因此時(shí)域特性逼近好。 缺點(diǎn)： 會(huì)產(chǎn)生頻譜混疊現(xiàn)象，只適合帶限濾波器雙線性變換法優(yōu)點(diǎn)： 克服多值映射得關(guān)系，可以消除頻率的混疊 缺點(diǎn)： 時(shí)域到頻域的變換是非線性的，在高頻處有較大的失真。實(shí)驗(yàn)要求：自己編寫用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的程序(調(diào)用impinvar，掌握方法即可)。若把上述濾波器改為高通濾波器，試分別用脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)對(duì)應(yīng)的數(shù)字濾波器，觀察設(shè)計(jì)結(jié)果，并說(shuō)明原因。答：用脈沖響應(yīng)不變法設(shè)計(jì)IIR

14、數(shù)字濾波器 附1雙線性變換法設(shè)計(jì)對(duì)應(yīng)的數(shù)字高通濾波器 附2脈沖響應(yīng)不變法不能用于設(shè)計(jì)數(shù)字高通濾波器和帶阻濾波器；附1、用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器T=2;  %設(shè)置采樣周期為2 fs=1/T;  %采樣頻率為周期倒數(shù) Wp=0.25*pi/T;  Ws=0.35*pi/T;  %設(shè)置歸一化通帶和阻帶截止頻率 Ap=20*log10(1/0.9); As=20*log10(1/0.18);  %設(shè)置通帶最大和最小衰減N,Wc=buttord(Wp,Ws,Ap,As,'s')

15、;  %調(diào)用butter函數(shù)確定巴特沃斯濾波器階數(shù) B,A=butter(N,Wc,'s');  %調(diào)用butter函數(shù)設(shè)計(jì)巴特沃斯濾波器 W=linspace(0,pi,400*pi);  %指定一段頻率值 hf=freqs(B,A,W);   %計(jì)算模擬濾波器的幅頻響應(yīng) subplot(2,1,1); plot(W/pi,abs(hf)/abs(hf(1);  %繪出巴特沃斯模擬濾波器的幅頻特性曲線grid on; title('巴特沃斯

16、模擬濾波器'); xlabel('Frequency/Hz'); ylabel('Magnitude'); D,C=impinvar(B,A,fs);  %調(diào)用脈沖響應(yīng)不變法 Hz=freqz(D,C,W);  %返回頻率響應(yīng) subplot(2,1,2); plot(W/pi,abs(Hz)/abs(Hz(1);  %繪出巴特沃斯數(shù)字低通濾波器的幅頻特性曲線grid on; title('巴特沃斯數(shù)字濾波器'); xl

17、abel('Frequency/Hz'); ylabel('Magnitude'); 附2、雙線性變換法設(shè)計(jì)對(duì)應(yīng)的數(shù)字高通濾波器Fs=10000; Ts=1/Fs; %采樣頻率和采樣周期wp=3000*2*pi;ws=2000*2*pi; %低通、高通數(shù)字濾波器截止頻率Rp=3;Rs=14; %通帶波紋和阻帶衰減 wpl=Fs*2*tan(ws*Ts/2); %按頻率轉(zhuǎn)換公式進(jìn)行轉(zhuǎn)換wsl=Fs*2*tan(wp*Ts/2); N,Wn=buttord(wpl,wsl,Rp

18、,Rs,'s'); %計(jì)算模擬濾波器的最小階數(shù)z,p,k=buttap(N); %設(shè)計(jì)模擬原型低通濾波器 b,a=zp2tf(z,p,k); %變零極點(diǎn)增益形式為系統(tǒng)傳遞函數(shù)形式H,w=freqs(b,a); %求模擬低通濾波器的頻率響應(yīng) figure; plot(w,abs(H);grid; %繪制低通頻率特性title('模擬低通原型'); xlabel('w');ylabel('幅度');e1,f1=lp2lp(b,a,2000*2*p

19、i)  %原型低通轉(zhuǎn)換為低通模擬濾波器 h1,w1=freqs(e1,f1); %求模擬低通濾波器的頻率響應(yīng)figure; plot(w1/(2*pi),abs(h1); %繪制高通頻率特性grid; title('模擬低通濾波器'); xlabel('頻率Hz');ylabel('幅度'); r1,q1=bilinear(e1,f1,10000); %模擬低通濾波器轉(zhuǎn)換為數(shù)字濾波 器 s1,w2=freqz(r1,q1); %求數(shù)字低通濾波器的頻率響應(yīng)figure; plot(w2/pi, 20*log