心率變異性HRV信號提取及時頻域(包含程序)

1、課程設計報告題 目： 心率變異性（HRV）信號的提取及時頻域分析 專 業(yè)： 生物醫(yī)學工程 班 級： XXXXXXX 學 號： XXXXXXX 姓 名： XXXXXXX 指導教師： XXXXXXX XXXXXX大學 XXXXX學院2016年 9月 29日一、 開題背景（一）HRV簡介傳統(tǒng)的醫(yī)學觀點認為，正常的心率為規(guī)則的竇性節(jié)律；后來發(fā)現(xiàn)在健康狀態(tài)下，許多生理系統(tǒng)中存在自然的變異性，人的心率正常情況下也是呈不規(guī)則性變化的，而心率變異就是指竇性心率的這種波動變化的程度。心率變異性（Heart Rate Variability，HRV）是指逐次心搏間期之間的微小變異特性。在生理條件下，HRV的產(chǎn)生主

2、要是由于心臟竇房結(jié)自律活動通過交感和迷走神經(jīng)，神經(jīng)中樞，壓力反射和呼吸活動等因素的調(diào)節(jié)作用，使得心臟每搏間期一般存在幾十毫秒的差異。 (二)HRV的研究現(xiàn)狀心率變異性(HRV)是近年來比較受關注的無創(chuàng)性心電監(jiān)測指標之一，對HRV的生理和病理意義進行了廣泛和深入的研究，其結(jié)果表明心率變異信號中蘊含著有關心血管調(diào)節(jié)的重要信息，對HRV進行分析可以間接地定量評價心肌交感、迷走神經(jīng)的緊張性和均衡性，而且還能分析自主神經(jīng)系統(tǒng)的活動情況。心率變異性還可以作為一個獨立的心源性猝死危險性的預測指標。同時心率變異性分析對多種惡性心律失常的預后判斷和藥物治療效果分析有指導作用。所以，對HRV的研究能夠極大的促進人

3、類對于心血管疾病的了解，從而在預防、治療心血管疾病等領域取得成果。（三）HRV的研究方法隨著對HRV研究的不斷深入，其蘊含的生理病理信息將進一步被揭示，使得HRV有更多的應用空間和應用價值。目前，心率變異性分析方法主要有時域分析法、頻域分析法、時頻分析法以及非線性分析法1。（四）HRV的臨床應用（1）心臟性猝死（SCD）預測：由于HRV是反映自主神經(jīng)張力的最敏感的指標，因此HRV降低是預測心臟性猝死最有價值的獨立指標。（2）急性心肌梗塞后患者危險性評估： HRV的降低是預測急性心肌梗塞后患者發(fā)生心臟性猝死和惡性心律失常危險的重要獨立指標。一般建議在梗塞后一周開始進行HRV的檢測。HRV在梗塞后

4、立即降低，并在幾周內(nèi)開始恢復（2周后逐漸回升），大約6-12個月恢復正常。因此，多次測定HRV可能比單次測定價值更大。梗塞后HRV恢復的快慢對以后死亡的危險性也有預測價值。（3）對糖尿病患者自主神經(jīng)系統(tǒng)損傷的評估：糖尿病患者不論病情輕重，均存在不同程度的自主神經(jīng)功能紊亂。HRV是判斷糖尿病患者是否伴有自主神經(jīng)系統(tǒng)損害最準確，最敏感的指標。（4）心力衰竭（CHF）患者危險性評估。（5）心率變異性生物反饋療法：對于不孕人群受孕幾率提高、懷孕人群孕期焦慮癥改善、產(chǎn)后人群產(chǎn)后抑郁癥情況緩解，起到很好的作用。（6）其它臨床應用范圍：心絞痛、高血壓、心肌病、非缺心臟病所致的慢性嚴重二尖瓣返流、二尖瓣脫垂、

5、心律失常、血管迷走性暈厥等心血管疾病。二、 課題目的（一） 基本掌握心電信號（ECG）的測量、數(shù)據(jù)采集的方法。（二） 學會使用MATLAB對ECG信號進行相關處理分析。主要包括從ECG信號中提取出所需的HRV信號，并分別對其進行時域、頻域、功率譜上的分析。（三） 掌握HRV信號的時頻域參數(shù)的意義，以及對其進行分析的基本方法。三、 課題研究的主要內(nèi)容（一）從網(wǎng)上下載正常人的心電信號以及各種病人的心電信號（ECG）數(shù)據(jù)。（二）首先HRV信號的提取，主要包括去除干擾、準確確定R波波峰位置、剔除異搏、確定R-R間期、線性內(nèi)插，并且繪出HRV信號曲線。（三）對HRV信號的時域分析，對HRV信號的頻譜圖和

