CN108814579B - 一种基于emd分解的心电、呼吸联合计算心率变异性的方法 - Google Patents

一种基于emd分解的心电、呼吸联合计算心率变异性的方法 Download PDF

Info

Publication number
CN108814579B
CN108814579B CN201810338619.1A CN201810338619A CN108814579B CN 108814579 B CN108814579 B CN 108814579B CN 201810338619 A CN201810338619 A CN 201810338619A CN 108814579 B CN108814579 B CN 108814579B
Authority
CN
China
Prior art keywords
heart rate
respiration
wavelet
frequency
hrv
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810338619.1A
Other languages
English (en)
Other versions
CN108814579A (zh
Inventor
张建保
刘彬彬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201810338619.1A priority Critical patent/CN108814579B/zh
Publication of CN108814579A publication Critical patent/CN108814579A/zh
Application granted granted Critical
Publication of CN108814579B publication Critical patent/CN108814579B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
    • A61B5/024—Detecting, measuring or recording pulse rate or heart rate
    • A61B5/02405—Determining heart rate variability
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
    • A61B5/024—Detecting, measuring or recording pulse rate or heart rate
    • A61B5/0245—Detecting, measuring or recording pulse rate or heart rate by using sensing means generating electric signals, i.e. ECG signals
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7264—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Cardiology (AREA)
  • Physics & Mathematics (AREA)
  • Biophysics (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Physiology (AREA)
  • General Health & Medical Sciences (AREA)
  • Animal Behavior & Ethology (AREA)
  • Surgery (AREA)
  • Molecular Biology (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Signal Processing (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Fuzzy Systems (AREA)
  • Psychiatry (AREA)
  • Mathematical Physics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
  • Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)

Abstract

本发明公开了一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，先将心率做EMD分解，将其分为几层固有模态，以心率EMD分解所得的各层信号与相应呼吸斜率之间的相关性作为依据，把心率的各层模态划分为与呼吸相关的模态，以及与呼吸不相关的模态，分别作为HRV的高频段(HF)、和低频段(LF)的表征，选用了连续小波变换，可以获得时频HRV变化结果，有别于传统HRV分析方法中HF、LF固定频率分割点(f＝0.15Hz)所存在的问题(如HF谱峰会随着呼吸频率平移，甚至低于0.15Hz)，本方法从生理机制出发，将呼吸对HRV的影响作为参考，避免HRV分析领域中长久存在的实验结果的不一致性、难以复现性、以及实验结果解释的不科学性的问题。

Description

一种基于EMD分解的心电、呼吸联合计算心率变异性的方法
技术领域
本发明属于生物医学信号处理技术领域，旨在提出一种基于EMD分解的心电、呼吸联合计算心率变异性的方法。
背景技术
心脏是人体中极为重要的器官，它通过节律性收缩与舒张，推动血液在血管中按照一定的方向不停地循环流动，从而保证了机体内环境的相对恒定以及新陈代谢的正常进行。研究证实，心脏除了按照窦房结的基本节律跳动，同时也受到体液以及神经调节的作用(其中后者占主导)，以适应机体及外界环境变化。神经调节主要通过自主神经系统(交感和副交感神经系统)实现。
心脏自主神经功能的失衡与心血管疾病和死亡有相当大的关系。失衡导致房颤、室性心律失常、心力衰竭等恶心心血管事件，进而导致患者死亡。因而，心脏自主神经功能的评测，对于预防和减少恶性心血管事件，检验相关药物的疗效，具有重要的临床应用价值。
自主神经功能的众多评测方法中，心率变异性(HRV)以其易计算、无创性和定量性，已成为使用最为频繁的指标。HRV有时域、频域、非线性域等多种研究手段，其中频域方法更为广泛使用，因为它能分别获得交感和副交感神经对心脏的调节指标。在大量使用和研究中，HRV频域方法所存在的问题也逐渐显现。最初研究认为，频域中的两个成分：低频LF(0.03Hz～0.15Hz)和高频HF(0.15Hz～0.4Hz)分别对应交感和副交感的调节作用。但近年来，有研究提出LF中也存在副交感调控信息，认为LF可能是表征压力感受器作用的指标。HRV的可重复性有待提高，不同研究者在使用同样的实验范式时，常得出不一致的HRV结果，这也导致相同实验，会出现各种相互冲突的生理解释。
为了对传统HRV算法进行改进，有研究者尝试用时频等新的信号处理方法对HRV进行频谱估计。但这只是技术上的改进，并未从生理机制上对问题进行探讨。最初，将HRV分为LF和HF是因为，研究者对心率做频谱分析时发现，其频谱存在两个明显地相互独立的谱峰，于是划定0.15Hz作为两个谱峰之间的“分割点”，由此划定了LF和HF的概念。此后，关于HRV频谱的解释和分析都以此为基础进行。
但众多研究发现，呼吸通过副交感神经对心率有着重大影响，吸气时心率升高，呼吸时心率下降。实际上，呼吸频率会直接影响HRV的频谱形态。其HF峰值点与呼吸频率十分一致。极端的例子是，当呼吸频率低于0.15Hz时，HRV中第二个谱峰，即HF峰值点会移向0.15Hz以下。这使得长期以来使用的固定频率分割点的方法，出现了其问题。此种情况下，若仍以0.15Hz作为分频点，势必造成LF大幅升高，而HF降低的计算结果，即交感增强副交感减弱，这显然与事实不符。此外，当呼吸频率很高时(如>0.4Hz)，亦能观察到相应的HRV频谱能量在大于0.4Hz处广泛分布。因此，传统HRV中，对HF的计算终止于0.4Hz也存在问题。
另外，传统的HRV算法通常使用基于FFT的方法进行功率谱估计，而生理信号属于非平稳时变信号，而基于FFT的算法假定信号局部平稳，用固定的窗长对信号进行变换。并且由于测不准原理，FFT方法无法在时域和频域同时有好的表征能力。因此传统的HRV分析方法，计算了一段时间内的心率变异性总体上包含了哪些成分，但是无法描述对各成分随着时间的变化情况，导致的结果是时域相差很大的两个信号，可能频谱图一样。而在具体分析一段时间的生理信号时，我们除了要了解其频率成分，也想知道各成分的出现时间和变化过程，各时间的瞬时频率及幅值。因此，在最后的HRV频谱计算中，本发明选用连续小波变换，因为小波变换采用可伸缩的小波基函数，具有同时在时、频域表征信号局部特征的能力，更适于表征心率变异性所含的细节变化。
发明内容
本发明的目的在于提供一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，以克服现有技术的不足。
为达到上述目的，本发明采用如下技术方案：
一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，包括以下步骤：
步骤1)、同步采集被测试者的心电信号ecg(t)与呼吸信号res(t)；
步骤2)、对心电信号ecg(t)进行预处理得到心率各模态，对呼吸信号res(t)进行预处理得到呼吸斜率集合s(n)；
步骤3)、对步骤2)得到的心率各模态与呼吸斜率集合s(n)进行相关性计算，分别获得与呼吸显著相关的心率成分HRr和与呼吸不相关的心率成分HRur；
步骤4)、进行HRV计算：将步骤3)得到的与呼吸显著相关的心率成分HRr和与呼吸不相关的心率成分HRur分别除以平均心率；
选用cmor3-1母小波，对除去平均心率的HRr进行小波时频分析，所得的总功率作为HRV的高频成分HF；
选用cmor3-1母小波，对除去平均心率的HRur进行小波时频分析，所得的总功率作为HRV的低频成分LF。
进一步的，步骤1)中，信号采样率为fs＝1000Hz。
进一步的，步骤2)中，对采集到的原始心电信号ecg(t)，首先利用小波分解去除基线漂移及高频噪声，得到去噪后的心电信号x(t)；然后利用小波模极大值的方法，选用的小波为支持紧支集且具有一阶消失矩db2小波，对R波峰进行检测，再对检测结果进行检视修改，获得一连串连续R波峰时间点的集合R(N)；对R波峰时间点的集合R(N)，利用后一个波峰时间点R(n+1)减去前一个波峰的时间点R(n)，获得逐拍RR间期RR(n)，利用HR＝60/RR得到逐拍心率HR(n)，其中(1≤n<N)；
将逐拍心率HR(n)重采样到f＝2Hz，得到心率间期HR(m)(1≤m≤f*(b-a))，其中a、b为重采样的起止时间，满足R(1)≤a且b≤R(N-1)；
然后对心率间期HR(m)进行EMD分解，得到心率各模态，根据信号的局部时间特征尺度，按频率由高到低把复杂的非线性、非平稳信号心率间期HR(m)分解为I个本征模态函数IMFci(m)(1<i<I)及一个残差r(m)，即
进一步的，利用小波模极大值的方法对R波峰进行检测，采用支持紧支集且具有一阶消失矩的db2小波，db2小波的等效滤波器系数如下：
h1＝0.3750,h2＝0.1250,h3＝0.0000
g1＝0.5798.g2＝0.0869,g3＝0.0061
hk＝h1-k,gk＝-g1-k
若k>3,hk＝gk＝0
对所采集到的心电信号用Mallat算法作小波分解。
进一步的，对采集到的原始心电信号ecg(t)，选择db2小波对信号进行9层分解，去除最后一层的基线漂移以及前三层的高频噪声，得到去噪后的心电信号x(t)。
进一步的，利用步骤2)中获得的R波峰时间点的集合R(N)，求相邻R(n)到R(n+1)时间内对应的呼吸信号的斜率，得到与心率对应的呼吸斜率集合s(n)，其中1≤n<N；
s(n)＝(y(R(n+1))-y(R(n)))/(R(n+1)-R(n))。
进一步的，对s(n)重采样时，选取了与心率重采样完全相同的起止时间a、b，以及相同的采样频率f＝2Hz，得到s(m)(1≤m≤f*(b-a))。
进一步的，将步骤2)中所得的心率EMD分解各模态，即ci(m)和r(m)，与步骤2)中所得的呼吸斜率集合s(n)，分别计算其相关性大小及相关显著性P值，将显著性P＜0.01的各模态进行标注，并进行叠加，获得与呼吸显著相关的心率成分HRr；而将剩余层叠加，获得与呼吸不相关的心率成分HRur。
与现有技术相比，本发明具有以下有益的技术效果：
本发明一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，先将心率做EMD分解，将其分为几层固有模态，然后考虑呼吸对心率的影响，主要体现在斜率与心率之间的关系上，以心率EMD分解所得的各层信号与相应呼吸斜率之间的相关性作为依据，把心率的各层模态划分为与呼吸相关的模态，以及与呼吸不相关的模态，分别作为HRV的高频段(HF)、和低频段(LF)的表征，选用了连续小波变换，可以获得时频HRV变化结果，为科研和实际应用提供了更为细节的分析手段，是一种更为准确和科学的HRV计算方法。有别于传统HRV分析方法中HF、LF固定频率分割点(f＝0.15Hz)所存在的问题(如HF谱峰会随着呼吸频率平移，甚至低于0.15Hz)，本方法从生理机制出发，将呼吸对HRV的影响作为参考，避免HRV分析领域中长久存在的实验结果的不一致性、难以复现性、以及实验结果解释的不科学性的问题。
附图说明
图1为一种基于EMD分解的心电、呼吸联合计算心率变异性的方法的系统框图；
图2为R波峰点检测算法流程图；
图3为R波峰点检测检测结果图；
图4为重采样心率的EMD分解结果图；
图5为重采样心率与对应时间段内的呼吸斜率结果图；
图6为心率EMD分解前4层与呼吸斜率相关性示意图；
图7为与呼吸显著相关的心率成分HRr及与呼吸不相关的心率成分HRur示意图；
图8为传统HRV分析方法及本实验HRV分析方法结果高频对照图；
图9为传统HRV分析方法及本实验HRV分析方法结果低频对照图。
具体实施方式
下面结合附图对本发明做进一步详细描述：
如图1所示，一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，包括以下步骤：
步骤1)、同步采集被测试者的心电信号ecg(t)与呼吸信号res(t)；
步骤2)、对采集到的原始心电信号ecg(t)，首先利用小波分解去除基线漂移及高频噪声，得到去噪后的心电信号x(t)；然后利用小波模极大值的方法，选用的小波为支持紧支集且具有一阶消失矩db2小波，对R波峰进行检测，再对检测结果进行检视修改，获得一连串连续R波峰时间点的集合R(N)；对R波峰时间点的集合R(N)，利用后一个波峰时间点R(n+1)减去前一个波峰的时间点R(n)，获得逐拍RR间期RR(n)，利用HR＝60/RR得到逐拍心率HR(n)，其中(1≤n<N)；
将逐拍心率HR(n)重采样到f＝2Hz，得到心率心率间期HR(m)(1≤m≤f*(b-a))，其中a、b为重采样的起止时间，满足R(1)≤a且b≤R(N-1)；
然后对心率HR(m)进行EMD分解，得到心率各模态，根据信号的局部时间特征尺度，按频率由高到低把复杂的非线性、非平稳信号心率间期HR(m)分解为I个本征模态函数IMFci(m)(1<i<I)及一个残差r(m)，即
Figure GDA0002411889640000071
步骤3)、对采集到的呼吸信号res(t)，首先利用小波分解去除高频噪声，得到y(t)；利用步骤2)中获得的R波峰时间点的集合R(N)，求相邻R(n)到R(n+1)时间内对应的呼吸信号的斜率，得到与心率对应的呼吸斜率集合s(n)，其中(1≤n<N)。
s(n)＝(y(R(n+1))-y(R(n)))/(R(n+1)-R(n))
为了与心率完全对应，对s(n)重采样时，选取了与心率重采样完全相同的起止时间a、b，以及相同的采样频率f＝2Hz，得到s(m)(1≤m≤f*(b-a))。
如图5所示，重采样心率间期HR(m)与对应时间内呼吸的斜率s(n)，心率间期HR(m)中有明显地随着呼吸斜率进行变化的成分；
步骤4)、对步骤2)和步骤3)得到的心率各模态与呼吸斜率s(m)进行相关性计算：将步骤2)中所得的心率EMD分解各模态，即ci(m)和r(m)，与步骤(3)中所得的呼吸斜率集合s(n)，分别计算其相关性大小及相关显著性P值，将显著性P＜0.01的各模态进行标注，并进行叠加，获得与呼吸显著相关的心率成分HRr；而将剩余层叠加，获得与呼吸不相关的心率成分HRur。如图6所示，分别计算心率分解所得的各层模态与呼吸斜率在延时为正负20s以内的相关性，统计最大相关系数的值，计算相应延时时的相关显著性P值。依据显著性P＜0.01为标准，将满足标准的各层进行标注，并叠加，如图6中的第1层和第2层满足P＜0.01，因此将心率EMD分解中的1、2层进行叠加，获得与呼吸显著相关的心率成分HRr；同理，心率EMD分解所得的其它层不满足P＜0.01，则为与呼吸信号无相关层，将其叠加，获得与呼吸不相关的心率成分HRur，结果如图7所示。可以看到HRr包含了原始心率中，随着呼吸进程而改变的高频变化部分。而HRur则显示出于原始心率中，与呼吸关联不大的低频变化部分，并且HRr与HRur相加，正是原始心率信号。
步骤5)、进行HRV计算：将得到的与呼吸显著相关的心率成分HRr和与呼吸不相关的心率成分HRur分别除以平均心率；
选用cmor3-1母小波，对除去平均心率的HRr进行小波时频分析，所得的总功率作为HRV的高频成分HF；
选用cmor3-1母小波，对除去平均心率的HRur进行小波时频分析，所得的总功率作为HRV的低频成分LF。
具体的，步骤1)中，信号采样率为fs＝1000Hz；
具体的，步骤2)中，对采集到的原始心电信号ecg(t)，选择db2小波对信号进行9层分解，去除最后一层的基线漂移以及前三层的高频噪声，得到去噪后的心电信号x(t)；
对步骤2)利用小波模极大值的方法对R波峰进行检测，采用支持紧支集且具有一阶消失矩的db2小波，db2小波的等效滤波器系数如下：
h1＝0.3750,h2＝0.1250,h3＝0.0000
g1＝0.5798.g2＝0.0869,g3＝0.0061
hk＝h1-k,gk＝-g1-k
若k>3,hk＝gk＝0
对所采集到的心电信号用Mallat算法作小波分解；
1)特征尺度的选择：
小波变换将信号分解为不同频段的成分，高频分量和噪声主要落在小尺度上，低频分量和噪声主要落在大尺度上。对于不同的人，ECG信号中QRS波的频谱稍有不同，但能量主要集中在尺度23和24上，且在23上能量最大。高频分量较多的QRS波，22上能量最大，而低频分量较多的QRS波，24上能量最大。因此，本论文中选用21-24四个尺度。采用db2小波，把心电信号进行四层分解，得到小波变换成分W2jf(j＝1,2,3,4)和光滑信号S24f。
2)R波模极大值列的确定：
R波在每个特征尺度上都会生成一对小波变换模极大值，即正极大值-负极小值对。高频噪声会在尺度21和22上生成模极大值，而低频的高P波或高T波会在尺度23和24上生成模极大值，因此，检测这4个尺度上的模极大值列，以减小噪声以及高P波和高T波对R波检测的影响。
通过以下步骤确定R波对应的模极大值点：
(1)从尺度24的小波变换结果中找出大于阈值ε4的模极大值点，得到这些点的位置集合
Figure GDA0002411889640000091
(2)从尺度23的小波变换结果中，在
Figure GDA0002411889640000092
邻域内找出大于阈值ε3且与
Figure GDA0002411889640000093
处模极大值点同符号的模极大值点，本文中“邻域”选为
Figure GDA0002411889640000094
左右各10个点，将其位置定为
Figure GDA0002411889640000095
若
Figure GDA0002411889640000096
附近有几个模极大值点，则选最大的一个。若这个点的幅值小于1.2倍其它几个模极大值点的幅值，则选最靠近
Figure GDA0002411889640000097
的点。若
Figure GDA0002411889640000098
邻域内没有与
Figure GDA0002411889640000099
处模极大值点同符号的模极大值点，则令
Figure GDA00024118896400000910
为0；得到所有这些点位置的集合
Figure GDA00024118896400000911
(3)重复步骤(2)，找到尺度22、21上的模极大值点位置的集合
Figure GDA00024118896400000912
和
Figure GDA0002411889640000101
此算法中，不同的尺度采用的不同的幅度阈值{ε4，ε3，ε2，ε1}是根据最新检测到的小波变换模极大值来刷新的，刷新公式为：
若，
Figure GDA0002411889640000102
则，
Figure GDA0002411889640000103
否则
Figure GDA0002411889640000104
而
Figure GDA0002411889640000105
其中
Figure GDA0002411889640000106
代表检测到的小波变换模极大值。这个方法可以保证，检测的QRS波幅值突然增大时不会影响以后的幅值判别。
使用的这种从大尺度向小尺度方向搜索小波变换模极大值点的方法可以去除小尺度上由噪声产生的模极大值点，准确地检测R波峰点，同时可以节约运算时间。
3)奇异点奇异度的计算：
令aj(nk)＝|W2jf(nk)|，假设a为Lipschitz指数的上限。令a≈log2 aj+1(nk)-log2aj(nk)。通过四层小波分解，可以得到a1、a2和a3。在R波峰点处必定有a1＞0，且通常情况下a2＞0，而且即使a2＜0时，a1+a2必定会大于0。对于大多数R波，通常有a3＜0而且a1+a2+a3≤0，而对于高频噪声和干扰剧烈，a1≤0，a2≤0，a3≤0，且a1+a2+a3≤0。因此，从a1+a2+a3的值不能分辨R波、高频噪声和干扰，而a1+a2有很好的分辨效果。所以，在检测R波时选用了a1、a2，并令a'＝(a1+a2)/2。若a'>0，则相应的模极大值点是R波峰值点所对应的；若a'突然减小变为负值，则相应的模极大值点一定是噪声或干扰所对应的，应删除相应的模极大值列。
4)去除孤立的和多余的模极大值列：
运动伪迹和肌电噪声的频带通常与QRS波的频带重叠。因此，应从前面所得到的模极大值列集合中，剔除由伪迹或肌电噪声引入的模极大值列。
(1)删除孤立的模极大值列。
R波在每个特征尺度上都对应于一对模极大值列，即正极大值-负极小值对。这两个模极大值点的间距在尺度21上比R波的宽度小。设
Figure GDA0002411889640000111
为尺度21上
Figure GDA0002411889640000112
的一个正极大值点，
Figure GDA0002411889640000113
为同一尺度上
Figure GDA0002411889640000114
的负极小值点，若
Figure GDA0002411889640000115
和
Figure GDA0002411889640000116
间距大于给定的阈值，设间距阈值为120ms，则称
Figure GDA0002411889640000117
为孤立模极大值点，应将其删除。
(2)删除多余的模极大值列。
一个R波在每个尺度上只产生一对模极大值点。但一些带噪声的R波，在一对正极大值-负极小值对的邻域(阈值为120ms)内，会产生两对或更多的模极大值列，而其中仅有一对是有用的。多余的模极大值列可以使用下面的准则来删除：
这里选择QRS波的能量主要集中的尺度23上的模极大值来判别。设两个负极小值点分别为Min1和Min2，其幅值分别为A1和A2，而它们与相应的正极大值点的距离分别为L1和L2。
准则1：若
Figure GDA0002411889640000118
则Min2是多余模极小值点；
准则2：若
Figure GDA0002411889640000119
则Min1是多余模极小值点；
准则3：否则，若Min1，Min2在该正极大值的同侧，那么离该正极大值远的是多余模极大值点；若Min1，Min2在该正极大值点的两侧，那么该正极大值点后的那一点为多余点。
5)R波峰点检测：
前面算法已经剔除了噪声和干扰以及P波T波对应的正极大值-负极小值对，并且删除了孤立的和多余的模极大值列，获得了尺度21上ECG信号R波峰点对应的的正极大值-负极小值对；检测出这些正极大值-负极大值对的过零点的位置，就得到了R波峰点的位置；
为了提高检测率，还运用了一下两条策略：
1)不应期：
一般人的心率小于300次/分钟，在一次心跳过后的一段时间内不会产生另一次心跳，也就是会有一段不应期；因此，在检测到一个R波之后，把其后200ms内的极值都忽略，可以避免因噪声干扰而造成的误检。
2)反向搜索：
在心率失常或其他情况下，R波幅度和频率会突然变小，导致模极大值点的幅度达不到阈值，导致漏检；在我们的算法中，先对前30秒内所检测的RR间期进行平均，得到最近一段时间的平均心动周期T，若本次检测的RR间期大于1.5T，则在此间期内在尺度23上用0.5ε3检测模极大值。若此区间内的一对正极大值-负极小值对点之间的间隔小于140ms，则认为有漏检，检测它们之间的零交叉点，并用3点的时移修正，得到重检的R波。采用这种方法可以减少绝大多数情况下的漏检。
检测算法流程如图2所示。
对以上方法检测结果，最后再进行检视并修改，保证无误后，获得一连串连续R波峰时间点的集合R(N)。其检测结果如图3所示：
对R波峰时间点的集合R(N)，利用后一个波峰时间点R(n+1)减去前一个波峰的时间点R(n)，获得逐拍RR间期RR(n)。HR＝60/RR得到逐拍心率HR(n)，其中(1≤n＜N)。
将心率HR(n)重采样到f＝2Hz，得到HR(m)(1≤m≤f*(b-a))，其中a、b为重采样的起止时间，必须满足：R(1)≤a且b≤R(N-1)。接着对HR(m)进行EMD分解，根据信号的局部时间特征尺度，按频率由高到低把复杂的非线性、非平稳信号心率间期HR(m)分解为I个经验模态函数(IMF)ci(m)(1＜i＜I)及一个残差r(m)。即：
Figure GDA0002411889640000131
其EMD分解结果如图4所示。可以看到EMD分解，根据信号本身的特性，既可以将信号按照其内部固有的频率成分做划分，也能避免误差的引入，是一种信号自适应型的算法。
计算HRV之前，首先要将HRr和HRur分别除以平均心率，这样可以消除平均心率对HRV计算结果的影响。
然后利用小波分解去除0.03Hz以下的频率成分。
在算法的选择上，不同传统的基于FFT的算法，本发明选用了基于连续小波变换时频分析方法。因为生理信号属于非平稳时变信号，而基于FFT的算法假定信号局部平稳，用固定的窗长对信号进行变换。因此传统的HRV分析方法，计算了一段时间内的心率变异性总体上包含了哪些成分，但是无法描述对各成分随着时间的变化情况，导致的结果是时域相差很大的两个信号，可能频谱图一样。而在具体分析一段时间的生理信号时，我们除了要了解其频率成分，也想知道各成分的出现时间和变化过程，各时间的瞬时频率及幅值。而小波变换选用的是可伸缩的小波基函数，具有同时在时、频域表征信号局部特征的能力，更适于表征心率变异性所含的细节变化。
对时间序列做分析时，期望能得到的小波振幅是平滑而连续的，因此选非正交小波函数比较合适。另外，要求得到振幅和相位两个信息，就要选择复小波。因为复小波偶虚部，可以很好的表达相位信息。Morlet小波不但有非正交性且是由Gaussian调节的指数复小波，因此本例做时频分析选取的是cmor3-1母小波。
母小波选定好以后，将信号x(t)进行MWT变换
对任意函数f(t)＝L2(R)，其连续小波变换为
Figure GDA0002411889640000141
a表示伸缩因子；而b则表示平移因子。得到时频解析信号y(t,fn)。则时频能量谱可以定义为：
ps(t,f)＝|y(t,f)|2
根据上述算法，对除去平均心率的HRr进行小波时频分析，所得的总功率作为HRV的高频成分HF；
对除去平均心率的HRur进行小波时频分析，所得的总功率作为HRV的低频成分LF；
如图8所示，(a)为对HR/mean(HR)做时频变换的结果，(b)和(c)分别是对HRr/mean(HR)和HRur/mean(HR)做时频变换的结果。可以看到，基于小波变换的时频HRV方法，在时间和频率上都有很好的分辨率。另外如果按传统方法f＝0.15Hz作为分割，图(a)中的大部分信息都会被分割到低频部分。而本发明的方法，如图(b)和(c)所示，分别对应LF和HF，由于以生理机制为划分依据，其划分结果由图可以明显看出，两张图的频率成分有明显区别，的确是更为科学和准确的方法。
对频谱图沿着频率方向积分，可获得总能量随时间的变化情况。举例说明，我们根据传统定义对时频结果沿着频率积分，低频段(0.03Hz～0.15Hz)和高频段(0.15Hz～0.04Hz)，即可获得低、高频段随时间变化情况。
Figure GDA0002411889640000151
Figure GDA0002411889640000152
同理，对频谱图沿着时间做积分，可得一段时间内，各频率点的能量变化情况。
为了与传统HRV计算结果做比较，如图9所示，对图8中的三个时频图，沿时间做叠加平均，绘制了传统HRV分析方法及本发明HRV分析方法结果对照图。
本次举例说明中，特别选用了一名呼吸频率低于0.15Hz的被试数据。可以看到，如果用传统方法，即图9(a)中的结果，直接将f＝0.15Hz作为频率划分点，那么与呼吸相关的副交感活动，即HRV中本应统计为HF成分的频谱成分，实则被划分到了LF部分。造成的统计结果是，副交感活性很小，而交感活性很大。这显然与事实不符。
但应用本发明所述的方法，如图9(b)、(c)所示，将呼吸信号作为参考，先将原始心率做EMD分解，在没有引入误差的情况下，根据心率自身特性将其进行自适应的分解。然后利用呼吸斜率与心率之间相关性，对分解结果进行划分，得出与呼吸相关HRr和不相关HRur的部分。最后把上述两个部分分别作为高(HF)和低(LF)的表征，进行频谱分析。看到，本算法中的LF(图9(b))和HF(图9(c))的频率成分存在极少部分的交叠。但其总和仍然与传统算法的结果一致，但其高低、频的划分更具有生理依据，是一种更为科学和有效的方法。应用此方法，可以避免HRV分析领域中长久存在的一些问题，诸如，实验结果的不一致性，难以复现性，以及实验结果解释的不科学性。
进一步，还可得到标准化低频(nLF)为nLF＝LFP/(LFP+HFP)，标准化高频(nHF)为nHF＝HFP/(LFP+HFP)，以及交感迷走平衡PR＝LFP/HFP。

