2017年4月30日 星期日

[機率論] 動差生成函數 的 常見應用 (1)

令 $\Omega$ 為 樣本空間,定義 $X : \Omega \to \mathbb{R}$ 為配備 機率分配函數 $f_X$ 的 隨機變數。

==================
Definition: k-th Moment
隨機變數 $X : \Omega \to \mathbb{R}$  的 k階 動差 (k-th moment) 定為 $E[X^k]$
==================


==================
Definition: Moment Generating Function
令 $X$ 為隨機變數,若存在 $\delta >0$ 使得對 $t \in (-\delta,\delta )$ 而言,期望值 $E[e^{tX}]$ 存在,則 $X$ 的 動差生成函數 (Moment Generating Function, mgf) 存在,且定義為
\[
M_X(t) := E[e^{tX}], \;\;\; t \in (-\delta,\delta )
\]除此之外,若上述條件成立則
\[
D_t^k E[e^{tX}] = E[D_t^k e^{tX}]
\]其中 $D_t^k $ 表示為對 $t$ 微分 $k$次 微分算子。
==================

Comments:
1. 動差生成函數 $M_X(t)$ 是一個以 $t$ 為變數 的函數,目的在於 "產生動差",至於如何產生我們會在下面進行討論。
2. 上述定義僅僅要求 mgf 在 $t = 0$ 附近開區間 $t \in (-\delta, \delta)$ 期望值存在,此條件也保證積分與微分互換性。
3. 若在開區間 $ t \in (-\delta, \delta)$ 期望值存在,立刻可得知 $M_X(0) = E[e^{0}] = 1$


FACT: 上述動差生成函數算是非常便利的工具,我們可以透過其產生各種具有常見分配的隨機變數之一,二階動差。假定某隨機變數之 mgf 存在,則此隨機變數的 k-th 動差表為
\[
D_t^k M_X (0) = E[X^k], \;\; k=1,2,...
\]

Comment:
由上述討論可知,若我們想求期望值則
\[
D_t M_X(0) = E[X]
\]若我們想求變異數則
\begin{align*}
Var(X) &= E[(X - E[X] )^2]\\
& = E[X^2] -(E[X])^2\\
& = D_t^2 M_X(0) - D_t M_X(0)
\end{align*}

==================
Theorem: MGF唯一決定分配函數且 兩 MGF 相等 保證 分配函數相等 
令 $X,Y$ 為兩隨機變數且其各自對應的 mgf $M_X, M_Y$ 在含 $0$ 開區間存在。則此兩隨機變數的分配函數 $$
F_X(z) = F_Y(z), \;\;\; \forall \;\; z \in \mathbb{R}
$$若且唯若 存在 $\delta >0$ 使得 對任意 $t \in (-\delta, \delta)$ 而言
$$M_X(t) = M_Y(t)$$
==================
Proof: omitted


Comments:
1. 上述定理是非常強大的結果,他說明了要決定某隨機變數的分配可以透過 mgf 來唯一決定。也就是說 mgf 與 分配函數為 1-1 對應,且如果當分配函數很困難取得,我們可以透過計算 mgf 來幫助我們確定分配。
2. 機率論中另外有一種函數稱作 特性函數 (characteristic function) 其作用與 mgf相仿,且永遠存在,定義為 $\varphi_X(t) := E[e^{i tX}]$,但是此式包含複數,一般在計算上會稍微比 mgf 複雜一些。


以下我們展示幾個例子來使用動差生成函數求得 期望值 與 變異數。首先我們看幾個離散隨機變數的例子:

==================
Example 1: Bernoulli Random Variable 的 期望值 與 變異數
令 $X$ 為 Bernoulli random variable 滿足 $P(X=1) = p$ 且 $P(X=0)=1-p$。
(a) 試問其 mgf 是否存在?若存在則求其 mgf
(b) 利用 part(a) 證明 $E[X] = p$ 與 $Var(X) = p (1-p)$。
==================

Proof (a): 以下我們使用 mgf 方法來求期望值 以及 變異數 ,首先計算 mgf ,由定義可知
\begin{align*}
  {M_X}(t) &= E\left[ {{e^{tX}}} \right] \hfill \\
   &= \sum\limits_{x \in \left\{ {0,1} \right\}}^{} {{e^{tx}}{f_X}\left( x \right)}  \hfill \\
   &= {e^{t1}}p + {e^{t0}}\left( {1 - p} \right) \hfill \\
   &= {e^t}p + \left( {1 - p} \right) \hfill \\