6、功率譜圖分析。四、 原理和方法（一） ECG信號的采集 本文主要使用100.hea、100dat、100.atr, 101.hea、101.dat、101.atr, 102.hea、102.dat、102.atr這三組數(shù)據(jù)來對HRV進行研究。實驗數(shù)據(jù)來源于PhysioNet。PhysioNet是一個基于Web的復雜生理和生物醫(yī)學信號的研究資源網(wǎng)站，其網(wǎng)址為http: /www.physionet. org。 PhysioNet由PhysioNet, PhysioBank和PhysioToolkit三個相互關聯(lián)的部分組成。數(shù)據(jù)庫中數(shù)據(jù)來源于正常人、各種病人（如心臟猝死、心力衰竭、心律失常、癲癇、

7、睡眠呼吸暫停綜合癥等）及運動、休息等不同狀態(tài)下的數(shù)據(jù)，樣本選取范圍廣泛，其中大部分數(shù)據(jù)都進行詳細的注釋，并將數(shù)據(jù)被劃分為3類，即Class l：專家已經(jīng)作出了標注；Class 2：原始數(shù)據(jù)；Class3：處于研究進展之中。因此，PhysioBank數(shù)據(jù)庫中的數(shù)據(jù)足已滿足生物醫(yī)學各領域研究者的需要。PhysioBank數(shù)據(jù)庫中的每一條數(shù)據(jù)記錄包括至少三類文件，頭文件（.hea）、數(shù)據(jù)文件（.dat）和注釋文件（.atr,.al,.aiM等）。頭文件是描述數(shù)據(jù)屬性的文本文件，其內(nèi)容包括記錄名、信號數(shù)目、貯存格式、信號數(shù)量和類型、采樣頻率、數(shù)字化特征、記錄的持續(xù)時間和起始時間等信息。一般可由Phys

8、ioToolkit軟件庫的WFDB庫函數(shù)的getinfo、putinfo函數(shù)讀和寫的字符。數(shù)據(jù)文件是定義了相應存貯格式的數(shù)字化采樣點的二進制存儲文件。數(shù)據(jù)存貯格式在頭文件中說明，一個數(shù)據(jù)組有相同的數(shù)據(jù)存貯格式，常用的是8位和16位格式。注釋文件是記錄了對信號特征的注釋信息2。（二） ECG信號的特征（1） 典型心電信號波形心臟搏動及其節(jié)律性是人體生命和生理狀態(tài)的重要標志之一。心電生理學的研究表明，心電信號來源于心肌細胞的生物電變化。心肌細胞的電激動稱為除極，心肌細胞恢復為靜息狀態(tài)稱為復極，心電信號的產(chǎn)生與心肌細胞的除極和復極過程密切相關。心臟電激動起源于竇房結(jié)，沿特化的心臟傳導系統(tǒng)下傳，其傳播

9、方向、途徑、次序及時間存在一定的規(guī)律。若心臟不能及時發(fā)出電激動，則心臟陷于停博。人體體液中充滿電解質(zhì)，具有導電性能，心臟電激動過程產(chǎn)生的有序生物電變化通過體液傳至身體表面使身體各部位出現(xiàn)有規(guī)律而各向異性的電位變化，通過測量電極采集體表特定點電位變化，并放大、顯示及記錄，即為體表心電信號，也即是通常的 ECG 信號3。一個心動周期正常心電信號波形如圖 4.1 所示。它是由特征波及其特征間期組成，每個心動周期包含一個 P 波，一個 QRS 波群和一個 T 波，有時還會出現(xiàn)一個小的 U 波。特征波及特征間期的含義如下：圖 4.1 典型的心電信號波形P 波：由左右心房的除極過程引起，其波形小而圓鈍，時

10、寬為 0.08s-0.11s，波幅不超過 0.25mV。QRS 波：反映左右心室除極產(chǎn)生的電位變化，在 P 波之后出現(xiàn)，為心電信號中最高大和最快速的波形。典型的 QRS 波群包括三個相連的波，第一個向下的波為 Q 波，緊接著為高而尖峭的向上的 R 波，最后是一個向下的 S 波。在體表不同位置（使用不同導聯(lián)記錄）時，三個波不一定都有，大小方向也會不同。T 波：代表心室復極時的電位變化，方向與 R 波方向相同，時寬為 0.05s-0.25s，波幅一般為 0.lmV-0.8mV。U 波：T 波之后可能出現(xiàn)的一個低而寬的波，與 T 波方向一致，其機理不十分清楚，可能反映普頃野纖維復極的電位變化。QRS

11、 間期：從 Q 波開始至 S 波結(jié)束之間的時程，反映心室除極時間。正常 QRS 間期為 0.04s-0.1s，若 QRS 間期>0.12s，反映室內(nèi)傳導阻滯。PR 間期：從 P 波開始到 QRS 波開始之間的時程，反映激動由竇房結(jié)產(chǎn)生經(jīng)由結(jié)間束、房室交界和左右束支抵達心室，并引起心室興奮所需要的時間，又稱為房室傳導時間。正常為 0.12s-0.2s。當發(fā)生房室傳導阻滯時，PR 間期增長。如當 PR 間期>0.21s。則為度房室傳導阻滯。QT 間期：從 QRS 波開始到 T 波終點的時程，反映心室除極和復極時間的總和。許多因素可影響 QT 間期，如心肌缺血、低血鉀、低血鈣等可使 QT

