離散傅立葉轉換

本文的閱讀等級:中級

考慮定義於區間 [-T/2,T/2]T-週期函數 f(t) 的指數傅立葉級數 (見“傅立葉級數 (下)”):

\displaystyle F(t)=\sum_{k=-\infty}^{\infty}c_ke^{2\pi ikt/T}

其中複傅立葉係數為

\displaystyle c_k=\frac{1}{T}\int_{-T/2}^{T/2}f(t)e^{-2\pi ikt/T}dt,~~k\in\mathbb{Z}

f(t) 滿足 Dirichlet 條件 (見“傅立葉級數 (上)”),可以證明 \Vert f-F\Vert^2=0,在此情況下,往後我們不再區分函數 f(t) 與其傅立葉級數 F(t),而一律以 f(t) 表示。本文將解除 f(t) 是週期函數的限制,以下僅假設 \int_{-\infty}^{\infty}\vert f(t)\vert dt 是有界的。當 T\to\infty,傅立葉級數可推廣為傅立葉轉換 (Fourier transform)。還有,若 f(t) 不再是連續函數而是一有限數列,傅立葉級數又可延伸為離散傅立葉轉換 (discrete Fourier transform,簡稱 DFT)。本文將介紹這兩種轉換的推導過程,並解說離散傅立葉轉換的線性代數性質。

 
傅立葉轉換是傅立葉級數於 T\to\infty 的表達形式。令 f_T(t) 為一 T-週期函數,對於 -T/2\le t\le T/2f_T(t)=f(t)。令 \nu_k=k/T,且 \Delta\nu=\nu_{k+1}-\nu_k=1/T,則 f_T(t) 的傅立葉級數為

\displaystyle f_T(t)=\sum_{k=-\infty}^{\infty}c_k(T)e^{2\pi i\nu_kt}

其中我們令

\displaystyle c_k(T)=\frac{1}{T}\int_{-T/2}^{T/2}f_T(t)e^{-2\pi i\nu_kt}dt

c_k(T) 代回傅立葉級數 f_T(t),可得

\displaystyle\begin{aligned}  f_T(t)&=\sum_{k=-\infty}^{\infty}\left[\frac{1}{T}\int_{-T/2}^{T/2}f_T(t)e^{-2\pi i\nu_kt}dt\right]e^{2\pi i\nu_kt}\\  &=\sum_{k=-\infty}^{\infty}\left[\int_{-T/2}^{T/2}f_T(t)e^{-2\pi i\nu_kt}dt\right]e^{2\pi i\nu_kt}\Delta\nu,\end{aligned}

注意,上式中括弧內的 t 是虛擬變數 (dummy variable)。當 T 逐漸增大時,發生兩件事:f(t)f_T(t) 於擴大區間 [-T/2,T/2] 內一致;\Delta\nu=1/T 變成一微小量 d\nu,使得 \nu_k 逼近連續量 \nu,總和因此趨於積分。令 T\to\infty,則 f_T(t)\to f(t),可得下列等式:

\displaystyle f(t)=\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty}f(t)e^{-2\pi i\nu t}dt\right]e^{2\pi i\nu t}d\nu

欲看清整個關係,將上式拆開為兩部分。令括弧內的函數為

\displaystyle \hat{f}(\nu)=\int_{-\infty}^{\infty}f(t)e^{-2\pi i\nu t}dt

稱作 f(t) 的傅立葉轉換,其中變數 \nu 代表頻率。因此得到

\displaystyle f(t)=\int_{-\infty}^{\infty}\hat{f}(\nu)e^{2\pi i\nu t}d\nu

稱為 \hat{f}(\nu) 的逆傅立葉轉換。傅立葉級數與傅立葉轉換之間的關係可以簡單解釋如下:將傅立葉轉換的微小面積 \hat{f}(\nu)d\nu 視同傅立葉係數 c_k,當 T\to\infty,傅立葉係數變成 f(t) 的傅立葉轉換:

\displaystyle c_k=\frac{1}{T}\int_{-T/2}^{T/2}f(t)e^{-2\pi ikt/T}dt\to\hat{f}(\nu)=\int_{-\infty}^{\infty}f(t)e^{-2\pi i\nu t}dt

傅立葉級數則變成 \hat{f}(\nu) 的逆傅立葉轉換:

