首先，介绍下FFT算法。
参考：
http://www.gatevin.moe/acm/fft%E7%AE%97%E6%B3%95%E5%AD%A6%E4%B9%A0%E7%AC%94%E8%AE%B0/
博客写的很好，也很基础。

简单介绍以下FFT：

为什么需要FFT

　　第一个问题是为什么要创造FFT，简单的说，为了速度。我们承认DFT很有用，但是我们发现他的速度不是很快，1D的DFT原始算法的时间复杂度是O(n^2)，这个可以通过公式观察出来，对于2D的DFT其时间复杂度是O（n^4），这个速度真的很难接受，也就是说，你计算一幅1024x768的图像时，你将等大概。。。大概。。。我也没试过，反正很久。不信的自己去试试。所以找到一种快速方法的方法计算FFT势在必行。
　　以下为DFT公式：
　　计算一个4点DFT。计算量如下：

如何得到FFT

　　通过观察DFT公式，我们发现DFT计算每一项X[u]的时候都遍历了完整的x[n]所以，我们的想法就是能不能通过其他的X[u+m]（m为一个整数，可正可负）得到X[u]，也就是充分利用之间的计算结构来构建现在的结果，这种方法就很容易表现成迭代的形式。本文介绍基2的FFT，及离散信号x[n]的个数为2的k次方，即如果完整的离散信号中有N个信分量{x1，x2，x3....xN}其中`(N=1< 　　数学基础
　　FFT的数学基础其实就是：


　　旋转因子具有以下性质：

　　这些性质使得我们可以利用前面的计算结果来完成出后续的计算。

2-FFT

　　观察：一个只有两个值的离散信号，假设N=2，利用性质2-对称性可以得到

4-FFT

　　与上面2点的一样，推广到4点，将会是这样，其中方框内的操作为上述2-FFT： 　　最终算法：

8-FFT

　　废话不多说，和上面一样：
　　计算过程为：

　　结果为：

2^n-FFT

　　同理，推广到2^n，可以得到类似的结果，不过我们发现，为了使输出结果为顺序结果，输入的顺序经过了一系列的调整，而且每一步WN的幂次参数也是变化，我们必须得到其变化规律才能更好的编写程序：
　　观察WN的变化规律为：
　　节点距离（设为z）就是从WN的0次幂开始连续加到WN的z次幂。然后间隔为z个式子，再次从0次幂加到z次。重复以上过程：
　　例如第二级，间隔z=2，节点为实心黑 {MOD}点，节点0,1，不做操作，节点2*W8^0,节点3*W8^2，节点4，5不做操作。。
　　其内在规律就是节点i是否乘以WN：
　　if（i%z==奇数）
　　节点i*WN^(step)；
　　step每次增加的数量由当前的所在的计算级决定
　　step=1<<（k（总级数）-i（当前级数））
　　输入参数重新排序
　　i行表示数组下标，蓝 {MOD}数字表示实际数据：其内在规律就是，下标为偶数的将被映射到自己本身下标的1/2处，下表为奇数时被映射到数组后半段（size_n/2）的（下标-1）1/2处，将排列后的数据分为前后两个部分，递归次过程，直到只有两个元素。则停止。过程如下：

　　观察算法
　　观察至此，我们已经基本把FFT蝶形算法的所有特征都搞定了，接下来就是使用代码来实现了。

fft.c #include "math.h" #include "fft.h" //精度0.0001弧度 void conjugate_complex(int n,complex in[],complex out[]) { int i = 0; for(i=0;iimag = -in[i].imag; out[i].real = in[i].real; } } void c_abs(complex f[],float out[],int n) { int i = 0; float t; for(i=0;ireal * f[i].real + f[i].imag * f[i].imag; out[i] = sqrt(t); } } void c_plus(complex a,complex b,complex *c) { c->real = a.real + b.real; c->imag = a.imag + b.imag; } void c_sub(complex a,complex b,complex *c) { c->real = a.real - b.real; c->imag = a.imag - b.imag; } void c_mul(complex a,complex b,complex *c) { c->real = a.real * b.real - a.imag * b.imag; c->imag = a.real * b.imag + a.imag * b.real; } void c_div(complex a,complex b,complex *c) { c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag); c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag); } #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void Wn_i(int n,int i,complex *Wn,char flag) { Wn->real = cos(2*PI*i/n); if(flag == 1) Wn->imag = -sin(2*PI*i/n); else if(flag == 0) Wn->imag = -sin(2*PI*i/n); } //傅里叶变化 void fft(int N,complex f[]) { complex t,wn;//中间变量 int i,j,k,m,n,l,r,M; int la,lb,lc; /*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----按照倒位序重新排列原信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i/2; while(k<=j) { j=j-k; k=k/2; } j=j+k; } /*----FFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离 /*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,M-m); for(n=l-1;n-1;n=n+la) //遍历每个分组，分组总数为N/la { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 Wn_i(N,r,&wn,1);//wn=Wnr c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算 c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr } } } } //傅里叶逆变换 void ifft(int N,complex f[]) { int i=0; conjugate_complex(N,f,f); fft(N,f); conjugate_complex(N,f,f); for(i=0;iimag = (f[i].imag)/N; f[i].real = (f[i].real)/N; } } fft.h #ifndef __FFT_H__ #define __FFT_H__ typedef struct complex //复数类型 { float real; //实部 float imag; //虚部 }complex; #define PI 3.1415926535897932384626433832795028841971 /////////////////////////////////////////// void conjugate_complex(int n,complex in[],complex out[]); void c_plus(complex a,complex b,complex *c);//复数加 void c_mul(complex a,complex b,complex *c) ;//复数乘 void c_sub(complex a,complex b,complex *c); //复数减法 void c_div(complex a,complex b,complex *c); //复数除法 void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中 void ifft(int N,complex f[]); // 傅里叶逆变换 void c_abs(complex f[],float out[],int n);//复数数组取模 //////////////////////////////////////////// #endif