fft2.cpp

fft变换程序

  //small   fft   program,using   complex.   
    
    
  #include   "complex"   
  #include   <iostream>   
  #include   <conio.h>   
  using   namespace   std;   //very   important   
  #define   PI   3.1415926   
  typedef   complex<double>   complexd;   //use   template   
  void   myfft(double   y[][2],   int   N)   //N   is   the   power   of   2   (as   8,64,   1024)   
  {   
      complexd   t1,t2,t;   
      int   u,mu,i,j,k,n,l,s;   
      complexd   *w   =     new   complexd[N/2];   
      complexd   *x   =     new   complexd[N];   
      mu=0;   u=N-1;   
      while   (u>0)   {mu++;   u/=2;};   
    
      for   (i=0;   i<N;   i++)   x[i   =complexd(y[i][0],y[i][1]);   
      for   (i=0;   i<N;   i++)   {   //sort   input   number   array   
          k=i;   l=0;   
          for   (j=2;   j<=N;   j=j*2)   {   
              l+=(N/j)*(k%2);   k=k/2;   
          }   
          if   (l>i)   {   
              t=x[i];   
              x[i   =x[l];   
              x[l   =t;   
          }   
      }   
      for   (i=0;   i<N/2;   i++)   w[i   =complexd(cos(2*PI/N*i),sin(2*PI/N*i));   
      s=1;   
      for   (i=0;   i<mu;   i++)   {   //fft   by   binary   dividing   
          for   (j=0;   j<N;   j+=s*2)   {   
              for   (k=0;   k<s;   k++)   {   
      t=x[j+k+s]*(w[k*N/s/2]);   
      t1=x[j+k]+t;   
      t2=x[j+k]-t;   
      x[j+k]   =t1;   
      x[j+k+s]   =t2;   
              }   
          }   
          s*=2;   
      }   
      for   (i=0;   i<N;   i++)   {y[i][0]   =x[i].real();   y[i][1]   =x[i].imag();}   
      delete   w;   //very   important.   If   not,   the   system   will   lost   memory   
      delete   x;   
  }   
  void   main()   
  {   
    int   i;   
    double   u[8][2];   
    for   (i=0;   i<8;   i++)   u[i][0   =i+1;   //the   input   
    myfft(u,8);   
    for   (i=0;   i<8;   i++)   cout<<'('<<u[i][0]<<','<<u[i][1]<<')'<<endl;   
    _getch();   
  }