1、課程設計報告題目：心率變異性（HRV信號的提取及時頻域分析專業：生物醫學工程班級：XXXXXXX學號：XXXXXXX姓名：XXXXXXX指導教師：_XXXXXXXXXXXXXX學XXXXX學院2016年9月29日一、開題背景（一）HR播介傳統的醫學觀點認為，正常的心率為規則的竇性節律；后來發現在健康狀態下，許多生理系統中存在自然的變異性，人的心率正常情況下也是呈不規則性變化的，而心率變異就是指竇性心率的這種波動變化的程度。心率變異性（HeartRateVariability,HRV是指逐次心搏間期之間的微小變異特性。在生理條件下，HRV勺產生主要是由于心臟竇房結自律活動通過交感和迷走神經，神經

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

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

4、，大約6-12個月恢復正常。因此，多次測定HRVT能比單次測定價值更大。梗塞后HRVS復的快慢對以后死亡的危險性也有預測價值。（3）對糖尿病患者自主神經系統損傷的評估：糖尿病患者不論病情輕重，均存在不同程度的自主神經功能紊亂。HRM判斷糖尿病患者是否伴有自主神經系統損害最準確，最敏感的指標。（4）心力衰竭（CHF患者危險性評估（5）心率變異性生物反饋療法：對于不孕人群受孕幾率提高、懷孕人群孕期焦慮癥改善、產后人群產后抑郁癥情況緩解，起到很好的作用。（6）其它臨床應用范圍：心絞痛、高血壓、心肌病、非缺心臟病所致的慢性嚴重二尖瓣返流、二尖瓣脫垂、心律失常、血管迷走性暈厥等心血管疾病。二、課題目的（

5、一）基本掌握心電信號（ECG的測量、數據采集的方法。（二）學會使用MATLABtECG®號進行相關處理分析。主要包括從ECGW號中提取出所需的HRVJ號，并分別對其進行時域、頻域、功率譜上的分析。（三）掌握HRVJ號的時頻域參數的意義，以及對其進行分析的基本方法。三、課題研究的主要內容（一）從網上下載正常人的心電信號以及各種病人的心電信號（ECG數據。（二）首先HRVIW號的提取，主要包括去除干擾、準確確定R波波峰位置、剔除異搏、確定R-R間期、線性內插，并且繪出HRV言號曲線。（三）對HRV言號的時域分析，對HRVIW號的頻譜圖和功率譜圖分析。四、原理和方法（一）ECG言號的采集本

6、文主要使用100.hea、100dat、100.atr,101.hea、101.dat、101.atr,102.hea、102.dat、102.atr這三組數據來對HRV!行研究。實驗數據來源于PhysioNet。PhysioNet是一個基于Web勺復雜生理和生物醫學信號的研究資源網站，其網址為。PhysioNet由PhysioNet,PhysioBank和PhysioToolkit三個相互關聯的部分組成。數據庫中數據來源于正常人、各種病人（如心臟猝死、心力衰竭、心律失常、癲癇、睡眠呼吸暫停綜合癥等）及運動、休息等不同狀態下的數據，樣本選取范圍廣泛

7、，其中大部分數據都進行詳細的注釋，并將數據被劃分為3類，即Classl:專家已經作出了標注；Class2：原始數據；Class3:處于研究進展之中。因此，PhysioBank數據庫中的數據足已滿足生物醫學各領域研究者的需要。PhysioBank數據庫中的每一條數據記錄包括至少三類文件，頭文件（.hea）、數據文件（.dat）和注釋文件（.atr,.al,.aiM等）。頭文件是描述數據屬性的文本文件，其內容包括記錄名、信號數目、貯存格式、信號數量和類型、采樣頻率、數字化特征、記錄的持續時間和起始時間等信息。一般可由PhysioToolkit軟件庫的WFD酢函數的getinfo、putinfo函數

8、讀和寫的字符。數據文件是定義了相應存貯格式的數字化采樣點的二進制存儲文件。數據存貯格式在頭文件中說明，一個數據組有相同的數據存貯格式，常用的是8位和16位格式。注釋文件是記錄了對信號特征的注釋信息2。（二）ECG言號的特征(1)典型心電信號波形心臟搏動及其節律性是人體生命和生理狀態的重要標志之一。心電生理學的研究表明，心電信號來源于心肌細胞的生物電變化。心肌細胞的電激動稱為除極，心肌細胞恢復為靜息狀態稱為復極，心電信號的產生與心肌細胞的除極和復極過程密切相關。心臟電激動起源于竇房結，沿特化的心臟傳導系統下傳，其傳播方向、途徑、次序及時間存在一定的規律。若心臟不能及時發出電激動，則心臟陷于停博。

