針對用於遠距醫療的量測裝置之心電圖特徵描繪演算法開發

全文

(1)國立台灣師範大學機電工程學系碩士論文. 指導教授：吳順德博士針對用於遠距醫療的量測裝置之心電圖特徵描繪演算法開發 Development of ECG delineation algorithm for measurement devices used in telemedicine. 研究生：林昆宏中. 華. 民. 國. 1. 0. 5. 撰年. 7. 月.

(2) 摘要隨著遠距醫療的發展，在過去的幾十年中，已有許多用於心電圖量測的可攜式裝置被研發出來。這些科技進步的成果使我們能夠在居家環境中進行心電圖量測，甚至讓操作變得更方便。然而，在這種模式下，我們無法保證訊號品質的良好。因為透過可攜式裝置量測而得的心電訊號，容易因為電極黏貼不穩造成嚴重的電極接觸雜訊干擾。此外，遠距醫療模式也使得資料可能蒐集自不同的國家，因而造成資料受到不同的電力線干擾。這些問題將導致許多心電圖特徵描繪演算法無法正常運作甚至失效。因此，本研究針對上述問題提出了一套新的特徵描繪演算法，其中包含： 1. 提出短時 Goertzel 演算法來偵測電力線頻率，並以零相位延遲的 Butterworth 陷波濾波器來處理電力線干擾。 2. 使用截止頻率為 0.8Hz 的零相位延遲 Butterworth 高通濾波器來處理心電訊號的基準線飄移。 3. 提出一個基於粗粒化流程，並配合微分、平方與移動平均的方法來處理電極接觸雜訊。 4. 改良 Pan-Tompkins 演算法並將其用於偵測心電訊號中的 R 波。 5. 改良 Elgendi 演算法並將其用來偵測心電訊號中的 P、T 波。最後，本研究提出的這套特徵描繪演算法已被用於偵測 CMATE® 心電圖量測裝置所量測的資料中的心電訊號特徵。而實驗結果也證明了本研究提出的演算法的可行性與不錯的表現。. 關鍵字：遠距醫療、心電圖、特徵描繪演算法. i.

(3) Abstract With the growth of telemedicine, kinds of portable devices for ECG measurement have been developed over the past few decades. It makes in-home ECG measurements possible and convenient. However, a good signal quality could not be always guaranteed. ECG signals measured by a portable device may be contaminated by serious electrode contact noise (ECN) due to misconnection of electrode. Besides, for ECG signals collecting from different countries, signals could be contaminated by power line interferences with different frequencies. These drawbacks make lots of conventional delineation algorithms fail to work. In this study, we propose a new delineation algorithm to overcome this difficulty. This algorithm consists of the following parts: 1.. A short-time Goertzel algorithm was proposed to detect the frequency of the power line interference, and then a zero phase Butterworth notch filter is designed to attenuate the power line interference.. 2.. The baseline wandering is filtered by using a zero phase Butterworth high pass filter with cutoff frequency at 0.8 Hz.. 3.. A novel algorithm based on coarse grain, derivative, squaring and moving averaging is developed to handle the electrode contact noise.. 4.. A modified Pan-Tompkins’ algorithm is used to detect R peaks of ECG signals. 5.. Finally, P and T waves are detected by modified Elgendi’s algorithm. The proposed algorithm has been implemented to delineate ECG signals. measured by CMATE® ECG devices successfully. Experimental results showed the feasibility and superiority of the proposed algorithm. Keyword: Telemedicine, electrocardiogram, delineation algorithm ii.

(4) 誌謝在碩士班的這段期間，非常感謝我的指導教授吳順德博士。由於在大學時期曾擔任過老師的課程助教，能明顯地感受到老師在教學與研究上的認真。而在進入碩班後，老師在研究方面多以指引方向為主，其餘皆放手讓我不斷地從實務中學習，使我著實體會到做研究該有的態度：勇於面對問題、接受磨練。另外也感謝系上的葉榮木教授、呂有勝教授、張國維教授與鄭淳護教授，在大學與碩班的這幾年間刺激了我的成長。在口試當天，特別感謝海洋大學林修國教授、台北科技大學劉益宏教授與逸奇科技王逸民博士，能夠抽空指導，提供建議與未來的研究方向，讓本研究可以更加完備。而本論文能夠有這樣的成果，還要感謝智創管理顧問公司賴冠吉總經理的協助。若沒有賴總經理提供的醫學知識、資料庫與技術支援，本研究是無法有如此的收穫的。而這兩年在訊號處理實驗室中，非常感謝我的好夥伴鈞翔。除了在課業與研究上的幫忙，你也在我面對一些重要抉擇時提供了另一套想法，讓我能夠重新審視。此外，也謝謝求文學長在這兩年裡的照顧，雖然不常見面，但您提供的人生經歷是寶貴的。還有，從大學時期就是同學的登元、彥廷、典勇、揚群、禹承、伯聰、姿穎、文乙及永宏、家顥學長，謝謝你們分享在自己領域上的研究想法與經驗，讓我接觸的領域不被侷限。最後感謝我的家人和異父異母的兄弟奕德，因為一路上有你們的支持與分擔，才能讓我即使跌跌撞撞地也能走到今天，萬分感激。 2016/07. iii.

(5) 目錄摘要 ...................................................................................................................... i Abstract ............................................................................................................... ii 誌謝 .................................................................................................................... iii 目錄 .................................................................................................................... iv 表目錄 ................................................................................................................ vi 圖目錄 ............................................................................................................... vii 第一章緒論 ....................................................................................................... 1 1-1 前言 ...................................................................................................... 1 1-2 研究動機與目的 .................................................................................. 2 1-3 論文架構 .............................................................................................. 4 1-4 文獻回顧 .............................................................................................. 5 1-4-1 電力線干擾處理 ....................................................................... 5 1-4-2 電極接觸雜訊處理 ................................................................... 6 1-4-3 QRS 波群偵測 ........................................................................... 6 1-4-4 P、T 波偵測 .............................................................................. 7 第二章演算法設計－前處理 ........................................................................... 8 2-1 電力線干擾頻率偵測 .......................................................................... 8 2-1-1 Goertzel 演算法 ......................................................................... 9 2-1-2 短時 Goertzel 演算法(Short-time Goertzel algorithm) .......... 14 2-2 電極接觸雜訊處理 ............................................................................ 17 2-2-1 電極接觸雜訊的特徵 ............................................................. 18 2-2-2 前處理 ..................................................................................... 19 2-2-3 粗粒化(Coarse-graining approach) ......................................... 20. iv.

(6) 2-2-4 切入點偵測(Cut-in point detection) ....................................... 22 2-2-5 正常訊號片段選擇 ................................................................. 25 2-3 濾波 .................................................................................................... 26 第三章演算法設計－特徵偵測..................................................................... 27 3-1 QRS 波群偵測 .................................................................................... 27 3-1-1 R 波偵測................................................................................... 28 3-1-2 誤判 R 波的處理 .................................................................... 43 3-1-3 Q、S 波偵測 ............................................................................ 45 3-2 P、T 波偵測 ....................................................................................... 46 3-2-1 Elgendi 演算法 ......................................................................... 46 3-2-2 P、T 波檢驗 ............................................................................ 53 第四章效能分析與偵測結果討論................................................................. 56 4-1 效能分析：短時 Goertzel 演算法 .................................................... 56 4-2 結果討論：QRS 波群偵測 ............................................................... 59 4-3 結果討論：P、T 波偵測 .................................................................. 63 第五章結論 ..................................................................................................... 66 5-1 結論 .................................................................................................... 66 5-2 本研究之具體貢獻 ............................................................................ 67 5-3 未來展望 ............................................................................................ 67 5-3-1 短時 Goertzel 演算法 ............................................................. 67 5-3-2 誤判 R 波的處理 .................................................................... 68 5-3-3 P、T 波檢驗 ............................................................................ 68 參考文獻 ........................................................................................................... 69. v.

(7) 表目錄表 2-1. Goertzel 演算法的虛擬程式碼 ........................................................... 13 表 2-2. 短時 Goertzel 演算法的虛擬程式碼 ................................................. 17 表 3-1. Pan-Tompkins 決策演算法中各符號之意義 ...................................... 36. vi.

(8) 圖目錄圖 1-1. 心電訊號與心臟組織活動之關係圖 ................................................... 2 圖 1-2. 心電圖處理上常見的干擾 ................................................................... 3 圖 1-2(續). 心電圖處理上常見的干擾............................................................. 4 圖 1-3. 特徵描繪演算法流程圖 ....................................................................... 4 圖 1-4. db3 小波 ................................................................................................. 6 圖 2-1. 電力線干擾 ........................................................................................... 9 圖 2-2. 供電頻率的變化 ................................................................................. 14 圖 2-3. 電力線干擾在頻譜上的分佈 ............................................................. 15 圖 2-4. 電極接觸雜訊的特徵 ......................................................................... 19 圖 2-5. 濾波後的電極接觸雜訊 ..................................................................... 20 圖 2-6. 粗粒化流程示意圖 ............................................................................. 21 圖 2-7. 切入點偵測（暫態峰） ..................................................................... 23 圖 2-8. 粗粒化前後的心電訊號振幅分佈直方圖 ......................................... 24 圖 2-9. 正常訊號片段選擇的處理流程 ......................................................... 25 圖 2-10. 心電訊號特徵波的頻率分佈 ........................................................... 26 圖 3-1. QRS 波群偵測的流程圖 ..................................................................... 28 圖 3-2. Pan-Tompkins 演算法的流程圖 .......................................................... 29 圖 3-3. Pan-Tompkins 演算法各步驟的處理結果 .......................................... 30 圖 3-3（續）. Pan-Tompkins 演算法各步驟的處理結果 .............................. 31 圖 3-4. 五點差分濾波器的頻率響應 ............................................................. 33 圖 3-5. QRS 波群與移動平均濾波器窗口的位置關係 ................................. 34 圖 3-6. Pan-Tompkins 演算法中閥值過濾的流程圖 ...................................... 37 圖 3-7. 含有 VPC 的心電訊號 ....................................................................... 44. vii.

