CN109635786A - 一种基于心音自相关函数的心率计算算法 - Google Patents
一种基于心音自相关函数的心率计算算法 Download PDFInfo
- Publication number
- CN109635786A CN109635786A CN201910057743.5A CN201910057743A CN109635786A CN 109635786 A CN109635786 A CN 109635786A CN 201910057743 A CN201910057743 A CN 201910057743A CN 109635786 A CN109635786 A CN 109635786A
- Authority
- CN
- China
- Prior art keywords
- autocorrelation function
- sequence
- heart sound
- heart
- envelope
- 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.)
- Granted
Links
- 238000005311 autocorrelation function Methods 0.000 title claims abstract description 79
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 34
- 238000001914 filtration Methods 0.000 claims abstract description 22
- 238000012545 processing Methods 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 34
- 238000005070 sampling Methods 0.000 claims description 14
- 230000009191 jumping Effects 0.000 claims description 11
- 238000000034 method Methods 0.000 claims description 10
- 238000010606 normalization Methods 0.000 claims description 3
- 238000002555 auscultation Methods 0.000 abstract description 8
- 238000003745 diagnosis Methods 0.000 abstract description 2
- 230000005236 sound signal Effects 0.000 description 20
- 208000037656 Respiratory Sounds Diseases 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 4
- 230000029058 respiratory gaseous exchange Effects 0.000 description 4
- 230000036387 respiratory rate Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000002349 favourable effect Effects 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 208000008203 tachypnea Diseases 0.000 description 2
- 206010043089 tachypnoea Diseases 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 208000019622 heart disease Diseases 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
- G16H50/20—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for computer-aided diagnosis, e.g. based on medical expert systems
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H80/00—ICT specially adapted for facilitating communication between medical practitioners or patients, e.g. for collaborative diagnosis, therapy or health monitoring
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Public Health (AREA)
- Biomedical Technology (AREA)
- Epidemiology (AREA)
- Pathology (AREA)
- Theoretical Computer Science (AREA)
- Primary Health Care (AREA)
- General Health & Medical Sciences (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Signal Processing (AREA)
- General Physics & Mathematics (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- Data Mining & Analysis (AREA)
- Databases & Information Systems (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Abstract
本发明提供了一种基于心音自相关函数的心率计算算法，其特征在于：获取心音序列，对心音序列作滤波处理和求包络；之后在搜索范围内计算自相关函数，搜索自相关函数的峰值点进而根据峰值点序号计算心率值。本发明计算算法可克服呼吸音、摩擦音以及第一、第二心音波形相似等因素的干扰，利用心音序列即可准确计算出心率值，便于在小型可穿戴设备和小型电子听诊器中实现，也适用于在医院和家庭用电子听诊辅助诊疗系统。
Description
技术领域
本发明涉及医疗器械与医学信号处理领域，更具体地说，涉及一种基于心音自相关函数的心率计算算法。
背景技术
心率是反映身体健康状况的重要指标之一，心率监测对于部分心脏病患者用药、家庭护理等均具有重要意义。
随着远程医疗的普及，电子听诊器将进入更多社区医院和家庭，心音的远程听诊也将日益成为包括心脏病患者家庭护理在内的众多远程医疗项目的重要手段之一。相对于基于脉搏测量心率等常用方法，利用心音信号计算心率准确度更高，使用方便且成本低廉。因此，基于心音计算心率势必成为今后心率监测的主要方式之一。
基于心音信号计算心率通常可通过确定心音包络峰值点间距离或心音自相关函数峰值点间距离来实现。但这些方法在实际应用中存在以下不利因素：
1)包络和自相关函数均会因呼吸音、摩擦音等因素的干扰而出现异常的波峰，导致心率计算错误；例如，即便是很浅的呼吸都可能因胸脯的轻微起伏使听诊器移位而出现摩擦音；而呼吸较急促时，这种干扰将更显著；
2)有的心音同时存在第一和第二心音，由于第一和第二心音波形相似，若直接应用传统的相关分析方法计算心率，容易得到虚假的周期值；虽然通过识别并分割第一和第二心音可克服这种因素的影响，但是相当部分心音第一和第二心音的分界并不明显，因此，不能确保通过去掉第二心音所致包络峰值或自相关函数峰值的形式来确定心率；
3)采样时间不可能过长；一般来说，听诊持续时间通常不超过半分钟，有时甚至更短，所以基于心音信号监测心率仅能根据数个心跳周期的心音信号来进行心率估算，此时现有的自相关函数的计算量极大，当采样频率较高时，为在短时间内完成自相关函数计算，往往需配备计算能力较强的硬件设备，导致硬件成本偏高。
故此，采用现有的方案和手段来实现基于心音估算心率，还存在很多困难。
发明内容
本发明的目的在于克服现有技术中的缺点与不足，提供一种基于心音自相关函数计算、可克服呼吸音等因素的干扰、只利用心音序列即可准确计算出心率值的心率计算算法。
为了达到上述目的，本发明通过下述技术方案予以实现：一种基于心音自相关函数的心率计算算法，其特征在于：包括如下步骤：
第一步，以采样频率fs获取持续时间为t、长度为N的心音序列x(k)，其中， k＝1，...，N，N为正整数，N＝tfs；
第二步，对心音序列x(k)作滤波处理得到序列
第三步，以频率对序列欠采样得到序列s(i)，其中，round()为四舍五入取整函数；
以LT为半窗长求序列s(i)的包络En(m)；其中，
第四步，设定移位序号为n；在移位序号的搜索范围内将包络直接进行自相关函数R(n)计算；或者是将包络进行处理后，再在移位序号的搜索范围内将处理后的包络进行自相关函数R(n)计算；其中，c1为上限系数，c2为下限系数；
第五步，初始化移位序号n、峰值点序号nP和峰值RP：令移位序号峰值点序号峰值RP＝0；
第六步，将移位序号n自加1；判断当前移位序号n和对应的自相关函数R(n) 数值大小：
若R(n)＞R(n-1)，且R(n)＞R(n+1)，且R(n)＞RP，且|n-ainP|＞nc，ai＝2，3，..，则将当前移位序号n作为峰值点序号nP，并将对应的自相关函数R(n)数值作为峰值RP；其中，nc为差值下限；之后跳至第七步；
否则直接跳至第七步；
第七步，判断当前移位序号n的大小：若则返回至第六步；否则跳转至第八步；
第八步，计算心率值：单位为次/分钟。
本发明计算算法的原理是：根据心音信号的特点，通过限定自相关函数峰值搜索范围、增加滤波等方式抑制呼吸音、摩擦音等因素引起的干扰。在搜索范围内计算自相关函数，通过搜索自相关函数的最大峰值点，并采用出现最大峰值点时的移位距离来确定心率周期。该方式的好处是：第一心音和第二心音产生机制相近，两者波形亦相似，导致心音自相关函数在单个心音周期内会出现因第一、第二心音之间的互相关而形成的多个峰值。但由于第一心音和第二心音波形存在一定差别，其因互相关所形成的峰值幅值较低。故此，本发明第六步中寻找R(n)的最大峰值点，并在第八步通过该峰值点的序号计算心率，而不是简单地通过自相关函数相邻峰值点间距离计算心率，以避免因第一、第二心音波形相似造成的心率计算错误。由于周期信号的自相关函数在基频的倍频处会出现峰值；若心音信号质量不佳，极可能出现自相关函数在倍频处的峰值高于基频处的峰值，而导致心率计算出现错误。故此，本发明第六步中采用 |n-ainP|＞nc，ai＝2，3，..作为峰值点选择的必要条件之一，可避免该错误的发生。此外本发明计算算法限定自相关函数峰值搜索范围，既可准确计算心率值，又可降低计算算法计算量，快速得出结果，有利于降低硬件实现要求和生产成本。
当采样频率较高时，直接计算实测听诊信号自相关函数的计算量极大；本发明在第三步通过欠采样和包络计算来降低所需计算自相关函数的序列长度。以采样频率为22050Hz并取且LT＝10为例，采用欠采样和包络计算，则处理前后所需计算自相关函数的序列的长度比高达220∶1；因此本发明计算算法计算量大大降低，从而降低对硬件实现的要求。
该计算算法既可以在小型可穿戴设备和小型电子听诊器中实现，又可以在大型电子听诊设备中实现。
优选地，所述第一步中，持续时间为t的取值范围为：2.5s＜t＜4s。
正常心率值范围为45＜HR＜180次/分钟，而呼吸频率一般不会超过40次/ 分钟，呼吸频率超过24次/分钟已属呼吸急促。故此，本发明第一步中选取持续时间落于2.5秒到4秒之间的心音序列来计算心率，可保证所选取的信号段至少包含两次以上心跳，满足自相关函数计算心率所需的心跳数；同时，也保证了在所选信号段内不会出现两次或两次以上呼吸，有效避免呼吸音对于心音周期计算的影响。
优选地，所述第二步中，对心音序列x(k)作滤波处理得到序列是指：对心音序列x(k)作截止频率为fc的低通滤波得到序列以抑制呼吸音和摩擦音对心音信号的干扰。
优选地，所述截止频率fc的取值范围为：fc∈[150，250]Hz；低通滤波是指二阶以上的低通滤波。相对于心音信号，呼吸音信号和摩擦音信号在高于250Hz 的频段中成份占比较高。而决定心音包络的成份则主要集中于心音信号的低频段[0，250]Hz，尤其是[0，150]Hz频段内。即仅根据心音信号在[0，150]Hz频段内的成份即可有效确定心率。故此，本发明第二步中对心音序列作截止频率为 fc∈[150，250]Hz的低通滤波，有助于降低呼吸音和摩擦音的幅值，抑制其对于心音包络计算以及对于心率计算的干扰。
优选地，所述第三步中，频率的取值范围为：
半窗长LT的取值范围为：
包络En(m)的计算方法为：
或者是
优选地，所述第四步中，上限系数c1和下限系数c2的取值范围分别为： c1∈[0.25，0.27]，c2∈[1.2，1.5]。
考虑到心率值的范围为45＜HR＜180次/分钟，本发明第四步仅计算搜索范围内的自相关函数R(n)，将和分别代入到心率值中，分别对应于心率下限和心率上限：取c2∈[1.2，1.5]相当于默认心率下限为40～50次/ 分钟，取c1∈[0.25，0.27]相当于默认心率上限为222～240次/分钟。该上下限的选择有利于排除呼吸音的影响，避免将呼吸频率误作为心率。相应地，第六步中，判断条件|n-ainP|＞nc的ai优选取值为ai＝2，3；即判断|n-2nP|＞nc，且 |n-3nP|＞nc。
其中一种方案是：所述第四步中，在移位序号的搜索范围内将包络直接进行自相关函数R(n)计算是指：计算包络En(m)的自相关函数R(n)：
其中，
另一种方案是：所述第四步中，将包络进行处理后，再在移位序号的搜索范围内将处理后的包络进行自相关函数R(n)计算，是指：
对包络En(m)作归一化处理：其中，max(En)为En(m) 的最大值；之后在移位序号的搜索范围内计算包络的自相关函数R(n)。
优选地，所述在移位序号的搜索范围内计算包络的自相关函数R(n)，是指：
其中，
与现有技术相比，本发明具有如下优点与有益效果：
1、本发明计算算法可克服呼吸音、摩擦音以及第一、第二心音波形相似等因素的干扰，只利用心音信号即可实现心率计算，准确程度高；
2、本发明计算算法计算量小，对硬件实现的要求低，非常便于在小型可穿戴设备和小型电子听诊器中实现，同时本发明计算算法也特别适用于在医院和家庭用电子听诊辅助诊疗系统中实现心率监测。
附图说明
图1是实施例一基于心音自相关函数的心率计算算法的流程图；
图2(a)～图2(e)是本发明实施例中计算算法各步骤所得结果；
图3(a)～图3(h)是本发明实施例中滤波前后呼吸音、心音信号曲线及相应频谱；
图4是实施例二基于心音自相关函数的心率计算算法的流程图。
具体实施方式
下面结合附图与具体实施方式对本发明作进一步详细的描述。
实施例一
本实施例基于心音自相关函数的心率计算算法，其工作流程如图1所示，包括如下步骤：
第一步，以采样频率fs获取持续时间为t、长度为N的心音序列x(k)，其中 k＝1，...，N，N为正整数，N＝tfs；持续时间为t的取值范围优选为： 2.5s＜t＜4s；
第二步，对心音序列x(k)作滤波处理得到序列优选是对心音序列x(k) 作截止频率为fc的低通滤波得到序列截止频率fc的取值范围优选为：fc∈[150，250]Hz；低通滤波可采用现有低通滤波方法，优选采用二阶以上低通滤波；采用其它可抑制呼吸音和摩擦音干扰的滤波处理也可；
第三步，以频率对序列欠采样得到序列s(i)，其中，以LT为半窗长求序列s(i)的包络 En(m)；其中，round()为四舍五入取整函数：
频率的取值范围优选为：
半窗长LT的取值范围优选为：
包络En(m)的计算方法优选为：
或者是
第四步，设定移位序号为n；将包络进行处理后，再在移位序号的搜索范围内将处理后的包络进行自相关函数R(n)计算；其中，c1为上限系数，c2为下限系数；上限系数c1和下限系数c2的取值范围分别优选为：c1∈[0.25，0.27]，c2∈[1.2，1.5]；
将包络进行处理是指：对包络En(m)作归一化处理：其中，max(En)为En(m)的最大值；之后在移位序号的搜索范围内计算包络的自相关函数R(n)：
实际应用中自相关函数也可以采用其它计算公式。
第五步，初始化移位序号n、峰值点序号nP和峰值RP：令移位序号峰值点序号峰值RP＝0；
第六步，将移位序号n自加1；判断当前移位序号n和对应的自相关函数R(n) 数值大小：
若R(n)＞R(n-1)，且R(n)＞R(n+1)，且R(n)＞RP，且|n-ainP|＞nc，ai＝2，3，..，则将当前移位序号n作为峰值点序号nP，并将对应的自相关函数R(n)数值作为峰值RP；其中，nc为差值下限；差值下限nc的取值范围优选为：nc∈[5，10]；之后跳至第七步；
否则直接跳至第七步；
第七步，判断移位序号n的大小：若则返回至第六步，否则跳转至第八步；
第八步，计算心率值：单位为次/分钟。
图2(a)是实测听诊信号，图2(b)是低通滤波后信号，图2(c)是欠采样所得信号，图2(d)是欠采样信号的包络，图2(e)包络的自相关函数。
图3(a)是兼有心音和呼吸音的信号，图3(b)是兼有心音和呼吸音的信号频谱，图3(c)是兼有心音和呼吸音的信号低通滤波后曲线，图3(d)是兼有心音和呼吸音的信号低通滤波后频谱，图3(e)是滤波前心音信号，图3(f)是滤波前心音频谱，图3(g)是低通滤波后心音信号，图3(h)是低通滤波后心音频谱。
本发明计算算法的原理是：
如图3(a)至图3(d)所示，相对于心音信号，呼吸音信号和摩擦音信号在高于250Hz的频段中成份占比较高。而决定心音包络的成份则主要集中于心音信号的低频段[0，250]Hz，尤其是[0，150]Hz频段内，如图3(e)至图3(h)所示。即仅根据心音信号在[0，150]Hz频段内的成份即可有效确定心率。故此，本发明第二步中对心音序列作截止频率为fc∈[150，250]Hz的低通滤波，有助于降低呼吸音和摩擦音的幅值(可对比图2(a)和图2(b)含呼吸音和摩擦音的片段在低通滤波前后曲线的幅值)，抑制其对于心音包络计算以及对于心率计算的干扰。
正常心率值范围为45＜HR＜180次/分钟，而呼吸频率一般不会超过40次/ 分钟，呼吸频率超过24次/分钟已属呼吸急促。故此，本发明第一步中选取持续时间落于2.5秒到4秒之间的心音序列来计算心率，可保证所选取的信号段至少包含两次以上心跳，满足自相关函数计算心率所需的心跳数；同时，也保证了在所选信号段内不会出现两次或两次以上呼吸，有效避免呼吸音对于心音周期计算的影响。
考虑到心率值的范围为45＜HR＜180次/分钟，本发明第四步仅计算搜索范围c1∈[0.25，0.27]，c2∈[1.2，1.5]，内的自相关函数R(n)，将和分别代入到心率值中，分别对应于心率下限和心率上限。由于故取c2∈[1.2，1.5]相当于默认心率下限为40～50次/分钟，取c1∈[0.25，0.27]相当于默认心率上限为222～240次/分钟(见图2(e)中所示峰值点搜索范围)。该上下限的选择有利于排除呼吸音的影响，避免将呼吸频率误作为心率。
周期信号的自相关函数在基频的倍频处会出现峰值。若心音信号质量不佳，极可能出现自相关函数在倍频处的峰值高于基频处的峰值(见图2(e)中所示心音包络自相关函数倍频处波峰)，而导致心率计算出现错误。故此，本发明第六步中采用|n-ainP|＞nc，ai＝2，3，..作为峰值点选择的必要条件之一，可避免该错误的发生。
第一心音和第二心音产生机制相近，两者波形亦相似，导致心音自相关函数在单个心音周期内会出现因第一、第二心音之间的互相关而形成的多个峰值(见图2(e)中所示因第一、第二心音波形相似所致心音包络自相关函数峰值点)。但由于第一心音和第二心音波形存在一定差别，其因互相关所形成的峰值幅值较低。故此，本发明第六步中寻找R(n)的最大峰值点，并在第八步通过该峰值点的序号计算心率，而不是简单地通过R(n)相邻峰值点间距离计算心率，以避免因第一、第二心音波形相似造成的心率计算错误。
当采样频率较高时，直接计算实测听诊信号自相关函数的计算量极大；本发明第三步通过包络计算或者是欠采样和包络计算来降低所需计算自相关函数的序列长度；其中，通过欠采样和包络计算可进一步降低所需计算自相关函数的序列长度。以采样频率为22050Hz并取且LT＝10为例，采用S1 步到S3步中的欠采样和包络计算，则处理前后所需计算自相关函数的序列的长度比高达220：1；因此本发明计算算法计算量大大降低，从而降低对硬件实现的要求。
下面以具体示例对本发明计算算法进行说明：
一种基于心音自相关函数的心率计算算法，包括如下步骤：
第一步：取长度为N＝77175点、采样频率为fs＝22050Hz的心音序列x(k)，其中k＝1，...，N，持续时间t＝N/fs，显然所取持续时间t满足条件2.5s＜t＜4s，即所分析的心音序列持续时间落于2.5秒到4秒之间。
应说明的是，心率计算对心音的保真程度要求不高，采样频率可行的取值范围可宽至fs∈[2000，48000]Hz。
第二步：对x(k)作截止频率为fc＝200Hz、二阶低通滤波得滤波后序列
其中，低通滤波参数采用巴特沃斯低通滤波器设计计算算法，由截止频率fc＝200Hz和采样频率fs＝22050Hz共同决定。
第三步，取频率则以对欠采样得其中，i＝1，...，3500。
取LT＝10，求s(i)的包络其中，
第四步，对En(m)作归一化处理：
分别取c1＝0.25，c2＝1.5，则 计算的自相关函数：
其中，
第五步，分别取峰值点序号峰值RP＝0。
第六步，n＝n+1；取nc＝10，若R(n)＞R(n-1)，且R(n)＞R(n+1)，且R(n)＞RP，且|n-ainP|＞nc，ai＝2，3，..，则令RP＝R(n)，且nP＝n；之后跳至第七步；否则直接跳至第七步。
第七步，若则返回到第六步，否则跳转至第八步。
第八步，计算心率得次/分钟。
实施例二
本实施例基于心音自相关函数的心率计算算法，其工作流程如图4所示，与实施例一的区别在于：本实施例中，第四步，在移位序号的搜索范围内将包络En(m)直接进行自相关函数R(n)计算。
具体地说，是指：计算包络En(m)的自相关函数R(n)：
其中，
本实施例的其余步骤与实施例一相同。
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
Claims (10)
1.一种基于心音自相关函数的心率计算算法，其特征在于：包括如下步骤：
第一步，以采样频率fs获取持续时间为t、长度为N的心音序列x(k)，其中，k＝1,…,N，N为正整数，N＝tfs；
第二步，对心音序列x(k)作滤波处理得到序列
第三步，以频率对序列欠采样得到序列s(i)，其中，round()为四舍五入取整函数；
以LT为半窗长求序列s(i)的包络En(m)；其中，
第四步，设定移位序号为n；在移位序号的搜索范围内将包络直接进行自相关函数R(n)计算；或者是将包络进行处理后，再在移位序号的搜索范围内将处理后的包络进行自相关函数R(n)计算；其中，c1为上限系数，c2为下限系数；
第五步，初始化移位序号n、峰值点序号nP和峰值RP：令移位序号峰值点序号峰值RP＝0；
第六步，将移位序号n自加1；判断当前移位序号n和对应的自相关函数R(n)数值大小：
若R(n)>R(n-1)，且R(n)>R(n+1)，且R(n)>RP，且|n-ainP|>nc，ai＝2,3,...，则将当前移位序号n作为峰值点序号nP，并将对应的自相关函数R(n)数值作为峰值RP；其中，nc为差值下限；之后跳至第七步；
否则直接跳至第七步；
第七步，判断当前移位序号n的大小：若则返回至第六步；否则跳转至第八步；
第八步，计算心率值：单位为次/分钟。
2.根据权利要求1所述的基于心音自相关函数的心率计算算法，其特征在于：所述第一步中，持续时间为t的取值范围为：2.5s＜t＜4s。
3.根据权利要求1所述的基于心音自相关函数的心率计算算法，其特征在于：所述第二步中，对心音序列x(k)作滤波处理得到序列是指：对心音序列x(k)作截止频率为fc的低通滤波得到序列以抑制呼吸音和摩擦音对心音信号的干扰。
4.根据权利要求3所述的基于心音自相关函数的心率计算算法，其特征在于：所述截止频率fc的取值范围为：fc∈[150,250]Hz；低通滤波是指二阶以上的低通滤波。
5.根据权利要求1所述的基于心音自相关函数的心率计算算法，其特征在于：所述第三步中，频率的取值范围为：
半窗长LT的取值范围为：
包络En(m)的计算方法为：
或者是
6.根据权利要求1所述的基于心音自相关函数的心率计算算法，其特征在于：所述第四步中，上限系数c1和下限系数c2的取值范围分别为：c1∈[0.25,0.27]，c2∈[1.2,1.5]。
7.根据权利要求1所述的基于心音自相关函数的心率计算算法，其特征在于：所述第四步中，在移位序号的搜索范围内将包络直接进行自相关函数R(n)计算是指：计算包络En(m)的自相关函数R(n)：
其中，
8.根据权利要求1所述的基于心音自相关函数的心率计算算法，其特征在于：所述第四步中，将包络进行处理后，再在移位序号的搜索范围内将处理后的包络进行自相关函数R(n)计算，是指：
对包络En(m)作归一化处理：其中，max(En)为En(m)的最大值；之后在移位序号的搜索范围内计算包络的自相关函数R(n)。
9.根据权利要求8所述的基于心音自相关函数的心率计算算法，其特征在于：所述在移位序号的搜索范围内计算包络的自相关函数R(n)，是指：
其中，
10.根据权利要求1至9中任一项所述的基于心音自相关函数的心率计算算法，其特征在于：所述第六步中，差值下限nc的取值范围为：nc∈[5,10]。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910057743.5A CN109635786B (zh) | 2019-01-22 | 2019-01-22 | 一种基于心音自相关函数的心率计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910057743.5A CN109635786B (zh) | 2019-01-22 | 2019-01-22 | 一种基于心音自相关函数的心率计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109635786A true CN109635786A (zh) | 2019-04-16 |
| CN109635786B CN109635786B (zh) | 2022-05-17 |
Family
ID=66063018
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910057743.5A Active CN109635786B (zh) | 2019-01-22 | 2019-01-22 | 一种基于心音自相关函数的心率计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109635786B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113171073A (zh) * | 2021-05-18 | 2021-07-27 | 南京润楠医疗电子研究院有限公司 | 一种基于检波器的无感式心率检测方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080051669A1 (en) * | 2006-08-28 | 2008-02-28 | Wolfgang Meyer | Diagnosis of sleep apnea |
| CN103845079A (zh) * | 2014-02-14 | 2014-06-11 | 广东工业大学 | 一种基于盲分离的多普勒胎心音瞬时心率的检测方法 |
| CN106923812A (zh) * | 2017-03-30 | 2017-07-07 | 华南理工大学 | 一种基于心音信号自相关分析的心率计算方法 |
| CN108652611A (zh) * | 2018-05-17 | 2018-10-16 | 广东工业大学 | 一种心率的计算方法、装置、设备及计算机可读存储介质 |
-
2019
- 2019-01-22 CN CN201910057743.5A patent/CN109635786B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080051669A1 (en) * | 2006-08-28 | 2008-02-28 | Wolfgang Meyer | Diagnosis of sleep apnea |
| CN103845079A (zh) * | 2014-02-14 | 2014-06-11 | 广东工业大学 | 一种基于盲分离的多普勒胎心音瞬时心率的检测方法 |
| CN106923812A (zh) * | 2017-03-30 | 2017-07-07 | 华南理工大学 | 一种基于心音信号自相关分析的心率计算方法 |
| CN108652611A (zh) * | 2018-05-17 | 2018-10-16 | 广东工业大学 | 一种心率的计算方法、装置、设备及计算机可读存储介质 |
Non-Patent Citations (1)
| Title |
|---|
| 蒋飞朋 等: "一种用于电子听诊器的心跳闪烁指示电路", 《自动化与信息工程》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113171073A (zh) * | 2021-05-18 | 2021-07-27 | 南京润楠医疗电子研究院有限公司 | 一种基于检波器的无感式心率检测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109635786B (zh) | 2022-05-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Arzeno et al. | Analysis of first-derivative based QRS detection algorithms | |
| CN105852851B (zh) | 心房颤动检测系统 | |
| Zhu et al. | Real-time monitoring of respiration rhythm and pulse rate during sleep | |
| EP2939592B1 (en) | Detection of parameters in cardiac output related waveforms | |
| JP6595706B2 (ja) | 血圧の推移を判定する装置 | |
| EP2375982B1 (en) | Method and apparatus for the analysis of ballistocardiogram signals | |
| US11234658B2 (en) | Photoplethysmogram data analysis and presentation | |
| Castro et al. | Heart sound segmentation of pediatric auscultations using wavelet analysis | |
| EP3479758B1 (en) | System and method for breathing pattern extraction from ppg signals | |
| WO2019191487A1 (en) | Photoplethysmogram data analysis and presentation | |
| CN105705079A (zh) | 用于处理生理信号的设备、方法和系统 | |
| JP4644268B2 (ja) | シャント状態検知装置 | |
| CN111031902B (zh) | 多传感器心搏出量监测系统和分析法 | |
| CN110680307A (zh) | 一种运动环境下基于脉搏波传导时间的动态血压监测方法 | |
| Kong et al. | Heart rate tracking using a wearable photoplethysmographic sensor during treadmill exercise | |
| CN109635786B (zh) | 一种基于心音自相关函数的心率计算方法 | |
| CN114027804A (zh) | 一种脉象诊断的方法、装置及可读存储介质 | |
| Garg et al. | An effective method to identify various factors for denoising wrist pulse signal using wavelet denoising algorithm | |
| Santo et al. | Respiration rate extraction from ECG signal via discrete wavelet transform | |
| CN107280650B (zh) | 生命体特征参数获取方法及装置 | |
| Grochowina et al. | Comparison of SVM and k-NN classifiers in the estimation of the state of the arteriovenous fistula problem | |
| Sundaram et al. | Discriminating normal phonocardiogram from artifact using a multiscale entropy technique | |
| JP2015213714A (ja) | 心電波形検出装置、心電波形検出方法、心電波形検出プログラム、及び撮像装置 | |
| Tsimenidis et al. | Reliability of clinical alarm detection in intensive care units | |
| Kozumplik et al. | Respiratory rate estimation using the photoplethysmogram: Towards the implementation in wearables |
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 |