以回饋式自動模板生成為基礎之正規化關聯值棘波偵測系統之設計及實現

全文

(1)國立臺灣師範大學資訊工程研究所碩士論文. 指導教授：黃文吉博士. 以回饋式自動模板生成為基礎之正規化關聯值棘波偵測系統之設計及實現 Spike Detection Based on Normalized Correlation with Automatic Template Generation. 研究生：王思淮撰中華民國. 一百零三年六. 月.

(2) 誌謝. 誌謝首先，我要先感謝黃文吉教授在研究所入學的時候肯收留我，並在這兩年間的不辭辛勞的指導。不只研究上的問題，連待人處事跟學習上的態度都獲益良多。使得我在這短短兩年碩士的時間中獲得了許多的成長與寶貴的經驗，在此獻上最誠摯的謝意。同時也感謝中央大學通訊工程學系張寶基教授與台灣師範大學光電科技系鄭超仁教授撥冗來參加我的碩士論文口詴，並給予論文上的建議，跟研究的改進方向。. 接著，很榮幸能夠來到國立臺灣師範大學資訊工程研究所就讀，並在多媒體通訊暨系統晶片實驗室這個大家庭和大家一起學習成長。我要先感謝已畢業的學長姐：建廷、清志、國璿、翰逸、任軒、聖穎，謝謝你們在課業上的教導，同時也感謝在我碩士求學生涯與我一起奮鬥的同學們：奇恩、光耀、煥元、皓棠、一修、雅姿、烝祺在研究上給予多方面的協助。. 最後，更要感謝家人支持我繼續出來升學，也感謝之前公司的包容。也要感謝所有關心我的朋友，有了你們的支持、關懷與鼓勵，讓我有勇氣去面對各種的挫折與壓力，並順利地完成學業。謹將此論文成果獻給所有關心我的人，希望大家與我分享這份喜悅與榮耀. i.

(3) 摘要. 中文摘要本論文提出了全新架構的回饋式棘波偵測演算法，主要是用來偵測一個未知棘波特色的棘波序列。此方法在初始階段使用了 Block energy 的棘波偵測法則，接著會把初始階段的結果輸出給 Osort 部份去進行分群並產生模板，最後再利用此模板來進行 Matched filter 的棘波偵測的動作。. 在偵測的過程中，閥值的訂定一直是我們非常困擾的問題，所以我們嘗詴了多種方式來制定出理想的閥值。一開始利用直接定義閥值的方式，給閥值一個訂值，但是此閥值無法適用於各種棘波序列。所以後來利用棘波序列的中間值來自動定義閥值，且在本系統的初始階段中使用它。同時我們也透過了將棘波序列、模板正規化來簡化系統中閥值的訂定，並提供了一個制訂閥值的依據。. 本論文還對棘波偵測系統進行加速的動作，使其不只在命中率上有更優異的表現，在產能上也能有所提升。最後也有將此棘波分類系統在 FPGA 上做實現更進一步的提升其棘波偵測的效能。. 關鍵字：棘波偵測、棘波排序、FPGA、Normalized Correlator ii.

(4) 目錄. 目錄第一章. 緒論 ............................................................................................................. 1. 1-1 研究背景與動機................................................................................................. 1 1-2 研究目的與方法................................................................................................. 2 第二章棘波偵測之研究背景與演算法則 ................................................................. 8 2-1 棘波分類的介紹................................................................................................. 8 2-2 棘波偵測演算法則............................................................................................. 9 2-2-1 Matched filter ............................................................................................. 9 2-2-2 Block Energy ............................................................................................. 11 2-2-3 LRT and GLRT tests ................................................................................ 13 2-2-4 Normalized Correlator............................................................................. 16 2-2-5 Advance Normalized Correlator ............................................................. 20 2-3 DETECTION SYSTEM.......................................................................................... 27 第三章棘波偵測系統架構 ....................................................................................... 29 3-1 BLOCK ENERGY COMPUTATION UNIT .............................................................. 30 3-2 CORRELATOR UNIT ........................................................................................... 31 第四章實驗結果與數據探討 ................................................................................... 32 4-1 4-2 4-3 4-4 4-5 4-6. 加速運算效能實驗........................................................................................... 34 回饋式演算法則的效能評測........................................................................... 36 不同閥值的偵測效果....................................................................................... 40 不同偵測法則的效果比較............................................................................... 45 回饋式演算法的分群效果............................................................................... 47 硬體電路所耗資源........................................................................................... 49. 第五章結論 ............................................................................................................... 52 REFERENCES........................................................................................................... 53. iii.

(5) 目錄. 附圖目錄. 圖 2-1 棘波分類流程圖 ..........................................................................8 圖 2-2 MATCHED FILTER 的運作流程圖 ...............................................9 圖 2-3 BLOCK ENERGY 的運作流程圖................................................. 11 圖 2-4 快速計算的流程圖 ....................................................................23 圖 2-5 OSORT 演算法的流程圖 ...........................................................26 圖 2-6 回饋式棘波偵測系統的架構圖 ................................................27 圖 3-1 棘波偵測系統架構 ....................................................................29 圖 3-2 BLOCK ENERGY COMPUTATION UNIT .......................................30 圖 3-3 CORRELATOR UNIT ....................................................................31 圖 4-1 SNR = 8 時的棘波序列 .............................................................33 圖 4-2 SNR = -2 時的棘波序列 ...........................................................33 圖 4-3 SNR = 8 時本系統的棘波圖形: ................................................38 圖 4-4 SNR = -2 時本系統的棘波圖形: ...............................................38 圖 4-5 模板跟雜訊的摺積值 ................................................................43 圖 4-6 模板跟棘波的摺積值 .................................................................43. iv.

(6) 目錄. 附表目錄. 表格 1-1 獲得棘波序列的方法 ................................................................2 表格 1-2 棘波偵測的兩大種類 ................................................................4 表格 2-1 比較不同類型 MATCHED FILTER 的計算複雜度。 ...............22 表格 4-1 不同 SNR 下各種 MATCHED FILTER 的運算速度 ................34 表格 4-2 當𝐓𝟏使用不同的秒數時，對第二階段偵測的影響 .............37 表格 4-3 我們提出的演算法在不同 SNR 下的效能............................39 表格 4-4 第二階段在不同閥值下的 TP RATE 跟 FA RATE 表現 ........40 表格 4-5 各種棘波偵測法則在不同雜訊等級的棘波序列下的 ..........45 表格 4-6 不同 SNR 值下，棘波序列偵測和分群的效果....................47 表格 4-7 三種棘波形狀的棘波序列偵測和分群的效果......................48 表格 4-8 電路簡化前後的硬體資源消耗 ............................................49 表格 4-9 整體電路的消耗 ......................................................................50 表格 4-10 硬體跟軟體執行時間比較 ....................................................51. v.

(7) 第一章. 第一章. 緒論. 緒論. 1-1 研究背景與動機大腦負責處理人類的一切感官資訊，是人類身體中最精密且複雜的器官。人類可以動作、記憶、學習、還有透過五感去感受這世界的一切，都是利用佈滿腦內的神經元(Neuron)所互相聯絡完成的。人類的神經元大概有十億個，因為有如此龐大數量的神經元，我們人類可以完成許多複雜的動作，還有產生多元的情感。即使是科技發達的現在，人類對於大腦的了解還是相當有限，有研究指出我們終其一生也只發揮了大腦 5%左右的能力而已。所以我們打算透過研究神經元訊號，進一步的了解大腦的運作，也許可以用來開發大腦的能力或是去重建且創造人們的感官。神經元，是一種能被電流訊號激發的神經系統核心單元，它能透過電流和化學能信號交替轉換的過程來傳遞信息。神經元所發出的訊號，其實就是細胞膜內鈉鉀離子變化所產生的電位差，我們稱之為動作電位(Action Potential)。而偵測棘波的方法一般有兩種方式(1) 透過植入體內的微電殛探針(2)體外的有線裝置來偵測訊號。將偵測到的這一連串動作電位組合起來，就形成了所謂的棘波序列(Spike Train)。棘波排序這項技術的目的就是把偵測到的混合訊號中將來自不同神經元細胞的訊號給區分出來。主要可分為棘波偵測、特徵擷取、棘波分類這三個部分。 1.

