離散傅立葉轉換

Posted on 04/19/2012 by ccjou

本文的閱讀等級：中級

考慮定義於區間 [-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/2，f_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/n，j=0,1,\ldots,n-1。其次，離散傅立葉轉換的輸出是複傅立葉係數構成的相同長度數列 c_k，k=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 n 階 F 稱為傅立葉矩陣。明顯地，F=[f_{pq}] 是一複對稱矩陣，f_{pq}=w^{pq}，我們標記各元列行指標 p,q\in\{0,1,\ldots,n-1\}。另外，F 也是一種特殊的 Vandermonde 矩陣 (見“特殊矩陣 (8)：Vandermonde 矩陣”)。離散傅立葉轉換要能成為真正有用的訊號處理工具，傅立葉矩陣 F 必須具備兩個性質：第一，不論 F 或其逆矩陣 F^{-1} 都有簡單形式。稍後我們會說明兩者形式極其類似，並可使用相同演算法運算。第二，F 和 F^{-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=4，w=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\} 似乎是自找麻煩。付出這個代價，能換來什麼好處？下一回我們將探討背後隱藏的細微道理。

繼續閱讀：
This entry was posted in 線性代數專欄, 應用之道 and tagged , , , , , . Bookmark the permalink.

1 Response to 離散傅立葉轉換

  1. 黃靖凱 says:
    11/12/2017 at 8:58 pm

    老師
    想問一個問題，當初在推導傅立葉級數時是以sin與cos作為基底的向量空間，
    但當T趨近於無窮大時所出現的傅立葉轉換，這樣的基底還存在嗎?
    即向量空間的基底可以是一個不可數的集合 B={exp{2pi i v t}, v為任意時數}
    另外
    可以把拉普拉斯轉換當成是週期從0到無窮大的週期函數的傅立葉轉換嗎?

    Reply

Leave a comment