跳轉到

6.3 FFT 演算法 1965#

本章導言#

1965 年 4 月，《Mathematics of Computation》刊登了一篇題為〈用於複雜傅立葉級數之機器計算的算法〉（An algorithm for the machine calculation of complex Fourier series）的短文，作者 James W. Cooley（IBM）與 John W. Tukey（普林斯頓/貝爾實驗室統計學家）。這篇 5 頁論文描述的快速傅立葉變換（Fast Fourier Transform，FFT）把 \(N\) 點 DFT 的運算量從 \(O(N^2)\) 降到 \(O(N\log N)\) —— \(N = 10^6\) 時快約 5 萬倍。FFT 不是新變換（變換就是 DFT），而是同一矩陣的因式分解；但正是它使傅立葉分析從「理論工具」變成「文明基礎設施」：數位通訊、醫學成像、音訊壓縮、氣象預報、引力波探測，全部站在 \(O(N\log N)\) 之上。本節講述 FFT 的歷史（包括冷戰諜報式的催生過程）、分治原理、以及它對 20 世紀後半科學的系統性衝擊。

歷史：從 Gauss 到冷戰的地震偵測#

前史（1805–1965）：Gauss 1805 手稿（6.2 節）、Runge 1903、Danielson–Lanczos 1942（X 光晶體學的 DFT 計算，已含分治思想）——同一想法被反覆發明至少五次，每次都因「沒有計算需求」而沉寂。Cooley–Tukey 的不同之處在於時代：電腦已存在，而需求正在爆炸。

催生的地緣政治：1960 年代美蘇核試談判（部分禁試條約 1963）的關鍵技術是地震波偵測——分辨地震與地下核試需要對地震儀訊號做譜分析。問題規模巨大（長序列、多測站）， \(O(N^2)\) 的 DFT 成為瓶頸。Tukey 在總統科學顧問委員會的會議上得知此需求（據傳他先以手算驗證了分治思路），Garwin（氫彈設計者、時任 IBM 顧問）把問題帶給 Cooley——1965 年論文因此被學界戲稱為「冷戰的副產品」。更有趣的一層：Tukey 後來回憶，推動的動機之一是驗證蘇聯是否違反禁試——FFT 一出生就參與了軍備控制。Cooley 的實作與傳播（1967 年 IEEE 論文、公開 FORTRAN 程式）使 FFT 迅速普及；Tukey 的統計學視野（他同時發明了 box plot、快速排序的 early ideas）保證了算法的「可傳授性」。

Gauss 優先權的發現：1984 年 Heideman–Johnson–Burrus 的文獻考古確認 Gauss 在 1805 年（Kronecker 1861 編輯其遺稿出版）已寫出等價算法——Cooley–Tukey 論文的第一句引言便值得重讀：他們當時並不知道 Runge 與 Danielson–Lanczos 的先例。科學史在此展示「發明 = 想法 × 時代需求 × 傳播系統」的乘積律：Gauss 有想法無需求，Danielson–Lanczos 有想法有需求（晶體學）但傳播不廣，Cooley–Tukey 三者俱全。

分治原理：矩陣的因式分解#

DFT 的困境： \(X = F_N x\) 直接計算需 \(N^2\) 次複數乘加—— \(F_N\) 表面上是稠密矩陣。

觀察（偶奇分解）：把 \(N = 2M\) 點輸入按奇偶下標拆成兩個 \(M\) 點序列 \(x^{(e)}_m = x_{2m}\) 、 \(x^{(o)}_m = x_{2m+1}\) ，則

\[ X_k = \sum_{m=0}^{M-1}x_{2m}e^{-2\pi i\frac{2mk}{N}} + \sum_{m=0}^{M-1}x_{2m+1}e^{-2\pi i(2m+1)k/N} = E_k + e^{-\frac{2\pi i k}{N}}\,O_k, \]

其中 \(E, O\) 是長 \(M\) 的 DFT。關鍵：利用 \(e^{-2\pi i(k+M)/N} = -e^{-2\pi ik/N}\) ，後半段輸出可重複使用同一對子 DFT：

\[ X_{k+M} = E_k - e^{-\frac{2\pi i k}{N}}\,O_k, \qquad k = 0, \ldots, M-1. \]

