IIR濾波器,kalman濾波器綜合設(shè)計(jì)

1、基于MATLAB的IIR數(shù)字帶通濾波器的設(shè)計(jì)設(shè)計(jì)利用沖擊響應(yīng)不變法設(shè)計(jì)一個(gè)切比雪夫帶通濾波器，設(shè)計(jì)的指標(biāo)為：通帶邊緣頻率 w1 = 0.4pi,w2 = 0.8pi阻帶邊緣頻率 w3 = 0.3pi,w4 = 1.0pi通帶波紋二0.5dB阻帶衰減二20 dB解:1.基本原理濾波器的傳遞函數(shù)的一般形式為：MMdzk(z zr)H(z)七k專一1HZ*(z zp)k 1k 1當(dāng) M<=N, N 階 IIR 系統(tǒng)；M>=N, H(z): N 階 IIR 系統(tǒng) +(M-N)階的FIR系統(tǒng)，以上兩種表示等價(jià)，部分分式形式和零極點(diǎn)增益形式IIR系統(tǒng)的逼近,就是找到濾波器的系數(shù)ak, bk，或

2、者是系統(tǒng)的零極點(diǎn)和增益(z,p,k)。切比雪夫I型的幅度平方函數(shù)為：H(j)|2112cN()CH(j )|的特點(diǎn)如下:(1)=0時(shí),N為偶數(shù)時(shí),Ha(j0)| =1,1+ 2當(dāng)N為奇數(shù)時(shí),Ha(jO)| =0.C時(shí),Ha(j)1=1.1+ 2即此時(shí)所有的幅度函數(shù)曲線都經(jīng)過治點(diǎn)，C即為切比雪夫?yàn)V波器的通帶截止頻率.(3) 在通帶內(nèi)，|Ha(j )|在1：干1亍之間等波紋地起伏。(4) 在通帶外，隨著 的增大，|Ha(j )|迅速單調(diào)的趨近于零。該濾波器在通帶內(nèi)具有等波紋起伏特性，在阻帶內(nèi)則單調(diào)下降且具有更大的衰減.相比于巴特沃斯濾波器，階數(shù)N較小.2沖擊響應(yīng)不變法沖擊響應(yīng)不變法是使數(shù)字濾波器的

3、單位沖擊響應(yīng)序列h(n)模仿模擬濾波器的單位沖擊響應(yīng)ha(t)。將模擬濾波器的單位沖擊響應(yīng)加以等間隔抽樣，使h(n)正好等于ha(t)的抽樣值，即滿足：h(n)= h a(nT)1 2既有:H(z)|z=eST = ZHa(S-j k)I k=-I沖擊響應(yīng)不變法是將模擬濾波器的 S平面變換成數(shù)字濾波器的Z平面。沖擊響應(yīng)不變法使得數(shù)字濾波器的單位脈沖響應(yīng)完全模仿模擬 濾波器的單位沖激響應(yīng)，也就是時(shí)域逼近良好。一個(gè)線性相位的模擬濾波器通過脈沖響應(yīng)不變法得到的仍然是一個(gè)線性相位的數(shù)字濾波3 / Tj|mZS平面Z平面但是，由于任何一個(gè)實(shí)際的模擬濾波器頻率響應(yīng)都不是嚴(yán)格限帶 的，變換后就會(huì)產(chǎn)生頻率響應(yīng)

4、的混疊失真。 這時(shí)數(shù)字濾波器的頻響就 不同于原模擬濾波器的頻響，而帶有一定的失真。當(dāng)模擬濾波器的頻 率響應(yīng)在折疊頻率以上處衰減越大、越快時(shí)，變換后頻率響應(yīng)混疊失 真就越小。這時(shí)，采用脈沖響應(yīng)不變法設(shè)計(jì)的數(shù)字濾波器才能得到良 好的效果。3 .設(shè)計(jì)流程：(1) 確定數(shù)字濾波器的技術(shù)指標(biāo)。(2) 將數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換成相應(yīng)的模擬濾波器的技術(shù)指標(biāo)。(3) 按照模擬低通濾波器的技術(shù)指標(biāo)設(shè)過渡模低通濾波器。(4) 用脈沖響應(yīng)不變法，將模擬濾波轉(zhuǎn)換成數(shù)字濾波器4. matlab 實(shí)現(xiàn)主要步驟 :(1) 先由數(shù)字頻率轉(zhuǎn)化為模擬的頻率 , 由于采用脈沖響應(yīng)不變法 , 所 以:Omega1 =fs*w1;O

