RAMA谱估计算法

1、RAMA谱估计算法 姓名：刘金强 专业：控制理论与控制工程 学号：2007255 实验目的：（1）、掌握RAMA谱估计算法的原理和实现方法（2）、了解Yule-Walker的原理（3）、熟练掌握Yule-Walker的实现方法实验要求： x(n)是一个满足差分方程： x(n)=0.2*x(n-1)+0.23*x(n-2)-0.06*x(n-3)+e(n)+0.7*e(n-1)+0.12*e(n-2)的RAMA过程。基本要求： 将系统视为黑箱(可以利用方程的结构信息，p=3,q=2),根据系统的输入输出关系确定，当e(n)为N(0,1)白噪声时AR参数。实验原理： 在Cadzow谱估计子和Kav

2、eh谱估计子中，都需要知道AR阶数和AR参数。在实际场合中，它们可以根据观测数据进行估计。为此，需要推导AR参数所服从的线性方程组。这一问题的答案由Gersch在1970年得到了解决。得到了我们现在的RA参数可辨识性的定理：若ARMA（p,q）模型的多项式A(z)和B(z)无可对消的公共因子，且0,该模型的AR参数可由p个修正Yule-Walker方程惟一确定或标识。由定理可以看出：只要已知了ARMA过程的结构参数p和q，同时已知自相关函数，AR参数可以由p个Yule-Walker方程唯一确定，从而，实现了AR参数的估计，并且节约大量的运算时间。实验方案流程：开始生成高斯白噪声序列结束生成x序

3、列的前三个序列估计AR参数生成完整的x序列构造Yule-Walker方程的系数矩阵和常数矩阵计算自相关函数Rx（1）、高斯白噪声的生成可以采用线性同余法、接收拒绝法等方法来实现，因为这一点不是本实验的重点内容，因此，此处不再赘述。（2）、x序列的生成由白噪声序列先生成x序列的前三个序列，即x0 = e0;x1 = 0.2 * x0 + e1 + 0.7 * e0;x2 = 0.2 * x1 + 0.23 * x0 + e2 + 0.7 * e1 + 0.12 * e0;然后，依据已有的白噪声序列和x序列（前三项）求出x序列的其他元素值，只需将已知的内容带入到关系式：中即可求出完整的x序列。（3

4、）、构造Yule-Walker方程的系数矩阵A和常数矩阵b由Yule-Walker方程和题目中的已知条件，即可以得到A和b如下：（4）、估计AR参数只需求解方程，即可求解出AR参数X序列。实验结果：AR的参数为：0.217062，0.216432，-0.0571144.实验分析实验正确估计出了AR参数，而且随着x序列长度的不断增长，估计的误差越来越小，但到达一定程度后，AR参数值基本保持稳定（不再变化或变化很小）。为了能够最优的估计出AR参数，同时兼顾程序的处理速度，序列长度的选择一定要适当。本程序得出的稳定结果是在序列长度为10000000时得出来的。程序代码：（1）主程序：/* 问题： x

5、(n)是一个满足差分方程 x(n)=0.2*x(n-1)+0.23*x(n-2)-0.06*x(n-3)+e(n)+0.7*e(n-1)+0.12*e(n-2)的RAMA过程。基本要求： 将系统视为黑箱(可以利用方程的结构信息，p=3,q=2),根据系统的输入输出关系确定，当e(n)为N(0,1)白噪声时AR参数。*设计作者：刘金强设计时间：2008年10月31日作者专业：控制理论与控制工程作者学号：2007255*/#include <iostream.h>#include "MathMatrix.h"#define NRANSI#define SWAP(a,

6、b) temp=(a);(a)=(b);(b)=temp;const int NUM=10000000;const int N =3;const int M =6;void GaussInverse(double *a, int n);void CMathMatrix:Mat_Vector(double *Arr, int Rows, int Cols, double B, double C);void main()int i,j,k;double *e = new double NUM; /高斯白噪声申请内存并初始化double *x = new double NUM; /x序列double

7、 *A = new double * N; /系数矩阵Adouble *B = new double N; /常数矩阵Bfor( i = 0; i < N; i+ )Ai = new double N;/生成高斯白噪声CMathMatrix MM;MM.random(NUM, e, 0, 1);/构造X序列x0 = e0;x1 = 0.2 * x0 + e1 +0.7 * e0;x2 = 0.2 * x1 + 0.23 * x0 + e2 + 0.7 * e1 + 0.12 * e0;/生成X序列for(i = 3; i < NUM; i+)xi =0.2 * xi-1 + 0.2

