[訊號與系統] LTI系統輸入輸出關係由 Convolution 決定 - 從離散時間觀點

以下我們討論 為何 線性非時變 (Linear Time-Invariant, LTI) 系統 輸入與輸出關係 由 所謂的 convolution 表示。為了避免過多繁雜的數學，以下僅討論離散時間的情況。首先我們需要一些定義的幫助：

========================
Definition: Unit Impulse Function in Discrete-Time (Kronecker Delta)
我們說函數 $\delta: \mathbb{N} \to \{0,1\}$ 為 unit impulse function in discrete-time time 若 $\delta$ 滿足
\[\delta \left[ n \right] = \left\{ \begin{gathered}
  1,\;\;\;\;n = 0 \hfill \\
  0,\;\;\;\; o.w. \hfill \\
\end{gathered}  \right.
\]========================

========================
Definition: Impulse Response
給定任意系統配備輸入 $x[n]$ 與輸出 $y[n]$ 關係為 $y[n] = T\{x[n]\}$其中 $T$ 視為 operator (定義在某函數空間)，若輸入為 $x[n]=\delta[n]$ 則 輸出
\[
h[n] := y[n] = T\{\delta[n]\}
\]稱為系統 $T$ 的 脈衝響應 (impulse response)
========================


========================
FACT: 任意離散訊號 $x[n]$ 可由 $\delta$ 做組合疊加，亦即
\[
x[n] = \sum_{k=-\infty}^{\infty}x[i] \delta[n-k]
\]========================
Proof: 證明顯然，在此不做贅述。$\square$



========================
Definition: Linear System
給定系統 $T$ 滿足以下輸入與輸出關係： $y_1[n]=T\{x_1[n]\}$ 且 $y_2[n] = T\{x_2[n]\}$。現在定義 $x[n] =ax_1[n] + bx_2[n]$ 其中 $a,b \in \mathbb{R}$則我們說系統 $T$ 為 linear 若下列條件成立
\[
y[n] = T\{x[n]\}
\]且 $y[n] = a y_1[n] + b y_2[n]$
========================

Remarks:
上述 系統 $T$  其實可看按作 泛函分析中的 operator。故線性系統 就是 泛函分析中的 線性算子 (linear operator)。



========================
Definition: Time-Invariant System
給定系統 $T$ 滿足以下輸入與輸出關係： $y[n]=T\{x[n]\}$。現在定義 $x[n] :=x[n-d]$  其中 $d \in \mathbb{N}$ 則我們說 $T$ 為 time-invariant 若下列條件成立
\[
y[n-d] = T\{x[n-d]\}
\]========================


有了以上定義的幫助，我們現在有辦法給出以下為何 LTI 系統的輸入輸出關係確實為 convolution。


========================
Theorem: LTI input-output Relationship is Governed by Convolution
給定 線性非時變 (LTI)系統 $L$ 滿足
\[
y[n] = L\{x[n]\}
\]其中 $L$ 為 linear operator，則
\[
y[n] = \sum_k h[k] x[n-k]
\]========================
Proof:
由於 任意 輸入 $x[n]$ 皆可被表為 脈衝函數的 組合：
\[
x[n] = \sum_{k = -\infty}^\infty x[k] \delta[n-k]
\]現在將此訊號作為LTI系統 $L$ 的輸入，則由於 $L$ 為 linear operator 我們有
\begin{align*}
  y[n] &= L\left\{ {x\left[ n \right]} \right\} \\
&= L\left\{ {\sum\limits_{k =  - \infty }^\infty  x [k]\delta [n - k]} \right\} \hfill \\
   &= \sum\limits_{k =  - \infty }^\infty  x[k] L\left\{ \delta [n - k] \right\}  \;\;\;\; (*)
\end{align*} 接著因為 $L$ 為 time-invariant 我們可進一步改寫
\[L\left\{ {\delta [n - k]} \right\} = h\left[ {n - k} \right] \;\;\;\; (**)
\]故由 $(*)$ 與 $(**)$ 可得
\[y[n] = \sum\limits_{k =  - \infty }^\infty  {x[k]h\left[ {n - k} \right]} \]即為所求。$\square$

Comments:
上述 \[
y[n] = \sum\limits_{k =  - \infty }^\infty  {x[k]h\left[ {n - k} \right]}
\]一般又稱作 $x[n]$ 與 $h[n]$ convolution sum。記作
\[
y[n] = x[n] * h[n]
\]

[訊號處理] 2d Convolution & 簡單的影像處理 (利用 MATLAB )

一般在處理影像的時候我們會把 影像(image) 視為 二維離散訊號， 更一般的說法是將其視為矩陣，亦即給定任一張 (灰階) 影像 我們可以將其分割成 $M \times N$ 矩陣 (像素 pixel)，並且將其記作 $x[m,n]$，其中 $m=0,...,M-1$ 與 $n = 0,...,N-1$ 。

Comments:
如果要處理的影像是彩色的，一個常用的做法是把該影像以 $R,G,B$ 三色分別存成 三個矩陣，最後在做疊合。在此不做贅述。

以下我們利用 MATLAB 來執行 (灰階)影像處理，我們使用的圖檔是 MATLAB內建的圖檔 cameraman.tif ，當然讀者可自行讀入任何自己想要的圖檔。在MATLAB 輸入

MATLAB Code For Loading the Image
x = imread('cameraman.tif' )
imagesc(x)

則會顯示

圖1a: cameraman.tif 原圖 ($256 \times 256$)

Comments:
1. 上述影像訊號 $x[m,n]$ 為 $M \times N = 256 \times 256 $ 矩陣。
2. 一般而言在影像處理領域，我們將圖片最左上角點視為座標原點，對應的影像訊號為 $x[0,0]$ 。
3. 上述影像透過 MATLAB 2017a 讀圖會顯示有一些藍綠等顏色，若要使此圖顯示為灰階(gray scale) 讀者可鍵入 colormap(gray) 來使其變成灰階影像，亦即

MATLAB Code For Gray Scale Colormap
x = imread('cameraman.tif' )
imagesc(x)
colormap(gray)

則我們會得到如下圖

圖1b: 灰階影像

現在我們可對上述影像進行一些常見的基本處理。令 $x[m,n]$ 為輸入影像訊號，且假設 影像 濾波器 為 線性非時變 (Linear Time Invariant, LTI) 故其對應的 impulse response $h[m,n]$ 可以用來完全描述我們的影像濾波器，最後 我們令 影像處理過後的輸出訊號 $y[m,n]$ 可表為 $x[m,n]$ 與 $h[m,n]$ 的 convolution ，差別僅在此時我們的 convolution運算為二維運算，故我們寫成
\[
y[m,n] = h[m,n]**x[m,n]: = \sum\limits_{k = 0}^{M - 1} {\sum\limits_{l = 0}^{N - 1} {h\left[ {k,l} \right]x\left[ {m - k,n - l} \right]} }
\] 其中 $**$ 表示二維 convolution，在 MATLAB 中我們可以使用 conv2 來執行此運算。下圖顯示了上述討論的觀點：

Comments:
1. 上述討論中提及的 脈衝響應 $h[m,n]$ 在影像處理中又被稱為 點擴散函數 (point spread function, PSF)
2. 給定任意 影像 $x[m,n]$ 其影像尺寸為 $M_1 \times N_1$ 且 影像濾波器  $h[m,n]$ 具有 影像尺寸為$M_2 \times N_2$ 則其經過 2d convolution之後的 輸出 $y[m,n]$ 之影像尺寸會變成
$$
(M_1 + M_2 -1) \times (N_1 + N_2 -1)
$$ 故上述 2d convolution 執行之後原圖尺寸會變成
$$
(256+2-1) \times (256 + 2 -1) = 257 \times 257
$$也就是說透過2d convolution之後 影像邊緣 會多跑出一些額外的像素。讀者可參考下方後續的討論。
3. 上述討論中我們假設 影像濾波器為 LTI ，故輸入與輸出關係為  (time-domain) convolution，那麼熟悉訊號與系統的讀者不難做出以下猜想，對輸入 $x$ 與 輸出 $y$ 與 影像濾波器 $h$ 取 Discrete Fourier Transform (DFT) 可得到 對應的 DFT係數 如 \[\begin{gathered}
  X[k,l] = DFT(x[m,n]); \hfill \\
  H[k,l] = DFT(h[m,n]); \hfill \\
  Y[k,l] = DFT(y[m,n]). \hfill \\
\end{gathered} \] 假設 $N,M$ 夠大，則其在 輸入輸出在頻域上對應的關係為 frequency domain  為相乘
$$
Y[k,l] = H[k,l] X[k,l]
$$
一但計算出 $Y[k,l]$ 在對其取 Inverse Digital Fourier Transform (IDFT) 即可得到 $y[m,n]$。在 MATLAB 中，上述討論可以使用 fft2 與 ifft2 來實現，在此不做贅述。


以下我們探討幾種常見的影像濾波器，亦即幾種常見的 $h[.]$：

Image Filtering:
以下為幾種常見的影像濾波器的例子：
1. 影像模糊 (Image blurring)
\[{h_b}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{1/10} \\
  {1/10}&{1/10}
\end{array}} \right];\]