9、人體體液中充滿電解質，具有導電性能，心臟電激動過程產生的有序生物電變化通過體液傳至身體表面使身體各部位出現有規律而各向異性的電位變化，通過測量電極采集體表特定點電位變化，并放大、顯示及記錄，即為體表心電信號，也即是通常的ECG信號3。一個心動周期正常心電信號波形如圖4,1所示。它是由特征波及其特征間期組成，每個心動周期包含一個P波，一個QRS波群和一個T波，有時還會出現一個小的U波。特征波及特征間期的含義如下：圖4.1典型的心電信號波形P波：由左右心房的除極過程引起，其波形小而圓鈍，時寬為0.08s-0.11s,波幅不超過0.25mVoQRS波：反映左右心室除極產生的電位變化，在P波之后出現，

10、為心電信號中最高大和最快速的波形。典型的QRS波群包括三個相連的波，第一個向下的波為Q波，緊接著為高而尖峭的向上的R波，最后是一個向下的S波。在體表不同位置(使用不同導聯記錄)時，三個波不一定都有，大小方向也會不同。T波：代表心室復極時的電位變化，方向與R波方向相同，時寬為0.05s-0.25s,波幅一般為0.lmV-0.8mV。U波：T波之后可能出現的一個低而寬的波，與T波方向一致，其機理不十分清楚，可能反映普頃野纖維復極的電位變化。QRS間期：從Q波開始至S波結束之間的時程，反映心室除極時間。正常QRS問期為0.04s-0.1s,若QRS問期0,12s,反映室內傳導阻滯。PR間期：從P波開

11、始到QRS波開始之間的時程，反映激動由竇房結產生經由結問束、房室交界和左右束支抵達心室，并引起心室興奮所需要的時間，又稱為房室傳導時問。正常為0.12s-0.2s。當發生房室傳導阻滯時，PR間期增長。如當PR問期0.21s。則為I度房室傳導阻滯。QT間期：從QRS波開始到T波終點的時程，反映心室除極和復極時間的總和。許多因素可影響QT間期，如心肌缺血、低血鉀、低血鈣等可使QT間期延長，QT問期延長使心室肌復極不均一，易誘發折返激動，導致嚴重室性心率失常。QT間期隨受心率變化的影響，心率越慢，QT間期越長；心率越快，QT間期越短。通常用QTc間期修正心率對QT問期影響，正常QTc問期小于0.43

12、0.44s。ST段：指從QRS波群終止點到T波起點之間的波形線段，反映心室部分己完全進入去極化狀態，正常時與基線平齊。PP間期：相鄰P波之間的間距稱為PP間期，反映心房率。正常情況下，PP問期與RR間期一致。在II度或II度以上房室傳導阻滯和某些心率失常，兩者可不一致。RR間期：相鄰QRS波群之間的間距稱為RR間期，反映心室率。正常情況下，RR間期與PP間期一致。在心電信號的測試中，對電極的放置部位和導聯的連接方式臨床有明確的規定。目前，國際公認的是標準12導聯，包括心電標準導聯（I、II、III）、加壓單極肢體導聯（aVRaVL、aVF）及胸導聯（VlV6）,共有12個導聯，具體可參考文獻4

13、。（三）典型心電信號的能量（頻譜圖）分布典型的心電信號的整個心動周期的頻譜估計圖如圖4.22所示，可以明顯看出心電信號各波的能量主要集中在低頻區域，且隨著頻率的增高，相應的能量逐漸降低。心電信號的整體頻譜范圍在0.05Hz100Hz,但能量主要集中在0.545Hz,能量的最高點在815Hz附近；QRS波群的頻譜帶寬為340Hz,積聚了將近99%勺能量，波峰能量集中在618Hz附近，P波的頻譜帶寬為018Hz,波峰能量集中在512Hz;T波的頻譜帶寬為08Hz,波峰能量集中在08Hz區間。1.0圖4.2典型的心電信號頻譜能量分布(四)ECG言號的噪聲分析在采集、放大及傳輸心電信號的過程中，由于受

