JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

庫利-圖基：從 O(N^2) 到 O(N log N)

第一篇老老實實地建立了離散傅立葉轉換——也讓我們看到它要花 O(N^2)。現在我們要看二十世紀最偉大的演算法之一，怎麼把問題對半切、再對半切，直到那道平方的高牆崩塌成 O(N log N)。這招比電腦還老，就藏在一個代數的拆分裡。

為什麼 O(N^2) 是一道牆，不只是個小麻煩

上一篇從頭建立了離散傅立葉轉換，也把它的成本攤開來看：要得到全部 N 個輸出頻率，你得算 N 個總和，而每個總和都碰到全部 N 個輸入，所以你大約做 N^2 次複數乘加。對一個玩具訊號 N = 1000 來說那是一百萬次運算——一瞬間就完了。但傅立葉轉換是在大尺度下生活的：一張一百萬像素影像的一條邊，N 就接近 10^6，這時 N^2 是 10^12，單單一次轉換就要一兆次運算。以每秒十億次運算來算，那是一刻鐘，只為了一張影像。音訊、雷達、磁振造影、天氣模型都要呼叫轉換好幾百萬次。在 O(N^2) 之下，整件事根本就是不可能的。

所以問題不是「我們能不能省下一個常數倍率？」——那沒有用。我們得改變成長速率本身，也就是 N 上的那個指數。這正是計算複雜度的核心：一個 O(N log N) 的演算法不是比 O(N^2) 快一點點，而是完全不同的世界。在 N = 10^6 時，N^2 是 10^12，但 N log N（以 2 為底）大約是 2 乘以 10^7——大約五萬倍的加速。那個一兆次運算的轉換變成兩千萬次運算，從幾分鐘變成幾毫秒。找到這樣一個演算法，正是讓快速傅立葉轉換成為人類寫進程式碼裡最具影響力的數學之一的原因。

那個搞定一切的代數拆分

整個祕密就在這裡，而且就是一個念頭。假設 N 是偶數。把你的 N 個輸入分成位在偶數位置的那些（x_0、x_2、x_4、……）和位在奇數位置的那些（x_1、x_3、x_5、……）。每一半有 N/2 個點。驚人的事實是：整個訊號那個大小為 N 的 DFT，可以免費地從偶數那半的大小 N/2 的 DFT、加上奇數那半的大小 N/2 的 DFT 重新拼出來。當你把 DFT 的總和按偶數、奇數索引項分開，再把一個共同的相位提出來，它就恰好重組成兩個較小的轉換黏在一起。我們把一個大小 N 的問題，變成了兩個大小 N/2 的問題——這正是分而治之的招牌動作。

Let w = exp(-2*pi*i / N)              # the basic twiddle (Nth root of unity)
For each output index k = 0 .. N/2 - 1 :
    E_k = DFT of the even-indexed samples, at frequency k   # size N/2
    O_k = DFT of the odd-indexed samples,  at frequency k   # size N/2
    X_k         = E_k +  w^k * O_k        # one butterfly...
    X_{k + N/2} = E_k -  w^k * O_k        # ...gives TWO outputs
庫利-圖基的重組：每一對輸出 k 與 k+N/2 都來自同一個 E_k 與 O_k，差別只在旋轉因子 w^k 的正負號。一次乘法，兩個輸出。

仔細看看為什麼一次乘法能生出兩個輸出。輸出 X_k 用的是 E_k + w^k O_k，而隔了半圈的輸出 X_{k+N/2} 用的是 E_k - w^k O_k——同樣那兩個數、同樣那個旋轉因子，只是翻了個正負號。這個小小的一加一減配對就是蝴蝶運算，因為畫出資料流向時那交叉的雙翼形狀而得名。由於每個蝴蝶運算用一次複數乘法產出兩個輸出，把兩半重組的成本只大約是 N/2 次乘法，而不是 N^2。昂貴的部分已經被外包給那兩個較小的轉換了。

一路遞迴到底：N log N 在哪裡誕生

