利用协方差法估计AR模型参数.docx

随机信号分析基础大作业

利用协方差法估计AR模型参数
进而估计功率谱

严 奎（学号：3222008008）
陈 韬（学号：3222008022）
朱燕豪（学号：3222008021）

2011年01月15日

[键入文字]

作业综述：
本作业中采用面向对象的程序设计方法，将用到的子程序封装在一个类中，防止其他函数的干扰，具有良好的信息内聚性。类中定义的有获得（0,1）分布随机数的函数uniform()，产生高斯分布随机数的函数gauss()，产生自回归滑动平均模型ARMA（p,q）数据的函数arma()，用乔布斯基（Cholesky）算法求解对称正定方程组的函数cholesky（），计算ARMA模型的功率谱密度的函数psd（），用协方差方法估计AR模型参数，进而实现功率谱估计的函数covar（）。采用的编程工具是VC++6.0以及VS2010，用MATLAB对生成的数据进行画图。
一．题目要求
给定一段信号数据及采样率，利用现代谱估计理论编程估计信号的功率谱。

二．基本原理及方法
现代谱估计是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱，主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的，应用最广的是AR参数模型。现代谱估计的参数模型有自回归滑动平均（ARMA）模型、自回归（AR）模型、滑动平均（MA）模型，Wold分解定理阐明了三者之间的关系：任何有限方差的ARMA或MA模型的平稳随机过程可以用无限阶的AR模型表示，任何有限方差的ARMA或MA模型的平稳随机过程可以用无限阶的AR模型表示。但是由于只有AR模型参数估计是一组线性方程，而实际的物理系统往往是全极点系统，因而AR应用最广。
我们用协方差法估计AR模型参数，进而实现功率谱估计。若已知平稳随机序列x（n）的AR模型为

其中a（i）是AR系数，w(n)是均值为零，均方差为σ的白噪声。
1. 计算协方差

2. 用乔布斯基算法解对称正定方程组
N阶对称正定方程组的矩阵形式为AX=B，即

矩阵A的乔布斯基分解

这里D是主对角元素都为正实数的对角阵，即D=diag(d1,d2,…,dn)，L为主对角元素是1的下三角矩阵。用乔布斯基算法解对称正定方程组的方法是，先用回代法求解方程组LY=B，得到Y之后，再用回代法求解方程

3.计算激励白噪声的方差

4．用AR模型参数的估计值，可以计算功率谱密度

三．算法设计与实现
1．程序流程图
采用协方差的方法进行功率谱估计。如下图所示
开始

输入有限序列AR模型系数

根据噪声均值、方差产生自高斯白噪声

产生自回归滑动平均模型ARMA（p,q）模型的数据

用协方差法估计无限序列AR模型参数

计算AR模型系数功率谱密度

根据已存储的数据用Matlab做图

结束

图1算法流程图1
2.主要模块的设计：
1. 产生随机序列的函数uniform（），
采用线性同余法由种子seed产生随机数。
2. 产生高斯白噪声的函数gauss（），
gauss(double mean,double sigma,long int * s)
{
int i;double x,y;
for(x=0,i=0;i<12;i++)
x+=uniform(0.0,1.0,s);
x=x-6.0;
y=mean+x*sigma;
return(y);
}
3. ARMA模型数据的生成函数为arma（）略
4. 乔里斯基算法解对称正定方程组的函数cholesky（）略
5. 由协方差函数covar（）求AR参数；
6. 再根据AR参数求出功率谱的函数psd()略；
7. 最终用MATLAB的 画图工具给出直观的功率谱图形，

四．结果分析
输入平稳随机序列x（n）的AR模型为

其中1，-2.76,3.809,-2.654,0.924为AR系数，
根据要求产生W(n)是均值为零，方差为1的白噪声。
根据均匀分布产生（0,1）分布的随机序列，再由均值和方差生成高斯白噪声如下图所示：

由图可知产生的随机序列近似于高斯分布，符合题目要求。
由白噪声求自回归滑动平均模型ARMA（p,q）模型的数据，
用协方差法估计AR模型参数，结果为：
a(0)= 1.0000000
a(1)=—2.7310949
a(2)= 3.7478402
a(3)=—2.5951549
a(4)= 0.9022404
可以看出估计出的AR模型参数与原AR模型系数基本接近，但是不相等，这是因为现代谱估计是由有限长序列估计无限长的随机序列AR模型参数，但是结果基本接近。
其中预测误差功率是Pe=1.0995336，与原方差1较接近。
计算AR模型系数功率谱密度
根据已存储在covar1.dat的数据，用Matlab做图
在归一化频率的基础上做的功率谱

