DSP實驗5 IIR濾波器的設計與信號濾波

實驗五IIR濾波器的設計與信號濾波一、實驗目的設計IIR濾波器，實現(xiàn)對存在加性干擾的時域離散信號進行濾波。二、實驗內容已知帶加性干擾的信號用x(n)表示，x(n)=xs(n)+η(n)，式中xs(n)是有用的信號，是一個0～0.2πrad的帶限信號。η(n)是一個干擾信號，其頻譜分布在0.3πrad以上。要求設計一個巴特沃斯IIR數字濾波器對信號x(n)進行濾波，將干擾η(n)濾除。要求在xs(n)所在的通帶內濾波器幅度平坦，在0.2πrad處幅度衰減不大于1dB，在噪聲所在的0.3πrad以上的頻帶內濾波器幅度衰減大于等于40dB。三、實驗步驟1.根據題目要求確定要設計的數字濾波器的技術指標(低通濾波器指標：wp=0.2πrad,ws=0.3πrad,αp=1dB,αs=40dB)；2.用雙線性變換法頻率轉換公式，將上面的技術指標轉換為巴特沃斯模擬濾波器的技術指標；3.調用MATLAB函數buttord和butter，設計該模擬濾波器；4.用雙線性變換法將設計的模擬濾波器轉換為數字濾波器，求出該濾波器的系統(tǒng)函數并畫出它的直接型網絡結構圖；5.畫出該濾波器的幅度特性曲線；6.根據1所確定的技術指標，調用MATLAB函數buttord和butter，直接設計數字濾波器，觀察設計結果與上面用雙線性變換法的設計結果是否相同。8.濾波仿真：調用MATLAB工具箱函數filter對下面給出的帶加性干擾的信號x(n)進行濾波，觀察濾波效果(繪制濾波前后信號的時域和頻域波形)。實驗信號為心電圖信號采樣序列，人體心電圖信號在測量過程中往往受到工業(yè)高頻干擾，所以必須經過低通濾波處理后，才能作為判斷心臟功能的有用信息。下面給出一實際心電圖信號采樣序列樣本x(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}四、實驗用MATLAB函數簡介filter功能：一維數字濾波器直接Ⅱ型實現(xiàn)函數。調用格式：yn=filter(B,A,xn)：按直接Ⅱ型實現(xiàn)結構對輸入信號向量xn濾波，返回濾波器輸出信號向量yn，調用參數B和A分別為濾波器系統(tǒng)函數的分子和分母多項式系數向量。其實質是求解差分方程：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時，則對系數關于A(1)歸一化后計算輸出信號y(n)。當A=1時，對應IIR-DF的直接Ⅱ型實現(xiàn)。五、實驗報告要求1.寫出所設計的數字濾波器的主要技術指標、系統(tǒng)函數、網絡結構；2.寫出兩種IIR數字濾波器的設計原理、步驟；3.寫出仿真濾波試驗的結果。附錄1：用雙線性變換法設計IIR數字濾波器的程序：%輸入信號及其波形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('心電圖信號采樣序列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('心電圖信號采樣序列的頻譜');%用雙線性變換法濾波器設計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');%計算巴特沃斯模擬濾波器階數N及%3dB截止頻率Wc[Z,P,K]=buttap(N);%求butterworth歸一化原型H(p)，零極點增益形式[b,a]=zp2tf(Z,P,K);%轉換成歸化系統(tǒng)函數形式，H(p)為分子分母多項式形式[bt,at]=lp2lp(b,a,Wc);%去歸一化轉換為實際低通Ha(s)%[b1,a1]=butter(N,Wc,'s');%可以代替以上三個步驟[bb,ab]=bilinear(bt,at,Fs);%雙線性變換得H(z)，F(xiàn)s為采樣頻率%求數字濾波器的幅頻特性[H,w]=freqz(bb,ab,100);%對應采樣間隔為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=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('進行濾波后的心電圖信號');附錄2：直接設計IIR巴特沃斯數字濾波器的程序：%直接設計巴特沃斯數字濾波器Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=40;wp=Wp/pi;ws=Ws/pi;%求數字濾波器的歸一化頻率[N,wc]=buttord(wp,ws,Rp,Rs);%計算巴特沃斯數字濾波器階數N及3dB截止頻率[b1,a1]=butter(N,wc);%設計巴特沃斯數字濾波器%求數字濾波器的幅頻特性[H,w]=freqz(b1,a1,100);%對應采樣間隔為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('濾波器的幅頻響應曲線');%輸入信號及其波形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('心電圖信號采樣序列x(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('進行濾波后的心電圖信號');比較直接設計IIR數字濾波器和用雙線性變換法設計IIR濾波器的結果是否相同？脈沖響應不變法和雙線性變換法在IIR數字濾波器的設計中起到什么作用？比較他們的優(yōu)缺點脈沖響應不變法的優(yōu)點：

1，模擬頻率到數字頻率的轉換時線性的；

2，數字濾波器單位脈沖響應的數字表示近似原型的模擬濾波器單位脈沖響應，因此時域特性逼近好。缺點：

會產生頻譜混疊現(xiàn)象，只適合帶限濾波器雙線性變換法優(yōu)點：

克服多值映射得關系，可以消除頻率的混疊缺點：

時域到頻域的變換是非線性的，在高頻處有較大的失真。實驗要求：自己編寫用脈沖響應不變法設計IIR數字濾波器的程序(調用impinvar，掌握方法即可)。若把上述濾波器改為高通濾波器，試分別用脈沖響應不變法和雙線性變換法設計對應的數字濾波器，觀察設計結果，并說明原因。答：用脈沖響應不變法設計IIR數字濾波器附1雙線性變換法設計對應的數字高通濾波器附2脈沖響應不變法不能用于設計數字高通濾波器和帶阻濾波器；附1、用脈沖響應不變法設計IIR數字濾波器T=2;

%設置采樣周期為2

fs=1/T;

%采樣頻率為周期倒數

Wp=0.25*pi/T;

Ws=0.35*pi/T;

%設置歸一化通帶和阻帶截止頻率

Ap=20*log10(1/0.9);

As=20*log10(1/0.18);

%設置通帶最大和最小衰減[N,Wc]=buttord(Wp,Ws,Ap,As,'s');

%調用butter函數確定巴特沃斯濾波器階數

[B,A]=butter(N,Wc,'s');

%調用butter函數設計巴特沃斯濾波器

W=linspace(0,pi,400*pi);

%指定一段頻率值

hf=freqs(B,A,W);

%計算模擬濾波器的幅頻響應

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);