\displaystyle f(t)=\sum_{k=-\infty}^{\infty}c_ke^{2\pi ikt/T}\to f(t)=\int_{-\infty}^{\infty}\hat{f}(\nu)e^{2\pi i\nu t}d\nu

 
在數值計算上,因為電腦僅能處理有限維向量,連續型態的傅立葉轉換必須換裝成離散傅立葉轉換。由指數傅立葉級數同樣可推導出離散傅立葉轉換。首先,離散傅立葉轉換的輸入不再是無限維函數 f(t),我們僅能運用一週期 [0,T] 內的 n 個等間距函數值 x_j\equiv f(t_j),其中 t_j=jT/nj=0,1,\ldots,n-1。其次,離散傅立葉轉換的輸出是複傅立葉係數構成的相同長度數列 c_kk=0,1,\ldots,n-1,而非無窮數列。因為傅立葉係數 c_k 與函數 f(t) 的關係是線性的 (積分是線性函數),我們也預期離散傅立葉轉換是數列 \{x_j\}\{c_k\} 之間的一可逆線性轉換,故可表示為 n\times n 階可逆矩陣。在連續情況下,將區間正規化,令 T=1,複傅立葉級數即為

\displaystyle c_k=\int_{0}^{1}f(t)e^{-2\pi ikt}dt

在離散情況下,以 t_j=j/n 取代 t,將連續函數的傅立葉係數積分改成計算和,可得對應的離散公式:對於 k=0,1,\ldots,n-1

\displaystyle y_k=\sum_{j=0}^{n-1}x_je^{-2\pi ikj/n}

此即離散傅立葉轉換。為了與傅立葉係數 c_k 區分,我們以 y_k 代表離散傅立葉轉換結果。使用歐拉公式 e^{i\theta}=\cos\theta+i\sin\theta,令

\displaystyle w=e^{-2\pi i/n}=\cos\left(\frac{2\pi}{n}\right)-i\sin\left(\frac{2\pi}{n}\right)

離散傅立葉轉換可表示成矩陣形式 \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 稱為傅立葉矩陣。明顯地,F=[f_{pq}] 是一複對稱矩陣,f_{pq}=w^{pq},我們標記各元列行指標 p,q\in\{0,1,\ldots,n-1\}。另外,F 也是一種特殊的 Vandermonde 矩陣 (見“特殊矩陣 (8):Vandermonde 矩陣”)。離散傅立葉轉換要能成為真正有用的訊號處理工具,傅立葉矩陣 F 必須具備兩個性質:第一,不論 F 或其逆矩陣 F^{-1} 都有簡單形式。稍後我們會說明兩者形式極其類似,並可使用相同演算法運算。第二,FF^{-1} 的矩陣乘法要能夠快速運算。這個算法稱為快速傅立葉轉換 (Fast Fourier transform,簡稱 FFT),他日將另文詳細介紹。

 
傅立葉矩陣 F 是一個型態特殊的複矩陣,因此有必要進一步了解其中性質。有些作者定義 w=e^{2\pi i/n},選擇共軛複數純屬個人偏好,在此我們採用負號的原因是為了與傅立葉轉換 \hat{f}(\nu) 一致。複數集 \{1,w,w^2,\ldots,w^{n-1}\} 為方程式 z^n=1 的根,稱為 n 次方單位根,這些根在複數平面上單位圓按順時針旋轉並以相同間距排列。若 k\ge n,則 w^k=w^{k\!\mod n},其中 k\!\mod n 代表 k 除以 n 的餘數。譬如,當 n=6,就有 w^6=1, w^7=w, w^8=w^2,\ldots。對於任一整數 k,由歐拉公式可確認

\displaystyle  \overline{w^k}=\overline{e^{-2\pi ik/n}}=\cos\left(\frac{2\pi k}{n}\right)+i\sin\left(\frac{2\pi k}{n}\right)=e^{2\pi ik/n}=w^{-k}

考慮下式

w^k(1+w^k+w^{2k}+\cdots+w^{(n-2)k}+w^{(n-1)k})=w^k+w^{2k}+\cdots+w^{(n-1)k}+1

即知 (w^k-1)(1+w^k+w^{2k}+\cdots+w^{(n-2)k}+w^{(n-1)k})=0。若 k\!\mod n\neq 0,則 w^k\neq 1,可推得

\displaystyle\sum_{j=0}^{n-1}w^{jk}=1+w^k+w^{2k}+\cdots+w^{(n-2)k}+w^{(n-1)k}=0