五．任务分工
三人合作进行了前期的资料查找，阅读文献，确定现代谱估计，分析算法。
严 奎 （学号：3222008008）完成了程序调试，绘图。
陈 韬 （学号：3222008022）完成答辩PPT的制作，以及负责主讲。
朱燕豪 （学号：3222008021）完成论文的撰写。
n 六.心得
通过这次的大作业提高了我们的合作能力，文献查取能力，编程能力，使我们掌握了书本上的知识，复习了前面的高斯分布，白噪声的产生，特别是掌握了功率谱的多种分析方法，了解了现代谱估计的方法与原理，极大地提高了我们的综合能力。
在选题时，我们以勇于专研问题的精神，选了现代谱分析。在做课题时，我们发现了很多问题，自己对谱分析的了解只停留在很基础的方面。特别是在完成算法分析时我们花了很多时间，开始我们只建立了AR模型，为了更加完善，我们加上了ARMA模型，最后在此基础上我们采用协方差分析使结果更趋于逼真。程序编写时，我们参照了大量的网上资源，但是调试过程中，变量的定义出了很多问题，很多地方都出了问题，我们只能一步一步调试改进。
虽然开始时我们遇到很多困难，编程能力太差，书本知识体系不完整缺少功率谱分析具体算法，上网条件差，图书馆资源有限等。但是怀着认真、踏实的态度我们完成了预期的任务，达到了一定的效果。总的来说，这次的课题我们都收获颇多。
七．参考文献
[1]殷福亮，宋爱军数字信号C语言程序集.辽宁科技出版社，1997
[2]张贤达，现代信号处理，清华大学出版社，2002
[3]常建军，李海林，随机信号分析，科学出版社，2006
八．附录
程序源代码
#include<iostream>
#include<cstddef>
#include"stdlib.h"
#include"stdio.h"
#include"math.h"
#include"malloc.h"

using namespace std;
class Power
{
public:
Power(){}
~Power(){}
double uniform(double a,double b,long int * seed);
double gauss(double mean,double sigma,long int * s);
void cholesky_1(double a[],double b[],int n);
void covar(double x[],int n,int p,double a[],double *v,int mode);
void arma(double a[],double b[],int p,int q,double mean,double sigma,long *seed,double x[],int n);
void psd(double b[],double a[],int q,int p,double sigma2,double fs,double x[],double freq[],int len,int sign);
};
double Power::uniform(double a,double b,long int *seed)
{
double t;
*seed=2045*(*seed)+1; //seed为种子
*seed=*seed-(*seed/1048576)*1048576;
t=(*seed)/1048576.0;
t=a+(b-a)*t;
return(t);
}
double Power::gauss(double mean,double sigma,long int * s)
{
int i;double x,y;
for(x=0,i=0;i<12;i++)
x+=uniform(0.0,1.0,s);
x=x-6.0;
y=mean+x*sigma;
return(y);
}