12、 間期延長，QT 間期延長使心室肌復極不均一，易誘發(fā)折返激動，導致嚴重室性心率失常。QT 間期隨受心率變化的影響，心率越慢，QT 間期越長；心率越快，QT 間期越短。通常用 QTc間期修正心率對 QT 間期影響，正常 QTc間期小于 0.430.44s。ST 段：指從 QRS 波群終止點到 T 波起點之間的波形線段，反映心室部分己完全進入去極化狀態(tài)，正常時與基線平齊。PP 間期：相鄰 P 波之間的間距稱為 PP 間期，反映心房率。正常情況下，PP 間期與 RR 間期一致。在度或度以上房室傳導阻滯和某些心率失常，兩者可不一致。RR 間期：相鄰 QRS 波群之間的間距稱為 RR 間期，反映心室率。

13、正常情況下，RR間期與 PP 間期一致2。在心電信號的測試中，對電極的放置部位和導聯(lián)的連接方式臨床有明確的規(guī)定。目前，國際公認的是標準 12 導聯(lián)，包括心電標準導聯(lián)（I、II、III）、加壓單極肢體導聯(lián)（aVR、aVL、aVF）及胸導聯(lián)（VlV6），共有 12 個導聯(lián)，具體可參考文獻4。(三)典型心電信號的能量（頻譜圖）分布典型的心電信號的整個心動周期的頻譜估計圖如圖4.22所示，可以明顯看出心電信號各波的能量主要集中在低頻區(qū)域，且隨著頻率的增高，相應的能量逐漸降低。心電信號的整體頻譜范圍在0.05Hz100Hz，但能量主要集中在0.545Hz，能量的最高點在815Hz附近；QRS 波群的頻譜

14、帶寬為340 Hz，積聚了將近99%的能量，波峰能量集中在618Hz附近, P波的頻譜帶寬為018Hz，波峰能量集中在512Hz；T波的頻譜帶寬為08Hz，波峰能量集中在08Hz區(qū)間5。圖 4.2 典型的心電信號頻譜能量分布(四)ECG信號的噪聲分析在采集、放大及傳輸心電信號的過程中，由于受人體、采集儀器、電磁環(huán)境、操作水平等的影響，不可避免會有許多干擾耦合到心電信號，主要干擾表現(xiàn)形式如下：（1）電源工頻干擾產(chǎn)生的原因主要由于電源磁場作用于心電圖儀的導聯(lián)與人體之間的環(huán)形電路所致，表現(xiàn)為心電信號上有明顯的正弦波或正弦波的疊加信號，其頻率為60Hz工頻及其諧波構(gòu)成，幅度較低。（2）基線漂移產(chǎn)生的原

15、因主要由于人體呼吸運動、電極接觸不良等因素所導致。表現(xiàn)為心電信號上疊加緩慢變化的信號，其頻率一般小于1Hz，幅度為ECG峰-峰的15%。（3）肌電干擾產(chǎn)生的原因主要由于人體活動，肌肉緊張所引起的干擾。表現(xiàn)為不規(guī)則的快速變化波形，其頻率范圍較寬，一般在52kHz之間，幅度為毫伏級。（4）運動偽跡產(chǎn)生的原因主要由于電極與人體間輕微移動或抖動而引入的干擾，表現(xiàn)為信號基線的短暫變化，但不是基線的躍變，其持續(xù)時間為100500ms，頻率一般在7Hz以下，幅度較大。（5）其他隨機噪聲心電信號還受到其他的隨機噪聲和環(huán)境干擾的影響,如加性白噪聲、極化噪聲、儀器內(nèi)部噪聲等。由此可見，噪聲信號基本覆蓋了有用的心電

16、信號的全頻率范圍。而其中60Hz及其倍頻附近的工頻干擾、1Hz以下的基線漂移以及肌電干擾噪聲是最主要的干擾源，在心電信號預處理中必須消除或抑制,以提高心電信號的信噪比2。（五）ECG信號去噪聲方法對于心電信號預處理一般要從硬件電路優(yōu)化設計和軟件數(shù)字濾波器的設計兩個方面考慮。根據(jù)心電信號的頻譜分布特點，在硬件方面來看，消除基線漂移的干擾應考慮分別設計下限頻率為0.5Hz的高通濾波器；消除肌電等高頻干擾應考慮設計上限頻率為100Hz的低通濾波器，同時還應考慮設計60Hz的陷波器來濾除工頻干擾。由于硬件濾波器的元器件精度、穩(wěn)定性要求以及難以實現(xiàn)嚴格線性相位等問題，使得僅采用硬件濾波不能滿足濾波性能的

