現(xiàn)代信號實驗ZY終

1、現(xiàn)代信號處理實驗報告專業(yè)：模式識別與智能系統(tǒng)姓名：曾勇學號：2006193實驗報告一 遞推最小二乘法1、問題的提出當由實驗提供了大量數(shù)據(jù)時，不能要求擬合函數(shù)在數(shù)據(jù)點處的偏差， 即(i=1,2,m) 嚴格為零，但為了使近似曲線盡量反映所給數(shù)據(jù)點的變化趨勢，需對偏差有所要求。通常要求偏差平方和最小,此即稱為最小二乘原理。2、方案設計（1）實驗要求已知，其中是均值為零方差相同的獨立隨機變量。觀測值如下表所示：xiyi0.1250.1246750.250.2474040.3750.3662730.50.4794260.6250.5850970.750.6816390.8750.76754410.841

2、471試用遞推最小二乘法估計cj ( j=1,3,5,7 )的值，并在計算機上實現(xiàn)該算法。（2）擬合模型的建立關于擬和模型必須能反映離散點分布基本特征。常選取j是線性擬和模型，既j所屬函數(shù)類為M=Spanj，j1， jn，其中 j 0，j1， jn 是線性無關的基函數(shù)，于是j(x)= c j jj(x)。通常選取每個jj是次數(shù)£j的簡單多項式，即M 是次數(shù)£ n 的n次多項式空間。取 jj(x)=x j ， j=0，1，n M =Span1 ,x , x2，x n，從而j(x)= C0 +C1 x1 + + C n x n =Pn(x)設離散數(shù)據(jù)模型j(x)= c j jj

3、(x) 則求解歸結為 n+1元函數(shù)S的 極值問題： S(c0，c1，c n)= wi y i -c j jj(xi) 2顯然S達最小值必要條件是=2wi y i -c j jj(xi) j k(x i)= 0,(k=0 ，1，n)這是關于 c0，c1，c n 的方程組，改寫成 (jj ，j k) c j =(y, j k ), (k=0,1,2,n)稱為正規(guī)方程組其中(jj ，j k )= wi jj(xi) j k(x i) 一般，n < m，函數(shù) j 0，j1，jn，線性無關能保證正規(guī)方程組的系數(shù)矩陣的行列式不為零。因此正規(guī)方程組有唯一解。設其解為c j =c j *，j=0，1，n

4、則所要求的離散點的擬合函數(shù)(最佳平方逼近)為j*(x)= c j *jj(x)3、數(shù)據(jù)處理本實驗的擬合函數(shù)為j(x)= c 2k+1 *x2k+1，其中j 0(x)=x, j 1(x)=x3, j 2(x)=x5, j 3(x)=x7。正規(guī)方程組為(jj ，j k) c j =(y, j k ), (k=0,1,2)，以及正規(guī)方程組的系數(shù)公式， 其中m=8，f即正規(guī)方程中的y，xi 和y i為已知的實驗數(shù)據(jù)。由此可得到相應的系數(shù)c 2k+1 (k=0,1,2)，VC+仿真的結果如圖1。圖1 遞推最小二乘法仿真結果4、結果分析VC+仿真的結果表明，擬合函數(shù)為j(x)= c 2k+1 *x2k+1

5、sinx，擬合函數(shù)近似為一個正弦函數(shù)，離散的數(shù)據(jù)基本在逼近正弦曲線。最小二乘法方曲線擬和是實驗數(shù)據(jù)處理的常用方法。最佳平方逼近可以在一個區(qū)間上比正弦較均勻的逼近函數(shù)且具有方法簡單易行，實效性大，應用廣泛等特點。但當正規(guī)方程階數(shù)較高時，往往出現(xiàn)病態(tài)。因此必須謹慎對待和加以巧妙處理。有效方法之一是引入正交多項式以改善其病態(tài)性。5、附錄矩陣相乘程序： /* 計算矩陣的乘積 */1） void brmul(a,b,m,n,k,c,col1,col2) /* m*n維矩陣A乘以n*k維矩陣B 結果放在矩陣C里 col1為A所在大矩陣的列數(shù) col2為B所在大矩陣的列數(shù)*/ int m,n,k,col1,

