3.2 快速傅立葉轉換FFT#

概念說明#

3.1 節的 DFT 定義要 \(O(N^2)\) 次運算，1965 年 Cooley–Tukey 提出的快速傅立葉轉換（FFT）把它降到 \(O(N \log N)\) ，被譽為二十世紀最重要的演算法之一。沒有 FFT，就沒有即時頻譜分析、MP3、語音辨識與現代通訊。

FFT 的關鍵觀察是分治：把 \(N\) 點 DFT 拆成兩個 \(N/2\) 點 DFT——一個只含偶數下標樣本，一個只含奇數下標樣本。利用旋轉因子 \(W_N = e^{-j2\pi/N}\) 的對稱性（ \(W_N^{k+N/2} = -W_N^k\) ），兩個半長 DFT 的結果可以組合成完整 DFT：

其中 \(E[k]\) 、 \(O[k]\) 是偶、奇子序列的 DFT。這個「蝶形運算（butterfly）」每層把問題減半，共 \(\log_2 N\) 層，每層 \(N\) 次乘加，總計 \(O(N \log N)\) 。

以 \(N = 4096\) 為例：直接 DFT 約需 1677 萬次乘加，FFT 約 4.9 萬次——快 340 倍。 \(N\) 越大差距越懸殊。

兩個實務重點：

1. 長度：基數-2（radix-2）FFT 要求 \(N = 2^m\) 。長度不為 2 的冪時，可補零（不影響頻譜形狀，只增加內插的頻點）或用混合基演算法。numpy 的 np.fft 兩者都支援。
2. 位址反轉：遞迴拆分會把輸入按二進位位址反轉重排，迭代實作先做位址反轉置換，再逐層做蝶形，即可原地計算、只需 \(O(1)\) 額外空間（numpy 內部已處理）。

演算法步驟：迭代式基數-2 FFT#

1. 輸入長度 \(N = 2^m\) 的複數序列。
2. 將輸入依位址反轉順序重排。
3. 對 \(s = 1, 2, \dots, \log_2 N\) （子段長 \(L = 2^s\) ）：
4. 對每個長度 \(L\) 的子段，取半段 \(L/2\) ；
5. 對每個 \(k \in [0, L/2)\) 計算旋轉因子 \(W = e^{-j2\pi k/L}\) ；
6. 執行蝶形： \(u = a[k]\) 、 \(v = W \cdot a[k+L/2]\) ； \(a[k] = u + v\) 、 \(a[k+L/2] = u - v\) 。
7. 輸出即為 DFT（下標依自然順序）。

Python 範例#

import numpy as np

def fft_radix2(x):
    x = np.asarray(x, dtype=complex)
    N = x.size
    assert N & (N - 1) == 0, "N 必須是 2 的冪"
    # 位址反轉重排
    idx = np.arange(N)
    rev = np.zeros_like(idx)
    bits = N.bit_length() - 1
    for b in range(bits):
        rev |= ((idx >> b) & 1) << (bits - 1 - b)
    a = x[rev].copy()
    # 逐層蝶形
    L = 2
    while L <= N:
        half = L // 2
        W = np.exp(-2j * np.pi * np.arange(half) / L)
        for start in range(0, N, L):
            u = a[start:start+half]
            v = a[start+half:start+L] * W
            a[start:start+half] = u + v
            a[start+half:start+L] = u - v
        L *= 2
    return a

# 正確性驗證：與 numpy 內建 FFT 比較
rng = np.random.default_rng(0)
for N in [8, 256, 1024]:
    x = rng.normal(size=N)
    err = np.max(np.abs(fft_radix2(x) - np.fft.fft(x)))
    print(f"N={N:5d}  與 np.fft 的最大誤差 = {err:.2e}")

# 效能比較
import time
N = 2**16
x = rng.normal(size=N)
t0 = time.perf_counter(); np.fft.fft(x); t_fft = time.perf_counter() - t0
print(f"\nN={N}: np.fft.fft 耗時 {t_fft*1000:.2f} ms（直接 DFT 需 ~{N**2/1e6:.0f} 百萬次乘加）")

輸出：

N=    8  與 np.fft 的最大誤差 = 1.78e-14
N=  256  與 np.fft 的最大誤差 = 2.84e-13
N= 1024  與 np.fft 的最大誤差 = 9.66e-13

N=65536: np.fft.fft 耗時 1.81 ms（直接 DFT 需 ~4295 百萬次乘加）

自製 FFT 與 numpy 內建結果一致（浮點誤差在 \(10^{-13}\) 量級），驗證了蝶形結構的正確性。實務上永遠使用 np.fft（以 C 實作、經過高度最佳化），自製版本用於理解原理。

複雜度分析#

時間： \(T(N) = 2T(N/2) + O(N)\) ，由主定理得 \(O(N \log N)\) 。乘法次數約 \(\frac{N}{2}\log_2 N\) ，加法約 \(N\log_2 N\) 。空間：原地迭代實作 \(O(1)\) 額外（不含輸出）；numpy 實作為 \(O(N)\) 。對比： \(N = 2^{20}\) （約 24 秒的 44.1 kHz 音訊）時， \(O(N^2)\) 約 \(10^{12}\) 次運算（數分鐘），FFT 約 \(2\times 10^7\) 次（數毫秒）——相差五個數量級。這就是為何頻譜圖（3.3 節）、卷積加速、以及第 5 章以後的特徵提取全部建構在 FFT 之上。

練習#

1. 將 fft_radix2 改為遞迴版本，並繪製（或文字描述） \(N=8\) 時蝶形運算的資料流圖，標出偶、奇子序列。
2. 量測 fft_radix2 與 np.fft.fft 在 \(N = 2^{10} \dots 2^{16}\) 的執行時間，驗證成長趨勢符合 \(N\log N\) （以雙對數座標的斜率判斷）。
3. 利用「循環卷積 = IFFT(FFT(a)·FFT(b))」實作長度 1024 的兩訊號卷積，與 np.convolve 的直接法比較結果與時間；說明卷積定理如何加速第 4 章的 FIR 濾波。