實驗報告-卡爾曼濾波

1、精選優(yōu)質(zhì)文檔-傾情為你奉上數(shù)字信號處理實驗報告姓名： 專業(yè): 通信與信息系統(tǒng) 學(xué)號: 日期：2015.11實驗內(nèi)容任務(wù)一：一連續(xù)平穩(wěn)的隨機信號，自相關(guān)函數(shù)，信號為加性噪聲所干擾，噪聲是白噪聲，測量值的離散值為已知，-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，自編卡爾曼濾波遞推程序

2、，估計信號的波形。任務(wù)二：設(shè)計一維納濾波器。（1） 產(chǎn)生三組觀測數(shù)據(jù)：首先根據(jù)產(chǎn)生信號，將其加噪（信噪比分別為20dB，10dB，6dB），得到觀測數(shù)據(jù)，。（2） 估計，的AR模型參數(shù)。假設(shè)信號長度為，AR模型階數(shù)為，分析實驗結(jié)果，并討論改變，對實驗結(jié)果的影響。實驗任務(wù)一1. 卡爾曼濾波原理1.1 卡爾曼濾波簡介早在20世紀(jì)40年代，開始有人用狀態(tài)變量模型來研究隨機過程，到60年代初，由于空間技術(shù)的發(fā)展，為了解決對非平穩(wěn)、多輸入輸出隨機序列的估計問題，卡爾曼提出了遞推最優(yōu)估計理論。它用狀態(tài)空間法描述系統(tǒng)，由狀態(tài)方程和量測方程所組成，即知道前一個狀態(tài)的估計值和最近一個觀測數(shù)據(jù)，采用遞推的算法估計

3、當(dāng)前的狀態(tài)值。由于卡爾曼濾波采用遞推法，適合于計算機處理，并且可以用來處理多維和非平穩(wěn)隨機信號，現(xiàn)已廣泛應(yīng)用于很多領(lǐng)域，并取得了很好的結(jié)果。卡爾曼濾波一經(jīng)出現(xiàn)，就受到人們的很大重視，并 在實踐中不斷豐富和完善，其中一個成功的應(yīng)用是設(shè)計運載體的高精度組合導(dǎo)航系統(tǒng)?？柭鼮V波具有以下的特點:（1） 算法是遞推的，且狀態(tài)空間法采用在時域內(nèi)設(shè)計濾波器的方法，因而適用于多維隨機過程的估計；離散型卡爾曼算法適用于計算機處理。（2） 用遞推法計算，不需要知道全部過去的值，用狀態(tài)方程描述狀態(tài)變量的動態(tài)變化規(guī)律，因此信號可以是平穩(wěn)的，也可以是非平穩(wěn)的，即卡爾曼濾波適用于非平穩(wěn)過程。（3） 卡爾曼濾波采取的誤差準(zhǔn)

4、則仍為估計誤差的均方值最小。1.2 卡爾曼濾波的狀態(tài)方程和測量方程假設(shè)某系統(tǒng)時刻的狀態(tài)變量為，狀態(tài)方程和量測方程（輸出方程）表示為其中，是狀態(tài)變量；表示輸入信號是白噪聲；是觀測噪聲；是觀測數(shù)據(jù)。為了推導(dǎo)簡單，假設(shè)狀態(tài)變量的增益矩陣不隨時間發(fā)生變化，都是均值為零的正態(tài)白噪聲，方差分別是和，并且初始狀態(tài)與，都不相關(guān)，表示相關(guān)系數(shù)。即：其中1.3 卡爾曼濾波的遞推算法卡爾曼濾波采用遞推算法來實現(xiàn)，其基本思想是先不考慮輸入信號和觀測噪聲的影響，得到狀態(tài)變量和輸出信號（即觀測數(shù)據(jù)）的估計值，再用輸出信號的估計誤差加權(quán)后校正狀態(tài)變量的估計值，使?fàn)顟B(tài)變量估計誤差的均方值最小。因此，卡爾曼濾波器的關(guān)鍵是計算出

5、加權(quán)矩陣的最佳值。當(dāng)不考慮觀測噪聲和輸入信號時，狀態(tài)方程和量測方程為顯然，由于不考慮觀測噪聲的影響，輸出信號的估計值與實際值是有誤差的，用表示為了提高狀態(tài)估計的質(zhì)量，用輸出信號的估計誤差來校正狀態(tài)變量其中，為增益矩陣，即加權(quán)矩陣。經(jīng)過校正后的狀態(tài)變量的估計誤差及其均方值分別用和表示，把未經(jīng)校正的狀態(tài)變量的估計誤差的均方值用表示卡爾曼濾波要求狀態(tài)變量的估計誤差的均方值為最小，因此卡爾曼濾波的關(guān)鍵即為通過選擇合適的，使得取得最小值。首先推導(dǎo)狀態(tài)變量的估計值和狀態(tài)變量的的估計誤差，然后計算的均方值，通過化簡，得到一組卡爾曼濾波的遞推公式：假設(shè)初始條件，已知，其中，那么遞推流程如下：，2. 卡爾曼濾波

