心率變異性HRV信號(hào)提取及時(shí)頻域分析(包含程序)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

25、的角度來(lái)分析心率變化的規(guī)律。它與時(shí)域分析既有相關(guān)性，又能揭示出心率的更復(fù)雜的變化規(guī)律。FFT是經(jīng)典譜估計(jì)方法，算法簡(jiǎn)單。輸入和輸出信號(hào)能量有線性關(guān)系，但對(duì)信號(hào)要作周期延拓假定，短數(shù)據(jù)譜分辨率較低，并有能量泄露現(xiàn)象。最近更仔細(xì)的研究發(fā)現(xiàn)，正常人基礎(chǔ)狀態(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)張交感神力的高低。五、 步驟（包括計(jì)算程序流程框圖）（一） 總程序流程圖。詳細(xì)程序間附件1。獲取心電信號(hào)去除工頻干擾和基線漂移確定波峰位置剔除異搏 確定R-R間期求出HRV信號(hào)進(jìn)行時(shí)域計(jì)算分析進(jìn)行線性?xún)?nèi)插進(jìn)行頻域計(jì)算分析進(jìn)行功率譜分析HRV信號(hào)NNVGRSDNNRMSSDSDSDpNN50頻譜圖TPHFLFLF/HF結(jié)論（二） 獲取心電信號(hào)將下載得到三組心電信號(hào)：100、101、102，通過(guò)Klaus Rheinberger編寫(xiě)的讀取心電的MATLAB程序ecg.load.m，加載原始心電信號(hào)，從而準(zhǔn)確提取相應(yīng)的心電信號(hào)。（三） 去除

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

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

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

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

31、的HRV圖像圖示 由上至下分別為來(lái)自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像（七）時(shí)域分析HRV的時(shí)域分析指標(biāo)NNVGR、SDNN、RMSSD、SDSD、NN50和pNN50計(jì)算結(jié)果如表6.1所示。表6.1 三組數(shù)據(jù)的時(shí)域指標(biāo)值數(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信號(hào)的

32、頻譜信號(hào)和功率譜信號(hào)。（1） 頻譜圖。三組心電數(shù)據(jù)對(duì)應(yīng)的HRV信號(hào)的頻譜圖如圖6.7所示。圖 6.7 消除基線后三組心電信號(hào)的頻譜圖圖示 由上至下分別為來(lái)自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像（2） 功率譜圖圖 6.8 HRV信號(hào)的功率譜圖圖示 由上至下分別為來(lái)自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像（3） 根據(jù)頻譜計(jì)算TP、HF、LF、LF/HF，結(jié)果如下表6.1所示。表6.1 三組數(shù)據(jù)的時(shí)域指標(biāo)值數(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é)論本文通過(guò)使用MATLAB軟件，實(shí)現(xiàn)了心電信號(hào)的讀取、采樣及頻譜分析，并且分離出心率變異性HRV信號(hào)，通過(guò)時(shí)域分析和頻域分析的方法進(jìn)行分析，對(duì)比資料發(fā)現(xiàn)：第101號(hào)數(shù)據(jù)組應(yīng)該來(lái)自正常人，第100號(hào)數(shù)據(jù)組應(yīng)該來(lái)自具有輕微患心血管疾病的患者，第102號(hào)數(shù)據(jù)組則最可能來(lái)自患有較為嚴(yán)重的心血管疾病的患者八、 討論心率變異信號(hào)中蘊(yùn)含著有關(guān)心血管調(diào)節(jié)的重要信息，對(duì)HRV進(jìn)行分析可以間接地定量評(píng)價(jià)心肌交感、迷走神經(jīng)的緊張性和均衡性，而且還能分析自主神經(jīng)系統(tǒng)的活動(dòng)情況。同時(shí)心率變異性分析對(duì)多種惡性心律失常

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

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

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

37、床應(yīng)用續(xù)一:心率變異性的生理研究.中國(guó)醫(yī)療器械信息,1998, 4(1):29-32.5 曹細(xì)武,鄧親愷.心電圖各波的頻率分析.中國(guó)醫(yī)學(xué)物理學(xué)雜志,2001,18(1):46-48.6 萬(wàn)相奎.心電信號(hào)分析與虛擬式心電自動(dòng)分析儀的開(kāi)發(fā):重慶大學(xué)博士學(xué)位論文.重慶:重慶大學(xué)機(jī)械電子工程,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 孫京霞, 白延強(qiáng), 楊玉星. 一種抑制心電信號(hào) 50Hz 工頻干擾的改進(jìn) Levkov 方法.航天醫(yī)學(xué)與醫(yī)學(xué)工程,2000,13(3):196-199.9 楊福生,呂揚(yáng)生.生物醫(yī)學(xué)信號(hào)的處理和識(shí)別.天津:天津科技翻譯出版司,1997,378-379.10 高鷹,謝勝利.一種變步長(zhǎng)LMS 自適應(yīng)濾波算法及分析. 電子學(xué)報(bào).2001,29(8):1094-1097.11 李延軍,嚴(yán) 洪,王增麗.心電基線漂移去除方法的比較研究.航天醫(yī)學(xué)與醫(yī)學(xué)工程.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）信號(hào)的提取及時(shí)頻域分析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