17、要求。隨著數(shù)字信號處理技術的發(fā)展，設計高精度、高可靠性、簡潔靈活的數(shù)字濾波器成為可能。目前，濾除心電信號干擾更多采用的是具有線性相位的數(shù)字濾波方法，并已逐漸顯示出取代硬件濾波器的趨勢。（1）工頻干擾對于工頻干擾，主要的消除方法有平均濾波器、梳狀濾波器、Levkov 濾波器、自適應工頻陷波器等方法。平均濾波器法是較早被應用的數(shù)字濾波方法，其特點是算法簡單，處理速度快，濾波效果較好。由于平均濾波算法的實時性好，被廣泛應用于實時心電監(jiān)護設備中或基于單片機心電數(shù)據(jù)采集系統(tǒng)中。從實際濾波效果看，濾波器能較好地濾除了60Hz 工頻，但對 ECG 中的 QRS 波有較大削峰，信號衰減較大。因此，在實際應用中

18、受到一定的限制。梳狀濾波器的名稱源于該濾波器的傳遞函數(shù)的幅頻特性形如梳狀。其特點是運算量小、形式靈活簡單，是一種快速的數(shù)字濾波器。梳狀濾波器被廣泛應用去除工頻及其高次諧波的干擾。選擇合適的梳狀濾波參數(shù)，可保證濾波器對有用信號基本無衰減和濾波后基本不產(chǎn)生相移6。改進的梳狀濾波器可用于消除 ECG 的工頻干擾且能有效地防止信號的失真。Levkov 濾波器是 Levkov 在 1984 年提出一種對 ECG 信號的線性段和非線性段采用不同處理方式的數(shù)字濾波方法。其基本原理是從原始 ECG 信號中直接減去在該線性段中確定的干擾信號的幅值7。改進的 Levkov 濾波方法中引入了判定線性段和非線性段的指

19、標 M。算法的濾波效果與選擇合適的 M 值有關。改進 Levkov 濾波法算法簡單，參數(shù)可調(diào)，運算量小，可實現(xiàn)實時心電信號濾波處理，但在噪聲較強時，將會造成 QRS 波群的削峰現(xiàn)象，濾波效果不能令人滿意8。自適應工頻陷波技術特點是只需要很少或根本不需要任何有關信號和噪聲的先驗知識的情況下，直接利用觀測數(shù)據(jù)不斷遞歸更新自身的參數(shù)，以逐步逼近某一最優(yōu)值，能自動跟蹤信號頻率漂移，具有非常好的適應性9。自適應的數(shù)字陷波器的權值調(diào)節(jié)主要有最小均方誤差（LMS）、最小遞推二乘法(RLS)算法10。一般來說，RLS 自適應算法與 LMS 自適應算法相比，權值有更快的收斂速度。缺點是每次迭代計算量較大，由于在

20、生物醫(yī)學應用中數(shù)據(jù)的計算量不是太大, RLS 自適應算法可獲得更好性能。（2）基線漂移的消除基線漂移可嚴重影響心電信號的分析與處理。濾波法和擬合基漂法是去除心電基線漂移的兩類主要算法。具體而言有形態(tài)濾波、中值濾波、線性相位FIR濾波、小波變換、形態(tài)和小波相結(jié)合、中值濾波和小波相結(jié)合等去除基漂的濾波方法，以及導數(shù)法與坐標法這兩類基漂擬合方法。研究表明濾波法實現(xiàn)簡單,但精度較低，分段擬合法去除基漂的效果更好,有自適應的特性而且失真較小。無論是實時性還是準確性,基漂擬合法都比濾波法更有優(yōu)勢但某些場合下,基漂擬合點的提取十分困難,此時只能使用濾波法11。對于基線擬合法，通常多采用簡潔、快速的分段擬合法

21、。其基本原理是通過對采樣的心電數(shù)據(jù)分段進行三次多項式或拋物線擬合，獲取基線漂移曲線，然后用原始心電曲線減去擬合曲線來消除基線漂移。算法的關鍵是插值點選擇，插值點應該是心電信號的零電位點，通常采用 TP 段。對于FIR 濾波器只有零點，不易獲得較好的通帶與阻帶衰減特性。要取得好的衰減特性，一般要求 H(z)的階次要高。對于小波變換，由于其提供了一個在時頻平面上可調(diào)的分析窗口。信號的時、頻分辨率可以隨分析任務的需要作出相應調(diào)整。在低頻時小波變換的時間分辨率較差，而頻率分辨率較高；在高頻時小波變換的時間分辨率較高，而頻率分辨率較低，自動滿足了信號分析。信號x(t)在不同尺度分解的逼近信號和細節(jié)信號具