(8) 第一章. 緒論. 1-2 研究目的與方法前面有提到棘波排序(Spike Sorting)[1]的技術主要分成棘波偵測、特徵擷取、棘波分類這三個部分。因為棘波偵測是棘波排序演算法中的第一個階段，所以棘波偵測效能的好壞就直接影響整個系統的優異了。所以本論文主要是提出一個新的回饋式棘波偵測演算法來提升棘波偵測的準確度，進一步的提升整個棘波排序的效能。而要決定方法之前，我們先從了解預計要偵測的資料來源開始，接著再介紹各部分所用到的技術。棘波排序的資料來源通常是透過腦機介面(Brain Machine Interface，BMI)[2]所獲得。腦機介面獲得棘波的方式主要可以分成兩種，(1) 侵入式 (2) 非侵入式。其優缺點如下表所示。. 侵入式. 非侵入式. 優點. 偵測正確率高. 安全性高. 缺點. 危險性較高. 正確率低. 表格 1-1 獲得棘波序列的方法. 2.

(9) 第一章. 緒論. (1) 侵入式所謂的侵入式腦機介面，就是指在顱內植入微電極探針，並利用探針偵測周圍的神經元細胞所發出的訊號，其好處是用這種方式獲得的訊號可靠度較高，比較不會受到雜訊的干擾。相對的壞處就是要將探針置入腦內的危險性較高，通常需要有專業的外科醫生在旁協助。 (2) 非侵入式所謂的非侵入式腦機介面，是指利用外部的有線裝置來獲得腦部的神經元活動資訊，這種方式的好處是對人體的危害程度低，也不需要專業的外科醫生在旁協助，只是相對來說，偵測出來的訊號受雜訊的影響就比較嚴重。. 不管是用哪種方式偵測，所獲取的棘波序列都是由不同的神經元細胞所發出的訊號的集合，且難免都會受到雜訊的干擾，所以一個好的棘波排序系統，其第一階段的棘波偵測的準確率相當重要，若第一階段的時候就能將受雜訊跟重疊干擾的棘波序列處理好，那後續的棘波分類階段就能更準確的分辨出來該棘波是不是屬於同個神經元細胞所發出的。. 而現行的棘波偵測方法中[3]，大致上可以分成兩種偵測方式(1) 振幅偵測 (Amplitude) 跟 (2)能量計算(Energy)。其比較如下表。. 3.

(10) 第一章. 方法. Amplitude. Energy. 運算複雜度. 低. 略高. 雜訊:低. 弱. 好. 雜訊:高. 弱. 弱. 偵測效果. 緒論. 表格 1-2 棘波偵測的兩大種類. (1) 振幅偵測(Amplitude) 最基本的振幅式棘波偵測方法，是訂定一個閥值，然後對棘波序列取絕對值，若棘波序列上有超過該閥值的部分，就判斷該部分為棘波。這個偵測法則的缺點就是準確率偏低。 (2)能量計算(Energy) 能量計算式的棘波偵測，主要是計算各點之間的關聯性，例如 NEO[4]，其演算法的核心為計算棘波序列中各點的平方與相鄰此點左右兩點的乘積之能量差，然後再判斷該能量是否大於閥值。這樣的優點在於偵測結果不會被突然出現的雜訊給影響。除了 NEO 之外還有 SWT(Stationary Wavelet Transform)[8,9]也都是屬於計算能量式的棘波偵測。以上提到的這些方法的運算複雜度都很低，效能也不錯，但是閥值難訂定，就算搭配了其他自動定義閥值的演算法[4,6,10]。而且當雜訊值逐漸增大的時候，偵測效能也會急速下降。. 4.

(11) 第一章. 緒論. 比較上面兩種偵測方法，我們新提出來的演算法基本上就使用了(2)能量計算式的棘波偵測法。但是為了改善(2)計算能量式的棘波偵測的缺點，我們也使用了另一種利用模板來偵測棘波的演算法。這種利用模板來偵測的演算法中，最典型的例子就屬於 Matched filter[5,11,12]了。. Matched filter 演算法是利用特定的棘波產生一組模板，並利用此模板與棘波序列做摺積後的值來進行偵測，若值大於閥值，就是偵測到棘波。這演算法可以視為是一種 correlator。 Matched filter 已經被廣泛運用在通訊系統中的雜訊偵測，而它的缺點在於效能太過依賴模板的選擇了，當模板與棘波序列中的波型相差太多的時候，偵測錯誤的機率就會大幅的提升。而且由於無法預測出棘波序列中所包含的棘波種類，所以也無法得知該準備幾種模板來偵測棘波序列。. 為了改善以上的缺點，有人提出了自動模板生成的 Matched filter[13]，該演算法的模板是利用高於閥值的棘波的中間值計算所獲得，而獲得的模板可以對應於不同的神經元發出的訊號。該演算法雖然解決了模板產生的問題，但是只會產生一個模板，而靠這單一個模板要偵測所有的棘波訊號似乎略顯不足。. 不管是振幅式或是能量式的棘波偵測，以上那些現有的棘波偵測方式的共同特點就是只在棘波偵測階段處理有雜訊的棘波序列，偵測完的結果就丟到特徵擷取階段，從此以後這些資訊就不會再被棘波偵測階段給使用。所以我們提出了一 5.

(12) 第一章. 緒論. 種新的回饋式棘波偵測系統。其系統會利用特徵擷取階段的結果產生出模板，再回饋給棘波偵測階段來提升棘波偵測的準確率。因為在特徵擷取階段的結果通常含有許多棘波序列的資訊。例如: 該棘波序列中所含的棘波個數和所含棘波的峰值，這些資訊有助於幫助模板的產生以及知道必頇產生的模板個數，並估計出各種波型的平均值。. 在棘波擷取的演算法[14-16]上，有許多是使用 Principal Component Analysis(PCA)[17]或是其相似法則。而在分群方面，則是 K-means、Fuzzy C-mean(FCM)[18]這兩種演算法較為廣泛使用，但是他們都需要經過長時間的離線訓練才可以達到一定的分群效果，所以很難用於即時的回饋式棘波偵測系統上。根據上述原因，我們新提出的回饋式棘波偵測系統在特徵擷取方面是使用 Osort 演算法。我們利用 Osort 演算法經過非監督式的棘波分類產生出模板，所以它並不需要經過特徵擷取跟分群的訓練，而在分群方面的效能也比 WaveClus[10] 和 KlustaKwik[14]這兩種演算法優異，所以我們的回饋式棘波偵測系統選擇使用了 Osort[19]法則。. 我們結合了 Osort 演算法解決了模板的產生問題後，接下來要解決棘波偵測系統中，閥值不好訂定的問題。我們詴著將模板和棘波序列做 Normalized correlator，並計算兩者的平方誤差(squared distance)，再經過一些簡單的推導來制定出適當的閥值。. 6.

(13) 第一章. 緒論. 到這邊為止，我們解決了模板的產生跟閥值的訂定這兩大問題。但是因為我們新提出來的架構，對模板跟棘波序列多做了 Normalized correlator 運算，所以會比傳統的 Matched filter 更加耗時。於是我們提出了一個加速計算的演算法，改善了這個缺點，此加速演算法是透過將能量過低的棘波序列過濾掉，使其不進行 Normalized correlator 的計算，進而降低整個系統的計算時間。. 統整一下上面所講的，我們提出的回饋式棘波偵測系統的主軸是使用由 Osort 演算法提供的模板來進行 Matched filter 棘波偵測。而在一開始的初始階段，因為還沒有棘波分類的結果回饋給 Matched filter 棘波偵測運用，所以我們是使用 generalized likelihhod ratio test(GLRT)[7]的動作來進行偵測。. 為了評估系統，我們使用由 Professor Leslie S. Smith 所開發的神經元訊號模擬器[21]，來產生各種雜訊干擾下的棘波序列。根據實驗結果顯示，我們提出的回饋式棘波偵測系統，可以有效的利用由 Osort 產生的模板來提升偵測的準確度，而且該演算法在高雜訊的影響下，跟其他基於能量計算的棘波偵測演算法比較下，還是有很好的表現。. 7.

(14) 第二章. 棘波偵測之研究背景與演算法則. 第二章棘波偵測之研究背景與演算法則 2-1 棘波分類的介紹傳統的棘波分類演算法主要包含了棘波偵測 (Spike Detection)、特徵擷取 (Feature Extraction) 與棘波分類 (Clustering) 三大步驟，詳見於圖 2-1。棘波分類演算法會先將輸入的棘波序列 (Spike Train) 經由棘波偵測階段來偵測出棘波序列中各種棘波發生的位置，然後將其截取出來並做對齊的動作。這些經過對齊過後的棘波會再經由特徵擷取的階段，擷取出所有棘波的特徵值，最後再根據這些特徵值進行棘波分類。棘波偵測階段最常碰到的問題是，一旦輸入的棘波序列受到大量雜訊的干擾，其原本的波形就會嚴重被破壞，使得偵測效能大幅下降。棘波偵測經常被用於醫學上或是大腦研究，一般來說我們會以植入式的微電極探針或是外接裝置來進行腦波的擷取，在擷取過程中都會受到雜訊的干擾，尤其是外接裝置受到的影響更是顯著，所以要將棘波從大量雜訊干擾的環境中擷取出來，是所有研究棘波偵測者共同面臨的最大問題。. 1. 2. Spike Detection. 3. Feature Extraction. 圖 2-1 棘波分類流程圖 8. Clustering.