MATLAB code for Image Blurring
h_b = 1/10 * [1 1; 1 1];
y_b = conv2(x, h_b, 'same')
imagesc( abs(y_b) )

圖2: 模糊效果

讀者可以比較此圖與之前的圖不難發現此圖較原圖為模糊。另外可注意到模糊之後的影像邊緣部分出現額外不屬於原圖的像素。


2. 邊緣偵測：圖像的邊緣可以透過特定的濾波器設計來偵測，比如說我們可以考慮
\[\begin{gathered}
  {h_h}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{1/10} \\
  { - 1/10}&{ - 1/10}
\end{array}} \right]; \hfill \\
  {h_v}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{ - 1/10} \\
  {1/10}&{ - 1/10}
\end{array}} \right]. \hfill \\
\end{gathered} \]其中 $h_h$ 用以偵測 影像的水平邊緣，$h_v$ 用以偵測影像的垂直邊緣。

MATLAB Code for Horizontal Edge Detection
h_h = 1/10 * [1 1; -1 -1];
y_h = conv2(x, h_h, 'same')
imagesc( abs(y_h) )

我們得到對於影像水平邊緣的偵測如下圖所示

圖3: 水平邊緣偵測

MATLAB Code for Vertical Edge Detection
h_v = 1/10 * [1 -1; 1 -1];
y_v = conv2(x, h_v, 'same')
imagesc( abs(y_v) )

