capon.cpp

信号处理中的谱分析

#include "iostream.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include <fstream.h>
#include <iomanip.h>
void burg(double x[],int n,int p,double a[],double *v);
void mlm(double x[],int n,int pm,double fs,double s[],double freq[],int len,int sdb);
void main()
{
	double freq[200],b[1];
    double x[200],fs,sigma2;
	double data_frist[228],a[14],v,data_array[1368];
	int data_size=1368,i,j,p,n,len;
    ifstream data_file("test.txt");
	ofstream data_file1("test1.txt");
	 for(i=0;i<data_size;i++)
	 {
		 data_file>>data_array[i];
	 }
	 for(i=0;i<data_size/6;i++)
	 {
		 j=6*i;
		 data_frist[i]=data_array[j];
	 }
	 b[0]=1.0;
	n=data_size/6;
	p=13;	
	burg(data_frist,n,p,a,&v);
	for(i=0;i<=p;i++)
	{
		cout<<a[i]<<endl;
	}
	len=200;
	sigma2=v;
	fs=1.0;
	mlm(data_frist,n,p,fs,x,freq,len,1);
	for(i=0;i<len;i++)
	{
		data_file1<<setprecision(13)<<freq[i]<<"    "<<x[i]<<endl;
	}

}
void burg(double x[],int n,int p,double a[],double *v)
{ 
	int i,k;
	double r,sumd,sumn,*b,*ef,*eb;
	b=(double *)malloc((p+1)*sizeof(double));
	ef=(double *)malloc(n*sizeof(double));
	eb=(double *)malloc(n*sizeof(double));
	a[0]=1.0;
	r=0.0;
	for(i=0;i<n;i++)
	{
		r+=x[i]*x[i]/n;
	}
	v[0]=r;          

	for(i=1;i<n;i++)
	{
		ef[i]=x[i];
		eb[i-1]=x[i-1];
	}
	
	for(k=1;k<=p;k++)
	{
		sumn=0.0;
		sumd=0.0;
		for(i=k;i<n;i++)
		{
			sumn+=ef[i]*eb[i-1];
			sumd+=ef[i]*ef[i]+eb[i-1]*eb[i-1];
		}
		a[k]=-2*sumn/sumd;
		for(i=1;i<k;i++)
		{
			b[i]=a[i]+a[k]*a[k-i];
		}
		for(i=1;i<k;i++)
		{
			a[i]=b[i];
		}
		v[0]=(1.0-a[k]*a[k])*v[0];
		for(i=(n-1);i>=(k+1);i--)
		{
			ef[i]=ef[i]+a[k]*eb[i-1];
			eb[i-1]=eb[i-2]+a[k]*ef[i-1];
		}
	}
	free(b);
	free(ef);
	free(eb);
}

void mlm(double x[],int n,int pm,double fs,double s[],double freq[],int len,int sdb)
{
	int i,k,p;
	double v,ai,ar,im,re,zi,zr,sar,ang,tpi,*a;
	a=(double*)malloc((pm+1)*sizeof(double));
	tpi=8.0*atan(1.0);
	for(i=0;i<len;i++)
	{
		s[i]=0.0;
	}
	for(p=0;p<=pm;p++)
	{
		burg(x,n,p,a,&v);
		for(k=0;k<len;k++)
		{
			ang=k*1.0/(len-1);
			freq[k]=ang*fs;
			zr=cos(-tpi*ang);
			zi=sin(-tpi*ang);
			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;
			}
			ar=ar+1.0;
			sar=v/(fs*(ar*ar+ai*ai));
			s[k]=s[k]+1.0/sar;
		}
	}
	for(k=0;k<len;k++)
	{
		s[k]=pm/s[k];
	}
	if(sdb==1)
	{
		for(k=0;k<len;k++)
		{
			s[k]=10.0*log10(s[k]);
		}
	}
	free(a);
}