fft编程
㈠ 一维实序列的快速傅里叶变换（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
FFT，快速傅里叶变换，是处理信号和数据分析的关键工具，其在数学、工程和科学领域发挥着重要作用。本文深入探讨了FFT的基本理论、物理意义以及工程应用，帮助读者全面理解FFT的内在逻辑和实际应用。
理论介绍中，我们首先回顾了傅里叶变换与离散傅里叶变换的基本概念。傅里叶变换基于连续信号，而实际应用中，数据是离散的，因此必须将傅里叶变换离散化。离散傅里叶变换（DFT）的公式体现了这一转换过程，旋转因子是其核心计算元素，用于将连续信号转换为频率域表示。
接着，我们介绍了快速傅里叶变换（FFT）的概念。FFT是DFT的高效算法，通过利用旋转因子的周期性、对称性等性质，大幅度减少了计算量。当采样点数为2的幂时，FFT的计算速度远超DFT，使得实际应用中能够快速获得信号的频谱信息。
FFT的物理意义在于揭示信号在频域的特性。任何信号都可以被分解为不同频率的正弦波组合，通过FFT处理，可以直观地观察到信号的频谱结构。在工程应用中，频谱分析是信号处理、通信、控制等领域中的重要工具。例如，在飞控设计中，对陀螺仪和加速度信号进行频谱分析，有助于滤除机架振动的噪声，优化飞机的控制性能。
工程应用中，FFT的编程实现依赖于对旋转因子的高效利用。通过递归方式，将原始序列拆分为蝴蝶形运算，每一层的运算都利用前一层的结果，最终得到所需的频谱信息。这一过程可以通过编程实现，具体细节可以参阅相关文献或教程。
结果分析部分，我们讨论了FFT的输出及其解释。使用FFT时，需要选择合适的采样点数，这决定了后续分析的精度和计算量。N个采样点经过FFT后，将得到N个频谱点，其中一半的点在实际应用中可能不直接相关，因此通常仅关注N/2+1个关键点。这些点对应于信号在不同频率上的幅度和相位信息，对于信号分析和滤波设计至关重要。
在实际应用中，通过FFT分析信号的频率特性，可以有效识别并处理噪声、干扰等，提高信号处理的准确性和可靠性。以一个示例信号为例，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），即为所求。