摺積與離散傅立葉轉換

本文的閱讀等級：中級

在電機工程領域中，特別是訊號處理、電路和控制系統，線性非時變系統 (linear time-invariant system，簡稱 LTI 系統) 是一種廣泛使用的系統模型。LTI 系統按照所處理的訊號分為兩大類：連續時間系統和離散時間系統。本文從線性代數觀點介紹離散時間 LTI 系統的基本運算：摺積 (卷積，convolution)，並說明如何利用離散傅立葉轉換計算摺積。

考慮兩個 $(n-1)$ 次實多項式

$\begin{aligned} p(t)&=a_0+a_1t+\cdots+a_{n-1}t^{n-1},\\ q(t)&=b_0+b_1t+\cdots+b_{n-1}t^{n-1},\end{aligned}$

乘積 $p(t)q(t)$ 為一 $(2n-2)$ 次多項式，稱為柯西乘積 (Cauchy product)：

$p(t)q(t)=c_0+c_1t+\cdots+c_{2n-2}t^{2n-2}$ ，

第 $k$ 次項由 $(a_jt^j)(b_{k-j}t^{k-j})$ 構造而成，其係數為

$\displaystyle c_k=\sum_{j=0}^ka_{j}b_{k-j},~~k=0,1,\ldots,2n-2$ 。

上式中， $a_j=b_j=0$ ，若 $j\ge n$ 。將多項式 $p(t)$ 和 $q(t)$ 的係數合併成 $n$ 維係數向量

$\mathbf{a}=\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix},~~\mathbf{b}=\begin{bmatrix} b_0\\ b_1\\ \vdots\\ b_{n-1} \end{bmatrix},$

這裡 $a_j$ 和 $b_j$ 的元指標從0開始算起，以下所有向量遵守同樣規則。同樣地，乘積 $p(t)q(t)$ 的係數 $c_k$ 亦可合併為係數向量，如下：

$\mathbf{a}\circledast\mathbf{b}\overset{\underset{\mathrm{def}}{}}{=}\begin{bmatrix} a_0b_0\\ a_0b_1+a_1b_0\\ a_0b_2+a_1b_1+a_2b_0\\ \vdots\\ a_{n-2}b_{n-1}+a_{n-1}b_{n-2}\\ a_{n-1}b_{n-1}\\ 0 \end{bmatrix}_{2n\times 1}$ ，

稱為 $\mathbf{a}$ 和 $\mathbf{b}$ 的摺積。我們加入最末0元的用意在使 $\mathbf{a}\circledast\mathbf{b}$ 成為 $2n$ 維向量，便於之後平衡公式。此外，設想 $\mathbf{a}$ 後面隱含 $n$ 個零元， $a_n=\cdots=a_{2n-1}=0$ ；同樣地， $\mathbf{b}$ 後面也隱含 $n$ 個零元， $b_n=\cdots=b_{2n-1}=0$ 。如此一來，摺積的每一元都有相同的表達式：

$\displaystyle (\mathbf{a}\circledast\mathbf{b})_k=\sum_{j=0}^ka_jb_{k-j},~~k=0,1,\ldots,2n-1$ 。

摺積各元的運算公式和一般實向量內積有兩個不同的地方：第一， $\mathbf{a}$ 和 $\mathbf{b}$ 有相反的指標序；第二，從第0元開始算起，摺積的第 $k$ 元加總 $k+1$ 個乘積 (其中涵蓋乘積0)。根據這兩項差異，我們可以設計一個視覺化摺積計算法：將向量 $\mathbf{a}$ 和 $\mathbf{b}$ 想像成兩列迎向對駛的火車，各自有 $n$ 節車廂，從車頭算起，編號為 $0,1,\ldots,n-1$ 。當兩列火車會車時， $a_0$ 對齊 $b_k$ 且 $b_0$ 對齊 $a_k$ ，加總位置對齊的 $k+1$ 節車廂乘積即得 $(\mathbf{a}\circledast\mathbf{b})_{k}$ ，如下所示。

$k=0:$

$\begin{matrix} ~&~&~&~&a_0&a_1&a_2&\cdots&a_{n-1}\\ ~&~&~&~&\times&~&~&~&~\\ b_{n-1}&\cdots&b_2&b_1&b_0&~&~&~&~ \end{matrix}$

$k=1:$

$\begin{matrix} ~&~&~&~&a_0&a_1&a_2&\cdots&a_{n-1}\\ ~&~&~&~&\times&\times&~&~&~\\ ~&b_{n-1}&\cdots&b_2&b_1&b_0&~&~&~ \end{matrix}$

