IIR濾波器,kalman濾波器綜合設(shè)計(jì)知識(shí)分享

1、Good is good, but better carries it.精益求精，善益求善。IIR濾波器,kalman濾波器綜合設(shè)計(jì)-基于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阻帶衰減=20dB解:1.基本原理濾波器的傳遞函數(shù)的一般形式為:當(dāng)M=N,H(z):N階IIR系統(tǒng)+(M-N)階的FIR系統(tǒng),以上兩種表示等價(jià)，部分分式形式和零極點(diǎn)增益形式IIR系統(tǒng)的逼近,就是找到濾波器的系數(shù)ak,bk，或者是系統(tǒng)的零極點(diǎn)和

2、增益(z,p,k)。切比雪夫I型的幅度平方函數(shù)為:SKIPIF10的特點(diǎn)如下:當(dāng)SKIPIF10=0時(shí),N為偶數(shù)時(shí),SKIPIF10=SKIPIF10,當(dāng)N為奇數(shù)時(shí)，SKIPIF10=0.當(dāng)SKIPIF10=SKIPIF10c時(shí)，SKIPIF10即此時(shí)所有的幅度函數(shù)曲線都經(jīng)過(guò)SKIPIF10點(diǎn),SKIPIF10c即為切比雪夫?yàn)V波器的通帶截止頻率.在通帶內(nèi),SKIPIF10在1SKIPIF10SKIPIF10之間等波紋地起伏。在通帶外,隨著SKIPIF10的增大,SKIPIF10迅速單調(diào)的趨近于零。該濾波器在通帶內(nèi)具有等波紋起伏特性,在阻帶內(nèi)則單調(diào)下降且具有更大的衰減.相比于巴特沃斯濾波器,階數(shù)

3、N較小.2沖擊響應(yīng)不變法沖擊響應(yīng)不變法是使數(shù)字濾波器的單位沖擊響應(yīng)序列h(n)模仿模擬濾波器的單位沖擊響應(yīng)ha(t)。將模擬濾波器的單位沖擊響應(yīng)加以等間隔抽樣,使h(n)正好等于ha(t)的抽樣值，即滿足:h(n)=ha(nT)既有:SKIPIF10沖擊響應(yīng)不變法是將模擬濾波器的s平面變換成數(shù)字濾波器的z平面。沖擊響應(yīng)不變法使得數(shù)字濾波器的單位脈沖響應(yīng)完全模仿模擬濾波器的單位沖激響應(yīng)，也就是時(shí)域逼近良好。一個(gè)線性相位的模擬濾波器通過(guò)脈沖響應(yīng)不變法得到的仍然是一個(gè)線性相位的數(shù)字濾波器。但是，由于任何一個(gè)實(shí)際的模擬濾波器頻率響應(yīng)都不是嚴(yán)格限帶的，變換后就會(huì)產(chǎn)生頻率響應(yīng)的混疊失真。這時(shí)數(shù)字濾波器的頻

4、響就不同于原模擬濾波器的頻響，而帶有一定的失真。當(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è)過(guò)渡模低通濾波器。(4)用脈沖響應(yīng)不變法，將模擬濾波轉(zhuǎn)換成數(shù)字濾波器4.matlab實(shí)現(xiàn)主要步驟:(1)先由數(shù)字頻率轉(zhuǎn)化為模擬的頻率,由于采用脈沖響應(yīng)不變法,所以:Omega1=fs*w1;Omega2=fs*w2;Omega3=fs*w3;Omega4=fs*w

5、4;采樣頻率fs=2000HZ由函數(shù)N,Wn=cheb1ord(Wp,Ws,Rp,Rs,s);z,p,k=cheb1ap(N,Rp);求出模擬低通濾波器。(2)再通過(guò)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,at,fs)將上述帶通模擬濾波器轉(zhuǎn)換為數(shù)字帶通濾波器.其中fs為采樣頻率(4):驗(yàn)證：預(yù)期輸出信號(hào):uzs=sin(2*pi*451*nn)加入噪聲

6、后的信號(hào)為:uz=0.5*cos(2*pi*21*nn)+uzs+2*sin(2*pi*1000*nn)觀察經(jīng)過(guò)該濾波器后的輸出波形。其中Fn=100;dn=1/fs;n=0:Fn-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)行頻率變換,變換成模擬的帶通濾波器,其幅頻,相頻曲線如圖2所示:圖二:模擬帶通濾波器(3)沖擊響應(yīng)不變法實(shí)現(xiàn)數(shù)字帶通濾波器的帶通曲線如圖三所示:其幅頻相頻曲線如圖四所示:圖四:帶通濾波器的幅頻和相頻曲線為了對(duì)設(shè)計(jì)的數(shù)字帶通濾波器的性能進(jìn)行驗(yàn)證,本

