實驗報告-卡爾曼濾波

-了-WQRD整式-：可編輯--專業(yè)資料-----__數(shù)字信號處理實驗報告姓名：專業(yè)：通信與信息系統(tǒng)學(xué)號：日期：2015.11實驗內(nèi)容任務(wù)一：一連續(xù)平穩(wěn)的隨機信號/），自相關(guān)函數(shù)rO=e-』，信號x。為加性噪聲所干擾，噪聲是白x噪聲，測量值的離散值旅）為已知，T=0.02s，-3.2，-0.8，-14，-16，-17，-18，-3.3，-2.4,-18，-0.3，-0.4，-0.8，-19，-2.0，-1.2，-11，-14，-0.9，-0.8，10，0.2，0.5，-0.5，2.4，-0.5，0.5，-13，0.5，10，-12，0.5，-0.6，-15，-0.7，15，0.5，-0.7，-2.0，-19，-17，-11，-14，自編卡爾曼濾波遞推程序，估計信號x（t）的波形。任務(wù)二：設(shè)計一維納濾波器。（1）產(chǎn)生三組觀測數(shù)據(jù)：首先根據(jù)s（n）=as（n-1）+w（n）產(chǎn)生信號s（n），將其加噪（信噪比分別為20dB，10dB，6dB），得到觀測數(shù)據(jù)氣（n），x2（n），x3（n）。（2）估計X,（n），i=1，2，3的AR模型參數(shù)。假設(shè)信號長度為L，AR模型階數(shù)為N，分析實驗結(jié)果，并討論改變L，N對實驗結(jié)果的影響。實驗任務(wù)一卡爾曼濾波原理1.1卡爾曼濾波簡介早在20世紀(jì)40年代，開始有人用狀態(tài)變量模型來研究隨機過程，到60年代初，由于空間技術(shù)的發(fā)展，為了解決對非平穩(wěn)、多輸入輸出隨機序列的估計問題，卡爾曼提出了遞推最優(yōu)估計理論。它用狀態(tài)空間法描述系統(tǒng)，由狀態(tài)方程和量測方程所組成，即知道前一個狀態(tài)的估計值和最近一個觀測數(shù)據(jù)，采用遞推的算法估計當(dāng)前的狀態(tài)值。由于卡爾曼濾波采用遞推法，適合于計算機處理，并且可以用來處理多維和非平穩(wěn)隨機信號，現(xiàn)已廣泛應(yīng)用于很多領(lǐng)域，并取得了很好的結(jié)果。卡爾曼濾波一經(jīng)出現(xiàn)，就受到人們的很大重視，并在實踐中不斷豐富和完善，其中一個成功的應(yīng)用是設(shè)計運載體的高精度組合導(dǎo)航系統(tǒng)。卡爾曼濾波具有以下的特點（1）算法是遞推的，且狀態(tài)空間法采用在時域內(nèi)設(shè)計濾波器的方法，因而適用于多維隨機過程的估計；離散型卡爾曼算法適用于計算機處理。（2）用遞推法計算，不需要知道全部過去的值，用狀態(tài)方程描述狀態(tài)變量的動態(tài)變化規(guī)律，因此信號可以是平穩(wěn)的，也可以是非平穩(wěn)的，即卡爾曼濾波適用于非平穩(wěn)過程。--完整版學(xué)習(xí)資料分享----（3）卡爾曼濾波采取的誤差準(zhǔn)則仍為估計誤差的均方值最小。1.2卡爾曼濾波的狀態(tài)方程和測量方程假設(shè)某系統(tǒng)頃刻的狀態(tài)變量為x，狀態(tài)方程和量測方程（輸出方程）表示為Lk-1X-1+w-1yCkX+Vk七是觀測數(shù)據(jù)。匕都是均值為零的正態(tài)白其中，x是狀態(tài)變量；wk-1表示輸入信號是白噪聲；Vk是觀測噪聲;為了推導(dǎo)簡單，假設(shè)狀態(tài)變量的增益矩陣A不隨時間發(fā)生變化，氣噪聲，方差分別是Qk和氣，并且初始狀態(tài)與七是觀測數(shù)據(jù)。匕都是均值為零的正態(tài)白噪聲，方差分別是Qk和氣，并且初始狀態(tài)與wk,

w:^Cw]=0,b2=Q

