FFT的算法
FFT是一种DFT的高效算法，称为快速傅立叶变换（fast Fourier transform），它根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。FFT算法可分为按时间抽取算法和按频率抽取算法，先简要介绍FFT的基本原理。从DFT运算开始，说明FFT的基本原理。DFT的运算为：式中 由这种方法计算DFT对于 的每个K值，需要进行4N次实数相乘和（4N-2）次相加，对于N个k值，共需4N*4N次实数相乘和（4N-2)(4N-2）次实数相加。改进DFT算法，减小它的运算量，利用DFT中 的周期性和对称性，使整个DFT的计算变成一系列迭代运算，可大幅度提高运算过程和运算量，这就是FFT的基本思想。FFT对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当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+（N^2）/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
FFT的公式是什么和算法是怎样实现
二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下：
先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩阵的第i行。注意这只是一个简单的记法，并不能完全照抄。还需要通过一些语句来生成各行的数据。同理，COL[i]是对矩阵的第i列的一种简单表示方法。
所以，关键是一维FFT算法的实现。下面讨论一维FFT的算法原理。
【1D-FFT的算法实现】
设序列h(n)长度为N，将其按下标的奇偶性分成两组，即he和ho序列，它们的长度都是N/2。这样，可以将h(n)的FFT计算公式改写如下 ：
(A)
由于
所以，(A)式可以改写成下面的形式：
按照FFT的定义，上面的式子实际上是：
其中，k的取值范围是 0~N-1。
我们注意到He(k)和Ho(k)是N/2点的DFT，其周期是N/2。因此，H(k)DFT的前N/2点和后N/2点都可以用He(k)和Ho(k)来表示
一维实序列的快速傅里叶变换（FFT）
通过前面的分析，我们认识到傅里叶变换本身是复数运算，地球物理获取的数据大多数是实数，对于实数的变换原则上可直接套用复序列的FFT算法，但那样是把实数序列当作虚部为零的复数对待，显然需要存储虚部的零并进行无功的运算，既浪费了一倍的计算内存，又降低了约一半的运算速度。为了不浪费不可不设的虚部内存和必然出现的复数运算，可否将一个实序列分为两个子实序列，分别作为实部与虚部构成一个复数序列，然后用复序列的FFT算法求其频谱，对合成的复序列频谱进行分离和加工得到原实序列的频谱呢？答案是肯定的，实现这一过程思路就是实序列FFT算法的基本思想。1.实序列的傅里叶变换性质对于一个N个样本的实序列x（k），其频谱为X（j），用Xr（j）和Xi（j）表示X（j）的实部和虚部， 表示X（j）的共轭，则 证明：已知 则地球物理数据处理基础上式两端取共轭，并注意到x（k）是实序列，则地球物理数据处理基础这就是实序列的傅里叶变换具有复共轭性。其同样具有周期性，即地球物理数据处理基础2.一维实序列的FFT算法（1）同时计算两个实序列的FFT算法已知两个实序列h（k），g（k）（k＝0，1，…，N－1），例如重磁异常平面数据中的两条剖面，或地震勘探中的两道地震记录，可以人为地构成一个复序列：地球物理数据处理基础设h（k）的频谱为H（j）＝Hr（j）＋iHi（j）g（k）的频谱为G（j）＝Gr（j）＋iGi（j）y（k）的频谱为Y（j）＝Yr（j）＋i Yi（j）利用上节的复序列FFT算法，求得Y（j），即Yr（j）和Yi（j）已知，来寻找Hr（j），Hi（j），Gr（j），Gi（j）与Yr（j），Yi（j）之间的关系。对式（8－22）作傅里叶变换：地球物理数据处理基础由于H（j），G（j）本身是复序列，所以不能仅从上式分离出H（j）和G（j）。应用Y（j）的周期性，容易得到Y（N－j）＝H（－j）＋iG（－j）上式取共轭：地球物理数据处理基础由于h（k），g（k）为实序列，对上式右端应用复共轭定理，得地球物理数据处理基础对式（8－23）展开，得地球物理数据处理基础对式（8－24）展开，并应用共轭关系，得地球物理数据处理基础把式（8－25）和式（8－26）与Y（j）＝Yr（j）＋iYi（j）进行对比，有地球物理数据处理基础整理得地球物理数据处理基础因此，对于两个实序列，通过构造一个复序列，应用复序列的FFT算法和式（8－28）的分离加工，即可得到两个实序列的频谱。（2）计算2 N个数据点的实序列FFT算法设有2N点的实序列u（k）（k＝0，1，…，2N－1），首先按k的偶、奇分成两个子实序列，并构成复序列，即地球物理数据处理基础通过调用复序列FFT算法，求得y（k）的频谱为Y（j）。另记h（k），g（k）的频谱为H（j）和G（j）。利用前面式（8－23）和式（8－24），容易求得地球物理数据处理基础下面分析用H（j），G（j）形成u（k）频谱的问题。记u（k）（k＝0，1，…，2 N－1）的频谱为V（j），分析V（j），H（j），G（j）之间的关系，根据定义地球物理数据处理基础利用式（8－31）和式（8－34）可换算出u（k）的前N个频谱V（j）（j＝0，1，…，N－1），还要设法求u（k）的后N个频谱V（N＋j）（j＝0，1，…，N－1）。利用实序列其频谱的复共轭和周期性：（1）H（N）＝H（0），G（N）＝G（0），WN1＝－1，得地球物理数据处理基础（2）由于u（k）（k＝0，1，…，2N－1）是实序列，同样利用实序列其频谱的复共轭和周期性，用已求出的前N个频谱V（j）表示出后面的N－1个频谱V（N＋j）：地球物理数据处理基础由于0<2N－j<N，所以可从V（j）（j＝0，1，…，N－1）中选出V（2N－j）（j＝N＋1，N＋2，…，2 N－1），并直接取其共轭 即可得到V（N＋1）～V（2 N－1），从而完成整个实序列频谱的计算。总结以上叙述，一维实序列u（k）（k＝0，1，…，2N－1）的FFT计算编程步骤如下：（1）按偶、奇拆分实序列u（k），并构造复序列：地球物理数据处理基础（2）调用复序列的FFT计算y（k）的频谱Y（j）（j＝0，1，…，N－1）；（3）用下式计算形成h（k），g（k）的频谱H（j）和G（j）；地球物理数据处理基础（4）用下式换算实序列u（k）的频谱V（j）（j＝0，1，…，2 N－1）：地球物理数据处理基础［例3］求实序列样本u（k）＝｛1，2，1，1，3，2，1，2｝（k＝0，1，…，7）的频谱。解：按偶、奇拆分实序列u（k），按式（8－37）构造复序列c（j）（j＝0，1，2，3），即c（0）＝1＋2i； c（1）＝1＋i； c（2）＝3＋2i； c（3）＝1＋2i。（1）调用复序列FFT求c（j）（j＝0，1，2，3）的频谱Z（k）（k＝0，1，2，3），得Z（0）＝6＋7i； Z（1）＝－3； Z（2）＝2＋i； Z（3）＝－1。地球物理数据处理基础（3）运用公式（8－38）计算H（j），G（j）：地球物理数据处理基础（4）根据式（8－39）求出u（k）（k＝0，1，…，7）的8个频谱V（j）（j＝0，1，…，7）：地球物理数据处理基础地球物理数据处理基础由上例可见，完成全部2 N个实序列频谱的计算只需做N次FFT计算，相比直接用复序列的FFT算法节省了约一半的计算量。
请用MATLAB、C语言或者其他语言编程实现8点序列的基2-DIT－FFT算法，并对结果进行分析验证。
typedef struct
{
float real;
float img;
}COMPLEX;
void reverse(COMPLEX *X,COMPLEX *x,int N) /*输入数组名称和元素个数，实现倒序 */
{
int LH = N/2;
int j = LH;
int N1 = N-2;
int i;
for(i = 1; i <= N1; i++)
{
int k = LH;
X[i] = (i<j ? x[j] : x[i]);
/*
if(i < j)
{
comp T = a[j];
a[j] = a[i];
a[i] = T;
*/
}
while(j >= k)
{
j = j-k;
k = k/2;
}
j = j+k;
}
}
/*
函数功能：对指定长度的采样数据进行快速傅立叶变换，并返回变换值, X = FFT(x)
入口参数：目标采样数据，保存变换值的指针，采样数据长度
出口参数：指向变换值的指针
注意事项：保存结果的空间是在函数外部申请，本函数不处理此内容
/*/
COMPLEX *FFT(COMPLEX *X,COMPLEX *x, unsigned N)
{
unsigned temp=N,L=0,M=0,LE=1,LE1=0,J=0,I=0,IP=0;
COMPLEX U={0,0},W={0,0},T={0,0};
while((temp>>=1)>0) /*获得阶码M*/
{
M++;
}
reverse(X,x,N);
for(L=0;L<M;L++)
{
LE=1<<L; /*LE=2^L,分组间隔*/
LE1 = LE>>1; /*LE1=LE/2，对偶跨距，实际也是一个分组中对偶点的对数*/
U.real = 1.0;
U.img = 1.0;
W.real = cos(PI / LE1);
W.img = sin(PI / LE1);
J = 0;
while(J<LE1)
{
I = J;
while(I < N)
{
IP = I + LE1;
ComMul(X + IP,&U,&T);
ComSub(X + I,&T,X + IP);
ComAdd(X + I,&T,X + I);
I += LE;
}
ComMul(&U,&W,&U);
J++;
}
}
return X;
}
FFT中有几个复数运算，自己实现，不想发给你，年轻人还是要自己动手做点东西。N为任意数，正常应该为2的幂次方。