(9) 圖 3-8. Elgendi 演算法的流程圖 ..................................................................... 47 圖 3-9. Elgendi 演算法的處理結果 ................................................................. 48 圖 3-9（續）. Elgendi 演算法的處理結果 ..................................................... 49 圖 3-10. P、T 波與 RR 區間的位置關係 ....................................................... 52 圖 3-11. 心律產生器的心跳的偵測結果 ....................................................... 53 圖 3-12. 心房顫動的偵測結果 ....................................................................... 54 圖 3-13. 倒置 P 波的偵測結果 ....................................................................... 54 圖 3-14. P、T 波檢驗示意圖 ........................................................................... 55 圖 4-1. 效能分析(int-32bits) ........................................................................... 57 圖 4-2. 效能分析(float-32bits) ........................................................................ 58 圖 4-3. 偵測結果比較（平坦線） ................................................................. 59 圖 4-4. 偵測結果比較（暫態峰） ................................................................. 60 圖 4-5. 受電極接觸雜訊干擾的訊號 ............................................................. 61 圖 4-6. 被過度拉高的閥值 ............................................................................. 61 圖 4-7. 因過高閥值無法恢復的偵測結果 ..................................................... 62 圖 4-8. 電極接觸雜訊處理與 Pan-Tompkins 演算法配合的結果 ................ 62 圖 4-9. 因 VPC 而產生的錯誤判斷 ............................................................... 63 圖 4-10. 經 k-means 分群演算法處理後 ........................................................ 63 圖 4-11. 心律產生器的心跳的偵測結果（對照圖 3-10） ........................... 63 圖 4-12. 心房顫動的偵測結果（對照圖 3-11） ........................................... 64 圖 4-13. 心房顫動資料中異常區間的濾波與曲線擬合結果 ....................... 64 圖 4-14. 倒置 P 波的偵測結果（對照圖 3-12）........................................... 65 圖 4-15. 對倒置 P 波定義感興趣區間(BOI)的結果（對照圖 3-12） ......... 65. viii.

(10) 第一章緒論 1-1 前言心電圖(Electrocardiography, ECG/EKG)是一項用來記錄心臟組織的電生理活動(Electrophysiological activity)的技術，其透過肢體上成對的電極來紀錄心臟組織因去極化(Depolarization)與再極化(Repolarization)的過程與節律（見圖 1-1）。因此我們能透過此技術來檢測異常的心臟節律與大部分因組織受損或瘀血造成的異常活動。過去人們有檢測心臟健康狀況的需求時，都是在醫院接受專業人士的處理與測量。而現在科技進步，為了因應健康監測的需求與即時性，不少新興的硬體與服務開始崛起，同時也促進了遠距醫療(Telemedicine)的發展。在這種情況下，因為便利性提昇與成本下降，長時間生理監測的需求也會有成長的趨勢。然而，這樣的趨勢也代表資料量增加。過去心電圖必須經由醫師或醫檢師判讀，而長時間的監測資料若沒有軟體的協助，僅由人工判讀會是一項重擔。雖然已有眾多學者在心電圖特徵的自動判讀演算法上努力，但是不見得適用於現今透過遠距醫療蒐集的心電圖資料。因為過去的資料庫組成較單純，通常只由同一家醫院提供，並有醫檢師把關品質與分類（如 Physionet 的 Physiobank 下各個資料庫[1]）。而在遠距醫療的模式下，使用者一方通常欠缺專業人士的輔助與乾淨穩定的量測場所，所以在資料上更容易見到干擾，許多特徵偵測演算法也容易因此失效。因此，過去在處理受干擾的資料時，由於通常監測時間短，所以可以透過干擾的嚴重性與影響時間長短來決定要丟棄資料重新量測或剪接，再使用演算法來找出特徵。但是對於長時間監測的資料，受干擾的片段可能會分散在許多時間點，此時不僅難以決定是否丟棄重新量測，連剪接出有用的片段 1.

(11) 也是一項重擔。所以，本研究將逐步探討在遠距醫療模式下，哪些訊號處理方面的問題是有別於以往在醫院量測時更需被注意的，並且提出一套實際可被應用的特徵偵測演算法。. 圖 1-1. 心電訊號與心臟組織活動之關係圖[1]. 1-2 研究動機與目的在心電圖的處理上，常見的干擾有基準線漂移(Baseline wandering)[3]、電力線干擾(Power line interference)[4]、動作干擾(Motion artifact)[5]、肌肉震顫(Muscle tremor)[6]及電極接觸雜訊(Electrode contact noise)…等。（見圖 12）上述的幾項干擾，至今已有許多相關訊號處理的研究被發表於世。然而在電力線干擾的部份，多數研究是在已知干擾頻率的情況下設計演算法。但在遠距醫療的模式下，資料有可能是從世界各地而來的，受到干擾的頻率不盡相同。除非倚賴額外的資料註記，否則許多演算法是無法使用的。. 2.

(12) 而在電極接觸雜訊的部份，歷年來主要以硬體方面的改良為主，鮮少有訊號處理相關的研究。本研究推測，由於導致電極接觸雜訊產生的原因是電極與皮膚的不良接觸，而且在實務上通常被視為操作不當而造成的干擾，因此解決方案是以硬體改良為主。但是，硬體的改良雖然可以減少電極接觸雜訊的影響，卻難以避免其發生，而且電極接觸雜訊往往是導致特徵偵測演算法失效的元兇。因此，在軟體上也應該要有應對的策略。此外，如 1-1 節所述，由於遠距醫療與過往醫療模式的差異，資料的長度會影響重新量測的成本與處理負擔，而且數量也會隨著需求上升而增加。所以面對遠距醫療的發展，除了要解決上述兩項問題，也應該要有一套完善的心電圖描繪演算法(ECG delineation algorithm)來協助醫師與醫檢師快速找到特徵波，以便於判斷是否有異常。而本研究將針對上述幾項問題，提出對應的前處理方法，並且分別對心電圖的 QRS 波群、P 波和 T 波以及 P 波起點和 T 波終點這幾項特徵提出偵測演算法，並組成一套完整的特徵描繪流程。. (a) 基準線漂移. (b) 電力線干擾. (c) 動作干擾 (d) 肌肉震顫圖 1-2. 心電圖處理上常見的干擾 3.

(13) (e) 電極接觸雜訊圖 1-2(續). 心電圖處理上常見的干擾. 1-3 論文架構本論文分為五個章節。第一章說明在遠距醫療模式下處理心電圖會遇到的困難，與回顧過去心電圖特徵偵測演算法；第二章將說明用來處理前述問題的訊號前處理(Preprocessing)演算法；第三章是說明組成特徵描繪流程的各個演算法流程；而在第四章將以效能分析與偵測結果來檢驗本研究的演算法成果；最後第五章將說明本論文之貢獻與未來需再進一步研究的部分。而本論文提出的特徵描繪演算法之流程如下：. 圖 1-3. 特徵描繪演算法流程圖. 4.

(14) 1-4 文獻回顧本研究的處理重點有電力線干擾、電極接觸雜訊、QRS 波群偵測、P 波與 T 波偵測。其中前兩項是為了因應在遠距醫療模式的發展下，可能對心電圖量測方面造成的問題；而後面兩項則是為了建立一套自動化的特徵偵測演算法，以因應受遠距醫療發展而帶動的心電圖量測需求。. 1-4-1 電力線干擾處理電力線干擾是指在擷取訊號的過程中，受到周圍的交流供電線路影響，在量測設備的電路或線路上因電磁感應而產生的電流干擾[4]。電力線干擾的特徵為其通常具有極強的基頻，視國家或地區的不同，目前主要可分為 50Hz 與 60Hz。而由於電力線為單一頻率的干擾，在許多電路設計或訊號處理上都會使用陷波濾波器(Notch filter)來處理。但是陷波濾波器的缺點是其對於那些頻率很靠近欲濾除的中心頻率的訊號也會有很大的影響，而且對於較高頻的訊號還會有嚴重的振鈴效應(Ringing effect)[7]。因此，有人提出適應性(Adaptive)演算法來做處理[4][8]。但適應性演算法通常需要良好的初始條件，而且濾波效果也受其收斂速度影響。而為了避免振鈴效應對訊號的影響與顧及訊號的處理速度，也有人提出透過線性準則(Linear criterion)判斷電力線干擾的相位，再利用減法流程 (Subtraction procedure)來消除電力線干擾[9]。但以上的方法都是建立於已知電力線干擾頻率的情況下，若無法得知頻率，則相關的處理參數便無法確定。. 5.

(15) 1-4-2 電極接觸雜訊處理關於電極接觸雜訊(Electrode contact noise)的文獻，情況已在 1-2 節敘述。由於與訊號處理有關的文獻僅以說明電極接觸雜訊的成因與特徵[10][11]，所以本論文在後續的演算法設計只能以這些文獻與實際處理心電訊號時的經驗來說明。. 1-4-3 QRS 波群偵測 QRS 波群的偵測，通常以 R 波為主要目標。因為 Q 波與 S 波分別為 R 波起點(Onset)與終點(Offset)，所以只要先找到 R 波，Q 波與 S 波的偵測就不會是太大的難題。在許多 QRS 波群偵測的演算法中，以 J. Pan 與 W.J. Tompkins 兩人所著的演算法[12]最為廣傳。此演算法會先利用濾波、微分、平方、移動平均這四個操作框選出 QRS 波群大概的位置，接著找出各個方框的最大值，最後以逐步迭代的閥值判斷與更新來找出 R 波的實際位置。除此之外，也有人利用模板比對(Template matching)的方式來尋找 QRS 波群。這類方法最為知名的是使用小波轉換(Wavelet transform)來進行，因為其中的 db3 小波(Daubechies 3 wavelet，見圖 1-4)外型與 QRS 波群非常相似，所以許多以小波轉換為基礎的偵測演算法會利用此小波來偵測[13]。. 圖 1-4. db3 小波[14] 6.

