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

./課程設(shè)計(jì)報(bào)告題目：心率變異性〔HRV信號(hào)的提取及時(shí)頻域分析專業(yè)：生物醫(yī)學(xué)工程班級(jí)：XXXXXXX學(xué)號(hào)：XXXXXXX姓名：XXXXXXX指導(dǎo)教師：XXXXXXXXXXXXX大學(xué)XXXXX學(xué)院2016年9月29日開題背景〔一HRV簡(jiǎn)介傳統(tǒng)的醫(yī)學(xué)觀點(diǎn)認(rèn)為,正常的心率為規(guī)則的竇性節(jié)律；后來發(fā)現(xiàn)在健康狀態(tài)下,許多生理系統(tǒng)中存在自然的變異性,人的心率正常情況下也是呈不規(guī)則性變化的,而心率變異就是指竇性心率的這種波動(dòng)變化的程度。心率變異性〔HeartRateVariability,HRV是指逐次心搏間期之間的微小變異特性。在生理?xiàng)l件下,HRV的產(chǎn)生主要是由于心臟竇房結(jié)自律活動(dòng)通過交感和迷走神經(jīng),神經(jīng)中樞,壓力反射和呼吸活動(dòng)等因素的調(diào)節(jié)作用,使得心臟每搏間期一般存在幾十毫秒的差異。<二>HRV的研究現(xiàn)狀心率變異性<HRV>是近年來比較受關(guān)注的無創(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)人類對(duì)于心血管疾病的了解,從而在預(yù)防、治療心血管疾病等領(lǐng)域取得成果?！踩鼿RV的研究方法隨著對(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)。一般建議在梗塞后一周開始進(jìn)行HRV的檢測(cè)。HRV在梗塞后立即降低,并在幾周開始恢復(fù)〔2周后逐漸回升,大約6-12個(gè)月恢復(fù)正常。因此,多次測(cè)定HRV可能比單次測(cè)定價(jià)值更大。梗塞后HRV恢復(fù)的快慢對(duì)以后死亡的危險(xiǎn)性也有預(yù)測(cè)價(jià)值?！?對(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)重二尖瓣返流、二尖瓣脫垂、心律失常、血管迷走性暈厥等心血管疾病。課題目的基本掌握心電信號(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)行分析的基本方法。課題研究的主要容〔一從網(wǎng)上下載正常人的心電信號(hào)以及各種病人的心電信號(hào)〔ECG數(shù)據(jù)?！捕紫菻RV信號(hào)的提取,主要包括去除干擾、準(zhǔn)確確定R波波峰位置、剔除異搏、確定R-R間期、線性插,并且繪出HRV信號(hào)曲線?！踩龑?duì)HRV信號(hào)的時(shí)域分析,對(duì)HRV信號(hào)的頻譜圖和功率譜圖分析。原理和方法ECG信號(hào)的采集本文主要使用100.hea、100dat、100.atr,101.hea、101.dat、101.atr,102.hea、102.dat、102.atr這三組數(shù)據(jù)來對(duì)HRV進(jìn)行研究。實(shí)驗(yàn)數(shù)據(jù)來源于PhysioNet。PhysioNet是一個(gè)基于Web的復(fù)雜生理和生物醫(yī)學(xué)信號(hào)的研究資源,其網(wǎng)址為http://.physionet.org。PhysioNet由PhysioNet,PhysioBank和PhysioToolkit三個(gè)相互關(guān)聯(lián)的部分組成。數(shù)據(jù)庫中數(shù)據(jù)來源于正常人、各種病人〔如心臟猝死、心力衰竭、心律失常、癲癇、睡眠呼吸暫停綜合癥等及運(yùn)動(dòng)、休息等不同狀態(tài)下的數(shù)據(jù),樣本選取圍廣泛,其部分?jǐn)?shù)據(jù)都進(jìn)行詳細(xì)的注釋,并將數(shù)據(jù)被劃分為3類,即Classl：專家已經(jīng)作出了標(biāo)注；Class2：原始數(shù)據(jù)；Class3：處于研究進(jìn)展之中。因此,PhysioBank數(shù)據(jù)庫中的數(shù)據(jù)足已滿足生物醫(yī)學(xué)各領(lǐng)域研究者的需要。PhysioBank數(shù)據(jù)庫中的每一條數(shù)據(jù)記錄包括至少三類文件,頭文件〔.hea、數(shù)據(jù)文件〔.dat和注釋文件〔.atr,.al,.aiM等。頭文件是描述數(shù)據(jù)屬性的文本文件,其容包括記錄名、信號(hào)數(shù)目、貯存格式、信號(hào)數(shù)量和類型、采樣頻率、數(shù)字化特征、記錄的持續(xù)時(shí)間和起始時(shí)間等信息。一般可由PhysioToolkit軟件庫的WFDB庫函數(shù)的getinfo、putinfo函數(shù)讀和寫的字符。數(shù)據(jù)文件是定義了相應(yīng)存貯格式的數(shù)字化采樣點(diǎn)的二進(jìn)制存儲(chǔ)文件。數(shù)據(jù)存貯格式在頭文件中說明,一個(gè)數(shù)據(jù)組有相同的數(shù)據(jù)存貯格式,常用的是8位和16位格式。注釋文件是記錄了對(duì)信號(hào)特征的注釋信息[2]。ECG信號(hào)的特征典型心電信號(hào)波形心臟搏動(dòng)及其節(jié)律性是人體生命和生理狀態(tài)的重要標(biāo)志之一。心電生理學(xué)的研究表明,心電信號(hào)來源于心肌細(xì)胞的生物電變化。心肌細(xì)胞的電激動(dòng)稱為除極,心肌細(xì)胞恢復(fù)為靜息狀態(tài)稱為復(fù)極,心電信號(hào)的產(chǎn)生與心肌細(xì)胞的除極和復(fù)極過程密切相關(guān)。心臟電激動(dòng)起源于竇房結(jié),沿特化的心臟傳導(dǎo)系統(tǒng)下傳,其傳播方向、途徑、次序及時(shí)間存在一定的規(guī)律。若心臟不能及時(shí)發(fā)出電激動(dòng),則心臟陷于停博。人體體液中充滿電解質(zhì),具有導(dǎo)電性能,心臟電激動(dòng)過程產(chǎn)生的有序生物電變化通過體液傳至身體表面使身體各部位出現(xiàn)有規(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波：由左右心房的除極過程引起,其波形小而圓鈍,時(shí)寬為0.08s-0.11s,波幅不超過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間期：從Q波開始至S波結(jié)束之間的時(shí)程,反映心室除極時(shí)間。正常QRS間期為0.04s-0.1s,若QRS間期>0.12s,反映室傳導(dǎo)阻滯。PR間期：從P波開始到QRS波開始之間的時(shí)程,反映激動(dòng)由竇房結(jié)產(chǎn)生經(jīng)由結(jié)間束、房室交界和左右束支抵達(dá)心室,并引起心室興奮所需要的時(shí)間,又稱為房室傳導(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波開始到T波終點(diǎn)的時(shí)程,反映心室除極和復(fù)極時(shí)間的總和。許多因素可影響QT間期,如心肌缺血、低血鉀、低血鈣等可使QT間期延長(zhǎng),QT間期延長(zhǎng)使心室肌復(fù)極不均一,易誘發(fā)折返激動(dòng),導(dǎo)致嚴(yán)重室性心率失常。QT間期隨受心率變化的影響,心率越慢,QT間期越長(zhǎng)；心率越快,QT間期越短。通常用QTc間期修正心率對(duì)QT間期影響,正常QTc間期小于0.43~0.44s。ST段：指從QRS波群終止點(diǎn)到T波起點(diǎn)之間的波形線段,反映心室部分己完全進(jìn)入去極化狀態(tài),正常時(shí)與基線平齊。PP間期：相鄰P波之間的間距稱為PP間期,反映心房率。正常情況下,PP間期與RR間期一致。在Ⅱ度或Ⅱ度以上房室傳導(dǎo)阻滯和某些心率失常,兩者可不一致。RR間期：相鄰QRS波群之間的間距稱為RR間期,反映心室率。正常情況下,RR間期與PP間期一致[2]。在心電信號(hào)的測(cè)試中,對(duì)電極的放置部位和導(dǎo)聯(lián)的連接方式臨床有明確的規(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)〔Vl～V6,共有12個(gè)導(dǎo)聯(lián),具體可參考文獻(xiàn)[4]。<三>典型心電信號(hào)的能量〔頻譜圖分布典型的心電信號(hào)的整個(gè)心動(dòng)周期的頻譜估計(jì)圖如圖4.2[2]所示,可以明顯看出心電信號(hào)各波的能量主要集中在低頻區(qū)域,且隨著頻率的增高,相應(yīng)的能量逐漸降低。心電信號(hào)的整體頻譜圍在0.05Hz～100Hz,但能量主要集中在0.5～45Hz,能量的最高點(diǎn)在8～15Hz附近；QRS波群的頻譜帶寬為3～40Hz,積聚了將近99%的能量,波峰能量集中在6～18Hz附近,P波的頻譜帶寬為0～18Hz,波峰能量集中在5～12Hz；T波的頻譜帶寬為0～8Hz,波峰能量集中在0～8Hz區(qū)間[5]。圖4.2典型的心電信號(hào)頻譜能量分布<四>ECG信號(hào)的噪聲分析在采集、放大及傳輸心電信號(hào)的過程中,由于受人體、采集儀器、電磁環(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)成,幅度較低?！?基線漂移產(chǎn)生的原因主要由于人體呼吸運(yùn)動(dòng)、電極接觸不良等因素所導(dǎo)致。表現(xiàn)為心電信號(hào)上疊加緩慢變化的信號(hào),其頻率一般小于1Hz,幅度為ECG峰-峰的15%?！?肌電干擾產(chǎn)生的原因主要由于人體活動(dòng),肌肉緊所引起的干擾。表現(xiàn)為不規(guī)則的快速變化波形,其頻率圍較寬,一般在5～2kHz之間,幅度為毫伏級(jí)?！?運(yùn)動(dòng)偽跡產(chǎn)生的原因主要由于電極與人體間輕微移動(dòng)或抖動(dòng)而引入的干擾,表現(xiàn)為信號(hào)基線的短暫變化,但不是基線的躍變,其持續(xù)時(shí)間為100～500ms,頻率一般在7Hz以下,幅度較大?！?其他隨機(jī)噪聲心電信號(hào)還受到其他的隨機(jī)噪聲和環(huán)境干擾的影響,如加性白噪聲、極化噪聲、儀器部噪聲等。由此可見,噪聲信號(hào)基本覆蓋了有用的心電信號(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),在硬件方面來看,消除基線漂移的干擾應(yīng)考慮分別設(shè)計(jì)下限頻率為0.5Hz的高通濾波器；消除肌電等高頻干擾應(yīng)考慮設(shè)計(jì)上限頻率為100Hz的低通濾波器,同時(shí)還應(yīng)考慮設(shè)計(jì)60Hz的陷波器來濾除工頻干擾。由于硬件濾波器的元器件精度、穩(wěn)定性要求以及難以實(shí)現(xiàn)嚴(yán)格線性相位等問題,使得僅采用硬件濾波不能滿足濾波性能的要求。隨著數(shù)字信號(hào)處理技術(shù)的發(fā)展,設(shè)計(jì)高精度、高可靠性、簡(jiǎn)潔靈活的數(shù)字濾波器成為可能。目前,濾除心電信號(hào)干擾更多采用的是具有線性相位的數(shù)字濾波方法,并已逐漸顯示出取代硬件濾波器的趨勢(shì)?！?工頻干擾對(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)用中受到一定的限制。梳狀濾波器的名稱源于該濾波器的傳遞函數(shù)的幅頻特性形如梳狀。其特點(diǎn)是運(yùn)算量小、形式靈活簡(jiǎn)單,是一種快速的數(shù)字濾波器。梳狀濾波器被廣泛應(yīng)用去除工頻及其高次諧波的干擾。選擇合適的梳狀濾波參數(shù),可保證濾波器對(duì)有用信號(hào)基本無衰減和濾波后基本不產(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濾波方法中引入了判定線性段和非線性段的指標(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)象,濾波效果不能令人滿意[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]。一般來說,RLS自適應(yīng)算法與LMS自適應(yīng)算法相比,權(quán)值有更快的收斂速度。缺點(diǎn)是每次迭代計(jì)算量較大,由于在生物醫(yī)學(xué)應(yīng)用中數(shù)據(jù)的計(jì)算量不是太大,RLS自適應(yīng)算法可獲得更好性能。〔2基線漂移的消除 基線漂移可嚴(yán)重影響心電信號(hào)的分析與處理。濾波法和擬合基漂法是去除心電基線漂移的兩類主要算法。具體而言有形態(tài)濾波、中值濾波、線性相位FIR濾波、小波變換、形態(tài)和小波相結(jié)合、中值濾波和小波相結(jié)合等去除基漂的濾波方法,以及導(dǎo)數(shù)法與坐標(biāo)法這兩類基漂擬合方法。研究表明濾波法實(shí)現(xiàn)簡(jiǎn)單,但精度較低,分段擬合法去除基漂的效果更好,有自適應(yīng)的特性而且失真較小。無論是實(shí)時(shí)性還是準(zhǔn)確性,基漂擬合法都比濾波法更有優(yōu)勢(shì)但某些場(chǎng)合下,基漂擬合點(diǎn)的提取十分困難,此時(shí)只能使用濾波法[11]。對(duì)于基線擬合法,通常多采用簡(jiǎn)潔、快速的分段擬合法。其基本原理是通過對(duì)采樣的心電數(shù)據(jù)分段進(jìn)行三次多項(xiàng)式或拋物線擬合,獲取基線漂移曲線,然后用原始心電曲線減去擬合曲線來消除基線漂移。算法的關(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)滿足了信號(hào)分析。信號(hào)x<t>在不同尺度分解的逼近信號(hào)和細(xì)節(jié)信號(hào)具有不同頻帶,而在不同的尺度上一定包含有x<t>的不同頻率信息。所以若對(duì)包含有不同頻率信息的細(xì)節(jié)信號(hào)進(jìn)行閾值的取舍,丟棄其中的干擾成分后,再對(duì)信號(hào)進(jìn)行重構(gòu),則重構(gòu)信號(hào)的干擾得到了非常好的抑制,即為小波變換消噪的基本原理[2]。通過對(duì)集中方法的實(shí)驗(yàn)效果來選擇合適的消除基線漂移的方法。〔六R波峰位置、剔除異搏、確定RR間期〔1波峰位置找出波峰位置的基本想法是尋找一個(gè)采樣點(diǎn),其電壓值既大于前一個(gè)采樣點(diǎn)的值也大于后一個(gè)采樣點(diǎn)的值。之后再將P波與T波中的偽峰值去除,剩下的就是QRS波的峰值,即R波峰。可以通過先屏蔽掉小于某一電壓的所有采樣點(diǎn),這樣就將P波和T波中的偽峰值提前去除,之后再尋找峰值。以此保證得到的峰值處于QRS波中,該峰值即為R波峰值。剔除異搏對(duì)得到R波的位置,檢查其中的異常點(diǎn),去掉異搏。確定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：全程相鄰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)。另一角度,即頻譜分析的角度來分析心率變化的規(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大小,LF/HF比值反映自主神經(jīng)系統(tǒng)的平衡狀態(tài),基本上代表經(jīng)交感神力的高低。步驟〔包括計(jì)算程序流程框圖總程序流程圖。詳細(xì)程序間附件1。獲取心電信號(hào)獲取心電信號(hào)去除工頻干擾和基線漂移去除工頻干擾和基線漂移確定波峰位置確定波峰位置剔除異搏剔除異搏確定R-R間期確定R-R間期求出HRV信號(hào)求出HRV信號(hào)進(jìn)行時(shí)域計(jì)算分析進(jìn)行線性內(nèi)插進(jìn)行時(shí)域計(jì)算分析進(jìn)行線性內(nèi)插進(jìn)行頻域計(jì)算分析進(jìn)行功率譜分析進(jìn)行頻域計(jì)算分析進(jìn)行功率譜分析HRV信號(hào)NNVGRSDNNHRV信號(hào)NNVGRSDNNRMSSDSDSDpNN50頻譜圖TP頻譜圖TPHFLFLF/HF結(jié)論結(jié)論獲取心電信號(hào)將下載得到三組心電信號(hào)：100、101、102,通過KlausRheinberger編寫的讀取心電的MATLAB程序ecg.load.m,加載原始心電信號(hào),從而準(zhǔn)確提取相應(yīng)的心電信號(hào)。去除工頻噪聲干擾對(duì)信號(hào)進(jìn)行頻域分析,發(fā)現(xiàn)信號(hào)主要集中在低頻端,并且在60Hz處有微弱信號(hào),通過判斷可得出該頻率信號(hào)為工頻噪聲。本文使用FIR濾波器將噪聲濾去。去基線漂移對(duì)比兩大類型發(fā)現(xiàn)：濾波法實(shí)現(xiàn)簡(jiǎn)單,但精度較低；分段擬合法去除基漂的效果更好,有自適應(yīng)的特性而且失真較小。無論是實(shí)時(shí)性還是準(zhǔn)確性,基漂擬合法都比濾波法更有優(yōu)勢(shì)；但某些場(chǎng)合下,基漂擬合點(diǎn)的提取十分困難,此時(shí)只能使用濾波法。通過進(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í)域分析,其中的技術(shù)指標(biāo)有NNVGR、SDNN、RNSSD、SDSD、pNN50等。<七>線性插對(duì)提取出來的HRV信號(hào)經(jīng)線性插獲得等間隔〔R-R間期均值R-R間期時(shí)間序列,然后再進(jìn)行頻譜分析。線性插的公式為。式中X<n>為插值,和分別為插值前、后的R-R間期序列值,和分別為X〔n所在位置與和所在位置的時(shí)間間距?！舶烁鶕?jù)快速傅里葉變換〔FFT得到HRV信號(hào)的頻譜信號(hào)和功率譜信號(hào)。主要有TP、HF、LF、LF/HF、頻譜圖、功率譜圖等。結(jié)果通過MATLAB程序分別得到100、101和102三組ECG信號(hào)數(shù)據(jù),然后對(duì)其進(jìn)行相關(guān)分析,具體結(jié)果以及圖像如下所述?！惨籈CG信號(hào)的提取通過MATLAB程序分別得到100、101和102三組ECG信號(hào)數(shù)據(jù)。其時(shí)序圖像具體如圖6.1所示。圖6.1三組原始心電信號(hào)圖示由上至下分別為來自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像〔二原始心電頻譜圖對(duì)三組信號(hào)通過皮譜圖進(jìn)行頻域分析。從圖6.2中可以看出,信號(hào)主要集中在20Hz以下的低頻端,并且在60Hz附近有一個(gè)明顯的幅度值上升,可以判斷其為工頻干擾。圖6.2三組心電信號(hào)的頻譜圖圖示由上至下分別為來自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像去除工頻干擾使用FIR濾波器法將工頻噪聲去除。去除工頻干擾之后頻譜圖如圖6.3所示圖6.3三組心電信號(hào)去除工頻干擾后的頻譜圖圖示由上至下分別為來自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)的頻譜圖圖示由上至下分別為來自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像確定波峰確定波峰結(jié)果如圖6.5所示。其中紅色點(diǎn)表示波峰。圖6.5確定三組心電信號(hào)的R波峰圖示由上至下分別為來自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)的HRV圖像圖示由上至下分別為來自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>NN50pNN50100.811.434.247.539.976.93%101.859.348.828.617.588.42%102.818.4143.7199.2187.51111.34%〔八頻域分析根據(jù)快速傅里葉變換〔FFT得到HRV信號(hào)的頻譜信號(hào)和功率譜信號(hào)。頻譜圖。三組心電數(shù)據(jù)對(duì)應(yīng)的HRV信號(hào)的頻譜圖如圖6.7所示。圖6.7消除基線后三組心電信號(hào)的頻譜圖圖示由上至下分別為來自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像功率譜圖圖6.8HRV信號(hào)的功率譜圖圖示由上至下分別為來自100號(hào)、101號(hào)、102號(hào)數(shù)據(jù)的對(duì)應(yīng)圖像根據(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.dat1900014005710.4150101.dat2060087111001.2063102.dat28300640037000.5812結(jié)論本文通過使用MATLAB軟件,實(shí)現(xiàn)了心電信號(hào)的讀取、采樣及頻譜分析,并且分離出心率變異性HRV信號(hào),通過時(shí)域分析和頻域分析的方法進(jìn)行分析,對(duì)比資料發(fā)現(xiàn)：第101號(hào)數(shù)據(jù)組應(yīng)該來自正常人,第100號(hào)數(shù)據(jù)組應(yīng)該來自具有輕微患心血管疾病的患者,第102號(hào)數(shù)據(jù)組則最可能來自患有較為嚴(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ì)多種惡性心律失常的預(yù)后判斷和藥物治療效果分析有指導(dǎo)作用。所以,對(duì)HRV的研究能夠極大的促進(jìn)人類對(duì)于心血管疾病的了解,從而在預(yù)防、治療心血管疾病等領(lǐng)域取得成果。本文從獲取的心電信號(hào)中,提取出有效的心電信號(hào),進(jìn)而提取出我們所研究的HRV信號(hào),然后對(duì)其進(jìn)行相關(guān)分析。在獲取心電信號(hào)的過程中,直接使用了來自美國PhysioNet的心電信號(hào)資源,在提取出有效心電信號(hào)的過程中,首先進(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間期；最后,得到我們所需要的HRV信號(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)過對(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)行了線性插。以上種種方法都是為了使數(shù)據(jù)具有更高的準(zhǔn)確性,基本確保研究沒有在研究對(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é)院,20XX.[2]董紅生.心電波形檢測(cè)與心率異變性分析研究:[博士論文].理工大學(xué),20XX[3]吳學(xué)勤.動(dòng)態(tài)心電圖技術(shù)與應(yīng)用.:中國科學(xué)技術(shù),1998,22-25.[4]胡大一,郭成軍,瑞杰.心率變異性一測(cè)量標(biāo)準(zhǔn),生理釋義與臨床應(yīng)用—續(xù)一:心率變異性的生理研究.中國醫(yī)療器械信息,1998,4<1>:29-32.[5]細(xì)武,鄧親愷.心電圖各波的頻率分析.中國醫(yī)學(xué)物理學(xué)雜志,2001,18<1>:46-48.[6]萬相奎.心電信號(hào)分析與虛擬式心電自動(dòng)分析儀的開發(fā):[大學(xué)博士學(xué)位論文].:大學(xué)機(jī)械電子工程,2005,16-17.[7]LevkovC,MichovG,IvanovRetal.Subtractionof50Hzinterferencefromtheelectrocardiogram.Med.&Biol.,Eng.&Comp.,1984,22<4>:371-373.[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é)工程.20XX第22卷第五期：381-386.附件1：%ThisprogrammreadsECGdatawhicharesavedinformat212.%<e.g.,100.datfromMIT-BIH-DB,cu01.datfromCU-DB,...>%Thedataaredisplayedinafiguretogetherwiththeannotations.%TheannotationsaresavedinthevectorANNOT,thecorresponding%times<inseconds>aresavedinthevectorATRTIME.%Theannotationsaresavedasnumbers,themeaningofthenumberscan%befoundinthecodetable"ecgcodes.h"availableat..%%ANNOTonlycontainsthemostimportantinformation,whichisdisplayed%withtheprogramrdann<availableon.>inthe3rdrow.%The4thto6throwarenotsavedinANNOT.%%%createdonFeb.27,2003by%RobertTratnig<VorarlbergUniversityofAppliedSciences>%<email:rtratniggmx.at>,%%algorithmisbasedonaprogramwrittenby%KlausRheinberger<UniversityofInnsbruck>%<email:klaus.rheinbergeruibk.ac.at>%%clc;clearall;%SPECIFYDATAPATH='C:\Users\asus你\Desktop\心率變異性〔HRV信號(hào)的提取及時(shí)頻域分析\ecg';%path,wheredataaresavedHEADERFILE='100.hea';%header-fileintextformatATRFILE='100.atr';%attributes-fileinbinaryformatDATAFILE='100.dat';%data-fileSAMPLES2READ=30000;%numberofsamplestoberead%incaseofmorethanonesignal:%2*SAMPLES2READsamplesareread%LOADHEADERDATAfprintf<1,'\\n$>WORKINGON%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>;%numberofsignalssfreq=A<2>;%samplerateofdataclearA;fork=1:nosigz=fgetl<fid1>;A=sscanf<z,'%*s%d%d%d%d%d',[1,5]>;dformat<k>=A<1>;%format;hereonly212isallowedgain<k>=A<2>;%numberofintegerspermVbitres<k>=A<3>;%bitresolutionzerovalue<k>=A<4>;%integervalueofECGzeropointfirstvalue<k>=A<5>;%firstintegervalueofsignal<totestforerrors>end;fclose<fid1>;clearA;%LOADBINARYDATAifdformat~=[212,212],error<'thisscriptdoesnotapplybinaryformatsdifferentto212.'>;end;signald=fullfile<PATH,DATAFILE>;%datainformat212fid2=fopen<signald,'r'>;A=fread<fid2,[3,SAMPLES2READ],'uint8'>';%matrixwith3rows,each8bitslong,=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<M2H,8>+A<:,3>-PRR;ifM<1,:>~=firstvalue,error<'inconsistencyinthefirstbitvalues'>;end;switchnosigcase2M<:,1>=<M<:,1>-zerovalue<1>>/gain<1>;M<:,2>=<M<:,2>-zerovalue<2>>/gain<2>;TIME=<0:<SAMPLES2READ-1>>/sfreq;case1M<:,1>=<M<:,1>-zerovalue<1>>;M<:,2>=<M<:,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%thiscasedidnotappearuptonow!%hereMhastobesorted!!!disp<'Sortingalgorithmformorethan2signalsnotprogrammedyet!'>;end;clearAM1HM2HPRRPRL;fprintf<1,'\\n$>LOADINGDATAFINISHED\n'>;%LOADATTRIBUTESDATAatrd=fullfile<PATH,ATRFILE>;%attributefilewithannotationdatafid3=fopen<atrd,'r'>;A=fread<fid3,[2,inf],'uint8'>';fclose<fid3>;ATRTIME=[];ANNOT=[];sa=size<A>;saa=sa<1>;i=1;whilei<=saaannoth=bitshift<A<i,2>,-2>;ifannoth==59ANNOT=[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;elseifannoth==60%nothingtodo!elseifannoth==61%nothingtodo!elseifannoth==62%nothingtodo!elseifannoth==63hilfe=bitshift<bitand<A<i,2>,3>,8>+A<i,1>;hilfe=hilfe+mod<hilfe,2>;i=i+hilfe/2;elseATRTIME=[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>>=[];%lastline=EOF<=0>ATRTIME<length<ATRTIME>>=[];%lastline=EOFclearA;ATRTIME=<cumsum<ATRTIME>>/sfreq;ind=find<ATRTIME<=TIME<end>>;ATRTIMED=ATRTIME<ind>;ANNOT=round<ANNOT>;ANNOTD=ANNOT<ind>;%DISPLAYDATAfigure<1>;clf,boxon,holdonplot<TIME,M<:,1>,'r'>;%修改,'LineWidth',2DATA1=M<:,1>';%ifnosig==2%plot<TIME,M<:,2>,'b'>;%end;%fork=1:length<ATRTIMED>%text<ATRTIMED<k>,0,num2str<ANNOTD<k>>>;%end;xlim<[TIME<1>,TIME<end>]>;xlabel<'Time/s'>;ylabel<'Voltage/mV'>;string=['ECGsignal',DATAFILE];title<string>;fprintf<1,'\\n$>DISPLAYINGDATAFINISHED\n'>;%fprintf<1,'\\n$>ALLFINISHED\n'>;%ecg頻域分析fs=length<TIME>/TIME<end>;[f1,y1]=fft_simple<DATA1,fs,0,180>;figure<2>;plot<f1,y1,'r','LineWidth',2>;xlabel<'Frequency/Hz'>;ylabel<'Amplitude'>;title<'ecg頻譜圖'>;%去除工頻噪聲f_p=[4080];f_s=[5961];R_p=3;R_s=20;Ws=f_s/<fs/2>;Wp=f_p/<fs/2>;[n,Wn]=buttord<Wp,Ws,R_p,R_s>;[b,a]=butter<n,Wn,'stop'>;DATA=filter<b,a,DATA1>;[f2,y2]=fft_simple<DATA,fs,0,180>;figure<3>;plot<f2,y2,'r','LineWidth',2>;xlabel<'Frequency/Hz'>;ylabel<'Amplitude'>;title<'去除工頻干擾后頻譜圖'>;%多項(xiàng)式擬合%p=polyfit<TIME,DATA,100>;%z=polyval<p,TIME>;%figure<4>;plot<TIME,e>;%消除基線漂移z=medfilt1<DATA,330>;%中值濾波法DATA=DATA-z;figure<4>;plot<TIME,DATA,'r'>;%修改,'LineWidth',2xlim<[TIME<1>,TIME<end>]>;xlabel<'Time/s'>;ylabel<'Voltage/mV'>;title<'消除基線后ecg時(shí)序圖'>;%確定波峰fori=1:length<TIME>ifDATA<1,i><0.8DATA<1,i>=0;endendfigure<5>;plot<TIME,DATA,'--o'>;xlim<[080]>;holdon;count=0;forj=2:length<TIME>-1ifDATA<1,j>~=0if<DATA<1,j>>=DATA<1,j-1>>&&<DATA<1,j>>DATA<1,j+1>>count=count+1;a<1,count>=TIME<1,j>;a<2,count>=DATA<1,j>;endendendfigure<5>;plot<a<1,:>,a<2,:>,'*r'>;xlim<[080]>;ylim<[02]>;title<'確定波峰位置'>;%確定R-R間距b=diff<a<1,:>>;figure<6>;plot<b,'r','LineWidth',3>;xlim<[1length<b>]>;ylim<[01.2]>;gridon;xlabel<'R-Ri'>;ylabel<'R-