6、col2; double a,b,c; int i,j,l,u; for (i=0; i<=m-1; i+) for (j=0; j<=k-1; j+) u=i*col2+j; cu=0.0; for (l=0; l<=n-1; l+) cu=cu+ai*col1+l*bl*col2+j; return; 2） 矩陣求逆 /* 計算矩陣的逆 */ #include "stdlib.h" #include "math.h" #include "stdio.h" int brinv(a,n) /*維數(shù)為n*n矩陣A求逆后

7、結果仍放在A中*/ int n; double a; int *is,*js,i,j,k,l,u,v; double d,p; is=malloc(n*sizeof(int); js=malloc(n*sizeof(int); for (k=0; k<=n-1; k+) d=0.0; for (i=k; i<=n-1; i+) for (j=k; j<=n-1; j+) l=i*n+j; p=fabs(al); if (p>d) d=p; isk=i; jsk=j; if (d+1.0=1.0) free(is); free(js); printf("err*

8、not invn"); return(0); if (isk!=k) for (j=0; j<=n-1; j+) u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if (jsk!=k) for (i=0; i<=n-1; i+) u=i*n+k; v=i*n+jsk; p=au; au=av; av=p; l=k*n+k; al=1.0/al; for (j=0; j<=n-1; j+) if (j!=k) u=k*n+j; au=au*al; for (i=0; i<=n-1; i+) if (i!=k) for (j=0; j

9、<=n-1; j+) if (j!=k) u=i*n+j; au=au-ai*n+k*ak*n+j; for (i=0; i<=n-1; i+) if (i!=k) u=i*n+k; au=-au*al; for (k=n-1; k>=0; k-) if (jsk!=k) for (j=0; j<=n-1; j+) u=k*n+j; v=jsk*n+j; p=au; au=av; av=p; if (isk!=k) for (i=0; i<=n-1; i+) u=i*n+k; v=i*n+isk; p=au; au=av; av=p; free(is); free

10、(js); return(1); 3）主程序 /*循環(huán)部分*/for(k=4;k<8;k+)a100 =xk;a110 = xxk;a120 =xxxk;a130 =xxxxk;/向量akfor(i=0;i<4;i+)a20i=a1i0; ;/ a2是向量ak的轉(zhuǎn)置brmul(a2,p,1,4,4,e); brmul(e,a1,1,4,1,f);f00=f00+1;brinv(f,1);g=f00;brmul(p,a1,4,4,1,ff);brmul(ff,a2,4,1,4,l);brmul(l,p,4,4,4,ll);for(i=0;i<4;i+)for(j=0;j<

11、4;j+)llij=g*llij;for(i=0;i<4;i+)for(j=0;j<4;j+)pij=pij-llij;/以上計算P(k+1）f00=yk;brmul(aaa,b,4,k,1,ff);brmul(a1,f,4,1,1,fff);for(i=0;i<4;i+) ffi0=ffi0+fffi0;brmul(p,ff,4,4,1,cc); /此處計算for(j=0;j<4;j+)printf("%f",ccj0) /輸出結果printf("n");實驗報告二 ARMA譜估計1、 問題的提出 參數(shù)不隨時間變化的系統(tǒng)稱為時不

12、變系統(tǒng)。相當多的平穩(wěn)隨機過程都可以通過用白色噪聲激勵一線性時不變系統(tǒng)產(chǎn)生，而線性系統(tǒng)又可以用線性差分方程進行描述，這種差分模型就是自回歸滑動平均（ARMA）模型。ARMA模型是描述平穩(wěn)隨機序列的最常用的一種模型，研究表明任何一個有理式的功率譜密度都是可以用一個ARMA隨機過程的功率譜密度精確逼近的。2、 方案設計（1）實驗要求x(n)是一個滿足差分方程的ARMA過程?；疽螅簩⑾到y(tǒng)視為灰箱（可以利用方程的結構信息，p=3, q=2），根據(jù)系統(tǒng)的輸入輸出關系確定當e(n)為N(0,1)白噪聲時輸出功率譜密度的有理分式表達式。發(fā)揮部分：將系統(tǒng)視為黑箱（p,q未知），辨識出AR和MA階數(shù)p,q。（

13、2）AR(p)模型參數(shù)的Yule-Walker估計(AR參數(shù)的可辨識性) 若ARMA(p,q)模型的多項式A(z)和B(z)無對消因子,且ap¹0，則該模型的參數(shù)可由以下p個修正Yule-Walker方程確定其中令，則修正Yule-Walker方程可以寫成由定理可知，當ARMA(p,q)過程x(n)的AR階數(shù)p和自相關函數(shù)Rx(t)已知時，只需求解p個Yule-Walker方程，便可辨識出AR參數(shù)。3、 數(shù)據(jù)處理本實驗已知p=3，q=2，MA參數(shù)：b1=0.7，b2=0.12。先由高斯噪聲函數(shù)帶入已知函數(shù)得到xi的數(shù)據(jù)。xi=ei+0.7*ei-1+0.17*ei-2+0.2*xi-

14、1+0.23*xi-2-0.06*xi-3其中。由Yule-Walker方程：及公式，仿真過程中N取值為100000。解得a1=0.214856, a2=0.241973, a3=-0.049168。圖1 AR(p)模型參數(shù)的Yule-Walker估計仿真結果4、 結果分析求解系數(shù)與已知系數(shù)接近。當數(shù)據(jù)xi很少時，即N值很小時，誤差比較大；當數(shù)據(jù)xi增多時，即N值增大時，誤差明顯減??；當N值到一定程度后，誤差不再減小。AR參數(shù)為：a1=0.214856, a2=0.241973, a3=-0.049168。自相關矩陣為5、 附錄1) /* Agaus.c為求解線性方程組 */ #include

15、 "stdlib.h" #include "math.h" #include "stdio.h" int agaus(a,b,n) /A為方程組的系數(shù)矩陣，B為結果矩陣，n為方程組的個數(shù) int n; double a,b; int *js,l,k,i,j,is,p,q; double d,t; js=malloc(n*sizeof(int); l=1; for (k=0;k<=n-2;k+) d=0.0; for (i=k;i<=n-1;i+) for (j=k;j<=n-1;j+) t=fabs(ai*n+j);

16、 if (t>d) d=t; jsk=j; is=i; if (d+1.0=1.0) l=0; else if (jsk!=k) for (i=0;i<=n-1;i+) p=i*n+k; q=i*n+jsk; t=ap; ap=aq; aq=t; if (is!=k) for (j=k;j<=n-1;j+) p=k*n+j; q=is*n+j; t=ap; ap=aq; aq=t; t=bk; bk=bis; bis=t; if (l=0) free(js); printf("failn"); return(0); d=ak*n+k; for (j=k+1

17、;j<=n-1;j+) p=k*n+j; ap=ap/d; bk=bk/d; for (i=k+1;i<=n-1;i+) for (j=k+1;j<=n-1;j+) p=i*n+j; ap=ap-ai*n+k*ak*n+j; bi=bi-ai*n+k*bk; d=a(n-1)*n+n-1; if (fabs(d)+1.0=1.0) free(js); printf("failn"); return(0); bn-1=bn-1/d; for (i=n-2;i>=0;i-) t=0.0; for (j=i+1;j<=n-1;j+) t=t+ai*n

18、+j*bj; bi=bi-t; jsn-1=n-1; for (k=n-1;k>=0;k-) if (jsk!=k) t=bk; bk=bjsk; bjsk=t; free(js); return(1); 2)高斯白噪聲程序void random1(long n,float e,float mean,float var)/以下是產(chǎn)生高斯白噪聲的程序long j;float a,b,tt;srand(unsigned)time(NULL);for(j=0;j<n;j+) tt=rand(); if(tt-0.000001)<0) j-; continue; a=sqrt(-2.

19、0*log(tt/32767.0); b=2*3.14159265*rand()/32767.0; ej=var*a*cos(b)+mean; 3)主程序main() double Rx6=0; double a33,b3; long i,j; random1(N,e,0,1);/產(chǎn)生N個均值為0，方差為1的白噪聲 x0=e0; x1=e1+0.7*e0+0.2*x0; x2=e2+0.7*e1+0.17*e0+0.23*x0+0.2*x1; for(i=3;i<N;i+) xi=ei+0.7*ei-1+0.17*ei-2+0.2*xi-1+0.23*xi-2-0.06*xi-3;/得出

