利用协方差法估计AR模型参数Word格式.docx
《利用协方差法估计AR模型参数Word格式.docx》由会员分享，可在线阅读，更多相关《利用协方差法估计AR模型参数Word格式.docx（20页珍藏版）》请在冰豆网上搜索。
如下图所示
图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）;
升级会员 