下面我們運用此等式計算逆傅立葉矩陣 F^{-1}

 
\mathbf{f}_p\mathbf{f}_q 代表 F 的第 p 行和第 q 行,若 p\neq q,則 (q-p)\!\mod n\neq 0,兩向量的內積為

\displaystyle\left\langle\mathbf{f}_p,\mathbf{f}_q\right\rangle\overset{\underset{\mathrm{def}}{}}{=}\mathbf{f}_p^{\ast}\mathbf{f}_q=\sum_{j=0}^{n-1}\overline{w^{pj}}w^{jq}=\sum_{j=0}^{n-1}w^{j(q-p)}=0

而且

\displaystyle\Vert\mathbf{f}_p\Vert^2=\mathbf{f}_p^{\ast}\mathbf{f}_p=\sum_{j=0}^{n-1}\overline{w^{pj}}w^{jp}=\sum_{j=0}^{n-1}1=n

得知 \Vert\mathbf{f}_p\Vert=\sqrt{n}。若將 F 的各行向量予以正規化,\frac{1}{\sqrt{n}}F 為一么正 (unitary) 矩陣 (見“特殊矩陣 (3):么正矩陣(酉矩陣)”)。因為 F^T=F,就有

\displaystyle\left(\frac{1}{\sqrt{n}}F\right)^{-1}=\left(\frac{1}{\sqrt{n}}F\right)^{\ast}=\frac{1}{\sqrt{n}}\overline{F}

F^{-1}=\overline{F}/n,即 \left(F^{-1}\right)_{pq}=w^{-pq}/n,或明確寫出

\displaystyle F^{-1}=\frac{1}{n}\begin{bmatrix}  1&1&1&\cdots&1\\  1&w^{-1}&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}

所以,逆離散傅立葉轉換 \mathbf{x}=F^{-1}\mathbf{y} 可表示如下:對於 j=0,1,\ldots,n-1

\displaystyle x_j=\frac{1}{n}\sum_{k=0}^{n-1}y_ke^{2\pi ikj/n}

 
離散傅立葉轉換 F\mathbf{x} 和逆轉換 F^{-1}\mathbf{y} 有相同的運算方式,因為

\displaystyle   F^{-1}\mathbf{y}=\frac{\overline{F}\mathbf{y}}{n}=\frac{\overline{F\overline{\mathbf{y}}}}{n}

假設吾人設計出一離散傅立葉轉換程序 \mathrm{DFT}(\cdot),則逆轉換 \mathbf{x}=F^{-1}\mathbf{y} 的算法如下:

  1. 計算共軛向量 \overline{\mathbf{y}}
  2. 計算離散傅立葉轉換 \mathbf{z}=\mathrm{DFT}(\overline{\mathbf{y}})
  3. 計算共軛向量的純量乘法 \mathbf{x}=(1/n)\overline{\mathbf{z}}

 
離散傅立葉轉換的計算以矩陣運算實現,其理論則根基於內積空間。見下例,當 n=4w=e^{-2\pi i/4}=\cos(\pi/2)-i\sin(\pi/2)=-i,傅立葉矩陣為

\displaystyle F=\begin{bmatrix}  1&1&1&1\\  1&(-i)&(-i)^2&(-i)^3\\  1&(-i)^2&(-i)^4&(-i)^6\\  1&(-i)^3&(-i)^6&(-i)^9  \end{bmatrix}=\left[\!\!\begin{array}{crrr}  1&1&1&1\\  1&-i&-1&i\\  1&-1&1&-1\\  1&i&-1&-i  \end{array}\!\!\right]

逆矩陣則為

\displaystyle  F^{-1}=\frac{1}{4}\begin{bmatrix}  1&1&1&1\\  1&(i)&(i)^2&(i)^3\\  1&(i)^2&(i)^4&(i)^6\\  1&(i)^3&(i)^6&(i)^9  \end{bmatrix}=\frac{1}{4}\left[\!\!\begin{array}{crrr}  1&1&1&1\\  1&i&-1&-i\\  1&-1&1&-1\\  1&-i&-1&i  \end{array}\!\!\right]

