心率变异性HRV信号提取及时频域分析(包含程序)

./课程设计报告题目：心率变异性〔HRV信号的提取及时频域分析专业：生物医学工程班级：XXXXXXX学号：XXXXXXX姓名：XXXXXXX指导教师：XXXXXXXXXXXXX大学XXXXX学院2016年9月29日开题背景〔一HRV简介传统的医学观点认为,正常的心率为规则的窦性节律；后来发现在健康状态下,许多生理系统中存在自然的变异性,人的心率正常情况下也是呈不规则性变化的,而心率变异就是指窦性心率的这种波动变化的程度。心率变异性〔HeartRateVariability,HRV是指逐次心搏间期之间的微小变异特性。在生理条件下,HRV的产生主要是由于心脏窦房结自律活动通过交感和迷走神经,神经中枢,压力反射和呼吸活动等因素的调节作用,使得心脏每搏间期一般存在几十毫秒的差异。<二>HRV的研究现状心率变异性<HRV>是近年来比较受关注的无创性心电监测指标之一,对HRV的生理和病理意义进行了广泛和深入的研究,其结果表明心率变异信号中蕴含着有关心血管调节的重要信息,对HRV进行分析可以间接地定量评价心肌交感、迷走神经的紧性和均衡性,而且还能分析自主神经系统的活动情况。心率变异性还可以作为一个独立的心源性猝死危险性的预测指标。同时心率变异性分析对多种恶性心律失常的预后判断和药物治疗效果分析有指导作用。所以,对HRV的研究能够极大的促进人类对于心血管疾病的了解,从而在预防、治疗心血管疾病等领域取得成果。〔三HRV的研究方法随着对HRV研究的不断深入,其蕴含的生理病理信息将进一步被揭示,使得HRV有更多的应用空间和应用价值。目前,心率变异性分析方法主要有时域分析法、频域分析法、时频分析法以及非线性分析法[1]。〔四HRV的临床应用〔1心脏性猝死〔SCD预测：由于HRV是反映自主神经力的最敏感的指标,因此HRV降低是预测心脏性猝死最有价值的独立指标。〔2急性心肌梗塞后患者危险性评估：HRV的降低是预测急性心肌梗塞后患者发生心脏性猝死和恶性心律失常危险的重要独立指标。一般建议在梗塞后一周开始进行HRV的检测。HRV在梗塞后立即降低,并在几周开始恢复〔2周后逐渐回升,大约6-12个月恢复正常。因此,多次测定HRV可能比单次测定价值更大。梗塞后HRV恢复的快慢对以后死亡的危险性也有预测价值。〔3对糖尿病患者自主神经系统损伤的评估：糖尿病患者不论病情轻重,均存在不同程度的自主神经功能紊乱。HRV是判断糖尿病患者是否伴有自主神经系统损害最准确,最敏感的指标。〔4心力衰竭〔CHF患者危险性评估。〔5心率变异性生物反馈疗法：对于不孕人群受孕几率提高、怀孕人群孕期焦虑症改善、产后人群产后抑郁症情况缓解,起到很好的作用。〔6其它临床应用围：心绞痛、高血压、心肌病、非缺心脏病所致的慢性严重二尖瓣返流、二尖瓣脱垂、心律失常、血管迷走性晕厥等心血管疾病。课题目的基本掌握心电信号〔ECG的测量、数据采集的方法。学会使用MATLAB对ECG信号进行相关处理分析。主要包括从ECG信号中提取出所需的HRV信号,并分别对其进行时域、频域、功率谱上的分析。掌握HRV信号的时频域参数的意义,以及对其进行分析的基本方法。课题研究的主要容〔一从网上下载正常人的心电信号以及各种病人的心电信号〔ECG数据。〔二首先HRV信号的提取,主要包括去除干扰、准确确定R波波峰位置、剔除异搏、确定R-R间期、线性插,并且绘出HRV信号曲线。〔三对HRV信号的时域分析,对HRV信号的频谱图和功率谱图分析。原理和方法ECG信号的采集本文主要使用100.hea、100dat、100.atr,101.hea、101.dat、101.atr,102.hea、102.dat、102.atr这三组数据来对HRV进行研究。实验数据来源于PhysioNet。PhysioNet是一个基于Web的复杂生理和生物医学信号的研究资源,其网址为http://.physionet.org。PhysioNet由PhysioNet,PhysioBank和PhysioToolkit三个相互关联的部分组成。数据库中数据来源于正常人、各种病人〔如心脏猝死、心力衰竭、心律失常、癫痫、睡眠呼吸暂停综合症等及运动、休息等不同状态下的数据,样本选取围广泛,其部分数据都进行详细的注释,并将数据被划分为3类,即Classl：专家已经作出了标注；Class2：原始数据；Class3：处于研究进展之中。因此,PhysioBank数据库中的数据足已满足生物医学各领域研究者的需要。PhysioBank数据库中的每一条数据记录包括至少三类文件,头文件〔.hea、数据文件〔.dat和注释文件〔.atr,.al,.aiM等。头文件是描述数据属性的文本文件,其容包括记录名、信号数目、贮存格式、信号数量和类型、采样频率、数字化特征、记录的持续时间和起始时间等信息。一般可由PhysioToolkit软件库的WFDB库函数的getinfo、putinfo函数读和写的字符。数据文件是定义了相应存贮格式的数字化采样点的二进制存储文件。数据存贮格式在头文件中说明,一个数据组有相同的数据存贮格式,常用的是8位和16位格式。注释文件是记录了对信号特征的注释信息[2]。ECG信号的特征典型心电信号波形心脏搏动及其节律性是人体生命和生理状态的重要标志之一。心电生理学的研究表明,心电信号来源于心肌细胞的生物电变化。心肌细胞的电激动称为除极,心肌细胞恢复为静息状态称为复极,心电信号的产生与心肌细胞的除极和复极过程密切相关。心脏电激动起源于窦房结,沿特化的心脏传导系统下传,其传播方向、途径、次序及时间存在一定的规律。若心脏不能及时发出电激动,则心脏陷于停博。人体体液中充满电解质,具有导电性能,心脏电激动过程产生的有序生物电变化通过体液传至身体表面使身体各部位出现有规律而各向异性的电位变化,通过测量电极采集体表特定点电位变化,并放大、显示及记录,即为体表心电信号,也即是通常的ECG信号[3]。一个心动周期正常心电信号波形如图4.1所示。它是由特征波及其特征间期组成,每个心动周期包含一个P波,一个QRS波群和一个T波,有时还会出现一个小的U波。特征波及特征间期的含义如下：图4.1典型的心电信号波形P波：由左右心房的除极过程引起,其波形小而圆钝,时宽为0.08s-0.11s,波幅不超过0.25mV。QRS波：反映左右心室除极产生的电位变化,在P波之后出现,为心电信号中最高大和最快速的波形。典型的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波开始到QRS波开始之间的时程,反映激动由窦房结产生经由结间束、房室交界和左右束支抵达心室,并引起心室兴奋所需要的时间,又称为房室传导时间。正常为0.12s-0.2s。当发生房室传导阻滞时,PR间期增长。如当PR间期>0.21s。则为Ⅰ度房室传导阻滞。QT间期：从QRS波开始到T波终点的时程,反映心室除极和复极时间的总和。许多因素可影响QT间期,如心肌缺血、低血钾、低血钙等可使QT间期延长,QT间期延长使心室肌复极不均一,易诱发折返激动,导致严重室性心率失常。QT间期随受心率变化的影响,心率越慢,QT间期越长；心率越快,QT间期越短。通常用QTc间期修正心率对QT间期影响,正常QTc间期小于0.43~0.44s。ST段：指从QRS波群终止点到T波起点之间的波形线段,反映心室部分己完全进入去极化状态,正常时与基线平齐。PP间期：相邻P波之间的间距称为PP间期,反映心房率。正常情况下,PP间期与RR间期一致。在Ⅱ度或Ⅱ度以上房室传导阻滞和某些心率失常,两者可不一致。RR间期：相邻QRS波群之间的间距称为RR间期,反映心室率。正常情况下,RR间期与PP间期一致[2]。在心电信号的测试中,对电极的放置部位和导联的连接方式临床有明确的规定。目前,国际公认的是标准12导联,包括心电标准导联〔I、II、III、加压单极肢体导联〔aVR、aVL、aVF及胸导联〔Vl～V6,共有12个导联,具体可参考文献[4]。<三>典型心电信号的能量〔频谱图分布典型的心电信号的整个心动周期的频谱估计图如图4.2[2]所示,可以明显看出心电信号各波的能量主要集中在低频区域,且随着频率的增高,相应的能量逐渐降低。心电信号的整体频谱围在0.05Hz～100Hz,但能量主要集中在0.5～45Hz,能量的最高点在8～15Hz附近；QRS波群的频谱带宽为3～40Hz,积聚了将近99%的能量,波峰能量集中在6～18Hz附近,P波的频谱带宽为0～18Hz,波峰能量集中在5～12Hz；T波的频谱带宽为0～8Hz,波峰能量集中在0～8Hz区间[5]。图4.2典型的心电信号频谱能量分布<四>ECG信号的噪声分析在采集、放大及传输心电信号的过程中,由于受人体、采集仪器、电磁环境、操作水平等的影响,不可避免会有许多干扰耦合到心电信号,主要干扰表现形式如下：〔1电源工频干扰产生的原因主要由于电源磁场作用于心电图仪的导联与人体之间的环形电路所致,表现为心电信号上有明显的正弦波或正弦波的叠加信号,其频率为60Hz工频及其谐波构成,幅度较低。〔2基线漂移产生的原因主要由于人体呼吸运动、电极接触不良等因素所导致。表现为心电信号上叠加缓慢变化的信号,其频率一般小于1Hz,幅度为ECG峰-峰的15%。〔3肌电干扰产生的原因主要由于人体活动,肌肉紧所引起的干扰。表现为不规则的快速变化波形,其频率围较宽,一般在5～2kHz之间,幅度为毫伏级。〔4运动伪迹产生的原因主要由于电极与人体间轻微移动或抖动而引入的干扰,表现为信号基线的短暂变化,但不是基线的跃变,其持续时间为100～500ms,频率一般在7Hz以下,幅度较大。〔5其他随机噪声心电信号还受到其他的随机噪声和环境干扰的影响,如加性白噪声、极化噪声、仪器部噪声等。由此可见,噪声信号基本覆盖了有用的心电信号的全频率围。而其中60Hz及其倍频附近的工频干扰、1Hz以下的基线漂移以及肌电干扰噪声是最主要的干扰源,在心电信号预处理中必须消除或抑制,以提高心电信号的信噪比[2]。〔五ECG信号去噪声方法对于心电信号预处理一般要从硬件电路优化设计和软件数字滤波器的设计两个方面考虑。根据心电信号的频谱分布特点,在硬件方面来看,消除基线漂移的干扰应考虑分别设计下限频率为0.5Hz的高通滤波器；消除肌电等高频干扰应考虑设计上限频率为100Hz的低通滤波器,同时还应考虑设计60Hz的陷波器来滤除工频干扰。由于硬件滤波器的元器件精度、稳定性要求以及难以实现严格线性相位等问题,使得仅采用硬件滤波不能满足滤波性能的要求。随着数字信号处理技术的发展,设计高精度、高可靠性、简洁灵活的数字滤波器成为可能。目前,滤除心电信号干扰更多采用的是具有线性相位的数字滤波方法,并已逐渐显示出取代硬件滤波器的趋势。〔1工频干扰对于工频干扰,主要的消除方法有平均滤波器、梳状滤波器、Levkov滤波器、自适应工频陷波器等方法。平均滤波器法是较早被应用的数字滤波方法,其特点是算法简单,处理速度快,滤波效果较好。由于平均滤波算法的实时性好,被广泛应用于实时心电监护设备中或基于单片机心电数据采集系统中。从实际滤波效果看,滤波器能较好地滤除了60Hz工频,但对ECG中的QRS波有较大削峰,信号衰减较大。因此,在实际应用中受到一定的限制。梳状滤波器的名称源于该滤波器的传递函数的幅频特性形如梳状。其特点是运算量小、形式灵活简单,是一种快速的数字滤波器。梳状滤波器被广泛应用去除工频及其高次谐波的干扰。选择合适的梳状滤波参数,可保证滤波器对有用信号基本无衰减和滤波后基本不产生相移[6]。改进的梳状滤波器可用于消除ECG的工频干扰且能有效地防止信号的失真。Levkov滤波器是Levkov在1984年提出一种对ECG信号的线性段和非线性段采用不同处理方式的数字滤波方法。其基本原理是从原始ECG信号中直接减去在该线性段中确定的干扰信号的幅值[7]。改进的Levkov滤波方法中引入了判定线性段和非线性段的指标M。算法的滤波效果与选择合适的M值有关。改进Levkov滤波法算法简单,参数可调,运算量小,可实现实时心电信号滤波处理,但在噪声较强时,将会造成QRS波群的削峰现象,滤波效果不能令人满意[8]。自适应工频陷波技术特点是只需要很少或根本不需要任何有关信号和噪声的先验知识的情况下,直接利用观测数据不断递归更新自身的参数,以逐步逼近某一最优值,能自动跟踪信号频率漂移,具有非常好的适应性[9]。自适应的数字陷波器的权值调节主要有最小均方误差〔LMS、最小递推二乘法<RLS>算法[10]。一般来说,RLS自适应算法与LMS自适应算法相比,权值有更快的收敛速度。缺点是每次迭代计算量较大,由于在生物医学应用中数据的计算量不是太大,RLS自适应算法可获得更好性能。〔2基线漂移的消除 基线漂移可严重影响心电信号的分析与处理。滤波法和拟合基漂法是去除心电基线漂移的两类主要算法。具体而言有形态滤波、中值滤波、线性相位FIR滤波、小波变换、形态和小波相结合、中值滤波和小波相结合等去除基漂的滤波方法,以及导数法与坐标法这两类基漂拟合方法。研究表明滤波法实现简单,但精度较低,分段拟合法去除基漂的效果更好,有自适应的特性而且失真较小。无论是实时性还是准确性,基漂拟合法都比滤波法更有优势但某些场合下,基漂拟合点的提取十分困难,此时只能使用滤波法[11]。对于基线拟合法,通常多采用简洁、快速的分段拟合法。其基本原理是通过对采样的心电数据分段进行三次多项式或抛物线拟合,获取基线漂移曲线,然后用原始心电曲线减去拟合曲线来消除基线漂移。算法的关键是插值点选择,插值点应该是心电信号的零电位点,通常采用TP段。对于FIR滤波器只有零点,不易获得较好的通带与阻带衰减特性。要取得好的衰减特性,一般要求H<z>的阶次要高。对于小波变换,由于其提供了一个在时频平面上可调的分析窗口。信号的时、频分辨率可以随分析任务的需要作出相应调整。在低频时小波变换的时间分辨率较差,而频率分辨率较高；在高频时小波变换的时间分辨率较高,而频率分辨率较低,自动满足了信号分析。信号x<t>在不同尺度分解的逼近信号和细节信号具有不同频带,而在不同的尺度上一定包含有x<t>的不同频率信息。所以若对包含有不同频率信息的细节信号进行阈值的取舍,丢弃其中的干扰成分后,再对信号进行重构,则重构信号的干扰得到了非常好的抑制,即为小波变换消噪的基本原理[2]。通过对集中方法的实验效果来选择合适的消除基线漂移的方法。〔六R波峰位置、剔除异搏、确定RR间期〔1波峰位置找出波峰位置的基本想法是寻找一个采样点,其电压值既大于前一个采样点的值也大于后一个采样点的值。之后再将P波与T波中的伪峰值去除,剩下的就是QRS波的峰值,即R波峰。可以通过先屏蔽掉小于某一电压的所有采样点,这样就将P波和T波中的伪峰值提前去除,之后再寻找峰值。以此保证得到的峰值处于QRS波中,该峰值即为R波峰值。剔除异搏对得到R波的位置,检查其中的异常点,去掉异搏。确定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：全程相邻NN间期长度之差的标准差,单位为ms。即其中,。pNN50：在全部NN间期的记录中,相邻的NN间期之差大于50ms的个数与总的NN间期的个数的比,以百分比表示。其中NNVGR用于评估心率总体变化水平,；SDNN用于评估心率总体变化的大小,即交感及迷走神经力大小,；SDANN用于评估心率变化中的长期慢变化成分,即交感神经力大小。RMSSD及pNN50反映心率快变化成分的大小,即副交感神经力的敏感指标。<八>HRV信号的频域分析：HRV时域分析的指标大多用于描述HRV整体的大小,不能仔细地分析交感神经和迷走神经各自的活动的情况。而HRV频谱分析则可以弥补这个缺点。另一角度,即频谱分析的角度来分析心率变化的规律。它与时域分析既有相关性,又能揭示出心率的更复杂的变化规律。FFT是经典谱估计方法,算法简单。输入和输出信号能量有线性关系,但对信号要作周期延拓假定,短数据谱分辨率较低,并有能量泄露现象。最近更仔细的研究发现,正常人基础状态下心率谱曲线在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反映HRV大小,LF/HF比值反映自主神经系统的平衡状态,基本上代表经交感神力的高低。步骤〔包括计算程序流程框图总程序流程图。详细程序间附件1。获取心电信号获取心电信号去除工频干扰和基线漂移去除工频干扰和基线漂移确定波峰位置确定波峰位置剔除异搏剔除异搏确定R-R间期确定R-R间期求出HRV信号求出HRV信号进行时域计算分析进行线性内插进行时域计算分析进行线性内插进行频域计算分析进行功率谱分析进行频域计算分析进行功率谱分析HRV信号NNVGRSDNNHRV信号NNVGRSDNNRMSSDSDSDpNN50频谱图TP频谱图TPHFLFLF/HF结论结论获取心电信号将下载得到三组心电信号：100、101、102,通过KlausRheinberger编写的读取心电的MATLAB程序ecg.load.m,加载原始心电信号,从而准确提取相应的心电信号。去除工频噪声干扰对信号进行频域分析,发现信号主要集中在低频端,并且在60Hz处有微弱信号,通过判断可得出该频率信号为工频噪声。本文使用FIR滤波器将噪声滤去。去基线漂移对比两大类型发现：滤波法实现简单,但精度较低；分段拟合法去除基漂的效果更好,有自适应的特性而且失真较小。无论是实时性还是准确性,基漂拟合法都比滤波法更有优势；但某些场合下,基漂拟合点的提取十分困难,此时只能使用滤波法。通过进行初步实验中比较各种方法的实际效果,发现使用中值滤波法的效果最好,故本文主要采用其进行去基线漂移处理。找出波峰、剔除异搏、确定RR间期对HRV信号进行时域分析对前面得到的HRV信号进行时域分析,其中的技术指标有NNVGR、SDNN、RNSSD、SDSD、pNN50等。<七>线性插对提取出来的HRV信号经线性插获得等间隔〔R-R间期均值R-R间期时间序列,然后再进行频谱分析。线性插的公式为。式中X<n>为插值,和分别为插值前、后的R-R间期序列值,和分别为X〔n所在位置与和所在位置的时间间距。〔八根据快速傅里叶变换〔FFT得到HRV信号的频谱信号和功率谱信号。主要有TP、HF、LF、LF/HF、频谱图、功率谱图等。结果通过MATLAB程序分别得到100、101和102三组ECG信号数据,然后对其进行相关分析,具体结果以及图像如下所述。〔一ECG信号的提取通过MATLAB程序分别得到100、101和102三组ECG信号数据。其时序图像具体如图6.1所示。图6.1三组原始心电信号图示由上至下分别为来自100号、101号、102号数据的对应图像〔二原始心电频谱图对三组信号通过皮谱图进行频域分析。从图6.2中可以看出,信号主要集中在20Hz以下的低频端,并且在60Hz附近有一个明显的幅度值上升,可以判断其为工频干扰。图6.2三组心电信号的频谱图图示由上至下分别为来自100号、101号、102号数据的对应图像去除工频干扰使用FIR滤波器法将工频噪声去除。去除工频干扰之后频谱图如图6.3所示图6.3三组心电信号去除工频干扰后的频谱图图示由上至下分别为来自100号、101号、102号数据的对应图像去除基线漂移从图6.1中可以看出这些心电信号均有明显的基线漂移。使用中值滤波法去除基线漂移。消除基线漂移后的ECG信号时序图如图6.4所示。图6.4消除基线后三组心电信号的频谱图图示由上至下分别为来自100号、101号、102号数据的对应图像确定波峰确定波峰结果如图6.5所示。其中红色点表示波峰。图6.5确定三组心电信号的R波峰图示由上至下分别为来自100号、101号、102号数据的对应图像确定R-R间期R-R间期与R-Ri的关系,亦即HRV的图像如图6.6所示。图6.6三组心电信号的HRV图像图示由上至下分别为来自100号、101号、102号数据的对应图像〔七时域分析HRV的时域分析指标NNVGR、SDNN、RMSSD、SDSD、NN50和pNN50计算结果如表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%〔八频域分析根据快速傅里叶变换〔FFT得到HRV信号的频谱信号和功率谱信号。频谱图。三组心电数据对应的HRV信号的频谱图如图6.7所示。图6.7消除基线后三组心电信号的频谱图图示由上至下分别为来自100号、101号、102号数据的对应图像功率谱图图6.8HRV信号的功率谱图图示由上至下分别为来自100号、101号、102号数据的对应图像根据频谱计算TP、HF、LF、LF/HF,结果如下表6.1所示。表6.1三组数据的时域指标值数据组TP<ms×ms>HF<ms×ms>LF<ms×ms>LF/HF100.dat1900014005710.4150101.dat2060087111001.2063102.dat28300640037000.5812结论本文通过使用MATLAB软件,实现了心电信号的读取、采样及频谱分析,并且分离出心率变异性HRV信号,通过时域分析和频域分析的方法进行分析,对比资料发现：第101号数据组应该来自正常人,第100号数据组应该来自具有轻微患心血管疾病的患者,第102号数据组则最可能来自患有较为严重的心血管疾病的患者讨论心率变异信号中蕴含着有关心血管调节的重要信息,对HRV进行分析可以间接地定量评价心肌交感、迷走神经的紧性和均衡性,而且还能分析自主神经系统的活动情况。同时心率变异性分析对多种恶性心律失常的预后判断和药物治疗效果分析有指导作用。所以,对HRV的研究能够极大的促进人类对于心血管疾病的了解,从而在预防、治疗心血管疾病等领域取得成果。本文从获取的心电信号中,提取出有效的心电信号,进而提取出我们所研究的HRV信号,然后对其进行相关分析。在获取心电信号的过程中,直接使用了来自美国PhysioNet的心电信号资源,在提取出有效心电信号的过程中,首先进行的是对原始心电信号,即三组每组三个相关文件进行相关处理,得出心电信号；其次是观察得到的初步心电信号进行时序图的观察,再进行去工频干扰和去基线漂移；再其次,对前面得到的信号确定波峰位置,并且剔除异搏,确定出R-R间期；最后,得到我们所需要的HRV信号。对于HRV信号的分析主要分为时域和频域两大部分,时域分析中,本文主要对HRV信号的五个时域指标NNVGR、SDNN、RMSSD、SDSD、pNN50进行了计算；对于频域部分,对HRV信号的频域四大指标TP、HF、LF、LF/HF进行了计算,然后作出了HRV信号的频谱图和功率谱图。在研究中使用多种方法进行对信号的处理,经过对比各种方法得到的数据对比,本文选取了效果相对而言比较好,算法相对而言比较简单的一些方法,比如心电数据的去工频干扰中使用的是60Hz陷波器滤波法,在心电信号的去基线漂移中使用的是中值滤波法,之后为了保证能对HRV信号进行正确的频域分析,对HRV信号进行了线性插。以上种种方法都是为了使数据具有更高的准确性,基本确保研究没有在研究对象的处理上的重大失误。在前人的基础之上,由于水平、能力有限,本文对心电信号的HRV进行了粗浅的研究。在心电数据的处理上和认识上都还有很多的东西都需要学习,希望以后能对此做进一步的研究。