上述code得到對於垂直邊緣的偵測如下圖

圖4: 垂直邊緣偵測

Comments:
1. 讀者可自行改變 $h_v, h_h, h_b$  的矩陣的係數來看看得到的圖有什麼不同。
2. 影像處理有非常多的有趣的問題可以進一步討論，比如說被 模糊後之後的影像是否可以把它還原？如果可以該怎麼做？ 直接用 conv2 或者用 fft2 何者運算較快？影像太大該怎麼壓縮？等等。

[訊號與系統] FIR 與 Finite Convolution

在 系統理論 或者 訊號處理 的領域中，我們一般將 濾波器 (filter) 視作可用以移除某特定頻段的訊號 並且僅讓 部分指定的頻段訊號 可以通過 的 系統(system)。

現在 給定 一組 標準 有限脈衝響應濾波器(Finite Impulse Response Filter,  FIR filter) ，其 輸入/輸出之關係可用下列表述
\[
y[n] = \sum_{k = 0}^M b_k x[n -k]
\] 其中 $b_k$ 稱為 FIR 濾波器的係數，$x[\cdot]$ 為輸入訊號，$y[\cdot]$ 為輸出訊號。

Comments:
1. 上述 FIR 的定義中並不要求 未來的輸入訊號，亦即我們僅需要現在與過去的輸入 $x[n], x[n-1],...,x[n-M]$，一般稱此 FIR 為 因果系統 (casual system)
2. 上述 FIR filter 的 輸入/輸出關係 可被視為 有限摺積(finite convolution) 運算
3. 令輸入$x[n]$為具有長度 $L_x$ 的 sequence 且 $h[n]$ 為具有 長度 $L_h$ 的sequence 則其輸出 $y[n]$ 在經過 convolution 運算之後會具有長度
$$
L_y =L_x + L_h -1 \text {( why ? ) }
$$ 4. 在 MATLAB 中 內建函數 filter() 可以用來建構 上述 FIR filter，舉例而言，考慮一組 FIR
\[y[n] = \sum\limits_{k = 0}^3 {\frac{1}{4}} x[n - k] = \frac{1}{4}\left( {x[n] + x[n - 1] + x[n - 2] + x[n-3]} \right)\] 且輸入為 $sin(0.1 \pi n), \;\;\; n=0,1,...,99$ 則 我們可用 MATLAB code 來計算對應的輸出 $y$ 如下：