(15) 第二章. 棘波偵測之研究背景與演算法則. 2-2 棘波偵測演算法則. 2-2-1 Matched filter 我們先從最基本型的Matched filter technique 討論起。Matched filter 是用來進行棘波偵測的一種技術，其核心的技術是利用多組事先儲存好的模板(template) 拿來循序的與棘波序列(Spike train)的區段(𝑥[m])做摺積運算。大致上的流程如下圖所示:. m=0 𝑥 =0. m=m+1. dot(𝑥 ,template). N miss. >Threshold Y. HIT. 圖 2-2 Matched filter 的運作流程圖. 9.

(16) 第二章. 棘波偵測之研究背景與演算法則. 為了方便解釋，我們先假設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 環境下做出有效的偵測。. 10.

(17) 第二章. 棘波偵測之研究背景與演算法則. 2-2-2 Block Energy 接下來介紹 Block Energy 演算法，這個演算法在進行棘波偵測的時候，也會將輸入進來的棘波序列先分成數個區段，然後再觀察每個區段中的採樣點彼此平方相加的值，若值超過了事先訂定好的閥值，則該值相對應的區段就會被判斷為棘波。流程如下圖:. 圖 2-3. Block energy 的運作流程圖. 11.

(18) 第二章. 棘波偵測之研究背景與演算法則. 假設𝑥[m]為棘波序列中第 m 個樣本(sample)。𝑥 = [x[𝑚], x[𝑚 − 1], … , x[𝑚 − N + 1]]𝑇 為棘波序列中第 m 個區段(segment)，其中 N 為區段長度。y[m]定義為 Block energy 在棘波序列中第 m 點的能量輸出值。. 𝑁. 𝑦[𝑚] = ∑ 𝑥[𝑚 − 𝑘]2 = 𝑥𝑇𝑚 x𝑚. (2 − 2). 𝑘=. 從Block energy的流程圖跟公式可以看的出來，雖然在偵測的觀念上略有不同，但是跟前面提到的Matched filter 的運作卻頗為相似，Block energy只是把 Matched filter 中模板的部分換成了自身的區段x[m]。 Block Energy的優點在於計算複雜度低，在雜訊等級低的時候棘波偵測效果好。缺點是當輸入訊號的雜訊太高的時候，因為區段內的能量會被雜訊影響而提高，所以就不好分辨出該區段是屬於棘波或是雜訊，導致偵測效果惡化。. 12.