\end{align*} 注意到上式對任意 $t \in \mathbb{R}$ 成立,故存在 $\delta >0$ 使得對 $t \in (-\delta,\delta )$ 而言,期望值 $E[e^{tX}]$ 存在。

Proof (b): 
以下我們計算 mgf 的微分
\[\left\{ \begin{align*}
  {D_t}{M_X}(t) &= {D_t}\left( {{e^t}p + \left( {1 - p} \right)} \right) = {e^t}p \hfill \\
  D_t^2{M_X}(t) &= {e^t}p \hfill \\
\end{align*}  \right.\] 故 $ E[X] = {D_t}{M_X}(0) = p $ 且
\begin{align*}
 Var(X) &= D_t^2{M_X}(0) - {({D_t}{M_X}(0))^2} \hfill \\
 &= p - {p^2} = p\left( {1 - p} \right) \hfill \\
 \end{align*}至此得證 $\square$





==================
Example 2: Poisson Random Variable 的 期望值 與 變異數
令 $X \sim Poisson(\lambda)$ 為 Poisson random variable  配備機率密度函數
\[
f_X(x) = \frac{\lambda^x e^{-\lambda}}{x!},\;\;\; x=0,1,2,...
\]
(a) 試問其 mgf 是否存在?若存在則求其 mgf
(b) 利用 part(a) 求 $E[X] = \lambda$ 與 $Var(X) = \lambda$。
==================

Proof (a):
由 mgr 定義,我們觀察 \begin{align*}
  {M_X}\left( t \right) &= E\left[ {{e^{tX}}} \right] \hfill \\
   &= \sum\limits_{x = 0}^\infty  {{e^{tx}}{f_X}(x)}  \hfill \\
   &= \sum\limits_{x = 0}^\infty  {{e^{tx}}\frac{{{\lambda ^x}{e^{ - \lambda }}}}{{x!}}}  \hfill \\
   &= {e^{ - \lambda }}\underbrace {\sum\limits_{x = 0}^\infty  {\frac{{{{\left( {{e^t}\lambda } \right)}^x}}}{{x!}}} }_{ = {e^{{e^t}\lambda }}} = {e^{\left( {{e^t} - 1} \right)\lambda }} \hfill \\
\end{align*}
上式最後一條等式成立因為
\[{e^z}: = \sum\limits_{k = 0}^\infty  {\frac{{{z^k}}}{{k!}}},\;\;\; \forall \; z \in \mathbb{C} \]注意到 $M_X(t)$ 對任意 $t \in \mathbb{R}$ 有定義,故存在 $\delta >0$ 使得對 $t \in (-\delta,\delta )$ 而言,期望值 $E[e^{tX}]$ 存在。

Proof (b):
首先對 mgf 求一階 與 二階導數,
\[\left\{ \begin{gathered}
  {D_t}{M_X}\left( t \right) = {D_t}\left( {{e^{\left( {{e^t} - 1} \right)\lambda }}} \right) = {e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t} \hfill \\
  D_t^2{M_X}\left( t \right) = D_t^2\left( {{e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t}} \right) = {e^{\left( {{e^t} - 1} \right)\lambda }}{\left( {\lambda {e^t}} \right)^2} + {e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t} \hfill \\
\end{gathered}  \right.\] 故期望值為 $E\left[ X \right] = {D_t}{M_X}\left( 0 \right) = \lambda $ 且變異數為
\begin{align*}
  Var\left( X \right) &= D_t^2{M_X}\left( 0 \right) - {\left( {{D_t}{M_X}\left( 0 \right)} \right)^2} \hfill \\
   &= {\left. {\left( {{e^{\left( {{e^t} - 1} \right)\lambda }}{{\left( {\lambda {e^t}} \right)}^2} + {e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t}} \right)} \right|_{t = 0}} - {\lambda ^2} = \lambda  \hfill \\
\end{align*} 至此得證 $\square$






下面例子是 連續隨機變數 的情況。

=======================
Example 3: Normal Distribution
令 $X \sim \mathcal{N}(\mu, \sigma^2)$ 配備機率密度函數
\[{f_X}(x): = \frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( { - \frac{{{{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}} \right), \;\;\; x \in \mathbb{R}
\]其中 $\mu \in \mathbb{R}, \sigma >0$
(a) 試求 $X$ 的 mgf
(b) 利用 part(a) 證明 $E[X] =\mu$ 與 $Var(X) =\sigma^2$
=======================

Proof: (a)
由定義出發,觀察
\begin{align*}
  {M_X}\left( t \right) &= E\left[ {{e^{tX}}} \right] \hfill \\
   &= \int_{ - \infty }^\infty  {{e^{tx}}{f_X}(x)} dx \hfill \\
   &= \int_{ - \infty }^\infty  {{e^{tx}}\frac{1}{{\sqrt {2\pi } \sigma }}{e^{ - \frac{{{{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= \frac{1}{{\sqrt {2\pi } \sigma }}\int_{ - \infty }^\infty  {{e^{\frac{{2{\sigma ^2}tx - {{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= \frac{1}{{\sqrt {2\pi } \sigma }}{e^{\frac{{ - {\mu ^2}}}{{2{\sigma ^2}}}}}\int_{ - \infty }^\infty  {{e^{\frac{{ - \left( {{x^2} - 2\left( {{\sigma ^2}t + \mu } \right)x + {{\left( {{\sigma ^2}t + \mu } \right)}^2}} \right)}}{{2{\sigma ^2}}}}}{e^{\frac{{{{\left( {{\sigma ^2}t + \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= \frac{1}{{\sqrt {2\pi } \sigma }}{e^{\frac{{ - {\mu ^2}}}{{2{\sigma ^2}}}}}{e^{\frac{{{{\left( {{\sigma ^2}t + \mu } \right)}^2}}}{{2{\sigma ^2}}}}}\int_{ - \infty }^\infty  {{e^{\frac{{ - {{\left( {x - \left( {{\sigma ^2}t + \mu } \right)} \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= {e^{\frac{{{\sigma ^2}{t^2} + 2\mu t}}{2}}}\underbrace {\frac{1}{{\sqrt {2\pi } \sigma }}\int_{ - \infty }^\infty  {{e^{\frac{{ - {{\left( {x - \left( {{\sigma ^2}t + \mu } \right)} \right)}^2}}}{{2{\sigma ^2}}}}}} dx}_{ = 1{\text{ }}\left( {{\text{pdf of normal }}\mathcal{N}\left( {\mu+\sigma^2 t ,{\sigma ^2}} \right)} \right)} \hfill \\
   &= {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}} \hfill \\
\end{align*}

Proof (b):
以下我們計算 mgf 的微分
\[\left\{ \begin{align*}
  {D_t}{M_X}\left( t \right) &= {D_t}\left( {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}} \right) = {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}\left( {\mu  + {\sigma ^2}t} \right) \hfill \\
  D_t^2{M_X}\left( t \right) &= {D_t}\left( {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}\left( {\mu  + {\sigma ^2}t} \right)} \right) = {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{\left( {\mu  + {\sigma ^2}t} \right)^2} + {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{\sigma ^2} \hfill \\
\end{align*}  \right.\]故對應的期望值為
\[E\left[ X \right] = {D_t}{M_X}\left( 0 \right) = {\left. {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}\left( {\mu  + {\sigma ^2}t} \right)} \right|_{t = 0}} = \mu \]
且 變異數 為
\begin{align*}
  Var\left( X \right) &= D_t^2{M_X}\left( 0 \right) - {\left( {D_t^1{M_X}\left( 0 \right)} \right)^2} \hfill \\
   &= {\left. {\left( {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{{\left( {\mu  + {\sigma ^2}t} \right)}^2} + {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{\sigma ^2}} \right)} \right|_{t = 0}} - {\mu ^2} \hfill \\
   &= {\sigma ^2} \hfill \\
\end{align*} 至此得證 $\square$

2017年4月26日 星期三

[訊號處理] 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 中,上述討論可以使用 fft2ifft2 來實現,在此不做贅述。


以下我們探討幾種常見的影像濾波器,亦即幾種常見的 $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 何者運算較快?影像太大該怎麼壓縮?等等。




2017年4月20日 星期四

[凸分析] 半正定對稱矩陣所成之集合為凸錐

首先定義 $S^n$ 為由所有 實係數對稱矩陣 所形成之集合,表為
\[
S^n := \{A \in \mathbb{R}^{n \times n} : A^T = A\}
\] 則不難證明 $S^n$ 為一個 subspace (why?),故 $S%n$ 必為 vector space 。

Comments:
由於  $S^n$ 為一個 subspace ,我們可以定義其維度,且值得一提的是 $\dim S^n = (n) (n+1)/2$,在此不做贅述。


現在我們收集所有實係數對稱 且 半正定 (positive semidefinite) 矩陣,定義
\[
S_+^n := \{A \in S^n: A  \succeq  0\}
\] 其中 $A \succeq 0$ 表示 $A$ 為 半正定矩陣:亦即給定任意 $x \in \mathbb{R}^n$ 我們有
\[
x^T A x \geq 0
\]

則我們聲稱上述 $S_+^n $為 凸錐 (convex cone),以下我們給出主要 FACT :
===================
FACT:
$S_+^n$ 為 convex cone。
===================

Proof:
要證明 $S_+^n$ 為 convex cone,我們要證明 $S_+^n$ 為一個 cone 且 $S_+^n$ 為 convex,故我們令 $A,B \in S_+^n$ 與 $\theta_1, \theta_2 \geq 0$ 且必須證明
\[
\theta_1 A + \theta_2 B \in S_+^n
\]
此等價證明 $ \theta_1 A + \theta_2 B $ 為對稱矩陣 且 半正定。現在我們首先證明對稱性:觀察
\[{({\theta _1}A + {\theta _2}B)^T} = {\theta _1}{A^T} + {\theta _2}{B^T} = {\theta _1}A + {\theta _2}B\]注意最後一條等式成立因為  $A,B \in S_+^n$ 故  $A=A^T,B=B^T$。

接著我們證明半正定性質:取 $x \in \mathbb{R}^n$ 觀察
\[{x^T}({\theta _1}A + {\theta _2}B)x = {\theta _1}{x^T}Ax + {\theta _2}{x^T}Bx \geqslant 0\]同樣最後一條等式成立因為   $A,B \in S_+^n$ 故  $x^TAx \geq 0$ 且 $x^T Bx \geq 0$。$\square$


Comments:
上述結果指出 線性代數中的 對稱半正定矩陣 可以與 凸分析 中的凸集拉上關係。

2017年4月9日 星期日

[凸分析] 仿射組合 仿射空間 與 線性方程關係

考慮 中學數學 中提及的 直線方程 (更嚴格的說法是 affine function 在此用 斜截式 表示):令 $x \in \mathbb{R}^1$,定義函數 $y: \mathbb{R}^1 \to \mathbb{R}^1$ 滿足
\[
y(x) = m x + b
\] 其中 $m$ 表示斜率, $b$ 表示截距。現在我們進一步觀察上式並將其改寫如下:
\[
y(x) = (m + b - b) x  + b
\]則讀者不難發現可得 $y(x) = (m + b) x  + b (1 - x) $ 現在若令 $a := m+b$ 則我們得到如下簡潔的形式
\[
y(x) = a x + b (1-x)
\]

Comments:
注意到 $ y(x):=y = a x + b (1-x)$ 一般稱 $y$ 為透過 $x, (1-x)$ 所成之 線性組合 (linear combination),若 $0 \le x \le 1$,則上式一般稱為 $a$ 與 $b$ 的 凸組合 (convex combination)


推廣到有限維度歐式空間:
上述結果可以推廣到 $\mathbb{R}^n$ 空間:考慮 $x_1 \neq x_2$ 為 $\mathbb{R}^n$ 中的兩(向量)點,則
\[
y := \theta x_1 + (1 - \theta) x_2  \;\;\;\; (*)
\] 其中 $\theta \in \mathbb{R}$ 形成  $\mathbb{R}^n$ 過點 $x_1$ 與 $x_2$ 之直線。讀者可觀察若 $\theta = 0$ 則 $y=x_2$。若 $\theta = 1$ 則 $y= x_1$。亦即當我們調整參數 $\theta \in [0,1]$ 可得到一條 $x_1$ 與 $x_2$ 的封閉線段 (line segment)。另外我們亦可將 $(*)$ 改寫如下
\begin{align*}
  &y = \theta {x_1} + (1 - \theta ){x_2} \hfill \\
   &\Rightarrow y = {x_2} + \theta \left( {{x_1} - {x_2}} \right) \hfill \\
\end{align*} 則此時我們可以用另一種觀點來看上述直線方程:亦即上述直線方程有 "基準點" $x_2$ (對應 $\theta = 0$) 與 透過 參數 $\theta$ 調整後的 "方向" $x_1 - x_2$ (從 $x_2$ 指向 $x_1$ )。有了上述觀念,我們可以進一步提出 仿射集(Affine Set) 的概念:


=====================
Definiton:  仿射集 (Affine Set)
我們說一個集合 $C \subset \mathbb{R}^n$ 為 affine 若下列條件成立:對任意兩點 $x_1, x_2 \in C$ 與 $\theta \in \mathbb{R}$ ,其兩點用參數 $\theta$ 所成之線段 滿足
\[
\theta x_1 + (1-\theta)x_2 \in C
\]=====================

Comments:
注意到上述定義要求 $x_1,x_2$ 的線性組合之係數和為 $\theta + (1 - \theta) = 1$。

事實上上述定義不必僅僅取兩點,我們可以取任意有限多點比如 $x_1,x_2,...,x_k$ 且我們建構
\[
\theta_1 x_1 + \theta_2 x_2 + ... + \theta_k x_k
\]其中 $\sum_{i=1}^k \theta_i = 1$。這種有額外要求  $x_1,x_2,...,x_k$  係數和 為 $1$ 的特殊線性組合又稱作 $x_1,x_2,...,x_k$ 的仿射組合 (affine combination)

=============
FACT: 令 $C$ 為 affine 且任取一點 $x_0 \in C$ ,則 集合
\[
V:= C - x_0 := \{x-x_0 : x \in C\}
\]為子空間 subspace (亦即滿足 向量加法封閉性 與 純量乘法封閉性)。換言之若 $V$ 為 subspace 且 $x_0$ 任取為 $C$ 中一點,則集合
\[
C = V + x_0
\]為 affine。
=============

Comments: 
關於(實數)子空間更嚴格的定義如下:我們說非空集合 $V$ 為 子空間 若且唯若 $V$ 滿足向量加法封閉性 與 純量乘法封閉性:
1. 向量加法封閉性:對任意 $v_1, v_2 \in V,$ $v_1+v_2 \in V$
2. 純量乘法封閉性:對任意 $v \in V$ 與 $c \in \mathbb{R}^1$,$c v \in V$。


FACT: 線性方程之解所成的集合為仿射集
事實上 仿射集合 離我們並不遙遠,比如說考慮 任意線性方程的解所成之集合
\[
C:= \{x\in \mathbb{R}^n: Ax = b\}
\]其中 $A \in \mathbb{R}^{m \times n}$ 與 $b \in \mathbb{R}^m$ 則此集合即為仿射集。

Proof : 要證明 $C$ 為 affine ,我們從定義出發:取 $x,y \in C$ 與 $\theta \in \mathbb{R}^1$ 我們要證明
\[
\theta x + (1-\theta) y \in C \;\;\;\; (**)
\]注意到  $x,y \in C$ ,故 $Ax = b$ 且 $Ay=b$,要證明 $(**)$成立,則等價證明
\[
A (\theta x + (1-\theta) y) = b
\]上述等式成立因為:
\begin{align*}
  A(\theta x + (1 - \theta )y) &= \theta Ax + (1 - \theta )Ay \hfill \\
   &= \theta b + (1 - \theta )b = b \hfill \\
\end{align*} 故此得證。

以下我們給出一些常見的 affine set 例子

Example:
1. 任意空集合 $\emptyset$ 為 affine
2. 任意單點集 $\{x\}$ 為 affine
3. 任意 subspace 為 affine

2017年4月6日 星期四

[訊號與系統] FIR系統的弦波響應

考慮 FIR 系統為 線性非時變(Linear Time-Invariant, LTI) 系統,且假設輸入為 離散 complex exponential 則其對應的輸出將非常容易計算:考慮 FIR 系統
\[
y[n] = \sum_{k=0}^M b_k x[n-k]
\]

假設輸入為 complex exponential 表為 $x[n] = x(nT_s) = A e^{j \varphi} e^{j \omega Ts n}$ 且 $-\infty < n < \infty$。則輸出為
\[\begin{array}{l}
y[n] = \sum\limits_{k = 0}^M {{b_k}} x[n - k]\\
 = \sum\limits_{k = 0}^M {{b_k}} A{e^{j\varphi }}{e^{j\omega Ts\left( {n - k} \right)}}\\
 = \left( {\sum\limits_{k = 0}^M {{b_k}} {e^{j\omega Ts\left( { - k} \right)}}} \right)A{e^{j\varphi }}{e^{j\omega Ts\left( n \right)}}\\
 := H\left( {{e^{ - j\widehat \omega }}} \right)\underbrace {A{e^{j\varphi }}{e^{j\omega Ts\left( n \right)}}}_{ = x\left[ n \right]}
\end{array}\]其中 $\widehat{\omega} := \omega T_s$ 且
\[H\left( {{e^{ - j\widehat \omega }}} \right) = \sum\limits_{k = 0}^M {{b_k}} {e^{ - j\widehat \omega k}} = \sum\limits_{k = 0}^M {h\left[ k \right]} {e^{ - j\widehat \omega k}}\]稱作 frequency-response function,一般而言我們簡稱為 frequency response。

回憶由於 FIR系統的脈衝響應 與 濾波器的係數相同,亦即 $b_k = h[k]$,我們可以將上述頻率響應改寫成
\[H\left( {{e^{ - j\widehat \omega k}}} \right) = \sum\limits_{k = 0}^M {{b_k}} {e^{ - j\widehat \omega k}} = \sum\limits_{k = 0}^M {h\left[ k \right]} {e^{ - j\widehat \omega k}}\]

Comments:
1. 輸入為離散時間訊號 $x[n] = x(nT_s) = A e^{j \varphi} e^{j \omega Ts n}$ 則 LTI FIR 系統亦為離散時間 complex exponential 乘上不同的複數大小,但頻率 $\widehat{\omega}$ 不變。

2. 上述以離散時間訊號為 complex exponential $x[n] = x(nT_s) = A e^{j \varphi} e^{j \omega Ts n}$ 為輸入,則 LTI FIR 系統可表為
\[y[n] = H\left( {{e^{ - j\widehat \omega k}}} \right)x\left[ n \right]\]

3. 注意到 $H\left( {{e^{ - j\widehat \omega k}}} \right)$ 為複數值函數,故我們可使用極座標表示法將其表示為
\[y[n] = \left( {\left| {H\left( {{e^{ - j\widehat \omega k}}} \right)} \right|{e^{j\angle H\left( {{e^{ - j\widehat \omega k}}} \right)}}} \right)A{e^{j\varphi }}{e^{j\omega Ts\left( n \right)}}\]
亦即,
\[\begin{array}{l}
y[n] = H\left( {{e^{ - j\widehat \omega k}}} \right)\underbrace {A{e^{j\varphi }}{e^{j\omega Ts\left( n \right)}}}_{ = x\left[ n \right]}\\
 = \left( {\left| {H\left( {{e^{ - j\widehat \omega k}}} \right)} \right|{e^{j\angle H\left( {{e^{ - j\widehat \omega k}}} \right)}}} \right)A{e^{j\varphi }}{e^{j\widehat \omega n}}\\
 = \left( {\left| {H\left( {{e^{ - j\widehat \omega k}}} \right)} \right|A} \right){e^{j\left( {\angle H\left( {{e^{ - j\widehat \omega k}}} \right) + \varphi } \right)}}{e^{j\widehat \omega n}}
\end{array}\]一般而言,我們會稱 ${\left| {H\left( {{e^{ - j\widehat \omega k}}} \right)} \right|}$ 為系統增益 (gain) 且 ${\angle H\left( {{e^{ - j\widehat \omega k}}} \right)}$ 稱為系統相位 (phase)。

[訊號與系統] 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$