Claims (8)

1.一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，包括以下步骤：
步骤1)、同步采集被测试者的心电信号ecg(t)与呼吸信号res(t)；
步骤2)、对心电信号ecg(t)进行预处理得到心率各模态，对呼吸信号res(t)进行预处理得到呼吸斜率集合s(n)；
步骤3)、对步骤2)得到的心率各模态与呼吸斜率集合s(n)进行相关性计算，分别获得与呼吸显著相关的心率成分HRr和与呼吸不相关的心率成分HRur；
步骤4)、进行HRV计算：将步骤3)得到的与呼吸显著相关的心率成分HRr和与呼吸不相关的心率成分HRur分别除以平均心率；
选用cmor3-1母小波，对除去平均心率的HRr进行小波时频分析，所得的总功率作为HRV的高频成分HF；
选用cmor3-1母小波，对除去平均心率的HRur进行小波时频分析，所得的总功率作为HRV的低频成分LF。
2.根据权利要求1所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，步骤1)中，信号采样率为fs＝1000Hz。
3.根据权利要求1所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，步骤2)中，对采集到的原始心电信号ecg(t)，首先利用小波分解去除基线漂移及高频噪声，得到去噪后的心电信号x(t)；然后利用小波模极大值的方法，选用的小波为支持紧支集且具有一阶消失矩db2小波，对R波峰进行检测，再对检测结果进行检视修改，获得一连串连续R波峰时间点的集合R(N)；对R波峰时间点的集合R(N)，利用后一个波峰时间点R(n+1)减去前一个波峰的时间点R(n)，获得逐拍RR间期RR(n)，利用HR＝60/RR得到逐拍心率HR(n)，其中1≤n<N；
将逐拍心率HR(n)重采样到f＝2Hz，得到心率间期HR(m)(1≤m≤f*(b-a))，其中a、b为重采样的起止时间，满足R(1)≤a且b≤R(N-1)；
然后对心率间期HR(m)进行经验模态分解EMD分解，得到心率各模态，根据信号的局部时间特征尺度，按频率由高到低把复杂的非线性、非平稳信号心率间期HR(m)分解为I个本征模态函数IMFci(m)(1<i<I)及一个残差r(m)，即
Figure FDA0002411889630000021
4.根据权利要求3所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，利用小波模极大值的方法对R波峰进行检测，采用支持紧支集且具有一阶消失矩的db2小波，db2小波的等效滤波器系数如下：
h1＝0.3750,h2＝0.1250,h3＝0.0000
g1＝0.5798.g2＝0.0869,g3＝0.0061
hk＝h1-k,gk＝-g1-k
若k>3,hk＝gk＝0
对所采集到的心电信号用Mallat算法作小波分解。
5.根据权利要求3所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，对采集到的原始心电信号ecg(t)，选择db2小波对信号进行9层分解，去除最后一层的基线漂移以及前三层的高频噪声，得到去噪后的心电信号x(t)。
6.根据权利要求3所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，利用步骤2)中获得的R波峰时间点的集合R(N)，求相邻R(n)到R(n+1)时间内对应的呼吸信号的斜率，得到与心率对应的呼吸斜率集合s(n)，其中1≤n<N；
s(n)＝(y(R(n+1))-y(R(n)))/(R(n+1)-R(n))。
7.根据权利要求6所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，对s(n)重采样时，选取了与心率重采样完全相同的起止时间a、b，以及相同的采样频率f＝2Hz，得到s(m)，1≤m≤f*(b-a)。
8.根据权利要求1所述的一种基于EMD分解的心电、呼吸联合计算心率变异性的方法，其特征在于，将步骤2)中所得的心率EMD分解各模态，即ci(m)和r(m)，与步骤2)中所得的呼吸斜率集合s(n)，分别计算其相关性大小及相关显著性P值，将显著性P＜0.01的各模态进行标注，并进行叠加，获得与呼吸显著相关的心率成分HRr；而将剩余层叠加，获得与呼吸不相关的心率成分HRur。
CN201810338619.1A 2018-04-16 2018-04-16 一种基于emd分解的心电、呼吸联合计算心率变异性的方法 Active CN108814579B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810338619.1A CN108814579B (zh) 2018-04-16 2018-04-16 一种基于emd分解的心电、呼吸联合计算心率变异性的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810338619.1A CN108814579B (zh) 2018-04-16 2018-04-16 一种基于emd分解的心电、呼吸联合计算心率变异性的方法

