课程设计报告
题 目： 心率变异性（HRV）信号的提取及时频域分析 专 业： 生物医学工程 班 级： XXXXXXX 学 号： XXXXXXX 姓 名： XXXXXXX 指导教师： XXXXXXX
XXXXXX大学 XXXXX学院
2016年 9月 29日
一、 开题背景
（一）HRV简介
传统的医学观点认为，正常的心率为规则的窦性节律；后来发现在健康状态下，许多生理系统中存在自然的变异性，人的心率正常情况下也是呈不规则性变化的，而心率变异就是指窦性心率的这种波动变化的程度。心率变异性（Heart Rate Variability，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: //www.physionet. org。 PhysioNet由PhysioNet, PhysioBank和PhysioToolkit三个相互关联的部分组成。数据库中数据来源于正常人、各种病人（如心脏猝死、心力衰竭、心律失常、癫痫、睡眠呼吸暂停综合症等）及运动、休息等不同状态下的数据，样本选取范围广泛，其中大部分数据都进行详细的注释，并将数据被划分为3类，即Class l：专家已经作出了标注；Class 2：原始数据；Class3：处于研究进展之中。因此，PhysioBank数据库中的数据足已满足生物医学各领域研究者的需要。PhysioBank数据库中的每一条数据记录包括至少三类文件，头文件（.hea）、数据文件（.dat）和注释文件（.atr,.al,.aiM等）。头文件是描述数据属性的文本文件，其内容包括记录名、信号数目、贮存格式、信号数量和类型、采样频率、数字化特征、记录的持续时间和起始时间等信息。一般可由PhysioToolkit软件库的WFDB库函数的getinfo、putinfo函数读和写的字符。数据文件是定义了相应存贮格式的数字化采样点的二进制存储文件。数据存贮格式在头文件中说明，一个数据组有相同的数据存贮格式，常用的是8位和16位格式。注释文件是记录了对信号特征的注释信息[2]。
（二） ECG信号的特征
（1） 典型心电信号波形
心脏搏动及其节律性是人体生命和生理状态的重要标志之一。心电生理学的研究表明，心电信号来源于心肌细胞的生物电变化。心肌细胞的电激动称为除极，心肌细胞恢复为静息状态称为复极，心电信号的产生与心肌细胞的除极和复极过程密切相关。心脏电激动起源于窦房结，沿特化的心脏传导系统下传，其传播方向、途径、次序及时间存在一定的规律。若心脏不能及时发出电激动，则心脏陷于停博。人体体液中充满电解质，具有导电性能，心脏电激动过程产生的有序生物电变化通过体液传至身体表面使身体各部位出现有规律而各向异性的电位变化，通过测量电极采集体表特定点电位变化，并放大、显示及记录，即为体表心电信号，也即是通常的 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～40 Hz，积聚了将近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波峰值。 （2） 剔除异搏
对得到R波的位置，检查其中的异常点，去掉异搏。 （3） 确定RR间期
用后一个峰值的采样时间减去前一个峰值的采样时间就可以得到RR间期的值。将信号以RR值为纵坐标，R-Ri值为横坐标，得到HRV的图像。
（七）HRV信号的时域分析
通常使用的HRV时域检测指标有5项：NNVGR、SDNN、RMSSD、SDSD和pNN50。 上述5个指标的定义分别为：
NNVGR：全部正常窦性心博间期（NN）的平均值，单位为ms。即NNVGR?1N?RR。
i?1NSDNN：标准差，即全部NN间期的标准差，单位为ms。即SDNN?1N?(RRi?1Ni?RR)2。
RNSSD：全程相邻NN间期之差的均方根值，单位为ms。即RNSSD?1N?1(RRi?1?RRi) 。 ?N?1i?1SDSD：全程相邻NN间期长度之差的标准差，单位为ms。即SDSD?其中RRi＇?RRi?RRi?1，RR＇?RR?RRi?1。
1N?1(RRi＇?RR＇)2 ?N?1i?1??2pNN50：在全部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。
HRV信号NNVGR SDNN RMSSD SDSD pNN50 进行时域计算分析 获取心电信号 去除工频干扰和基线漂移 确定波峰位置 剔除异搏 确定R-R间期 求出HRV信号 进行线性内插 进行功率谱分析 进行频域计算分析 结论 频谱图 TP HF LF LF/HF （二） 获取心电信号 将下载得到三组心电信号：100、101、102，通过Klaus Rheinberger编写的读取心电的MATLAB程序ecg.load.m，加载原始心电信号，从而准确提取相应的心电信号。 （三） 去除工频噪声干扰
对信号进行频域分析，发现信号主要集中在低频端，并且在60Hz处有微弱信号，通过判断可得出该频率信号为工频噪声。本文使用FIR滤波器将噪声滤去。 （四） 去基线漂移
对比两大类型发现：滤波法实现简单，但精度较低；分段拟合法去除基漂的效果更好,有自适应的特性而且失真较小。无论是实时性还是准确性,基漂拟合法都比滤波法更有优势；但某些场合下,基漂拟合点的提取十分困难,此时只能使用滤波法。通过进行初步实验中比较各种方法的实际效果，发现使用中值滤波法的效果最好，故本文主要采用其进行去基线漂移处理。
（五） 找出波峰、剔除异搏、确定RR间期 （六） 对HRV信号进行时域分析
对前面得到的HRV信号进行时域分析，其中的技术指标有NNVGR、SDNN、RNSSD、SDSD、pNN50等。 (七) 线性内插
对提取出来的HRV信号经线性内插获得等间隔（R-R间期均值）R-R间期时间序列，然后再进行频谱分析。线性内插的公式为X(n)?d2RRidRRi?1。式中X(n)为?1d1?d2d1?d2插值，RRi和RRi?1分别为插值前、后的R-R间期序列值，d1和d2分别为X（n）所在位置与RRi和RRi?1所在位置的时间间距。
（八）根据快速傅里叶变换（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) NN50 pNN50 100. 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%
（八）频域分析
根据快速傅里叶变换（FFT）得到HRV信号的频谱信号和功率谱信号。 （1） 频谱图。三组心电数据对应的HRV信号的频谱图如图6.7所示。
图 6.7 消除基线后三组心电信号的频谱图 图示 由上至下分别为来自100号、101号、102号数据的对应图像 （2） 功率谱图
图 6.8 HRV信号的功率谱图 图示 由上至下分别为来自100号、101号、102号数据的对应图像
（3） 根据频谱计算TP、HF、LF、LF/HF，结果如下表6.1所示。
表6.1 三组数据的时域指标值 数据组 TP(ms×ms) HF(ms×ms) LF(ms×ms) LF/HF 100.dat 19000 1400 571 0.4150 101.dat 20600 871 1100 1.2063 102.dat 28300 6400 3700 0.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)在运动性疲劳诊断中应用的实验研究:[硕士论文].江西
南昌：江西师范大学体育学院，2004年.
[2] 董红生.心电波形检测与心率异变性分析研究:[博士论文].兰州理工大学，2012年 [3] 吴学勤.动态心电图技术与应用.合肥:中国科学技术出版社,1998,22-25. [4] 胡大一,郭成军,李瑞杰.心率变异性一测量标准,生理释义与临床应用—续一:心率
变异性的生理研究.中国医疗器械信息,1998, 4(1):29-32.
[5] 曹细武,邓亲恺.心电图各波的频率分析.中国医学物理学杂志,2001,18(1):46-48. [6] 万相奎.心电信号分析与虚拟式心电自动分析仪的开发:[重庆大学博士学位论文].
重庆:重庆大学机械电子工程,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.
[8] 孙京霞, 白延强, 杨玉星. 一种抑制心电信号 50Hz 工频干扰的改进 Levkov 方
法.航天医学与医学工程,2000,13(3):196-199.
[9] 杨福生,吕扬生.生物医学信号的处理和识别.天津:天津科技翻译出版司,1997,
378-379.
[10] 高鹰,谢胜利.一种变步长LMS 自适应滤波算法及分析. 电子学报.2001,29(8):
1094-1097.
[11] 李延军,严 洪,王增丽.心电基线漂移去除方法的比较研究.航天医学与医学工
程.2009年第22卷第五期：381-386.
附件1：
% This programm reads ECG data which are saved in 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 meaning of the numbers can % be found in the codetable \%
% ANNOT only contains the most important information, which is displayed % with the program rdann (available on www.physionet.org) in the 3rd row. % The 4th to 6th row are not saved in ANNOT. % %
% created on Feb. 27, 2003 by
% Robert Tratnig (Vorarlberg University of Applied Sciences) % (email: rtratnig@gmx.at), %
% algorithm is based on a program written by % Klaus Rheinberger (University of Innsbruck) % (email: klaus.rheinberger@uibk.ac.at) %
% ------------------------------------------------------------------------- clc; clear all;
%------ SPECIFY DATA ------------------------------------------------------ PATH= 'C:\\Users\\asus你\\Desktop\\心率变异性（HRV）信号的提取及时频域分析
\\ecg'; % path, where data are saved
HEADERFILE= '100.hea'; % header-file in text format ATRFILE= '100.atr'; % attributes-file in binary format DATAFILE='100.dat'; % data-file
SAMPLES2READ=30000; % number of samples to be read % in case of more than one signal: % 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 signals sfreq=A(2); % sample rate of data clear 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); % first 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 212 fid2=fopen(signald,'r');
A= fread(fid2, [3, SAMPLES2READ], 'uint8')'; % matrix with 3 rows, each 8 bits
long, = 2*12bit fclose(fid2);
M2H= bitshift(A(:,2), -4); M1H= bitand(A(:,2), 15);
PRL=bitshift(bitand(A(:,2),8),9); % sign-bit PRR=bitshift(bitand(A(:,2),128),5); % sign-bit M( : , 1)= bitshift(M1H,8)+ A(:,1)-PRL;
M( : , 2)= bitshift(M2H,8)+ A(:,3)-PRR;
if M(1,:) ~= firstvalue, error('inconsistency in the first bit values'); end; switch nosig case 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( : , 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 PRR PRL;
fprintf(1,'\\\\n$> LOADING DATA FINISHED \\n');
%------ LOAD ATTRIBUTES DATA ---------------------------------------------- atrd= fullfile(PATH, ATRFILE); % attribute file with annotation data fid3=fopen(atrd,'r');
A= fread(fid3, [2, inf], 'uint8')'; fclose(fid3); ATRTIME=[]; ANNOT=[]; sa=size(A); saa=sa(1); i=1;
while i<=saa
annoth=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=bitshift(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 = EOF clear A;
ATRTIME= (cumsum(ATRTIME))/sfreq; ind= find(ATRTIME <= TIME(end)); ATRTIMED= ATRTIME(ind); ANNOT=round(ANNOT); ANNOTD= ANNOT(ind);
%------ DISPLAY DATA ------------------------------------------------- figure(1); clf, box on, hold on
plot(TIME, M(:,1),'r'); %修改,'LineWidth',2
DATA1=M(:,1)'; % if nosig==2
% plot(TIME, M(:,2),'b'); % end;
% for k=1:length(ATRTIMED)
% text(ATRTIMED(k),0,num2str(ANNOTD(k))); % end;
xlim([TIME(1), TIME(end)]);
xlabel('Time / s'); ylabel('Voltage / mV'); string=['ECG signal ',DATAFILE]; title(string);
fprintf(1,'\\\\n$> DISPLAYING DATA FINISHED \\n');
% ---------------------------------------------------------------- fprintf(1,'\\\\n$> ALL FINISHED \\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=[40 80]; f_s=[59 61]; 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',2 xlim([TIME(1), TIME(end)]);
xlabel('Time / s'); ylabel('Voltage / mV'); title('消除基线后ecg时序图');
%------ 确定波峰 -------------------------------------------- for i=1:length(TIME) if DATA(1,i)<0.8 DATA(1,i)=0; end end
figure(5);plot(TIME,DATA,'--o');xlim([0 80]);hold on; count=0;
for j=2:length(TIME)-1 if DATA(1,j)~=0
if (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); end end end
figure(5);plot(a(1,:),a(2,:),'*r');xlim([0 80]);ylim([0 2]); title('确定波峰位置');
%------ 确定R-R间距 ------------------------------------------------- b=diff(a(1,:));
figure(6);plot(b,'r','LineWidth',3);xlim([1 length(b)]);ylim([0 1.2]); grid on;
xlabel('R-Ri'); ylabel('R-R'); title('HRV图像');
%------ 时域分析 ---------------------------------------------------- NNVGR=mean(b) SDNN=std(b) c=abs(diff(b));
RNSSD=sqrt(sum(c.^2)/length(c)) SDSD=std(c) d=find(c>0.05); NN50=length(d) PNN50=NN50/length(c)
%------ 频域分析 ------------------------------------------------ b1=interp(b,100); %插值计算 fs=100;
[f,y]=fft_simple(b1,fs,0,50); %频谱 figure(7);
plot(f,y,'r','LineWidth',3); xlim([0 0.5]);
xlabel('Frequency / Hz'); ylabel('Amplitude'); title('HRV频谱图');
%------ 功率谱分析 -------------------------------------------- [Pxx,g]=periodogram(b1,window,128,fs); %功率谱 figure(9);
plot(g,10*log10(Pxx),'r','LineWidth',3)
xlabel('Frequency / Hz'); ylabel('Power / dB'); title('HRV功率谱图');
%------ 求总功率 ------------------------------------------------ ct0=0;
for l=1:length(f)
if f(1,l)<=0.4 ct0=ct0+1; Y1(1,ct0)=f(1,l); Y1(2,ct0)=y(1,l); end end
for l0=1:length(Y1(1,:)) T1(1,l0)=Y1(2,l0)*0.0104; end TP=sum(T1)
%------ 求出高频功率 ------------------------------------- ct1=0;
for l2=1:length(f)
if (f(1,l2)>=0.15)&&(f(1,l2)<0.4) ct1=ct1+1;
Y2(1,ct1)=f(1,l2); Y2(2,ct1)=y(1,l2); end end
for l3=1:length(Y2(1,:)) T2(1,l3)=Y2(2,l3)*0.0104; end HF=sum(T2)
%------ 求出低频功率 ----------------------------------------- ct2=0;
for l4=1:length(f)
if (f(1,l4)>=0.04)&&(f(1,l4)<0.15) ct2=ct2+1;
Y3(1,ct2)=f(1,l4); Y3(2,ct2)=y(1,l4); end end
for l5=1:length(Y3(1,:)) T3(1,l5)=Y3(2,l5)*0.0104; end
LF=sum(T3)
%------ 低频功率与高频功率之比 --------------------------------------- A=LF/HF
%====================================================================== % 文件名：fft_simple.m【新建脚本文件于上脚本同一路径】 % 功 能：函数，做短时傅里叶变换，计算频谱
% 输 入：一维数据x;采样频率srate;频率下限fl;频率上限fh % 输 出：频率上下限之间的频率序列及所对应频谱序列
% 备 注：返回值只包含频率上下限之间的值，需根据需求设置fl和fh %====================================================================== %% function
function [f,y] = fft_simple(x,srate,fl,fh) if fl<0
error('MATLAB:fft_simple:频率下限不能小于零') end
if fh>srate/2
error('MATLAB:fft_simple:频率上限不能大于采样频率的一半') end
N = length(x); %求输入数据的长度 z = abs(fft(x,N))*2/N; %短时傅里叶变换， p_s = 0:1:N/2; %频率序列对应的点序列 df = srate/N; %频率梯度，连续两点的频率差值
f_s = p_s*df; %产生频率序列，由点序列乘频率梯度产生 f = f_s(fl/df+1:1:fh/df+1); %由fl和fh产生频率序列 y = z(fl/df+1:1:fh/df+1); %由fl和fh产生频谱序列 end