5、mega2 = fs*w2;Omega3 =fs*w3;Omega4 = fs*w4;采樣頻率 fs = 2000HZ由函數(shù)N,Wn=cheb1ord(Wp,Ws,Rp,Rs,'s');z,p,k = cheb1ap(N,Rp); 求出模擬低通濾波器。 再通過bt,at = lp2bp(bb,aa,Wo,Bw)函數(shù)實(shí)現(xiàn)由模擬低通濾波 器向帶通濾波器的轉(zhuǎn)換 .其中 :Wo = sqrt(Omega1*Omega2);Bw = Omega2 - Omega1; 進(jìn)行頻率變換后實(shí)現(xiàn)了帶通模擬濾波器(3) : 采用脈沖響應(yīng)不變法 : 利用函數(shù) :bz, az=impinvar(bt,a

6、t,fs)將上述帶通模擬濾波器轉(zhuǎn)換為數(shù)字帶通濾波器 . 其中 fs 為采樣頻率(4) : 驗(yàn)證：預(yù)期輸出信號(hào) :uzs=sin(2*pi*451*nn)加入噪聲后的信號(hào)為：uz = 0.5*cos(2*pi*21* nn) + uzs+ 2*si n(2*pi*1000* nn)觀察經(jīng)過該濾波器后的輸出波形。其中 Fn = 100;dn二 1/fs;n=0:F n-1; nn二 n*dn;5. matlab實(shí)現(xiàn)結(jié)果：(1)根據(jù)濾波器的設(shè)計(jì)指標(biāo)，首先實(shí)現(xiàn)一個(gè)模擬的低通濾波器，其曲線如圖一所示：圖一：模擬低通濾波器(2)將設(shè)計(jì)的模擬低通濾波器進(jìn)行頻率變換，變換成模擬的帶通濾波器,其幅頻,相頻曲線如