(19) 第二章. 棘波偵測之研究背景與演算法則. 2-2-3 LRT and GLRT tests 假設模板 t 是已知的棘波，而且輸入系統的棘波序列中的雜訊 n 是屬於 additive zero-mean white Gaussian noise 的話，則對區段𝐱 和模板 t 做 Match filter (𝐱 𝑇 t ≥ η) 的意義等於對區段𝐱 做 Likelihood ratio test (LRT)。所以我們現在觀察區段𝐱 ，LRT test 的目的是要根據區段𝐱 判斷以下兩個假說(分別設為 H0、H1) 何者為真。 H0 : xm = n, H1 : xm = t + n.. 因為一開始就假設雜訊 n 是屬於高斯分布，所以 n 的發生機率為 N(0,Σ)， Σ = σ2 𝐈，σ2 是雜訊的變異數，I 是單位矩陣。 𝐱 在 H0 假說成立下發生的機率為 P(xm/H0)，也等於 N(0,Σ)，同理，𝐱 在 H1 假說成立下發生的機率為 P(xm/H1)= N(t,Σ)。. P(x |H ) =. 𝑥𝑇 𝛴. 𝑡. √2𝜋|Σ|. e. H > < H0. 13. 1 (x 2 𝑚. 𝐭)𝑇 Σ−1 (x𝑚 𝐭). (2). 𝑃(𝐻 )(𝑐. 𝑐00 ). 1. 𝑐11 ). log 𝑃(𝐻0 )(𝑐10 01. ≡𝛾. (3).

(20) 第二章. 棘波偵測之研究背景與演算法則. 設 cij 是當假說H𝑗 為真，但我們卻接受假說H𝑖 的機率， P(H𝑖 ) 是x 屬於在 𝑃(𝐻 )(𝑐. 𝑐00 ). 1. 𝑐11 ). 假說H𝑖 的機率。根據 optimal test[22]，我們設 log 𝑃(𝐻0 )(𝑐10 並給予符號 γ，並檢驗之，若 𝑥 𝑇 𝛴. 01. 為檢定統計量，. 𝑡 > γ 則x 符合H 假說，反之則符合H0 假. 說。觀察(式 3)，因為Σ = σ2 𝐈，故我們可以把其代換成:. x𝑇 𝑡. H > < H0. 𝜎 2𝛾 ≡ η. (2-3). 觀察(式 2-3)與(式 2-1)，可以發現若把閥值η定為𝜎 2 𝛾，則對x 做 LRT test 等價於對x 做 matched filter。. 若我們一開始能知道棘波的資訊，則對於 matched filter 演算法的準確性就能大幅的提升，但是這在實際的過程中可能無法實現，我們從欲偵測的區段x 中先估計出一個模板 t 來解決這個問題。假設x 是符合假說H ，我們使用 maximum likelihood (ML) 來估計出模板 t，當模板𝐭 = x 時，模板 t 使得P(𝑥 |𝐻 )的機率最大，故我們選擇此情況，並新定義一個符號𝐭̂。. 𝐭̂ = x. = maxP(x |H ). 14. (4).

(21) 第二章. 棘波偵測之研究背景與演算法則. 在這個情況下，假說H 不用依賴模板 t 的值，故問題從 LRT test 變成 GLRT ， (式 4)可改寫為 :. H > < H0. 𝑥𝑇 Σ x. log. 𝑃(𝐻0 )(𝑐10 𝑐00 ) 𝑃(𝐻1 )(𝑐01 𝑐11 ). 為了方便起見，我們假設代價為 uniform cost (c00 = c. ≡𝛾. (5). = 0，c0 = c. 0. = 1)，. 另外因為不失一般性，我們又假設沒有 Spike 的機率比有 Spike 產生的機率高 (P(𝐻0 )>P(𝐻 ))，這時我們觀察(式 2-3)，得知 γ>0。若Σ = σ2 𝐈，GRLT tset 等效於 block energy 𝑥 𝑇 𝑥 ，且閥值η = 𝜎 2 𝛾。. H > < H0. 𝑥𝑇 x. σ2 𝛾 ≡ 𝜂. (6). 總而言之，GLRT 是一種 block energy detector。假設x 的平均值是 0，我們可以用以下方法估計出區段x 的標準差 σ[10]:. σ = median {. |𝑥[𝑘]| 0.6745. , 𝑘 = 𝑚 − 𝑀 + 1, … , 𝑚}. (7). 我們可以根據上式，透過對輸入訊號採樣得到的中間值估計出 σ，符號 M 代表: 對輸入訊號的採樣個數。. 15.

(22) 第二章. 棘波偵測之研究背景與演算法則. 2-2-4 Normalized Correlator 回顧上述提到的棘波偵測演算法，Matched filter 和 LRT test 其實都有一個很大的缺點，就是制定理想的閥值是一個很大的挑戰。觀察(式 3)可知，如果我們沒有辦法知道𝑃(𝐻0 )、𝑃(𝐻 )和𝜎 2 的資訊，就很難估算出理想的 η。. 因此，我們選擇另外一種方式制定閥值，透過量測模板 t 跟區段𝐱 的誤差上限來決定 η。我們事先訂定其誤差上限的值，當輸入的區段𝐱 小於誤差上限，就判斷為棘波。這樣一種 correlator 的設計方式，可以對應到另一個僅用平方誤差的簡單比對技術。. 我們利用平方誤差 (squared distance) 來觀察區段 𝐱. 與模板 t 之間的關. 係，我們首先可以觀察下面所列出的平方誤差 (squared distance)公式：. d(𝐱 , t) = ‖ 𝐱 ‖2 + ‖ 𝐭 ‖2 - 2𝐱 T t.. (2-4). d(𝐱 , t) ≤ ‖ 𝐱 ‖2 + ‖ 𝐭 ‖2 - 2 η.. (2-5). 因此, 當 𝐱 𝑇 t > η 時,. 觀察上式得知，當偵測到 𝐱. 為棘波時 (i.e., 𝐱 𝑇 t > η)， d(𝐱 , t)的上限值取決於. ‖ 𝐱 ‖2、‖ 𝐭 ‖2 及 η，而‖𝐱 ‖2 大小取決於輸入棘波序列。所以當 ‖ 𝐱 ‖2 變大時， d(𝐱 , t)的值就會跟著變大，即使 𝐱 𝑇 t > η。在這種情況下，就可能會把雜訊判斷成棘波，而產生 false alarm 的情形。 16.

(23) 第二章. 一種解決上述問題的方法，就是讓 𝐱. 棘波偵測之研究背景與演算法則. 及 t 在做 correlation 計算前先做正. 規化(Normalization)的動作。接著我們將做完正規化的 𝐱. 及 t 分別定義為 𝐱̅. 與 𝐭̅ ，如下 𝐱𝑚. 𝐱 =. ‖𝐱 𝑚. , 𝐭̅ = ‖. 𝐭 ‖𝐭‖. .. (8). 而 Normalized correlator 之輸出為 𝑁 𝑇. 𝑦̅[𝑚] = ∑ 𝑥̅ [𝑘]𝑡[𝑘] = 𝐱 𝐭̅. (2 − 6). 𝑘=. 其中 𝑥̅ [k] 及 𝑡[𝑘] 分別為 𝐱. 及 𝐭̅ 之第k個element，. 因此，(式 2-5)可以改寫成. d(𝐱̅ , 𝐭)̅ = 2 - 2𝐱̅ 𝑇 𝐭.̅. (9). 因為 d(𝐱 𝑛 , t) ≥ 0, 我們可以得到. 𝐱̅ 𝑇 𝐭̅ ≤ 1.. 所以我們的 Normalized correlator 是指 𝐱̅ 們就可以說此時對應到的 𝐱. (10). 與 𝐭̅ 上的計算。當 𝐱̅ 𝑇 𝐭̅ ≥ η 時，我. 就是我們偵測到的棘波。從(式 10) 我們可以推出. η ≤ 1. 17. (11).

(24) 第二章. 棘波偵測之研究背景與演算法則. 另外，當 𝐱 𝑇 t > η ，根據(式 9) 我們可以得到. d(x̅ , 𝐭)̅ ≤ 2(1 − 𝜂) ≡ 𝐷. (12). D可以表示為HIT時的平方誤差上界，經過Normalized correlator後得到的唯一閥值，可以當成使用在template matching時的平方誤差上界。此外d(x̅ , 𝐭)̅ 的值取決於閥值η的大小，當閥值η變大時d(x̅ , 𝐭)̅ 會變小，而閥值η的最大值是1，且並不會受到輸入的棘波序列所影響。. 事實上，用來判斷輸入的棘波區段xm 跟模板t是否匹配的平方誤差最大值為 D = 2(1 − 𝜂)，當閥值η=1時。Normalized correlator這方法對於閥值η提供了一種新的解釋。如果我們設定閥值η=1，也就是當區段xm 跟模板t 做運算後其值超過 η(x̅ 𝑇 𝐭̅ ≥ 1)才算是偵測到棘波的話，就代表該區段xm 跟模板t有完全相關性(full correlation)，這時他們之間的平方誤差為0 。換句話說，如果當閥值η=0.5，就表示了區段xm 跟模板t 有半相關性(half correlation)，只要區段xm 跟模板t運算後的值大於η(x̅ 𝑇 𝐭̅ ≥ 0.5)就算偵測到棘波，此時它們之間的平方誤差上限值為1。然而如果當閥值η=0，就表示了區段xm 跟模板t 毫無相關性(no correlation)，代表只要有區段xm 進來就算偵測到棘波，此時它們之間的平方誤差上限值為2。. 18.

(25) 第二章. 棘波偵測之研究背景與演算法則. 不過在真實情況下因為有雜訊的干擾，我們無法使用完全相關性(full correlation)當偵測標準(η=1)，根據我們實驗的結果，即使在雜訊很高的環境下我們只需要 70% 的相關性 (i.e. η = 0.7)，就可以讓我們 Normalized correlator 偵測達到高True Positive,低 False Positive的效果。. 19.

(26) 第二章. 棘波偵測之研究背景與演算法則. 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 主要的運作原理是利用𝐱. 跟𝐱 這兩者的相關性。. 觀察下式:. ‖𝑥 ‖2 = ‖𝑥. 因為‖𝑥. ‖2在運算𝐱̅ 𝑇. ‖2 − 𝑥 2 [𝑚 − 𝑁] + 𝑥 2 [𝑚]. (13). ̅ 就已算過了，所以當我們想得到‖𝑥 ‖2 時，僅需要再 𝐭時. 使用 2 個乘法、2 個加法和 1 個平方根運算。. ̅ 順序，減少運 Post-correlation normalization 技術主要透過改變了計算𝐱̅ 𝑇 𝐭的算次數，進而增快運算時間。方法如下式:. 𝐱̅ 𝑇 𝐭̅ = ‖𝑥. 𝑚‖. 20. 𝑥 𝑇 𝐭̅. (14).

(27) 第二章. 棘波偵測之研究背景與演算法則. 改善前的運算流程，都是先算出𝐱̅ 𝑇 後再跟 t 做內積，現在我們先計算𝑥 𝑇 𝐭̅ ̅ 純量)。再對它做正規化。如此一來，在這個階段我們就只需要 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 的技術核心是在比較‖𝑥 ‖2 跟‖𝐭‖2 的大小，假如‖𝑥 ‖2 <λ‖𝐭‖2 ， 0 < λ < 1，‖𝑥 ‖2 就會先因為能量過低而被判斷成不是棘波。這個方法可以有效地降低正規化的次數。. 21.

(28) 第二章. Matched Filter Basic Normalized Correlator Fast Normalized Correlator. 棘波偵測之研究背景與演算法則. Addition. Multiplication. Division. Squared Root. N-1. N-. 0. 0. 2N-2. 2N. N. 1. N+1. N+2. 1. 1. p×(N−1) + 2. p×N+2. p. 1. (No Pre-Screening) Fast Normalized Correlator (Pre-Screening). 表格 2-1 比較不同類型 Matched filter 的計算複雜度。. 從實驗結果可以發現，結合了 Pre-Screening 的 Fast normalized correlator 可以再更一步的增加其運算速度，甚至比最基本型態的 Matched filter 更有效率。令 p (0<p<1) 是‖𝑥 ‖2 <η‖𝐭‖2 成立的機率，如(表 2-1)所示，Pre-Screening 的 Fast normalized correlator 的平均的加法數量為 p(N-1)+2、乘法除法數量為 p*N+2、除法的數量是 p 次。. (圖 2-4)是棘波偵測部分的流程圖，流程圖裡詳細的表達了在一個棘波序列進來後，我們先利用 Fast energy computation 簡化了獲得‖𝑥 ‖2 的運算，再利用 Pre-Screening 過濾掉能量太低的區段，減少 correlation 的次數，最後再利用 Post-correlation normalization 降低內積的運算成本。. 22.

(29) 第二章. 圖 2-4 快速計算的流程圖. 23. 棘波偵測之研究背景與演算法則.

(30) 第二章. 棘波偵測之研究背景與演算法則. 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 演算法會把這新群組的平均值𝐭 𝑖 ∗ 跟其他群組做比較找出最相近的群組，設這最相近群組的平均值為𝐭𝑗 ∗ , 𝑗 ∗ ≠ 𝑖 ∗。然後會檢查𝑑(𝐭 𝑖 ∗ , 𝐭𝑗 ∗ ) 是否小於事先訂定好的第二個閥值𝜏2，若𝑑(𝐭 𝑖 ∗ , 𝐭𝑗 ∗ ) < 𝜏2，則合併這兩個群組𝐂𝑖 ∗ 、 𝐂𝑗 ∗ 。合併後總群組數量減少 1 群。. 24.

(31) 第二章. 棘波偵測之研究背景與演算法則. 在經過𝐓 秒後，我們就可以把這些經由 Osort 演算法產生的棘波群組的平均值(𝐭 𝑖 , 𝑖 = 1, … , c)，拿來當作模板回饋給系統使用，這𝐓 的時間可以事先訂定。在這𝐓 的時間內，我們還可以做其他調整，例如:當遇到某個群組內的棘波數量太少時，我們就刪除掉該群組，因為該群組的產生可能是因為雜訊的誤判造成的; 或是經由計算𝐭 𝑖 中最大值跟最小值之間的距離，若太接近或是太遠，我們便判斷該群組為雜訊，並刪除該群。. 25.