8、3 * xi-2 -0.06 * xi-3 + ei+0.7 * ei-1 + 0.12 * ei-2;/生成自相关函数double *R = new doubleM;for(i = 0; i < M; i+)Ri = 0;k = 0;for(i = 0; i < M; i+)for(j = 0; j < NUM - M - 1; j+)Ri = Ri + xj * xj + k; k+;Ri = Ri / (NUM - M);/生成系数矩阵A和常数矩阵Bk = 0;for(i = M - 1; i > M - N - 1; i-)Bk = -Ri;for (j = 1

9、; j < M - N + 1; j+)AM - i -1j - 1 = Ri - j; k+;/计算估计值double *X = new double N;GaussInverse(A, N); /计算A的逆A_1MM.Mat_Vector(A, N, N, B, X);/输出估计值for(i = 0; i < N; i+)cout<<Xi<<endl;cout<<endl;/释放内存for(i = 0; i < N; i+) delete Ai;delete A;delete e;delete x;delete X;delete B;d

10、elete R;（2）CMathMatrix文件的源文件CMathMatrix.cpp：/ mathMatrix.cpp: implementation of the CmathMatrix class./#include "mathMatrix.h"#include <time.h>#include <math.h>#include <stdio.h>#include <stdlib.h>#include <iostream.h>#define SWAP(a,b) temp=(a);(a)=(b);(b)=tem

11、p;/ Construction/Destruction/CMathMatrix:CMathMatrix()CMathMatrix:CMathMatrix()/矩阵求逆int *ivector(long nl, long nh);void nrerror(char error_text);void free_ivector(int *v, long nl, long nh);void GaussInverse(double *a, int n)int *indxc, *indxr, *ipiv;int i, icol, irow, j, k, l, ll;double big, dum, pi

12、vinv, temp;indxc = ivector(1,n);indxr = ivector(1,n);ipiv = ivector(1,n);for ( j = 1; j <= n; j+ ) ipivj = 0;for ( i = 1; i <= n; i+ ) big = 0.0;for (j = 1; j <= n; j+)if (ipivj != 1)for (k = 1;k <= n; k+) if (ipivk = 0) if (fabs(aj-1k-1) >= big)big = fabs(aj-1k-1);irow = j;icol = k;

13、else if ( ipivk > 1) nrerror("gaussj: Singular Matrix-1");+(ipivicol);if ( irow != icol ) for ( l = 1; l <= n; l+ )SWAP(airow-1l-1,aicol-1l-1); indxri = irow;indxci = icol;if (aicol-1icol-1 = 0.0) nrerror("gaussj: Singular Matrix-2");pivinv = 1.0/aicol-1icol-1;aicol-1icol-1

14、 = 1.0;for (l = 1; l <= n; l+) aicol-1l-1 *= pivinv; for (ll = 1; ll <= n; ll+)if ( ll != icol )dum = all-1icol-1;all-1icol-1 = 0.0;for ( l = 1; l <= n; l+) all-1l-1 -= aicol-1l-1 * dum; for (l = n; l >= 1; l-) if (indxrl != indxcl)for (k = 1; k <= n; k+)SWAP(ak-1indxrl-1, ak-1indxcl-

15、1);free_ivector(ipiv,1,n);free_ivector(indxr,1,n);free_ivector(indxc,1,n);#undef SWAP#undef NRANSI/生成高斯白噪声void CMathMatrix:random(int n, double *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

16、(-2.0 * log(tt / 32767.0);b = 2 * 3.14159265 * rand() / 32767.0;ej = var * a * cos(b) + mean;/矩阵与矢量相乘void CMathMatrix:Mat_Vector(double *Arr, int Rows, int Cols, double B, double C)int i, j;double Sum = 0.0;for (i = 0; i < Rows; i+)for (j = 0; j < Cols; j+)Sum += Arrij * Bj;Ci = Sum;Sum = 0.0;

17、(3)、CMathMatrix文件的头文件CMathMatrix.h:/ MathMatrix.h: interface for the CMathMatrix class./#if !defined(AFX_MATHMATRIX_H_718FA3C9_4408_4B30_BA72_F9C28742A8F5_INCLUDED_)#define AFX_MATHMATRIX_H_718FA3C9_4408_4B30_BA72_F9C28742A8F5_INCLUDED_#if _MSC_VER > 1000#pragma once#endif / _MSC_VER > 1000class CMathMatrix public:/生成高斯白噪声void random(int n, double *e, float mean, float var);/矩阵与矢量相乘void Mat_Vector(doubl