%調用脈沖響應不變法

Hz=freqz(D,C,W);

%返回頻率響應

subplot(2,1,2);

plot(W/pi,abs(Hz)/abs(Hz(1)));

%繪出巴特沃斯數字低通濾波器的幅頻特性曲線grid

on;

title('巴特沃斯數字濾波器');

xlabel('Frequency/Hz');

ylabel('Magnitude');

附2、雙線性變換法設計對應的數字高通濾波器Fs=10000;

Ts=1/Fs;

%采樣頻率和采樣周期wp=3000*2*pi;ws=2000*2*pi;

%低通、高通數字濾波器截止頻率Rp=3;Rs=14;

%通帶波紋和阻帶衰減

wpl=Fs*2*tan(ws*Ts/2);

%按頻率轉換公式進行轉換wsl=Fs*2*tan(wp*Ts/2);

[N,Wn]=buttord(wpl,wsl,Rp,Rs,'s');

%計算模擬濾波器的最小階數[z,p,k]=buttap(N);

%設計模擬原型低通濾波器

[b,a]=zp2tf(z,p,k);

%變零極點增益形式為系統(tǒng)傳遞函數形式[H,w]=freqs(b,a);

%求模擬低通濾波器的頻率響應

figure;

plot(w,abs(H));grid;

%繪制低通頻率特性title('模擬低通原型');

xlabel('w');ylabel('幅度');[e1,f1]=lp2lp(b,a,2000*2*pi)

%原型低通轉換為低通模擬濾波器

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

%求模擬低通濾波器的頻率響應figu