(16) 1-4-4 P、T 波偵測 P、T 波相較於 QRS 波群是屬於較低頻的心電圖特徵，因此多數偵測演算法會選擇先將 QRS 波群找出，再進行 P、T 波的偵測。以 S.S. Mehta 與 N.S. Lingayat 提出的方法[15]為例，其在移除 QRS 波群後，算出 12 個導程訊號中的各點之間的斜率並正規化 (Normalize) ，以之作為支持向量機 (Support vector machine)的訓練資料，再配合一個滑動窗口(Sliding window) 找出 P、T 波的位置。而因為 P、T 波在一個 RR 區間(RR interval，指兩個相鄰 R 波之間的區間)內是有相對的位置關係，因此 Chao Lin 等人提出一個以貝氏階層模型 (Bayesian hierarchical model)為基礎、使用 Gibbs 取樣程序(Gibbs sampling)的方法，在 RR 區間中分別定出 P、T 波可能出現的位置，再從其中找到 P、T 波[16]。近年來，有一位學者 Mohamed Elgendi 則是提出以移動平均交叉(Moving average crossover)的方式來尋找 P 波與 T 波[32]。這項方法與過往方法不同的地方在於，其能夠同時找出一個 RR 區間內可能有 P 波與 T 波的位置，因此只要再配上一些生理上的相關參數，便能判斷該 RR 區間是否有正常的 P 波與 T 波。. 7.

(17) 第二章演算法設計－前處理在訊號前處理的部分，本研究將針對電力線干擾 (Power line interference/Mains hum)與電極接觸雜訊(Electrode contact noise)作討論。原因如 1-2 節所述，電力線干擾處理的部分已有許多演算法被提出來，但本研究的重點在於偵測電力線干擾的頻率；而電極接觸雜訊則是因為在訊號處理領域幾乎沒有相關的研究發表，所以本研究提出一套演算法來處理，且可以配合後續特徵偵測的演算法使用。. 2-1 電力線干擾頻率偵測電力線干擾是生醫訊號處理(Biomedical signal processing)中的一個重要議題。在現代的都市生活模式下，許多測量設備都是電氣化的，而且受測者通常也都處在有電力線佈置的空間中，所以透過設備擷取的訊號通常都會受到電力線的干擾。因此，標準的訊號擷取設備都會加上處理這種干擾的電路與屏蔽設計，以免強度過大的干擾蓋過我們所要擷取的訊號。但在實際測量時，仍可能會有其他外在的因素，造成電力線干擾在經過訊號擷取設備的處理後，依然在擷取的訊號中佔有不小的成份。因此，通常要對生醫訊號的原始資料(Raw data)進行分析前，都會進行數位濾波，讓電力線造成的干擾降到最低，以利後續處理。要處理電力線干擾，可以使用陷波濾波器(Notch filter)來針對特定的窄頻寬進行濾波，但前提是必須知道干擾的頻率為何。以一般居家環境來說，電力線干擾的頻率主要有 50Hz 與 60Hz 兩種，其與各個國家或地區的發電機組運轉的頻率有關。而要判斷訊號所受的干擾頻率為何，通常我們可以使用快速傅立葉轉換 (Fast Fourier Transform, FFT)來做頻率分析。根據 Nyquist-Shannon 取樣定理， 8.

(18) 透過 FFT 分析的有效頻率項是從 0Hz 到. f s / 2 Hz。但我們所關注的只有. 50Hz. 與 60Hz 這兩項，在這樣的目的下，使用 FFT 來分辨電力線干擾頻率並不符合效益。. 圖 2-1. 電力線干擾. 2-1-1 Goertzel 演算法針對這樣的議題，我們可以使用 Goertzel 演算法[17]來分析，這個演算法在只需要求出少數離散傅立葉轉換(Discrete Fourier Transform, DFT)項時特別有效率。以 DFT 的概念來說，要對一訊號作頻率分析，就是求出訊號在各種不同頻率的基底上的投影量。若訊號 x[n]長度為 N，則以數學式表示 DFT 的結果 X[k]為 X [k ] . N -1.  x[n]e. n0. i. 2 nk N. , k  0, 1, 2 ... f s. 將(2-1)式改寫為摺積(Convolution)的形式[18]：. 9. (2-1).

(19) X [k ]. . N -1. . x[n]WNnk. , WN  e. i. 2 N. n0. . N -1.  x[n]WN k ( N  n) ,. WNkN  1. n 0. . N -1.  x[r ]WN k ( N  r ). n 0.   n N  x[n]  WN nk u[n]  n N.  x[n]  WNnk. (2-2). 由上方(2-2)式可以看出，DFT 即為訊號 x[n]與一個線性非時變的(Linear time-invariant, LTI)濾波器 WNnk u[n] 的摺積。而若將這個 LTI 濾波器以 z 轉換的型式改寫則為 . WN mk z  m. H ( z) . (2-3). m0. 根據(2-3)式，此濾波器的轉移函數為 H (z ). . . WN mk z  m. m0. . 1  WN k z 1. . . WN mk z  m  k 1 1  WN z m  0 . . . . WN mk z  m m0. WN0 z 0. 1  WN k z 1. . . WN(m 1)k z (m 1). m0 1  WN k z 1. . 1 1  WN k z 1. 因此，輸入訊號 x[n]與輸出訊號 yk[n]的關係可以表示為. 10. (2-4).

(20) yk [n]  x[n]  WNk y[n  1]. . .  x[n]  WN k x[n  1]  WN k y[n  2]. . .  x[n]  WN k x[n  1]  WN 2k x[n  2]  WN k y[n  3] . . (expand and simplify ) n.  x[k ]WN(n  k ). k  .  WN n. n.  x[k ]WNk. (2-5). k  . 利用(2-5)式，則可算出一長度為 N 的訊號 x[n]的第 k 項頻率分量 yk [n]  WN N. N. . x[n]WNk m0.  1. N.  x[n]WNk. (2-6). m0. (2-6)式中，等號右側的式子與 DFT 的形式相同（見(2-1)式），但差別在於(2-6)式的加總上界為 N。這表示，使用這個形式的算法時，訊號必須要有 N+1 個點。但實際上訊號只有 N 個點，因此額外需要的那一點資料必須自行補上。然而，資料是不可以隨意更動的，隨意補上的點很可能因為與實際的訊號邊界數值相差過多而造成分析的誤差。因此，另一個較好的做法是透過尤拉公式(Euler’s equation)將(2-4)式改寫如下 H (z ). . . 1  WNk z 1. 1. 1  WN k z 1 1  WNk z 1 1  WNk z 1. 1  (WNk  WN k ) z 1  z  2.  1  (e. i. 2 k N. 1  WNk z 1 e. i. 2 k N ) z 1  z  2. 11.

(21) 1  WNk z 1  2k 1 1  2 cos( )z  z 2 N. (2-7). 並且利用連鎖規則(Chain rule)引入暫時的變數 S(z)來重新表示輸入與輸出的關係. H ( z) . 1  WNk z 1 S ( z) Y ( z)  X ( z ) S ( z ) 1  2 cos(2k ) z 1  z  2 N. (2-8). 由(2-8)式，可以知道. 2k sk [n]  x[n]  2 cos( ) x[n  1]  x[n  2] N yk [n] . (2-9). sk [n]  WNk sk [n  1]. 從(2-9-2)式可以知道，對於一長度為 N 的訊號 x[n]，要求得 yk [N ] 的話，需要算到 sk [N ] 為止。而由(2-9-1)式可知，欲求得 sk [N ] ，需要算到 x[N ] 。綜合以上兩式的結論與(6)式相比，確實不需要額外補點來計算。不過另外要注意的是，雖然當 n  0 時， x[n  1] 與 x[n  2] 是沒有定義的，但我們可以將其設為 0，使這個濾波器的響應為零態響應(Zero state response)。綜合上述結論與數學式，我們便能實作 Goertzel 演算法，以下為虛擬程式碼(Pseudocode)：. 12.

(22) 表 2-1. Goertzel 演算法的虛擬程式碼 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21. // --- definition --sig: signal series n: signal length fs: sampling frequency ft: target frequency k = (int)(0.5+n*ft/fs) omega = 2*pi*k/N cosine = cos(omega) sine = sin(omega) coeff = 2*cosine s0, s1, s2 = 0 // temporary space for saving previous data for i from 0 to n s0 = coeff*s1-s2+sig[i] s2 = s1 s1 = s0 real = (s1-s2*cosine) imag = s2*sine mag = sqrt(real*real+imag*imag). // result. 至於計算複雜度方面，在給定 K 項欲求得的 DFT 項，對於一長度為 N 的訊號的情況下，根據(2-9)式，Goertzel 演算法的複雜度為 OKNM  ，其中 M 表示(2-9-1)式每一次迭代所花費的運算量。而 FFT 則必須視使用的演算法而定，若以 2 基底-快速傅立葉轉換(radix-2 FFT)來說，則 FFT 的計算複雜度為 OMN log 2 N  。回到本節一開始提到「使用 FFT 來分辨電力線干擾頻率並不符合效益」這點敘述。在不考慮上段所述 M 項的前提下，比較推得的 Goertzel 演算法與 radix-2 FFT 的計算複雜度；若 Goertzel 演算法要比 FFT 更有效率，則須滿足以下條件. 13.

