dft计算

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算法节省了约一半的计算量。

二维实序列的快速傅里叶变换（FFT）

在地球物理数据处理中，经常遇到处理二维实数据的情况。例如在地震勘探中，对面波勘探数据作频散分析解释时，要将时间－空间域的信息转换为频率－波数域频谱；在重磁异常的滤波或转换中，要将空间域的异常f（x，y）转换为波数域F（ω，υ）等。这些分析都需要进行二维的傅里叶变换（FFT）。根据傅里叶变换的定义，对于连续二维函数f（x，y），其傅里叶变换对为地球物理数据处理基础对于离散的二维序列fjk（j＝0，1，…，M－1；k＝0，1，…，N－1），其傅里叶变换为地球物理数据处理基础1.二维复序列的FFT算法对于M条测线，每条测线N个测点，构成复序列yjk（j＝0，1，…，M－1；k＝0，1，…，N－1），根据离散傅里叶公式（8－41），其傅里叶变换为地球物理数据处理基础于是，可以分两步套用一维复FFT完成二维复FFT的计算。（1）沿测线方向计算对于j＝0，1，…，M－1逐测线套用一维复FFT，执行式（8－43）。定义复数组 则算法为1）对于j＝0，1，…，M－1，作2）～7）；2）将yjk输入A1（k），即A1（k）＝yjk（k＝0，1，…，N－1）；3）计算Wr，存入W（r），即 4）q＝1，2，…，p（p＝log2N），若q为偶数执行6），否则执行5）；5）k＝0，1，2，…，（2p－q－1）和n＝0，1，2，…，（2q－1－1）循环，作A2（k2q＋n）＝A1（k2q－1＋n）＋A1（k2q－1＋n＋2p－1）A2（k2q＋n＋2q－1）＝［A1（k2q－1＋n）－A1（k2q－1＋n＋2p－1）］·W（k2q－1）至k，n循环结束；6）k＝0，1，2，…，（2p－q－1）和n＝0，1，2，…，（2q－1－1）循环，作A1（k2q＋n）＝A2（k2q－1＋n）＋A2（k2q－1＋n＋2p－1）A1（k2q＋n＋2q－1）＝［A2（k2q－1＋n）－A2（k2q－1＋n＋2p－1）］·W（k2q－1）至k，n循环结束；7）q循环结束，若p为偶数，将A1（n）输入到Yjn，否则将A2（n）输入到Yjn（n＝0，1，…，N－1）；8）j循环结束，得到Yjn（j＝0，1，…，M－1；n＝0，1，…，N－1）。（2）垂直测线方向计算对于n＝0，1，…，N－1逐一套用一维复FFT，执行式（8－44）。即1）对于n＝0，1，…，N－1，作2）～7）；2）将Yjn输入A1（j），即A1（j）＝Yjn（j＝0，1，…，M－1）；3）计算Wr存入W（r），即 4）q＝1，2，…，p（p＝log2M），若q为偶数执行6），否则执行5）；5）j＝0，1，2，…，（2p－q－1）和m＝0，1，2，…，（2q－1－1）循环，作A2（j2q＋m）＝A1（j2q－1＋m）＋A1（j2q－1＋m＋2p－1）A2（j2q＋m＋2q－1）＝［A1（j2q－1＋m）－A1（j2q－1＋m＋2p－1）］·W（j2q－1）至j，m循环结束；6）j＝0，1，2，…，（2p－q－1）和m＝0，1，2，…，（2q－1－1）循环，作 A1（j2q＋m）＝A2（j2q－1＋m）＋A2（j2q－1＋m＋2p－1） A1（j2q＋m＋2q－1）＝［A2（j2q－1＋m）－A2（j2q－1＋m＋2p－1）］·W（j2q－1）至j，m循环结束；7）q循环结束，若p为偶数，将A1（m）输入到Ymn，否则将A2（m）输入到Ymn（m＝0，1，…，M－1）；8）n循环结束，得到二维复序列的傅氏变换Ymn（m＝0，1，…，M－1；n＝0，1，…，N－1），所求得的Ymn是复数值，可以写为Ymn＝Rmn＋iImn （m＝0，1，…，M－1；n＝0，1，…，N－1）其中，Rmn，Imn的值也是已知的。2.二维实序列的FFT算法对于二维的实序列，我们把其看作是虚部为零的复序列，套用上述的二维复序列FFT方法来求其频谱算法上也是可行的，但势必会增加大量的无功运算。因此，有必要研究二维实序列FFT的实用算法，同一维实序列FFT的实现思路一样，同样把二维实序列按一定的规律构造成二维复序列，调用二维复序列FFT，然后通过分离和加工得到原实序列的频谱。例如采样区域有2 M条测线，每条测线有N个点，并且M，N都是2的整数幂，需要计算实样本序列xjk（j＝0，1，2，…，2 M－1；k＝0，1，2，…，N－1）的傅氏变换：地球物理数据处理基础类似于一维实序列FFT的思想，直接建立下面的二维实序列FFT算法：（1）将一个二维实序列按偶、奇线号分为两个二维子实序列，分别作为实部和虚部组合为一个二维复序列。即令地球物理数据处理基础（2）调用二维复FFT过程，求出yjk的二维傅氏变换Ymn的复数值：地球物理数据处理基础式中：Rmn，Imn是Ymn的实部和虚部。（3）利用Rmn，Imn换算Xmn的值。前两步容易实现，下面分析第（3）步的实现。记hjk，gjk的傅氏变换为Hmn，Gmn。根据傅里叶变换的定义，我们导出Xmn与Hmn，Gmn的关系式：地球物理数据处理基础式中，Hmn，Gmn为复数，我们用上标r和i表示其实部和虚部，将上式右端实部、虚部分离地球物理数据处理基础其中：地球物理数据处理基础下面的任务是将Hmn，Gmn各分量与通过二维复FFT求出的Rmn，Imn值联系起来。为此先给出奇、偶分解性质和类似于一维情况的三个二维傅氏变换性质：（1）奇偶分解性任何一个正负对称区间定义的函数，均可唯一地分解为如下偶（even）、奇（odd）函数之和：地球物理数据处理基础（2）周期性地球物理数据处理基础（3）复共轭性地球物理数据处理基础现在我们来建立Rmn，Imn与Hmn，Gmn的关系。对Ymn作奇偶分解：地球物理数据处理基础根据线性性质地球物理数据处理基础对照式（8－54）和式（8－55），得地球物理数据处理基础由于hjk，gjk是实函数，根据复共轭性质，上面两式对应的奇偶函数相等。即地球物理数据处理基础再由奇偶分解性和周期性，得地球物理数据处理基础将式（8－57）代入式（8－50），得地球物理数据处理基础再利用Hmn，Gmn周期性及复共轭性，可以得到m＝M/2＋1，…，M－1；n＝0，1，…，N－1的傅氏变换，即地球物理数据处理基础将式（8－50）中M，N改为M－m，N－n，并将上式代入，得地球物理数据处理基础由式（8－58）、式（8－59）和式（8－61）即可得到原始序列xjk（j＝0，1，…，2M－1；n＝0，1，…，N－1）在m＝0，1，…，M－1；n＝0，1，…，N－1区间的傅氏变换Xmn。具体二维实序列的FFT算法如下：（1）令hjk＝x2j，k，gjk＝x2j＋1，k，形成yjk＝hjk＋igjk （j＝0，1，…，2 M－1；n＝0，1，…，N－1）（2）调用二维复序列FFT过程，即从两个方向先后调用一维复FFT算法式（8－43）和式（8－44），求得yjk的二维傅氏变换Ymn的复数值：Ymn＝Rmn＋iImn （m＝0，1，…，M－1；n＝0，1，…，N－1）（3）用下列公式由Rmn，Imn的值换算Xmn的值：地球物理数据处理基础地球物理数据处理基础