$k=2:$

$\begin{matrix} ~&~&~&~&a_0&a_1&a_2&\cdots&a_{n-1}\\ ~&~&~&~&\times&\times&\times&~&~\\ ~&~&b_{n-1}&\cdots&b_2&b_1&b_0&~&~ \end{matrix}$

$\vdots$
$k=2n-2:$

$\begin{matrix} a_0&a_1&a_2&\cdots&a_{n-1}&~&~&~&~\\ ~&~&~&~&\times&~&~&~&~\\ ~&~&~&~&b_{n-1}&\cdots&b_2&b_1&b_0 \end{matrix}$

對於任意 $n$ 維向量 $\mathbf{a},\mathbf{b},\mathbf{c}$ ，純量 $\alpha$ ，運用視覺化算法不難推知摺積具有下列性質：

(1) 交換律： $\mathbf{a}\circledast\mathbf{b}=\mathbf{b}\circledast\mathbf{a}$

(2) 分配律： $(\mathbf{a}+\mathbf{b})\circledast\mathbf{c}=\mathbf{a}\circledast\mathbf{c}+\mathbf{b}\circledast\mathbf{c}$

(3) 純量結合律： $\alpha(\mathbf{a}\circledast\mathbf{b})=(\alpha\mathbf{a})\circledast\mathbf{b}=\mathbf{a}\circledast(\alpha\mathbf{b})$

備妥摺積的基本性質，現在是線性代數上場的時候。首先，我們將摺積視為一變換，令 $T(\mathbf{x})=\mathbf{x}\circledast\mathbf{h}$ ，其中 $\mathbf{x}=[x_i]$ 和 $\mathbf{h}=[h_i]$ 同是 $n$ 維實向量。利用分配律可得

$T(\mathbf{x}+\mathbf{y})=(\mathbf{x}+\mathbf{y})\circledast\mathbf{h}=(\mathbf{x}\circledast\mathbf{h})+(\mathbf{y}\circledast\mathbf{h})=T(\mathbf{x})+T(\mathbf{y})$ ，

利用純量結合律可得

$T(\alpha\mathbf{x})=(\alpha\mathbf{x})\circledast\mathbf{h}=\alpha(\mathbf{x}\circledast\mathbf{h})=\alpha T(\mathbf{x})$ ，

因此證明 $T(\mathbf{x})=\mathbf{x}\circledast\mathbf{h}$ 為一從 $\mathbb{R}^n$ 映至 $\mathbb{R}^{2n}$ 的線性變換。接著我們寫出摺積 $T$ 的矩陣表達式。考慮 $\mathbb{R}^n$ 的標準基底 $\{\mathbf{e}_0,\mathbf{e}_1,\ldots,\mathbf{e}_{n-1}\}$ ，其中 $\mathbf{e}_i$ 的第 $i$ 元等於1，其餘各元為零。對於 $i=0,1,\ldots,n-1$ ，運用視覺化算法可知標準單位向量 $\mathbf{e}_i$ 的像為

$T(\mathbf{e}_i)=\mathbf{e}_i\circledast\mathbf{h}=\begin{bmatrix} 0\\ \vdots\\ h_0\\ h_1\\ \vdots\\ h_{n-1}\\ \vdots\\ 0 \end{bmatrix}_{2n\times 1},$

其中 $h_0$ 出現在 $T(\mathbf{e}_i)$ 的第 $i$ 元位置。將向量 $\mathbf{x}$ 寫成標準基底向量 $\{\mathbf{e}_i\}$ 的線性組合： $\mathbf{x}=x_0\mathbf{e}_0+x_1\mathbf{e}_1+\cdots+x_{n-1}\mathbf{e}_{n-1}$ ，代入 $T(\cdot)$ ，利用線性變換基本性質，可得 $T(\mathbf{x})=\mathbf{x}\circledast\mathbf{h}$ 的矩陣乘法運算表達式，如下：

$\begin{aligned} T(\mathbf{x})&=T(x_0\mathbf{e}_0+x_1\mathbf{e}_1+\cdots+x_{n-1}\mathbf{e}_{n-1})\\ &=x_0T(\mathbf{e}_0)+x_1T(\mathbf{e}_1)\cdots+x_{n-1}T(\mathbf{e}_{n-1})\\ &=\begin{bmatrix} T(\mathbf{e}_0)&T(\mathbf{e}_1)&\cdots&T(\mathbf{e}_{n-1}) \end{bmatrix}\begin{bmatrix} x_0\\ x_1\\ \vdots\\ x_{n-1} \end{bmatrix}\\ &=H\mathbf{x},\end{aligned}$