参考文献[1]钟运健.心率变异性<HRV>在运动性疲劳诊断中应用的实验研究:[硕士论文].：师大学体育学院,20XX.[2]董红生.心电波形检测与心率异变性分析研究:[博士论文].理工大学,20XX[3]吴学勤.动态心电图技术与应用.:中国科学技术,1998,22-25.[4]胡大一,郭成军,瑞杰.心率变异性一测量标准,生理释义与临床应用—续一:心率变异性的生理研究.中国医疗器械信息,1998,4<1>:29-32.[5]细武,邓亲恺.心电图各波的频率分析.中国医学物理学杂志,2001,18<1>:46-48.[6]万相奎.心电信号分析与虚拟式心电自动分析仪的开发:[大学博士学位论文].:大学机械电子工程,2005,16-17.[7]LevkovC,MichovG,IvanovRetal.Subtractionof50Hzinterferencefromtheelectrocardiogram.Med.&Biol.,Eng.&Comp.,1984,22<4>:371-373.[8]京霞,白延强,玉星.一种抑制心电信号50Hz工频干扰的改进Levkov方法.航天医学与医学工程,2000,13<3>:196-199.[9]福生,吕扬生.生物医学信号的处理和识别.:科技翻译出版司,1997,378-379.[10]高鹰,胜利.一种变步长LMS自适应滤波算法及分析.电子学报.2001,29<8>:1094-1097.[11]延军,严洪,王增丽.心电基线漂移去除方法的比较研究.航天医学与医学工程.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信号的提取及时频域分析\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<'去除工频干扰后频谱图'>;%多项式拟合%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时序图'>;%确定波峰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-