利用协方差法估计AR模型参数(2)

{

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) {

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

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]; // } } }

cholesky_1(c,a,p); // 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(j,0) 解得a(i) 计算Cxx(0,k)

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

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++)

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

{

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模型参数进而估计功率谱

}

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];

分母有理化后

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

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); }