CN110558959B - 一种基于极值能量分解法的冥想训练的hrv信号分析方法 - Google Patents

一种基于极值能量分解法的冥想训练的hrv信号分析方法 Download PDF

Info

Publication number
CN110558959B
CN110558959B CN201910841609.4A CN201910841609A CN110558959B CN 110558959 B CN110558959 B CN 110558959B CN 201910841609 A CN201910841609 A CN 201910841609A CN 110558959 B CN110558959 B CN 110558959B
Authority
CN
China
Prior art keywords
signal
extreme
meditation
training
energy
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
CN201910841609.4A
Other languages
English (en)
Other versions
CN110558959A (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.)
Jiangsu Huakang Information Technology Co ltd
Original Assignee
Jiangsu Huakang Information Technology Co ltd
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 Jiangsu Huakang Information Technology Co ltd filed Critical Jiangsu Huakang Information Technology Co ltd
Priority to CN201910841609.4A priority Critical patent/CN110558959B/zh
Priority to PCT/CN2019/121417 priority patent/WO2021042592A1/zh
Publication of CN110558959A publication Critical patent/CN110558959A/zh
Application granted granted Critical
Publication of CN110558959B publication Critical patent/CN110558959B/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/0205—Simultaneously evaluating both cardiovascular conditions and different types of body conditions, e.g. heart and respiratory condition
    • 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/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
    • 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/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • 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
    • 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/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • Cardiology (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Physiology (AREA)
  • Signal Processing (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Psychiatry (AREA)
  • Pulmonology (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明公开了一种基于极值能量分解法的冥想训练的HRV信号分析方法，包括获取给定时间和给定采样频率下的冥想训练中的ECG信号，去噪后得到RRI信号；将RRI信号x(t)作为原始信号，将原始信号x(t)分解为n个极值模态函数分量和一个余量，将原始信号x(t)分解得的n个极值模态函数分量，代表了原始信号不同频段的分量，根据n个极值模态函数分量判定该RRI信号是否进入冥想状态。本发明采用极值能量分解方法分析RRI信号，将原始信号分解为多个分量，也就是极值分量函数，计算每一个分量的能量，得到其能量分布。

Description

一种基于极值能量分解法的冥想训练的HRV信号分析方法
技术领域
本发明涉及一种心电图信号分析，尤其涉及一种基于极值能量分解法的冥想训练的HRV信号分析方法。
背景技术
生理信号是由生命体多个系统相互作用产生的，不同系统作用的时间和强度不同，导致生理信号具有时间和空间上的复杂性。心率变异性(HRV)是指测量连续心动周期之间的时间变异数，准确地说，应该是测量连续出现的正常P-P间期之间的差异的变异数。但由于P波不如R波明显或P波顶端有时宽钝，所以对心率变异性信号的研究通常用与P-P间期相等的R-R间期信号(RRI)来代替。研究表明，HRV可做为植物神经系统活动的无创性检测指标，尤其在判断某些心血管疾病的预后方面有着重要意义。
现有技术中研究心脏的长时调节节律(<1Hz))常用心律变异性信号(HRV信号)作为分析对象。大量的研究表明，人体HRV信号具有长时相关性和非线性动力学复杂性，并且年龄和疾病会导致动力学复杂性降低。对心率变异性(Heart Rate Varibility,HRV)信号的研究常用的是RR间期(Interbeat Intervals，RRI)信号，即连续RRI信号R波之间的时间间隔信号。
研究HRV信号的能量改变最常用的方法是功率谱分析(PSD)。PSD通过傅里叶变换将HRV信号的功率转化成频率的函数，研究不同频域范围的功率大小，通常HRV频谱分为高频(HF)、低频(LF)和极低频(VLF)等几个频段,频域分割的典型方式为：HF(0.15～0.4Hz)，LF(0.04～0.15Hz)，VLF(0.0033～0.04Hz)。LF/HF比值有重要的临床价值。心脏疾病会引起HRV功率谱的改变，比如心衰和心肌梗死引起归一化HF增高、LF和VLF降低。然而，PSD不是一种基于数据驱动的方法，并且对频域的分段比较粗糙，导致细节缺失，分割也不够灵活。
大量的研究表明，人体心率变异性(Heart Rate Varibility,HRV)信号具有长时相关性和非线性动力学复杂性，并且年龄和疾病会导致动力学复杂性降低，另外心理因素也会影响HRV的波动特征。冥想训练是一种通过集中注意力来调节呼吸、心跳等的训练方法，通过冥想训练，人们可以调节交感和副交感神经系统和心血管系统，达到放松身心、减轻心理压力的目的。但冥想训练对于HRV信号的影响
因此，亟待解决上述问题。
发明内容
发明目的：本发明的目的是提供一种可采用较少数据即可从HRV的极值能量分布角度分析其在不同分解层次上的波动强度，根据冥想训练对于HRV信号的影响判断是否进入冥想状态的基于极值能量分解法的冥想训练的HRV信号分析方法。
技术方案：为实现以上目的，本发明公开了一种基于极值能量分解法的冥想训练的HRV信号分析方法，包括如下步骤：
(1)、获取给定时间和给定采样频率下的冥想训练中的ECG信号，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练中RRI信号x(t)；
(2)、将RRI信号x(t)作为原始信号，求出原始信号的所有局部极值点，然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin，得到上、下包络线的包络均值信号m(t)＝(emax+emin)/2；
(3)、将原始信号x(t)减去包络均值信号m(t)，得到h(t)＝x(t)-m(t)；然后判断h(t)是否满足极值模态函数的判定条件，如果不满足，将h(t)作为原始信号返回至步骤(2)，直到hk(t)满足极值模态函数的判定条件，则记c1(t)＝hk(t)，作为第一个极值模态函数分量；
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t)，得到余量r1(t)＝x(t)-c1(t)，然后判断hk(t)是否满足停止准则，如果不满足，将r1(t)作为新的原始序列x(t)，返回至步骤(2)和(3)；如果hk(t)满足停止准则但n＜8时，返回步骤(1)重新获取原始信号；如果hk(t)满足停止准则且n≥8时，得到第2、3、…、n个极值模态函数分量及余量rn(t)，于是将原始信号x(t)分解为n个极值模态函数分量和一个余量，即
Figure GDA0002235838320000021
(5)、将原始信号x(t)分解得的n个极值模态函数分量，代表了原始信号不同频段的分量，然后计算其各个分量的能量
Ei＝∫|ci(t)|2dt,i＝1，2，…，n
将每一个能量值归一化，得到归一化的能量分布向量
pi＝Ei/E,i＝1,2,…,n
其中，
Figure GDA0002235838320000022
第一个分量p1表示最高频段的能量，代表了信号在最高频段范围内能量分布的比例，最后一个分量pn表示信号在最低频段范围内能量分布的比例；根据归一化的能量分布向量绘制归一化能量分布图，其中横坐标表示分量层次，纵坐标表示归一化的能量分布向量值，曲线表示平均值，误差棒表示标准差；
(6)、分别选取第二个分量p2至第七个分量p7，计算冥想训练中RRI信号的能量差异值EDV，EDV＝(p2+p3+p4)-(p5+p6+p7)，当冥想训练中RRI信号的EDV大于EDV标准值，则判定该冥想训练中的RRI信号为已进入冥想状态的RRI信号。
其中，所述原始信号x(t)所需最少数据量N＝2n+1，其中n为分解出的极值模态函数分量的数量。
再者，所述步骤(1)中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移。
进一步，所述步骤(3)中极值模态函数的判定条件为：(a)、在整个数据序列中，极值点的数量与过零点的数量相等或者相差一个；(b)、在任意时刻，上下包络线对于时间轴对称。
优选的，所述步骤(4)中hk(t)满足停止准则公式为：
Figure GDA0002235838320000031
ε表示筛选门限，取0.2～0.3之间。
再者，所述冥想训练为太极训练，其步骤(6)中的标准值为0.26。
进一步，所述冥想训练为瑜伽训练，其步骤(6)中的标准值为0.68。
有益效果：与现有技术相比，本发明具有以下显著优点：
本发明采用极值能量分解方法(Extremum Energy Decomposition,EED)分析RRI信号，将原始信号分解为多个分量，也就是极值分量函数，计算每一个分量的能量，得到其能量分布；本发明的可依据RRI信号自身的波动规律将信号分解为从高频到低频的不同时间层次信号，对频段的分割较为细致；冥想训练对于RRI信号的影响很敏感，通过实验与研究，冥想训练，打破了其训练前的能量分布，显示出与随机序列一样的特征，说明冥想训练导致RRI信号的长时相关性消失；本发明极值分解在所有层次上得到的数据长度相同，因而不会导致数据长度减小，从而使其可以用于短时间数据分析，即需要很少数据量即可分析得到准确结果；EED对于不同层次分量能量分析不容易受噪声的影响。
附图说明
图1为本发明中原始信号的示意图；
图2为本发明中原始信号求取包络线的示意图；
图3为本发明中原始信号的减去包络均值信号的示意图；
图4为本发明中得到第一个极值模态函数分量的示意图；
图5为本发明中极值能量分解法的流程示意图；
图6为本发明实施例1中太极冥想训练前和训练中的RRI信号的示意图；
图7为本发明实施例1中瑜伽冥想训练前和训练中的RRI信号的示意图；
图8为本发明实施例2中太极冥想训练的RRI信号的EED分解示意图；
图9为本发明实施例2中太极冥想训练的RRI信号的EED分解示意图。
具体实施方式
下面结合附图对本发明的技术方案作进一步说明。
如图1、图2、图3、图4和图5所示，本发明一种基于极值能量分解法的冥想训练的HRV信号分析方法，包括如下步骤：
(1)、获取给定时间和给定采样频率下的冥想训练中的ECG信号，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练中RRI信号x(t)；其中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移；
(2)、将RRI信号x(t)作为原始信号，原始信号x(t)所需最少数据量N＝2n+1，其中n为分解出的极值模态函数分量的数量，求出原始信号的所有局部极值点，然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin，得到上、下包络线的包络均值信号m(t)＝(emax+emin)/2；
(3)、将原始信号x(t)减去包络均值信号m(t)，得到h(t)＝x(t)-m(t)；然后判断h(t)是否满足极值模态函数的判定条件，如果不满足，将h(t)作为原始信号返回至步骤(2)，直到hk(t)满足极值模态函数的判定条件，则记c1(t)＝hk(t)，作为第一个极值模态函数分量；其中极值模态函数的判定条件为：(a)、在整个数据序列中，极值点的数量与过零点的数量相等或者相差一个；(b)、在任意时刻，上下包络线对于时间轴对称；
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t)，得到余量r1(t)＝x(t)-c1(t)，然后判断hk(t)是否满足停止准则，如果不满足，将r1(t)作为新的原始序列x(t)，返回至步骤(2)和(3)；如果hk(t)满足停止准则但n＜8时，返回步骤(1)重新获取原始信号；如果hk(t)满足停止准则且n≥8时，得到第2、3、…、n个极值模态函数分量及余量rn(t)，于是将原始信号x(t)分解为n个极值模态函数分量和一个余量，即
Figure GDA0002235838320000051
其中hk(t)满足停止准则公式为：
Figure GDA0002235838320000052
ε表示筛选门限，取0.2～0.3之间；
(5)、将原始信号x(t)分解得的n个极值模态函数分量，代表了原始信号不同频段的分量，然后计算其各个分量的能量
Ei＝∫|ci(t)|2dt,i＝1，2，…，n
将每一个能量值归一化，得到归一化的能量分布向量
pi＝Ei/E,i＝1,2,…,n
其中，
Figure GDA0002235838320000053
第一个分量p1表示最高频段的能量，代表了信号在最高频段范围内能量分布的比例，最后一个分量pn表示信号在最低频段范围内能量分布的比例；根据归一化的能量分布向量绘制归一化能量分布图，其中横坐标表示分量层次，纵坐标表示归一化的能量分布向量值，曲线表示平均值，误差棒表示标准差；
(6)、分别选取第二个分量p2至第七个分量p7，计算冥想训练中RRI信号的能量差异值EDV，EDV＝(p2+p3+p4)-(p5+p6+p7)，当冥想训练中RRI信号的EDV大于EDV标准值，则判定该冥想训练中的RRI信号为已进入冥想状态的RRI信号；当冥想训练为太极训练，其步骤(6)中的标准值为0.26；当冥想训练为瑜伽训练，其步骤(6)中的标准值为0.68。
本发明采用的极值能量分解方法(Extremum Energy Decomposition,EED)，其是基于极值模态函数的概念的方法，极值模态函数是同时满足下面两个条件的具有单一频率的一类信号，两个条件为：
(a)、在整个数据序列中，极值点(包括极大值和极小值)的数量与过零点的数量必须相等或者最多相差一个；
(b)在任意时刻，局部极大值点形成的上包络线与局部极小值点形成的下包络线的均值为零，也就是说局部上下包络线对于时间轴局部对称；
上面两个条件，条件(a)类似于高斯正态平稳过程对于传统窄带的要求，条件(b)保证了由极值模态函数计算得到的瞬时频率有物理意义。
本发明的极值模态函数分解终止的标准选择要适中，条件太严格，会导致最后几个极值模态函数分量失去意义；条件太宽松，会导致有用的分量丢失；在实际应用中，也可以根据需求设定需要分解的极值模态函数分量层数，当满足分解层数时即终止计算。
实施例1
太极(Chi)和瑜伽(Yoga)分别是流行于中国和印度的冥想训练方法，本发明将分别采用太极和瑜伽的冥想训练数据进行HRV信号分析方法，探寻心理调节改变HRV能量分布的规律。
太极训练受试者的一种基于极值能量分解法的冥想训练的HRV信号分析方法，包括如下步骤：
(1)、如图6所示，从physionet的meditation数据库获取给定时间和给定采样频率下的太极训练受试者冥想训练前的ECG信号，该数据库包含8个太极训练受试者(年龄26～35,平均29)，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练前RRI信号x(t)；其中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移；
(2)、将去噪后的RRI信号x(t)作为原始信号，其中原始信号x(t)所需最少数据量N＝2n+1，其中n为分解出的极值模态函数分量的数量；求出原始信号的所有局部极值点，然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin，得到上、下包络线的包络均值信号m(t)＝(emax+emin)/2；
(3)、将原始信号x(t)减去包络均值信号m(t)，得到h(t)＝x(t)-m(t)；然后判断h(t)是否满足极值模态函数的判定条件，如果不满足，将h(t)作为原始信号返回至步骤(2)，直到hk(t)满足极值模态函数的判定条件，则记c1(t)＝hk(t)，作为第一个极值模态函数分量；其中极值模态函数的判定条件为：(a)、在整个数据序列中，极值点的数量与过零点的数量相等或者相差一个；(b)、在任意时刻，上下包络线对于时间轴对称；
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t)，得到余量r1(t)＝x(t)-c1(t)，然后判断hk(t)是否满足停止准则，如果不满足，将r1(t)作为新的原始序列x(t)，返回至步骤(2)和(3)；如果hk(t)满足停止准则但n＜8时，返回步骤(1)重新获取原始信号；如果hk(t)满足停止准则且n≥8时，得到第2、3、…、n个极值模态函数分量及余量rn(t)，于是将原始信号x(t)分解为n个极值模态函数分量和一个余量，即
Figure GDA0002235838320000071
其中hk(t)满足停止准则公式为：
Figure GDA0002235838320000072
ε表示筛选门限，取0.2～0.3之间；
(5)、将原始信号x(t)分解得的n个极值模态函数分量，代表了原始信号不同频段的分量，然后计算其各个分量的能量
Ei＝∫|ci(t)|2dt,i＝1，2，…，n
将每一个能量值归一化，得到归一化的能量分布向量
pi＝Ei/E,i＝1,2,…,n
其中，
Figure GDA0002235838320000073
第一个分量p1表示最高频段的能量，代表了信号在最高频段范围内能量分布的比例，最后一个分量pn表示信号在最低频段范围内能量分布的比例；根据归一化的能量分布向量绘制归一化能量分布图，其中横坐标表示分量层次，纵坐标表示归一化的能量分布向量值，曲线表示平均值，误差棒表示标准差；
(6)、从physionet的meditation数据库获取给定时间和给定采样频率下的太极训练受试者的冥想训练中的ECG信号，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练中的RRI信号，返回步骤(2)至步骤(5)得到冥想训练中的ECG信号的归一化的能量分布向量；其中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移；
(7)、将从physionet的meditation数据库获取给定时间和给定采样频率下的太极训练受试者的冥想训练前RRI信号数据随机打乱形成训练前随机信号，返回步骤(2)至步骤(5)得到冥想训练前的ECG随机信号的归一化的能量分布向量；
(8)、将从physionet的meditation数据库获取给定时间和给定采样频率下的太极训练受试者的冥想训练中RRI信号数据随机打乱形成训练中随机信号，返回步骤(2)至步骤(5)得到冥想训练中的ECG随机信号的归一化的能量分布向量；
本发明将太极的HRV信号以及其随机信号分解得的极值模态函数分量Ci(t)，计算获得归一化能量分布向量，并作EED曲线图，如图8所示。图8为太极训练RR间期信号的EED分析，数据长度为3000，曲线上方的*表示训练前和训练中能量值T检验p<0.01。
如图8所示，太极训练前的EED曲线比较平滑，随着分解层次增加呈递增趋势；而在训练中表现出由高到低的变化趋势；太极训练中的能量值在低层次(层次2、3)高于训练前，在较高层次(层次4、5、6)低于训练前，并且在层次2、3、5、6具有显著性差异，另外可以看出冥想训练中的EED曲线与替代数据EED曲线分布非常相似。
为了评估冥想训练前与冥想训练中的EED曲线在低层次和高层次分解上的能量分布差异，将分量低层次能量和分量高层次能量相减，得到参数能量差值(EnergyDifferntial Value，EDV)EDV＝(p2+p3+p4)-(p5+p6+p7)，高的EDV值表示RRI信号更高的分量低层次能量分布和更低的分量高层次能量分布，太极训练的EDV表如表1所示。
表1太极训练EDV值
Figure GDA0002235838320000081
*表示太极训练前和训练中T检验p<0.0001。
**表示替代数据与其原始数据T检验p<0.0001。
从表1可得到，太极训练中的EDV值高于训练前的七倍多，差异显著；太极训练前EDV值小于其替代数据九倍多，可见明显小于；太极训练中的EDV值与其替代数据并没有表现出显著性差异；即太极冥想训练打破了原有的能量分布，显示出与随机序列一样的特征，说明冥想训练导致HRV信号的长时相关性消失。
瑜伽训练受试者的一种基于极值能量分解法的冥想训练的HRV信号分析方法，包括如下步骤：
(1)、如图7所示，从physionet的meditation数据库获取给定时间和给定采样频率下的太极训练受试者冥想训练前的ECG信号，该数据库包含4个瑜伽训练受试者(年龄20～52,平均33)，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练前RRI信号x(t)；其中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移；
(2)、将去噪后的RRI信号x(t)作为原始信号，其中原始信号x(t)所需最少数据量N＝2n+1，其中n为分解出的极值模态函数分量的数量；求出原始信号的所有局部极值点，然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin，得到上、下包络线的包络均值信号m(t)＝(emax+emin)/2；
(3)、将原始信号x(t)减去包络均值信号m(t)，得到h(t)＝x(t)-m(t)；然后判断h(t)是否满足极值模态函数的判定条件，如果不满足，将h(t)作为原始信号返回至步骤(2)，直到hk(t)满足极值模态函数的判定条件，则记c1(t)＝hk(t)，作为第一个极值模态函数分量；其中极值模态函数的判定条件为：(a)、在整个数据序列中，极值点的数量与过零点的数量相等或者相差一个；(b)、在任意时刻，上下包络线对于时间轴对称；
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t)，得到余量r1(t)＝x(t)-c1(t)，然后判断hk(t)是否满足停止准则，如果不满足，将r1(t)作为新的原始序列x(t)，返回至步骤(2)和(3)；如果hk(t)满足停止准则但n＜8时，返回步骤(1)重新获取原始信号；如果hk(t)满足停止准则且n≥8时，得到第2、3、…、n个极值模态函数分量及余量rn(t)，于是将原始信号x(t)分解为n个极值模态函数分量和一个余量，即
Figure GDA0002235838320000091
其中hk(t)满足停止准则公式为：
Figure GDA0002235838320000092
ε表示筛选门限，取0.2～0.3之间；
(5)、将原始信号x(t)分解得的n个极值模态函数分量，代表了原始信号不同频段的分量，然后计算其各个分量的能量
Ei＝∫|ci(t)|2dt,i＝1，2，…，n
将每一个能量值归一化，得到归一化的能量分布向量
pi＝Ei/E,i＝1,2,…,n
其中，
Figure GDA0002235838320000093
第一个分量p1表示最高频段的能量，代表了信号在最高频段范围内能量分布的比例，最后一个分量pn表示信号在最低频段范围内能量分布的比例；根据归一化的能量分布向量绘制归一化能量分布图，其中横坐标表示分量层次，纵坐标表示归一化的能量分布向量值，曲线表示平均值，误差棒表示标准差；
(6)、从physionet的meditation数据库获取给定时间和给定采样频率下的瑜伽训练受试者的冥想训练中的ECG信号，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练中的RRI信号，返回步骤(2)至步骤(6)得到冥想训练中的ECG信号的归一化的能量分布向量；其中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移；
(7)、将从physionet的meditation数据库获取给定时间和给定采样频率下的瑜伽训练受试者的冥想训练前RRI信号数据随机打乱形成训练前随机信号，返回步骤(2)至步骤(5)得到冥想训练前的ECG随机信号的归一化的能量分布向量；
(8)、将从physionet的meditation数据库获取给定时间和给定采样频率下的瑜伽训练受试者的冥想训练中RRI信号数据随机打乱形成训练中随机信号，返回步骤(2)至步骤(5)得到冥想训练中的ECG随机信号的归一化的能量分布向量；
本发明将瑜伽的HRV信号以及其随机信号分解得的极值模态函数分量Ci(t)，计算获得归一化能量分布向量，并作EED曲线图，如图9所示，图9为瑜伽训练RR间期信号的EED分析，数据长度为437到1126(平均805)，曲线上方的*表示训练前和训练中能量值T检验p<0.01，△表示p<0.05；由于瑜伽数据长度较短，所有数据长度介于437到1126之间(平均805)，为了提高计算的准确性，将所有数据用于EED分析。
如图9所示，瑜伽训练前EED曲线比较平滑，随着分解层次增加呈递增趋势；而在训练中表现出由高到低的变化趋势。太极训练中的能量值在低层次(层次2、3)高于训练前，在较高层次(层次4、5、6)低于训练前，并且在层次2、5、6具有显著性差异；另外，可以看出冥想训练中的EED曲线与替代数据EED曲线分布非常相似。
为了评估冥想训练前与冥想训练中的EED曲线在低层次和高层次分解上的能量分布差异，将分量低层次能量和分量高层次能量相减，得到参数能量差值(EnergyDifferntial Value，EDV)EDV＝(p2+p3+p4)-(p5+p6+p7)，高的EDV值表示RRI信号更高的分量低层次能量分布和更低的分量高层次能量分布，瑜伽训练的EDV表如表2所示。
表2瑜伽训练EDV值
Figure GDA0002235838320000101
△表示瑜伽训练前和训练中T检验p<0.0005。
△△表示替代数据与其原始数据T检验p<0.0005。
从表2可知，瑜伽冥想训练中的EDV值高于训练前近20倍，差异显著；冥想训练前的EDV值小于其替代数据近20倍，明显小于；训练中EDV值与其替代数据并没有表现出显著性差异，即表明瑜伽冥想训练打破了原有的能量分布，显示出与随机序列一样的特征，说明冥想训练导致HRV信号的长时相关性消失。

