棘波偵測之研究背景與演算法則 - 以回饋式自動模板生成為基礎之正規化關聯值棘波偵測系統之設計及實現

第二章棘波偵測之研究背景與演算法則

2-1 棘波分類的介紹

傳統的棘波分類演算法主要包含了棘波偵測 (Spike Detection)、特徵擷取

(Feature Extraction) 與棘波分類 (Clustering) 三大步驟，詳見於圖 2-1。棘波分類演算法會先將輸入的棘波序列 (Spike Train) 經由棘波偵測階段來偵測出棘波序

Extraction Clustering

2-2 棘波偵測演算法則

2-2-1 Matched filter

我們先從最基本型的Matched filter technique 討論起。Matched filter 是用來進行棘波偵測的一種技術，其核心的技術是利用多組事先儲存好的模板(template) 拿來循序的與棘波序列(Spike train)的區段(𝑥[m])做摺積運算。大致上的流程如下

圖所示:

m=0 𝑥=0

dot(𝑥 ,template)

>Threshold

HIT miss

m=m+1

Y N

圖 2-2 Matched filter 的運作流程圖

為了方便解釋，我們先假設Matched filter目前只存在一組模板。以下我們對一些符號做定義; 𝑥[m]為棘波序列中第m個樣本(sample)。𝑥 = [x[𝑚], x[𝑚 − 1], … , x[𝑚 − N + 1]]^𝑇為棘波序列中第m個區段(segment)，其中N為區段長度。而

用於 Matched filter 的模板同樣擁有 N 個元素 (element)，我們將其定義為：

𝐭 = [t[1], t[2], … , t[N − 1]]^𝑇。y[m]定義為 Matched filter 在棘波序列中第 m 點的 摺積輸出值。

𝑦[𝑚] = ∑ 𝑥[𝑚 − 𝑘]𝑡[𝑘] = 𝐱

^𝑇

𝐭 (2 − 1)

𝑁 𝑘=

值得注意的是，摺積的計算相當於區段 𝐱 與模板

t

做內積的計算，這也表示摺積算出來的值代表區段 𝐱 與模板

t 這兩向量之間的相關性。當

y[m] 比事先訂定的閥值η 大時(y[m]> η)，區段𝐱 則被判斷為棘波。然而 Matched filter 的缺點是單一的 η 值無法在不同 SNR 環境下做出有效的偵測。

2-2-2 Block Energy

接下來介紹 Block Energy 演算法，這個演算法在進行棘波偵測的時候，也會將輸入進來的棘波序列先分成數個區段，然後再觀察每個區段中的採樣點彼此平方相加的值，若值超過了事先訂定好的閥值，則該值相對應的區段就會被判斷為棘波。流程如下圖:

圖 2-3 Block energy 的運作流程圖

假設𝑥[m]為棘波序列中第 m 個樣本(sample)。𝑥 = [x[𝑚], x[𝑚 − 1], … , x[𝑚 − N + 1]]^𝑇為棘波序列中第 m 個區段(segment)，其中 N 為區段長度。y[m]定義為 Block energy 在棘波序列中第 m 點的能量輸出值。

𝑦[𝑚] = ∑ 𝑥[𝑚 − 𝑘]

=

^𝑥_𝑚^𝑇^x_𝑚

(2 − 2)

𝑁 𝑘=

從Block energy的流程圖跟公式可以看的出來，雖然在偵測的觀念上略有不同，但是跟前面提到的Matched filter 的運作卻頗為相似，Block energy只是把

Matched filter 中模板的部分換成了自身的區段x[m]。 Block Energy的優點在於計算複雜度低，在雜訊等級低的時候棘波偵測效果好。缺點是當輸入訊號的雜訊太高的時候，因為區段內的能量會被雜訊影響而提高，所以就不好分辨出該區段是屬於棘波或是雜訊，導致偵測效果惡化。

2-2-3 LRT and GLRT tests

