FFT简介（来源：百度百科）：

快速傅立叶变换（FFT）是离散傅立叶(DFT)的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。

 

 

Iterative_FFT

参数 n ——a的元素个数

        a[0,1,…,(n-1)]——系数

局部变量

        wd——单位元的n次原根

        w——当前点的估值多项式

        y——变换结果

y = BitReverse（a）

for j = 1 to logn

       d =

       wd =

      w = 1

       for k = 0 to d/2-1

              for m = k to n-1 step d

                     t = w×y[m + d/2]

            x = y[m]

            y[m] = x + t

            y[m + d/2] = x – t

        endfor

        w = w×wd

    endfor

endfor

return y

 

 

C语言实现：

 

void Iterative_FFT(double *ir,double *ii,double *or,double *oi,int len)

{

       int i,j,n,d,k,m;

    double t;

       double wdr,wdi,wr,wi;

       double wtr,wti;

       double tr,ti;

       double xr,xi;

      

       n=pow2(len);

 

       for(i=0;i<n;i++)

              {           

                     j = BitReverse(i,len);

                     or[j] = ir[i];

                     oi[j] = ii[i];

              }

      

       for(j=1;j<=len;j++)

       {

              d  = pow2(j);

 

              t = 2*pi/d;

              wdr = cos(t); wdi = sin(t);

             

              wr = 1; wi = 0;

 

              for (k=0;k<=d/2-1;k++)

              {

                     for (m=k;m<=n-1;m+=d)

                     {

                            tr = wr*or[m + d/2] - wi*oi[m + d/2];

                            ti = wr*oi[m + d/2] + wi*or[m + d/2];

 

                            xr = or[m]; xi = oi[m];

 

                            or[m] = xr + tr; oi[m] = xi + ti;

                            or[m+d/2] = xr - tr; oi[m+d/2] = xi - ti;                           

                     }

 

                     wtr = wr*wdr - wi*wdi;

                     wti = wr*wdi + wi*wdr;

 

                     wr = wtr;

                     wi = wti;

              }    

       }

}

 

FFT的串行递归实现伪代码：

Recursive_FFT

参数 n ——a的元素个数

        a[0,1,…,(n-1)]——系数

局部变量

        wn——单位元的n次原根

        w——当前点的估值多项式

        a0——偶数项系数

        a1——奇数项系数

        y——变换结果

        y0——a0的FFT结果

        y1——a1的FFT结果

 

if n = =1 then return a

else

wn =

w = 1

a0 = (a[0],a[2],…a[n-2])

a1 = (a[2],a[3],…a[n-1])

y0 = Recursive_FFT(a0,n/2)

y1 = Recursive_FFT(a1,n/2)

for k = 0 to n/2-1

       y[k] = y0[k] + w×y1[k]

    y[k + n/2] = y0[k]－w×y1[k]

    w = w×wn

endfor

return y

endif

 

C语言实现：

 

void Recursive_FFT(double *ir,double *ii,double *or,double *oi,int len)

{

       int i,k;

       int n;

       double wnr,wni;

       double wr,wi;

       double wtr,wti;

       double temp0_ir[MAX_N], temp0_ii[MAX_N], temp1_ir[MAX_N], temp1_ii[MAX_N];

       double temp0_or[MAX_N], temp0_oi[MAX_N], temp1_or[MAX_N], temp1_oi[MAX_N];

 

       memset(temp0_ir, 0 ,sizeof(temp0_ir));

       memset(temp0_ii, 0 ,sizeof(temp0_ii));

       memset(temp1_ir, 0 ,sizeof(temp1_ir));

       memset(temp1_ii, 0 ,sizeof(temp1_ii));

       memset(temp0_or, 0 ,sizeof(temp0_oi));

       memset(temp0_or, 0 ,sizeof(temp0_oi));

 

       n=pow2(len);

      

       if(len == 0)

       {

              or = ir;

              oi = ii;

       }

 

       else

       {

              wnr = cos(2*pi/n);   

              wni = sin(2*pi/n);

 

              //printf("wnr:%.2f wni:%.2f\n",wnr,wni);

             

              wr = 1;

              wi = 0;

 

              for(i=0; i<n; i+=2)

              {

                     temp0_ir[i/2] = ir[i];

                     temp0_ii[i/2] = ii[i];

 

                     temp1_ir[i/2] = ir[i+1];

                     temp1_ii[i/2] = ii[i+1];

              }

 

              //for(i=0; i<2; i++)

              //{

              //     printf("%.2f %.2f %.2f %.2f\n",temp0_ir[i],temp0_ii[i],temp1_ir[i],temp1_ii[i]);

              //}

             

              Recursive_FFT(temp0_ir, temp0_ii, temp0_or, temp0_oi, len-1);

              Recursive_FFT(temp1_ir, temp1_ii, temp1_or, temp1_oi, len-1);

 

              for(k=0; k<=n/2-1; k++)

              {

                     or[k] = temp0_or[k] + (wr*temp1_or[k] - wi*temp1_oi[k]);

                     oi[k] = temp0_oi[k] + (wr*temp1_oi[k] + wi*temp1_or[k]);

 

                     or[k+n/2] = temp0_or[k] - (wr*temp1_or[k] - wi*temp1_oi[k]);

                     oi[k+n/2] = temp0_oi[k] - (wr*temp1_oi[k] + wi*temp1_or[k]);

                    

                     wtr = wr*wnr - wi*wni;

                     wti = wr*wni + wi*wnr;

 

                     wr = wtr;

                     wi = wti;

                    

              }

       }

 

}