Publications (2)

Publication Number Publication Date
CN108814579A CN108814579A (zh) 2018-11-16
CN108814579B true CN108814579B (zh) 2020-05-22

Family

ID=64154616

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810338619.1A Active CN108814579B (zh) 2018-04-16 2018-04-16 一种基于emd分解的心电、呼吸联合计算心率变异性的方法

Country Status (1)

Country Link
CN (1) CN108814579B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109620179A (zh) * 2018-12-26 2019-04-16 南京茂森电子技术有限公司 一种呼吸对心率调制的定量分析方法及相关系统
CN110051325A (zh) * 2019-03-29 2019-07-26 重庆邮电大学 基于小波变换及改进eemd的心电信号综合滤波方法
CN110464337B (zh) * 2019-09-06 2020-09-01 江苏华康信息技术有限公司 一种基于极值能量分解法的心率变异性信号分析方法
CN110558959B (zh) * 2019-09-06 2020-09-01 江苏华康信息技术有限公司 一种基于极值能量分解法的冥想训练的hrv信号分析方法
CN110558974B (zh) * 2019-09-06 2020-11-03 江苏华康信息技术有限公司 一种基于极值能量分解法的心电图信号分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103876733A (zh) * 2014-03-12 2014-06-25 西安交通大学 用于心肺系统相位同步分析的系统和方法
CN104545888A (zh) * 2014-12-27 2015-04-29 迪姆软件(北京)有限公司 一种基于动态心电与呼吸波采集的睡眠呼吸暂停采集分析系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6993378B2 (en) * 2001-06-25 2006-01-31 Science Applications International Corporation Identification by analysis of physiometric variation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103876733A (zh) * 2014-03-12 2014-06-25 西安交通大学 用于心肺系统相位同步分析的系统和方法
CN104545888A (zh) * 2014-12-27 2015-04-29 迪姆软件(北京)有限公司 一种基于动态心电与呼吸波采集的睡眠呼吸暂停采集分析系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于RQA方法的心肺系统非线性耦合测量;刘一辉等;《心脏杂志》;20141231;第16卷(第2期);118-120 *