(23) KNM  MN log 2 N. (2-10).  K  log 2 N. 因此，根據(2-10)式，在理論上只要使用 Goertzel 演算法時欲求得的 DFT 項數小於訊號長度取 2 的對數，計算速度就會比直接使用 FFT 來的快。. 2-1-2 短時 Goertzel 演算法(Short-time Goertzel algorithm) 然而，在實際處理訊號時，直接使用 Goertzel 演算法來分辨電力線干擾的頻率會遇到兩個問題： 1.. 頻率漂移(Frequency drift)：電力線的頻率並非固定不變的，通常會隨著電廠供電與外界用電之. 間的供需不平衡而變化[19]。當電力的供給大於需求時，供電的頻率會變高；反之則變低。而漂移的範圍依各國的規定與發電機組效能而有所不同，通常在±0.5Hz 內。. 圖 2-2. 供電頻率的變化[20] 2.. 頻率解析度過高： Goertzel 演算法和 FFT 的頻率解析度皆會受到訊號長度的影響。如. (2-1)式所示，若訊號長度 N 越大，其中 e. 14. i. 2 nk N 所代表的頻率基底就會.

(24) 被分割成越多項，而頻率解析度就越高。由於存在頻率飄移的問題，而 Goertzel 演算法一次只能算出一項 DFT 項，若指定 50Hz 與 60Hz 作為 Goertzel 演算法要分析的頻率，則算出的數值將可能不是實際的電力線干擾的成份。如圖 2-3，在大尺度的觀察下，電力線的干擾頻率約在 60Hz；但是再更進一步檢視時，發現實際上電力線干擾頻率最強的分佈並不是準確的落在 60Hz 上。. (a). (b). 圖 2-3. 電力線干擾在頻譜上的分佈 (a) 檢視範圍：0~100Hz；(b) 檢視範圍：56~640Hz 對於此問題，有一個簡易的解法是：用 Goertzel 演算法在 50Hz 與 60Hz 附近多分析數個點，再以個別區域的最大值當作分析結果。雖然電力線的漂移範圍小，不過當訊號長度變長，頻率解析度就隨之增高，而在同樣的頻率範圍內需要算出的 DFT 項數就越多，使用 Goertzel 演算法花費的時間就越多。然而，這個問題的關鍵在於頻率解析度，原因如下： 15.

(25) (1) 欲求的 DFT 項數的多寡會影響計算時間的長短。在一固定的範圍內，DFT 各項間距越小，須求的 DFT 項數就越多，而 DFT 的各項間距即為 DFT 的頻率解析度。 (2) 對於現今大部分生理訊號擷取設備的取樣頻率來說，頻率漂移的範圍不大（根據美國心臟學會(American Heart Association, AHA)的建議，以監測為主要用途的心電訊號量測設備，取樣頻率應在 150Hz 以上[21]；因此，以±0.5Hz 的頻率漂移範圍來說，在取樣頻率 150Hz 時，造成的誤差約 0.3%）。因此只要能降低頻率解析度，便能降低 Goertzel 演算法耗費的時間。但是，直接對訊號降取樣(Downsampling)很容易受混疊現象(Aliasing)影響。所以，為了達到此目的，應以調整分析的訊號長度為手段。而本研究根據短時傅立葉轉換(Short-time Fourier Transform)的概念，提出短時 Goertzel 演算法(Short-time Goertzel algorithm)。此演算法的流程為： (1) 將訊號以每秒取樣的點數為單位分割，若最後的片段不足一單位，可以選擇補足至一單位或是捨棄。 (2) 用 Goertzel 演算法分析每單位長度的訊號片段。 (3) 將每一段的分析結果加總。在流程(1)中，因為是以取樣頻率(1000Hz)的值為單位作分割，所以每一段訊號頻譜的頻率解析度會固定為 1Hz。也因為此解析度完全涵蓋了頻率漂移的範圍(±0.5Hz)，因此分析的結果不會受頻率漂移影響，而且 Goertzel 演算法也只需分析 50Hz 與 60Hz 這兩個目標頻率即可。. 16.

(26) 表 2-2. 短時 Goertzel 演算法的虛擬程式碼 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18. // --- definition --sig: signal series n: signal length fs: sampling frequency ft: target frequency w: desired length of a signal segment goertzel(sig, fs, ft, w): a function implementing Goertzel algorithm rem = n%w // remaining dlen = n-rem val = 0 cnt = 0 for i from 0 to dlen val += goertzel(sig[i:i+width], fs, ft, w) cnt += 1 i += width val /= cnt. // result. 2-2 電極接觸雜訊處理電極接觸雜訊是指因電極與皮膚接觸不良而造成的雜訊。由於目前生醫訊號的擷取設備多採用接觸式電極，因此訊號的品質難免會受到電極與皮膚接面耦合穩定度的影響。要降低電極接觸雜訊的影響，在硬體方面，主要是以改良電極的材質、幾何形狀為主。過去用於量測心電訊號的電極多採用溼式電極，是於皮膚與電極之間加入導電膠的方式，降低膚電阻抗(Skin-electrode impedance)與增加接觸面積，以減少電極與皮膚接觸不穩而造成影響。但由於導電膠的導電性會隨著時間增加而流失[22]，無法長時間維持良好的訊號品質，因此近年十來逐漸有乾式電極的應用。乾式電極通常透過加工表面的產生微結構來增加與皮膚的接觸面積及 17.

(27) 穩定度，在使用上確實比濕式電極較方便，也除去了上述的困擾。但也因為沒有導電膠的關係，不僅膚電阻抗大為提高，接面的穩定度也容易受皮膚上毛髮的影響而降低[23]。為維持訊號品質的良好，首要仍為電極的改良。但是在長時間的監控需求下，要確保訊號不受到電極接觸雜訊的影響仍是個難題。因此在軟體方面，也需有對應的措施來處理電極接觸雜訊。. 2-2-1 電極接觸雜訊的特徵觀察未經濾波處理的心電圖，電極接觸雜訊以突波 (Spike) 或斷點 (Breakpoint)為起始，在其後通常會接著一段平坦線(Flat line)。從實際電極與皮膚接觸的情況來解釋，突波的出現代表該瞬間電極鬆脫，原因是電極導線受拉扯或者其他外力因素影響到電極與皮膚的耦合狀況；平坦線則代表電極與皮膚耦合未穩定，其持續時間與干擾的強弱和訊號擷取電路的穩定時間 (Settling time)有關[23]，且電壓數值通常為訊號擷取設備的訊號飽和值；而斷點的成因與突波類似，只是時間點通常發生於訊號擷取設備正處於要恢復正常交流耦合(AC coupling)時，或是電極完全脫離皮膚接觸的那一瞬間開始，所以從圖上可觀察到斷點處即為訊號振幅有極大落差的時間點（如圖 2-3 所示）。為方便描述，後文對於突波與斷點的情況將只以突波兩字來代表。而因為電極接觸雜訊的干擾模式不固定，而突波又代表極高頻的干擾，若使用一般的線性濾波器來處理，會造成嚴重的失真(Distortion)，影響後續的處理。因此，為了避免濾波造成的訊號失真，本研究採用的策略為：找出沒有受到電極接觸雜訊干擾的訊號片段，提供其兩端點位置給後續的特徵偵測演算法使用。. 18.

(28) 圖 2-4. 電極接觸雜訊的特徵. 2-2-2 前處理由於正常的心電訊號特徵在頻域上分佈各有不同，因此要達到上述的策略，必須改為利用電極接觸雜訊的突波特徵來反向操作。然而訊號可能同時受到多種不同的干擾影響，因此第一步必須先設法讓區隔電極接觸雜訊與其他的干擾。所以，訊號會先經過截止頻率為 0.8Hz 與 50Hz 的 4 階 Butterworth 帶通濾波器(Bandpass filter)處理，藉此濾除低頻的基準線漂移干擾與電力線干擾等高頻的干擾。經過帶通濾波器的處理後，電極接觸雜訊的特徵會稍作改變（如圖 2-5 所示），而我們後續將根據下述的兩種情況來做處理： . 情況 1：暫態峰. 暫態峰(Transient peak)為線性濾波器對於突波與斷點此類變化極快的訊號所產生的暫態響應。而暫態峰的持續時間與突波的高低和濾波器本身的響應速度有關。 . 情況 2：平坦線 19.

(29) 未經處理訊號中，平坦線的電壓數值通常為於訊號飽和的數值，可以視為一段單純的直流偏壓(DC bias)訊號。而因為經過低頻截止頻率為 0.8Hz 的帶通濾波器，直流訊號的直流成份會被濾除，因此平坦線的電壓值會被拉近至 0。注意，電壓值只是接近 0，並不能保證準確為 0。這是因為心電訊號中各個特徵波的起伏高低不同，而經過濾除直流的處理後，只能確定整段訊號的平均值為零，表示 0 電位以上與以下的所有訊號面積總合為 0。. 圖 2-5. 濾波後的電極接觸雜訊. 2-2-3 粗粒化(Coarse-graining approach) 我們知道經過濾波的電極接觸雜訊有暫態峰這個明顯的特徵可以利用，因為其峰值通常遠比正常心電訊號中 R 波的峰值還大，所以直覺上可以利用搜索極值的方式來找到電極接觸雜訊發生的時間點。但因為電極接觸雜訊發生的次數不固定且不一定連續，而且暫態峰屬於一個連續的訊號，所以直接作極值搜索並不恰當。然而，在看過許多受電極接觸雜訊干擾的資料後，本研究發現人眼可以 20.

