本帖最后由 look_w 于 2017-9-23 15:35 编辑
- /*
- * 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 function
- of 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 is
- stored in r(0) to r(lag-1).
- in Chapter 1 and 11
- 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.real-r[k-i].imag*a.imag;
- sum.imag+=r[k-i].real*a.imag+r[k-i].imag*a.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.real=a.real+a[k-i].real*a[k].real+
- a[k-i].imag*a[k].imag;
- x.imag=a.imag+a[k-i].real*a[k].imag-
- a[k-i].imag*a[k].real;
- }
- for(i=1;i<k;i++) {
- a.real=x.real;
- a.imag=x.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-1
- N
 ata 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-1
- Note: 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 _DEBUG
- printf(" N is not a power of 2 ! \n");
- #endif
- return;
- }
- 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/=n;
- xi/=n;
- }
- }
- /*---------------------------------------------------------------------
|