Also Published As

Publication number Publication date
CN108814579A (zh) 2018-11-16

Similar Documents

Publication Publication Date Title
CN108814579B (zh) 一种基于emd分解的心电、呼吸联合计算心率变异性的方法
Lin et al. Comparison of heart rate variability from PPG with that from ECG
US7809433B2 (en) Method and system for limiting interference in electroencephalographic signals
US7623912B2 (en) Method, apparatus and system for characterizing sleep
Alcaraz et al. Adaptive singular value cancelation of ventricular activity in single-lead atrial fibrillation electrocardiograms
Pueyo et al. QRS slopes for detection and characterization of myocardial ischemia
US9125577B2 (en) Extraction of fetal cardiac signals
Willson et al. Relationship between detrended fluctuation analysis and spectral analysis of heart-rate variability
US10912479B2 (en) Method for accurately extracting abnormal potential within QRS
Jeyarani et al. Analysis of noise reduction techniques on QRS ECG waveform-by applying different filters
Singh et al. Effects of RR segment duration on HRV spectrum estimation
Ranjan et al. Cardiac artifact noise removal from sleep EEG signals using hybrid denoising model
Burri et al. Wavelet transform for analysis of heart rate variability preceding ventricular arrhythmias in patients with ischemic heart disease
Colak Preprocessing effects in time–frequency distributions and spectral analysis of heart rate variability
Mendez et al. Automatic detection of sleep macrostructure based on bed sensors
Bojorges-Valdez et al. Scaling patterns of heart rate variability data
Liu et al. Detection of myocardial infarction from multi-lead ECG using dual-Q tunable Q-factor wavelet transform
Lemay et al. Phase-rectified signal averaging used to estimate the dominant frequencies in ECG signals during atrial fibrillation
CN109044335B (zh) 一种基于瞬时声音刺激的心脏功能评价方法
Chen et al. A fast ECG diagnosis using frequency-based compressive neural network
Song et al. A Robust and Efficient Algorithm for St–T Complex Detection in Electrocardiograms
Joshi et al. Arterial pulse rate variability analysis for diagnoses
Antsiperov et al. A new PVC/SPB detection method: Based on analytical spectra technique
CN115886742B (zh) 一种基于变分心肺耦合的睡眠动力学分析方法
Zakaria et al. Heart rate variability (HRV) analysis using DSP for the detection of myocardial infarction

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant