fft1.cpp
上传用户:alisonmail
上传日期:2013-02-28
资源大小:500k
文件大小:1k
- #include "stdafx.h"
- #include "complex"
- 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;
- }