6、遞推程序編程思想題目分析（1） 由于信號為加性噪聲所干擾，可知，所以（2） 又因為噪聲為白噪聲，所以（3） 因為，所以由此可知，即，可得到：，因為抽樣間隔，所以：。 （4）因此，所以因此編程分析由上面的分析可知初始條件，已知，在仿真中假設(shè)，則，由以上參數(shù)可得卡爾曼實際遞推公式將得到的公式代入前面分析的遞推公式，即可進行迭代得到結(jié)果。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

7、);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.9 0.8 10 0.2 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;%數(shù)據(jù)長度N=length(Y);for k=1:N if k=1 %k=1時由初值開始計算 P_(k)=Ak*P0*Ak'+Qk; H(k)=P_(k)*Ck

8、9;*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

9、*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.

10、 維納濾波器原理維納-霍夫方程當(dāng)是一個長度為的因果序列（即一個長度為的濾波器）時，維納-霍夫方程表述為定義則可寫成矩陣的形式，即對上式求逆，得到由以上式子可知：若已知期望信號與觀測數(shù)據(jù)的互相關(guān)函數(shù)及觀測數(shù)據(jù)的自相關(guān)函數(shù)，則可以通過矩陣求逆運算，得到維納濾波器的最佳解。同時可以看到，直接從時域求解因果的維納濾波器，當(dāng)選擇的濾波器的長度較大時，計算工作量很大，并且需要計算的逆矩陣，從而要求的存儲量也很大預(yù)測是根據(jù)觀測到對的過去數(shù)據(jù)來估計當(dāng)前或?qū)淼男盘栔?。維納預(yù)測是已知以前時刻的個數(shù)據(jù)，估計當(dāng)前時刻，或者未來時刻的信號值，即估計，估計得到的結(jié)果仍然要求滿足均方誤差最小的準(zhǔn)則。信號可以預(yù)測是由于信號

11、內(nèi)部存在著關(guān)聯(lián)性。預(yù)測是利用數(shù)據(jù)前后的關(guān)聯(lián)性，根據(jù)其中一部分推知其余部分。一步線性預(yù)測的時域解已知，預(yù)測，假設(shè)噪聲，這樣的預(yù)測成為一步線性預(yù)測。設(shè)定系統(tǒng)的單位脈沖響應(yīng)為。根據(jù)現(xiàn)行系統(tǒng)的基本理論，輸出信號令，則預(yù)測誤差其中要使均方誤差為最小值，要求，.，又因為，我們可以得到，.，所以，.， （1）由于預(yù)測器的輸出是輸入信號的線性組合，所以可得：以上說明誤差信號與輸入信號滿足正交性原理，預(yù)測誤差與預(yù)測信號值同樣滿足正交性原理。預(yù)測誤差的最小均方值 （2）由（1）（2）聯(lián)立方程組，寫成矩陣形式可得這就是有名的Yule-Walker（維納-霍夫）方程。2. 實驗編程思想在本實驗中，首先根據(jù)要求產(chǎn)生加噪

12、不同的觀測數(shù)據(jù)，然后可利用已知條件代入Yule-Walker方程，即可求解AR模型參數(shù)。在本實驗中，假設(shè)，信號的初值。3. MATLAB代碼function Wiener_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;for n = 2 : L S(n) = a * S(n-1) + W(n); end%產(chǎn)生觀測信號Am = sum(abs(S).2) / L;P1 = Am / (10(SN1/20);

13、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);for n=1:L X1(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('信號S(n)

14、9;);ylabel('幅度');grid on;subplot(2,2,2);plot(X1,'b');title('觀測信號X1(n)');ylabel('幅度');grid on;subplot(2,2,3);plot(X2,'b');title('觀測信號X2(n)');ylabel('幅度');grid on;subplot(2,2,4);plot(X3,'b');title('觀測信號X3(n)');ylabel('幅度')

15、;grid 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); function AR(X,N) L = length(X);rx = zeros(1, N + 1);R = zeros(N + 1, N + 1);for i = 1 : (N + 1) sum = 0; for j = 1 : (L - i + 1); su