20、 X的值 for(i=0;i<(N-5);i+) Rx0=Rx0+xi*xi+0; Rx1=Rx1+xi*xi+1; Rx2=Rx2+xi*xi+2; Rx3=Rx3+xi*xi+3; Rx4=Rx4+xi*xi+4; Rx5=Rx5+xi*xi+5; Rx0=Rx0/(N-5); Rx1=Rx1/(N-5); Rx2=Rx2/(N-5); Rx3=Rx3/(N-5); Rx4=Rx4/(N-5); Rx5=Rx5/(N-5);/以上計算自相關函數(shù)值 b0=Rx3*(-1); b1=Rx4*(-1); b2=Rx5*(-1);/得出B的矩陣 for(i=0;i<3;i+) a0i=

21、Rx2-i; a1i=Rx3-i; a2i=Rx4-i;得出系數(shù)矩陣 agaus(a,b,3);/計算出結果 printf("系數(shù)A為:n"); for(i=0;i<=2;i+) printf("%10.6fn",bi); getchar();實驗報告三 Kalman濾波1、問題的提出卡爾曼濾波器是一個最優(yōu)化自回歸數(shù)據(jù)處理算法，實質(zhì)上就是線性、無偏、最小均方誤差的估計。這樣估計的實際意義為：對于一輸入未知的體系，可根據(jù)輸出Y、體系特征和統(tǒng)計性質(zhì)來最佳地估計輸入。2、方案設計（1）實驗要求,：均值為零方差為1的均勻分布白噪聲；：均值為零方差為0.2的