(32) 第二章. 圖 2-5. 棘波偵測之研究背景與演算法則. Osort 演算法的流程圖. 26.

(33) 第二章. 2-3. 棘波偵測之研究背景與演算法則. Detection System. 圖 2-6 回饋式棘波偵測系統的架構圖. 上圖(圖2-6)是我們新提出來的回饋式棘波系統的架構圖，從圖中可以看出該架構是由GLRT test為基礎的棘波偵測系統、正規關聯性棘波偵測系統跟Osort algorithm所組成的。GLRT test為基礎的棘波偵測系統跟正規關聯性棘波偵測系統負責處理棘波偵測的部分，而Osort algorithm負責分群的部分。而此新架構的特點在於用來跟棘波序列對照的模板是自動產生的，當Osort經過一定的時間後，會產生數個模板並回饋給系統使用。. 27.

(34) 第二章. 棘波偵測之研究背景與演算法則. 該系統在一開始的時候沒有模板可以用來做棘波偵測，所以在初始階段，我們利用以GLRT test為基礎的block energy 棘波偵測系統來進行棘波偵測，Block energy棘波偵測系統所需要用到的閥值是由系統自動產生的(透過式6和式7)。偵測到的棘波再經過Osort algorithm做分群跟模板的產生，當Osort經過一段時間(𝐓 ) 的運算後，就會把產生好的模板回饋給系統使用。. 系統在使用模板前，會先將模板做Normalized correlation，這是因為透過這個動作可以簡化閥值的選擇(式11)，我們也可以根據自己的需求去訂定閥值，每個閥值都對應到不同的關聯性(式12)，在這個系統裡面，我們為了確保棘波偵測的準確度，在模板比對的時候選擇了固定的閥值。並且在每經過T2 的時間後，會再透過Osort演算法自動更新所需要用的模板。適時的更新模板對於長時間的棘波偵測效能會更加穩定。. 28.

(35) 第四章. 實驗結果與數據探討. 第三章棘波偵測系統架構本章節會介紹一下我們所提出來的棘波偵測系統的電路架構，本架構中包含四大單元，如下圖(3-1)所示，分別為 Filter Unit、Block Energy Computation Unit、 Correlator Unit 以及 Thresholding Unit。Filter Unit 會將輸入的棘波序列 s[m] 的直流成分與高頻雜訊的部分去除，並得到輸出 x[m] 。接著 x[m]會被 Block Energy Computation Unit 拿去做平方相加的計算得到‖𝐱 ‖2 。之後 Correlator Unit 會將 x[m] 及‖𝐱 ‖2 跟事先準備好的模板 𝐭̅ 與 𝐭̅2 做 correlation 計算，得到 ̅ [m] 與 ̅2 [m]。最後 Thresholding Unit 將 ̅ [m] 與 ̅ 2 [m] 的值與閥值 η 做比較，若大於閥值就判斷為棘波。接下來會描述一下 Block Energy Computation Unit 跟 Computation Unit 的架構。. s[m]. x[m]. Filter Unit. Block Energy Computation Unit ‖𝐱 𝑚 ‖2. 𝑦̅ [m]. Correlator Unit. 𝐭̅. 𝑦̅2 [m]. Thresholding Unit. η. 𝐭̅ 2. 圖 3-1 棘波偵測系統架構. 29. Detection Output.

(36) 第四章. 實驗結果與數據探討. 3-1 Block Energy Computation Unit 下圖(圖 3-2)可以看出 Block Energy Computation Unit 內部架構的運算過程，在此架構中，我們有使用前面所提到的加速演算法 Fast energy computation 來降低我們的 area cost。當中 x[m]是資料寬度為 8-bits 的棘波採樣點， Reg1 到 Reg64 是寬度 8-bits 的暫存器，而 Reg0 是 16 個 bits 的暫存器，儲存 Reg64 中值的平方。此單元會先將 Filter Unit 輸出的訊號先存入一個 64 階的位移站存器，再將 x[m] 平方後的值累加，同時減掉位移暫存器最後一階(i.e. Reg64)的值的平方後再存入 Reg0，透過不斷加上最新的 x[m] 與減掉延遲 64 個時脈週期的 x[m - 64] 之計算，就可以達到以 sliding window 的方式計算出 ‖𝐱 ‖2 的值。. x[m]. -. +. x. Reg1. 圖 3-2. ‖𝒙𝑚 ‖2. .... Reg0. Reg64. Block Energy Computation Unit. 30. x.

(37) 第四章. 實驗結果與數據探討. 3-2 Correlator Unit 圖 3-3 是 Correlator Unit 內部的架構，這電路架構已經有將之前提到的加速流程 Post-correlation normalization 作實現了，若要比較此單元簡化前後之硬體資源消耗的話，後續會在實驗結果上再作探討。我們可以看到 Correlator Unit 會把 ‖𝐱 ‖2 經由 SQRT 電路與 INVERSE 的計算後變成. ‖𝐱𝑚 ‖. 。同時，也會把輸入. 的棘波序列區段 x[m]，跟事先準備好的兩個模板𝐭̅ = [𝑡̅ [1], . . . , 𝑡̅ [64]]𝑇 、 𝐭̅2 = [𝑡2̅ [1], . . . , 𝑡2̅ [64]]𝑇 做內積。最後會把 x[m]與 𝐭̅ 及 𝐭̅2 內積的值再乘以得到𝑦̅ [m]、𝑦̅2 [m]。. ‖𝐱 𝑚 ‖2. SQRT. x[m]. 𝒕̅ [1]. INVERS. …. Reg1. 1 ‖𝐱 𝑛 ‖. 𝒕̅ [64]. Reg64. …. 𝒕̅2 [64]. … +. +. x. x. 𝑦̅ [m]. 圖 3-3. 𝑦̅2 [m]. Correlator Unit. 31. x. x. x. x. 𝒕̅2 [1]. ‖𝐱𝑚 ‖. ，.

(38) 第四章. 實驗結果與數據探討. 第四章實驗結果與數據探討為了評估我們新提出來的演算法，我們實驗用的輸入資料是使用由[21]所產生的棘波序列，這模擬器不只可以產生出類似真實情況的棘波，也提供了完整的棘波序列資訊，例如:棘波個數、棘波種類、棘波出現的時間。所以我們可以很有效的利用這個模擬器來評估系統的效能。在實驗的過程中，我們模擬了各種不同的 SNR(signal-to-noise)值情況下的棘波序列，而採樣頻率(Sampling rate)是維持在一秒 24,000 個採樣點，因為每個棘波訊號的長度是 2.67 毫秒，所以我們每個棘波訊號的長度為 64 個採樣點(N=64)。. 所謂的 SNR(signal-to-noise)，又稱為訊噪比（Signal/ Noise），它是指訊號強度與同時發出雜訊強度之間的比率，單位為分貝(dB)。SNR 值越高，代表這訊號的品質越好，換句話說 SNR 值越低，代表這個訊號受雜訊的干擾非常嚴重。在我們的實驗中，SNR=100 為沒有雜訊干擾的棘波序列。透過以下我們擷取的 SNR=-2 (圖 4-1)還有 SNR=8 (圖 4-1)時候的兩個棘波序列，可以更清楚的了解到訊噪比的高低對於訊號的影響。下圖中標記的部分代表該區段有棘波產生。. 觀察(圖 4-1)跟(圖 4-2)，可以發現隨著 SNR 的下降，我們越來越不容易判別棘波序列中的雜訊跟棘波，這也代表著我們棘波偵測的難度也隨之提高了。. 32.

(39) 第四章. 實驗結果與數據探討. 圖 4-1 SNR = 8 時的棘波序列. 圖 4-2. SNR = -2 時的棘波序列. 實驗步驟的部分，我們設計了五個實驗來評估我們新提出的演算法: (1)確認本系統的加速計算的效能，(2)檢視本系統的棘波偵測效能，(3)證明我們能有效的簡化閥值的選擇，(4)本系統跟其他演算法做比較，(5)本系統產生模板的效能。. 33.

