6.cpp

用传统的方法进行双谱估计与用传统的方法进行功率谱估计类似

 #include <math.h>
  #define M_PI 3.14159265358979323846
  typedef struct{
   float r,i;
   }complex;
   
 typedef struct{
  double r,i;
  
  }doublecomplex;
  void complex_exp(complex *r,complex *z)
  {double expx;
   expx=exp((double)z->r);
   r->r=(float)expx*cos((double)z->i);
   r->i=(float)expx*sin((double)z->i);
   
   }
 void fftc(complex *x,long *n,long *isign)  
 {/* Local variables*/
  long i,l,m,mr,tmp_int;
  complex t,tmp_complex,tmp;
  float pisign;
  pisign= (float)((double )*isign *M_PI);
  mr=0;
  for (m=1;m<*n;++m)
   {l=*n;
   l/=2;
   while (mr+1>=*n)
    {l/=2;
    }mr=mr%l+l;
    if(mr>m)
      {t.r=x[m].r;
       t.i=x[m].i;
       x[m].r=x[mr].r; 
       x[m].i=x[mr].i;
       x[mr].r=t.r;
       x[mr].i=t.i;
       
      }
   }
   l=1;
   while(l<*n)
    {for(m=0;m<l;++m)
     {tmp_int=l*2;
      for(i=m;tmp_int<0?i>=(*n-1):i<*n;i+=tmp_int)
       {tmp.r=0.0;
        tmp.i=(float)m*pisign/(float)l;
        complex_exp(&tmp_complex,&tmp);
        t.r=x[i+l].r*tmp_complex.r-x[i+l].i*tmp_complex.i;
        t.i=x[i+l].r*tmp_complex.i+x[i+l].i*tmp_complex.r;
        x[i+l].r=x[i].r-t.r;
        x[i+l].i=x[i].i-t.i;
        x[i].r=x[i].r+t.r;
        x[i].i=x[i].i+t.i;
       
        }
      }l*=2;
    }return;
  
  }