22、有不同頻帶，而在不同的尺度上一定包含有x(t)的不同頻率信息。所以若對包含有不同頻率信息的細節(jié)信號進行閾值的取舍，丟棄其中的干擾成分后，再對信號進行重構(gòu)，則重構(gòu)信號的干擾得到了非常好的抑制，即為小波變換消噪的基本原理2。通過對集中方法的實驗效果來選擇合適的消除基線漂移的方法。（六）R波峰位置、剔除異搏、確定RR間期（1）波峰位置找出波峰位置的基本想法是尋找一個采樣點，其電壓值既大于前一個采樣點的值也大于后一個采樣點的值。之后再將P波與T波中的偽峰值去除，剩下的就是QRS波的峰值，即R波峰?？梢酝ㄟ^先屏蔽掉小于某一電壓的所有采樣點，這樣就將P波和T波中的偽峰值提前去除，之后再尋找峰值。以此保證得

23、到的峰值處于QRS波中，該峰值即為R波峰值。（2） 剔除異搏對得到R波的位置，檢查其中的異常點，去掉異搏。（3） 確定RR間期用后一個峰值的采樣時間減去前一個峰值的采樣時間就可以得到RR間期的值。將信號以RR值為縱坐標，R-Ri值為橫坐標，得到HRV的圖像。（七）HRV信號的時域分析通常使用的HRV時域檢測指標有5項：NNVGR、SDNN、RMSSD、SDSD和pNN50。上述5個指標的定義分別為：NNVGR：全部正常竇性心博間期（NN）的平均值，單位為ms。即。SDNN：標準差，即全部NN間期的標準差，單位為ms。即。RNSSD：全程相鄰NN間期之差的均方根值，單位為ms。即 。SDSD：全

24、程相鄰NN間期長度之差的標準差，單位為ms。即其中，。pNN50：在全部NN間期的記錄中，相鄰的NN間期之差大于50ms的個數(shù)與總的NN間期的個數(shù)的比，以百分比表示。其中NNVGR用于評估心率總體變化水平,；SDNN用于評估心率總體變化的大小，即交感及迷走神經(jīng)張力大小,；SDANN用于評估心率變化中的長期慢變化成分，即交感神經(jīng)張力大小。RMSSD及pNN50反映心率快變化成分的大小，即副交感神經(jīng)張力的敏感指標。(八) HRV信號的頻域分析：HRV時域分析的指標大多用于描述HRV整體的大小，不能仔細地分析交感神經(jīng)和迷走神經(jīng)各自的活動的情況。而HRV頻譜分析則可以彌補這個缺點。另一角度，即頻譜分析

25、的角度來分析心率變化的規(guī)律。它與時域分析既有相關性，又能揭示出心率的更復雜的變化規(guī)律。FFT是經(jīng)典譜估計方法，算法簡單。輸入和輸出信號能量有線性關系，但對信號要作周期延拓假定，短數(shù)據(jù)譜分辨率較低，并有能量泄露現(xiàn)象。最近更仔細的研究發(fā)現(xiàn)，正常人基礎狀態(tài)下心率譜曲線在0-0.4Hz之間，0.003-0.04Hz為極低頻段（VLF），0.04-0.15Hz為低頻段（LF），0.15-0.4Hz高頻段（HF），0-0.40Hz為總功率譜（TP）。研究表明，VLF反映心率變化受熱調(diào)節(jié)（體溫），血管舒縮張力和腎血管緊張素系統(tǒng)的影響；LF反映交感和迷走神經(jīng)的雙重調(diào)節(jié)；HF只反映迷走神經(jīng)的調(diào)節(jié)；TP反映HRV

26、大小，LF/HF比值反映自主神經(jīng)系統(tǒng)的平衡狀態(tài)，基本上代表經(jīng)張交感神力的高低。五、 步驟（包括計算程序流程框圖）（一） 總程序流程圖。詳細程序間附件1。獲取心電信號去除工頻干擾和基線漂移確定波峰位置剔除異搏 確定R-R間期求出HRV信號進行時域計算分析進行線性內(nèi)插進行頻域計算分析進行功率譜分析HRV信號NNVGRSDNNRMSSDSDSDpNN50頻譜圖TPHFLFLF/HF結(jié)論（二） 獲取心電信號將下載得到三組心電信號：100、101、102，通過Klaus Rheinberger編寫的讀取心電的MATLAB程序ecg.load.m，加載原始心電信號，從而準確提取相應的心電信號。（三） 去除