兩個 \(M\) 點 DFT（ \(2M^2\) 次乘加）+ \(M\) 次蝶形合成 = 約一半的代價。遞迴展開： \(T(N) = 2T(N/2) + O(N) \Rightarrow T(N) = O(N\log N)\) 。矩陣語言： \(F_N = P\begin{pmatrix}F_{N/2} & \\ & F_{N/2}\end{pmatrix}\cdot(\text{蝶形因子})\cdots\) ——稠密矩陣的隱藏稀疏因式分解（ \(\log_2 N\) 個因子、每個因子僅 \(N/2\) 個非平凡元素「旋轉因子」）。 \(\log_2 N\) 個因子的乘積，這就是「快」的全部秘密。

蝶形圖： \(X_k = E_k + w^k O_k\) 、 \(X_{k+M} = E_k - w^k O_k\) 的形狀如蝶（兩進兩出、一個旋轉因子 \(w = e^{-2\pi i/N}\) ）——蝶形是數位訊號處理的圖騰，硬體 FFT 晶片就是蝶形單元的流水線。

變體景觀：（1）Radix-2/4（ \(N = 2^m\) 最常見）；（2）混合基（ \(N\) = 高度合成數）；（3）質數長度用 Rader/Winograd 或 Bluestein（chirp-z）；（4）實數輸入用 \(N/2\) 點複數 FFT 合成（省一半）；（5）原地計算（位反轉排列 \(+\) 同址蝶形，記憶體 \(O(N)\) ）。工程上最常用的實現（如 FFTW 庫）是自適應組合：編譯時按機器特性自動選擇分解路徑——「最快的 FFT」是個參數族，不是單一算法。

計算複雜度的意義#

\(N\) 直接 \(N^2\) FFT \(N\log_2 N\) 加速比
1024 \(10^6\) \(10^4\) ×100
\(10^6\) \(10^{12}\) \(2\times10^7\) ×50,000
\(10^9\) \(10^{18}\) \(3\times10^{10}\) ×3×10⁷

複雜度從 \(N^2\) 到 \(N\log N\) 是「可行/不可行」的質變： \(N = 10^9\) （今日標準解析度）下，直接法需數十年單核時間，FFT 秒級完成。漸進複雜度不是學術裝飾，是文明能力的邊界。附帶的數值性質：FFT 的舍入誤差增長 \(O(\sqrt{\log N})\) （隨機模型下）——不僅快，還異常穩定（ \(F_N\) 接近酉矩陣，誤差不被放大）。

一個具體例題#

例題 6.3.1：手算 \(N = 4\) 的 radix-2 FFT。（1）寫出偶奇分解與蝶形步驟；（2）統計乘法次數對比直接法；（3）驗證與 6.2 節練習 1（ \(x = (1,2,3,4)\) 的 DFT \(= (10, -2+2i, -2, -2-2i)\) ）一致。

解：

（1）偶序列 \(x^{(e)} = (x_0, x_2) = (1, 3)\) 、奇序列 \(x^{(o)} = (x_1, x_3) = (2, 4)\) 。兩個 2 點 DFT（2 點 DFT = 減法， \(F_2 = \begin{pmatrix}1&1\\1&-1\end{pmatrix}\) ）： \(E = (1+3,\ 1-3) = (4, -2)\) 、 \(O = (2+4,\ 2-4) = (6, -2)\) 。旋轉因子 \(w = e^{-2\pi i/4} = -i\) 。蝶形：

\[ X_0 = E_0 + w^0 O_0 = 4 + 6 = 10,\qquad X_2 = E_0 - O_0 = -2, $$ $$ X_1 = E_1 + w^1 O_1 = -2 + (-i)(-2) = -2 + 2i,\qquad X_3 = E_1 - w^1 O_1 = -2 - 2i. \]

（2）乘法次數：兩個 2 點 DFT（0 乘法，僅加減）+ 2 次蝶形乘法（ \(w^0\) 可省）= 2 次乘法（含 trivial 因子），直接法 \(N^2 = 16\) 次。 \(N = 4\) 的理論省率 \(\frac{N\log_2N}{N^2} = \frac{8}{16} = 50\%\) ✓。

（3）結果與直接計算完全一致 ✓。注意 \(X_2 = -2\) 實數、 \(X_1, X_3\) 共軛——實訊號的共軛對稱（6.2 節）自然湧現。

歷史啟示#

  • 算法即科學的放大器：FFT 之後，「運算量」成為科學問題可行性的判據之一。20 世紀後半科學的「計算轉向」，其加速器就是少數幾個 \(O(N\log N)\) 算法（FFT、快速排序、快速多極）。Tukey 的統計學家身份提醒我們：算法思維是跨學科的通用能力。
  • 重複發明的教訓：同一算法五次發明、四次湮沒。科學傳播的斷層（跨領域不讀彼此文獻、語言壁壘、無檢索系統）是真實的歷史成本；今日的 arXiv 與演算法庫是社會性修正。
  • 軍事需求與民用紅利：地震偵測 → FFT → MP3 與 MRI。冷戰科技樹的民用果實，再次印證基礎技術的不可預期外溢——資助基礎研究的理由從來不只是眼前應用。

