现代信号实验ZY终

1、现代信号处理实验报告专业：模式识别与智能系统姓名：曾勇学号：2006193实验报告一 递推最小二乘法1、问题的提出当由实验提供了大量数据时，不能要求拟合函数在数据点处的偏差， 即(i=1,2,m) 严格为零，但为了使近似曲线尽量反映所给数据点的变化趋势，需对偏差有所要求。通常要求偏差平方和最小,此即称为最小二乘原理。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 )的值，并在计算机上实现该算法。（2）拟合模型的建立关于拟和模型必须能反映离散点分布基本特征。常选取j是线性拟和模型，既j所属函数类为M=Spanj，j1， jn，其中 j 0，j1， jn 是线性无关的基函数，于是j(x)= c j jj(x)。通常选取每个jj是次数£j的简单多项式，即M 是次数£ 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)设离散数据模型j(x)= c j jj

3、(x) 则求解归结为 n+1元函数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)称为正规方程组其中(jj ，j k )= wi jj(xi) j k(x i) 一般，n < m，函数 j 0，j1，jn，线性无关能保证正规方程组的系数矩阵的行列式不为零。因此正规方程组有唯一解。设其解为c j =c j *，j=0，1，n

4、则所要求的离散点的拟合函数(最佳平方逼近)为j*(x)= c j *jj(x)3、数据处理本实验的拟合函数为j(x)= c 2k+1 *x2k+1，其中j 0(x)=x, j 1(x)=x3, j 2(x)=x5, j 3(x)=x7。正规方程组为(jj ，j k) c j =(y, j k ), (k=0,1,2)，以及正规方程组的系数公式， 其中m=8，f即正规方程中的y，xi 和y i为已知的实验数据。由此可得到相应的系数c 2k+1 (k=0,1,2)，VC+仿真的结果如图1。图1 递推最小二乘法仿真结果4、结果分析VC+仿真的结果表明，拟合函数为j(x)= c 2k+1 *x2k+1

5、sinx，拟合函数近似为一个正弦函数，离散的数据基本在逼近正弦曲线。最小二乘法方曲线拟和是实验数据处理的常用方法。最佳平方逼近可以在一个区间上比正弦较均匀的逼近函数且具有方法简单易行，实效性大，应用广泛等特点。但当正规方程阶数较高时，往往出现病态。因此必须谨慎对待和加以巧妙处理。有效方法之一是引入正交多项式以改善其病态性。5、附录矩阵相乘程序： /* 计算矩阵的乘积 */1） void brmul(a,b,m,n,k,c,col1,col2) /* m*n维矩阵A乘以n*k维矩阵B 结果放在矩阵C里 col1为A所在大矩阵的列数 col2为B所在大矩阵的列数*/ 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) /*维数为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）主程序 /*循环部分*/for(k=4;k<8;k+)a100 =xk;a110 = xxk;a120 =xxxk;a130 =xxxxk;/向量akfor(i=0;i<4;i+)a20i=a1i0; ;/ a2是向量ak的转置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、 问题的提出 参数不随时间变化的系统称为时不

12、变系统。相当多的平稳随机过程都可以通过用白色噪声激励一线性时不变系统产生，而线性系统又可以用线性差分方程进行描述，这种差分模型就是自回归滑动平均（ARMA）模型。ARMA模型是描述平稳随机序列的最常用的一种模型，研究表明任何一个有理式的功率谱密度都是可以用一个ARMA随机过程的功率谱密度精确逼近的。2、 方案设计（1）实验要求x(n)是一个满足差分方程的ARMA过程。基本要求：将系统视为灰箱（可以利用方程的结构信息，p=3, q=2），根据系统的输入输出关系确定当e(n)为N(0,1)白噪声时输出功率谱密度的有理分式表达式。发挥部分：将系统视为黑箱（p,q未知），辨识出AR和MA阶数p,q。（

13、2）AR(p)模型参数的Yule-Walker估计(AR参数的可辨识性) 若ARMA(p,q)模型的多项式A(z)和B(z)无对消因子,且ap¹0，则该模型的参数可由以下p个修正Yule-Walker方程确定其中令，则修正Yule-Walker方程可以写成由定理可知，当ARMA(p,q)过程x(n)的AR阶数p和自相关函数Rx(t)已知时，只需求解p个Yule-Walker方程，便可辨识出AR参数。3、 数据处理本实验已知p=3，q=2，MA参数：b1=0.7，b2=0.12。先由高斯噪声函数带入已知函数得到xi的数据。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)模型参数的Yule-Walker估计仿真结果4、 结果分析求解系数与已知系数接近。当数据xi很少时，即N值很小时，误差比较大；当数据xi增多时，即N值增大时，误差明显减小；当N值到一定程度后，误差不再减小。AR参数为：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为方程组的系数矩阵，B为结果矩阵，n为方程组的个数 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)/以下是产生高斯白噪声的程序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);/产生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);/以上计算自相关函数值 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;得出系数矩阵 agaus(a,b,3);/计算出结果 printf("系数A为:n"); for(i=0;i<=2;i+) printf("%10.6fn",bi); getchar();实验报告三 Kalman滤波1、问题的提出卡尔曼滤波器是一个最优化自回归数据处理算法，实质上就是线性、无偏、最小均方误差的估计。这样估计的实际意义为：对于一输入未知的体系，可根据输出Y、体系特征和统计性质来最佳地估计输入。2、方案设计（1）实验要求,：均值为零方差为1的均匀分布白噪声；：均值为零方差为0.2的

22、均匀分布白噪声；相互独立。试用卡尔曼滤波算法估计。（2）基于状态空间模型的卡尔曼滤波设随机信号的模型方程为，它的测量方程为。利用含噪的观测数据，对时的状态向量各分量的最小均方误差估计。当时，Kalman算法为滤波问题。给定观测值的信息定义为其中为的一步最小均方误差预测，卡尔曼算法已知参数：状态转移矩阵观测矩阵过程噪声向量的相关矩阵观测噪声向量的相关矩阵，其中计算：初值：, ,n<N,n=n+1NY图1 kalman算法流程图3、数据处理由题目已知条件：,计算得，,又因：均值为零方差为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