/** ar_model.c** Created on: 2013-8-11* Author: monkeyzx*/#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdlib.h>
//#include "msp.h"
#include "ar_model.h"
#include "time.h"float mabs(complex a)
{float m;m=a.real*a.real+a.imag*a.imag;m=sqrt(m);return(m);
}/*---------------------------------------------------------------------Routine MCORRE1:To estimate the biased cross-correlation functionof complex arrays x and y. If y=x,then it is auto-correlation.input parameters:x :n dimensioned complex array.y :n dimensioned complex array.n :the dimension of x and y.lag:point numbers of correlation.output parameters:r :lag dimensioned complex array, the correlation function isstored in r(0) to r(lag-1).in Chapter 1 and 11
---------------------------------------------------------------------*/
void mcorre1(complex x[],complex y[],complex r[],int n,int lag)
{int m,j,k;for(k=0;k<lag;k++) {m=n-1-k;r[k].real=0.0f;r[k].imag=0.0f;for(j=0;j<=m;j++) {r[k].real+=y[j+k].real*x[j].real+y[j+k].imag*x[j].imag;r[k].imag+=y[j+k].imag*x[j].real-y[j+k].real*x[j].imag;}r[k].real=r[k].real/n;r[k].imag=r[k].imag/n;}return;
}/*---------------------------------------------------------------------Routine maryuwa: To determine the autoregressive coefficients bysolving Yule-Walker equation with Levinson algorithm.Input Parameters:n : Number of data samples (integer)ip : Order of autoregressive modelx : Array of complex data values, x(0) to x(n-1)Output Parameters:ep : Driving noise variance (real)a : Array of complex autoregressive coefficients, a(0) toa(ip)ierror=0 : No error=1 : ep<=0 .r : complex work array, auto-correlationin chapter 12
--------------------------------------------------------------------*/
void maryuwa(complex x[],complex a[],complex r[],int n,int ip,
float *ep,int *ierror)
{complex sum;int i,k;float r0;*ierror=1;mcorre1(x,x,r,n,ip+1);a[0].real=1.0;a[0].imag=0.0;r0=r[0].real;a[1].real=-r[1].real/r0;a[1].imag=-r[1].imag/r0;*ep=r0*(1.0f-pow(mabs(a[1]),2));for(k=2;k<=ip;k++) {sum.real=0.;sum.imag=0.;for(i=1;i<k;i++) {sum.real+=r[k-i].real*a[i].real-r[k-i].imag*a[i].imag;sum.imag+=r[k-i].real*a[i].imag+r[k-i].imag*a[i].real;}sum.real+=r[k].real;sum.imag+=r[k].imag;a[k].real=-sum.real/(*ep);a[k].imag=-sum.imag/(*ep);(*ep)*=1.-pow(mabs(a[k]),2);if(*ep<=0.0)return;for(i=1;i<k;i++) {x[i].real=a[i].real+a[k-i].real*a[k].real+a[k-i].imag*a[k].imag;x[i].imag=a[i].imag+a[k-i].real*a[k].imag-a[k-i].imag*a[k].real;}for(i=1;i<k;i++) {a[i].real=x[i].real;a[i].imag=x[i].imag;}}*ierror=0;
}/*----------------------------------------------------------------------routinue mrelfft:To perform split-radix DIF fft algorithm.input parameters:xr,xi:real and image part of complex data for DFT/IDFT,n=0,...,N-1N :Data point number of DFT compute .isign:Transform direction disignator ,isign=-1: For Forward Transform.isign=+1: For Inverse Transform.output parameters:xr,xi:real and image part of complex result of DFT/IDFT,n=0,...,N-1Note: N must be a power of 2 .in chapter 5
---------------------------------------------------------------------*/
void mrelfft(float xr[],float xi[],int n,int isign)
{float e,es,cc1,ss1,cc3,ss3,r1,s1,r2,s2,s3,xtr,xti,a,a3;int m,n2,n4,j,k,is,id,i0,i1,i2,i3,n1,i,nn;for(m=1;m<=16;m++) {nn=pow(2,m);if(n==nn)break;}if(m>16) {
#ifdef _DEBUGprintf(" N is not a power of 2 ! \n");
#endifreturn;}n2=n*2;es=-isign*atan(1.0)*8.0;for(k=1;k<m;k++) {n2=n2/2;n4=n2/4;e=es/n2;a=0.0;for(j=0;j<n4;j++) {a3=3*a;cc1=cos(a);ss1=sin(a);cc3=cos(a3);ss3=sin(a3);a=(j+1)*e;is=j;id=2*n2;do {for(i0=is;i0<n;i0+=id) {i1=i0+n4;i2=i1+n4;i3=i2+n4;r1=xr[i0]-xr[i2];s1=xi[i0]-xi[i2];r2=xr[i1]-xr[i3];s2=xi[i1]-xi[i3];xr[i0]+=xr[i2];xi[i0]+=xi[i2];xr[i1]+=xr[i3];xi[i1]+=xi[i3];if(isign!=1) {s3=r1-s2;r1=r1+s2;s2=r2-s1;r2=r2+s1;} else {s3=r1+s2;r1=r1-s2;s2=-r2-s1;r2=-r2+s1;}xr[i2]=r1*cc1-s2*ss1;xi[i2]=-s2*cc1-r1*ss1;xr[i3]=s3*cc3+r2*ss3;xi[i3]=r2*cc3-s3*ss3;}is=2*id-n2+j;id=4*id;}while(is<n-1);}}
/* ------------ special last stage -------------------------*/is=0;id=4;do {for(i0=is;i0<n;i0+=id) {i1=i0+1;xtr=xr[i0];xti=xi[i0];xr[i0]=xtr+xr[i1];xi[i0]=xti+xi[i1];xr[i1]=xtr-xr[i1];xi[i1]=xti-xi[i1];}is=2*id-2;id=4*id;} while(is<n-1);j=1;n1=n-1;for(i=1;i<=n1;i++) {if(i<j) {xtr=xr[j-1];xti=xi[j-1];xr[j-1]=xr[i-1];xi[j-1]=xi[i-1];xr[i-1]=xtr;xi[i-1]=xti;}k=n/2;while(1) {if(k>=j)break;j=j-k;k=k/2;}j=j+k;}if(isign==-1) return;for(i=0;i<n;i++) {xr[i]/=n;xi[i]/=n;}
}/*---------------------------------------------------------------------Routine mpsplot: To plot the normalized power spectum curve on thenormalized frequency axis from -.5 to +.5 .mfre : Points in frequency axis and must be the power of 2.ts : Sample interval in seconds (real).psdr : Real array of power spectral density values.psdi : Real work array.in chapter 11,12
--------------------------------------------------------------------*/
void mpsplot(float psdr[],float psdi[],int mfre,float ts)
{FILE *fp;char filename[30];int k,m2;float pmax,fs,faxis;m2=mfre/2;for(k=0;k<m2;k++){psdi[k]=psdr[k];psdr[k]=psdr[k+m2];psdr[k+m2]=psdi[k];}pmax=psdr[0];for(k=1;k<mfre;k++)if(psdr[k]>pmax)pmax=psdr[k];for(k=0;k<mfre;k++) {psdr[k]=psdr[k]/pmax;if(psdr[k]<=0.0)psdr[k]=.000001;}fs=1./ts;fs=fs/(float)(mfre);printf("Please input filename:\n");scanf("%s",filename);if((fp=fopen(filename,"w"))==NULL) {printf("cannot open file\n");exit(0);}for(k=0;k<mfre;k++) {faxis=fs*(k-m2);fprintf(fp,"%f,%f\n",faxis,10.*log10(psdr[k]));}fclose(fp);return;
}/*----------------------------------------------------------------------Routine mar1psd: To compute the power spectum by AR-model parameters.Input parameters:ip : AR model order (integer)ep : White noise variance of model input (real)ts : Sample interval in seconds (real)a : Complex array of AR parameters a(0) to a(ip)Output parameters:psdr : Real array of power spectral density valuespsdi : Real work arrayin chapter 12
---------------------------------------------------------------------*/
void mar1psd(complex a[],int ip,int mfre,float *ep,float ts)
{static float psdr[4096];static float psdi[4096];int k;float p;for(k=0;k<=ip;k++) {psdr[k]=a[k].real;psdi[k]=a[k].imag;}for(k=ip+1;k<mfre;k++) {psdr[k]=0.;psdi[k]=0.;}mrelfft(psdr,psdi,mfre,-1);for(k=0;k<mfre;k++) {p=pow(psdr[k],2)+pow(psdi[k],2);psdr[k]=(*ep)*ts/p;}mpsplot(psdr,psdi,mfre,ts);return;
}/** Below are examples for using @maryuwa and @mar1psd*/
#define PI (3.1415926)
#define N (1024)
#define AN (10)
complex x[N];
complex r[N];
complex a[AN];/** generate random number which satify guass distribution*/
double guass_rand(void)
{static double V1, V2, S;static int phase = 0;double X;if ( phase == 0 ) {do {double U1 = (double)rand() / RAND_MAX;double U2 = (double)rand() / RAND_MAX;V1 = 2 * U1 - 1;V2 = 2 * U2 - 1;S = V1 * V1 + V2 * V2;} while(S >= 1 || S == 0);X = V1 * sqrt(-2 * log(S) / S);} else {X = V2 * sqrt(-2 * log(S) / S);}phase = 1 - phase;return X;
}void zx_ar_model(void)
{int i=0;float ep = 0;int ierror = 0;/** generate x[N]*/srand(time(NULL));for (i=0; i<N; i++) {x[i].real = sin(2*PI*i/N) + guass_rand();x[i].imag = 0;}/* Find parameters for AR model */maryuwa(x, a, r, N, AN, &ep, &ierror);/* Calculate power spectum using parameters of AR model */mar1psd(a, AN, N, &ep, 1);
}
[cpp] view plaincopy print?
/*?
?*?main.c?
?*?
?*??Created?on:?2013-8-11?
?*??????Author:?monkeyzx?
?*/??
#include?"ar_model.h"??
??
int?main(void)??
{??
????zx_ar_model();??
??
????return?0;??
}??
/** main.c** Created on: 2013-8-11* Author: monkeyzx*/
#include "ar_model.h"int main(void)
{zx_ar_model();return 0;
}
上面的實例中給定輸入信號為余弦信號，采樣點數為1024個點，通過計算后的功率譜通過mpsplot函數保存到文本文件output.txt中，保存格式如下：