與其他章節及領域的關聯#

  • 6.1–6.2 節：取樣定理 + DFT + FFT = 數位訊號處理的完整迴路。
  • 4.3 節：頻域濾波的工程可行性由 FFT 決定（重疊相加法）。
  • 7.1–7.2 節：LTI 系統的數值分析（頻率響應測量）靠 FFT。
  • 8.2 節：短時傅立葉轉換 = 重疊窗 FFT；「FFT 的分形複用」。
  • 9.1 節：譜方法解 PDE——FFT 加速使譜方法可與有限差分競爭。
  • 10.2 節：壓縮感知的重建迭代（OMP、ISTA）核心運算是 FFT。
  • 大整數乘法：Schönhage–Strassen（1971）用 FFT 把 \(n\) 位乘法降到 \(O(n\log n\log\log n)\) ——FFT 連乘法本身都加速了。
  • 引力波：LIGO 的匹配濾波以 FFT 為底層運算（每秒百萬次相關運算）。

練習#

  1. 手算 \(N = 8\) 的 radix-2 FFT（輸入 \(x = (1, 1, 0, 0, 0, 0, 0, 0)\) ），並與例題 6.2.1 的 DFT 結果對照。
  2. 證明 radix-2 FFT 的乘法次數上界 \(\frac{N}{2}\log_2 N\) （每級 \(\frac{N}{2}\) 個蝶形、 \(\log_2 N\) 級、每蝶形一次乘法），並計算 \(N = 4096\) 時對直接法的加速比。
  3. 位反轉排序：寫出 \(N = 8\) 時下標 0–7 的位反轉排列，解釋為何「偶奇遞迴拆分」等價於「輸入按位反轉重排」。
  4. 實現（或推導）重疊相加法（overlap-add）：長訊號與 FIR 核的線性捲積用 FFT 分段實現，寫出分段長度 \(L\) 與 FFT 長度 \(M \ge L + K - 1\) 的關係（ \(K\) 為核長），估計總計算量。
  5. 思考題：FFT 的輸出順序（位反轉）與輸入順序（自然）的取捨，「同址計算」省記憶體但打亂順序。查閱「流水線 FFT 晶片」如何用延遲線記憶體（delay-line commutator）解決即時串流的順序問題——演算法與硬體的協同設計如何改變演算法形態？

進階研讀：FFT 的五個深化場景#

場景一：Tukey 的另一半遺產——探索性數據分析。Tukey 同年（1962）發表〈The Future of Data Analysis〉，倡導「探索性數據分析」（EDA）：先看數據再建模型。FFT 是 EDA 的旗艦工具（看頻譜就是探索）；今日的頻譜圖（spectrogram，8.2 節）是 EDA 的標配視覺。「先畫圖後建模」與「先變換後運算」是同一種務實哲學的兩面。

場景二：FFT 的親族——變換域算法家族。數論變換（NTT，用質數模運算避免浮點誤差，密碼學/大整數乘法）、快速 Walsh–Hadamard（方形波基， \(O(N\log N)\) ，壓縮感測測量矩陣）、快速餘弦/正弦變換（JPEG/MP3 底層）、快速 Hartley 變換、快速多項式變換——「分治 + 對稱群結構」的模板可以套用到任何具有循環群對稱的基。FFT 的深層身份：循環群 \(\mathbb{Z}/N\) 表示的快速計算；推廣到非交換群得「快速群傅立葉」（Clifford 代數、圖形訊號處理）。

場景三：離散與連續的計算橋樑誤差。用 DFT 近似連續 FT 的誤差源：混疊（取樣太疏）+ 洩漏（截斷太短）+ 柵欄（箱距太粗）。「Nyquist 準則 + 窗長準則 + 補零準則」三重檢查是數值譜分析的 SOP。9.1 節譜方法會反方向使用這座橋（譜收斂定理）——同一座橋，分析方向相反。

場景四：並行與量子。FFT 的蝶形圖是天然的並行拓撲（蝴蝶網路）：GPU 上 FFT 達到記憶體頻寬上限；量子傅立葉變換（QFT） 用 \(O(n^2)\) 個閘完成 \(2^n\) 點「變換」（指數加速！），是 Shor 質因數分解算法（1994）的核心——但 QFT 只能在量子疊加態上操作，輸入輸出的讀寫仍受測量限制。QFT 是「FFT 再加速」的量子續章，也提醒我們：複雜度的聖杯永遠伴隨新的 I/O 瓶頸。