27、工頻噪聲干擾對信號進行頻域分析，發(fā)現(xiàn)信號主要集中在低頻端，并且在60Hz處有微弱信號，通過判斷可得出該頻率信號為工頻噪聲。本文使用FIR濾波器將噪聲濾去。（四） 去基線漂移對比兩大類型發(fā)現(xiàn)：濾波法實現(xiàn)簡單，但精度較低；分段擬合法去除基漂的效果更好,有自適應的特性而且失真較小。無論是實時性還是準確性,基漂擬合法都比濾波法更有優(yōu)勢；但某些場合下,基漂擬合點的提取十分困難,此時只能使用濾波法。通過進行初步實驗中比較各種方法的實際效果，發(fā)現(xiàn)使用中值濾波法的效果最好，故本文主要采用其進行去基線漂移處理。（五） 找出波峰、剔除異搏、確定RR間期（六） 對HRV信號進行時域分析對前面得到的HRV信號進行時域

28、分析，其中的技術指標有NNVGR、SDNN、RNSSD、SDSD、pNN50等。(七) 線性內(nèi)插對提取出來的HRV信號經(jīng)線性內(nèi)插獲得等間隔（R-R間期均值）R-R間期時間序列，然后再進行頻譜分析。線性內(nèi)插的公式為。式中X(n)為插值，和分別為插值前、后的R-R間期序列值，和分別為X（n）所在位置與和所在位置的時間間距。（八）根據(jù)快速傅里葉變換（FFT）得到HRV信號的頻譜信號和功率譜信號。主要有TP、HF、LF、LF/HF、頻譜圖、功率譜圖等。六、 結(jié)果通過MATLAB程序分別得到100、101和102三組ECG信號數(shù)據(jù)，然后對其進行相關分析，具體結(jié)果以及圖像如下所述。（一）ECG信號的提取通

29、過MATLAB程序分別得到100、101和102三組ECG信號數(shù)據(jù)。其時序圖像具體如圖6.1所示。 圖 6.1 三組原始心電信號 圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（二）原始心電頻譜圖對三組信號通過皮譜圖進行頻域分析。從圖6.2中可以看出，信號主要集中在20Hz以下的低頻端，并且在60Hz附近有一個明顯的幅度值上升，可以判斷其為工頻干擾。 圖 6.2 三組心電信號的頻譜圖圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（三） 去除工頻干擾使用FIR濾波器法將工頻噪聲去除。去除工頻干擾之后頻譜圖如圖6.3所示圖 6.3 三組心電信號去除工頻干擾后

30、的頻譜圖圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（四） 去除基線漂移從圖6.1中可以看出這些心電信號均有明顯的基線漂移。使用中值濾波法去除基線漂移。消除基線漂移后的ECG信號時序圖如圖6.4所示。圖 6.4 消除基線后三組心電信號的頻譜圖圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（五） 確定波峰確定波峰結(jié)果如圖6.5所示。其中紅色點表示波峰。圖 6.5確定三組心電信號的R波峰圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（六） 確定R-R間期R-R間期與R-Ri的關系，亦即HRV的圖像如圖6.6所示。圖 6.6 三組心電信號

31、的HRV圖像圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（七）時域分析HRV的時域分析指標NNVGR、SDNN、RMSSD、SDSD、NN50和pNN50計算結(jié)果如表6.1所示。表6.1 三組數(shù)據(jù)的時域指標值數(shù)據(jù)組 NNVGR(ms) SDNN(ms) RMSSD(ms) SDSD(ms) NN50 pNN50100. 811.4 34.2 47.5 39.9 7 6.93%101. 859.3 48.8 28.6 17.5 8 8.42%102. 818.4 143.7 199.2 187.5 11 11.34%（八）頻域分析根據(jù)快速傅里葉變換（FFT）得到HRV信號的

32、頻譜信號和功率譜信號。（1） 頻譜圖。三組心電數(shù)據(jù)對應的HRV信號的頻譜圖如圖6.7所示。圖 6.7 消除基線后三組心電信號的頻譜圖圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（2） 功率譜圖圖 6.8 HRV信號的功率譜圖圖示 由上至下分別為來自100號、101號、102號數(shù)據(jù)的對應圖像（3） 根據(jù)頻譜計算TP、HF、LF、LF/HF，結(jié)果如下表6.1所示。表6.1 三組數(shù)據(jù)的時域指標值數(shù)據(jù)組 TP(ms×ms) HF(ms×ms) LF(ms×ms) LF/HF100.dat 19000 1400 571 0.4150101.dat 20

33、600 871 1100 1.2063102.dat 28300 6400 3700 0.5812七、 結(jié)論本文通過使用MATLAB軟件，實現(xiàn)了心電信號的讀取、采樣及頻譜分析，并且分離出心率變異性HRV信號，通過時域分析和頻域分析的方法進行分析，對比資料發(fā)現(xiàn)：第101號數(shù)據(jù)組應該來自正常人，第100號數(shù)據(jù)組應該來自具有輕微患心血管疾病的患者，第102號數(shù)據(jù)組則最可能來自患有較為嚴重的心血管疾病的患者八、 討論心率變異信號中蘊含著有關心血管調(diào)節(jié)的重要信息，對HRV進行分析可以間接地定量評價心肌交感、迷走神經(jīng)的緊張性和均衡性，而且還能分析自主神經(jīng)系統(tǒng)的活動情況。同時心率變異性分析對多種惡性心律失常

