數(shù)字信號處理實驗報告姓名：專業(yè)：通信與信息系統(tǒng)學(xué)號：日期：2015。11實驗內(nèi)容任務(wù)一：一連續(xù)平穩(wěn)的隨機(jī)信號x(t)，自相關(guān)函數(shù)r()=g-H，信號x(t)為加性噪聲所干擾，噪聲是白噪聲，測量x值的離散值Z4)為已知,T=0.02S,—3.2，—0。8，—14,-16，—17，-18，—3。3,—2。4，-18,-0。3，S—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è)計一維納濾波器。產(chǎn)生三組觀測數(shù)據(jù)：首先根據(jù)s(n)=as(n-1)+w(n)產(chǎn)生信號K)，將其加噪(信噪比分別為20dB，10dB,6dB)，得到觀測數(shù)據(jù)x(n)，x(n)，x(n).123估計x(n)，i=1,2,3的AR模型參數(shù)。假設(shè)信號長度為L,AR模型階數(shù)為N，分析實驗結(jié)果，并討論改變L，N對實驗結(jié)果的影響。實驗任務(wù)-1.卡爾曼濾波原理1。1卡爾曼濾波簡介早在20世紀(jì)40年代，開始有人用狀態(tài)變量模型來研究隨機(jī)過程，到60年代初，由于空間技術(shù)的發(fā)展，為了解決對非平穩(wěn)、多輸入輸出隨機(jī)序列的估計問題，卡爾曼提出了遞推最優(yōu)估計理論。它用狀態(tài)空間法描述系統(tǒng)，由狀態(tài)方程和量測方程所組成，即知道前一個狀態(tài)的估計值和最近一個觀測數(shù)據(jù)，采用遞推的算法估計當(dāng)前的狀態(tài)值。由于卡爾曼濾波采用遞推法，適合于計算機(jī)處理，并且可以用來處理多維和非平穩(wěn)隨機(jī)信號，現(xiàn)已廣泛應(yīng)用于很多領(lǐng)域并取得了很好的結(jié)果?？柭鼮V波一經(jīng)出現(xiàn)，就受到人們的很大重視并在實踐中不斷豐富和完善，其中一個成功的應(yīng)用是設(shè)計運(yùn)載體的高精度組合導(dǎo)航系統(tǒng)。卡爾曼濾波具有以下的特點，八*、?算法是遞推的，且狀態(tài)空間法采用在時域內(nèi)設(shè)計濾波器的方法，因而適用于多維隨機(jī)過程的估計;離散型卡爾曼算法適用于計算機(jī)處理。用遞推法計算，不需要知道全部過去的值，用狀態(tài)方程描述狀態(tài)變量的動態(tài)變化規(guī)律，因此信號可以是平穩(wěn)的，也可以是非平穩(wěn)的，即卡爾曼濾波適用于非平穩(wěn)過程.卡爾曼濾波采取的誤差準(zhǔn)則仍為估計誤差的均方值最小。1.2卡爾曼濾波的狀態(tài)方程和測量方程假設(shè)某系統(tǒng)k時刻的狀態(tài)變量為七，狀態(tài)方程和量測方程(輸出方程)表示為、=A_i氣_「％一1+匕其中，x是狀態(tài)變量;w表示輸入信號是白噪聲；V是觀測噪聲;J是觀測數(shù)據(jù).kk—1kk為了推導(dǎo)簡單，假設(shè)狀態(tài)變量的增益矩陣A不隨時間發(fā)生變化，wk，匕都是均值為零的正態(tài)白噪聲，方差分別是Q和人，并且初始狀態(tài)與w，V都不相關(guān)，Y表示相關(guān)系數(shù)。即:kkkkwkVk:^Iw]=0q2=Q,ywkVkkwkw^,w,kkj:e\v]=0,b2=R,y=R5kvkVk,v.kkj其中1。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)方程和量測方程為A4Axk=Akxk—1J'=Cx=CAx

kkkkkk—1顯然，由于不考慮觀測噪聲的影響，輸出信號的估計值與實際值是有誤差的，用Jk表示為了提高狀態(tài)估計的質(zhì)量，用輸出信號的估計誤差Jk來校正狀態(tài)變量fJ-CAxkkkk—1/V/V/VjTOC\o"1-5"\h\zf，\J—Jkk\J其中，H為增益矩陣，即加權(quán)矩陣。經(jīng)過校正后的狀態(tài)變量的估計誤差及其均方值分別用可和P表示，kkk把未經(jīng)校正的狀態(tài)變量的估計誤差的均方值用fJ-CAxkkkk—1/V/V/VjPk=Etx八Pk=Eik八Pk=E卡爾曼濾波要求狀態(tài)變量的估計誤差的均方值p為最小，因此卡爾曼濾波的關(guān)鍵即為通過選擇合適的kH，使得P取得最小值。首先推導(dǎo)狀態(tài)變量的估計值X和狀態(tài)變量的的估計誤差X，然后計算擊的均方值kkkP，通過化簡Pk，得到一組卡爾曼濾波的遞推公式：L入(AAx=Ax+Hy-CAxkkk-1ktkkkk-17<H=PCTCP'Ct+R)1kkkkkkkP=APAt+Qk(kk-1k、k-1、P=Q-hrp假設(shè)初始條件A,C，Q，R，y，X-，P已知，其中X。=E\x］，P=varL［那么遞推流程如下:f，，，，，kkkkkk-1k-10000x.，p—pgh—n，p

k-1k-1kkkk2.卡爾曼濾波遞推程序編程思想題目分析由于信號x(t)為加性噪聲所干擾，可知y=x+v，所以ck=1又因為噪聲為白噪聲，所以Rk2=s(o)=1⑶因為,所以XIj(AIIS\37=r\tn)z~m==em^-mxxxx由此可知,8=一—一，即K）-g-ixC-1）=代2-1）,可得到：A=w,因為抽樣間隔r=O.O2S,所以：z1-e-iz-isA=e—Tq=e-o.02。（4）因此x（n+1）_e-rsx（n）,所以=。2=rE)]=eLvG)ivG)]=1-e-2rwww因此Q=1—e-o.04k由上面的分析可知初始條件A,k編程分析Ck，",七已知，在仿真中假設(shè)丁。，則n二°，匕=由上面的分析可知初始條件A,kTOC\o"1-5"\h\zrc\x=e-o.02x+Hy—e-0.02xkk—1k/k—1/V\/VJ<H=pG+1)ikkkP=e-o.o4尸+1—e-o.04k,0P=\I-H)Pkk將得到的公式代入前面分析的遞推公式，即可進(jìn)行迭代得到結(jié)果1。k3.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..?！?1—14-0。90。8-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時由初值開始計算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);endend(洲矚wordMHLNHH0.02*M-%奮匝、畫EEx(i)SM?figure(1)POUT、Y、-:、-LinewaihLl)；x-abe-(d-y-abe-(-y3二2-_e(-7H少WSM—渲MmffiF血y(i)??二grid;figure(2)POUT、X、-b':unewaih-、1二x-abe-(T)-y-abe-(-x3)；2-_e(7H少WSM僉豐而血x(UM?)-grid;4.實驗結(jié)果卡爾曼濾波-估計信號x(t)波形t1510-5-10-15-20卡爾曼濾波測量信號z（t）

估計信號x（t）0實驗任務(wù)二1.維納濾波器原理維納-霍夫方程rd?)=尤h(m)r(k-m)=h(k)*r(k)m=0當(dāng)h（n）是一個長度為M的因果序列（即一個長度為M的FIR濾波器）時，維納-霍夫方程表述為r(k)=Mhm)r(k-m)=h(k)*r(k)k=0,1,2,m=0定義