(30) 很容易地從連續且規律的訊號中找出突然改變的部份。而這是因為人類視覺系統(Human visual system)具有低通濾波的特性[24]，使得映入眼中的邊界（即高頻變化）有被強化的效果，如馬賀帶(Mach bands)錯覺。所以，雖然心電訊號中的 QRS 波群(QRS complex)屬於較高頻的特徵，但通常其數量遠多於電極接觸雜訊的特徵，而且相較之下是有規律的訊號，因此人眼能輕易地區分其與電極接觸雜訊的差別。但是，電腦並沒有這種視覺功能。因此，為了模擬人眼的這種功能，本研究利用了粗粒化(Coarse graining)[25]的方法來實現。粗粒化可以視為一種降取樣(Down-Sampling)的方法，其概念為：將一個訊號序列以每 n 個資料點視為一個顆粒向量(Grain)，最後用各個顆粒的一個代表性數值（平均值、最大值、最小值…等）取代原本的 n 個資料點。以下見圖 2-6 來說，第一行為原始的訊號序列，若將顆粒向量的大小設為 2，表示將每兩個訊號點視為一個顆粒，而粗粒化的結果以 x2,1, x2,2…表示。（下標第一個數字為顆粒向量的大小，第二個數字為該顆粒向量的編號）. 圖 2-6. 粗粒化流程示意圖所以，經由將訊號粗粒化後，不僅使資料量減少，更因為連續的暫態峰 21.

(31) 被離散化(Discretized)，讓搜索操作可以被用上。. 2-2-4 切入點偵測(Cut-in point detection) 雖然透過粗粒化解決了暫態峰的問題，可以比較容易地找到電極接觸雜訊的發生位置，但這還不足以用來定位出正常訊號的片段。因為我們還要考慮從電極接觸雜訊發生後，到皮膚與電極耦合情況恢復正常，使訊號擷取設備正常運作的時間，也就是前面所述的情況 2：平坦線。暫態峰與平坦線最明顯的差異就在於訊號準位(Signal level)：暫態峰的數值通常都很大，而平坦線是接近 0。因此我們可以透過描繪出訊號的上包線(Upper envelope)，再配合適當的閥值(Threshold value)來區隔暫態峰與平坦線。而要找到訊號的上包線，本研究使用的粗粒化方法是以每個顆粒中數個資料點的最大值來表示。以下是這兩種情況的處理流程： . 情況 1：暫態峰. 步驟 1. 由於暫態峰的變化速度相較於原始的突波或斷點來的慢，表示頻率分佈自高頻降低。為了與心電訊號中的 R 波有更明顯的區隔，首先會將訊號透過 Pan-Tompkins[12]的前處理流程：帶通濾波、微分、平方、移動平均，這四個方法來強化差異。步驟 2. 將訊號粗粒化找到上包線( UE A )，其中顆粒寬度( wA )為. wA  0.08 * f s (= 80 ms). (2-11). 之所以顆粒寬度為如此，是因為正常的 QRS 波群最短寬度為 80 毫秒。因此粗粒化後的訊號不會因為顆粒寬度過大而讓 QRS 波群與其他特徵波群難以分辨。步驟 3. 為了定義閥值，先暫時忽略其他特徵波，我們假設每一秒至少能夠擷取到一個 QRS 波群（取正常人心跳為 60 bpm 的情況）。而因為. 22.

(32) 在步驟 2 中，顆粒寬度( wA )為 80 毫秒，所以若要讓一個 QRS 波群能夠代表寬度為一秒的顆粒，我們必須乘上一個補償係數( rcomp )： rcomp . fs wA. (2-12). 而閥值( TA )則定義為： TA  rcomp  Mean(UE A ). (2-13). 步驟 4. 找出所有上包線( UE A )上超過閥值( TA )的點。而這些點就是電極接觸雜訊出現的位置。（見下圖 2-7）. 圖 2-7. 切入點偵測（暫態峰） . 情況 2：平坦線. 步驟 1. 在考慮處理的訊號只有受到平坦線的影響下，由於平坦線的訊號準位接近於 0，與正常的 QRS 波群明顯可分辨，所以直接將訊號粗粒化。而構成此上包線( UE B )的顆粒寬度( wB )為. wB  f s (= 1 sec) 23. (2-14).

(33) 而此處顆粒寬度( wB )為 1 秒，是因為要區別平坦線與 QRS 波群只需要從訊號準位作區分。因此在取正常人每分鐘心跳次數是 60 下的情況下，一個顆粒寬度只需要為 1 秒，便能確保每次心跳至少能被一個顆粒捕捉。由於經過粗粒化後，訊號的上包線( UE B )只會有平坦線與代表 QRS 波群的顆粒存在。如下圖 2-8 所示，(a)圖為訊號未經粗粒化時，即便 R 波波峰振幅較高，但是由於訊號已經經過濾波，直流成份被濾除，因此在振幅分佈的直方圖中難以區別平坦線與 QRS 波群；而(b)圖為訊號經粗粒化的結果，可以明顯看到 R 波的特徵被凸顯出來，剩下接近 0 的分佈則為平坦線。. (a). (b). 圖 2-8. 粗粒化前後的心電訊號振幅分佈直方圖 24.

(34) 步驟 2. 由於粗粒化後，QRS 波群的特徵與平坦線已經可以很容易地被區別，因此直接定義閥值( TB )為. TB  Mean(UE B ). (2-15). 步驟 3. 因為已知平坦線只會出現在暫態峰之後，而前述的情況 1 已經可以用來找出電極接觸雜訊的位置，所以可以確定平坦線只會出現在每一段正常訊號的前方。因此在此步驟中只要找出上包線( UE B )上第一個超過閥值( TB )的點，即可確定為偵測的切入點。. 2-2-5 正常訊號片段選擇有了上一小節對平坦線與暫態峰的處理方式後，配合平坦線只會出現在暫態峰之後的特性，要找到正常訊號片段的處理流程如下. 圖 2-9. 正常訊號片段選擇的處理流程 25.

(35) 2-3 濾波由於後續為心電訊號特徵波的偵測，所以在這個階段的目的是設法突顯特徵波並減少其他干擾對訊號的影響。然而，根據圖 2-10 的特徵波頻率分佈圖，QRS 波群、P 波與 T 波在頻域上分佈頗為分散，因此在本階段的處理是以與各個特徵波頻段較無重疊的干擾為主。而後要偵測各種特徵波時，再根據其頻率分佈來做帶通濾波才會是比較好的做法。而目前前處理的部分已經討論完電力線干擾與電極接觸雜訊這兩項。電力線干擾較重要的問題在於干擾頻率的偵測，而我們可以透過前面所提的短時 Goertzel 演算法來判斷，再根據其頻率使用 3 階的 Notch filter 來降低電力線干擾的影響。至於電極接觸雜訊，由於本研究是透過其特徵來找到沒有被干擾的訊號片段，因此後續沒有對其做額外的濾波處理。. 圖 2-10. 心電訊號特徵波的頻率分佈[26]. 26.

(36) 第三章演算法設計－特徵偵測心電圖特徵描繪演算法的目標是找出 QRS 波群(QRS complex)、P 波與 T 波。由於這些特徵波分別代表心臟的電生理活動資訊（如圖 1-1 所示），因此藉由判斷特徵波的異常與否便能快速檢測心臟的健康狀況。 P 波為心臟電生理活動的起始，表示心房去極化(Depolarization)，也就是組織收縮。QRS 波群代表心室的去極化，此時左右兩個心室分別要將血液輸出至肺部與全身。而因心室的纖維較粗、收縮力道較大，因此在心電圖上是最明顯的特徵。T 波代表心室的再極化(Repolarization)，即心室的舒張，其波形的振幅通常較 P 波來的大。其中 QRS 波群中的 R 波持續時間短，振幅通常是所有特徵波中最大的。因此許多心電圖的特徵描繪演算法都會以先找到 R 波為基礎，再從各個 RR 區間內找其他的特徵波。而本研究所提出的特徵描繪演算法也是如此，且分析的訊號是屬於第二導程(Lead II)，但分析的訊號片段是由第二章所提到的電極接觸雜訊處理演算法所決定。. 3-1 QRS 波群偵測本研究在 R 波偵測的部分，是以 Jiapu Pan 與 Willis J. Tompkins 於 1985 年所提出的演算法[12]為基礎（以下稱其為 Pan-Tompkins 演算法），並校正 R 波波峰的定位與閥值過濾，且透過 k-means 分群演算法來區別異常的 R 波或偵測錯誤的點。而 Q 波與 S 波則是在找到 R 波波峰後，從其左右兩側向外搜尋最低點。QRS 波群完整的偵測流程如下圖 3-1 所示. 27.

(37) 圖 3-1. QRS 波群偵測的流程圖. 3-1-1 R 波偵測 Pan-Tompkins 演算法的架構可分為三個階段：線性濾波(Linear filtering)、非線性轉換(Nonlinear transformation)、決策演算法(Decision rule algorithm)。在經過線性濾波後，可以降低心電訊號受到的干擾與相對的凸顯出 QRS 波群的特徵；而非線性轉換的部分則是為了定出 QRS 波群的位置；最後再透過決策演算法來尋找正確的 R 波波峰。而由於決策的部分採用迭代的方式處理，因此更容易實現為即時(Real time)處理的演算法而 Pan-Tompkins 演算法的處理流程圖如下圖 3-2 所示. 28.

(38) 圖 3-2. Pan-Tompkins 演算法的流程圖接下來本論文將更詳細的說明 Pan-Tompkins 的處理流程與本論文在實作上修改的部分。然而 Pan-Tompkins 的論文並沒有提供原始碼，為方便比對，本論文是以 Hooman Sedghamiz 所寫的 Matlab 程式碼[27]為基礎，改由 Python 語言實作出自己的版本。而圖 3-2 中的流程從 Filtering 至 Fiducial mark 可以與圖 3-3 的處理結果對照。. 29.

(39) (a). (b). (c). (d). 圖 3-3. Pan-Tompkins 演算法各步驟的處理結果 (a) 原始訊號；(b) 經濾波處理；(c) 經微分處理；(d) 經平方處理. 30.