void Power::cholesky_1(double a[],double b[],int n)
{
int i,j,k,m;
double *d,*y,*xl,eps;
d=(double *)malloc(n*sizeof(double));
y=(double *)malloc(n*sizeof(double));
xl=(double *)malloc(n*n*sizeof(double));
eps=1.0e-15;
m=0;
d[0]=a[m];
for(i=1;i<n;i++)
{
for(j=0;j<i;j++)
{
m=m+1;
xl[i*n+j]=a[m]/d[j];
if(j==0)continue;
for(k=0;k<j;k++)
{xl[i*n+j]=xl[i*n+k]*xl[j*n+k]*d[k]/d[j];}
}
m=m+1;
d[i]=a[m];
for(k=0;k<i;k++)
{d[i]=d[i]-d[k]*xl[i*n+k]*xl[i*n+k];}
if(fabs(d[i])<eps)
{
printf("\nill-conditioned! \n");
return;
}
}
y[0]=b[0];
for(k=1;k<n;k++)
{
y[k]=b[k];
for(j=0;j<k;j++)
{y[k]=y[k]-xl[k*n+j]*y[j];}
}
b[n-1]=y[n-1]/d[n-1];
for(k=(n-2);k>=0;k--)
{
b[k]=y[k]/d[k];
for(j=(k+1);j<n;j++)
{b[k]=b[k]-xl[j*n+k]*b[j];}
}
free(d);
free(y);
free(xl);
}
void Power::covar(double x[],int n,int p,double a[],double *v,int mode)
{
int i,j,k,m;
double cc,sum,*c;
c=(double *)malloc((p*(p+1)/2)*sizeof(double));
m=0;
for(k=1;k<=p;k++)
{
for(j=1;j<=k;j++)
{
c[m]=0.0;
for(i=p;i<n;i++)
{
c[m]+=x[i-j]*x[i-k];
}
if(mode==1)
{
for(i=0;i<(n-p);i++)
{
c[m]+=x[i+j]*x[i+k];// 计算Cxx(i,k)
}
}
m=m+1;
}
}
for(j=1;j<=p;j++)
{
a[j-1]=0.0;
for(i=p;i<n;i++)
{
a[j-1]-=x[i-j]*x[i];
}
if(mode==1)
{
for(i=0;i<(n-p);i++)
{
a[j-1]-=x[i+j]*x[i]; //计算Cxx(j,0)
}
}
}
cholesky_1(c,a,p); //解得a(i)
for(k=(p-1);k>=0;k--)
{
a[k+1]=a[k];
}
a[0]=1.0;
sum=0.0;
for(k=0;k<=p;k++)
{
cc=0.0;
for(i=p;i<n;i++)
{
cc+=x[i]*x[i-k];
}
if(mode==1)
{
for(i=0;i<(n-p);i++)
{
cc+=x[i]*x[i+k]; //计算Cxx(0,k)
}
}
if(k==0)
{
sum+=cc;
}
else
{
sum+=cc*a[k]; //计算a(k)*Cxx(0,k)

}
}
if(mode==1)
{
v[0]=sum/(2*(n-p));
}
else
{
v[0]=sum/(n-p); //计算sigma2
}
free(c);
}
void Power::arma(double a[],double b[],int p,int q,double mean,double sigma,long *seed,double x[],int n)
{
int i,k,m;
double s,*w;
w=(double *)malloc(n*sizeof(double));
for(k=0;k<n;k++)
w[k]=gauss(mean,sigma,seed); //得到高斯白噪声
x[0]=b[0]*w[0];
for(k=1;k<=p;k++) //得到前p个数据
{
s=0.0;
for(i=1;i<=k;i++)
{
s+=a[i]*x[k-i];
}
s=b[0]*w[k]-s;
if(q==0)
{
x[k]=s;
continue;
}
m=(k>q)?q:k;
for(i=1;i<=m;i++)
{
s+=b[i]*w[k-i];
}
x[k]=s;
}
for(k=(p+1);k<n;k++) //得到p+1到n的数据
{
s=0.0;
for(i=1;i<=p;i++)
{
s+=a[i]*x[k-i];
}
s=b[0]*w[k]-s;
if(q==0)
{
x[k]=s;
continue;
}
for(i=1;i<=q;i++)
{s+=b[i]*w[k-i];}
x[k]=s;
}
free(w);
}
void Power::psd(double b[],double a[],int q,int p,double sigma2,double fs,double x[],double freq[],int len,int sign)
{
int i,k;
double ar,ai,br,bi,zr,zi,im,re,xre,xim;
double ang,den,numr,numi,temp;
for(k=0;k<len;k++)
{
ang=k*0.5/(len-1);
freq[k]=ang*fs;
zr=cos(-8.0*atan(1.0)*ang);
zi=sin(-8.0*atan(1.0)*ang);
br=0.0;
bi=0.0;
for(i=q;i>0;i--)
{
re=br;
im=bi;
br=(re+b[i])*zr-im*zi;
bi=(re+b[i])*zi+im*zr; //分子的傅里叶变换
}
ar=0.0;
ai=0.0;
for(i=p;i>0;i--)
{
re=ar;
im=ai;
ar=(re+a[i])*zr-im*zi; //分母的傅里叶变换
ai=(re+a[i])*zi+im*zr;
}
br=br+b[0];
ar=ar+1.0;
numr=ar*br+ai*bi; //分母有理化后分子的实部
numi=ar*bi-ai*br;
den=ar*ar+ai*ai;
xre=numr/den;
xim=numi/den;
switch(sign)
{
case 0:
{
x[k]=xre*xre+xim*xim;
x[k]=sigma2*x[k]/fs;
break;

}
case 1:
{
temp=xre*xre+xim*xim;
temp=sigma2*temp/fs;
if(temp==0.0)
temp=1.0e-20;
x[k]=10.0*log10(temp);
}
}
}
}
void main()
{
Power P;
int i,n,p,q,len;
long seed;
double v,mean,var,c[10],x[500],freq[200];
double fs,sigma2;
static double a[5]={1.0,-2.76,3.809,-2.645,0.924};
static double b[1]={1.0};
FILE *fp;
p=4;
q=0;
seed=135791;
mean=0.0;
var=1.0;
n=500;
P.arma(a,b,p,q,mean,var,&seed,x,n);
for(i=0;i<300;i++)
x[i]=x[i+200];
n=300;
P.covar(x,n,p,c,&v,0);
printf("The coefficient of AR model\n");
for(i=0;i<=p;i++)
{
printf("a(%d)=%10.7lf\n",i,c[i]);
}
printf("The reflet coefficient of AR model\n");
printf("Pe=%10.7lf\n",v);
fs=1.0;
sigma2=v;
len=200;
P.psd(b,c,q,p,sigma2,fs,x,freq,len,1);
fp=fopen("covar1.dat","w");
for(i=0;i<len;i++)
fprintf(fp,"%lf %lf\n",freq[i],x[i]);
fclose(fp);
}