7、設(shè)計(jì)加入了噪聲信號(hào),如圖五所示即為原來(lái)理想的信號(hào)波形uzs=sin(2*pi*451*nn)如圖六所示即為加入噪聲信號(hào)后的波形曲線。噪聲信號(hào)的曲線方程為:uz=0.5*cos(2*pi*21*nn)+uzs+2*sin(2*pi*1000*nn)(7)如圖七所示即為通過(guò)濾波器后輸出的波形圖五:理想波形圖六:噪聲信號(hào)圖七:濾波后輸出波形Matlab程序代碼:w1=0.4*pi,w2=0.8*pi;%digitalsignalfreqRp=0.5;Rs=20;w3=0.3*pi,w4=1*pi;fs=2000;Omega1=fs*w1;Omega2=fs*w2;Omega3=fs*w3;Omega

8、4=fs*w4;f1=Omega1/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);%lowpassfiltergrid;bb=k*real(poly(z);aa=real(poly(p);Wo=sqrt(Omega1*Omeg

9、a2);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).2);axis(100150001);grid;subplot(2,1,2);plot(w/2/pi,angle(h)*180/pi);grid;H=tf(bt,at);bz,az=impinvar(bt,at,fs);hz,wz=freqz(bz,az,fs);figure;plot(wz/pi,abs(hz);grid;figure;freqz(bz,az);Fn=100;dn=1/

10、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(00.0502);ylabel(uz);Y=filter(bz,az,uz);figure;plot(nn,uzs);axis(00.0501.5);xlabel(預(yù)輸出信號(hào)波形);grid;figure;plot(nn,Y);axis(00.0501.5);xlabel(實(shí)際輸出波形);grid;設(shè)計(jì)二:設(shè)計(jì)一Kalman濾波器的設(shè)計(jì)1.概要:Kalman濾波理

11、論是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)方程(stateequation)狀態(tài)向量：狀態(tài)轉(zhuǎn)移矩陣：狀態(tài)噪聲輸入矩陣：系統(tǒng)狀態(tài)噪聲：其中噪聲為均值為零的高斯白噪聲(2)觀測(cè)方程(measurementequation)觀測(cè)向量:觀測(cè)矩陣:觀測(cè)噪聲:其中噪聲為均值為零的高斯白噪聲2.應(yīng)用:假設(shè)一理想質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡是一理想的正弦波.如:y=5*sin(pi*t*0.5)但是當(dāng)質(zhì)點(diǎn)受到外部環(huán)境的干擾后就

12、會(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)過(guò)的位移為s(n),速度為v(n),所以有:對(duì)于質(zhì)點(diǎn)的干擾可以視為質(zhì)點(diǎn)產(chǎn)生的加速度,從而有:其中:a表示質(zhì)點(diǎn)受到的加速度T表示采樣的間隔從而狀態(tài)方程可以表示為:令x(n)=(s(n),v(n)T,v1(n-1)=a(n-1)狀態(tài)方程化為:其中狀態(tài)轉(zhuǎn)移矩陣F(n,n-1)=觀測(cè)方程z(n)=C(n)x(n)+V2(n)其中C(n)=(10);V2(n)=noise3.matlab實(shí)現(xiàn):

13、(1)理想質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡如圖一所示:圖一:理想質(zhì)點(diǎn)運(yùn)動(dòng)軌跡(2)實(shí)際的質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡,在matlab中用加上均值為零的高斯白噪聲進(jìn)行模擬，波形如圖二所示:圖二：實(shí)際的質(zhì)點(diǎn)運(yùn)動(dòng)軌跡(3)經(jīng)過(guò)Kalman濾波器后的輸出的波形如圖三所示:圖三:輸出質(zhì)點(diǎn)的位移4.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);nois

14、e=randn(1,200);noise=noise-mean(noise);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)2T;%V=TT;howtogetthevalueofmatixV%V=11;C=1,0;P=0,0;0,1;%Q1=(0.25)2;%Q2=(0.25)2;Q1=cov(noise);Q2=cov(noise);%-卡爾曼算法-%forn=1:2