14、人體、采集儀器、電磁環境、操作水平等的影響，不可避免會有許多干擾耦合到心電信號，主要干擾表現形式如下：(1)電源工頻干擾產生的原因主要由于電源磁場作用于心電圖儀的導聯與人體之間的環形電路所致，表現為心電信號上有明顯的正弦波或正弦波的疊加信號，其頻率為60Hz工頻及其諧波構成，幅度較低。(2)基線漂移產生的原因主要由于人體呼吸運動、電極接觸不良等因素所導致。表現為心電信號上疊加緩慢變化的信號，其頻率一般小于1Hz,幅度為ECG1-峰的15%(3)肌電干擾產生的原因主要由于人體活動，肌肉緊張所引起的干擾。表現為不規則的快速變化波形，其頻率范圍較寬，一般在52kHz之間，幅度為毫伏級。(4)運動偽跡

15、產生的原因主要由于電極與人體間輕微移動或抖動而引入的干擾，表現為信號基線的短暫變化，但不是基線的躍變，其持續時間為100500m6頻率一般在7Hz以下，幅度較大。(5)其他隨機噪聲心電信號還受到其他的隨機噪聲和環境干擾的影響，如加性白噪聲、極化噪聲、儀器內部噪聲等。由此可見，噪聲信號基本覆蓋了有用的心電信號的全頻率范圍。而其中60Hz及其倍頻附近的工頻干擾、1Hz以下的基線漂移以及肌電干擾噪聲是最主要的干擾源，在心電信號預處理中必須消除或抑制，以提高心電信號的信噪比2（五）ECG言號去噪聲方法對于心電信號預處理一般要從硬件電路優化設計和軟件數字濾波器的設計兩個方面考慮。根據心電信號的頻譜分布特

16、點，在硬件方面來看，消除基線漂移的干擾應考慮分別設計下限頻率為0.5Hz的高通濾波器；消除肌電等高頻干擾應考慮設計上限頻率為100Hz的低通濾波器，同時還應考慮設計60Hz的陷波器來濾除工頻干擾。由于硬件濾波器的元器件精度、穩定性要求以及難以實現嚴格線性相位等問題，使得僅采用硬件濾波不能滿足濾波性能的要求。隨著數字信號處理技術的發展，設計高精度、高可靠性、簡潔靈活的數字濾波器成為可能。目前，濾除心電信號干擾更多采用的是具有線性相位的數字濾波方法，并已逐漸顯示出取代硬件濾波器的趨勢。（1）工頻干擾對于工頻干擾，主要的消除方法有平均濾波器、梳狀濾波器、Levkov濾波器、自適應工頻陷波器等方法。平

17、均濾波器法是較早被應用的數字濾波方法，其特點是算法簡單，處理速度快，濾波效果較好。由于平均濾波算法的實時性好，被廣泛應用于實時心電監護設備中或基于單片機心電數據采集系統中。從實際濾波效果看，濾波器能較好地濾除了60Hz工頻，但對ECG中的QRS波有較大削峰，信號衰減較大。因此，在實際應用中受到一定的限制。梳狀濾波器的名稱源于該濾波器的傳遞函數的幅頻特性形如梳狀。其特點是運算量小、形式靈活簡單，是一種快速的數字濾波器。梳狀濾波器被廣泛應用去除工頻及其高次諧波的干擾。選擇合適的梳狀濾波參數，可保證濾波器對有用信號基本無衰減和濾波后基本不產生相移6。改進的梳狀濾波器可用于消除ECG的工頻干擾且能有效

18、地防止信號的失真。Levkov濾波器是Levkov在1984年提出一種對ECG信號的線性段和非線性段采用不同處理方式的數字濾波方法。其基本原理是從原始ECG信號中直接減去在該線性段中確定的干擾信號的幅值7o改進的Levkov濾波方法中引入了判定線性段和非線性段的指標Mo算法的濾波效果與選擇合適的M值有關。改進Levkov濾波法算法簡單，參數可調，運算量小，可實現實時心電信號濾波處理，但在噪聲較強時，將會造成QRS波群的削峰現象，濾波效果不能令人滿意網。自適應工頻陷波技術特點是只需要很少或根本不需要任何有關信號和噪聲的先驗知識的情況下，直接利用觀測數據不斷遞歸更新自身的參數，以逐步逼近某一最優值