n = 0:99
x = sin( 0.1 * pi * n);
b = [1/4 1/4 1/4  1/4]
y = filter(b, 1, x)


另外注意到前述我們提及 輸入訊號與其 脈衝響應的和：
\[
y[n] = \sum_{k = 0}^M h[k] x[n -k]\;\;\;\;\;\;\; (*)
\] 上述結果稱作 finite convolution sum 且我們稱輸出 透過 $x[n]$ 與 $h[n]$ 做 convolution 而得。一個 更一般的 輸入/輸出 關係可表為
\[
y[n] = \sum_{k=-\infty}^\infty h[k] x[n-k] \;\;\;\;\; (**)
\]讀者不難發現若我們取 $h[n]=0$ 當 $n<0$ 以及 $n > M$ 時， $(**)$ 退化回 finite convolution 形式 $(*)$。



Example: Impulse Response For 3-Points Averaging System
令 $x[n] := \delta[n]$ 則 脈衝響應
\[
y[n] = \sum_{k = 0}^2 \frac{1}{3} x[n -k]
\]

在 MATLAB 中，convolution 運算可以透過 conv() 來實現：比如說
xx = sin( 0.1*pi*(0:50) )
hh = ones(11,1)/11
yy = conv(hh, xx);


線性非時變 FIR 濾波器

==================
FACT: 考慮 FIR 濾波器
\[
y[n] = \sum_{k = 0}^M b_k x[n -k]
\]則此 FIR 為 線性非時變( linear and time-invariant)
==================

Proof: 首先證明線性。令 $x_1[n]$ 與 $x_2[n]$ 分別為兩輸入，且其對應的輸出分別為
\[\begin{array}{l}
{y_1}[n] = \sum\limits_{k = 0}^M {{b_k}} {x_1}[n - k];\\
{y_2}[n] = \sum\limits_{k = 0}^M {{b_k}} {x_2}[n - k];
\end{array}\]現在考慮前述兩輸入的線性組合，亦即對任意 $\alpha, \beta \in \mathbb{R}$ 定義
\[
x[n] := \alpha x_1[n] + \beta x_2[n]
\]我們觀察
\begin{align*}
y[n] &= \sum\limits_{k = 0}^M {{b_k}} x[n - k]\\
 &= \sum\limits_{k = 0}^M {{b_k}} \left( {\alpha {x_1}[n - k] + \beta {x_2}\left[ {n - k} \right]} \right)\\
 &= \sum\limits_{k = 0}^M {{b_k}} \left( {\alpha {x_1}[n - k]} \right) + \sum\limits_{k = 0}^M {{b_k}} \left( {\beta {x_2}\left[ {n - k} \right]} \right)\\
 &= \alpha \sum\limits_{k = 0}^M {{b_k}} {x_1}[n - k] + \beta \sum\limits_{k = 0}^M {{b_k}} {x_2}\left[ {n - k} \right]\\
 &= \alpha {y_1}\left[ n \right] + \beta {y_2}\left[ n \right]