22、均勻分布白噪聲；相互獨立。試用卡爾曼濾波算法估計。（2）基于狀態(tài)空間模型的卡爾曼濾波設隨機信號的模型方程為，它的測量方程為。利用含噪的觀測數(shù)據(jù)，對時的狀態(tài)向量各分量的最小均方誤差估計。當時，Kalman算法為濾波問題。給定觀測值的信息定義為其中為的一步最小均方誤差預測，卡爾曼算法已知參數(shù)：狀態(tài)轉(zhuǎn)移矩陣觀測矩陣過程噪聲向量的相關矩陣觀測噪聲向量的相關矩陣，其中計算：初值：, ,n<N,n=n+1NY圖1 kalman算法流程圖3、數(shù)據(jù)處理由題目已知條件：,計算得，,又因：均值為零方差為1的均勻分布白噪聲；：均值為零方差為0.2的均勻分布白噪聲；相互獨立。則當時，和, 時，按方案設計中得遞推

23、公式計算。4、結果分析用VC+仿真Kalman濾波結果如圖2所示。圖中為理論均值，為估計均值。應用卡爾曼濾波可以達到很好得預測效果，濾波器是收斂的，理論上隨時間的推移，估計誤差會越來越小。圖2 kalman算法仿真結果5、附錄1）矩陣相乘程序：void brmul(a,b,m,n,k,c) /* m*n維矩陣A乘以n*k維矩陣B 結果放在矩陣C里 */ int m,n,k; double a,b,c; int i,j,l,u; for (i=0; i<=m-1; i+) for (j=0; j<=k-1; j+) u=i*k+j; cu=0.0; for (l=0; l<=n

24、-1; l+) cu=cu+ai*n+l*bl*k+j; return; 2）矩陣求逆程序程序同實驗一中的矩陣求逆程序。3）高斯白噪聲程序void random1(int nn,float n,float mean,float var)int j;float a,b,tt;srand(unsigned)time(NULL);for(j=0;j<nn;j+) tt=rand(); if(tt-0.000001)<0) j-; continue; a=sqrt(-2.0*log(tt/32767.0); b=2*3.14159265*rand()/32767.0; nj=var*a*c

25、os(b)+mean; 4）主程序 for(i=0;i<2;i+) for(j=0;j<3;j+) Ctji=Cij; for(i=0;i<3;i+) for(j=0;j<3;j+) F2ij=Fij; brinv(F2,3); for(i=0;i<3;i+) for(j=0;j<3;j+) Ftji=Fij; brmul(F,X,3,3,1,Xg); for(k=1;k<M;k+) random1(10,n,0,1); v100=0; v110=0; v120=n0*sin(0.1*k); random1(10,n,0,0.2); v200=n1; v210=n2; Q100=0;Q101=0;Q102=0; Q110=0;Q111=0;Q112=0; Q120=0;Q121=0;Q122=pow(sin(0.1*k),2); Q200=0.2;Q201=0; Q210=0;Q211=0.2; brmul(F,K,3,3,3,G1); brmul(G1,Ct,3,3,2,G2); brmul(C,K,2,3,3,G3); brmul(G3,Ct,2,3,2,G4); for(i=0;i<2;i+) for(j=0;j<2;j+) G4