我們還沒做完——我們只是切了一半。但同樣的拆分對每個 N/2 的轉換都適用：把它切成自己的偶數與奇數，做出兩個大小 N/4 的轉換，依此類推。如果 N 是 2 的次方，你就能一直切半，直到抵達大小為 1 的轉換，而單一個數的 DFT 就是那個數本身——沒有東西要算。於是遞迴乾淨地觸底。按層數來數工作量：每一層切半要花大約 N 次乘加去重組，而總共有 log_2(N) 層，因為那正是你能把 N 一路切半到 1 的次數。總工作量：大約 N 乘以 log_2(N)，也就是那赫赫有名的 O(N log N)。

  1. 基本情形：若 N = 1，單一樣本的轉換就是那個樣本。直接回傳。
  2. 把 N 個樣本分成偶數索引那一半與奇數索引那一半，各為大小 N/2。
  3. 對每一半遞迴地做 FFT——這就是演算法在更小的問題上呼叫它自己。
  4. 合併：對 k = 0 .. N/2 - 1，做蝴蝶運算 X_k = E_k + w^k O_k 與 X_{k+N/2} = E_k - w^k O_k。
  5. 回傳這 N 個合併後的輸出。重活是遞迴做的；這一層只負責黏合。

這就是庫利-圖基演算法（庫利-圖基 FFT），發表於 1965 年——不過高斯早在大約 1805 年就在手算裡悄悄用過同樣的切半念頭，那時候連傅立葉自己的著作都還沒出版。上面這個版本是 radix-2（基 2）形式，最容易理解。它跟第一篇的 DFT 是完全一樣的——同樣的輸出、同樣的數學——只是用更聰明的順序去求值。快速傅立葉轉換不是 DFT 的近似；它是一模一樣的答案，只是重新組織過，讓共用的子計算只做一次，而不是做 N 次。

在真正的機器上，到底是什麼讓它快

浮點運算次數是頭條，但它不是故事的全部，而誠實地交代很重要。把資料重排成偶數與奇數、在每一層重複這件事，最後會把輸入排列成位元反轉的順序——索引 6（二進位 110）和索引 3（二進位 011）對調，依此類推。實用的 FFT 會就地做這個排列，然後完全不用遞迴、掃過 log_2(N) 趟蝴蝶運算。遞迴是理解它的乾淨方式；迴圈是執行它的快速方式。旋轉因子 w^k 通常會預先算進一張小表裡，因為在內層迴圈裡重算正弦與餘弦會反過來主宰成本。

那反向轉換呢？反離散傅立葉轉換是同樣形狀的總和，只是指數的正負號翻過來、再加一個 1/N 的縮放，所以同樣那批蝴蝶運算、配上共軛的旋轉因子，就能以同樣的 O(N log N) 成本跑它。一個演算法、正著跑與反著跑，就給了你轉換的兩個方向——而這恰好是接下來兩篇要利用的，去把卷積變快、去用驚人的準確度對函數做微分。

誠實的但書：它很快，但不是魔法，也不是精確的

兩個誠實的限定能讓你不至於過度信任 FFT。第一，乾淨的 radix-2 故事需要 N 是 2 的次方。真實訊號沒那麼乖——你可能拿到 N = 1000，或一個質數。函式庫用混合基切分來應付（任何因式分解 N = p 乘以 q 都同樣行得通），而對於難搞的 N，則用 Bluestein 的把戲把轉換重新表達成一個卷積。結論是：一個好的 FFT 對幾乎每個 N 都跑在 O(N log N)，但質數長度會比最接近的 2 的次方慢，這就是為什麼人們在能補零時常常把長度補到一個方便的大小。

第二——而這正是整條學習階梯一再教我們的一課——FFT 是在浮點算術裡計算的，所以它的輸出是近似的，不是精確的 DFT。不過這裡有個真正的好消息：FFT 比那個天真的 O(N^2) 總和更準，而不是更不準。直接求和會在 N 次加法中累積捨入誤差，所以它的誤差像 O(N) 乘以單位捨入那樣成長。FFT 只在一棵深度 log_2(N) 的平衡樹裡相加，所以它的捨入誤差像 O(log N) 那樣成長——讓它省時間的同一招分而治之，也限制了誤差的累積。更少、組織得更好的運算，意味著又快又準。這是一場難得而美麗的雙贏。