(40) 第四章. 實驗結果與數據探討. 4-1 加速運算效能實驗首先，我們為了確認前面提到的加速系統的效能，而用 C code 測量了不同 template matching 的棘波偵測方法的計算時間。此實驗輸入的棘波序列長度為 100 秒，共 2,400,000 個採樣點，棘波序列內含有兩種不同類型的棘波，其結果如下:. SNR(dB). -2. 0. 2. 4. 6. 8. Matched Filter. 0.67 sec. 0.67 sec. 0.67 sec. 0.67 sec. 0.67 sec. 0.67 sec. 1.58 sec. 1.58 sec. 1.58 sec. 1.58 sec. 1.58 sec. 1.58 sec. 0.44 sec. 0.40 sec. 0.37 sec. 0.28 sec. 0.27 sec. 0.26 sec. Basic Normalized Correlator Fast Normalized Correlator (Pre-Screening). 表格 4-1 不同 SNR 下各種 Matched filter 的運算速度. 我們觀察表 4-1. 可以發現，Basic Normalized Correalator 所花時間明顯比其他兩種方式高，因為此方法是我們對 matched filter 又多做了正規化的動作，所以其時間花費理所當然比 matched filter 高。. 34.

(41) 第四章. 實驗結果與數據探討. 而 Fast Normalized Coeerlator 是從 Basic Normalized Correalator 添加了幾種減少運算的技術而得來，其中用到的技術有: 為了減少計算能量複雜度的 fast energy computation(式 13)跟 post-correlation normalization(式 14)，還有為了降低 Normalized Correalator 運算次數而使用的 pre-screening，這裡的 pre-screening 預設為λ=0.5。. 以上方法所減少的計算量可以由(表 2- 1)看出。然而，就算 SNR 降至 0 或是 0 以下的時候，我們的加速法則還是可以降低運算時間，這結果顯示出了，就算我們新提出的演算法多做了一個 Normalized Correalator 的運算，其運算時間也不會比其他法則慢。. 35.

(42) 第四章. 實驗結果與數據探討. 4-2 回饋式演算法則的效能評測在評測我們的系統效能前，先來討論 true positive (TP) rate 和 false alarm (FA) rate，這是用來評量棘波偵測演算法的標準之一，TP rate 的定義是偵測正確的棘波個數除以實際上的棘波總數，棘波總數是用模擬器[21]產生棘波序列時就會提供的資訊。FA rate 的定義是誤判的訊號總數除以偵測到的訊號總數(包括誤判的棘波個數跟正確的棘波個數)。數學式如下:. TP rate:. 偵測正確的棘波個數. FA rate:. 棘波總數. 誤判的棘波個數偵測到的訊號個數. 在我們規劃的實驗中，我們會測詴在不同 SNR 情況下棘波序列的 TP rate 和 FA rate。此外，在系統的第一階段跟第二階段的 TP rate 跟 FA rate 是各自獨立的。在第一階段偵測的時候是使用 GLRT detection ，偵測時間是T =2 秒；第二階段是使用 Matched filter detection，偵測時間是T2 =20 秒。. 36.

(43) 第四章. 實驗結果與數據探討. 第一階段 GLRT detection 的執行時間T 訂為 2 秒，是經過仔細評估的。雖然對整個系統的效能而言，系統第一階段的時間越短越好，因為系統越快產生模板就越能增加棘波偵測的效能。但是當T 小於某個時間的話，後續的 Osort 法則很容易因為偵測到的雜訊比棘波還多，而誤把雜訊當成模板回饋給系統使用，造成棘波偵測效能不佳的情況。其不同T 秒數的實驗情況如下表:. T(s). 0.5. 1. 1.5. 2. 2.5. 3. Second. TP Rate. 48.85. 82.12. 88.71. 92.92. 91.24. 93.46. stage. FA Rate. 67.47. 15.35. 6.49. 0.55. 0.72. 0.67. 表格 4-2 當𝐓𝟏 使用不同的秒數時，對第二階段偵測的影響. 從上表可以看出當T = 2以上時，其偵測效果並不沒有明顯提升，而在 T = 2以下時卻會產生 FA rate 大幅提升的問題，其主要原因應該是歸因於當T 太短的時候，Osort 演算法則會產生不符合棘波的模板。Osort 演算法產生的模板效果，如(圖 4-3)、(圖 4-4)所示。從圖中可以看出，雖然經過 Osort 分群後的棘波形狀跟原始的棘波形狀有些出入，但是該模板還是可以提供給系統做有效的棘波偵測。. 37.

(44) 第四章. (a-1) 圖 4-3. (a-2) SNR = 8 時本系統的棘波圖形:. 實驗結果與數據探討. (a-3). (a-1) 實際棘波的平均值。 (a-2) Osort 分類的結果。 (a-3). (b-1). Osort 提供給系統的模板。. (b-2) 圖 4-4. (b-3). SNR = -2 時本系統的棘波圖形:. (b-1) 實際棘波的平均值。 (b-2) Osort 分類的結果。 (b-3) Osort 提供給系統的模板。而關於第一階段的閥值的選擇，我們是根據(6)和(7)這兩個公式訂定了系統 38.

(45) 第四章. 實驗結果與數據探討. 最初的閥值(η=σ2 𝛾)，這裡的 σ 是根據(7)由輸入的棘波序列動態產生，而 γ 的選擇是由[10]提供的γ = 1.2N = 76.8。至於第二階段的主要模板來源是根據 Osort 演算法產生的，此階段因為有經過 Normalized correlator，所以在閥值的選擇上較為直觀，我們選擇 η=0.7(輸入訊號跟模板有七成相似就算棘波)。實驗用的輸入訊號是由模擬器[21]產生內含兩個棘波形狀的棘波序列。其實驗結果如表 3 所示。. SNR(dB). -2. 0. 2. 4. 6. 8. First. TP Rate. 83.37. 87.75. 86.64. 87.03. 90.18. 89.42. Stage. FA Rate. 20.25. 15.17. 9.45. 9.46. 7.26. 5.12. Second. TP Rate. 82.39. 88.22. 89.47. 90.74. 92.92. 93.04. Stage. FA Rate. 0.41. 0.59. 1.06. 0.91. 0.55. 0.36. 表格 4-3 我們提出的演算法在不同 SNR 下的效能. 可以從上表看出，因為第一階段的 GLRT 偵測法則是屬於沒有使用模板的非關聯式的能量偵測，所以當 SNR 低的時候會產生低 TP Rate 和高 FA Rate 這種比較不理想的結果。然而，隨著模板的輔助下，該偵測法則在高雜訊(SNR=-2) 的情況下實現了高 TP rate 和低 FArate，這也證明了模板有助於增進偵測的效能。. 39.

(46) 第四章. 實驗結果與數據探討. 4-3 不同閥值的偵測效果關於我們新提出來的棘波偵測演算法，不只改善了在高雜訊時偵測效能不佳的問題，其優點還有簡化了閥值的選擇。在第二階段偵測的時候因為對輸入訊號做了 normalized correlator，所以閥值的選擇不需要根據不同的輸入源不同而做更改。而是根據模板跟棘波的最大平方距離 D 來判斷是否命中[12]。根據表 3-4 我們可以看出在不同閥值下棘波偵測的效果。. SNR(dB). -2. 0. 2. 4. 6. 8. TP Rate. 76.47 %. 84.63 %. 86.47 %. 87.34 %. 89.18 %. 88.10 %. FA Rate. 0.15 %. 0.52 %. 0.52 %. 0.90 %. 0.23 %. 0.12 %. TP Rate. 79.41 %. 90.66 %. 88.37 %. 90.65 %. 90.09 %. 97.32 %. FA Rate. 0.22 %. 1.92 %. 0.43 %. 0.47 %. 0.26 %. 0.45 %. TP Rate. 82.39 %. 88.22 %. 89.47 %. 90.74 %. 92.92 %. 93.04 %. FA Rate. 0.41 %. 0.59 %. 1.06 %. 0.91 %. 0.55 %. 0.36 %. TP Rate. 90.51 %. 92.18 %. 95.67 %. 96.34 %. 95.03 %. 94.33 %. FA Rate. 22.78 %. 40.26 %. 20.51 %. 14.28 %. 18.72 %. 20.54 %. η=0.8. η=0.75. η=0.7. η=0.65. 表格 4-4 第二階段在不同閥值下的 TP rate 跟 FA rate 表現. 40.

