快速傅立葉轉換

本文的閱讀等級:中級

給定一序列 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 nF 稱為傅立葉矩陣。若採用一般矩陣乘法運算,離散傅立葉轉換的計算複雜度為 \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}=1n/2 次方單位根。同樣道理,對於 m=2^kk 是正整數,\{1,w^m,w^{2m},\ldots,w^{n-m}\}z^{n/m}=1n/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 nk 是正整數,因為 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

About these ads
This entry was posted in 應用之道, 主題專欄 and tagged , , , , . Bookmark the permalink.

6 則回應給 快速傅立葉轉換

  1. Watt Lin 說:

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

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

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

    • ccjou 說:

      據我所知,寫過科普文章的大學老師不多 (自然數理可能多些,工程則少之又少)。一方面因為寫科普文章不會提昇學術地位,另一方面可能是在學校久了,不習慣將社會大眾當作聽眾。

      若想要說服他們改變…這個…有點難。我的一位同事問我:「恐怖份子與大學老師有何不同?」一時無語。他接著又說:「恐怖份子是negotiable,大學老師則不。」聞罷,再度無語。

      • Watt Lin 說:

        寫給一般民眾看,真的不容易。
        如果文章對象是高中生,FFT寫淺一點,過程不必詳盡,主要講歷史背景,FT的源流,new idea如何發生?高中教材有講的式子,可多著墨,例如n次方單位根在複數平面上,專家巧妙運用於FFT,快速計算答案。
        假設我當年唸高中時,曾經看到FFT的粗略介紹,那麼,看見 cos \theta + i sin \theta  的時候,會增添一些學習樂趣,而非僅為了應付考試,生硬記公式。

      • Watt Lin 說:

        若在電機工程的課堂上,請大學生各自回憶高中時期
        cos \theta + i sin \theta
        是否很無聊?不知道有什麼用處?沒興趣學習?

  2. Chenlogy 說:

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

    • ccjou 說:

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

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 變更 )

Twitter picture

You are commenting using your Twitter account. Log Out / 變更 )

Facebook照片

You are commenting using your Facebook account. Log Out / 變更 )

Google+ photo

You are commenting using your Google+ account. Log Out / 變更 )

連結到 %s