(40) (e). (f). (g). (h). 圖 3-3（續）. Pan-Tompkins 演算法各步驟的處理結果 (e) 經移動平均處理；(f) 對位標記；(g) 決策流程；(h) 偵測結果 . 線性濾波：由於 QRS 波群在頻域上的分佈以 5~15Hz 最明顯，因此可以利用帶通濾 31.

(41) 波器使此頻帶的訊號凸顯出來。而原論文為了利用微處理機實現即時處理的演算法，需要將浮點數運算(Floating-point arithmetic)轉為整數運算(Integer arithmetic)，因此採用了透過在 z 平面上作極零對消(Pole-zero cancellation)設計而成的遞迴式濾波器。但由於原論文使用 200Hz 的取樣頻率，無法直接設計出通帶為 5~15Hz 的帶通濾波器，因此分別設計出低通濾波器與高通濾波器，再串接後使得通帶達到約 5~11Hz。但是經過濾波器的處理後，訊號難免會受到影響產生相位延遲(Phase delay)。對於即時處理的系統來說，這是無法避免的問題。然而本研究的演算法是用於離線(Offline)處理的，因此可以透過正反向濾波來達到零相位延遲(Zero-phase)的結果。（見圖 3-3(b)） . 非線性轉換：. 1. 微分(Derivative) 經過濾波後，會再透過微分的操作算出訊號點與點之間斜率。由於 QRS 波群的頻率分佈是心電訊號幾個特徵中較高的，所以微分後的數值會較高；而其他不是位於 5~15Hz 的訊號也已被上一階段的濾波器衰減，因此現在訊號經過微分後會以 QRS 波群的特徵最為明顯。（見圖 3-3(c)）而這裡的微分方式是五點差分(Five-point derivative)，其轉移函式為 1 H ( z )  ( T )(  z  2  2 z 1  2 z  z 2 ) 8. (3-1). 由上方(3-1)式轉為差分方程式為 1 y(nT )  ( T )[  x(nT  2T )  2 x(nT  T )  2 x(nT  T )  x(nT  2T )] 8. (3-2). 如此，由上方(3-2)式便能作為一個差分濾波器，其頻率響應如下圖所示. 32.

(42) 圖 3-4. 五點差分濾波器的頻率響應 2. 平方(Squaring) 平方的作用不僅是將所有負數轉為正數，也使得訊號被非線性的放大。由於訊號經過微分後，高頻的成份會被放大，使得我們能夠從斜率變化去分辨出高頻與低頻成份的差異。而經過平方後，則是能直接強化振幅方面的差異，使得 QRS 波群與其他特徵更容易分辨。（見圖 3-3(d)） 3. 移動平均(Moving average) 經過上述的微分與平方的處理後，QRS 波群的位置已經凸顯出來。然而，因為後續需要作極值的搜索，為減少區域極值造成搜索的困難，此時還須透過移動平均濾波器將經上述幾項步驟處理過的訊號作平滑化。但是移動平均濾波器的窗口大小必須適當，理想的結果是在訊號經過移動平均濾波器後只有 QRS 波群的位置會被圈選出來。如果窗口太大，有可能使得兩個 QRS 波群被圈選在同一個區塊，使得後續偵測的流程無法辨識；若窗口太小，則無法良好地做到平滑化的效果，使得做極值搜索時可能錯誤地找到區域極值。在原論文中，作者在取樣頻率為 200Hz 的情況下，使用的窗口大小為 150ms。而本論文根據 Borys Surawicz 等人在美國心臟協會 33.

(43) (American Heart Association, AHA)上發表的文章[28]判斷，健康成人的平均 QRS 波群的最大寬度可達 110ms 以上，因此也將移動平均濾波器的窗口大小設為 150ms。處理結果請見圖 3-3(e)。 . 決策演算法：. 1. 對位標記(Fiducial mark) 觀察圖 3-3(b)與(e)，可以發現訊號經過移動平均濾波器處理後產生的輪廓確實能夠將 QRS 波群圈選出來。然而，我們還需要從這些輪廓裡找到實際的 R 波波峰。由於訊號在經過平方的處理後，所有數值皆大於或等於 0，因此使用移動平均濾波器時，只要其窗口內涵蓋 QRS 波群的範圍增加，移動平均濾波器的輸出便有增加的趨勢；反之，則輸出結果會有下降的趨勢；而當窗口完全涵蓋 QRS 波群但沒有觸及其他非 QRS 波群的訊號，則移動平均的輸出會保持固定。這樣的現象如下圖 3-5 所示. 圖 3-5. QRS 波群與移動平均濾波器窗口的位置關係[12] 因此，我們可以確定 QRS 波群會位於移動平均濾波器輸出結果的上升邊緣至平坦區之間。所以，要找到 R 波波峰，可以先尋找經過移動平均濾波器處理後的結果的最大值，再進一步做定位即可。. 34.

(44) 要找到各個 QRS 波群外框（即圖 3-3(e)中各個凸起的波形）的最大值，並不能直接搜索整段訊號。因為這項步驟要處理的問題為：在一個具有多個數值接近但不盡相同的最大值的數列中，找出所有最大值的位置。在 PanTompkins 的論文中沒有詳細敘述此步驟的實作；Hooman Sedghamiz 則是利用 Matlab 中的 findpeaks 函數來逐步搜索；而本論文的實作架構與 Hooman Sedghamiz 的版本類似，僅有修部份處理流程。 findpeaks 函數的原理是：逐步比對資料點，當資料點數值比前一個小時，表示找到極值。但由於訊號有許多起伏，直接利用這種概念處理的話，會找到非常多不是我們要找的結果。因此還需配合其他手段來處理，在 findpeaks 中有兩個主要的處理方式： 1. 僅搜尋數值大於給定閥值以上的資料 2. 搜索到一個極值後，此點位置至少和上一個被搜索到的極值有一定的距離，否則不計而本論文採用的方式與 Hooman Sedghamiz 相同，由於每個人心室收縮的力道不盡相同，使得 R 波振幅也不一樣，所以僅使用上述的第二個方法來處理。但考慮到心臟有異常活動的心率(Heart rate)範圍，目前以房室結迴旋性頻脈 (AV nodal re-entry tachycardia, AVNRT) 的情況最高可能達到 280bpm(Beats per minutes)的速率，因此設定以 300bpm 為可能偵測到心率的最高上限，也就是每個被搜索到的最大值之間距離至少要大於 200ms。（如下式所示）. f  300 (bpm). T  1/ f.  300 / 60(beat / sec)  5(beat / sec)  1/ 5(sec/ beat)  200 (ms /beat). (3-3). 如此，便能大幅減少搜尋到非 QRS 波群外框的最大值。（見圖 3-3(f)）但後續仍要過濾此步驟找到的結果，才能找到真正對應 R 波波峰的標記。 2. 閥值過濾演算法 35.

(45) 此處理流程是將上一步驟所找到的標記一個個地依序處理，每一次的處理迴圈內容如下方圖 3-6 所示。由於 Pan-Tompkins 原文僅有流程介紹，所以實作方面會因人而異。而本論文在這部分的演算法是前述 Hooman Sedghamiz 所著的程式碼改良的版本，圖 3-6 的流程即分析自其程式。接下來會說明各個步驟與本論文修改的部份，其中符號與對應的意義如下表所示表 3-1. Pan-Tompkins 決策演算法中各符號之意義符號. 意義. SIG_BPF 經過帶通濾波器處理的訊號(圖 3-3(b)) SIG_MAI 經過帶通濾波、微分、平方與移動平均處理的訊號(圖 3-3(e)) PKF. SIG_BPF 中可能為 R 波的標記，位置與振幅表示為(tf , yf). PKF’. 同 PKF，但是在往回搜索階段才需找的值. SPKF. 當 PKF 被標記為 R 波時的訊號準位估計值. NPKF. 當 PKF 被標記為雜訊時的訊號準位估計值. PKI. SIG_MAI 中可能為 R 波的標記，位置與振幅表示為(ti , yi). PKI’. 同 PKI，但是在往回搜索(Search-back)階段才需找的值. SPKI. 當 PKI 被標記為 R 波時的訊號準位估計值. NPKI. 當 PKI 被標記為雜訊時的訊號準位估計值. RRn. 距離當前處理迴圈的最近的 RR 區間（由先前的 PKF 算得）. RRnew. 目前處理迴圈的 RR 區間. RRavg1. 最近的 8 個 RR 區間之平均. RRavg2. 同上，但作為暫存空間，用於往回搜索階段. PKI_ARY 儲存 PKI 的陣列，可被用來算 RRavg1 PKIprev. PKI_ARY 中最新的點. THR_I1. 用於判斷 PKI 是否為 R 波的第一閥值. THR_I2. 用於判斷 PKI 是否為 R 波的第二閥值. THR_F1. 用於判斷 PKF 是否為 R 波的第一閥值. THR_F2. 用於判斷 PKF 是否為 R 波的第二閥值 36.

(46) 圖 3-6. Pan-Tompkins 演算法中閥值過濾的流程圖 37.