(47) 第四章. 實驗結果與數據探討. 當我們選擇了較大閥值的時候可以降低 FA rate，因為這代表了要當輸入訊號跟模板有更高的關聯度才能判斷為命中，這同時也降低了雜訊被判斷成棘波的機率。反之，當選擇較低閥值的話，雖然在雜訊高的時候 TP rate 會提高，但是也提高了雜訊被判斷成棘波的機率。所以我們選擇了η=0.7，這是為了在各個 SNR 下都能有高 TP rate 跟低 FA rate 的表現。. 理想中的閥值會使得大部分的雜訊被濾掉而僅讓棘波通過。為了說明我們選擇η=0.7 的原因，首先，我們先回顧一下 FA rate 的定義。FA rate 的定義是誤判的訊號總數除以偵測到的訊號總數。分母其實就等校於誤判的棘波個數跟正確判斷的棘波個數的總和(如式 3-1)，所以 FA rate 的值會被誤判的棘波個數或是正確判斷的棘波個數所影響。. FA rate:. 誤判的棘波個數誤判的棘波個數 + 正確判斷的棘波個數. (式 3 − 1). 下面我們假設 g(𝑥 ) 是正規化後的棘波序列區段跟模板內積的最大值，如下式所示:. g(𝐱 ) = max 𝐱̅ 𝑇 𝐭̅𝐢 ≤i≤c. 41. (式 3 − 2).

(48) 第四章. 實驗結果與數據探討. 令 𝑓 是在 SNR=-2 下，雜訊跟模板的 g(𝑥 ) 分布(如圖 3-4)，𝑓2 是在 SNR=-2 下棘波跟模板的 g(𝑥 ) 分布(如圖 3-5)。我們詴著從(圖 3-4)跟(圖 3-5)說明 FA rate 跟閥值η的關係，所以我們令𝐹 (𝜂)是𝑓 在𝜂以上的面積，0<η<1。這時𝐹 (𝜂)是代表當閥值= 𝜂時，所有被誤判為棘波的區段個數，表示如下式:. 𝐹 (𝜂) = ∫ 𝑓 (𝑎)𝑑𝑎. (式 3 − 3). 𝜂. 同理，令𝐹2 (𝜂)是𝑓2 在 𝜂以上的面積，0<η<1，𝐹2 (𝜂)代表當閥值= 𝜂時，所有被正確判斷為棘波的區段個數。所以 FA rate 的公式可以換成以下的表示法:. FA rate:. 𝐹 (𝜂) 𝐹 (𝜂) + 𝐹2 (𝜂). (式 3 − 4). 以 η=0.7 時為例，這時𝐹 (𝜂) = 10、 𝐹2 (𝜂) = 1400、 FA rate = 0.41%。當 η=0.65 時，𝐹 (𝜂) = 430、𝐹2 (𝜂) = 1450、 FA rate = 22%。所以這也解釋了為何在(表 3- 4)中當 η=0.65、η=0.7 時，FA rate 會有巨大的變化的原因。. 42.

(49) 第四章. 圖 4-5 模板跟雜訊的摺積值. 圖 4-6 模板跟棘波的摺積值. 43. 實驗結果與數據探討.

(50) 第四章. 實驗結果與數據探討. (圖 4-5)、(圖 4-6)會有這樣的趨勢的原因，其實主要還是都依賴模板的選擇。若我們選用了一個較不好的模板的話，會導致模板跟棘波的摺積值變低，使得整個值都往 0 的方向靠近，這樣的話就無法有效地分辨棘波跟雜訊了。. 44.

(51) 第四章. 實驗結果與數據探討. 4-4 不同偵測法則的效果比較為了更進一步的評估我們所提出來的演算法，(表 4-5)中比較了其他現行的棘波偵測法則在各種 SNR 下的 TP rate 和 FA rate。實驗用的棘波序列都是用內含兩種棘波形狀的棘波序列，棘波序列長度為 20 秒。. SNR(dB). -3. -2. 0. 1. 8. 10. Proposed. TP Rate. 82.71 %. 82.39 %. 88.22 %. 90.04 %. 93.04 %. 93.64 %. Algorithm. FA Rate. 1.06 %. 0.41 %. 0.59 %. 0.92 %. 0.36 %. 0.40 %. NEO. TP Rate. 80.53 %. 82.05 %. 83.47 %. 87.21 %. 92.24 %. 93.10 %. FA Rate. 57.87 %. 56.16 %. 39.62 %. 22.49 %. 7.75 %. 3.57 %. TP Rate. 86.66 %. 87.38 %. 91.37 %. 92.43 %. 94.26 %. 94.82 %. FA Rate. 82.43 %. 81.80 %. 79.00 %. 79.36 %. 19.00 %. 6.77 %. TP Rate. 22.22 %. 51.08%. 54.65 %. 67.81 %. 81.74 %. 90.90 %. FA Rate. 3.70 %. 0.84 %. 1.52 %. 1.25 %. 12.44 %. 17.96%. TP Rate. 80.31 %. 80.83 %. 82.20 %. 82.90 %. 86.32 %. 89.65 %. FA Rate. 8.92 %. 7.69 %. 4.62 %. 3.02 %. 2.94 %. 2.80 %. Algorithm [4] SWT Algorithm [8]. Absolute Value [10] Matched Algorithm [11]. 表格 4-5 各種棘波偵測法則在不同雜訊等級的棘波序列下的 TP rate 跟 FA rate表現. 45.

(52) 第四章. 實驗結果與數據探討. 在這邊要注意的是，因為我們提出來的演算法在第一個階段是為了產生模板，主要的偵測效能是要觀察第二階段的表現，所以在(表 4-5)中是顯示第二階段的結果。而 matched filter[11] 則是假設我們已經事先準備好需要的模板了。. 從(表 4-5)中可以看出，我們的回饋式演算法透過自適應的模板產生跟簡化的閥值選擇，在各個 SNR 情況下的表現都比其他演算法來的好。就算是跟已經假設有準確模板的 Matched filter 演算法相比，仍然略勝一籌。這是因為 Matched filter 的閥值需要隨著 SNR 跟棘波峰值的變化而做更動，所以要選擇到一個理想的閥值更為困難，而單一的閥值更不可能適用於各種雜訊下的棘波序列。. 雖然以能量偵測為基礎的 NEO[4]和 SWT[8]在高 SNR 的時候有良好的表現，但是隨著 SNR 的降低，雜訊的提高使得這兩個演算法的效能急速惡化。在 SNR-3 下 NEO 和 SWT 的 FA rate 高達 50%。而與此相對的，我們提出的演算法在 SNR=-3 的時候 FA rate 只有 1%而已。這結果證明了我們提出來的演算法，在高雜訊的情況下依然能有不錯的表現。. 46.

(53) 第四章. 實驗結果與數據探討. 4-5 回饋式演算法的分群效果 (表 4-6)和(表 4- 7)顯示了我們提出的演算法中 Osort 演算法的效能，表中紀錄了沒被偵測到的棘波個數和偵測到但是卻分類錯誤的棘波個數，測量用的輸入訊號是各種 SNR 下包含兩種和三種棘波形狀的棘波序列，棘波序列長度為 20 秒。我們可以從(表 4-6) 跟(表 4-7)觀察到，在不同的 SNR 下，棘波序列經過我們提出的演算法後，沒有被偵測到的棘波個數跟誤判的棘波個數都十分的低，這可以有效地證明該系統抵抗雜訊的性能十分強大。. SNR(dB). Total number. Number of. of spikes. spikes Not detected. Number of. Number of spikes. spikes detected& detected& mis-classified correctly-classified. -2. 2378. 406 (17.07 %). 10 (0.42 %). 1962 (82.51 %). 0. 2308. 260 (11.27 %). 11 (0.48 %). 2035 (88.17 %). 2. 2375. 240 (10.10 %). 10 (0.42 %). 2115 (89.05 %). 4. 2387. 197 (8.25 %). 10 (0.42 %). 2180 (91.32 %). 6. 2359. 194 (8.22 %). 9 (0.38 %). 2156 (91.38 %). 8. 2332. 165(7.08%). 7 (0.30 %). 2160 (92.62 %). 表格 4-6 不同 SNR 值下，棘波序列偵測和分群的效果. 47.

(54) 第四章. 實驗結果與數據探討. SNR(dB). Total number of spikes. Number of spikes Not detected. Number of Number of spikes spikes detected& detected& mis-classified correctly-classified. -2. 2014. 354 (17.57 %). 206 (10.23 %). 1454 (72.19 %). 0. 1990. 327 (16.43 %). 194 (9.75 %). 1469 (73.81 %). 2. 1983. 217 (10.94 %). 153 (7.71 %). 1613 (81.34 %). 4. 1988. 227 (11.42 %). 81 (4.07 %). 1680 (84.51 %). 6. 2008. 219 (10.91 %). 166 (8.27 %). 1623 (80.83 %). 8. 1987. 179 (9.00 %). 123 (6.34 %). 1682 (84.65 %). 表格 4-7 三種棘波形狀的棘波序列偵測和分群的效果. 48.