3一「r(03一「r(0)]h=h2R=roXd:xd:hLM_r（M-1）J心）廣（0）=XXxx::XX??r（M-1）r（M-2）…r（r（M-1）rX（M-2）:〔（0）RX「5對上式求逆，得到h=R:七由以上式子可知:若已知期望信號與觀測數(shù)據(jù)的互相關(guān)函數(shù)及觀測數(shù)據(jù)的自相關(guān)函數(shù)則可以通過矩陣求逆運(yùn)算，得到維納濾波器的最佳解。同時可以看到，直接從時域求解因果的維納濾波器，當(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時刻的信號值，即估計?(n+N),N>0，估計得到的結(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)為龍0。根據(jù)現(xiàn)行系統(tǒng)的基本理論，輸出信號y（n）=x（n）=^jh（k）x（n-k）k=lpk，則SCi）=-如ox（n-k）pkpk

k=l預(yù)測誤差e（n）=xGi）-x（n）=x（n）+Xax（n—k）=x（n—k）pkk=lpkk=l其中a=1po要使均方誤差為最小值，要求又因為，我們可以得到要使均方誤差為最小值，要求又因為，我們可以得到所以r（/）+Y'ar（k-L）=0I=1,2,.o.,p（1）XXpkXXk=l

由于預(yù)測器的輸出X（n）是輸入信號的線性組合，所以可得:Ee*(n')x(n)以上說明誤差信號與輸入信號滿足正交性原理，預(yù)測誤差與預(yù)測信號值同樣滿足正交性原理。預(yù)測誤差的最小均方值EU)2]min=Ee*Cn[x(n)一x(n)=Ex*(EU)2]min=Ee*Cn[x(n)一x(n)=Ex*(n)+如ax*pkk=1=EL(n)x(n)](n-k)1x(n)（2）由（1）（2）聯(lián)立方程組，寫成矩陣形式可得「,(0)r0...r(p)]「1]Ee(n*」Exxr(0)...rip-1)a.min0xx...xx......xx(0)p1...=***r(p)xxr(p-1)...rxx(0)a1-pp」..._0_這就是有名的Yule-Walker維納一霍夫）方程。2.實驗編程思想在本實驗中，首先根據(jù)要求產(chǎn)生加噪不同的觀測數(shù)據(jù)X1（n）,x2（n），x（n），然后可利用已知條件代入Yule—Walker方程，即可求解AR模型參數(shù)。在本實驗中，假設(shè)a=0.2，信號s（n）的初值sG）=0.3.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/(10A(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('信號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;c。p_E(tes)而邪_>((fu)sn^e￥_E」(、q"s)3_d二ICXJCXJModqns鑒|1咿皿壬?%-(」$u=、-■'PEnssluLeNS、S)u6希HmxK」$u=」PEns$lu」以NS、S)U6MPHrxJx-(」$u=」PEns$lu"2NS、S)U6希HIX咿aw會k卅柜ft嗤呂％PUCD=U)M+(Tu)s%H(u)s蕾管平施蚩你真P」OM翩探)s+N+n)so」cdzH+N、I)SO」CDZHxmx)￡6ucd-h_i(N、XM<UOEUE二n)v<(c<n8?mx侵U/WUK』-(ncnxm<(c<7rik<r刪?巨贏*affi妄N痞?咿erxjx侵u/)tuK』_(N、IX注<(c<RIK<R刪?巨贏*affi妄N痞?咿皿艾侵u/)tuK』<￡%%-uopy何(?s)Qq受」(fu)mx咿-lllns會K￥_E二-qXX)lock(寸、rxJCXJModqnsCOP&二財?、)而qp_>」(、(5浸咿伽￡誨￥_-2二-qCXJx)lodsCN、rxj)lodqnscoP&二tes)而qp兵-(，艾咿空菸冥￥_專(q、IX)lo_d二rxJCNCN)lodqnsfori=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)]);

4.實驗結(jié)果與分析觀測數(shù)據(jù)產(chǎn)生度幅原始信號S(n)度幅圖1.原始信號與觀測信號(L=50)模型階數(shù)N對實驗結(jié)果的影響N=1對X1信號來說N階模型參數(shù)和誤差分別為：度幅原始信號S(n)度幅對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=3Ln8I17rn.?！?奪寸1。0ICNCNCNIoo—I一轍燎止<6uLn6oo-測購貝aLn8LnLnz。?！猺nI0"9.0ILnlzLnrn。?TH-\6尋。-c-測眼矣a179rn6Ioo—8CN60Ln0oo仁9。1。?I一轍燎止<蕾管平施蚩你真P」OM翩探)保Lnrn.IliwnH<矣a8oLnCNo——68CNS.0寸寸68"0。0RCNOLOrnlzrnRoo—TH一轍燎止<Lnolz"。1-測眼矣aRrnILnO0'寸LnCNw。0LnooizLn.。——6I99Ln.oLnILnLnrnoo—TH-gnN旨rn66.0-測眼矣aAR系數(shù)：1-0。365940O41414—0.416650.66894—0.60712均方誤差：1。1025分析：由以上實驗結(jié)果可知：在數(shù)據(jù)的長度一定的條件下改變AR模型的階數(shù)，均方誤差會改變，當(dāng)階數(shù)在某個值時，均方誤差的值最小，因此濾波器的階數(shù)對實驗結(jié)果有很大影響。在本次實驗中，仿真情況有限，在以上仿真中我們可以看到當(dāng)模型階數(shù)N為某一固定值時，均方誤差明顯較小。信號長度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=10dB觀測信號X3(n),SNR=6dBX1信號: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-N-srnx617rn0。1-88ILnI。0寸目。LolI-^<-n-8CNX寸"886.。迪貝岑「褊glwnH<g隱rnuLnCN.o8Ln9oCNo——I-^<-n-sxx"nN卷管平施蚩你真p」OMa{IR)6S868。0-依6。寸。.0rnwCN-ro寸1。61.?"孚云.0zrnLnCNOOI【巔㈱止<-n-s"X69