其中

$H=\begin{bmatrix} h_0&0&\cdots&0\\ h_1&h_0&\cdots&0\\ h_2&h_1&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ h_{n-1}&h_{n-2}&\cdots&h_0\\ 0&h_{n-1}&\cdots&h_1\\ 0&0&\cdots&h_2\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&h_{n-1}\\ 0&0&\cdots&0 \end{bmatrix}_{2n\times n}$

即為摺積 $T$ 參考標準基底的表示矩陣 (見“線性變換表示矩陣”)。

摺積是一個重要的數學分析運算，最常見的應用在於描述離散時間 LTI 系統 $T$ 的輸入與輸出關係，圖示如下：

$\mathbf{x}\longrightarrow\boxed{~~\mathbf{h}~~}\longrightarrow\mathbf{x}\circledast\mathbf{h}$

從摺積的表示矩陣 $H$ 很容易理解離散時間 LTI 系統的特性，考慮等價圖像：

$\mathbf{x}\longrightarrow\boxed{~~H~~}\longrightarrow H\mathbf{x}$

我們將向量 $\mathbf{x}=\begin{bmatrix} x_0&x_1&\cdots&x_{n-1} \end{bmatrix}^T$ 看成是一有限序列訊號， $x_i$ 代表在時間座標 $i\ge 0$ 的輸入值。顧名思義，線性非時變系統 $T$ 有兩個特性：「線性」與「非時變」。前者說明 $T$ 是一線性變換，後者的意思是 $T$ 不隨時間改變。考慮特殊輸入訊號 $\mathbf{x}=\mathbf{e}_0=\begin{bmatrix} 1&0&0&\cdots&0 \end{bmatrix}^T$ ，稱為單位脈衝 (unit impulse)，輸出 $T(\mathbf{e}_0)$ 稱為脈衝響應 (impulse response)，即 LTI 系統 $T$ 的表示矩陣 $H$ 的第一行： $\begin{bmatrix} h_0&h_1&\cdots&h_{n-1}&0&\cdots&0 \end{bmatrix}^T$ 。脈衝響應的時間座標也從0開始，我們說它滿足因果律 (causality)，也就是說系統的輸出只與當前以及過去的輸入有關。如果延遲一單位時間後輸入單位脈衝，即 $\mathbf{x}=\mathbf{e}_1=\begin{bmatrix} 0&1&0&\cdots&0 \end{bmatrix}^T$ ，輸出 $T(\mathbf{e}_1)$ 為 $H$ 的第二行，脈衝響應內容不變，但往後推遲一單位時間。從以上討論可知線性非時變系統 $T$ 完全由脈衝響應 $\mathbf{h}$ 決定，「線性」指出 $T$ 可表示成矩陣乘法運算 $T(\mathbf{x})=H\mathbf{x}$ ，「非時變」則說明 $H$ 的每一行皆由同一脈衝響應 $\mathbf{h}$ 構成，惟各行的脈衝響應依序往下平移。

見下例，設一離散時間 LTI 系統的脈衝響應為 $\mathbf{h}=\begin{bmatrix} 2&1&1&1 \end{bmatrix}^T$ ，對於輸入訊號 $\mathbf{x}=\begin{bmatrix} 1&2&0&0 \end{bmatrix}^T$ ，其輸出訊號可由摺積定義 $\mathbf{x}\circledast\mathbf{h}$ 或矩陣乘法 $H\mathbf{x}$ 算得：

$T(\mathbf{x})=H\mathbf{x}=\begin{bmatrix} 2&0&0&0\\ 1&2&0&0\\ 1&1&2&0\\ 1&1&1&2\\ 0&1&1&1\\ 0&0&1&1\\ 0&0&0&1\\ 0&0&0&0 \end{bmatrix}\begin{bmatrix} 1\\ 2\\ 0\\ 0 \end{bmatrix}=1\begin{bmatrix} 2\\ 1\\ 1\\ 1\\ 0\\ 0\\ 0\\ 0 \end{bmatrix}+2\begin{bmatrix} 0\\ 2\\ 1\\ 1\\ 1\\ 0\\ 0\\ 0 \end{bmatrix}=\begin{bmatrix} 2\\ 5\\ 3\\ 3\\ 2\\ 0\\ 0\\ 0 \end{bmatrix}$ 。

下圖清楚解釋了摺積與 LTI 系統的運作。將輸入訊號分解成單位脈衝的線性組合，分別送入 LTI 系統，再重組 (加總) 各自的脈衝響應即得到輸出訊號。