假設模板 t 是已知的棘波，而且輸入系統的棘波序列中的雜訊 n 是屬於 additive zero-mean white Gaussian noise 的話，則對區段𝐱 和模板 t 做 Match filter

(𝐱^𝑇t ≥ η) 的意義等於對區段𝐱 做 Likelihood ratio test (LRT)。所以我們現在觀察

設 c_ij 是當假說H_𝑗為真，但我們卻接受假說H_𝑖 的機率， P(H_𝑖) 是x 屬於在假說H_𝑖的機率。根據 optimal test[22]，我們設 log^𝑃(𝐻_𝑃(𝐻⁰^)(𝑐¹⁰^𝑐⁰⁰⁾

1)(𝑐₀₁ 𝑐₁₁) 為檢定統計量，

並給予符號γ，並檢驗之，若 𝑥^𝑇𝛴𝑡 > γ 則x 符合H 假說，反之則符合H₀假說。觀察(式 3)，因為Σ = σ²𝐈，故我們可以把其代換成:

觀察(式 2-3)與(式 2-1)，可以發現若把閥值η定為𝜎²𝛾，則對x 做 LRT test 等價於對x 做 matched filter。

若我們一開始能知道棘波的資訊，則對於 matched filter 演算法的準確性就能大幅的提升，但是這在實際的過程中可能無法實現，我們從欲偵測的區段x 中 先估計出一個模板 t 來解決這個問題。假設x 是符合假說H，我們使用 maximum likelihood (ML) 來估計出模板 t，當模板𝐭 = x 時，模板 t 使得P(𝑥 |𝐻 )的機率最大，故我們選擇此情況，並新定義一個符號𝐭̂。

𝐭̂ = x = maxP(x |H ) (4) x^𝑇𝑡

H H₀

< 𝜎²𝛾 ≡ η (2-3)

總而言之，GLRT 是一種 block energy detector。假設x 的平均值是 0，我們可以用以下方法估計出區段x 的標準差σ[10]:

2-2-4 Normalized Correlator

回顧上述提到的棘波偵測演算法，Matched filter 和 LRT test 其實都有一個很大的缺點，就是制定理想的閥值是一個很大的挑戰。觀察(式 3)可知，如果我們沒有辦法知道𝑃(𝐻₀)、𝑃(𝐻 )和𝜎²的資訊，就很難估算出理想的η。

因此，我們選擇另外一種方式制定閥值，透過量測模板 t 跟區段

𝐱

^的誤差

上限來決定η。我們事先訂定其誤差上限的值，當輸入的區段

𝐱

^{小於誤差上限，}

就判斷為棘波。這樣一種 correlator 的設計方式，可以對應到另一個僅用平方誤差的簡單比對技術。

我們利用平方誤差 (squared distance) 來觀察區段 𝐱 與模板

t

之間的關係，我們首先可以觀察下面所列出的平方誤差 (squared distance)公式：

d(𝐱 , t) = ‖ 𝐱 ‖

+ ‖ 𝐭 ‖

- 2𝐱

t. (2-4)

因此, 當 𝐱^𝑇

t

> η 時,

d(𝐱 , t) ≤ ‖ 𝐱 ‖

+ ‖ 𝐭 ‖

- 2 η. (2-5)

觀察上式得知，當偵測到 𝐱 為棘波時 (i.e., 𝐱^𝑇t > η)， d(𝐱 , t)的上限值取決於

‖ 𝐱 ‖²、‖ 𝐭 ‖²及η，而‖𝐱 ‖²大小取決於輸入棘波序列。所以當 ‖ 𝐱 ‖² 變大時，

d(𝐱 , t)的值就會跟著變大，即使 𝐱^𝑇t > η。在這種情況下，就可能會把雜訊判斷 成棘波，而產生 false alarm 的情形。

規化(Normalization)的動作。接著我們將做完正規化的 𝐱 及 t 分別定義為 𝐱̅

與 𝐭̅ ，如下

𝐱 =

_‖𝐱^𝐱^𝑚

𝑚‖

, 𝐭̅ =

_‖𝐭‖^𝐭

. (8)

而 Normalized correlator 之輸出為

𝑦̅[𝑚] = ∑ 𝑥̅ [𝑘]𝑡[𝑘] = 𝐱

^𝑇

𝐭̅ (2 − 6)

𝑁 𝑘=

其中 𝑥̅ [k] 及 𝑡[𝑘] 分別為 𝐱 及 𝐭̅ 之第k個element，

因此，(式 2-5)可以改寫成

d(𝐱̅ , 𝐭̅) = 2 - 2𝐱̅

^𝑇

𝐭̅. (9)

因為 d(𝐱_𝑛, t) ≥ 0, 我們可以得到

𝐱̅

^𝑇

𝐭̅ ≤ 1. (10)

所以我們的 Normalized correlator 是指 𝐱̅ 與 𝐭̅ 上的計算。當 𝐱̅^𝑇𝐭̅ ≥ η 時，我 們就可以說此時對應到的 𝐱 就是我們偵測到的棘波。從(式 10) 我們可以推出

η ≤ 1. (11)

另外，當 𝐱^𝑇t > η ，根據(式 9) 我們可以得到

d(x̅ , 𝐭̅) ≤ 2(1 − 𝜂) ≡ 𝐷

(12)

D可以表示為HIT時的平方誤差上界，經過Normalized correlator後得到的唯一閥值，可以當成使用在template matching時的平方誤差上界。此外d(x̅ , 𝐭̅)的值取決於閥值η的大小，當閥值η變大時d(x̅ , 𝐭̅)會變小，而閥值η的最大值是1，且

並不會受到輸入的棘波序列所影響。

事實上，用來判斷輸入的棘波區段x_m跟模板t是否匹配的平方誤差最大值為 D = 2(1 − 𝜂)，當閥值η=1時。Normalized correlator這方法對於閥值η提供了一種新的解釋。如果我們設定閥值η=1，也就是當區段x_m跟模板t 做運算後其值超過 η(x̅^𝑇𝐭̅ ≥ 1)才算是偵測到棘波的話，就代表該區段x_m跟模板t有完全相關性(full

correlation)，這時他們之間的平方誤差為0 。換句話說，如果當閥值η=0.5，就表示了區段x_m跟模板t 有半相關性(half correlation)，只要區段x_m跟模板t運算後的值大於η(x̅^𝑇𝐭̅ ≥ 0.5)就算偵測到棘波，此時它們之間的平方誤差上限值為1。然而如果當閥值η=0，就表示了區段x_m跟模板t 毫無相關性(no correlation)，代表只要有區段x_m進來就算偵測到棘波，此時它們之間的平方誤差上限值為2。

不過在真實情況下因為有雜訊的干擾，我們無法使用完全相關性(full

correlation)當偵測標準(η=1)，根據我們實驗的結果，即使在雜訊很高的環境下我 們只需要 70% 的相關性 (i.e. η = 0.7)，就可以讓我們 Normalized correlator 偵 測達到高True Positive,低 False Positive的效果。

2-2-5 Advance Normalized Correlator

上述的 Normalized Correlator 若是直接實現它的話，會花費非常多的計算時間，所以我們提出了一個改進方法去節省它的運算時間。整個系統最花時間的地方是對輸入的棘波序列的區段(𝐱 )做正規化的動作(𝐱 = ^𝐱^𝑚

‖𝐱𝑚‖)。因為一個區段 𝐱 中含有 N 個樣本，所以我們為了計算𝐱 ，需要 N 個乘法、N 個除法、N-1 個加法、還有一個平方根運算。然後計算完的𝐱 還要跟𝐭̅做內積(𝐱̅^𝑇𝐭̅)，所以還需要 N 個乘法跟 N-1 個加法。根據上面分析的結果，一次完整的 Normalized Correlator 總共需要 2N 個乘法、(2N-2)個加法、N 個除法、1 個平方根運算。

我們為了加速 Normalized Correlator 的計算速度，使用了 Fast energy

computation 和 Post-correlation normalization 這兩種技術。Fast energy computation 主要的運作原理是利用𝐱跟𝐱 這兩者的相關性。

觀察下式:

‖𝑥 ‖

= ‖𝑥

‖

− 𝑥

[𝑚 − 𝑁] + 𝑥

[𝑚] (13)

因為‖𝑥‖²在運算𝐱̅^𝑇 𝐭̅時就已算過了，所以當我們想得到‖𝑥 ‖²時，僅需要再使用 2 個乘法、2 個加法和 1 個平方根運算。

Post-correlation normalization 技術主要透過改變了計算𝐱̅^𝑇𝐭̅的順序，減少運算次數，進而增快運算時間。方法如下式:

𝐱̅

^𝑇

𝐭̅ =

_‖𝑥

𝑚‖

𝑥

^𝑇

𝐭̅ (14)

再對它做正規化。如此一來，在這個階段我們就只需要 1 個除法(因為𝑥^𝑇𝐭̅是純量)。

透過以上兩個優化 Normalized Correlator 的方式，我們的 Advance Normalized Correlator 總共只要花費 N+2 個乘法、N+1 個加法、1 個除法、1 個平方根計算。

(表 2-1)中比較了各種模板匹配方法的計算複雜度，Basic normalized correlator 代表了加入正規化運算的 Matched filter 演算法，而 Fast normalized correlator 代表了利用 Post-correlation normalization 和 Fast energy computation 技術的 Matched filter。

從表上可以觀察到，傳統的 Matched filter 在加上正規化運算之後所需要的加法跟乘法運算都上升了 2 倍，但是跟 Fast normalized correlator 比較後可以發現後者在計算兩個 N 維向量之間的相關性上，運算量只需要花費前者一半的加法與減法，而除法也從 N 次降至 1 次。這證明了我們提出的加速流程應該可以節省近一半的運算時間。

但是我們也觀察到 Fast normalized correlator 計算量就算經過優化了，也還是比傳統的 Matched filter 略高一點。為了克服這個問題，我們再進一步的加入 Pre-Screening 的技術，使它可以更接近 Matched filter 的效能。Pre-Screening 的技術核心是在比較‖𝑥 ‖²跟‖𝐭‖²的大小，假如‖𝑥 ‖²<λ‖𝐭‖²， 0 < λ < 1，‖𝑥 ‖²就會先因為能量過低而被判斷成不是棘波。這個方法可以有效地降低正規化的次數。

Addition Multiplication Division Squared Root

從實驗結果可以發現，結合了 Pre-Screening 的 Fast normalized correlator 可以再更一步的增加其運算速度，甚至比最基本型態的 Matched filter 更有效率。

令 p (0<p<1) 是‖𝑥 ‖²<η‖𝐭‖²成立的機率，如(表 2-1)所示，Pre-Screening 的 Fast normalized correlator 的平均的加法數量為 p(N-1)+2、乘法除法數量為 p*N+2、

除法的數量是 p 次。

(圖 2-4)是棘波偵測部分的流程圖，流程圖裡詳細的表達了在一個棘波序列進來後，我們先利用 Fast energy computation 簡化了獲得‖𝑥 ‖²的運算，再利用 Pre-Screening 過濾掉能量太低的區段，減少 correlation 的次數，最後再利用 Post-correlation normalization 降低內積的運算成本。

圖 2-4 快速計算的流程圖

2.2.6 Osort Algorithm

Osort Algorithm 是一種非監督式的模板分群演算法，這個演算法的特色是不需要經過提取棘波的特徵就能進行分類，而且棘波的群數會由演算法自己增減。

Osort Algorithm 的運作動作如下:

設 s 是目前正要進行分群的棘波。令C_𝑖, 𝑖 = 1, … , c，是代表 Osort algorithm 目前的棘波群集編號， c 是棘波群集的個數。每個群集的平均值記為

𝐭_𝑖, 𝑖 = 1, … , c，也就是說C_𝑖這個群集中所有的棘波平均值為𝐭_𝑖，Osort 演算法一開 始會先計算目前的棘波 s 跟目前所有棘波的平均值𝐭_𝑖的距離，記為𝑑_𝑖 =

𝑑(𝐬, 𝐭_𝑖), i = 1, … , 𝑡。

接著我們找出其中最小的𝑑_𝑖相對的群集編號，令其為𝑖^∗，(𝑖^∗ = arg min_𝑖𝑑_𝑖)。

假如這個最小距離𝑑_𝑖^∗小於我們事先訂定好閥值的𝜏，那我們就把這個當前的棘波 s 歸類到群集𝐶_𝑖^∗，然後我們重新計算該群的平均值𝐭_𝑖^∗並更新它。不然的話，就創造出一個新的群集𝐂_𝑐+，再把當前的棘波 s 分類到這新的群集𝐂_𝑐+，所以該系統的總群數增加了 1 群。

新增了一群後，Osort 演算法會把這新群組的平均值𝐭_𝑖^∗跟其他群組做比較找出最相近的群組，設這最相近群組的平均值為𝐭_𝑗^∗, 𝑗^∗≠ 𝑖^∗。然後會檢查𝑑(𝐭_𝑖^∗, 𝐭_𝑗^∗) 是否小於事先訂定好的第二個閥值𝜏₂，若𝑑(𝐭_𝑖^∗, 𝐭_𝑗^∗) < 𝜏₂，則合併這兩個群組𝐂_𝑖^∗、 𝐂_𝑗^∗。合併後總群組數量減少 1 群。

值(𝐭_𝑖, 𝑖 = 1, … , c)，拿來當作模板回饋給系統使用，這𝐓 的時間可以事先訂定。在這𝐓 的時間內，我們還可以做其他調整，例如:當遇到某個群組內的棘波數量太

在文檔中以回饋式自動模板生成為基礎之正規化關聯值棘波偵測系統之設計及實現 (頁 14-35)

棘波偵測之研究背景與演算法則

第二章 棘波偵測之研究背景與演算法則

2-1 棘波分類的介紹

2-2 棘波偵測演算法則

2-2-1 Matched filter

𝑦[𝑚] = ∑ 𝑥[𝑚 − 𝑘]𝑡[𝑘] = 𝐱

𝐭 (2 − 1)

t

t 這兩向量之間的相關性。 當

2-2-2 Block Energy

𝑦[𝑚] = ∑ 𝑥[𝑚 − 𝑘]

=

(2 − 2)

2-2-3 LRT and GLRT tests

2-2-4 Normalized Correlator

𝐱

𝐱

t

d(𝐱 , t) = ‖ 𝐱 ‖

+ ‖ 𝐭 ‖

- 2𝐱

t. (2-4)

t

d(𝐱 , t) ≤ ‖ 𝐱 ‖

+ ‖ 𝐭 ‖

- 2 η. (2-5)

𝐱 =

, 𝐭̅ =

. (8)

𝑦̅[𝑚] = ∑ 𝑥̅ [𝑘]𝑡[𝑘] = 𝐱

𝐭̅ (2 − 6)

d(𝐱̅ , 𝐭̅) = 2 - 2𝐱̅

𝐭̅. (9)

𝐱̅

𝐭̅ ≤ 1. (10)

η ≤ 1. (11)

d(x̅ , 𝐭̅) ≤ 2(1 − 𝜂) ≡ 𝐷

(12)

2-2-5 Advance Normalized Correlator

‖𝑥 ‖

= ‖𝑥

‖

− 𝑥

[𝑚 − 𝑁] + 𝑥

[𝑚] (13)

𝐱̅

𝐭̅ =

𝑥

𝐭̅ (14)

2.2.6 Osort Algorithm

第二章棘波偵測之研究背景與演算法則

t 這兩向量之間的相關性。當