場景五：FFT 的社會學。1965 年至今，FFT 相關論文數十萬篇、專利數千件、每台手機每天執行億次 FFT。一個 5 頁的算法論文成為百年基礎設施——科學傳播史家稱 FFT 為「20 世紀最有影響力的算法」之一（與快速排序、PageRank 並列）。它的啟示是雙向的：理論數學（群表示）可以突然變成產業；產業需求（軍備控制）可以突然激活沉睡理論。數學與應用的相互拯救，是本書第 10 章的主題預告。

補充：Cooley–Tukey 論文的開頭與結尾#

論文開頭一句話就定位問題：「A finite Fourier series method of complex series and their practical computation…」——它沒有宣稱發明新變換，只宣稱「機器計算的算法」。結尾感謝了 Garwin。論文的自謙與時代的轟動形成反差：作者以為在寫一篇工程筆記，歷史把它讀成里程碑。同時代的 Danielson–Lanczos（1942）與 Runge（1903）無人引用的命運，襯托出「正確時間 + 正確渠道」的傳播學意義——科學史不是思想史的簡單疊加，而是思想 × 載體的卷積。

名詞速查表（本節出現的核心術語）#

  • FFT：DFT 的 \(O(N\log N)\) 算法總稱。
  • 分治 / 偶奇分解：radix-2 的核心遞迴。
  • 蝶形運算： \(a \pm wb\) ；FFT 的基本單元。
  • 旋轉因子 \(w^k = e^{-2\pi ik/N}\) ：蝶形的複數乘數。
  • 位反轉排序：遞迴拆分的等價輸入重排。
  • 同址計算： \(O(N)\) 記憶體的原地算法。
  • 重疊相加 / 重疊儲存：長訊號分塊捲積的兩法。
  • Rader / Bluestein / Winograd：非 2 冪長度的補救算法。
  • QFT：量子傅立葉變換；Shor 算法核心。
  • NTT：數論變換（模質數的精確 FFT）。
  • FFTW：自適應 FFT 庫（1998，Gordon Bell 獎）。

年表速查（本節涉及的關鍵年份）#

  • 1805：Gauss 手稿（失傳的 FFT）。
  • 1903：Runge 的快速方法。
  • 1942：Danielson–Lanczos（晶體學中的分治 DFT）。
  • 1958：Goertzel 算法。
  • 1963：部分禁試條約——地震偵測需求爆發。
  • 1965：Cooley–Tukey 論文發表（4 月 15 日收稿）。
  • 1967–69：FFT 傳播（IEEE 論文、公開程式、應用爆發）。
  • 1971：Schönhage–Strassen 大整數 FFT 乘法。
  • 1976–78：Winograd 與 Rader 系列算法。
  • 1994：Shor 算法（QFT）。
  • 1998：FFTW 自適應庫。
  • 1984：Gauss 優先權考古確認。

深化問答三則#

問：FFT 沒有「新數學」，為什麼被視為偉大成就？ 答：因為它的價值在複雜度質變而非概念新穎—— \(O(N^2)\to O(N\log N)\) 改變了「什麼問題可算」。科學的偉大成就分兩種：新概念（分配理論）與新能力（FFT）。後者對文明的直接槓桿往往更大；本書的立場是兩者等價地值得敬禮。

問：為什麼 \(N\log N\) 是「極限」了嗎？ 答：對一般輸入的 DFT，在標準計算模型下 \(\Omega(N\log N)\) 是普遍相信的下界（未嚴格證明！這是理論計算機科學的懸案之一）。但特定結構的輸入可以更快：稀疏譜（6.2 節、10.2 節壓縮感知的 sparse FFT， \(O(s\log N)\) ， \(s\) 為非零譜線數）與整數數據有亞線性算法——下界屬於最壞情況，好數據永遠有可乘之機。

問：FFT 對「傅立葉分析」的理論地位有何改變？ 答：微妙而深刻：FFT 使「變換到頻域」從概念操作變成預設動作——工程師不再問「要不要變換」而問「用哪個窗、多長段」。理論上，它使 \(L^2\) 理論、時頻分析（第 8 章）、譜方法（第 9 章）的全部構想獲得計算載體。理論提出問題，算法使問題成為日常——FFT 是「理論—算法—文化」三級火箭的教科書案例。

（本節完）