\end{align*}由於輸出為 $y_1[n],y_2[n]$ 之線性組合，亦即 $y[n] = \alpha y_1[n] + \beta y_2[n]$ ，故可知此 FIR 濾波器為線性。

接著我們證明非時變性質，對任意 $n_0 \in \mathbb{N}$ 考慮時延輸入訊號 $x[n-n_0]$，則不難發現
\[\sum\limits_{k = 0}^M {{b_k}} x[n - {n_0} - k] = y\left[ {n - {n_0}} \right]\]
亦即等量的時延輸入 導致 等量時延輸出，故此系統為非時變。$\square$

[數學分析] Weierstrass Theorem (1) - 先備概念

回憶在數學分析的內容中，我們試圖利用 $\mathbb{Q}$ 在 $\mathbb{R}$ 中 dense的想法，指出 任意在 $\mathbb{R}$ 上的實數 $r$，皆可透過 一組 sequence $\{q_n\} \in \mathbb{Q}$ 逼近 。也就是說 $q_n \rightarrow r$ 當 $n \rightarrow \infty$。那麼我們想問在函數上是否也有類似的概念? Weierstrass Approximation Theorem 便是試圖回答這個問題。

Weierstrass  Approximation Theorem 主要想法： 利用多項式 均勻收斂 連續函數!!

不過在介紹之前，我們需要一些先備知識。
首先看個 算子 (operator) 的概念：
定義 $A: \text{one function} \rightarrow \text{different function}$ 為一個算子(operator)

我們看個例子：

------------

Example

Fourier transform of 函數 $f$ 為一個算子 (將函數 $f$ 映射到另一個函數 $F$)

\[F(j\omega ) = \int_{ - \infty }^\infty  f (t){e^{ - j\omega }}dt
\]-----------

那麼算子何其多? 哪一種算子適合我們?? 以下我們介紹一個即為有用的特殊算子：摺積(Convolution)

===================

Definition: Convolution (Integral)

給定兩可積函數 $f,g$ on $\mathbb{R}$，則其折積(convolution) 定義為

\[(f*g)\left( x \right): = \int f (x - y)g(y)dy = \int g (x - y)f(y)dy
\]===================

Example

