實(shí)驗(yàn)報(bào)告卡爾曼濾波

數(shù)字信號(hào)解決實(shí)驗(yàn)報(bào)告姓名：專業(yè):通信與信息系統(tǒng)學(xué)號(hào):日期：.11實(shí)驗(yàn)內(nèi)容任務(wù)一：一持續(xù)平穩(wěn)的隨機(jī)信號(hào)，自有關(guān)函數(shù)，信號(hào)為加性噪聲所干擾，噪聲是白噪聲，測(cè)量值的離散值為已知，，-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，自編卡爾曼濾波遞推程序，預(yù)計(jì)信號(hào)的波形。任務(wù)二：設(shè)計(jì)一維納濾波器。產(chǎn)生三組觀察數(shù)據(jù)：首先根據(jù)產(chǎn)生信號(hào)，將其加噪（信噪比分別為20dB，10dB，6dB），得到觀察數(shù)據(jù)，，。預(yù)計(jì)，，，的AR模型參數(shù)。假設(shè)信號(hào)長(zhǎng)度為，AR模型階數(shù)為，分析實(shí)驗(yàn)成果，并討論變化，對(duì)實(shí)驗(yàn)成果的影響。實(shí)驗(yàn)任務(wù)一卡爾曼濾波原理1.1卡爾曼濾波介紹早在20世紀(jì)40年代，開(kāi)始有人用狀態(tài)變量模型來(lái)研究隨機(jī)過(guò)程，到60年代初，由于空間技術(shù)的發(fā)展，為理解決對(duì)非平穩(wěn)、多輸入輸出隨機(jī)序列的預(yù)計(jì)問(wèn)題，卡爾曼提出了遞推最優(yōu)預(yù)計(jì)理論。它用狀態(tài)空間法描述系統(tǒng)，由狀態(tài)方程和量測(cè)方程所構(gòu)成，即懂得前一種狀態(tài)的預(yù)計(jì)值和近來(lái)一種觀察數(shù)據(jù)，采用遞推的算法預(yù)計(jì)現(xiàn)在的狀態(tài)值。由于卡爾曼濾波采用遞推法，適合于計(jì)算機(jī)解決，并且能夠用來(lái)解決多維和非平穩(wěn)隨機(jī)信號(hào)，現(xiàn)已廣泛應(yīng)用于諸多領(lǐng)域，并獲得了較好的成果?？柭鼮V波一經(jīng)出現(xiàn)，就受到人們的很大重視，并在實(shí)踐中不停豐富和完善，其中一種成功的應(yīng)用是設(shè)計(jì)運(yùn)載體的高精度組合導(dǎo)航系統(tǒng)。卡爾曼濾波含有下列的特點(diǎn):算法是遞推的，且狀態(tài)空間法采用在時(shí)域內(nèi)設(shè)計(jì)濾波器的辦法，因而合用于多維隨機(jī)過(guò)程的預(yù)計(jì)；離散型卡爾曼算法合用于計(jì)算機(jī)解決。用遞推法計(jì)算，不需要懂得全部過(guò)去的值，用狀態(tài)方程描述狀態(tài)變量的動(dòng)態(tài)變化規(guī)律，因此信號(hào)能夠是平穩(wěn)的，也能夠是非平穩(wěn)的，即卡爾曼濾波合用于非平穩(wěn)過(guò)程。卡爾曼濾波采用的誤差準(zhǔn)則仍為預(yù)計(jì)誤差的均方值最小。1.2卡爾曼濾波的狀態(tài)方程和測(cè)量方程假設(shè)某系統(tǒng)時(shí)刻的狀態(tài)變量為，狀態(tài)方程和量測(cè)方程（輸出方程）表達(dá)為其中，是狀態(tài)變量；表達(dá)輸入信號(hào)是白噪聲；是觀察噪聲；是觀察數(shù)據(jù)。為了推導(dǎo)簡(jiǎn)樸，假設(shè)狀態(tài)變量的增益矩陣不隨時(shí)間發(fā)生變化，，都是均值為零的正態(tài)白噪聲，方差分別是和，并且初始狀態(tài)與，都不有關(guān)，表達(dá)有關(guān)系數(shù)。即：其中1.3卡爾曼濾波的遞推算法卡爾曼濾波采用遞推算法來(lái)實(shí)現(xiàn)，其基本思想是先不考慮輸入信號(hào)和觀察噪聲的影響，得到狀態(tài)變量和輸出信號(hào)（即觀察數(shù)據(jù)）的預(yù)計(jì)值，再用輸出信號(hào)的預(yù)計(jì)誤差加權(quán)后校正狀態(tài)變量的預(yù)計(jì)值，使?fàn)顟B(tài)變量預(yù)計(jì)誤差的均方值最小。因此，卡爾曼濾波器的核心是計(jì)算出加權(quán)矩陣的最佳值。當(dāng)不考慮觀察噪聲和輸入信號(hào)時(shí)，狀態(tài)方程和量測(cè)方程為顯然，由于不考慮觀察噪聲的影響，輸出信號(hào)的預(yù)計(jì)值與實(shí)際值是有誤差的，用表達(dá)為了提高狀態(tài)預(yù)計(jì)的質(zhì)量，用輸出信號(hào)的預(yù)計(jì)誤差來(lái)校正狀態(tài)變量其中，為增益矩陣，即加權(quán)矩陣。通過(guò)校正后的狀態(tài)變量的預(yù)計(jì)誤差及其均方值分別用和表達(dá)，把未經(jīng)校正的狀態(tài)變量的預(yù)計(jì)誤差的均方值用表達(dá)卡爾曼濾波規(guī)定狀態(tài)變量的預(yù)計(jì)誤差的均方值為最小，因此卡爾曼濾波的核心即為通過(guò)選擇適宜的，使得獲得最小值。首先推導(dǎo)狀態(tài)變量的預(yù)計(jì)值和狀態(tài)變量的的預(yù)計(jì)誤差，然后計(jì)算的均方值，通過(guò)化簡(jiǎn)，得到一組卡爾曼濾波的遞推公式：假設(shè)初始條件，，，，，，已知，其中，，那么遞推流程以下：，，卡爾曼濾波遞推程序編程思想題目分析由于信號(hào)為加性噪聲所干擾，可知，因此又由于噪聲為白噪聲，因此由于，因此由此可知，，即，可得到：，由于抽樣間隔，因此：。（4）因此，因此因此編程分析由上面的分析可知初始條件，，，，已知，在仿真中假設(shè)，則，，由以上參數(shù)可得卡爾曼實(shí)際遞推公式將得到的公式代入前面分析的遞推公式，即可進(jìn)行迭代得到成果。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ù)長(zhǎng)度N=length(Y);fork=1:Nifk==1%k=1時(shí)由初值開(kāi)始計(jì)算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時(shí)，開(kāi)始遞推%遞推公式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;%作圖，畫(huà)出x(t)的波形figure(1)plot(T,Y,'r','LineWidth',1);xlabel('t');ylabel('y(t)');title('卡爾曼濾波-測(cè)量信號(hào)y(t)波形');grid;figure(2)plot(T,X,'b','LineWidth',1);xlabel('t');ylabel('x(t)');title('卡爾曼濾波-預(yù)計(jì)信號(hào)x(t)波形');grid;實(shí)驗(yàn)成果實(shí)驗(yàn)任務(wù)二維納濾波器原理維納-霍夫方程當(dāng)是一種長(zhǎng)度為的因果序列（即一種長(zhǎng)度為的濾波器）時(shí)，維納-霍夫方程表述為定義則可寫成矩陣的形式，即對(duì)上式求逆，得到由以上式子可知：若已知盼望信號(hào)與觀察數(shù)據(jù)的互有關(guān)函數(shù)及觀察數(shù)據(jù)的自有關(guān)函數(shù)，則能夠通過(guò)矩陣求逆運(yùn)算，得到維納濾波器的最佳解。同時(shí)能夠看到，直接從時(shí)域求解因果的維納濾波器，當(dāng)選擇的濾波器的長(zhǎng)度較大時(shí)，計(jì)算工作量很大，并且需要計(jì)算的逆矩陣，從而規(guī)定的存儲(chǔ)量也很大預(yù)測(cè)是根據(jù)觀察到對(duì)的過(guò)去數(shù)據(jù)來(lái)預(yù)計(jì)現(xiàn)在或?qū)?lái)的信號(hào)值。維納預(yù)測(cè)是已知以前時(shí)刻的個(gè)數(shù)據(jù)，預(yù)計(jì)現(xiàn)在時(shí)刻，或者將來(lái)時(shí)刻的信號(hào)值，即預(yù)計(jì)，預(yù)計(jì)得到的成果仍然規(guī)定滿足均方誤差最小的準(zhǔn)則。信號(hào)能夠預(yù)測(cè)是由于信號(hào)內(nèi)部存在著關(guān)聯(lián)性。預(yù)測(cè)是運(yùn)用數(shù)據(jù)前后的關(guān)聯(lián)性，根據(jù)其中一部分推知其它部分。一步線性預(yù)測(cè)的時(shí)域解已知，，…，，預(yù)測(cè)，假設(shè)噪聲，這樣的預(yù)測(cè)成為一步線性預(yù)測(cè)。設(shè)定系統(tǒng)的單位脈沖響應(yīng)為。根據(jù)現(xiàn)行系統(tǒng)的基本理論，輸出信號(hào)令，則預(yù)測(cè)誤差其中要使均方誤差為最小值，規(guī)定，，...，又由于，我們能夠得到，，...，因此，，...，（1）由于預(yù)測(cè)器的輸出是輸入信號(hào)的線性組合，因此可得：以上闡明誤差信號(hào)與輸入信號(hào)滿足正交性原理，預(yù)測(cè)誤差與預(yù)測(cè)信號(hào)值同樣滿足正交性原理。預(yù)測(cè)誤差的最小均方值（2）由（1）（2）聯(lián)立方程組，寫成矩陣形式可得這就是有名的Yule-Walker（維納-霍夫）方程。實(shí)驗(yàn)編程思想在本實(shí)驗(yàn)中，首先根據(jù)規(guī)定產(chǎn)生加噪不同的觀察數(shù)據(jù)，，，然后可運(yùn)用已知條件代入Yule-Walker方程，即可求解AR模型參數(shù)。在本實(shí)驗(yàn)中，假設(shè)，信號(hào)的初值。MATLAB代碼functionWiener_predict(L,N)%clc;clear;%信噪比SN1=6;SN2=10;SN3=20;%產(chǎn)生信號(hào)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)生觀察信號(hào)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('norm',0,P3,L,1);forn=1:LX1(n)=S(n)+V1(n);X2(n)=S(n)+V2(n);X3(n)=S(n)+V3(n);endsubplot(2,2,1);plot(S,'b');title('信號(hào)S(n)');ylabel('幅度');gridon;subplot(2,2,2);plot(X1,'b');title('觀察信號(hào)X1(n)');ylabel('幅度');gridon;subplot(2,2,3);plot(X2,'b');title('觀察信號(hào)X2(n)');ylabel('幅度');gridon;subplot(2,2,4);plot(X3,'b');title('觀察信號(hào)X3(n)');ylabel('幅度');gridon;fprintf('\n對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：\n')AR(X1,N);fprintf('\n對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：\n')AR(X2,N);fprintf('\n對(duì)X3信號(hào)來(lái)說(shuō)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)生信號(hào)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)生觀察限號(hào)X1=awgn(S,SNR1,'measured','linear');X2=awgn(S,SNR2,'measured','linear');X3=awgn(S,SNR3,'measured','linear');%畫(huà)出信號(hào)圖像subplot(2,2,1);plot(S,'b');title('信號(hào)S(n)');ylabel('幅度');gridon;subplot(2,2,2);plot(X1,'b');title('觀察信號(hào)X1(n)');ylabel('幅度');gridon;subplot(2,2,3);plot(X2,'b');title('觀察信號(hào)X2(n)');ylabel('幅度');gridon;subplot(2,2,4);plot(X3,'b');title('觀察信號(hào)X3(n)');ylabel('幅度');gridon;%%預(yù)計(jì)ＡＲ模型參數(shù)fprintf('\n對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：\n');AR(X1,N);fprintf('\n對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：\n');AR(X2,N);fprintf('\n對(duì)X3信號(hào)來(lái)說(shuō)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)]);實(shí)驗(yàn)成果與分析觀察數(shù)據(jù)產(chǎn)生圖1.原始信號(hào)與觀察信號(hào)（L=50）模型階數(shù)N對(duì)實(shí)驗(yàn)成果的影響N=1對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.27766均方誤差：1.1289對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.29326均方誤差：0.97283對(duì)X3信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.26441均方誤差：1.0531N=2對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.344940.2854均方誤差：1.0958對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.16960.10742均方誤差：1.1639對(duì)X3信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.195320.17033均方誤差：0.92331N=3對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.106730.050928-0.19364均方誤差：1.4197對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.354510.6-0.75585均方誤差：0.95739對(duì)X3信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.122210.14428-0.34185均方誤差：0.99317N=5對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.355150.56619-0.540050.65254-0.51327均方誤差：1.2405對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.273430.102270.0289440.21289-0.2508均方誤差：1.3557對(duì)X3信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.365940.41414-0.416650.66894-0.60712均方誤差：1.1025分析：由以上實(shí)驗(yàn)成果可知：在數(shù)據(jù)的長(zhǎng)度一定的條件下，變化AR模型的階數(shù)，均方誤差會(huì)變化，當(dāng)階數(shù)在某個(gè)值時(shí)，均方誤差的值最小，因此濾波器的階數(shù)對(duì)實(shí)驗(yàn)成果有很大影響。在本次實(shí)驗(yàn)中，仿真狀況有限，在以上仿真中我們能夠看到當(dāng)模型階數(shù)N為某一固定值時(shí)，均方誤差明顯較小。信號(hào)長(zhǎng)度L對(duì)實(shí)驗(yàn)成果的影響L=100對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：10.0229140.0078698均方誤差：1.2033對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：10.0174140.19629均方誤差：1.1607對(duì)X3信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.0127890.13086均方誤差：1.1483L=200對(duì)X1信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.176790.073726均方誤差：1.3371對(duì)X2信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.26490.16751均方誤差：0.99844對(duì)X3信號(hào)來(lái)說(shuō)N階模型參數(shù)和誤差分別為：AR系數(shù)：1-0.271450.17666均方誤差：0.99289分析：由以上仿真成果可知，實(shí)驗(yàn)中存在誤差，但仍然能夠看出，隨著信號(hào)長(zhǎng)度的增加，均方誤差減小，預(yù)測(cè)更精確。L=100，N=1X1信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.15954預(yù)測(cè)誤差的最小均方值：1.0612X2信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.1682預(yù)測(cè)誤差的最小均方值：1.1551X3信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.1161預(yù)測(cè)誤差的最小均方值：1.2883N=2X1信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.206580.25733預(yù)測(cè)誤差的最小均方值：0.98824X2信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.105740.15188預(yù)測(cè)誤差的最小均方值：1.0349X3信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.173760.23089預(yù)測(cè)誤差的最小均方值：1.2323N=5X1信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.215870.22397-0.243060.24469-0.13453預(yù)測(cè)誤差的最小均方值：0.88869X2信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.25370.31482-0.190140.122430.040983預(yù)測(cè)誤差的最小均方值：0.89859X3信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.278120.33384-0.278810.25752-0.11447預(yù)測(cè)誤差的最小均方值：1.0384N=10X1信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.0857460.17042-0.248870.3049-0.290130.34871-0.331290.48704-0.399420.37265預(yù)測(cè)誤差的最小均方值：0.96799X2信號(hào):N階模型參數(shù)和預(yù)測(cè)誤差的最小均方值分別為：AR系數(shù)：1-0.101830.20689-0.18967