(47) 在 Pan-Tompkins 原論文中，此閥值過濾演算法有三個階段，分為：學習階段 1、學習階段 2 與偵測階段。學習階段 1 會利用到訊號前兩秒的片段作論文沒有詳細說明初始化的數值設定，因此本論文使用與 Hooman Sedghamiz 版本相同的初始化設定，細節如下： THR_I1  max( SIG _ MAI 2 sec ) / 3 THR_I2  mean (SIG _ MAI 2 sec ) / 2 (3-4) THR_F1  max( SIG _ BPF2 sec ) / 3 THR_F2  mean (SIG _ BPF2 sec ) / 2 其中 SIG_MAI2sec 為 SIG_MAI 的前兩秒片段，而 SIG_BPF2sec 為 SIG_BPF 前兩秒片段。但此設定在遇到前兩秒訊號的起伏不大時，可能使得閥值設定的結果不合理，詳細情形將在第四章說明。然而根據圖 3-1 所示，由於本研究在偵測 R 波前，會先由電極接觸雜訊處理流程決定要分析的訊號片段。而此問題的情況與其處理流程中的描述平坦線類似，因此本研究的實作可以避免此問題的影響。而學習階段 2 則是初始化用來判斷 RR 區間長度的閥值 RRavg1 與 RRavg2。但因為未行偵測前並無法確定哪些標記為實際的 R 波，因此在初期會先配合振幅高度的閥值過濾來逐步地記錄可能為 R 波的標記，直到數量大於 8 個後才會用來計算 RRavg1 與 RRavg2。接下來的偵測流程會利用 SIG_BPF 與 SIG_MAI 兩種訊號來判斷輸入點是否為 R 波。因為移動平均濾波器的遮罩寬度為 150ms，雖然能適當的把各個 QRS 波群獨立地框選出來，但如果訊號受到較高頻的干擾，可能使得濾波器的行為無法如圖 3-5 那樣正常地表現。以下配合圖 3-6 說明： (1) 新的待評估資料輸入：由於可能為 QRS 波群的訊號位置已如圖 33(e)的情況，因此會先以 PKI 為基礎來作搜索。 38.

(48) (2) 搜尋 PKF：根據上一小節關於對位標記的敘述，由於 PKI 是從 SIG_MAI 中搜索而得的，實際的 R 波位置會落在 QRS 波群外框的上升邊緣(見圖 3-5)。而因為移動平均濾波器始訊號產生 150ms 的延遲，所以搜尋 PKF 的範圍訂為[ti -150ms, ti]。 (3) 計算 RRn 與 RRavg1：由儲存於 RR_ARY 裡最接近目前時間點的 8 個 PKI 算出 RRn 與 RRavg1，用以判斷心率的變化。而 RRn 與 RRavg1 的關係為. RRavg1  ( RRn 7  RRn 6  ...  RRn ) / 8. (3-5). (4) 判斷 RRn 是否在範圍內：正常心率變化的範圍為. 0.92 RRavg 2  RRn  1.16 RRavg 2. (3-6). 演算法會藉由步驟(3)、(4)定義接下來要判斷標記的範圍，直到找到下一個被判斷是 R 波的 PKI 為止。 (5) 更新 RRavg2：若 RRn 落在範圍內，表示上一個 RR 區間內沒有遺漏該被偵測到的 R 波。 (6) 判斷 yi 是否大於 THR_I1：若是，則此 PKI 可能是 R 波。接下來必須確認 PKF 是否也有符合條件。 (7) RRnew 是否大於 360ms：若否，表示目前偵測的點與前一個 R 波距離太短，必須確認是否誤把 T 波當作 R 波。. RRnew  ti  ti prev. (3-7). 其中 xi' 為 PKI_ARY 中最新點 PKIprev 的位置。 (8) 判斷 yf 是否大於 THR_F1：若是，表示 PKF 為實際的 R 波；否則跳過步驟(9)與(10)。 (9) 將 PKF 標記為 R 波：經過步驟(6)與(8)的檢驗，此 PKF 可以被視為 R 波，因此存入輸出結果。但實際上可能並非 R 波，此問題將於下 39.

(49) 一節 k-means 分群演算法討論。 (10) 更新 SPKF：SPKF 是用來校正 THR_F1 與 THR_F2 的訊號準位，當 PKF 被視為 R 波時，SPKF 就會被更新。其關係式為. SPKF  0.125 y f  0.875 SPKF. (3-8). (11) 將 PKI 存入 PKI_ARY 並更新 SPKI：由於 PKI 在 SIG_MAI 上被視為 R 波，所以更新 SPKI 並存入 PKI_ARY 用於下一迴圈時更新 RRn 與 RRavg。而 SPKI 的更新式為. SPKI  0.125 yi  0.875 SPKI. (3-9). (12) 檢查 yi 是否介於 THR_I1 與 THR_I2 間：由於在步驟(6)判斷不符合條件，因此 PKI 不會是 R 波。在 Hooman Sedghamiz 的版本中，會利用此步驟決定是否將 PKI 標記為雜訊。 (13) 將 PKI 標記為雜訊：在 Hooman Sedghamiz 的版本中，並沒有再額外處理被標記為雜訊的 PKI。而本研究在比對 Pan-Tompkins 原論文後，也確定沒有對非 R 波的 PKI 做其他處理，因此本研究的實作將此步驟與步驟(12)刪去。 (14) 更新 NPKI 與 NPKF：經過先前的判斷，如步驟(6)或(15)，PKI 被判定不是 R 波，因此視為雜訊，所以更新此兩數值。. NPKI  0.125 yi  0.875 NPKI NPKF  0.125 y f  0.875 NPKF. (3-10). (15) 檢測是否為 T 波：當步驟(7)檢測不合條件時，接著由此步驟利用該 PKI 與 PKI_ARY 中最新點 PKIPrev 的斜率差異判斷其是否為 T 波。原論文是使用最大斜率來比較，但為了避免雜訊影響而誤判，本論文採用與 Hooman Sedghamiz 相同的方法，改以平均斜率來判斷。以下定義：. 40.

(50) Slope1: 在 SIG_MAI 上，PKI 前 75ms 各點之間斜率的平均值 Slope2: 在 SIG_MAI 上，PKIprev 前 75ms 各點之間斜率的平均值若 Slope1<0.5Slope2，則將此 PKI 視為 T 波。 (16) 降低 THR_I1 與 THR_F1：由於在步驟(4)發現 RRn 不在心率變化的正常範圍內，因此降低閥值 THR_I1 與 THR_I2，讓後續的偵測流程確認是否有遺漏偵測。. NPKI  0.125 yi  0.875 NPKI NPKF  0.125 y f  0.875 NPKF. (3-11). (17) 檢查 RRnew 是否大於 1.66 倍的 RRavg2：若是，表示可能遺漏偵測。 (18) 往回搜索，找出 PKI’：為了確認是否真的有遺漏偵測，所以在當前 PKI 與前一個點 PKIprev 之間搜尋最高點 PKI’，並由後續的步驟確認 PKI’是否為 R 波。在 Hooman Sedghamiz 的版本中，根據前述對位標記處理部分的(3-3)式，由於目前已知最快的心率不會超過 300bpm，因此搜索範圍由頭尾向內縮短 400ms，即[ti + 200ms, ti prev - 200ms]，藉此減少電腦處理步驟。然而，這步額外的處理可能使得程式出錯。如果該筆資料的心率偏快，便有可能使得步驟(17)中的 1.66RRavg2 小於 400ms，進而讓向內縮減 400ms 的指令出錯。. 1.66 RRavg 2  400 ms  0.4s  RRavg 2  240 .96 ms  0.241s. (3-12).  HRavg 2  60 / RRavg 2  60 / 0.241  248 .96 (bpm) 根據上方的推算，在任一時間點，其前 8 個點的平均心率必須小於 248.96bpm，否則就會造成錯誤。雖然 248.96bpm 是一個非常快的心率，正常人不太可能會達到這種情況，但因為我們無法保證接受量測的人心臟都正常，所以在處理流程上仍須將這種問題考慮進去。. 41.

(51) 而本研究的解決方式是在步驟(17)加上判斷[ti , ti prev]區間是否大於 400ms，若否，則放棄後續的處理，從步驟(1)開始處理下一個 PKI。 (19) 檢查 yi’是否大於 THR_I2：觸發往回搜尋的原因可能是先前的閥值過高，導致新的 R 波無法被偵測到。因此，在往回搜尋的流程中，會以第二閥值作為判斷條件。若 yi’大於 THR_I2，表示其可能為 R 波；若否，則判斷是可能是受其他原因影響造成 RR 區間突然變長，因此便沒有做其他處理，回歸到步驟(6)繼續檢查 PKI。 (20) 搜尋 PKF’：原因同步驟(2)，因為要從 SIG_MAI 轉到 SIG_BPF 上搜尋，因此搜尋範圍訂為[ti’ -150ms, ti’]。 (21) 檢查 yf’是否大於 THR_F2：若是，表示 PKF’為先前偵測遺漏的 R 波；若否則忽略該 PKF’。 (22) 將 PKF’標記為 R 波：經過步驟(19)與(21)的檢驗，可以確定 PKF’為 R 波，因此將其加入輸出結果。 (23) 更新 SPKF：與步驟(10)的作用相同，但注意此步驟更新 SPKF 的係數與其不同。. SPKF  0.25 y f  0.75SPKI. (3-13). (24) 將 PKI’存入 PKI_ARY，並更新 SPKI：與步驟(11)的作用相同，但注意此步驟更新 SPKI 的係數與其不同。. SPKI  0.25 yi  0.75SPKI. (3-14). (25) 更新閥值 THR_I1、THR_I2、THR_F1 與 THR_F2：經過前述步驟處理，SPKI、SPKF、NPKI 與 NPKF 可能有被更新，因此重新校正閥值。公式如下. THR_I1  NPKI  0.25 | SPKI  NPKI |. 42. (3-15).

