第二章 相關研究
此章節為本篇論文之相關研究。在 2.1 章節介紹偵測動作電位 (action potential)及閥值 (threshold)的設計；2.2 章節介紹分析 MER 訊號特徵分析方法；2.3 章節介紹動作電位分群及 對齊方法。
2.1 動作電位偵測
在過去的文獻中，區分 action potential 及背景噪音 (background noise)的方法大多以判斷 訊號的 spike 能量峰值是否超過所設定之 threshold。Threshold 計算方式多數將 MER 訊號切成 等大區段計算區段本身的均方根 (root mean square, RMS)[11]值並乘上一個適當倍率。此方法 雖然簡易且直觀，但 RMS 計算過程中含有平均概念，加上各神經部位放電量程度不相同，易 受到高峰值 spike 影響造成 RMS 值大小變化急遽，使得 threshold 值起伏較大，造成某些 action potential 無法被順利偵測，而某些 background noise 被錯誤偵測。因此希望找到一個穩定且適 當的 threshold，盡量偵測出所有的 action potential 與避免偵測出過多的 background noise。
Thakur et. al [12]提出 cap-fitting 演算法，如圖 2.1 所示。Cap-fitting 方法假設此段訊號中 純噪音(pure noise)能量分佈為高斯分佈 (Gaussian distribution)，利用此概念，將 MER 訊號切 分成數個連續相等長度 frame，計算每一個 frame 的能量分布，找出落在 99.87%分佈區間之 能量值，同時計算此段 frame 訊號之標準差 (standard deviation)。若某段訊號能量為 Gaussian distribution (pure noise)，其 99.87%分佈區間中的能量應約等同 3 倍標準差，根據以上原則可 判斷各 frame 是否為 pure noise，若訊號不屬於 background noise，將此段訊號的 action potential 移除並重複上述過程直到此段訊號滿足 Gaussian distribution。
Harrison et. al [13]提出一套適合嵌入式硬體偵測法 duty-cycle，如圖 2.2 為演算法之流程。
此篇將預設工作週期設定為 15.85%，接著計算每一段訊號的 RMS 值並乘上一個適當倍率來
做為判斷工作週期之閥值，將計算得到的工作週期與其預設工作週期相比較，若此段訊號工 作週期未超過其所設定之工作週期，則認定此段訊號為背景噪音。若此段訊號為背景噪音，
則將重新調整判斷工作週期閥值之倍率，並重新計算此段訊號工作週期直到符合預設之工作 週期。反之若訊號工作週期符合預設值，則計算此段訊號之標準差且搭配一個適當的倍率即 可做為偵測動作電位閥值。
圖 2.1. Cap-fitting 法之流程圖
圖 2.2. Duty-cycle 演算法之流程圖
Chan et. al [14-15]提出了 (max-min spread, MMS)方法，希望能降低運算時間，如圖 2.3 所 示。將訊號切分成相等長度前後 50%重疊的 frame，將每個 frame 之訊號能量最大值(max value) 與最小值(min value)相減，即可獲得該段 frame 之 MMS 值。將連續 200 個 MMS 值存入 buffer 中，並且由小到大排序並給予一個適當的α值作為擷取 MMS 值的範圍，預設的α值採用 0.6，
最後將排序過第α 5⁄ 至1 − α (即第 24 個至 80 個) MMS 值加總平均並乘上適當的倍率即為 threshold。，結果部分此方法可以得到一個最平滑的閥值並在利用高斯噪音與顏色噪音所模擬 的訊號中偵測率可以達到 80%以上。
俞仲麒 [16]提出混合 RMS 與 duty-cycle 的方法，希望能得到較穩定的 threshold。此方法 將訊號切分成連續相等長度的 frame 並將各個 frame RMS*2.5 為 threshold，計算各 frame 的 duty-cycle 若 duty-cycle 小於 1.5%則標記此 frame 為背景噪音。最後串聯所有標記為背景噪音 的 frame，計算其 RMS 做為偵測 action potential 之 threshold，利用此方法可以得到一個固定 的閥值此閥值不會隨著不同深度中訊號的變化而改變，本研究計算 threshold 方式與此方法相
圖 2.3. MMS 演算法之流程圖
2.2 特徵分析
MER 訊號特徵方析大致上可分成兩類，第一類針對 MER 時域部分進行特徵分析，第二 類針對 MER 頻域部分進行特徵分析，最終目標都是希望透過不同的分析方法量化數據及結 果，用來協助外科醫師在手術中判斷 MER 訊號來自大腦中哪個區域，以便將電極植入預定目 標。
Simona et. al [17]提出在 STN、SNr 或其他區域中不同的位置特徵數值會有差異，並分析 在不同性別特徵數值是否依然存有差異。此方法在對誤置進行分析時可以得知在 STN 中 number of spike 與 average nonlinear energy 較高、rising time 與 duration time 較低。在對性別 分析時可得知男生有較高的 average nonlinear energy 與 rising time，女生有較高的 number of spike。以下為六種特徵值
(1) Number of spike：平均每秒鐘出現的 action potential 次數。
(2) Rising time：如圖 2.4 (a)所示，為某 action potential 起始點至 amplitude 最高點經過時間(a 點至 b 點)。
(3) Duration time：如圖 2.4 (b)所示，為某 action potential 起始點至結束點經過時間(c 點至 d 點)。
(4) Average nonlinear energy：定義為可將振幅及頻率轉變為單一變量：
E(n)= 1
N-2∑ xi2-xi-1*xi+1
N-1
i=2
(5) Peaks
1
2∑ max{0, 𝑠𝑔𝑛[𝑥𝑖+2− 𝑥𝑖+1] − 𝑠𝑔𝑛[𝑥𝑖+1− 𝑥𝑖]}
𝑛−2
𝑖=1
(6) Zero crossings
1
2∑|𝑠𝑔𝑛(𝑥𝑖+1) − 𝑠𝑔𝑛(𝑥𝑖)|
𝑛−1
𝑖=1
黃麒璋 [18]沿襲上段所使用的相同特徵值並先利用經驗模態分解法(empirical mode decomposition, EMD)分解 MER 訊號，接著利用本質函數(intrinsic mode function, IMF)觀察訊 號特徵，而結果 number of spike 在 STN 中與其他區域有較明顯的提升而 rising time 在 STN 中 較明顯的穩定，訊號利用經驗模態分解法降躁後計算的特徵值其特徵值隨不同位置變化更為 明顯，因此提高辨識的正確率。
Ciecierski et. al [19]與 Novak et. al [20-21]相較於以上文獻皆利用時域進行特徵值計算，
提出頻域特徵值計算方法。Ciecierski et. al 利用 MER 訊號在經過雜訊抑制號後 1 ~ 100 Hz 的 各個深度的頻譜能量加總，利用特徵值在 STN 中會出現明顯的提高進行 STN 位置的判斷，
利用此方法不僅可以判斷 STN 位置也可以判斷較好植入路徑，此方法輔助進行路徑挑選其結 果在 20 位病人中，有 71%可以挑選到更好的路徑。
(a) (b)
圖 2.4. 時域特徵值示意圖(a)動作電位之 rising time 及 rising time area 示意圖 (b) 動作電位之 duration 及 duration area 示意圖
Novak et. al [20-21]提出利用 MER 訊號中 multiunit activity (MER 訊號中 500 Hz ~ 2000 Hz 部分)的頻譜能量密度在 STN 中會出現明顯的提高，利用此方法與醫生判斷進行比較其結果 在 10 位病人 50 段訊號中判斷進入 STN 深度誤差為 0.31±0.84mm 離開 STN 深度誤差為 0.44±0.47 mm，彼此之間差異性只有 5.6%相關係數為 0.79，與醫生在術中臨床判斷相當吻合。
Wong et. al [22]提出為先分別對於微電極訊號以及動作電位進行特徵計算再利用多種特 徵在多維空間中 k-mean clustering 的歐基里德距離，再將其繪製出 activity map，利用 activity map 不僅可以判斷出整段訊號中 STN 的位置，更可區分不同深度錄製的訊號在何處神經元區 域中，如圖 2.5 所示，此結果在 39 位病人中判斷 STN 上緣邊界深度誤差為 0.0422 mm ± 0.3164 mm，STN 下緣邊界深度誤差為 0.0422 mm ± 0.3164 mm，有 95%以上誤差小於 1mm。
Gemmar et. al [23]提出利用三個層級對 STN 進行判斷，每個層級都會用不同的判斷方法 來判斷 STN，第一層先利用各 frame 點中中位數與平均的比例判斷是否有背景活動的增加並 標記出疑似 STN 的區域，第二層在將訊號進行小波分解計計算出典型 STN 會出現的不規則 變化進行疑似 STN 的標記，最後在第三層對於前兩層的結論進行整合得到確切的 STN 位置，
此方法與手術紀錄相比有 95%的正確性。
圖 2.5. Activity map 判斷之結果
2.3 動作電位分群對齊
MER 訊號中在錄製時可能會因為時間偏移與噪音等外在因素，導致同一神經元的動作電 位波形有差異。因此 Chan et. al[24]提出將動作電位(action potential)轉換至相空間(phase space) 的方式，將原本一維空間的 action potential 轉換至二維空間，再從對應的 phase space 取出其 輪廓能量值 (major portrait radius, MPR)來代表其轉換過後的 action potential，希望可以降低 MER 訊號在錄製過程中可能因為 noise 或其他因素造成 action potential 的偏移，得到更好的 分群結果。如圖 2.6 所示，利用此方法可以放大動作電位波形的差異性，讓原本僅有些微差異 的波形有較明顯的區隔，使 k-mean clustering 分群後歸類在同一群集內的動作電位波形相似 度更高。使用 k-mean clustering 演算法之前會用主成份分析(principal component analysis, PCA) 降低 MPR 降維，藉此降低分群動作的時間複雜度。
圖 2.6. 一維 action potential 轉換至二維與輪廓能量值示意圖
Chan et.al [25]為了更方便分析動作電位波形，提出將所有擷取出的動作電位依固定參考 點對齊 (alignment) 的方式。其他文獻也曾提出同樣的概念，分別計算動作電位波形極值、微 分或斜率等方式決定對齊的參考點。然而本篇作者認為不同型態的動作電位應該使用不同的 對齊參考點，因此提出一個選擇的流程，在每個動作電位的起點 (P1)、終點 (P2) 或最高點 (D) 之中選擇一個做為對齊參考點。選擇流程如圖 2.7 所示，對於每個擷取出的動作電位分別 計算 P1、P2 及 D 三個點的 IMTE (instantaneous multi-resolution teager energy, IMTE) 值，依序 標記為 AP1、AP2及 AD，再依 AP1、AP2及 AD值的大小判斷動位電位型態並選擇合適的對齊參 考點。利用此方法可將不明確的對齊比率降低至約為 0.08~0.11%，而對齊錯誤的比率降低至 約為 0~0.02%。
圖 2.7. 對齊參考點選擇流程圖
第三章 多維特徵分析
本章節將詳細介紹多維特徵分析法(analysis of multi-feature)。並利用多維特徵分析法建立 出 activity map 當此段訊號特徵值與 STN 的特徵值出現較為一致性的變化時，其顏色比例會 最大，因此就可以利用顏色的比例對 MER 訊號進行 STN 的判斷。多維特徵分析可同時分析 多個 MER 訊號中 action potential 的特徵值，特徵值分析利用在不同植入位置特徵值的變化進 行位置的判斷。3.1 章節介紹 MER 訊號擷取組成方式；3.2 章節介紹特徵值計算方式及閥值 偵測設計；3.3 章節介紹多維特徵分析及其流程；3.4 章節則為多維特徵值分析結果說明及討 論。
3.1 MER 訊號之組成方法
目前針對帕金森氏症 DBS 手術是公認有效治療的方法之一，在 DBS 手術中須將 lead 植 入目標位置，而目前 lead 主要植入目標位置為 STN，在進行 lead 植入時需要對目標位置的精 準定位，因此會錄製 MER 訊號供醫生做手術中的判斷，因此本研究希望分析 MER 訊號中的 特徵進行 STN 區域的判別。
本研究使用的 MER 訊號由中國醫藥大學提供，於 DBS 手術中採用 Medtronic Leadpoint 系統錄製，取樣頻率為 24KHz。DBS 手術前會以 stereotaxic instrument 軟體決定初始植入 trajectory 及頭骨上的插入點、深度及角度，之後將 probe 沿初始植入 trajectory，植入至預計 深度，沿途測量錄製 MER 訊號。
本研究選用 5 位病患區分左右部分共 10 段 trajectory 的 MER 訊號，去除部分醫生術中測 試反應訊號後再依照深度將錄製的訊號串聯組合成，再以低通濾波器去除低頻雜訊如圖 3.1，
本研究選用 5 位病患區分左右部分共 10 段 trajectory 的 MER 訊號，去除部分醫生術中測 試反應訊號後再依照深度將錄製的訊號串聯組合成，再以低通濾波器去除低頻雜訊如圖 3.1，