C语言实现FFT(非递归蝶形运算版)

蝶形算法可参考链接
代码如下：

#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点FFT
    complex<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;
}