(52) THR_I2  0.5 THR_I1 THR_F1  NPKF  0.25 | SPKF  NPKF | THR_F2  0.5 THR_F1 以上即為本研究所採用的 Pan-Tompkins 演算法。但因為此演算法是根據電極接觸雜訊處理流程給定的訊號片段做分析，因此當全部訊號分析完畢後還須注意到每一段分析結果的銜接。. 3-1-2 誤判 R 波的處理 R 波偵測的演算法容易受到其他因心臟活動異常的特徵波影響，但 R 波通常是用來偵測其他特徵波的重要基礎，因此提高演算法偵測的準確率是非常重要的。在上一節有提到，某些標記雖然已經通過雙重的檢測，但仍有可能不是真正的 R 波。這個問題與 Pan-Tompkins 演算法中閥值過濾的設計有關，其傾向於保留更大的數值當作 R 波，詳細的內容會在第四章探討。而因為 Pan-Tompkins 演算法的這個特性，使其無法分辨出正常 R 波與部分異常的高振幅特徵，其中 VPC(Ventricular premature contraction, 心室早期收縮)是常見造成錯誤判斷的原因之一。根據早期的統計資料[29]，VPC 即使在健康的人身上也偶爾會見到，而全年齡層都有相關的紀錄，且隨著年紀增加有越高的發生機率。因此，區分 VPC 與正常 R 波是一項重要的課題。要解決這項問題，必須先了解 VPC 的特徵。VPC 發生的時間點在 T 波之後，而 T 波代表一個正常的心跳周期要進到結尾(P→Q→R→S→T)，此時心室纖維逐漸舒張，但普金氏纖維(Purkinje fiber)卻異常地重新刺激心室纖維收縮，使原本應由竇房結(Sinoatrial node, SA node)引發的自主性心跳跳脫，即為所謂的室性逸搏 (Ventricular escape)。因此 VPC 在心電圖上的特徵與 R 波類似，為 T 波之後出現一個持續時間短但比 R 波稍長、振幅比 R 波大的波峰，跟著其後的第一個 QRS 波群前不一定有 P 波。（如圖 3-7） 43.

(53) 圖 3-7. 含有 VPC 的心電訊號由於振幅是 R 波與 VPC 最明顯的差異，因此本研究利用 k-means 分群演算法(k-means clustering algorithm)[30]來區別這兩者。 . k-means 分群演算法 k-means 分群演算法中所謂的 k 即代表要將資料分群的數目，其目的是. 將Ｉ維度中的 N 個資料點分為 K 個群。而其分群的概念是藉由反覆的尋找新的群中心並重新分群，直到各個群中心與自己群內的資料點的歐氏距離 (Euclidean distance)的平方（即平方誤差）達到最小，便結束演算法。因此，對於一個具有 N 點的資料群 x，k-means 演算法流程如下[31]： (1) 初始化：隨機設定 K 個中心點的位置，第 k 群的中心以 m(k ) 表示。 (2) 分類：將所有資料點分配到離自己最近的中心，以 kˆ ( n) 表示各點 x(n) 所屬的群。 x(n) 為 N 個資料點內的一點，而 kˆ ( n) 的定義為 kˆ ( n)  arg min {d (m( k ) , x( n) )}. (3-16). k. 其中 d (m( k ) , x(n) ) 為中心 m (k ) 與資料點 x(n) 的平方誤差，而 I. d (m( k ) , x( n) )   mi( k )  xi( n). 2. (3-17). i 1. 也就是說只要判斷某點 x(n) 與某群中心 m (k ) 的平方誤差是否為最小，. 44.

(54) 即可以用一個標籤(Label) rk(n) 以 0 或 1 表示該點是否屬於該群，即. 1, if kˆ ( n)  k rk( n)   0, if kˆ ( n)  k. (3-18). (3) 更新中心點：各群根據自己群內的資料點更新，以第 k 群的新中心點 m(k ) 來說. m(k ) .  rk(n) xi(n) n. i. (3-19). R(k ). 其中 R ( k )   rk( n) n. (3-20). (4) 重複步驟(2)與(3)，直到資料點的分群不再變動。而因為重複步驟(2)與(3)時，每次都是挑選最小值，因此可以確定結果會收斂，但並不能保證是全域的最佳解。由於本研究使用 k-means 演算法處理問題的維度只有 1，而且明確知道分群的結果是兩個振幅高低不同的群，因此初始化時不使用隨機產生的中心，而是直接根據 R 波與 VPC 的振幅設定兩個中心（以 CR 與 CVPC 表示）。而因為這樣的設定，初始的中心點與待分群的資料已有較小的平方誤差，因此運算可以更快地收斂。然而，並不是每筆資料都會有 VPC。但根據上述的初始化設定，在這種情況下，第一次進行步驟(2)時幾乎已經把所有 R 波都歸類到 CR 這群，因此不會陷入硬是把 R 波資料點拆成兩群的情況。. 3-1-3 Q、S 波偵測處理完 R 波後，Q、S 波的偵測就變得相對容易。由於 Q、S 波分別為 R 波的起始(Onset)與結束(Offset)，因此只要從 R 波前後搜尋最靠近 R 波的 45.

(55) 低點便能找到 Q、S 波。而考量到一般成人的 QRS 波群寬度，其最大可達 110ms 以上，因此本研究將 Q、S 波的搜尋範圍分別訂為. Q : x R  60 ms ~ x R S : x R ~ xR  60 ms. (3-21). 其中 xR 為 R 波的位置。. 3-2 P、T 波偵測本研究在 P、T 波偵測的部分，是改良自 Mohamed Elgendi 於 2013 年提出的演算法[32]（以下稱其為 Elgendi 演算法）。Elgendi 演算法除了可以同時分析出心電訊號中的 P、T 波，還能夠分辨出心房早期收縮(Premature Atrial complex, PAC/APC)或部分類似心室性頻脈(Ventricular tachycardia, VT)造成 T 波與 P 波看似疊合的波形。但是因為其演算法的設計，使得部分無明顯 P、 T 波的資料也會被判斷成有正常 P、T 波。因此，本研究除了加上流程來處理這個問題，也有針對其中部分流程做修正。. 3-2-1 Elgendi 演算法 Elgendi 演算法主要分為三個階段：前處理(Preprocessing)、特徵萃取 (Feature extraction)、分類(Classification)。演算法的流程如圖 3-7 所示，在前處理階段主要是為了凸顯 P、T 波的特徵，因此本研究會根據這兩種波的頻率分佈進行濾波，並且使用前一節提到的 QRS 偵測演算法的結果來移除 QRS 波群，以利後續的處理；而特徵萃取的部分是利用移動平均線交叉 (Moving average crossover)來找到感興趣的區間(Blocks of interest, BOI)；最後再利用一般 P、T 波在 RR 區間內的位置來找到實際的 P、T 波。. 46.

(56) 圖 3-8. Elgendi 演算法的流程圖[32] 而因為原論文是針對 MIT-BIH 的心律不整(Arrhythmia)資料庫處理，因此在設計上，本研究與其會有部分差異。所以接下來會根據上圖 3-8 來說明本研究實作此演算法的細節還有與原論文的差別。而下圖 3-9 是對應各個步驟的處理結果. 47.

(57) (a). (b). (c). (d). 圖 3-9. Elgendi 演算法的處理結果 (a) 原始訊號；(b) 經帶通濾波；(c)移除 QRS 波群；(d) 移動平均交叉. 48.

(58) (e). (f). 圖 3-9（續）. Elgendi 演算法的處理結果 (e) 定義感興趣的區間；(f) 偵測結果. . 前處理. 1.. 帶通濾波(Bandpass filtering) 由於 Elgendi 演算法原始設計也並非用於即時處理的系統上，因此本論. 文與原論文都是使用零相位延遲的帶通濾波器，並根據圖 2-9 與原論文，將帶通濾波器的頻寬設為 0.5~10Hz，而階數為 3。（見圖 3-9(b)） 2.. 移除 QRS 波群(QRS removal) 由於 QRS 波群通常是 ECG 中最明顯的特徵，因此將其移除可以使得 P、. T 波相對明顯。由於原論文使用的資料庫中，R 波的位置已有提供，所以沒有再使用其他的 QRS 偵測演算法。而移除 QRS 波群則是以資料庫的 R 波位置，分別再往前 83ms 與往後 166ms 的範圍內尋找第一個小於 0 的點，再將介於這兩點之內的數值全部設為 0。. 49.

(59) 而因為本研究使用的資料庫僅有原始資料，因此是使用前一章節所提的演算法來找出 QRS 波群的位置，再將 Q、S 點內的各點以線性內插(Linear interpolation)的方法補起。即. y[i] . y s  yq ts  tq. (t s  t[i]). (3-22). 其中 ts 與 ys 分別為 S 波的位置與高度，tq 與 yq 分別為 Q 波的位置與高度， t[i]與 y[i]則是介於 Q 波與 S 波內某點的位置與高度，處理結果為圖 3-8(c)。 . 特徵萃取. 1.. 移動平均交叉移動平均交叉這項技巧，是透過兩個不同遮罩大小的移動平均濾波器處. 理訊號後，再將兩個訊號相疊比對，找出交叉點。由於 P 波與 T 波兩者相較於 R 波都是較平緩的特徵波，因此使用移動平均交叉的技巧時，必須注意遮罩的大小。根據原論文，對於正常成人而言，在心率 60bpm 的情況下，平均 P 波的寬度為 90~130ms。而 P 波又比 T 波的寬度小，因此在此會以 P 波的資訊為基準。兩個移動平均的結果為. MApeak [n] . W 1 W 1 1 ( y[n  1 ]  ...  y[n  1 ]) W1 2 2. W 1 W 1 1 MAp _ wave [n]  ( y[n  2 ]  ...  y[n  2 ]) W2 2 2. (3-23). 其中 W1 為 0.055fs，即上述平均 P 波寬度的一半（55ms），使得 MApeak[n]的響應速度較快；而 W2 為 0.11fs，即上述平均 P 波寬度（110ms），因此 MAp_wave[n] 會是響應速度較慢的結果。（見圖 3-9(d)） 2.. 尋找感興趣的區間(BOI) 根據上一段的(3-23)式，將點對點的比較，定義出 BOI Blocks[n]，規則. 如下所示，處理結果見圖 3-8(e). 50.