Claims (7)

1.一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于，包括如下步骤：
(1)、获取给定时间和给定采样频率下的冥想训练中的ECG信号，然后对ECG信号进行去噪预处理，从中提取RRI信号，得到冥想训练中RRI信号x(t)；
(2)、将RRI信号x(t)作为原始信号，求出原始信号的所有局部极值点，然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin，得到上、下包络线的包络均值信号m(t)＝(emax+emin)/2；
(3)、将原始信号x(t)减去包络均值信号m(t)，得到h(t)＝x(t)-m(t)；然后判断h(t)是否满足极值模态函数的判定条件，如果不满足，将h(t)作为原始信号返回至步骤(2)，直到hk(t)满足极值模态函数的判定条件，则记c1(t)＝hk(t)，作为第一个极值模态函数分量；
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t)，得到余量r1(t)＝x(t)-c1(t)，然后判断hk(t)是否满足停止准则，如果不满足，将r1(t)作为新的原始序列x(t)，返回至步骤(2)和(3)；如果hk(t)满足停止准则但n＜8时，返回步骤(1)重新获取原始信号；如果hk(t)满足停止准则且n≥8时，得到第2、3、…、n个极值模态函数分量及余量rn(t)，于是将原始信号x(t)分解为n个极值模态函数分量和一个余量，即
Figure FDA0002193891170000011
(5)、将原始信号x(t)分解得的n个极值模态函数分量，代表了原始信号不同频段的分量，然后计算其各个分量的能量
Ei＝∫|ci(t)|2dt,i＝1，2，…，n
将每一个能量值归一化，得到归一化的能量分布向量
pi＝Ei/E,i＝1,2,…,n
其中，
Figure FDA0002193891170000012
第一个分量p1表示最高频段的能量，代表了信号在最高频段范围内能量分布的比例，最后一个分量pn表示信号在最低频段范围内能量分布的比例；根据归一化的能量分布向量绘制归一化能量分布图，其中横坐标表示分量层次，纵坐标表示归一化的能量分布向量值，曲线表示平均值，误差棒表示标准差；
(6)、分别选取第二个分量p2至第七个分量p7，计算冥想训练中RRI信号的能量差异值EDV，EDV＝(p2+p3+p4)-(p5+p6+p7)，当冥想训练中RRI信号的EDV大于EDV标准值，则判定该冥想训练中的RRI信号为已进入冥想状态的RRI信号。
2.根据权利要求1的一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于：所述原始信号x(t)所需最少数据量N＝2n+1，其中n为分解出的极值模态函数分量的数量。
3.根据权利要求1的一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于：所述步骤(1)中去噪预处理的具体方法为：将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声，然后经过中值滤波器去除基线漂移。
4.根据权利要求1的一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于：所述步骤(3)中极值模态函数的判定条件为：(a)、在整个数据序列中，极值点的数量与过零点的数量相等或者相差一个；(b)、在任意时刻，上下包络线对于时间轴对称。
5.根据权利要求1的一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于：所述步骤(4)中hk(t)满足停止准则公式为：
Figure FDA0002193891170000021
ε表示筛选门限，取0.2～0.3之间。
6.根据权利要求1的一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于：所述冥想训练为太极训练，其步骤(6)中的标准值为0.26。
7.根据权利要求1的一种基于极值能量分解法的冥想训练的HRV信号分析方法，其特征在于：所述冥想训练为瑜伽训练，其步骤(6)中的标准值为0.68。
CN201910841609.4A 2019-09-06 2019-09-06 一种基于极值能量分解法的冥想训练的hrv信号分析方法 Active CN110558959B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201910841609.4A CN110558959B (zh) 2019-09-06 2019-09-06 一种基于极值能量分解法的冥想训练的hrv信号分析方法
PCT/CN2019/121417 WO2021042592A1 (zh) 2019-09-06 2019-11-28 一种基于极值能量分解法的冥想训练的hrv信号分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910841609.4A CN110558959B (zh) 2019-09-06 2019-09-06 一种基于极值能量分解法的冥想训练的hrv信号分析方法

