fft算法

// FFT_t.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include "complex.h"
#include "cmath"
#include "iostream"
#include "stdio.h"
#include "fstream"
using namespace std;
#include "stdio.h"
#define NUM  128

int changeadd(int k,int m)//倒序方法一，直接算出其值
{
//m是二进制数的位数
//这个地主要注意的是当m很大时，K1可能会超出范围
int k1=0;
for(int i=0;i<m;i++)
{
  if(k&(int)pow(2.0,i))
   k1+=pow(2.0,m-i-1);
}
return k1;
}
void FFT(complex *in,complex * out,int M)
{
ofstream file("fft.dat");
int N=pow(2.0,M);
out[0]=in[0];
out[N-1]=in[N-1];
complex t;
for(int i=1;i<N-1;i++)
{
  if(i<changeadd(i,M))
  {
  t=in[changeadd(i,M)];
  out[changeadd(i,M)]=in[i];
  out[i]=t;
  }
}
file<<"in"<<endl;
for(int i=0;i<N;i++)
  file<<in[i];
//cout<<"out"<<endl;
//for(int i=0;i<N;i++)
// cout<<out[i];
for(int L=1;L<=M;L++)//L第几层蝶形，
{
  int B=pow(2.0,L-1);//B同一个蝶形的两个点间的距离
 
  for(int J=0;J<=B-1;J++)
  {
   int P=pow(2.0,M-L)*J;
  
   complex temp(cos(2*PI*P/N),-sin(2*PI*P/N));
   file<<temp;
   for(int K=J;K<=N-1;K+=pow(2.0,L))
   {
    //if(L==1&&J==0&&K==2){__asm int 3}
    out[K]=out[K]+out[K+B]*temp;
    out[K+B]=out[K]-out[K+B]*temp*2;
   }
  }
  //观察蝶形的具体样子
  //file<<"第"<<L<<"层"<<endl;
  //for(int i=0;i<N;i++)
  //{
  // file<<out[i];
  //}
  //file<<endl;
}
file.close();
}
//////////////////////////////////////////////////////
////在网上找到的一种算法
void _FFT(complex   Input[],int Length,int isign)
{
//isign代表是正还是逆变换 -1表示正变换，1表示逆变换
ofstream _file("_fft.dat");
int   l,i,m,mr=0;  
complex   t;  
float   tm,pisign=isign*PI;  
for(m=1;m<Length;m++)//倒序
{  
  l=Length>>1;  
  while(mr+l>=Length)   l>>=1;  
  mr=mr%l+l;  
  if(mr>m)
  {  
   t=Input[m];  
   Input[m]=Input[mr];  
   Input[mr]=t;  
  }  
}  
_file<<"in"<<endl;
for(int i=0;i<Length;i++)
_file<<Input[i];
/////蝶形计算
l=1;  
while(l<Length)
{  
  for(m=0;m<l;m++)
  {  
   tm=pisign/l*m;  
   _file<<complex(cos(tm),sin(tm));
   for(i=m;i<Length;i+=l<<1)
   {  
    //if(l==1&&m==0&&i==2){__asm int 3}
    t=Input[i+l]*complex(cos(tm),sin(tm));  
    Input[i+l]=Input[i]-t;  
    Input[i]+=t;  
   }  
  }  
  l<<=1;
// _file<<"第"<<l<<"层"<<endl;
// for(int i=0;i<Length;i++)
// {
// 
//  _file<<Input[i];
// }
// _file<<endl;
}  
if(isign==1)   for(l=0;l<Length;l++)   Input[l]/=Length;  
_file.close();
}
//////////////////////////////////////
int _tmain(int argc, _TCHAR* argv[])
{
complex in[NUM];
complex _in[NUM];
double t;
int x=log(NUM*1.0)/log(2.0);
for(int i=0;i<NUM;i++)
{
  t=sin(0.01*2*PI*i)+sin(0.02*2*PI*i)+sin(0.4*2*PI*i);
  _in[i]=in[i]=complex(t,0);
}
FFT(in,in,x);
_FFT(_in,NUM,-1);
//for(int i=0;i<NUM;i++)
// cout<<in[i];
//int x;
//while(1)
//{
// cin>>x;
// if(x==0)break;
// //changeadd;
// printf("%x\n",x);
  //printf("%x\n",changeadd(2,7));
//}
return 0;
}
里边的调试代码我也没有删除，能帮助你理解