19、，能自動跟蹤信號頻率漂移，具有非常好的適應性陰。自適應的數字陷波器的權值調節主要有最小均方誤差（LMS、最小遞推二乘法（RLS）算法1°。一般來說，RLS自適應算法與LMS自適應算法相比，權值有更快的收斂速度。缺點是每次迭代計算量較大，由于在生物醫學應用中數據的計算量不是太大，RLS自適應算法可獲得更好性能。（2）基線漂移的消除基線漂移可嚴重影響心電信號的分析與處理。濾波法和擬合基漂法是去除心電基線漂移的兩類主要算法。具體而言有形態濾波、中值濾波、線性相位FIR濾波、小波變換、形態和小波相結合、中值濾波和小波相結合等去除基漂的濾波方法，以及導數法與坐標法這兩類基漂擬合方法。研究表明濾

20、波法實現簡單，但精度較低，分段擬合法去除基漂的效果更好，有自適應的特性而且失真較小。無論是實時性還是準確性，基漂擬合法都比濾波法更有優勢但某些場合下，基漂擬合點的提取十分困難，此時只能使用濾波法11。對于基線擬合法，通常多采用簡潔、快速的分段擬合法。其基本原理是通過對采樣的心電數據分段進行三次多項式或拋物線擬合，獲取基線漂移曲線，然后用原始心電曲線減去擬合曲線來消除基線漂移。算法的關鍵是插值點選擇，插值點應該是心電信號的零電位點，通常采用TP段。對于FIR濾波器只有零點，不易獲得較好的通帶與阻帶衰減特性。要取得好的衰減特性，一般要求H(z)的階次要高。對于小波變換，由于其提供了一個在時頻平面上

21、可調的分析窗口。信號的時、頻分辨率可以隨分析任務的需要作出相應調整。在低頻時小波變換的時間分辨率較差，而頻率分辨率較高；在高頻時小波變換的時間分辨率較高，而頻率分辨率較低，自動滿足了信號分析。信號x(t)在不同尺度分解的逼近信號和細節信號具有不同頻帶，而在不同的尺度上一定包含有x(t)的不同頻率信息。所以若對包含有不同頻率信息的細節信號進行閾值的取舍，丟棄其中的干擾成分后，再對信號進行重構，則重構信號的干擾得到了非常好的抑制，即為小波變換消噪的基本原理。通過對集中方法的實驗效果來選擇合適的消除基線漂移的方法。(六)R波峰位置、剔除異搏、確定RR間期(1)波峰位置找出波峰位置的基本想法是尋找一個

22、采樣點，其電壓值既大于前一個采樣點的值也大于后一個采樣點的值。之后再將P波與T波中的偽峰值去除，剩下的就是QR徽的峰值，即R波峰。可以通過先屏蔽掉小于某一電壓的所有采樣點，這樣就將P波和T波中的偽峰值提前去除，之后再尋找峰值。以此保證得到的峰值處于QR皺中，該峰值即為R波峰值。(2)剔除異搏對得到R波的位置，檢查其中的異常點，去掉異搏。(3)確定RR間期用后一個峰值的采樣時間減去前一個峰值的采樣時間就可以得到RR間期的值。將信號以RR值為縱坐標，R-Ri值為橫坐標，得到HRV勺圖像。(七)hrV言號的時域分析通常使用的HRV寸域檢測指標有5項：NNVGRSDNNRMSSDSDStJ口pNN50

