第三節 - DFT 的加速運算（FFT）

講解完 DFT 的常見性質之後，這節要來介紹它的快速演算法，後面那三個小節的應用基本上都是因為有了這些快速演算法才能凸顯出它們的意義。FFT 的全名為 fast Fourier transform，它其實只是一個統稱，這世界上有太多太多能加速 DFT 的演算法 [3] 了，只要是企圖達成這個目標的方法都可以算是一種 FFT。

一般來說最知名的 FFT 是 Cooley 和 Tukey 在 1965 年所提出的以「因數分解」為核心思想的方法，[4] 有簡單俐落的推導和時間複雜度分析 (可以先弄懂這一節最後的範例之後再來看這篇)，而 [5] 則是當年發表的 paper，只是大部分剛入門 FFT 的初心者一開始會學習的僅是該算法的其中一種特例。

若土法煉鋼計算 DFT 頻譜，按照原本的公式，一條頻譜需要$N$點時域訊號，全部$N$條頻譜計算完畢所須時間複雜度為$\mathcal{O}(N^2)$。我們會先介紹上一段所說的特例，在$N$是 2 的冪次方的條件下，如果好好利用旋轉因子 (twiddle factor) 在複數平面上的單位圓旋轉的特性，便能發現某些運算結果其實可以重複使用，進而透過分治法 (divide and conquer) 策略將整個時間的複雜度降至 $\mathcal O(N\log N)$；理解這個特例之後再進一步觀察怎麼將點數按照原本的模式推廣，最後竟然會發現又多了一層 DFT，才得到因數分解的結論。初學者可以僅先理解第一段的特例即可，第二段的推廣版本不用急著弄懂，有空再回來補即可。

# Radix-2 Decimation in Time

先回顧一下 DFT 的公式，如果我們先定義旋轉因子 $W_N = e^{j2\pi(\frac{-1}{N})}$且有符號慣例 $W_N^{kn} = (W_N)^{kn}$，那麼： 　　　　　　　$\displaystyle X[k] = \sum _{n=0}^{N-1}x[n]\ e^{-j2\pi{k\over N}n} = \sum _{n=0}^{N-1}x[n]\ \Big(e^{j{2\pi}(\frac{-1}N)}\Big)^{kn} = \sum _{n=0}^{N-1}x[n]\ W_N^{kn}$

夠細心的讀者應會發現$W_N$其實就是在複數平面單位圓上切 $N$等分、從$0$度開始往順時針方向取第一格的值，如下圖所示。

有了旋轉因子的概念之後，便可以把 DFT 公式重新改寫如下： \displaystyle \begin{alignat*}{2}X[k] &= \sum_{n=0}^{N-1}x[n]\ W_N^{kn} \\&= \sum_{m=0}^{\frac N2-1}x[2m]\ W_N^{k(2m)} + \sum_{m=0}^{\frac N2-1}x[2m+1]\ W_N^{k(2m+1)} \quad\quad//\ 分成偶數時點與奇數時點\\&=\sum_{m=0}^{\frac N2-1}x[2m]\ W_N^{k(2m)} + W_N^k\sum_{m=0}^{\frac N2-1}x[2m+1]\ W_N^{k(2m)} \quad\quad//\ 把奇數邊的一個旋轉因子提出來 \\&=\sum_{m=0}^{\frac N2-1}x_e[m]\ W_N^{k(2m)} + W_N^k\sum_{m=0}^{\frac N2-1}x_o[m]\ W_N^{k(2m)}\quad\quad//\ 把偶訊號與奇訊號重新命名 \\&=\sum_{m=0}^{\frac N2-1}x_e[m]\ W_{\frac N2}^{km} + W_N^k\sum_{m=0}^{\frac N2-1}x_o[m]\ W_{\frac N2}^{km}\quad\quad//\ N\ 等分的第\ 2km\ 個點就相當於\ \frac N2\ 等分的第\ km\ 個點\ (很直觀的想法) \\&=X_e[k] + W_N^k\cdot X_o[k]\quad\quad\ //\ 把偶訊號和奇訊號分別寫成各自的\text{ DFT }形式\end{alignat*} 到這裡，可以發現到原本的 DFT 竟然能拆成偶數時點的 DFT 再加上奇數時點的 DFT 乘上一個旋轉因子的$k$次方。那麼要如何利用這個「旋轉」的特性很快的去算出另外一條頻譜的值呢？讀者應可以理解到這個式子主要的目的在於把問題的大小 (problem size) 減半，原先是 $N$個點的難度被拆解成兩組 $N/2$ 個點的難度，那麼自然 $X_e$和 $X_o$ 的週期就都會是 $N/2$ 個點，而且和奇訊號 DFT 相乘的旋轉因子也會因為多轉半圈的關係而多了個負號：

　　　　　　　　　　　　　　　　　　$X_e[k+\dfrac N2] = X_e[k]$　　　　　　　　　　　　　　　　　　

　　　　　　　　　　　　　　　　　　$X_o[k+\dfrac N2] = X_o[k]$

　　　　　　　　　　　　　　　　　　　　 $W_N^{k+\frac N2} = -W_N^k$

有了這三個關鍵條件便足以說明 $X[k]$ 和 $X[k+\dfrac N2]$ 幾乎可以同時算出。

　　　　　　$X[k] 　　　　　　　　　　　　　　　　　　=\ \ \ X_e[k] + W_N^{k}\cdot X_o[k]$

