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

實(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è)0～0.2πrad的帶限信號(hào)。η(n)是一個(gè)干擾信號(hào)，其頻譜分布在0.3πrad以上。要求設(shè)計(jì)一個(gè)巴特沃斯IIR數(shù)字濾波器對(duì)信號(hào)x(n)進(jìn)行濾波，將干擾η(n)濾除。要求在xs(n)所在的通帶內(nèi)濾波器幅度平坦，在0.2πrad處幅度衰減不大于1dB，在噪聲所在的0.3πrad以上的頻帶內(nèi)濾波器幅度衰減大于等于40dB。三、實(shí)驗(yàn)步驟1.根據(jù)題目要求確定要設(shè)計(jì)的數(shù)字濾波器的技術(shù)指標(biāo)(低通濾波器指標(biāo)：wp=0.2πrad,ws=0.3πrad,α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ì)下面給出的帶加性干擾的信號(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,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ù)簡介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)不等于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,-2,-2,0,0,-2,-2,-2,-2,0];n=0:55;subplot(2,2,1);stem(n,x,'.');axis([060-10050]);holdon;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)');title('心電圖信號(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截止頻率Wc[Z,P,K]=buttap(N);%求butterworth歸一化原型H(p)，零極點(diǎn)增益形式[b,a]=zp2tf(Z,P,K);%轉(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]);gridontitle('濾波器的幅頻響應(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([060-10050]);holdon;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=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]);gridontitle('濾波器的幅頻響應(yīng)曲線');%輸入信號(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,-2,-2,0,0,-2,-2,-2,-2,0];n=0:55;subplot(3,1,2);stem(n,x,'.');axis([060-10050]);holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');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([060-10050]);holdon;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ù)字濾波器的設(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é)果，并說明原因。答：用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(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');

%調(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('巴特沃斯模擬濾波器');

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ù)字濾波器');

xlabel('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,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*pi)

%原型低通轉(zhuǎn)換為低通模擬濾波器

[h1,w1]=freqs(e1,f1);

%求模擬低通濾波器的頻率響應(yīng)figu