23、上述5個指標的定義分別為：1NNNVGR全部正常竇性心博間期(NN的平均值，單位為ms即NNVGRRR。NidSDNN標準差，即全部NN問期的標準差，單位為ms即SDNN=口g（RR_職2。一NiJ7RNSSD全程相令口NN間期之差的均方根值，單位為ms即RNSSDN_1N1L(RRi_RR)i1N1SDSD全程相鄰NN|i期長度N差的標準差，單包為ms即SDSD=1'（RR'_R2一N_1iJ一2其中RR=RR-RR1,RR=（RR_RR手）。pNN50在全部NN間期的記錄中，相鄰的NN間期之差大于50ms的個數與總的NN間期的個數的比，以百分比表示。其中NNVGRI于評估心

24、率總體變化水平,;SDNNW于評估心率總體變化的大小，即交感及迷走神經張力大小，；SDANW于評估心率變化中的長期慢變化成分，即交感神經張力大小。RMSSDtpNN50反映心率快變化成分的大小，即副交感神經張力的敏感指標。（八）HRV信號的頻域分析：HRV寸域分析的指標大多用于描述HRV6體的大小，不能仔細地分析交感神經和迷走神經各自的活動的情況。而HRV»譜分析則可以彌補這個缺點。另一角度，即頻譜分析的角度來分析心率變化的規律。它與時域分析既有相關性，又能揭示出心率的更復雜的變化規律。FFT是經典譜估計方法，算法簡單。輸入和輸出信號能量有線性關系，但對信號要作周期延拓假定，短數據譜

25、分辨率較低，并有能量泄露現象。最近更仔細的研究發現，正常人基礎狀態下心率譜曲線在0-0.4HZ之間，0.003-0.04HZ為極低頻段（VLF）,0.04-0.15HZ為低頻段（LF）,0.15-0.4HZ高頻段（HF,0-0.40HZ為總功率譜（TP）。研究表明，VLF反映心率變化受熱調節（體溫），血管舒縮張力和腎血管緊張素系統的影響；LF反映交感和迷走神經的雙重調節；HF只反映迷走神經的調節；TP反映HRS小，LF/HF比值反映自主神經系統的平衡狀態，基本上代表經張交感神力的高低。五、步驟（包括計算程序流程框圖）（一）總程序流程圖。詳細程序間附件1（二）獲取心電信號將下載得到三組心電信號：

26、100、101、102,通過KlausRheinberger編寫的讀取心電的MATLAB?序ecg.load.m,加載原始心電信號，從而準確提取相應的心電信號。（三）去除工頻噪聲干擾對信號進行頻域分析，發現信號主要集中在低頻端，并且在60Hz處有微弱信號，通過判斷可得出該頻率信號為工頻噪聲。本文使用FIR濾波器將噪聲濾去。（四）去基線漂移對比兩大類型發現：濾波法實現簡單，但精度較低；分段擬合法去除基漂的效果更好，有自適應的特性而且失真較小。無論是實時性還是準確性，基漂擬合法都比濾波法更有優勢；但某些場合下，基漂擬合點的提取十分困難，此時只能使用濾波法。通過進行初步實驗中比較各種方法的實際效果，

27、發現使用中值濾波法的效果最好，故本文主要采用其進行去基線漂移處理。（五）找出波峰、剔除異搏、確定RR可期（六）XtHR%號進行時域分析對前面得到的HR%號進行時域分析，其中的技術指標有NNVGRSDNNRNSSDSDSDpNN5第。（七）線性內插對提取出來的HRV言號經線性內插獲得等間隔（R-R間期均值）R-R間期時間序列，然后再進行頻譜分析。線性內插的公式為x（n）=d2RRi+dlRRf。式中X（n）為d1d2d1d2插值，RR和RR4分別為插值前、后的R-R間期序列值，d1和d2分別為X（n）所在位置與RRi和RRi書所在位置的時間問距。（八）根據快速傅里葉變換（FFD得到HRVJ號的頻

28、譜信號和功率譜信號。主要有TRHF、LF、LF/HF、頻譜圖、功率譜圖等。六、結果通過MATLA毓序分別得至IJ100、101和102三組ECG®號數據，然后對其進行相關分析，具體結果以及圖像如下所述。（一）ECG言號的提取通過MATLA毓序分別得至IJ100、101和102三組ECG®號數據。其時序圖像具體如圖6.1所示。Ti力立/信圖6.1三組原始心電信號圖示由上至下分別為來自100號、101號、102號數據的對應圖像（二）原始心電頻譜圖對三組信號通過皮譜圖進行頻域分析。從圖6.2中可以看出，信號主要集中在20Hz以下的低頻端，并且在60Hz附近有一個明顯的幅度值上升，

29、可以判斷其為工頻干擾-liO<b-IAOa口nrhaoloona-n圖6.2三組心電信號的頻譜圖圖示由上至下分別為來自100號、101號、102號數據的對應圖像（三）去除工頻干擾使用FIR濾波器法將工頻噪聲去除。去除工頻干擾之后頻譜圖如圖6.3所示543cMdmo-o-去除工壯心抗潛圖Q204060SOIDO12Q1401601SOFrequency/Hz圖6.3三組心電信號去除工頻干擾后的頻譜圖圖示由上至下分別為來自100號、101號、102號數據的對應圖像（四）去除基線漂移從圖6.1中可以看出這些心電信號均有明顯的基線漂移。使用中值濾波法去除基線漂移。消除基線漂移后的ECGJ號時序圖

30、如圖6.4所示«1謂隈睛很精如中叫211H舞十且看I-it-m圖6.4消除基線后三組心電信號的頻譜圖圖示由上至下分別為來自100號、101號、102號數據的對應圖像（五）確定波峰確定波峰結果如圖6.5所示。其中紅色點表示波峰。i_r2-r.>JtaC！i1IlII圖6.5確定三組心電信號的R波峰圖示由上至下分別為來自100號、101號、102號數據的對應圖像（六）確定R-R間期R-R間期與R-Ri的關系，亦即HRV勺圖像如圖6.6所示圖6.6三組心電信號的HR像圖示由上至下分別為來自100號、101號、102號數據的對應圖像（七）時域分析HRV勺時域分析指標NNVGRSDNNR

31、MSSDSDSDNN5解口pNN5（#算結果如表6.1所示。表6.1三組數據的時域指標值數據組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%（八）頻域分析根據快速傅里葉變換（FFD得到HRV®號的頻譜信號和功率譜信號（1）頻譜圖。三組心電數據對應的HRVJ號的頻譜圖如圖6.7所示。石.4.21.S.61110;0;魯登g2Q3S口口4g«Fra-qu&ncy

32、/HzHRvrJjfi/TH15a3Pm冒<DD.aS0.1D15O2J503035D.4口用丐D.SPnQi.JnciyJil-fcLJILl_2D25303S4口45Frequency/HzHRVMilftHI圖6.7消除基線后三組心電信號的頻譜圖圖示由上至下分別為來自100號、101號、102號數據的對應圖像(2)功率譜圖HRVlJjF川I制20BPLZ至£IO15ZQZ53DFrequerKy/Hzhrv.Uj-4aoM30Fraquc-nciy/Hz圖示由上至下分別為來自100號、101號、102號數據的對應圖像(3)根據頻譜計算TP、HRLF、LF/HF,結果如下表

33、6.1所示。表6.1三組數據的時域指標值數據組TP(msms)HF(msms)LF(msms)LF/HF100.dat1900014005710.4150101.dat2060087111001.2063102.dat28300640037000.5812七、結論本文通過使用MATLA歆件，實現了心電信號的讀取、采樣及頻譜分析，并且分離出心率變異性HRVJ號，通過時域分析和頻域分析的方法進行分析，對比資料發現：第101號數據組應該來自正常人，第100號數據組應該來自具有輕微患心血管疾病的患者，第102號數據組則最可能來自患有較為嚴重的心血管疾病的患者八、討論心率變異信號中蘊含著有關心血管調節的

34、重要信息，對HRV!行分析可以間接地定量評價心肌交感、迷走神經的緊張性和均衡性，而且還能分析自主神經系統的活動情況。同時心率變異性分析對多種惡性心律失常的預后判斷和藥物治療效果分析有指導作用。所以，對HRV勺研究能夠極大的促進人類對于心血管疾病的了解，從而在預防、治療心血管疾病等領域取得成果。本文從獲取的心電信號中，提取出有效的心電信號，進而提取出我們所研究的HRV信號，然后對其進行相關分析。在獲取心電信號的過程中，直接使用了來自美國PhysioNet的心電信號資源，在提取出有效心電信號的過程中，首先進行的是對原始心電信號，即三組每組三個相關文件進行相關處理，得出心電信號；其次是觀察得到的初步

35、心電信號進行時序圖的觀察，再進行去工頻干擾和去基線漂移；再其次，對前面得到的信號確定波峰位置，并且剔除異搏，確定出R-R間期；最后，得到我們所需要的HRV信號。對于HRV®號的分析主要分為時域和頻域兩大部分，時域分析中，本文主要對HRV言號的五個時域指標NNVGRSDNNRMSSDSDSDpNN501行了計算；對于頻域部分，對HRV言號的頻域四大指標TP、HF、LF、LF/HF進行了計算，然后作出了HRVB號的頻譜圖和功率譜圖。在研究中使用多種方法進行對信號的處理，經過對比各種方法得到的數據對比，本文選取了效果相對而言比較好，算法相對而言比較簡單的一些方法，比如心電數據的去工頻干擾中

36、使用的是60Hz陷波器濾波法，在心電信號的去基線漂移中使用的是中值濾波法，之后為了保證能對HRV®號進行正確的頻域分析，對HRVJ號進行了線性內插。以上種種方法都是為了使數據具有更高的準確性，基本確保研究沒有在研究對象的處理上的重大失誤。在前人的基礎之上，由于水平、能力有限，本文對心電信號的HRVa行了粗淺的研究。在心電數據的處理上和認識上都還有很多的東西都需要學習，希望以后能對此做進一步的研究。參考文獻1鐘運健.心率變異性(HRV心運動性疲勞診斷中應用的實驗研究：碩士論文.江西南昌：江西師范大學體育學院，2004年.2董紅生.心電波形檢測與心率異變性分析研究：博士論文.蘭州理工大學

37、，2012年3吳學勤.動態心電圖技術與應用.合肥：中國科學技術出版社,1998,22-25.4胡大一，郭成軍，李瑞杰.心率變異性一測量標準，生理釋義與臨床應用一續一：心率變異性的生理研究.中國醫療器械信息，1998,4(1):29-32.5曹細武，鄧親愷.心電圖各波的頻率分析.中國醫學物理學雜志,2001,18(1):46-48.6萬相奎.心電信號分析與虛擬式心電自動分析儀的開發：重慶大學博士學位論文.重慶：重慶大學機械電子工程,2005,16-17.7LevkovC,MichovG,IvanovRetal.Subtractionof50Hzinterferencefromtheelectro

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

39、isprogrammreadsECGdatawhicharesavedinformat212.%(e.g.,100.datfromMIT-BIH-DB,cu01.datfromCU-DB,.)%Thedataaredisplayedinafiguretogetherwiththeannotations.%TheannotationsaresavedinthevectorANNOT,thecorresponding%times(inseconds)aresavedinthevectorATRTIME.%Theannotationsaresavedasnumbers,themeaningofthe

40、numberscan%befoundinthecodetable"ecgcodes.h".%ANNOTonlycontainsthemostimportantinformation,whichisdisplayed%withtheprogramrdann()inthe3rdrow.%The4thto6throwarenotsavedinANNOT.%createdonFeb.27,2003by%RobertTratnig(VorarlbergUniversityo

41、fAppliedSciences)%(email:rtratniggmx.at),%algorithmisbasedonaprogramwrittenby%KlausRheinberger(UniversityofInnsbruck)%(email:klaus.rheinbergeruibk.ac.at)%clc;clearall;%SPECIFYDATAPATH='C:Usersasus你心率變異性(HRV信號的提取及時頻域分析ecg'%path,wheredataaresavedHEADERFILE='100.hea'%header-fileintextfo

42、rmatATRFILE='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=

43、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);%

44、bitresolutionzerovalue(k)=A(4);%integervalueofECGzeropointfirstvalue(k)=A(5);%firstintegervalueofsignal(totestforerrors)end;fclose(fid1);clearA;%LOADBINARYDATAifdformat-=212,212,error('thisscriptdoesnotapplybinaryformatsdifferentto212.');end;signald=fullfile(PATH,DATAFILE);%datainformat212fi

45、d2=fopen(signald,'r');A=fread(fid2,3,SAMPLES2READuint8')'%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)=bi

46、tshift(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)；M(:,2)=(M(:,2)-zerovalue)；M=M'M(1)二口；sM=size(M)

47、;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$>LOADINGDATAFINISHEDn');%LOADATTRIBUTESDATAatrd=full

48、file(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,

49、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=ANN

50、OT;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);%DISPLAYDATA修改,'LineWidth',2figure(1);clf,boxon,hol

51、donplot(TIME,M(:,1),'r');DATA1=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,

52、'n$>DISPLAYINGDATAFINISHEDn');%fprintf(1,'n$>ALLFINISHEDn');%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=

53、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(,去除工頻干擾后頻

54、譜圖');%-多項式擬合%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時序圖'

55、);%-確定波峰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,:),'