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），即為所求。