蝶形算法可参考链接
代码如下：
#include<iostream>#include<complex.h>#include<math.h>#define PI 3.14159#define N 32#define STAGE ((int)(log(N)/log(2)))using namespace std;typedef complex<double> data_t;void FFT(data_t Xin[N],data_t Xout[N]){//旋转因子计算data_t W[N/2];for(int i=0;i<N/2;i++){W[i].real(cos(2*PI*i/N));W[i].imag(-sin(2*PI*i/N));}//中间变量data_t Xtmp[STAGE+1][N];for(int i=0;i<N;i++){Xtmp[0][i]=Xin[i];}//计算for(int i=0;i<STAGE;i++) //共i个stage{//第i个stage,根据Xtmp[i]计算Xtmp[i+1]int span=1<<i;//跨度=蝶形单元大小的一半int group_num=(int)N/(2*span); //蝶形单元数目int group_size=(int)N/group_num; //蝶形单元大小for(int k=0;k<group_num;k++){for(int t=0;t<group_size/2;t++){//需要group_size次单位根的t次方data_t w={cos(2*PI*t/group_size),-sin(2*PI*t/group_size)};Xtmp[i+1][k*group_size+t] = Xtmp[i][k*group_size+t]+w*Xtmp[i][k*group_size+t+span];Xtmp[i+1][k*group_size+t+group_size/2] = Xtmp[i][k*group_size+t]-w*Xtmp[i][k*group_size+t+span];}}}for(int i=0;i<N;i++){Xout[i]=Xtmp[STAGE][i];}}void bit_reverse(data_t* Xin,int n){int fft8table[8]={0,4,2,6,1,5,3,7};int fft16table[16]={0,8,4,12,2,10,6,14,1,9,5,13,3,11,7,15};int fft32table[32]={0,16,8,24,4,20,12,28,2,18,10,26,6,22,14,30,1,17,9,25,5,21,13,29,3,19,11,27,7,23,15,31};data_t *Tmp=new data_t[n];switch(n){case 8:{for(int i=0;i<8;i++)Tmp[i]=Xin[fft8table[i]];for(int i=0;i<8;i++)Xin[i]=Tmp[i];break;}case 16:{for(int i=0;i<16;i++)Tmp[i]=Xin[fft16table[i]];for(int i=0;i<16;i++)Xin[i]=Tmp[i];break;}case 32:{for(int i=0;i<32;i++)Tmp[i]=Xin[fft32table[i]];for(int i=0;i<32;i++)Xin[i]=Tmp[i];break;}default:cout<<"error value of n"<<endl;}}int main(){//支持8,16,32点FFTcomplex<double> X[N];complex<double> Y[N];for(int i=0;i<N;i++){X[i].real(i);X[i].imag(0);}bit_reverse(X,N);FFT(X,Y);for(int i=0;i<N;i++)cout<<Y[i]<<endl;return 0;}