Discrete convolution From http://cnx.org/content/m10087/2.17/dtconv4e.png

最後說明如何利用離散傅立葉轉換計算摺積。我們需要一個新的向量乘法運算，稱為 Hadamard 乘積或分元 (pairwise) 乘積，定義如下：

$\mathbf{a}\circ\mathbf{b}=\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix}\circ\begin{bmatrix} b_0\\ b_1\\ \vdots\\ b_{n-1} \end{bmatrix}\overset{\underset{\mathrm{def}}{}}{=}\begin{bmatrix} a_0b_0\\ a_1b_1\\ \vdots\\ a_{n-1}b_{n-1} \end{bmatrix}$ 。

令 $\hat{\mathbf{x}}$ 和 $\hat{\mathbf{h}}$ 表示加入 $n$ 個零元的 $2n$ 維增廣向量：

$\hat{\mathbf{x}}=\begin{bmatrix} x_0\\ \vdots\\ x_{n-1}\\ 0\\ \vdots\\ 0 \end{bmatrix}_{2n\times 1},~~\hat{\mathbf{h}}=\begin{bmatrix} h_0\\ \vdots\\ h_{n-1}\\ 0\\ \vdots\\ 0 \end{bmatrix}_{2n\times 1}$ ，

且 $F$ 為 $2n\times 2n$ 階傅立葉矩陣 (見“離散傅立葉轉換”)：