$f,g$ 為在 $[-1,1]$ 上的週期函數，且 $|\delta| <1$ \[f\left( x \right): = \left\{ \begin{array}{l} 1/2\delta ,\begin{array}{*{20}{c}} {}&{} \end{array}x \in \left[ { - \delta ,\delta } \right]\\ 0,\begin{array}{*{20}{c}} {}&{}&{} \end{array}o.w \end{array} \right.
\]試求其 convolution $(f * g)(x)=?$
Proof:
\[\begin{array}{l}
f*g: = \int_{ - 1}^1 {f\left( y \right)g\left( {x - y} \right)dy}  = \int_{ - \delta }^\delta  {\frac{\delta }{2}g\left( {x - y} \right)dy} \\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = \frac{\delta }{2}\int_{ - \delta }^\delta  {g\left( {x - y} \right)dy}. \ \ \ \ \square
\end{array}
\]

Convolution 的好處：可以保留原函數的本身的優點!!

===================
FACT: 令 $K$ 為 compact interval，若 $f \in \mathcal{C}^{\infty}(\mathbb{R})$ (smooth function) 且 $g$ 為 可積函數，則
\[
(f * g)(x) = \int_K f(y)g(x-y)dy
\]亦為在 $K$ 上 smooth function
==================
Proof: omitted.

現在我們看個函數
===================
Definition: A Specific Smooth Function
令 $\phi(x)$ 為 $\mathbb{R}$ 上的 smooth function 且滿足下列條件
在 $(-1,1)$， $\phi>0$ ；且在 $(-1,1)^c$， $\phi=0$ ；另外我們限定此函數  $\phi$ 必須滿足下列積分
\[
\int_{-1}^{1} \phi(t) dt =1
\]===================

有了上述函數，我們可以定義算子 Operator $A$ 如下：
對 $s>0$，定義 $\phi_s$ 為\[{\phi _s}(t): = \frac{1}{s}\phi (\frac{t}{s});\;\;\;\int_{-1}^1 \phi  dt = 1\]且

\[
A_s f(t) :=( \phi_s * f)(t) = \int \phi_s(t) f(x-t) dt
\]其中 $\phi \ge 0$ on $[-1,1]$ 且 $\phi =0$ on $[-1,1]^c$ 。
============
Theorem: Preliminary Lemma for Weierstrass Approximation Theorem
令 $f$ 為 有界 連續 函數 on $\mathbb{R}$ ($f \in \mathcal{C}(\mathbb{R})$) 且 令 $J$ 為任意 compact interval in $\mathbb{R}$。則 當 $s \rightarrow 0$，我們有 $A_s f \rightarrow f$ 均勻收斂 on $J$。
============

Proof:
令 $\varepsilon>0$ ，我們要證明 $A_s f \rightarrow f$ 均勻收斂；(注意到上述定理陳述是當 $s \rightarrow 0$，故令 $n := 1/s$) 故要證 存在 $N>0$ 使得 $n > N$ 我們有
\[
|A_sf(x) - f(x) | < \varepsilon
\]對任意 $x \in J$ 。

觀察
\[\small
\begin{array}{l}
|{A_s}f(x) - f(x)| = \left| {\int_{ - \infty }^\infty  {{\phi _s}\left( t \right)f\left( {x - t} \right)dt}  - f\left( x \right)} \right|\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}
\end{array} = \left| {\int_{ - \infty }^\infty  {{\phi _s}\left( t \right)f\left( {x - t} \right)dt}  - f\left( x \right)\underbrace {\int_{ - \infty }^\infty  {{\phi _s}\left( t \right)dt} }_{ = 1}} \right|\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}
\end{array} = \left| {\int_{ - \infty }^\infty  {\left[ {f\left( {x - t} \right) - f\left( x \right)} \right]{\phi _s}\left( t \right)dt} } \right| \le \int_{ - \infty }^\infty  {\left| {f\left( {x - t} \right) - f\left( x \right)} \right|{\phi _s}\left( t \right)dt}  < \varepsilon
\end{array}\]注意到上式中當 $|t| >s$時， $\phi_s =0$故上式積分範圍為
\[|{A_s}f(x) - f(x)| \le \int_{ - s}^s {\left| {f\left( {x - t} \right) - f\left( x \right)} \right|{\phi _s}\left( t \right)dt}  < \varepsilon \]
故如果可以讓上式 $|f(x-t) - f(x)|$ 對任意 $x \in J$ 都可使其任意小便完成證明。故現在透過 $J$ 為任意 compact interval on $\mathbb{R}$，我們可推知 $f$ 在 $J$ 上為 uniform continuous。亦即我們可選 $\delta >0$ 使得對任意 $u,v \in J$，
\[
|u-v| < \delta \Rightarrow |f(u) - f(v)| < \varepsilon
\]取 $u :=x, v := x+t \in J$ 則我們可知
\[
|u-v| = |x - (x+t)| < \delta \Rightarrow  |f(x) - f(x-t)| < \varepsilon
\]也就是說只要能找到 $N$ 使得 $n >N$ 且滿足 $|t| < \delta$ 則 我們便會有
\[|{A_s}f(x) - f(x)| \le \int_{ - s}^s {\left| {f\left( {x - t} \right) - f\left( x \right)} \right|{\phi _s}\left( t \right)dt}  < \varepsilon \]由於 $n =1/s; \; n >N \Rightarrow s < 1/N$ 且又由 $\phi_s(t)$ 定義可知 $t \in [-s,s]$故取 $N = 1/\delta$ 則上述自動滿足。$\square$