kkVk:E[pk]=0q2=R,7Vk都不相關(guān)，7表示相關(guān)系數(shù)。即:wk'7w,w=Qk5kjk，j'=R5kkj其中kj1.3卡爾曼濾波的遞推算法卡爾曼濾波采用遞推算法來實現(xiàn)，其基本思想是先不考慮輸入信號wk和觀測噪聲Vk的影響，得到狀態(tài)變量和輸出信號（即觀測數(shù)據(jù)）的估計值，再用輸出信號的估計誤差加權(quán)后校正狀態(tài)變量的估計值，使?fàn)顟B(tài)變量估計誤差的均方值最小。因此，卡爾曼濾波器的關(guān)鍵是計算出加權(quán)矩陣的最佳值。當(dāng)不考慮觀測噪聲和輸入信號時，狀態(tài)方程和量測方程為AAxk=Akx-1y=Cx=CAx

kkkkkk-1.....................~.顯然，由于不考慮觀測噪聲的影響，輸出信號的估計值與實際值是有誤差的，用y表示k..-...........~....為了提高狀態(tài)估計的質(zhì)量，用輸出信號的估計誤差yk來校正狀態(tài)變量x=x=Ax,+Hkkk-1k=5-1+H..——.....-.-........■-....?其中,Hk為增益矩陣，即加權(quán)矩陣。經(jīng)過校正后的狀態(tài)變量的估計誤差及其均方值分別用xk和1表示’把未經(jīng)校正的狀態(tài)變量的估計誤差的均方值用P:表示?Axk=xk-xk八、Tx-xPk=E「raw八、Tx-xPk=EXk-xPk=E卡爾曼濾波要求狀態(tài)變量的估計誤差的均方值Pk為最小’因此卡爾曼濾波的關(guān)鍵即為通過選擇..一A一^?合適的H，使得P取得最小值。首先推導(dǎo)狀態(tài)變量的估計值xk和狀態(tài)變量的的估計誤差xk，然后計算X的均方值Pk通過化簡Pk，得到一組卡爾曼濾波的遞推公式:(y-CAX1然后計算X的均方值Pk通過化簡Pk，得到一組卡爾曼濾波的遞推公式:(y-CAX1kkk\H=P'CtCP'Ct+R)】kkkkkkkP'=APAt+Qk(kk-1k、k-1P=(I-HC)PkX=AX+Hkkk-1kkkk-1假設(shè)初始條件Ak么遞推流程如下:QkAxk-1Pk-1已知,其中X,=EL]，P=varix]o，那k-1Pk-1■=>P^=：>H^>Xk卡爾曼濾波遞推程序編程思想題目分析（1）由于信號x\）為加性噪聲所干擾，可知yk=xk+vk,（1）由于信號x\）為加性噪聲所干擾，可知yk=xk+vk,所以Ck=1（2）（3）因為（2）（3）因為rO=e-H，所以X/\m=^r)m=^im=om=^\z)=Jrm女-m=Je-mz-m=Jemz-m+e-mz-mSXXXXm=-sm=-sm=-s1+5=__L):2__1-e-1z1-e-1z-1i一e-11-e-1z-1(ezXTm=0又因為噪聲為白噪聲，所以Rk=。2=Sv（0）=1由此可知，B=1：1z］，即x（n）-e-1x（n-1）=w（n—1），可得到：A^=e-1，因為抽樣間隔T=0.02s，所以：A=e-ts=e-0.02。（4）因此w（n）=x（n+1）一e-Tsx（n），所以—rt)]=E\w(nk(n)]=1-e-2t因此Q=1-e-0.04k編程分析TOC\o"1-5"\h\z由上面的分析可知初始條件A，C，Q，R，y已知，在仿真中假設(shè)x—0，則x—0，P—1，kkkkk000由以上參數(shù)可得卡爾曼實際遞推公式fAA了了八x—e-0.02x+Hy一e-0.02xkk-1k】kk-1kk/L—p房）kkkP'—e-0.04P+1一e-0.04、P=G-H^p將得到的公式代入前面分析的遞推公式，即可進(jìn)行迭代得到結(jié)果xk。3.MATLAB源代碼根據(jù)以上分析，編寫matlab程序如下：%%%卡爾曼濾波%說明%X(k+1)=Ak*X(k)+W(k);%Y(k)=Ck*X(k)+V(k)%%clear;clc;%基本參數(shù)值A(chǔ)k=exp(-0.02);Ck=1;Qk=1-exp(-0.04);Rk=1;%初始值設(shè)置X0=0;P0=1;%觀測值y(k)Y=[-3.2-0.8-14-16-17-18-3.3-2.4-18-0.3-0.4-0.8-19-2.0-1.2...-11-14-0.90.8100.20.52.4-0.50.5-130.510-120.5-0.6-15-0.715...0.5-0.7-2.0-19-17-11-14];%數(shù)據(jù)長度N=length(Y);fork=1:Nifk==1%k=1時由初值開始計算--完整版學(xué)習(xí)資料分享----P_(k)=Ak*P0火Ak'+Qk;H(k)=P_(k)*Ck'*inv(Ck*P_(k)*Ck'+Rk);X(k)=Ak*X0+H(k)*(Y(k)-Ck*Ak*X0);I=eye(size(H(k)));P(k)=(I-H(k)*Ck)*P_(k);else%k>1時，開始遞推%遞推公式P_(k)=Ak*P(k-1)*Ak'+Qk;H(k)=P_(k)*Ck'*inv(Ck*P_(k)*Ck'+Rk);X(k)=Ak*X(k-1)+H(k)*(Y(k)-Ck*Ak*X(k-1));I=eye(size(H(k)));P(k)=(I-H(k)*Ck)*P_(k);endendM=1:N;T=0.02火M;%作圖，畫出x(t)的波形figure(1)plot(T,Y,'r','LineWidth',1);xlabel('t');ylabel('y(t)');title('卡爾曼濾波-測量信號y(t)波形')；grid;figure(2)plot(T,X,'b','LineWidth',1);xlabel('t');ylabel('x(t)');title('卡爾曼濾波-估計信號x(t)波形')；grid;4.實驗結(jié)果實驗任務(wù)二1.維納濾波器原理維納-霍夫方程rd?）=無h（m）r（k-m）=h（k）*r（k）m=0當(dāng)h（n）是一個長度為M的因果序列（即一個長度為M的FIR濾波器）時，維納-霍夫方程表述r（k）=￡h（m）r（k一m）=h（k）*r（k）k=0,1，2,…定義xdR=XXr0XXZ\r\定義R=XXr0XXZ\r\1）r1油）r（M-1）rXX（M-2）…r（0）~h一[r（0）]h2R=（1）xd,h1-MJxdr（M-1）xd則可寫成矩陣的形式，即Rxd=Rxh對上式求逆，得到h=R-Rd由以上式子可知：若已知期望信號與觀測數(shù)據(jù)的互相關(guān)函數(shù)及觀測數(shù)據(jù)的自相關(guān)函數(shù)，則可以通過矩陣求逆運算，得到維納濾波器的最佳解。同時可以看到，直接從時域求解因果的維納濾波器，當(dāng)選擇的濾波器的長度較大時，計算工作量很大，并且需要計算R的逆矩陣，從而要求的存儲量也很大預(yù)測是根據(jù)觀測到對的過去數(shù)據(jù)來估計當(dāng)前或?qū)淼男盘栔?。維納預(yù)測是已知以前時刻的p個數(shù)據(jù)x（n-1）x（n-2）…，x（n-p），估計當(dāng)前時刻n，或者未來n+N時刻的信號值，即估計1（n+N），N＞。，估計得到的結(jié)果仍然要求滿足均方誤差最小的準(zhǔn)則。信號可以預(yù)測是由于信號內(nèi)部存在著關(guān)聯(lián)性。預(yù)測是利用數(shù)據(jù)前后的關(guān)聯(lián)性，根據(jù)其中一部分推知其余部分。一步線性預(yù)測的時域解已知x（n-1），x（n-2），…，x（n-p），預(yù)測x（n），假設(shè)噪聲v（n）=0，這樣的預(yù)測成為一步線性預(yù)測。設(shè)定系統(tǒng)的單位脈沖響應(yīng)為h（n）。根據(jù)現(xiàn)行系統(tǒng)的基本理論，輸出信號y（n）=x（n）=￡ph（k）x（n-k）k=1令a=-h（k），則x(n)=-習(xí)ax(n-k)pkk=1預(yù)測誤差pkk=0e(n)=x(n)-x(n)=x(n)+習(xí)ax(n-k)=iax(npkk=0k=1其中a=1EE\x(n)+ipapkk=1pl要使均方誤差為最小值l=1，2，...，p又因為，我們可以得到eL*（h'）x（n-1止01=1，2，...，p所以rG)+iar(k-1)=0l=1，2，…，p(1)k=1由于預(yù)測器的輸出氐）是輸入信號的線性組合，所以可得:e*(n')x(n)以上說明誤差信號與輸入信號滿足正交性原理，預(yù)測誤差與預(yù)測信號值同樣滿足正交性原理。預(yù)測誤差的最小均方值Ek[Jmin=Ee*(nIx(n)-x(n)=E=E^e*(n')x(nx*(n)+明ax*(n-k)x(n)pkk=1（2）=r^(0)+ia人(k)k=1由（1）（2）聯(lián)立方程組，寫成矩陣形式可得「,（0）rG)...r(p)「1]e\e(n)2」「r(0)...rip-1)a10minxx.........r(0)P1...=一?)r(p-1)...rPa1-pp」]...0」這就是有名的Yule-Walker(維納-霍夫)方程。實驗編程思想在本實驗中，首先根據(jù)要求產(chǎn)生加噪不同的觀測數(shù)據(jù)氣。)，x2。)，與。)，然后可利用已知條件代入Yule-Walker方程，即可求解AR模型參數(shù)。在本實驗中，假設(shè)a=0.2，信號s(n)的初值0=0。MATLAB代碼functionWiener_predict(L,N)%clc;clear;%信噪比SN1=6;SN2=10;SN3=20;%產(chǎn)生信號s(n)a=0.2;W=random('norm',0,1,L,1);S(1)=0;forn=2:LS(n)=a火S(n-1)+W(n);end%產(chǎn)生觀測信號Am=sum(abs(S).八2)/L;P1=Am/(10八(SN1/20));P2=Am/(10八(SN2/20));P3=Am/(10八(SN3/20));V1=random('norm',0,P1,L,1);V2=random('norm',0,P2,L,1);V3=random('normrandom('normforn=1:LX1(n)X2(n)X3(n)0,P3,L,1);S(n)+V1(n);S(n)+V2(n);S(n)+V3(n);endsubplot(2,2,1);plot(S,'b');title('信號S(n)');ylabel('幅度');gridon;subplot(2,2,2);plot(X1,'b');title('觀測信號X1(n)');ylabel('幅度');gridon;subplot(2,2,3);plot(X2,'b');title('觀測信號X2(n)');ylabel('幅度');gridon;subplot(2,2,4);plot(X3,'b');title('觀測信號X3(n)');ylabel('幅度');gridon;fprintf('\n對X1信號來說N階模型參數(shù)和誤差分別為：\n')AR(X1,N);fprintf('\n對X2信號來說N階模型參數(shù)和誤差分別為：\n')AR(X2,N);fprintf('\n對X3信號來說N階模型參數(shù)和誤差分別為：\n')AR(X3,N);functionAR(X,N)L=length(X);rx=zeros(1,N+1);R=zeros(N+1,N+1);fori=1:(N+1)sum=0;forj=1:(L-i+1);sum=sum+X(j)*X(j+i-1);endrx(i)=sum/(L-i+1);endfori=1:N+1R(i,1:(i-1))=rx((i-1):-1:1);R(i,i:(N+1))=rx(1:(N-i+2));endzx=rx(2:(N+1));ap=inv(R(1:N,1:N))火(-zx)';a=[1,ap'];e=rx(1)+zx*ap;disp(['AR系數(shù)：'，num2str(a)]);disp(['均方誤差：',num2str(e)]);functionWiener_new1(L,N)%%產(chǎn)生三組觀測數(shù)據(jù)%信噪比(dB)SNR1=20;SNR2=10;SNR3=6;%產(chǎn)生信號s(n)a=0.2;W=random('norm',0,1,L,1);S(1)=0;forn=2:LS(n)=a火S(n-1)+W(n);end%加噪聲產(chǎn)生觀測限號X1=awgn(S,SNR1,'measured','linear');X2=awgn(S,SNR2,'measured','linear');X3=awgn(S,SNR3,'measured','linear');%畫出信號圖像on;on;on;subplot(2,2,1);plot(S,'b');title('信號S(n)');ylabel('幅度');gridon;subplot(2,2,2);plot(X1,'b');title('觀測信號X1(n)');ylabel('幅度');gridsubplot(2,2,3);plot(X2,'b');title('觀測信號X2(n)');ylabel('幅度');gridsubplot(2,2,4);plot(X3,'b');title('觀測信號X3(n)');ylabel('幅度'on;on;on;fprintf('\n對X1信號來說N階模型參數(shù)和誤差分別為：\n');AR(X1,N);fprintf('\n對X2信號來說N階模型參數(shù)和誤差分別為：\n');AR(X2,N);fprintf('\n對X3信號來說N階模型參數(shù)和誤差分別為：\n');AR(X3,N);functionAR(X,N)L=length(X);rx=zeros(1,N+1);R=zeros(N+1,N+1);fori=1:(N+1)sum=0;forj=1:(L-i+1);endrx(i)sum=sum+X(j)*X(j+i-endrx(i)sum/(L-i+1);--完整版學(xué)習(xí)資料分享---------WORD格式--可編輯--專業(yè)資料-----endfori=1:N+1R(i,1:(i-1))=rx((i-1)：-1:1);R(i,i:(N+1))=rx(1:(N-i+2));endzx=rx(2:(N+1));ap=inv(R(1:N,1:N))火(-zx)';a=[1,ap'];e=rx(1)+zx*ap;disp(['AR系數(shù)：'，num2str(a)]);disp(['均方誤差：',num2str(e)]);4.實驗結(jié)果與分析原始信號S（n）度幅觀測數(shù)據(jù)產(chǎn)生20-2-4度幅020406020-20204060-4圖1.原始信號與觀測信號（L=50）原始信號S（n）度幅度幅0204060（2）模型階數(shù)N對實驗結(jié)果的影響N=1對X1信號來說N階模型參數(shù)和誤差分別為：對X1信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.27766均方誤差：1.1289對X2信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.29326均方誤差：0.97283對X3信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.26441均方誤差：1.0531N=2對X1信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：1-0.344940.2854均方誤差：1.0958對X2信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：1-0.16960.10742均方誤差：1.1639對X3信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：1-0.195320.17033均方誤差：0.92331N=3對X1信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：1-0.106730.050928-0.19364均方誤差：1.4197對X2信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.354510.62013-0.75585均方誤差：0.95739對X3信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.122210.14428-0.34185均方誤差：0.99317N=5對X1信號來說N對X1信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.355150.56619-0.54005均方誤差：1.2405對X2信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.273430.102270.028944均方誤差：1.3557對X3信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.365940.41414-0.416650.65254-0.513270.21289-0.25080.66894-0.60712均方誤差：1.1025分析：由以上實驗結(jié)果可知：在數(shù)據(jù)的長度一定的條件下，改變由以上實驗結(jié)果可知：在數(shù)據(jù)的長度一定的條件下，改變AR模型的階數(shù)，均方誤差會改變，當(dāng)階數(shù)在某個值時，均方誤差的值最小，因此濾波器的階數(shù)對實驗結(jié)果有很大影響。在本次實驗中，仿真情況有限，在以上仿真中我們可以看到當(dāng)模型階數(shù)N為某一固定值時，均方誤差明顯較小。（3）信號長度L對實驗結(jié)果的影響L=100對X1信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：10.0229140.0078698均方誤差：1.2033對X2信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：10.0174140.19629均方誤差：1.1607對X3信號來說N階模型參數(shù)和誤差分別為:AR系數(shù)：1-0.0127890.13086

均方誤差：1.1483L=200對X1信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.176790.073726均方誤差：1.3371對X2信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.26490.16751均方誤差：0.99844對X3信號來說N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.271450.17666均方誤差：0.99289分析：由以上仿真結(jié)果可知，實驗中存在誤差，但仍然可以看出，隨著信號長度的增加，均方誤差減小，預(yù)測更準(zhǔn)確。L=100，N=1度幅觀測信號X1(n),SNR=20dB度幅觀測信號X2(n),SNR=10dBX1信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為:AR系數(shù)：1-0.15954預(yù)測誤差的最小均方值：1.0612X2信號:度幅觀測信號X1(n),SNR=20dB度幅觀測信號X2(n),SNR=10dBX1信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為:AR系數(shù)：1-0.15954預(yù)測誤差的最小均方值：1.0612X2信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為:AR系數(shù)：1-0.1682預(yù)測誤差的最小均方值：1.1551X3信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為:AR系數(shù)：1-0.1161預(yù)測誤差的最小均方值：1.2883觀測信號X3(n),SNR=6dBAR系數(shù)：1-0.206580.25733預(yù)測誤差的最小均方值：0.98824X2信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)：1-0.105740.15188預(yù)測誤差的最小均方值：1.0349X3信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)：1-0.173760.23089預(yù)測誤差的最小均方值：1.2323N=5X1信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)：1-0.215870.22397-0.243060.24469-0.13453預(yù)測誤差的最小均方值：0.88869X2信號:N階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)：1-0.25370.31482-0.190140.122430.040983預(yù)測誤差的最小均方值：