快速傅立葉轉換

本文的閱讀等級：中級

給定一序列 $x_0,x_1,\ldots,x_{n-1}$ ，離散傅立葉轉換的計算公式為 (見“離散傅立葉轉換”)

$\displaystyle y_k=\sum_{j=0}^{n-1}x_je^{-2\pi ikj/n},~~~k=0,1,\ldots,n-1$ 。

令 $\displaystyle w=e^{-2\pi i/n}=\cos(2\pi/n)-i\sin(2\pi/n)$ 。離散傅立葉轉換可表示成矩陣形式 $\mathbf{y}=F\mathbf{x}$ ，如下：

$\begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{n-1} \end{bmatrix}=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&w&w^2&\cdots&w^{n-1}\\ 1&w^2&w^4&\cdots&w^{2(n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&w^{n-1}&w^{2(n-1)}&\cdots&w^{(n-1)^2} \end{bmatrix}\begin{bmatrix} x_0\\ x_1\\ x_2\\ \vdots\\ x_{n-1} \end{bmatrix}$ ，

其中 $n\times n$ 階 $F$ 稱為傅立葉矩陣。若採用一般矩陣乘法運算，離散傅立葉轉換的計算複雜度為 $\mathcal{O}(n^2)$ 。公元1965年，美國學者庫利 (James Cooley) 與圖基 (John Tukey) 提出了一個複雜度為 $\mathcal{O}(n\log_2n)$ 的演算法，稱為快速傅立葉轉換 (fast Fourier transform，簡稱 FFT)，後來人們發現原來高斯 (Carl Friedrich Gauss) 早在1805年就已經提出同樣的演算法。本文採用別於傳統訊號處理的推導方法^[1]，我們先檢視傅立葉矩陣 $F$ 擁有的特殊性質，然後運用分塊矩陣技巧推導 Cooley-Tukey 演算法。除了本文介紹的 Cooley-Tukey 演算法，還有其他多種形式的快速傅立葉轉換演算法，請讀者參閱維基百科^[2]。

Cooley-Tukey 演算法的主要構想是將序列長為 $n$ 的離散傅立葉轉換分割成兩個長為 $n/2$ 的子序列的離散傅立葉轉換。為簡化分析，以下僅考慮 $n=2^r$ 的情況。令 $F_n=[f_{pq}]$ 表示 $n\times n$ 階傅立葉矩陣，列行指標設定為 $p,q\in\{0,1,\ldots,n-1\}$ ，其中 $f_{pq}=w^{pq}$ ，組成單元 $w=e^{-2\pi i/n}$ 具有下列性質。

(1) 複數集 $\{1,w,w^2,\ldots,w^{n-1}\}$ 為方程式 $z^n=1$ 的根，稱為 $n$ 次方單位根，這些根在複數平面上單位圓按順時針旋轉並以相同間距排列 (下圖是 $n=8$ 的例子)。此外， $\{1,w^2,w^4,\ldots,w^{n-2}\}$ 是 $z^{n/2}=1$ 的 $n/2$ 次方單位根。同樣道理，對於 $m=2^k$ ， $k$ 是正整數， $\{1,w^m,w^{2m},\ldots,w^{n-m}\}$ 是 $z^{n/m}=1$ 的 $n/m$ 次方單位根，快速傅立葉轉換即建立於此性質之上。

(2) $w$ 具有週期性。若 $k\ge n$ ，則 $w^k=w^{k\!\mod n}$ ，其中 $k\!\mod n$ 代表 $k$ 除以 $n$ 的餘數。譬如，若 $n=8$ ，則 $w^8=1, w^9=w, w^{10}=w^2,\ldots$ 。

(3) $w$ 具有對稱性，即 $w^{k+\frac{n}{2}}=-w^k$ ，因為 $w^p=e^{2\pi ip/n}=\cos(2\pi\frac{p}{n})+i\sin(2\pi\frac{p}{n})$ ，且

$\displaystyle \cos\left(2\pi\frac{k+\frac{n}{2}}{n}\right)=\cos\left(2\pi\frac{k}{n}+\pi\right)=-\cos\left(2\pi\frac{k}{n}\right)$ ，

$\displaystyle \sin\left(2\pi\frac{k+\frac{n}{2}}{n}\right)=\sin\left(2\pi\frac{k}{n}+\pi\right)=-\sin\left(2\pi\frac{k}{n}\right)$ 。

若 $n=8$ ，就有 $w^4=-1, w^5=-w, w^6=-w^2, w^7=-w^3$ 。

8th roots of unity

為便利說明，我們以 $n=8$ 為例推導 Cooley-Tukey 演算法。利用性質 (2) 和 (3)， $8\times 8$ 階傅立葉矩陣可化簡為

$\begin{aligned} F_8&=\begin{bmatrix} 1&1&1&1&1&1&1&1\\ 1&\omega&\omega^2&\omega^3&\omega^4&\omega^5&\omega^6&\omega^7\\ 1&\omega^2&\omega^4&\omega^6&1&\omega^2&\omega^4&\omega^6\\ 1&\omega^3&\omega^6&\omega&\omega^4&\omega^7&\omega^2&\omega^5\\ 1&\omega^4&1&\omega^4&1&\omega^4&1&\omega^4\\ 1&\omega^5&\omega^2&\omega^7&\omega^4&\omega&\omega^6&\omega^3\\ 1&\omega^6&\omega^4&\omega^2&1&\omega^6&\omega^4&\omega^2\\ 1&\omega^7&\omega^6&\omega^5&\omega^4&\omega^3&\omega^2&\omega \end{bmatrix}\\ &=\left[\!\!\begin{array}{rrrrrrrr} 1&1&1&1&1&1&1&1\\ 1&\omega&\omega^2&\omega^3&-1&-\omega&-\omega^2&-\omega^3\\ 1&\omega^2&-1&-\omega^2&1&\omega^2&-1&-\omega^2\\ 1&\omega^3&-\omega^2&\omega&-1&-\omega^3&\omega^2&-\omega\\ 1&-1&1&-1&1&-1&1&-1\\ 1&-\omega&\omega^2&-\omega^3&-1&\omega&-\omega^2&\omega^3\\ 1&-\omega^2&-1&\omega^2&1&-\omega^2&-1&\omega^2\\ 1&-\omega^3&-\omega^2&-\omega&-1&\omega^3&\omega^2&\omega \end{array}\!\!\right],\end{aligned}$

再利用性質 (1)， $4\times 4$ 階傅立葉矩陣為

$F_4=\begin{bmatrix} 1&1&1&1\\ 1&\omega^2&\omega^4&\omega^6\\ 1&\omega^4&1&\omega^4\\ 1&\omega^6&\omega^4&\omega^2 \end{bmatrix}=\left[\!\!\begin{array}{rrrr} 1&1&1&1\\ 1&\omega^2&-1&-\omega^2\\ 1&-1&1&-1\\ 1&-\omega^2&-1&\omega^2 \end{array}\!\!\right]$ 。

觀察後發現 $F_8$ 的偶數行，即第0,2,4,6行，完全由 $F_4$ 構造而成，於是我們將 $F_8$ 各行重新置換，偶數行排在前，奇數行排在後。設定下列排列矩陣：

$P_8^T=\begin{bmatrix} \mathbf{e}_0&\mathbf{e}_2&\mathbf{e}_4&\mathbf{e}_6&\mathbf{e}_1&\mathbf{e}_3&\mathbf{e}_5&\mathbf{e}_7 \end{bmatrix}=\begin{bmatrix} 1&0&0&0&0&0&0&0\\ 0&0&0&0&1&0&0&0\\ 0&1&0&0&0&0&0&0\\ 0&0&0&0&0&1&0&0\\ 0&0&1&0&0&0&0&0\\ 0&0&0&0&0&0&1&0\\ 0&0&0&1&0&0&0&0\\ 0&0&0&0&0&0&0&1 \end{bmatrix}$ ，

即得

$F_8P_8^T=\left[\!\!\begin{array}{rrrrrrrr} 1&1&1&1&1&1&1&1\\ 1&\omega^2&-1&-\omega^2&\omega&\omega^3&-\omega&-\omega^3\\ 1&-1&1&-1&\omega^2&-\omega^2&\omega^2&-\omega^2\\ 1&-\omega^2&-1&\omega^2&\omega^3&1&-\omega^3&-\omega\\ 1&1&1&1&-1&-1&-1&-1\\ 1&\omega^2&-1&-\omega^2&-\omega&-\omega^3&\omega&\omega^3\\ 1&-1&1&-1&-\omega^2&\omega^2&-\omega^2&\omega^2\\ 1&-\omega^2&-1&\omega^2&-\omega^3&-\omega&\omega^3&\omega \end{array}\!\!\right]=\begin{bmatrix} F_4&D_4F_4\\ F_4&-D_4F_4 \end{bmatrix}$ ，

上式中， $D_4=\mathrm{diag}(1,w,w^2,w^3)$ ，注意 $\{1,w,w^2,w^3\}$ 是 $z^8=1$ 的前四個根。一般情況可歸納如下：同樣令 $w=e^{-2\pi i/n}$ ，對於 $m=2^k\le n$ ， $k$ 是正整數，因為 $P_m^{-1}=P_m^T$ ，故 $m\times m$ 階傅立葉矩陣 $F_m$ 有分塊矩陣表達式

$F_m=\begin{bmatrix} F_{\frac{m}{2}}&D_{\frac{m}{2}}F_{\frac{m}{2}}\\ F_{\frac{m}{2}}&-D_{\frac{m}{2}}F_{\frac{m}{2}} \end{bmatrix}P_{m}$ ，

其中 $F_{\frac{m}{2}}$ 是 $m/2\times m/2$ 階傅立葉矩陣， $D_{\frac{m}{2}}=\mathrm{diag}(1,w^{n/m},w^{2n/m},\ldots,w^{n/2-n/m})$ ，主對角元為 $z^m=1$ 的前 $m/2$ 個根，而 $P_m^T=\begin{bmatrix} \mathbf{e}_0&\mathbf{e}_2&\cdots\mathbf{e}_{m-2}&\mathbf{e}_1&\mathbf{e}_3&\cdots&\mathbf{e}_{m-1} \end{bmatrix}$ 。

下面以 $n=8$ 為例展示快速傅立葉轉換的演算程序。給定待變換向量 $\mathbf{x}_8$ 及其置換結果如下：

$\mathbf{x}_8=\begin{bmatrix} x_0\\ x_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7 \end{bmatrix},~P_8\mathbf{x}_8=\begin{bmatrix} x_0\\ x_2\\ x_4\\ x_6\\ x_1\\ x_3\\ x_5\\ x_7 \end{bmatrix}=\begin{bmatrix} \mathbf{x}_4^{(0)}\\ --\\ \mathbf{x}_4^{(1)} \end{bmatrix}$ ，

其中 $\mathbf{x}_4^{(0)}$ 和 $\mathbf{x}_4^{(1)}$ 分別由 $\mathbf{x}_8$ 的偶數元和奇數元組成。運用分塊矩陣乘法運算，可得

$F_8\mathbf{x}_8=\begin{bmatrix} F_4&D_4F_4\\ F_4&-D_4F_4 \end{bmatrix}P_8\mathbf{x}_8=\begin{bmatrix} F_4&D_4F_4\\ F_4&-D_4F_4 \end{bmatrix}\begin{bmatrix} \mathbf{x}_4^{(0)}\\ \mathbf{x}_4^{(1)} \end{bmatrix}=\begin{bmatrix} F_4\mathbf{x}_4^{(0)}+D_4F_4\mathbf{x}_4^{(1)}\\ F_4\mathbf{x}_4^{(0)}-D_4F_4\mathbf{x}_4^{(1)} \end{bmatrix}$ 。

上式清楚顯示快速傅立葉轉換的分治 (divide-and-conquer) 過程：將8階乘法 $F_8\mathbf{x}_8$ 替換為兩個4階乘法 $F_4\mathbf{x}_4^{(0)}$ 和 $F_4\mathbf{x}_4^{(1)}$ 。使用遞歸策略繼續分解，令

$P_4\mathbf{x}_4^{(0)}=\begin{bmatrix} x_0\\ x_4\\ x_2\\ x_6\end{bmatrix}=\begin{bmatrix} \mathbf{x}_2^{(0)}\\ --\\ \mathbf{x}_2^{(1)} \end{bmatrix},~~P_4\mathbf{x}_4^{(1)}=\begin{bmatrix} x_1\\ x_5\\ x_3\\ x_7\end{bmatrix}=\begin{bmatrix} \mathbf{x}_2^{(2)}\\ --\\ \mathbf{x}_2^{(3)} \end{bmatrix}$ ，

就有

$F_4\mathbf{x}_4^{(0)}=\begin{bmatrix} F_2&D_2F_2\\ F_2&-D_2F_2 \end{bmatrix}P_4\mathbf{x}_4^{(0)}=\begin{bmatrix} F_2&D_2F_2\\ F_2&-D_2F_2 \end{bmatrix}\begin{bmatrix} \mathbf{x}_2^{(0)}\\ \mathbf{x}_2^{(1)} \end{bmatrix}=\begin{bmatrix} F_2\mathbf{x}_2^{(0)}+D_2F_2\mathbf{x}_2^{(1)}\\ F_2\mathbf{x}_2^{(0)}-D_2F_2\mathbf{x}_2^{(1)} \end{bmatrix}$ ，

$F_4\mathbf{x}_4^{(1)}=\begin{bmatrix} F_2&D_2F_2\\ F_2&-D_2F_2 \end{bmatrix}P_4\mathbf{x}_4^{(1)}=\begin{bmatrix} F_2&D_2F_2\\ F_2&-D_2F_2 \end{bmatrix}\begin{bmatrix} \mathbf{x}_2^{(2)}\\ \mathbf{x}_2^{(3)} \end{bmatrix}=\begin{bmatrix} F_2\mathbf{x}_2^{(2)}+D_2F_2\mathbf{x}_2^{(3)}\\ F_2\mathbf{x}_2^{(2)}-D_2F_2\mathbf{x}_2^{(3)} \end{bmatrix}$ ，

其中 $D_2=\begin{bmatrix} 1&0\\ 0&w^2 \end{bmatrix}$ 。因為 $F_2=\begin{bmatrix} 1&1\\ 1&w^4 \end{bmatrix}=\begin{bmatrix} 1&1\\ 1&-1 \end{bmatrix}$ ，很容易算得 $F_2\mathbf{x}_2^{(j)}$ ， $j=0,1,2,3$ 。

綜合以上討論，快速傅立葉轉換的分治運算可表示成下列樹狀圖：

$\begin{matrix} &&& F_8\mathbf{x}_8 &&&\\ && \nearrow && \nwarrow &&\\ &F_4\mathbf{x}_4^{(0)}&&&&F_4\mathbf{x}_4^{(1)}&\\ \nearrow&&\nwarrow&&\nearrow&&\nwarrow\\ F_2\mathbf{x}_2^{(0)}&&F_2\mathbf{x}_2^{(1)}&&F_2\mathbf{x}_2^{(2)}&&F_2\mathbf{x}_2^{(3)} \end{matrix}$

在實際計算時，有兩點需要注意。第一，各層的輸入向量如下表所示，括弧內是各標記的二進位表示：

$\begin{matrix} \underline{\mathrm{~Natural~Order~}} & \underline{\mathrm{~1st~Level~}}&\underline{\mathrm{~2nd~Level~}}\\ 0(000)&0(000)&0(000)\\ 1(001)&2(010)&\underline{4(100)}\\ 2(010)&4(100)&2(010)\\ 3(011)&\underline{6(110)}&\underline{6(110)}\\ 4(100)&1(001)&1(001)\\ 5(101)&3(011)&\underline{5(101)}\\ 6(110)&5(101)&3(011)\\ \underline{7(111)}&\underline{7(111)}&\underline{7(111)}\end{matrix}$

第二層的輸入向量即為初始順序，仔細觀察後發現它們不過就是將自然排序 (natural order) 的二進位數字顛倒後的結果。第一層的輸入向量順序則是固定第二層的最低有效位 (least significant bit, LSB)，將其餘數字顛倒後的結果。換句話說，按照上述規律即可得到各層所需的輸入向量，所以無須執行排列矩陣乘法。第二， $D_{\frac{m}{2}}$ 是對角矩陣，故不必逐項計算乘法。以 $D_4$ 為例，令 $\mathbf{d}_4=(1,w,w^2,w^3)^T$ ，使用 Hadamard 乘積 (或稱分元乘積) 即可，如下：

$D_4\mathbf{z}=\begin{bmatrix} 1&&&\\ &w&&\\ &&w^2\\ &&&w^3 \end{bmatrix}\begin{bmatrix} z_0\\ z_1\\ z_2\\ z_3 \end{bmatrix}=\mathbf{d}_4\circ\mathbf{z}=\begin{bmatrix} z_0\\ wz_1\\ w^2z_2\\ w^3z_3 \end{bmatrix}$ 。

下圖顯示 $n=8$ 的 Cooley-Tukey 演算法流程。若 $n=2^r$ ，快速傅立葉轉換可分解成 $r=\log_2n$ 層，每一層使用 $n/2$ 個乘法 (圖中未標記乘數的直線表示「乘以1」，因此無須計算)，所以總乘法運算量為 $\frac{n}{2}\log_2n$ 。

FFT

參考來源：
[1] 維基百科：Cooley-Tukey FFT algorithm
[2] 維基百科：FFT

6 Responses to 快速傅立葉轉換

Watt Lin says:

05/26/2012 at 9:07 am

2005年7月，台中工業區發生大火。
行政院環保署中區環境督導大隊等單位獲報趕往，初步研判刺鼻濃煙的成分應是廠房內亞硝酸鈉導致，緊急調派「傅立葉轉換紅外線光譜儀」到下風處，檢測空氣中是否含有毒物質，再據此決定已飄向市區的濃煙是否會危及市民健康，視狀況對市民提出警告，或甚至進一步的疏散動作。
(以上網頁版新聞，摘錄自 http://www.libertytimes.com.tw/2005/new/jul/4/today-t1.htm )

記得當年，我看紙本新聞，有寫「快速傅立葉轉換光譜儀」，3秒鐘，立即分析，知道空氣中所含的化學物質是什麼，它對人體不會引起重大危害。短時間可以決策，不必緊急疏散居民，可防止過度恐慌。

一般民眾，很少注意到「快速傅立葉轉換」，也不太會去思考「快速」兩字的含義。若有人寫「科普」文章，簡介「快速傅立葉轉換」，可讓大家知道，數學的進展，在生活中，有重要應用。

- ccjou says:
  
  05/28/2012 at 12:28 pm
  
  據我所知，寫過科普文章的大學老師不多 (自然數理可能多些，工程則少之又少)。一方面因為寫科普文章不會提昇學術地位，另一方面可能是在學校久了，不習慣將社會大眾當作聽眾。
  
  若想要說服他們改變…這個…有點難。我的一位同事問我：「恐怖份子與大學老師有何不同？」一時無語。他接著又說：「恐怖份子是negotiable，大學老師則不。」聞罷，再度無語。
  
  - Watt Lin says:
    
    05/28/2012 at 3:46 pm
    
    寫給一般民眾看，真的不容易。
    如果文章對象是高中生，FFT寫淺一點，過程不必詳盡，主要講歷史背景，FT的源流，new idea如何發生？高中教材有講的式子，可多著墨，例如n次方單位根在複數平面上，專家巧妙運用於FFT，快速計算答案。
    假設我當年唸高中時，曾經看到FFT的粗略介紹，那麼，看見 $cos \theta + i sin \theta$ 的時候，會增添一些學習樂趣，而非僅為了應付考試，生硬記公式。
    
  - Watt Lin says:
    
    05/28/2012 at 4:08 pm
    
    若在電機工程的課堂上，請大學生各自回憶高中時期
    $cos \theta + i sin \theta$
    是否很無聊？不知道有什麼用處？沒興趣學習？
    
Chenlogy says:

05/27/2012 at 9:35 pm

教授 ~ 印象中 ~ 快速FT 有兩種~ 一種是分時! 一種分頻! 都能以矩陣表達~ 但是結構樹略有不同~有空我在詳細與您討論! ( 不管到哪處處都可見到數學 ^_^ )

- ccjou says:
  
  05/28/2012 at 9:40 am
  
  FFT 的種類非常多，還有不少書專門討論它。我不是這方面專家，所以僅介紹最常見的Cooley-Tukey算法。主要目的是利用線性代數解構傅立葉矩陣。

	jianglong on Strassen 演算法──分治矩陣乘法
	jianglong on Strassen 演算法──分治矩陣乘法
	xmj on 內積的定義
	Ning ChingSan on 線性代數的第一堂課──矩陣乘法的定義
	momo on 兩岸線性代數用詞參照
	訪客 on 克拉瑪公式的簡易幾何證明