34、的預后判斷和藥物治療效果分析有指導作用。所以，對HRV的研究能夠極大的促進人類對于心血管疾病的了解，從而在預防、治療心血管疾病等領域取得成果。本文從獲取的心電信號中，提取出有效的心電信號，進而提取出我們所研究的HRV信號，然后對其進行相關分析。在獲取心電信號的過程中，直接使用了來自美國PhysioNet的心電信號資源，在提取出有效心電信號的過程中，首先進行的是對原始心電信號，即三組每組三個相關文件進行相關處理，得出心電信號；其次是觀察得到的初步心電信號進行時序圖的觀察，再進行去工頻干擾和去基線漂移；再其次，對前面得到的信號確定波峰位置，并且剔除異搏，確定出R-R間期；最后，得到我們所需要的HR

35、V信號。對于HRV信號的分析主要分為時域和頻域兩大部分，時域分析中，本文主要對HRV信號的五個時域指標NNVGR、SDNN、RMSSD、SDSD、pNN50進行了計算；對于頻域部分，對HRV信號的頻域四大指標TP、HF、LF、LF/HF進行了計算，然后作出了HRV信號的頻譜圖和功率譜圖。在研究中使用多種方法進行對信號的處理，經(jīng)過對比各種方法得到的數(shù)據(jù)對比，本文選取了效果相對而言比較好，算法相對而言比較簡單的一些方法，比如心電數(shù)據(jù)的去工頻干擾中使用的是60Hz陷波器濾波法，在心電信號的去基線漂移中使用的是中值濾波法，之后為了保證能對HRV信號進行正確的頻域分析，對HRV信號進行了線性內(nèi)插。以上種

36、種方法都是為了使數(shù)據(jù)具有更高的準確性，基本確保研究沒有在研究對象的處理上的重大失誤。在前人的基礎之上，由于水平、能力有限，本文對心電信號的HRV進行了粗淺的研究。在心電數(shù)據(jù)的處理上和認識上都還有很多的東西都需要學習，希望以后能對此做進一步的研究。參考文獻1 鐘運健. 心率變異性(HRV)在運動性疲勞診斷中應用的實驗研究:碩士論文.江西南昌：江西師范大學體育學院，2004年.2 董紅生.心電波形檢測與心率異變性分析研究:博士論文.蘭州理工大學，2012年3 吳學勤.動態(tài)心電圖技術與應用.合肥:中國科學技術出版社,1998,22-25.4 胡大一,郭成軍,李瑞杰.心率變異性一測量標準,生理釋義與臨

37、床應用續(xù)一:心率變異性的生理研究.中國醫(yī)療器械信息,1998, 4(1):29-32.5 曹細武,鄧親愷.心電圖各波的頻率分析.中國醫(yī)學物理學雜志,2001,18(1):46-48.6 萬相奎.心電信號分析與虛擬式心電自動分析儀的開發(fā):重慶大學博士學位論文.重慶:重慶大學機械電子工程,2005,16-17.7 Levkov C,Michov G, Ivanov Ret al. Subtraction of 50Hz interference from the electrocardiogram.Med.&Biol., Eng.&Comp.,1984,22(4):371-373.

38、8 孫京霞, 白延強, 楊玉星. 一種抑制心電信號 50Hz 工頻干擾的改進 Levkov 方法.航天醫(yī)學與醫(yī)學工程,2000,13(3):196-199.9 楊福生,呂揚生.生物醫(yī)學信號的處理和識別.天津:天津科技翻譯出版司,1997,378-379.10 高鷹,謝勝利.一種變步長LMS 自適應濾波算法及分析. 電子學報.2001,29(8):1094-1097.11 李延軍,嚴 洪,王增麗.心電基線漂移去除方法的比較研究.航天醫(yī)學與醫(yī)學工程.2009年第22卷第五期：381-386.附件1：% This programm reads ECG data which are saved in

39、format 212.% (e.g., 100.dat from MIT-BIH-DB, cu01.dat from CU-DB,.)% The data are displayed in a figure together with the annotations.% The annotations are saved in the vector ANNOT, the corresponding% times (in seconds) are saved in the vector ATRTIME.% The annotations are saved as numbers, the mea

40、ning of the numbers can% be found in the codetable "ecgcodes.h" available at .% ANNOT only contains the most important information, which is displayed% with the program rdann (available on ) in the 3rd row.% The 4th to 6th row are not saved in ANNOT.% crea

41、ted on Feb. 27, 2003 by% Robert Tratnig (Vorarlberg University of Applied Sciences)% (email: rtratniggmx.at),% algorithm is based on a program written by% Klaus Rheinberger (University of Innsbruck)% (email: klaus.rheinbergeruibk.ac.at)% -clc; clear all; %- SPECIFY DATA -PATH= 'C:Usersasus你Deskt