Publications (2)

Publication Number Publication Date
CN110558959A CN110558959A (zh) 2019-12-13
CN110558959B true CN110558959B (zh) 2020-09-01

Family

ID=68778220

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910841609.4A Active CN110558959B (zh) 2019-09-06 2019-09-06 一种基于极值能量分解法的冥想训练的hrv信号分析方法

Country Status (2)

Country Link
CN (1) CN110558959B (zh)
WO (1) WO2021042592A1 (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114159077B (zh) * 2022-02-09 2022-05-31 浙江强脑科技有限公司 基于脑电信号的冥想评分方法、装置、终端及存储介质
CN114652330B (zh) * 2022-02-11 2023-03-24 北京赋思强脑科技有限公司 一种基于历史脑电信号评估冥想训练的方法、装置和设备
CN116548928B (zh) * 2023-07-11 2023-09-08 西安浩阳志德医疗科技有限公司 一种基于互联网的护理服务系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101496716A (zh) * 2009-02-26 2009-08-05 周洪建 利用ecg信号检测睡眠呼吸暂停的测量方法
CN101642369A (zh) * 2008-08-04 2010-02-10 南京大学 自主神经功能生物反馈方法和系统
CN107085651A (zh) * 2016-02-12 2017-08-22 飞比特公司 可穿戴装置及操作所述可穿戴装置的方法
CN108814579A (zh) * 2018-04-16 2018-11-16 西安交通大学 一种基于emd分解的心电、呼吸联合计算心率变异性的方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080269628A1 (en) * 2007-04-25 2008-10-30 Siemens Medical Solutions Usa, Inc. Denoising and Artifact Rejection for Cardiac Signal in a Sensis System
US9451898B2 (en) * 2014-07-02 2016-09-27 National Central University Method and system for extracting ventricular fibrillation signals in electrocardiogram using spline interpolation with uniform phase ensembles
CN109567743A (zh) * 2018-10-24 2019-04-05 加康康健有限公司 一种基于emd的信号重构方法、装置、终端设备及存储介质
CN109998527A (zh) * 2019-04-09 2019-07-12 湖北工业大学 一种基于多尺度熵的心脏疾病检测方法
CN110141206A (zh) * 2019-05-31 2019-08-20 四川长虹电器股份有限公司 一种基于集合经验模态分解的生理电信号分析方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101642369A (zh) * 2008-08-04 2010-02-10 南京大学 自主神经功能生物反馈方法和系统
CN101496716A (zh) * 2009-02-26 2009-08-05 周洪建 利用ecg信号检测睡眠呼吸暂停的测量方法
CN107085651A (zh) * 2016-02-12 2017-08-22 飞比特公司 可穿戴装置及操作所述可穿戴装置的方法
CN108814579A (zh) * 2018-04-16 2018-11-16 西安交通大学 一种基于emd分解的心电、呼吸联合计算心率变异性的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
一种基于改进经验模态分解的癫痫脑电识别新方法;庞春颖等;《中国生物医学工程学报》;20131231;第32卷(第6期);第663-669页 *
基于Hilbert谱的心率变异信号时频分析方法;董红生等;《仪器仪表学报》;20110228;第32卷(第2期);第271-278页 *

Also Published As

Publication number Publication date
CN110558959A (zh) 2019-12-13
WO2021042592A1 (zh) 2021-03-11

Similar Documents

Publication Publication Date Title
CN110558973B (zh) 一种执行基于极值能量分解法的心电图信号量化分析方法的计算机设备
CA2979135C (en) Systems, apparatus and methods for sensing fetal activity
CN110558959B (zh) 一种基于极值能量分解法的冥想训练的hrv信号分析方法
CN110464337B (zh) 一种基于极值能量分解法的心率变异性信号分析方法
US7809433B2 (en) Method and system for limiting interference in electroencephalographic signals
Castiglioni et al. Local scale exponents of blood pressure and heart rate variability by detrended fluctuation analysis: effects of posture, exercise, and aging
Costin et al. Atrial fibrillation onset prediction using variability of ECG signals
CN106539580B (zh) 一种自主神经系统动态变化的连续监测方法
CN109620206B (zh) 包含异位心跳判断的房颤人工智能识别方法及装置
CN112057087B (zh) 精神分裂症高风险人群自主神经功能数据处理方法及装置
CN108814579A (zh) 一种基于emd分解的心电、呼吸联合计算心率变异性的方法
CN110558974B (zh) 一种基于极值能量分解法的心电图信号分析方法
Hugeng et al. Development of the ‘Healthcor’system as a cardiac disorders symptoms detector using an expert system based on arduino uno
Gajare et al. MATLAB-based ECG R-peak Detection and Signal Classification using Deep Learning Approach
Khavas et al. Robust heartbeat detection using multimodal recordings and ECG quality assessment with signal amplitudes dispersion
Ghosh et al. Pre-ictal epileptic seizure prediction based on ECG signal analysis
Mashin et al. Analysis of the heart rate variability in negative functional states in the course of psychological relaxation sessions
Patil et al. Discrimination between atrial fibrillation (AF) & normal sinus rhythm (NSR) using linear parameters
Georgieva-Tsaneva An Information System for Cardiological Data Studying
Thalange et al. HRV analysis of arrhythmias using linear-nonlinear parameters
Tobón et al. Improved heart rate variability measurement based on modulation spectral processing of noisy electrocardiogram signals
Bouziane et al. The ANS sympathovagal balance using a hybrid method based on the wavelet packet and the KS-segmentation algorithm
Jabloun et al. Phase-rectified signal averaging method applied to heart rate variability signals for assessment of the changes in sympathovagal balance during rest and tilt
Zhuravlev et al. The Modern Information System of Sleep Structure Detection in Polysomnographic Data Systems Approaches
Zhang et al. Study of pulse rate variability signals using bispectrum analysis

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