16、m = sum + X(j) * X(j + i - 1); end rx(i) = sum / (L - i + 1);endfor i = 1 : N + 1 R(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);funct

17、ion Wiener_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;for n = 2 : L S(n) = a * S(n-1) + W(n); end%加噪聲產(chǎn)生觀測限號X1= awgn(S,SNR1,'measured','linear'); X2= awgn(S,SNR2,'measured','linear'); X3= awg

18、n(S,SNR3,'measured','linear'); %畫出信號圖像subplot(2,2,1);plot(S,'b'); title('信號S(n)'); ylabel('幅度');grid on;subplot(2,2,2);plot(X1,'b');title('觀測信號X1(n)');ylabel('幅度');grid on;subplot(2,2,3);plot(X2,'b');title('觀測信號X2(n)');y

19、label('幅度');grid on;subplot(2,2,4);plot(X3,'b');title('觀測信號X3(n)');ylabel('幅度');grid on;% 估計模型參數(shù)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); function A

20、R(X,N) L = length(X);rx = zeros(1, N + 1);R = zeros(N + 1, N + 1);for i = 1 : (N + 1) sum = 0; for j = 1 : (L - i + 1); sum = sum + X(j) * X(j + i - 1); end rx(i) = sum / (L - i + 1);endfor i = 1 : N + 1 R(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,

21、 1:N) * (-zx)'a = 1, ap'e = rx(1) + zx * ap;disp('AR系數(shù)： ',num2str(a);disp('均方誤差：',num2str(e);4. 實驗結(jié)果與分析（1） 觀測數(shù)據(jù)產(chǎn)生圖1. 原始信號與觀測信號（L=50）（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

22、信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.26441均方誤差：1.0531N=2對X1信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.34494 0.2854均方誤差：1.0958對X2信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.1696 0.10742均方誤差：1.1639對X3信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.19532 0.17033均方誤差：0.92331N=3對X1信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.10673 0. -0.19364均方誤差：1.4197對X2信號來說 N 階模型參數(shù)和誤

23、差分別為：AR系數(shù)： 1 -0.35451 0.62013 -0.75585均方誤差：0.95739對X3信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.12221 0.14428 -0.34185均方誤差：0.99317N=5對X1信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.35515 0.56619 -0.54005 0.65254 -0.51327均方誤差：1.2405對X2信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.27343 0.10227 0. 0.21289 -0.2508均方誤差：1.3557對X3信號來說 N 階模型參數(shù)和誤差分別為

24、：AR系數(shù)： 1 -0.36594 0.41414 -0.41665 0.66894 -0.60712均方誤差：1.1025分析：由以上實驗結(jié)果可知：在數(shù)據(jù)的長度一定的條件下，改變AR模型的階數(shù)，均方誤差會改變，當(dāng)階數(shù)在某個值時，均方誤差的值最小，因此濾波器的階數(shù)對實驗結(jié)果有很大影響。在本次實驗中，仿真情況有限，在以上仿真中我們可以看到當(dāng)模型階數(shù)N為某一固定值時，均方誤差明顯較小。（3） 信號長度L對實驗結(jié)果的影響L=100對X1信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 0. 0.均方誤差：1.2033對X2信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 0. 0.19629

25、均方誤差：1.1607對X3信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0. 0.13086均方誤差：1.1483L=200對X1信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.17679 0.均方誤差：1.3371對X2信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.2649 0.16751均方誤差：0.99844對X3信號來說 N 階模型參數(shù)和誤差分別為：AR系數(shù)： 1 -0.27145 0.17666均方誤差：0.99289分析：由以上仿真結(jié)果可知，實驗中存在誤差，但仍然可以看出，隨著信號長度的增加，均方誤差減小，預(yù)測更準(zhǔn)確。L=100，N=1X1信號

26、: 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.2883N=2X1信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0.20658 0.25733預(yù)測誤差的最小均方值：0.98824X2信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0.10574 0.151

27、88預(yù)測誤差的最小均方值：1.0349X3信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0.17376 0.23089預(yù)測誤差的最小均方值：1.2323N=5X1信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0.21587 0.22397 -0.24306 0.24469 -0.13453預(yù)測誤差的最小均方值：0.88869X2信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0.2537 0.31482 -0.19014 0.12243 0.預(yù)測誤差的最小均方值：0.89859X3信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0.27812 0.33384 -0.27881 0.25752 -0.11447預(yù)測誤差的最小均方值：1.0384N=10X1信號: N 階模型參數(shù)和預(yù)測誤差的最小均方值分別為：AR系數(shù)： 1 -0. 0.17042 -0.24887 0.3