傅立葉矩陣 F 的共軛行向量(或列向量,因為 F^T=F

\overline{\mathbf{f}_0}=\begin{bmatrix}  1\\  1\\  1\\  1  \end{bmatrix},~\overline{\mathbf{f}_1}=\left[\!\!\begin{array}{r}  1\\  i\\  -1\\  -i  \end{array}\!\!\right],~\overline{\mathbf{f}_2}=\left[\!\!\begin{array}{r}  1\\  -1\\  1\\  -1  \end{array}\!\!\right],~\overline{\mathbf{f}_3}=\left[\!\!\begin{array}{r}  1\\  -i\\  -1\\  i  \end{array}\!\!\right]

構成一組 C^4 的正交基底 \mathfrak{B},輸入數列 \mathbf{x} 可分解為這些基底向量的線性組合:

\displaystyle\begin{aligned}  \mathbf{x}&=\frac{\left\langle\overline{\mathbf{f}_0},\mathbf{x}\right\rangle}{\Vert\overline{\mathbf{f}_0}\Vert^2}\overline{\mathbf{f}_0}+\frac{\left\langle\overline{\mathbf{f}_1},\mathbf{x}\right\rangle}{\Vert\overline{\mathbf{f}_1}\Vert^2}\overline{\mathbf{f}_1}+\frac{\left\langle\overline{\mathbf{f}_2},\mathbf{x}\right\rangle}{\Vert\overline{\mathbf{f}_2}\Vert^2}\overline{\mathbf{f}_2}+\frac{\left\langle\overline{\mathbf{f}_3},\mathbf{x}\right\rangle}{\Vert\overline{\mathbf{f}_3}\Vert^2}\overline{\mathbf{f}_3}\\  &=\frac{1}{4}\left((\mathbf{f}_0^T\mathbf{x})\overline{\mathbf{f}_0}+(\mathbf{f}_1^T\mathbf{x})\overline{\mathbf{f}_1}+(\mathbf{f}_2^T\mathbf{x})\overline{\mathbf{f}_2}+(\mathbf{f}_3^T\mathbf{x})\overline{\mathbf{f}_3}\right).\end{aligned}

因為

\mathbf{y}=\begin{bmatrix}  y_0\\  y_1\\  y_2\\  y_3  \end{bmatrix}=F\mathbf{x}=\begin{bmatrix}  \mathbf{f}^T_0\\  \mathbf{f}^T_1\\  \mathbf{f}^T_2\\  \mathbf{f}^T_3  \end{bmatrix}\mathbf{x}=\begin{bmatrix}  \mathbf{f}^T_0\mathbf{x}\\  \mathbf{f}^T_1\mathbf{x}\\  \mathbf{f}^T_2\mathbf{x}\\  \mathbf{f}^T_3\mathbf{x}  \end{bmatrix}

可得

\displaystyle\begin{bmatrix}  x_0\\  x_1\\  x_2\\  x_3  \end{bmatrix}=\frac{1}{4}\left(y_0\begin{bmatrix}  1\\  1\\  1\\  1  \end{bmatrix}+y_1\left[\!\!\begin{array}{r}  1\\  i\\  -1\\  -i  \end{array}\!\!\right]+y_2\left[\!\!\begin{array}{r}  1\\  -1\\  1\\  -1  \end{array}\!\!\right]+y_3\left[\!\!\begin{array}{r}  1\\  -i\\  -1\\  i  \end{array}\!\!\right]\right)

從線性代數觀點,離散傅立葉轉換的輸出數列 y_0,y_1,y_2,y_3 除以一常數 n=4,就是輸入數列 x_0,x_1,x_2,x_3 參考基底 \mathfrak{B} 的座標,或者說,離散傅立葉轉換得到的數列向量 \mathbf{y}/4 即為輸入 \mathbf{x} 參考基底 \mathfrak{B} 的座標向量 [\mathbf{x}]_{\mathfrak{B}}。考慮下列轉換對:

\mathbf{x}=\begin{bmatrix}  1\\  1\\  0\\  0  \end{bmatrix}~~\rightleftharpoons~~\mathbf{y}=\begin{bmatrix}  2\\  1-i\\  0\\  1+i  \end{bmatrix}

離散傅立葉轉換將實數列 \{x_j\} 轉換成複數列 \{y_k\} 似乎是自找麻煩。付出這個代價,能換來什麼好處?下一回我們將探討背後隱藏的細微道理。

繼續閱讀:
Advertisements
本篇發表於 線性代數專欄, 應用之道 並標籤為 , , , , , 。將永久鏈結加入書籤。

發表迴響

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

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s