利用协方差法估计AR模型参数Word格式文档下载.docx
- 文档编号:18267626
- 上传时间:2022-12-14
- 格式:DOCX
- 页数:17
- 大小:94.88KB
利用协方差法估计AR模型参数Word格式文档下载.docx
《利用协方差法估计AR模型参数Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《利用协方差法估计AR模型参数Word格式文档下载.docx(17页珍藏版)》请在冰豆网上搜索。
如下图所示
图1算法流程图1
2.主要模块的设计:
1.产生随机序列的函数uniform(),
采用线性同余法由种子seed产生随机数。
2.产生高斯白噪声的函数gauss(),
gauss(doublemean,doublesigma,longint*s)
{
inti;
doublex,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)完成论文的撰写。
⏹六.心得
通过这次的大作业提高了我们的合作能力,文献查取能力,编程能力,使我们掌握了书本上的知识,复习了前面的高斯分布,白噪声的产生,特别是掌握了功率谱的多种分析方法,了解了现代谱估计的方法与原理,极大地提高了我们的综合能力。
在选题时,我们以勇于专研问题的精神,选了现代谱分析。
在做课题时,我们发现了很多问题,自己对谱分析的了解只停留在很基础的方面。
特别是在完成算法分析时我们花了很多时间,开始我们只建立了AR模型,为了更加完善,我们加上了ARMA模型,最后在此基础上我们采用协方差分析使结果更趋于逼真。
程序编写时,我们参照了大量的网上资源,但是调试过程中,变量的定义出了很多问题,很多地方都出了问题,我们只能一步一步调试改进。
虽然开始时我们遇到很多困难,编程能力太差,书本知识体系不完整缺少功率谱分析具体算法,上网条件差,图书馆资源有限等。
但是怀着认真、踏实的态度我们完成了预期的任务,达到了一定的效果。
总的来说,这次的课题我们都收获颇多。
七.参考文献
[1]殷福亮,宋爱军数字信号C语言程序集.辽宁科技出版社,1997
[2]张贤达,现代信号处理,清华大学出版社,2002
[3]常建军,李海林,随机信号分析,科学出版社,2006
八.附录
程序源代码
#include<
iostream>
cstddef>
#include"
stdlib.h"
stdio.h"
math.h"
malloc.h"
usingnamespacestd;
classPower
public:
Power(){}
~Power(){}
doubleuniform(doublea,doubleb,longint*seed);
doublegauss(doublemean,doublesigma,longint*s);
voidcholesky_1(doublea[],doubleb[],intn);
voidcovar(doublex[],intn,intp,doublea[],double*v,intmode);
voidarma(doublea[],doubleb[],intp,intq,doublemean,doublesigma,long*seed,doublex[],intn);
voidpsd(doubleb[],doublea[],intq,intp,doublesigma2,doublefs,doublex[],doublefreq[],intlen,intsign);
};
doublePower:
:
uniform(doublea,doubleb,longint*seed)
doublet;
*seed=2045*(*seed)+1;
//seed为种子
*seed=*seed-(*seed/1048576)*1048576;
t=(*seed)/1048576.0;
t=a+(b-a)*t;
return(t);
voidPower:
cholesky_1(doublea[],doubleb[],intn)
inti,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;
n;
{
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;
{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;
y[k]=b[k];
k;
{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);
{b[k]=b[k]-xl[j*n+k]*b[j];
free(d);
free(y);
free(xl);
covar(doublex[],intn,intp,doublea[],double*v,intmode)
doublecc,sum,*c;
c=(double*)malloc((p*(p+1)/2)*sizeof(double));
=p;
for(j=1;
=k;
c[m]=0.0;
for(i=p;
{
c[m]+=x[i-j]*x[i-k];
}
if(mode==1)
for(i=0;
(n-p);
{
c[m]+=x[i+j]*x[i+k];
//计算Cxx(i,k)
}
for(j=1;
a[j-1]=0.0;
for(i=p;
a[j-1]-=x[i-j]*x[i];
if(mode==1)
for(i=0;
a[j-1]-=x[i+j]*x[i];
//计算Cxx(j,0)
cholesky_1(c,a,p);
//解得a(i)
for(k=(p-1);
a[k+1]=a[k];
a[0]=1.0;
sum=0.0;
cc=0.0;
cc+=x[i]*x[i-k];
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);
arma(doublea[],doubleb[],intp,intq,doublemean,doublesigma,long*seed,doublex[],intn)
inti,k,m;
doubles,*w;
w=(double*)malloc(n*sizeof(double));
w[k]=gauss(mean,sigma,seed);
//得到高斯白噪声
x[0]=b[0]*w[0];
k++)//得到前p个数据
s=0.0;
for(i=1;
s+=a[i]*x[k-i];
s=b[0]*w[k]-s;
if(q==0)
x[k]=s;
continue;
m=(k>
q)?
q:
=m;
s+=b[i]*w[k-i];
x[k]=s;
}
for(k=(p+1);
k++)//得到p+1到n的数据
s=0.0;
for(i=1;
s+=a[i]*x[k-i];
=q;
{s+=b[i]*w[k-i];
}
free(w);
psd(doubleb[],doublea[],intq,intp,doublesigma2,doublefs,doublex[],doublefreq[],intlen,intsign)
inti,k;
doublear,ai,br,bi,zr,zi,im,re,xre,xim;
doubleang,den,numr,numi,temp;
for(k=0;
len;
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;
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)
case0:
{
x[k]=xre*xre+xim*xim;
x[k]=sigma2*x[k]/fs;
break;
}
case1:
temp=xre*xre+xim*xim;
temp=sigma2*temp/fs;
if(temp==0.0)
temp=1.0e-20;
x[k]=10.0*log10(temp);
voidmain()
PowerP;
inti,n,p,q,len;
longseed;
doublev,mean,var,c[10],x[500],freq[200];
doublefs,sigma2;
staticdoublea[5]={1.0,-2.76,3.809,-2.645,0.924};
staticdoubleb[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;
300;
x[i]=x[i+200];
n=300;
P.covar(x,n,p,c,&
v,0);
printf("
ThecoefficientofARmodel\n"
printf("
a(%d)=%10.7lf\n"
i,c[i]);
TherefletcoefficientofARmodel\n"
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"
fprintf(fp,"
%lf%lf\n"
freq[i],x[i]);
fclose(fp);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 利用 协方差 估计 AR 模型 参数