42、op心率變異性（HRV）信號的提取及時頻域分析ecg' % path, where data are savedHEADERFILE= '100.hea' % header-file in text formatATRFILE= '100.atr' % attributes-file in binary formatDATAFILE='100.dat' % data-fileSAMPLES2READ=30000; % number of samples to be read % in case of more than one signa

43、l: % 2*SAMPLES2READ samples are read %- LOAD HEADER DATA -fprintf(1,'n$> WORKING ON %s .n', HEADERFILE);signalh= fullfile(PATH, HEADERFILE);fid1=fopen(signalh,'r');z= fgetl(fid1);A= sscanf(z, '%*s %d %d %d',1,3);nosig= A(1); % number of signalssfreq=A(2); % sample rate of

44、dataclear A;for k=1:nosig z= fgetl(fid1); A= sscanf(z, '%*s %d %d %d %d %d',1,5); dformat(k)= A(1); % format; here only 212 is allowed gain(k)= A(2); % number of integers per mV bitres(k)= A(3); % bitresolution zerovalue(k)= A(4); % integer value of ECG zero point firstvalue(k)= A(5); % firs

45、t integer value of signal (to test for errors)end;fclose(fid1);clear A; %- LOAD BINARY DATA -if dformat= 212,212, error('this script does not apply binary formats different to 212.'); end;signald= fullfile(PATH, DATAFILE); % data in format 212fid2=fopen(signald,'r');A= fread(fid2, 3,

46、 SAMPLES2READ, 'uint8')' % matrix with 3 rows, each 8 bits long, = 2*12bitfclose(fid2);M2H= bitshift(A(:,2), -4);M1H= bitand(A(:,2), 15);PRL=bitshift(bitand(A(:,2),8),9); % sign-bitPRR=bitshift(bitand(A(:,2),128),5); % sign-bitM( : , 1)= bitshift(M1H,8)+ A(:,1)-PRL;M( : , 2)= bitshift(M2

47、H,8)+ A(:,3)-PRR;if M(1,:) = firstvalue, error('inconsistency in the first bit values'); end;switch nosigcase 2 M( : , 1)= (M( : , 1)- zerovalue(1)/gain(1); M( : , 2)= (M( : , 2)- zerovalue(2)/gain(2); TIME=(0:(SAMPLES2READ-1)/sfreq;case 1 M( : , 1)= (M( : , 1)- zerovalue(1); M( : , 2)= (M(

48、: , 2)- zerovalue(1); M=M' M(1)=; sM=size(M); sM=sM(2)+1; M(sM)=0; M=M' M=M/gain(1); TIME=(0:2*(SAMPLES2READ)-1)/sfreq;otherwise % this case did not appear up to now! % here M has to be sorted! disp('Sorting algorithm for more than 2 signals not programmed yet!');end;clear A M1H M2H

49、PRR PRL;fprintf(1,'n$> LOADING DATA FINISHED n'); %- LOAD ATTRIBUTES DATA -atrd= fullfile(PATH, ATRFILE); % attribute file with annotation datafid3=fopen(atrd,'r');A= fread(fid3, 2, inf, 'uint8')'fclose(fid3);ATRTIME=;ANNOT=;sa=size(A);saa=sa(1);i=1;while i<=saa ann

50、oth=bitshift(A(i,2),-2); if annoth=59 ANNOT=ANNOT;bitshift(A(i+3,2),-2); ATRTIME=ATRTIME;A(i+2,1)+bitshift(A(i+2,2),8)+. bitshift(A(i+1,1),16)+bitshift(A(i+1,2),24); i=i+3; elseif annoth=60 % nothing to do! elseif annoth=61 % nothing to do! elseif annoth=62 % nothing to do! elseif annoth=63 hilfe=bi

51、tshift(bitand(A(i,2),3),8)+A(i,1); hilfe=hilfe+mod(hilfe,2); i=i+hilfe/2; else ATRTIME=ATRTIME;bitshift(bitand(A(i,2),3),8)+A(i,1); ANNOT=ANNOT;bitshift(A(i,2),-2); end; i=i+1;end;ANNOT(length(ANNOT)=; % last line = EOF (=0)ATRTIME(length(ATRTIME)=; % last line = EOFclear A;ATRTIME= (cumsum(ATRTIME)

52、/sfreq;ind= find(ATRTIME <= TIME(end);ATRTIMED= ATRTIME(ind);ANNOT=round(ANNOT);ANNOTD= ANNOT(ind); %- DISPLAY DATA -figure(1); clf, box on, hold onplot(TIME, M(:,1),'r'); %修改,'LineWidth',2DATA1=M(:,1)'% if nosig=2% plot(TIME, M(:,2),'b');% end;% for k=1:length(ATRTIMED)% text(AT