$F=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&w&w^2&\cdots&w^{2n-1}\\ 1&w^2&w^4&\cdots&w^{2(2n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&w^{2n-1}&w^{2(2n-1)}&\cdots&w^{(2n-1)^2} \end{bmatrix}$

其中 $w=e^{-2\pi i/(2n)}$ 。摺積定理指出兩向量摺積的離散傅立葉轉換即為兩向量的離散傅立葉轉換的 Hadamard 乘積：

$F(\mathbf{x}\circledast\mathbf{h})=(F\hat{\mathbf{x}})\circ(F\hat{\mathbf{h}})$ 。

證明於下。令 $\mathbf{f}_p$ 和 $\mathbf{f}_q$ 表示傅立葉矩陣 $F$ 的第 $p$ 行和第 $q$ 行，不難驗證

$\mathbf{f}_p\circ\mathbf{f}_q=\mathbf{f}_{p+q},~~p,q=0,1,\ldots,n-1$ 。

寫出 $F\hat{\mathbf{x}}$ 和 $F\hat{\mathbf{h}}$ 的線性組合表達式，利用上述性質，並使用類似多項式乘積的計算方式，可推得

$\displaystyle \begin{aligned} (F\hat{\mathbf{x}})\circ(F\hat{\mathbf{h}})&=\left(x_0\mathbf{f}_0+x_1\mathbf{f}_1+\cdots+x_{n-1}\mathbf{f}_{n-1}\right)\circ\left(h_0\mathbf{f}_0+h_1\mathbf{f}_1+\cdots+h_{n-1}\mathbf{f}_{n-1}\right)\\ &=\sum_{k=0}^{2n-2}\left(\sum_{j=0}^kx_jh_{k-j}\right)\mathbf{f}_k\\ &=\sum_{k=0}^{2n-2}(\mathbf{x}\circledast\mathbf{h})_k\mathbf{f}_k\\ &=F(\mathbf{x}\circledast\mathbf{h})\end{aligned}$

等號兩邊同乘 $F^{-1}$ ，即有

$\displaystyle\mathbf{x}\circledast\mathbf{h}=F^{-1}\left((F\hat{\mathbf{x}})\circ(F\hat{\mathbf{h}})\right)$ 。

上式指出兩向量的摺積可透過離散傅立葉轉換而得 (逆轉換 $F^{-1}$ 亦可由離散傅立葉轉換實現)，這個結果令人費解：為何將簡單的 $2n\times n$ 階矩陣乘法運算 $H\mathbf{x}$ 轉換成三次 $2n$ 階離散傅立葉轉換呢？答案於公元1965年揭曉。美國學者庫利 (James Cooley) 與圖基 (John Tukey) 提出了一個複雜度為 $\mathcal{O}(n\log_2 n)$ 的傅立葉矩陣演算法，後人稱為快速傅立葉轉換 (fast Fourier transform, 簡稱 FFT)。以一般矩陣乘法計算 $n$ 維向量摺積必須執行 $n^2$ 次乘法運算 ( $H$ 矩陣有半數為 $0$ ，故實際運算量從 $2n^2$ 縮減為 $n^2$ )，但利用快速傅立葉轉換僅需大約 $3n\log_2n$ 次乘法運算。當 $n$ 增大時，快速傅立葉轉換更顯現其優勢。究竟快速傅立葉轉換使用了甚麼秘訣，能有如此驚人的效能？這個主題將留待日後討論。

繼續閱讀：

快速傅立葉轉換

12 Responses to 摺積與離散傅立葉轉換

Watt Lin says:

06/10/2012 at 5:27 pm

請問老師：
聲音編輯軟體，淡入(Fade In)、淡出(Fade Out)效果，是採用Convolution處理嗎？
能否簡短說明，或需另外寫個長篇文章？
謝謝！

- ccjou says:
  
  06/10/2012 at 9:21 pm
  
  我以為fader只是將訊號衰減(震幅減小)，不須經過特殊處理。
  好的，我再找時間寫一篇離散傅立葉轉換與頻域分析的基本原理簡介。
  
  - Watt Lin says:
    
    06/10/2012 at 9:55 pm
    
    感謝老師的回答。
    我是外行人，好像也問得不恰當，
    發問後，看見老師答案之前，自己思考，
    淡入淡出，在時域應該可以處理，不必轉換到頻域作處理。
    我所好奇的，頻域分析，有沒有日常生活經常遇到的應用，能否作些初步簡介？
    
    - ccjou says:
      
      06/11/2012 at 8:22 am
      
      好的，但我不是訊號處理專家，還要花點時間研究一下。
      
    - 秦紀維 says:
      
      09/28/2013 at 4:17 pm
      
      醫療影像學。將影像轉換到頻域後，取高頻域再做逆轉換會凸顯細節，相對的取低頻域做逆轉換會忽略掉細節。後者用處不大，畢竟近視眼、注意力短缺、阿茲海默症就有類似功能（不正經）。
      
      - ccjou says:
        
        09/28/2013 at 8:02 pm
        
        JPEG 壓縮的量化歩驟就是保留低頻忽略高頻細節。
        
Watt Lin says:

11/04/2013 at 3:58 pm

請問老師：摺積運算，有沒有卡通動畫影片?
我在wiki看到簡單的幾張圖與動畫， http://en.wikipedia.org/wiki/Convolution
不知外國的其他網站，有沒有更仔細的動畫方式解說？

Watt Lin says:

11/04/2013 at 4:01 pm

http://en.wikipedia.org/wiki/Convolution
在wiki，摺積有幾張圖可供想像，還有個簡單動畫。
請問老師：國外其他網站，對於摺積，有沒有更仔細的卡通動畫？

- ccjou says:
  
  11/04/2013 at 6:43 pm
  
  謝謝你提供的動畫訊息。我所知道的動畫模擬也都很類似，例如，
  http://mathworld.wolfram.com/Convolution.html
  基本原理就是上文所述，兩列迎向對駛的火車內積。
  
  不知道為甚麼你的迴響被系統當作垃圾。原本我選擇的設定是「如果一則迴響中超過 3條鏈結，即判斷為垃圾，必須審核通過才可以發表。」
  
  - Watt Lin says:
    
    11/04/2013 at 7:12 pm
    
    老師提供MathWorld這個網頁連結，也很有幫助。
    我從好幾年前，偶爾看一點點有關傅立葉轉換的介紹，一直看不懂摺積。
    到了2012年，購買交通大學出版社的「訊號與系統」光碟，
    謝世福教授採取幽默的解說，我才開始懂摺積。
    他說：「吃一粒米，第一天很有力氣，第二天還剩餘一些力氣，第三天力氣微弱，並且在黑板畫上三條不同高度的線條。然後又說，第二天也吃一粒米，第二天的力氣，是第二粒米的力氣，加上第一天剩餘的力氣。第三天又吃一粒米，便有三天個別力氣的總和。」
    這種比喻法，聽起來很誇張，但我很喜歡這樣的幽默，幫助理解！
    
    - ccjou says:
      
      11/04/2013 at 8:06 pm
      
      哈哈，謝老師的每天吃米粒的譬喻比我的兩列迎向對駛的火車更傳神。
      
  - Watt Lin says:
    
    11/04/2013 at 7:17 pm
    
    如果有人講，吃一顆藥，第一天血液中藥物濃度高，第二天降至一半左右，第三天仍有一點點藥物濃度。
    第二天又吃相同的一顆藥，血中濃度，不僅受第二天藥物影響，也要加上第一天的剩餘影響。
    這樣的比喻法，沒有誇張，而且聽起來近乎學理。
    但是，藥理學的研究，好像不會用到摺積這種觀念。

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