一维实序列的快速傅里叶变换（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）
设x（N）为N点有限长离散序列，代入式（8－3）、式（8－4），并令 其傅里叶变换（DFT）为地球物理数据处理基础反变换（IDFT）为地球物理数据处理基础两者的差异只在于W的指数符号不同，以及差一个常数1/N，因此下面我们只讨论DFT正变换式（8－5）的运算量，其反变换式（8－6）的运算是完全相同的。一般来说，W是复数，因此，X（j）也是复数，对于式（8－5）的傅里叶变换（DFT），计算一个X（j）值需要N次复数乘法和N－1次复数加法。而X（j）一共有N个值（j＝0，1，…，N－1），所以完成整个DFT运算总共需要N2次复数乘法和N（N－1）次复数加法。直接计算DFT，乘法次数和加法次数都是与N2成正比的，当N很大时，运算量会很大，例如，当N＝8时，DFT需64次复数乘法；而当N＝1024时，DFT所需乘法为1048576次，即一百多万次的复数乘法运算，对运算速度要求高。所以需要改进DFT的计算方法，以减少运算次数。分析Wjk，表面上有N2个数值，由于其周期性，实际上仅有N个不同的值W0，W1，…，WN－1。对于N＝2m时，由于其对称性，只有N/2个不同的值W0，W1，…，地球物理数据处理基础因此可以把长序列的DFT分解为短序列DFT，而前面已经分析DFT与N2成正比，所以N越小越有利。同时，利用ab＋ac＝a（b＋c）结合律法则，可以将同一个Wr对应的系数x（k）相加后再乘以Wr，就能大大减少运算次数。这就是快速傅里叶变换（FFT）的算法思路。下面，我们来分析N＝2m情况下的FFT算法。1.N＝4的FFT算法对于m＝2，N＝4，式（8－5）傅里叶变换为地球物理数据处理基础将式（8－7）写成矩阵形式地球物理数据处理基础为了便于分析，将上式中的j，k写成二进制形式，即地球物理数据处理基础代入式（8－7），得地球物理数据处理基础分析Wjk的周期性来减少乘法次数地球物理数据处理基础则 代回式（8－9），整理得地球物理数据处理基础上式可分层计算，先计算内层，再计算外层时就利用内层计算的结果，可避免重复计算。写成分层形式地球物理数据处理基础则X（j1 j0）＝X2（j1 j0）。上式表明对于N＝4的FFT，利用Wr的周期关系可分为m＝2步计算。实际上，利用Wr的对称性，仍可以对式（8－11）进行简化计算。考虑到地球物理数据处理基础式（8－11）可以简化为地球物理数据处理基础令j＝j0；k＝k0，并把上式表示为十进制，得地球物理数据处理基础可以看到，完成上式N＝4的FFT计算（表8－1）需要N·（m－1）/2＝2次复数乘法和N·m＝8次复数加法，比N＝4的DFT算法的N2＝16次复数乘法和N·（N－1）＝12次复数加法要少得多。表8－1 N＝4的FFT算法计算过程注：W0＝1；W1＝－i。［例1］求N＝4样本序列1，3，3，1的频谱（表8－2）。表8－2 N＝4样本序列2.N＝8的FFT算法类似N＝4的情况，用二进制形式表示，有地球物理数据处理基础写成分层计算的形式：地球物理数据处理基础则X（j2 j1 j0）＝X3（j2 j1 j0）。对式（8－14）的X1（k1 k0 j0）进行展开，有地球物理数据处理基础还原成十进制，并令k＝2k1＋k0，即k＝0，1，2，3，有地球物理数据处理基础用类似的方法对式（8－14）的X2（k0 j1 j0），X3（j2 j1 j0）进行展开，整理得地球物理数据处理基础用式（8－16）、式（8－17）逐次计算到X3（j）＝X（j）（j＝0，1，…，7），即完成N＝23＝8的FFT计算，其详细过程见表8－3。表8－3 N＝8的FFT算法计算过程注：对于正变换 对于反变换 所 ［例2］求N＝8样本序列（表8－4）x（k）＝1，2，1，1，3，2，1，2的频谱。表8－4 N＝8样本序列3.任意N＝2m的FFT算法列出N＝4，N＝8的FFT计算公式，进行对比地球物理数据处理基础观察式（8－18）、式（8－19），不难看出，遵循如下规律：（1）等式左边的下标由1递增到m，可用q＝1，2，…，m代替，则等式右边为q－1；（2）k的上限为奇数且随q的增大而减小，至q＝m时为0，所以其取值范围为k＝0，1，2，…，（2m－q－1）；（3）j的上限为奇数且随q的增大而增大，且q＝1时为0，其取值范围为j＝0，1，2，…，（2q－1－1）；（4）k的系数，在等式左边为2q，等式右边为2q－1（包括W的幂指数）；（5）等式左边序号中的常数是2的乘方形式，且幂指数比下标q小1，即2q－1；等式右边m对式子序号中的常数都是定值2m－1。归纳上述规则，写出对于任意正整数m，N＝2m的FFT算法如下：由X0（p）＝x（p）（p＝0，1，…，N－1）开始：（1）对q＝1，2，…，m，执行（2）～（3）步；（2）对k＝0，1，2，…，（2m－q－1）及j＝0，1，2，…，（2q－1－1），执行地球物理数据处理基础（3）j，k循环结束；（4）q循环结束；由Xm（p）（p＝0，1，…，N－1）输出原始序列x（p）的频谱X（p）。在计算机上很容易实现上述FFT算法程序，仅需要三个复数数组，编程步骤如下：（1）设置复数数组X1（N－1），X2（N－1）和 （数组下界都从0开始）；（2）把样本序列x赋给X1，即X1（k）＝x（k）（k＝0，1，…，N－1）；（3）计算W，即正变换 反变换 （4）q＝1，2，…，m，若q为偶数，执行（6），否则执行第（5）步；（5）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作X2（2qk＋j）＝X1（2q－1k＋j）＋X1（2q－1k＋j＋2m－1）X2（2qk＋j＋2q－1）＝［X1（2q－1k＋j）－X1（2q－1k＋j＋2m－1）］W（2q－1k）至k，j循环结束；（6）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作X1（2qk＋j）＝X2（2q－1k＋j）＋X2（2q－1k＋j＋2m－1）X1（2qk＋j＋2q－1）＝［X2（2q－1k＋j）－X2（2q－1k＋j＋2m－1）］W（2q－1k）至k，j循环结束；（7）q循环结束，若m为偶数，输出X1（j），否则输出X2（j）（j＝0，1，…，N－1），即为所求。
关于二维傅里叶变换频率分布示意图
信号经过一维傅里叶变换后，其对应的数字频率域的范围为【0，2*pi）；
一般情况下，我们所观测的信号都是实信号，并不是解析信号，同时由于数字频率的周期性，所以数字频率域的范围必须是【0，pi】并（-pi，0），因此最开始和结束都在零频率附近，所以都认为是低频。
二维情况下基本相同，
左下角是x方向的正低频和y方向的正低频；
左上角是x方向的正低频和y方向的负低频；
右下角是x方向的负低频和y方向的正低频；
右上角是x方向的负低频和y方向的负低频；
所以，这四个角全部是低频~
要想看的舒服的话，可以用Matlab中的fftshift函数~ 这样中心就是低频，外围全部是高频了
