fft1.cpp
上传用户:alisonmail
上传日期:2013-02-28
资源大小:500k
文件大小:1k
源码类别:

图片显示

开发平台:

Visual C++

  1. #include "stdafx.h"
  2. #include "complex"
  3. using namespace std; //very important
  4. #define PI 3.1415926
  5. typedef complex<double> complexd; //use template 
  6. void myfft(double y[][2], int N) //N is the power of 2 (as 8,64, 1024)
  7. {
  8.   complexd t1,t2,t;
  9.   int u,mu,i,j,k,n,l,s;
  10.   complexd *w =  new complexd[N/2];
  11.   complexd *x =  new complexd[N];
  12.   mu=0; u=N-1;
  13.   while (u>0) {mu++; u/=2;};
  14.   for (i=0; i<N; i++) x[i]=complexd(y[i][0],y[i][1]);
  15.   for (i=0; i<N; i++) { //sort input number array
  16.     k=i; l=0;
  17.     for (j=2; j<=N; j=j*2) {
  18.       l+=(N/j)*(k%2); k=k/2;
  19.     }
  20.     if (l>i) {
  21.       t=x[i];
  22.       x[i]=x[l];
  23.       x[l]=t;
  24.     }
  25.   }
  26.   for (i=0; i<N/2; i++) w[i]=complexd(cos(2*PI/N*i),sin(2*PI/N*i));
  27.   s=1;
  28.   for (i=0; i<mu; i++) { //fft by binary dividing
  29.     for (j=0; j<N; j+=s*2) {
  30.       for (k=0; k<s; k++) {
  31. t=x[j+k+s]*(w[k*N/s/2]);
  32. t1=x[j+k]+t;
  33. t2=x[j+k]-t;
  34. x[j+k]=t1;
  35. x[j+k+s]=t2;
  36.       }
  37.     }
  38.     s*=2;
  39.   }
  40.   for (i=0; i<N; i++) {y[i][0]=x[i].real(); y[i][1]=x[i].imag();}
  41.   delete w; //very important. If not, the system will lost memory
  42.   delete x;
  43. }