(55) 第四章. 實驗結果與數據探討. 4-6 硬體電路所耗資源本論文提出的棘波偵測法則之所以能夠在不同的 SNR 環境下，有高準確率、低誤判率的優點，主要都是透過計算 Nomalized 及 Correlation ion 來達到此效果。而完成這兩種運算需要大量的加法器及乘法器，這會消耗大量的硬體資源。為了降低硬體資源的消耗，所以我們使用了在前些章節所提到的加速技術，希望可以降低硬體的資源。透過(表 4-8)可以得知電路簡化前後的硬體消耗，可以發現使用了加速技巧簡化後的電路，在硬體資源的消耗上有明顯的減少。. Embedded Logic Elements. Logic Registers. Multipliers. 12262. 4624. 525. 10603. 2257. 106. Normalized Correlator Without Post Normalized Normalized Correlator With Post Normalized. 表格 4-8 電路簡化前後的硬體資源消耗. 以上數據所使用的開發版為 Altera Quartus II 13.0 的 Cyclone IV E : EP4CES115F29C7。因為是在開發版上做實現，我們需要用到開發板上面的 Nios 49.

(56) 第四章. 實驗結果與數據探討. CPU 來驅動我們所設計的電路，所以不可避免的會多些額外的硬體資源消耗。再加上我們的回饋式棘波偵測系統，除了透過上述兩個計算單元之外，還需要提供一個 buffer 給後續的 Osort 去做產生模板的動作。基於以上兩個理由，我們可再透過下表，得到較完整的棘波偵測電路的硬體消耗。. Logic Elements. Logic Registers. Embedded Multipliers. Buffer. 57064. 49595. 489. Nios cpu + on-chip ram + jtag. 3250. 1982. 4. 10603. 2257. 106. 114,480. 114,480. 532. Normalized Correlator Total resources. 表格 4-9 整體電路的消耗. 最後，測量硬體的執行時間，我們期望可以透過把棘波偵測演算法做成硬體，並增加其偵測的效能。假設輸入資料長度為100秒的棘波序列，也就是240萬個採樣點，測詴時的硬體時脈頻率分別是50Mhz、100Mhz、500Mhz、1Ghz、3Ghz。不同時脈頻率所需的計算時間如(表4-10)所示。. 50.

(57) 第四章. Clock Rate. 實驗結果與數據探討. 50Mhz. 100Mhz. 500Mhz. 1G. 3G. 52 ms. 26 ms. 5.2 ms. 2.6 ms. 0.87 ms. Proposed Algorithm. Implemented On FPGA. 表格 4-10 硬體跟軟體執行時間比較. 接著我們可以把硬體執行時間(表4-10)跟軟體執行時間(表4-1)做個比較。從比較兩個表格後可以看出，當輸入的棘波序列SNR=8的時候，就算再用上了 Pre-Screening的加速技巧，軟體的計算時間還是只能達到0.26 sec，還是遠不及硬體的處理速度(無使用Pre-Screening)。. 51.

(58) 第五章. 結論. 第五章結論根據前面的實驗結果，可以證明我們新提出的回饋式棘波偵測系統的偵測效能在任何 SNR 下都是可靠的。而且在我們新提出的演算法中，有對傳統的 Matched filter 加入了 Normalized correlator 的動作且又進一步的使用加速技巧，有效的縮短了計算的時間同時也增加了單位時間的生產量。而也是因為使用了 Normalized correlator 的緣故，我們克服了以往閥值在選擇上的困難處，我們從模板跟棘波序列的平方距離的角度去詮釋了閥值的另一種意義，使得只要決定一個閥值就能對應到各種雜訊程度的棘波序列，特別是在 SNR=-2dB 的時候，我們可以達到高達 82.39%的正確率和只有 0.41%的錯誤率，這明顯優於其他也是利用能量偵測棘波的演算法。而在棘波分類方面，即使在不同的 SNR 下我們提出的演算法也有很高的準確率。以上種種，都證明了回饋式自動模板生成為基礎之正規化關聯值棘波偵測系統可以取代目前所使用的棘波偵測系統或是棘波分類系統。. 52.

(59) 參考. References 1. S. Gibson, J. W. Judy, and D. Markovic, Spike sorting: the first step in decoding the brain, IEEE Signal Processing Magazine, pp.124-143, 2012. 2. M.A. Lebedev and M.A.L. Nicolelis, Brainmachine interfaces: past, present and future, Trends in Neurosciences, Vol. 29, pp.536-546, 2006. 3. M.S., Lewicki, A review of methods for spike sorting: The detection and classification of neural action potentials. Netw. Comput. Neural Syst., Vol. 9, pp. R53R78, 1998. 4. S. Mukhopadhyay and G. C. Ray, A new interpretation of nonlinear energy operator and its efficacy in spike detection, IEEE Trans. Biomed. Eng., Vol. 45, pp. 180187, 1998. 5. I. Obeid and P. D. Wolf, Evaluation of Spike-Detection Algorithms for a Brain-Machine Interface Application, IEEE Trans. Biomed. Eng., Vol. 51, pp. 905-911, 2004. 6. S. Gibson, J. W. Judy, and D. Markovic, Technology-Aware Algorithm Design for Neural Spike Detection, Feature Extraction, and Dimensionality Reduction, IEEE Trans. Neural Systems and Rehabilitation Engineering, Vol. 18, pp.469-478, 2010. 7. K. Oweiss and M. Aghagolzadeh, Detection and classification of extracellular action potential recordings, Chapter 2 of Statistical Signal Processing for Neuroscience, pp.15-74, 2010. 8. K. Kim and S. Kim, A wavelet-based method for action potential detection from extracellular neural signal recording with low signal-to-noise ratio, IEEE Trans. Biomed. Eng., Vol. 50, pp. 999-l011,2003. 9. R. J. Brychta, S. Tuntrakool, M. Appalsamy, N. R. Keller, D. Robertson, R. G. Shiavi, Wavelet Methods for Spike Detection in Mouse Renal Sympathetic Nerve Activity, IEEE Trans. Biomed.Eng., Vol. 54, pp. 82-93, 2007.. 53.

(60) 參考. 10. R. Q. Quiroga, Z. Nadasdy, and Y. Ben-Shaul, Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering, Neural Comp., Vol. 16, pp. 16611687, 2004. 11. N.Mtetwa, L. S. Smith, Smoothing and thresholding in neuronal spike detection, Neurocomputing,Vol. 69, pp. 1366-1370, 2006. 12. T. Sato, T. Suzuki, and K. Mabuchi, Fast Template Matching for Spike Sorting, Electronics and Communications in Japan, Vol. 92, pp.57-63, 2009. 13. S. Kim and J. McNames, Automatic spike detection based on adaptive template matching for extracellular neural recordings, J. Neurosci. Methods, Vol.165, pp.165-174, 2007. 14. K.D. Harris, Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurophysiol., Vol. 84, pp. 401-414, 2000. 15. A. Oliynyk, C. Bonifazzi1, F. Montani, L. Fadiga1, Automatic online spike sorting with singular value decomposition and fuzzy C-mean clustering. BMC Neural Sci., Vol. 13, 2012,doi:10.1186/1471-2202-13-96. 16. W. J. Hwang, W. H. Lee, S. J. Lin and S. Y. Lai, Efficient Architecture for Spike Sorting in Reconfigurable Hardware, Sensors, Vol. 13, pp.14860-14887, 2013. 17. I.T. Jolliffe, Principal Component Analysis, 2nd ed.; Springer: Berlin, Heidelberg, Germany, 2002. 18. S. Miyamoto,H. Ichihashi, and K. Honda, Algorithms for Fuzzy Clustering, Springer:Berlin/Heidelberg, Germany, 2010. 19. U. Rutishauser, Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo, J. Neurosci. Methods, Vol.154, pp.204-224, 2006. 20. J. Wild, Z. Prekopcsak, T. Sieger, D. Novak, and R. Jech, Performance comparison of extracellular spike sorting algorithms for single-channel recordings, J. Neurosci. Methods, Vol.203, pp.369-376,2012. 54.

(61) 參考. 21. L.S. Smith and N. Mtetwa, A tool for synthesizing spike trains with realistic interference. J. Neurosci. Methods Vol. 159, pp. 170-180, 2007. 22. H. Vicent Poor, An Introduction to Signal Detection and Estimation, Springer-Verlag, 1988 23. Kou-Hsuan Wu, Efficient VLSI Architecture for Spike Detection Based on Normalized Correlators,pp,10-17, 2013. 55.

(62)

以回饋式自動模板生成為基礎 之 正規化關聯值棘波偵測系統 之設計及實現

以回饋式自動模板生成為基礎之正規化關聯值棘波偵測系統之設計及實現