7、圖2所示：圖二：模擬帶通濾波器(3)沖擊響應(yīng)不變法實(shí)現(xiàn)數(shù)字帶通濾波器的帶通曲線如圖三所示0.1.0.00.0.421J86)4-2/匸1000.0.1°.0.2 °0.3 00.4 °，0.5 °%.6 00.70.0.8 0.0.91°.0其幅頻相頻曲線如圖四所示32- 4- 6 - - - kDdr eanhnaa M2esanpJseslasHr escocch01 11-11 00 0M0 0.10.80.910.20.3NNmmazelFFqqncH0.40.5(圖四：帶通濾波器的幅頻和相頻曲線(5) 為了對(duì)設(shè)計(jì)的數(shù)字帶通濾波器的性能

8、進(jìn)行驗(yàn)證，本設(shè)計(jì)加入了噪聲信號(hào)，如圖五所示即為原來理想的信號(hào)波形uzs=si n(2*pi*451* nn)(6) 如圖六所示即為加入噪聲信號(hào)后的波形曲線。噪聲信號(hào)的曲線方程為：uz=0.5*cos(2*pi*21* nn )+uzs+2*si n(2*pi*1000* nn)(7) 如圖七所示即為通過濾波器后輸出的波形1.5預(yù)輸出信號(hào)波形圖五：理想波形圖七：濾波后輸出波形Matlab程序代碼:w1 = 0.4*pi,w2 = 0.8*pi; %digital signal freqRp=0.5;Rs = 20;w3 = 0.3*pi,w4=1*pi;fs = 2000;Omega1 =fs*

9、w1;Omega2 = fs*w2;Omega3 =fs*w3;Omega4 = fs*w4;f1 =0mega1/2/pi;f2 = Omega2/2/pi;f3 =Omega3/2/pi;f4 = Omega4/2/pi;Wp = Omega1,Omega2;Ws = Omega3,Omega4;N,Wn=cheb1ord(Wp,Ws,Rp,Rs,'s');z,p,k = cheb1ap(N,Rp);bs,as = zp2tf(z,p,k);n = 0:0.01:2;hs,ws=freqs(bs,as,n);figure;plot(ws,abs(hs).八2); %low

10、pass filtergrid;bb=k*real(poly(z);aa=real(poly(p);Wo = sqrt(Omega1*Omega2);Bw = Omega2 - Omega1;bt,at = lp2bp(bb,aa,Wo,Bw);h,w = freqs(bt,at);figure;subplot(2,1,1);plot(w/2/pi,(abs(h).A2);axis(100 1500 0 1);grid; subplot(2,1,2);plot(w/2/pi,angle(h)*180/pi);grid;H = tf(bt,at);bz, az=impinvar(bt,at,fs

11、);hz,wz=freqz(bz,az,fs);figure;plot(wz/pi,abs(hz);grid;figure;freqz(bz,az);Fn = 100;dn= 1/fs;n=0:Fn-1;nn= n*dn;uzs=sin(2*pi*451*nn);uz = 0.5*cos(2*pi*21*nn) + uzs+ 2*sin(2*pi*1000*nn); figure;plot(nn,uz);grid;axis(0 0.05 0 2);ylabel('uz');Y = filter(bz,az,uz);figure;plot(nn,uzs);axis(0 0.05

12、0 1.5);xlabel(' 預(yù)輸出信號(hào)波形 '); grid;figure;plot(nn,Y);axis(0 0.05 0 1.5);xlabel(' 實(shí)際輸出波形 '); grid;設(shè)計(jì)二 :設(shè)計(jì) 一 Kalman 濾波器的設(shè)計(jì)1. 概要 :Kalman 濾波理論是 Wiener 濾波理論的發(fā)展， 具有如下特點(diǎn) (1): 數(shù) 學(xué)公式用狀態(tài)空間概念描述； (2): 它的解釋遞推計(jì)算的，與 Wiener 濾波器不同，Kalman濾波器是一種自適應(yīng)濾波器。對(duì)于一般的時(shí)變系統(tǒng) , 系統(tǒng)狀態(tài)方程和觀測(cè)方程可表示為一下形式 :(1) 狀態(tài)方程 (state equa

13、tion)x(n)= F (n,n- 1)x(n- 1)+G(n,n- 1)v1(n- 1)N ' 1狀態(tài)向量：x(n)?C狀態(tài)轉(zhuǎn)移矩陣：F (n,n- 1)? C狀態(tài)噪聲輸入矩陣：G(n,n- 1)? CN，S系統(tǒng)狀態(tài)噪聲：w(n- 1)?CS'1其中噪聲w (n)為均值為零的高斯白噪聲(2) 觀測(cè)方程(measurement equation)z(n)= C (n)x (n)+ v2 (n)觀測(cè)向量：z(n)? CM r觀測(cè)矩陣：C(n)? CM，N觀測(cè)噪聲：v2 (n)? CM 1其中噪聲v2(n為均值為零的高斯白噪聲2. 應(yīng)用 ：假設(shè)一理想質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡是一理想的正弦波

14、.如:y=5*sin (pi*t*0.5)但是當(dāng)質(zhì)點(diǎn)受到外部環(huán)境的干擾后就會(huì)產(chǎn)生一定的偏移 ，利用函 數(shù)noise二randn(1,200) 產(chǎn)生隨機(jī)的信號(hào)進(jìn)行模擬，所以實(shí)際的質(zhì)點(diǎn) 偏離中心位移的曲線為y_true= y + noise.可令質(zhì)點(diǎn)的運(yùn)動(dòng)模型為 y二sin (x),質(zhì)點(diǎn)做勻速直線運(yùn)動(dòng)。令其經(jīng)過的位移為s(n),速度為v(n),所以有：s(n)=s(n-1)+Tv(n-1) v(n)=v(n-1)對(duì)于質(zhì)點(diǎn)的干擾可以視為質(zhì)點(diǎn)產(chǎn)生的加速度，從而有：1 2s( n)=s( n-1)+Tv( n-1)+ T2a( n-1)v( n)=v( n-1)+ Ta( n-1)其中:a表示質(zhì)點(diǎn)受到的

15、加速度T表示米樣的間隔從而狀態(tài)方程可以表示為s(n)1 Ts(n-1)T2= +2a(n-1)v(n)0 1v(n-1)T令 x(n)=(s(n),v(n)T廠zv獰 G (n, n -二家± 1 (n-1)=a(n-1)狀態(tài)方程化為：x(n)= F(n,n- 1)x( n- 1)+G (n,n- 1)v1 (n- 1)其中狀態(tài)轉(zhuǎn)移矩陣F(n,n-1)=0 1觀測(cè)方程 z(n)=C(n)x(n) + V ?(n)其中 C(n)=(1 0);V2(n)二no ise3. matlab 實(shí)現(xiàn)：(1)理想質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡如圖一所示-80-4-62468101214161820時(shí)間/s圖一：理

16、想質(zhì)點(diǎn)運(yùn)動(dòng)軌跡理想質(zhì)點(diǎn)運(yùn)動(dòng)軌跡6 42 O-2m卩移位點(diǎn)質(zhì)(2)實(shí)際的質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡，在matlab中用加上均值為零的咼斯白噪聲進(jìn)行模擬，波形如圖二所示mm/移位點(diǎn)質(zhì)1i；：ii.爲(wèi)1 1 AA仍:j詁"；Jf'l1i jni| J：=l'l11i 1jIJf 1i11-i11rI1Ik11J:Pl I! I?* 1LI j/11.1 t1 *h 1 L1'Ah 1 7叫nJ1實(shí)際的質(zhì)點(diǎn)運(yùn)動(dòng)軌跡-4-6-8186 42 O-20246810121416時(shí)間/s圖二：實(shí)際的質(zhì)點(diǎn)運(yùn)動(dòng)軌跡(3) 經(jīng)過Kalman濾波器后的輸出的波形如圖三所示卡爾曼濾波后輸岀質(zhì)點(diǎn)最優(yōu)

17、化運(yùn)動(dòng)位移86420-2-4-6-818m嚴(yán)移位點(diǎn)質(zhì)02468101214162020時(shí)間/s4. matlab 程序代碼% 產(chǎn)生正弦波信號(hào) %T=0.1;t=0:T:20-0.1; y=5*sin(pi*t*0.5);figure;plot(t,y);grid;axis(0,20,-8,8);xlabel(' 時(shí)間 /s');ylabel(' 質(zhì)點(diǎn)位移 /mm'); title(' 理想質(zhì)點(diǎn)運(yùn)動(dòng)軌跡 ');% 添加噪聲 %noise=randn(1,200);noise = randn(1,200); noise=noise-mean(nois

18、e); noise=sqrt(var(noise)*noise/sqrt(var(noise);Z=y+noise; figure; plot(t,Z);grid;axis(0,20,-8,8);xlabel(' 時(shí)間 /s ');ylabel(' 質(zhì)點(diǎn)位移 /mm');title(' 實(shí)際的質(zhì)點(diǎn)運(yùn)動(dòng)軌跡 ');% 建立系統(tǒng)模型 %X=zeros(2,200);X(:,1)=0,1;F=1,T;0,1;V=1/2*(T)A2 T'%V=T T; how to get the value of matix V %V=1 1'C=1,0;P=0,0;0,1;%Q1=(0.25)A2;%Q2=(0.25)A2;Q1=cov(noise);Q2=cov(noise);% 卡爾曼算法 %