　　　$X[k+\dfrac N2]\ \ \ =\ \ \ X_e[k+\dfrac N2] + W_N^{k+\frac N2}\cdot X_o[k+\dfrac N2]\ \ \ =\ \ \ X_e[k] - W_N^{k}\cdot X_o[k]$　　　

因此時間複雜度 $T(n) = 2T(n/2) + \mathcal{O}(n) = \mathcal{O}(n\log n)$，確實比 $\mathcal O(n^2)$還要快的！我們這邊不會解釋如何推導複雜度，有興趣的讀者可自行修習演算法相關課程。

下圖是展示如何用分治法的策略去拿「計算 4 點 DFT」的這個小模組來計算 8 點 DFT 的架構圖。

MATLAB 遞迴版本程式碼如下：

function Xdft = fftrecur(x)
    N = length(x); % N should be a power of 2
    if N == 1
        Xdft = x;
    else
        m = N/2;
        XE = fftrecur(x(1:2:N));
        XO = fftrecur(x(2:2:N));
        W = exp(-2*pi*sqrt(-1)/N).^(0:m-1).';
        temp = W.*XO;
        Xdft = [ XE+temp ; XE-temp ];
    end
end

# General Version (DIT + DIF)

上一段介紹的是針對「$N$是 2 的冪次方」的情況下以 2 為基底的時間抽取法，那如果$N$不是 2 的冪次方又該怎麼辦呢？我們這邊以 $N=15=5\times3$ 作為例子，把上一段以除以 2 得到的餘數為依據的分類法推廣到除以 3 的情況，看看會有什麼有趣的事情發生。 \displaystyle \begin{alignat*}{2}X[k] &= \sum_{n=0}^{15-1}x[n]\ W_{15}^{kn}\quad\quad\ //\ 令\ N=5\times3 \\&= \sum_{m=0}^{5-1}x[3m]\ W_{15}^{k(3m)} + \sum_{m=0}^{5-1}x[3m+1]\ W_{15}^{k(3m+1)} + \sum_{m=0}^{5-1}x[3m+2]\ W_{15}^{k(3m+2)} \quad\quad//\ 按照\ N\ 除以\ 3\ 的餘數分成三堆\\&=\sum_{m=0}^{5-1}x[3m]\ W_{15}^{k(3m)} + W_{15}^k\sum_{m=0}^{5-1}x[3m+1]\ W_{15}^{k(3m)} + W_{15}^{2k}\sum_{m=0}^{5-1}x[3m+2]\ W_{15}^{k(3m)} \quad\quad//\ 按照餘數提出適當的旋轉因子 \\&=\sum_{m=0}^{5-1}x_0[m]\ W_{15}^{k(3m)} + W_{15}^k\sum_{m=0}^{5-1}x_1[m]\ W_{15}^{k(3m)} + W_{15}^{2k}\sum_{m=0}^{5-1}x_2[m]\ W_{15}^{k(3m)}\quad\quad//\ 把分類好的訊號重新命名 \\&=\sum_{m=0}^{5-1}x_0[m]\ W_{5}^{km} + W_{15}^k\sum_{m=0}^{5-1}x_1[m]\ W_{5}^{km} + W_{15}^{2k}\sum_{m=0}^{5-1}x_2[m]\ W_{5}^{km}\quad\quad//\ 15\ 等分的第\ 3km\ 個點就相當於\ 5\ 等分的第\ km\ 個點\ (很直觀的想法) \\&=\text{DFT}(x_0)[k] + W_{15}^k\cdot\text{DFT}(x_1)[k] + W_{15}^{2k}\cdot\text{DFT}(x_2)[k]\quad\quad\ //\ 把分類好的訊號分別寫成各自的\text{ DFT }形式 \\&=\text{DFT}(x_0)[5q+r] + W_{15}^{5q+r}\cdot\text{DFT}(x_1)[5q+r] + W_{15}^{2(5q+r)}\cdot\text{DFT}(x_2)[5q+r]\quad\quad\ //\ 把\ k\ 除以\ 5\ 得到商數\ q\ 和餘數\ r \\&=\text{DFT}(x_0)[r] + W_{3}^{q}\Big(W_{15}^r\cdot\text{DFT}(x_1)[r]\Big) + W_{3}^{2q}\Big(W_{15}^{2r}\cdot\text{DFT}(x_2)[r]\Big)\quad\quad\ //\ 利用圓的週期性化簡 \\&=\sum_{i=0}^{3-1}\Big(W_{N}^{ir}\cdot\text{DFT}(x_i)[r]\Big)\cdot W_{3}^{qi}\quad\quad\ //\ 再度表成以小的\text{ DFT }為訊號的大\text{ DFT }形式\end{alignat*} 你會發現一開始都和上一段差不多，到了後面寫成各自的 DFT 形式那邊就開始不太一樣了，旋轉因子會有多的餘數要處理，最後又會發現另一層 DFT 的存在，這就是它奇妙的地方。

歸納一下： $\displaystyle X[k]=X[5q+r]=X_r[q]=\sum_{i=0}^{3-1}\Big(W_{N}^{ir}\cdot\text{DFT}(x_i)[r]\Big)\cdot W_{3}^{qi}$。

下圖是上面 $N=5\times3$ 的推導過程的架構化，應該還蠻清楚的，讀者可以自行撰寫 MATLAB 程式碼驗證這個架構是正確的。

(如果覺得太小的話可以直接用滑鼠對圖片按一下，會自動放大)