CN106923812A - 一种基于心音信号自相关分析的心率计算方法 - Google Patents

Abstract

本发明公开了一种基于心音信号自相关分析的心率计算方法，所述方法包括以下步骤：首先利用录音设备获取心音信号；然后对心音信号进行抽取及死区非线性处理，得到待分析处理的信号；之后对抽取得到的信号进行自相关运算获取自相关序列；接着对自相关序列再次进行死区非线性处理，从处理后得到的序列中提取出第二个周期的位置，根据第二个周期的峰值，计算出心音信号的周期，获得心率。本发明通过自相关处理以及死区非线性处理得到自相关序列，从自相关序列中计算出心率，只利用心音信号即可准确实现心率的计算，计算方法简单，结果精准，对硬件实现的要求低，非常便于在可穿戴设备中实现。

Description

一种基于心音信号自相关分析的心率计算方法

技术领域

本发明涉及医疗器械与医学信号处理领域，具体涉及一种基于心音信号自相关分析的心率计算方法。

背景技术

静息心率是反映身体健康状况的重要指标之一，心率监测对于心脏病患者的用药、家庭护理等均具有重要意义。

目前市面上常用的心率检测仪涉及到的检测方法有以下几种：

1.血氧定量法：血管中的氧含量随心脏搏动呈消耗——心脏泵血增加——再消耗的周期变化过程，且搏动周期与心率一致。血氧定量法使用波长660nm的红光和940nm的近红外光作为射入光源。血管中携氧的血红蛋白和不携氧的血红蛋白对这两种光的吸收率不同，故通过测定血管的光传导强度随时间变化的周期即可确定心率。这种测量方法多用于临床监测。

2.光电容积法：通过追踪可见光(绿光)在人体组织中的反射强度变化周期以确定心率。毛细血管和动、静脉的容积随脉搏而变化，该变化将周期性改变光的反射强度。测定反射光强度变化的频率可确定脉搏频率，而脉搏频率通常与心率一致。目前市面上运动腕表都采用这种测量方法。

3.心电信号法：心脏的搏动伴随心电的变化，通过测量体表心电可精确计算心率，医院用的心电仪、市面上最精确的可穿戴心率测量仪器、心率带均采用这种方法。

4.动脉血压法：通过压力传感器在手腕或者颈部两侧测量动脉压力有规律的涨落，可以确定脉搏频率，进而估计心率。其为最传统的心率测量方法，但由于需要长时间压迫血管造成使用不便，目前在商用中是最不成熟的，所以一般只在医院中对术中和术后的静息病人使用。

比较上述几种方法可知：血氧定量法、动脉血压法使用不便；对部分心脏病患者而言，脉搏频率不一定能准确反映心率值，此时动脉血压法、光电容积法测量结果存在误差。心电信号法测量的结果最准确，但是所需的设备价格也较高。

随着家用电子听诊器的普及，心音远程听诊将日益成为心脏病患者家庭护理的重要手段之一。而利用心音信号测量心率准确度高，使用方便且成本低廉。受成本、功耗和体积等限制，电子听诊器的数字处理器运算速度和存储空间有限；因此，通过心音信号计算心率所用的算法必须具有较低的计算复杂度，且对可靠性要求高，能适用于多种类型的心脏病患者。

相关分析是计算信号周期的常用方法，并已成功应用于语音信号基音周期的计算。但心音信号中第一和第二心音信号波形相似，若直接应用传统的相关分析方法计算心率，容易得到虚假的周期值。为在家用电子听诊器中实现心率计算功能，需要针对心音信号的特点对传统相关分析方法作适当的调整。

发明内容

本发明的目的是针对上述现有技术的不足，提供了一种基于心音信号自相关分析的心率计算方法，该方法根据心音信号的特点，引入特定的非线性变换改进传统的相关分析方法，进而基于电子听诊器所采集的心音信号计算心率，计算量低、可靠性高。

本发明的目的可以通过如下技术方案实现：

一种基于心音信号自相关分析的心率计算方法，所述方法包括以下步骤：

第一步、取长度为N点、采样频率为f_s的心音信号x(k)，其中k∈{0,…,N-1}，N为正整数，且N/f_s不小于3s，以保证所采集的心音信号x(k)至少包含三个完整的心跳周期；

第二步、确定第一步所得的心音信号x(k)的最大值x_max，以及对心音信号x(k)进行死区非线性变换的阈值C_l；

第三步、对第一步所得的心音信号x(k)作死区非线性变换，即将心音信号x(k)中小于阈值C_l的部分置为0，得心音信号

第四步、对第三步所得的心音信号进行自相关计算，得到自相关序列其中m∈{-N+1,…,0,…,N-1}；

第五步、确定对第四步得到的自相关序列作死区非线性变换的阈值C₂；

第六步、对第四步所得的自相关序列作死区非线性变换，即将自相关序列小于阈值C₂的部分置为0，得序列其中，为所得序列的一个峰值；

第七步、取第六步所得序列的一子序列将该子序列I按帧长L_w划分，并计算各帧绝对值均值，其中第i帧的绝对值均值为：

根据自相关函数的特性知：取峰值阈值C_p＝r_pEA(0)，其中，峰值阈值系数r_p∈[0.2,0.4]，自i＝1至各帧依次寻找高于峰值阈值C_p的峰值点，即满足条件：EA(i-1)<EA(i)、EA(i)>EA(i+1)且C_p<EA(i)，记与峰值点(0,EA(0))相邻的、满足上述条件的峰值点为(i_p,EA(i_p))，则对应于序列中的峰值点为其中N_T＝L_w×i_p；

第八步、求出心率为：每分钟60f_s/N_T次。

进一步地，第二步中，所述死区非线性变换的阈值C_l＝r_lx_max，r_l为预设比例系数，r_l∈[0.05,0.2]，x_max为心音信号x(k)的最大值。

进一步地，所述预设比例系数r_l值的大小能够根据心音信号x(k)中噪声幅值均方根的大小进行调整。

进一步地，第三步中，所述死区非线性变换的公式为：

进一步地，第四步中，所述自相关序列的计算方法为：

当m≥0时，否则其中，m∈{-N+1,…,0,…,N-1}。

进一步地，第五步中，所述死区非线性变换的阈值C₂＝aR_max，a为预设比例系数，a∈[0.6,0.9]，

进一步地，第六步中，所述死区非线性变换的公式为：

进一步地，第七步中，所述帧长L_w＝[r_wN/f_s]，其中帧长系数r_w∈[3,8]。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明针对心音信号，通过自相关处理以及非线性处理得到自相关序列，从自相关序列中计算出心率，只利用心音信号即可准确实现心率的计算，特别适用于医院和家庭电子听诊辅助诊疗系统中心率的监测。

2、本发明的计算方法代码简短且计算量小，对硬件实现的要求低，非常便于在可穿戴设备中实现。

附图说明

图1为本发明实施例一种基于心音信号自相关分析的心率计算方法的流程图。

图2为本发明实施例获取的一段心音信号示意图。

图3为本发明实施例获取的一段心音信号进行非线性变换预处理后得到的心音信号的示意图。

图4为本发明实施例对心音信号进行自相关处理后得到的序列的示意图。

图5为本发明实施例对序列作非线性变换后得到的序列的示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例：

本实施例提供了一种基于心音信号自相关分析的心率计算方法，流程如图1所示，包括以下步骤：

第一步、利用录音硬件设备获取心音信号其中N_s为正整数；

通常采样心音信号时硬件设备的采样频率一般共分为22.05KHz、44.1KHz、48KHz三个等级，因为心率计算对心音的保真要求不高，本实施例选取最低的采样率22.05KHz。为准确计算心率，需要处理至少两个完整周期的心音信号，同时为降低计算量，选取的信号不宜过长，本实施例选取三个完整周期的心音信号进行计算，按心率60次/分计算，可取心音信号的长度为3～5秒，若取采样频率为22.05KHz，并采集4秒钟心音信号，则N_s＝88200，所获得的心音信号如图2所示。

第二步、对第一步所得心音信号进行每隔c_us-1点抽取一点的欠采样，得到心音信号k∈{1,…,N}，其中即对作四舍五入的取整结果；当采样频率为22.05KHz时，可取c_us＝10，即 此时，本实施例所指的采样频率即为f_s＝2.205KHz；

作该步骤处理的原因在于：通常心音频率为1～1200Hz，对22.05KHz等高采样频率所得心音信号进一步欠采样至2000Hz左右亦不影响本发明的心率计算结果，但可极大地降低运算量。

第三步、求得第二步所得心音信号x(k)的最大值x_max＝max{x(k)}，从而确定非线性变换的阈值C_l＝r_lx_max，选取r_l＝0.1，对心音信号x(k)进行非线性处理，所得心音信号的示意图如图3所示。

第四步、对第三步所得的心音信号进行自相关计算，得到自相关序列序列的示意图如图4所示，其中m∈{-N+1,…,0,…,N-1}；

给定一个无限长序列s(n)，其自相关函数通常定义为：

若按该定义计算自相关函数，由于序列为有限长序列，且其下标最大值为N，故需要将下标n+m大于N的点置为0，计算得到的自相关函数将存在偏差，为此，本实施例的自相关计算采样如下方法：当m≥0时，否则其中，m∈{-N+1,…,0,…,N-1}。

第五步、确定对第四步得到的自相关序列作死区非线性变换的阈值C₂，C₂＝aR_max，a为预设比例系数，a∈[0.6,0.9]，

第六步、对第四步所得的自相关序列作死区非线性变换：

即将自相关序列小于阈值C₂的部分置为0，得序列序列的示意图如图5所示，其中，为所得序列的一个峰值，选取合适的阈值C₂即可消除自相关序列中影响心率判断的部分，即削弱甚至消除自相关函数中对应于真实心脏搏动周期的两个波峰之间的、由于第一第二心音波形相关而引起的虚假峰值。

第七步、取第六步所得序列的一子序列将该子序列I按帧长10划分，并计算各帧绝对值均值，其中第i帧的绝对值均值为：

根据自相关函数的特性知：取峰值阈值C_p＝0.3EA(0)，自i＝1至i<[N/10]各帧依次寻找第二周期的起始点，即满足大于阈值即为起始点，然后小于阈值的即可视为该周期的终点，记与峰值点(0,EA(0))相邻的、满足上述条件的峰值点为(i_p,EA(i_p))，则对应于序列中的峰值点为其中N_T＝10×i_p；本步骤中，以帧为单位进行计算的优点在于：自相关序列会因高频噪声的影响而存在毛刺，按帧计算可以减小因毛刺存在而对峰值点作出误判的可能。

第八步、求出心率为：每分钟60f_s/N_T次。

本实施例计算方法的原理是：心音为周期信号，故可通过其自相关序列的峰值间距确定心音信号的周期，即确定心率值。

第一心音是由于房室瓣关闭和室内血液冲击房室瓣，以及心室射出的血液撞击动脉壁引起振动而产生的。第二心音是由于主动脉瓣和肺动脉瓣关闭，血流冲击大动脉根部及心室内壁引起的振动而产生的。由于两种心音信号产生机制相近，两者对应的波形亦相似，导致心音信号的自相关序列在单个心音周期内会出现因第一心音第二心音之间的互相关而形成的多个峰值，单纯取自相关序列中相邻峰值间距计算心率将得到错误的结果。而另一方面，第一心音和第二心音波形尽管相似，但依然存在一定差别，故两者互相关所形成的峰值幅值较低，通过合理地设置阈值将可去掉。鉴于此，本实施例先后对心音信号和自相关序列作死区非线性处理，极大地减小因第一心音与第二心音两者互相关形成的峰值，提高心率计算准确性。

以上所述，仅为本发明专利较佳的实施例，但本发明专利的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明专利所公开的范围内，根据本发明专利的技术方案及其发明专利构思加以等同替换或改变，都属于本发明专利的保护范围。

Claims (8)

1.一种基于心音信号自相关分析的心率计算方法，其特征在于，所述方法包括以下步骤：

第一步、取长度为N点、采样频率为f_s的心音信号x(k)，其中k∈{0,…,N-1}，N为正整数，且N/f_s不小于3s，以保证所采集的心音信号x(k)至少包含三个完整的心跳周期；

第二步、确定第一步所得的心音信号x(k)的最大值x_max，以及对心音信号x(k)进行死区非线性变换的阈值C_l；

第三步、对第一步所得的心音信号x(k)作死区非线性变换，即将心音信号x(k)中小于阈值C_l的部分置为0，得心音信号

第四步、对第三步所得的心音信号进行自相关计算，得到自相关序列其中m∈{-N+1,…,0,…,N-1}；

第五步、确定对第四步得到的自相关序列作死区非线性变换的阈值C₂；

第六步、对第四步所得的自相关序列作死区非线性变换，即将自相关序列小于阈值C₂的部分置为0，得序列其中，为所得序列的一个峰值；

第七步、取第六步所得序列的一子序列将该子序列I按帧长L_w划分，并计算各帧绝对值均值，其中第i帧的绝对值均值为：

$EA\left(i\right)=\frac{1}{L_w}\underset{n=L_{w}i}{\overset{L_w(i+1)}{\&Sigma;}}\left|I\left(n\right)\right|,i\&Element;\{0,...,\&lsqb;\frac{N}{L_w}\&rsqb;\}$

根据自相关函数的特性知：取峰值阈值C_p＝r_pEA(0)，其中，峰值阈值系数rp∈[0.2,0.4]，自i＝1至各帧依次寻找高于峰值阈值C_p的峰值点，即满足条件：EA(i-1)<EA(i)、EA(i)>EA(i+1)且C_p<EA(i)，记与峰值点(0,EA(0))相邻的、满足上述条件的峰值点为(i_p,EA(i_p))，则对应于序列中的峰值点为其中N_T＝L_w×i_p；

第八步、求出心率为：每分钟60f_s/N_T次。

2.根据权利要求1所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：第二步中，所述死区非线性变换的阈值C_l＝r_lx_max，r_l为预设比例系数，r_l∈[0.05,0.2]，x_max为心音信号x(k)的最大值。

3.根据权利要求2所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：所述预设比例系数r_l值的大小能够根据心音信号x(k)中噪声幅值均方根的大小进行调整。

4.根据权利要求1所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：第三步中，所述死区非线性变换的公式为：

$\stackrel{\&OverBar;}{x}\left(k\right)=\left\{\begin{array}{ccc}1,& if& x\left(k\right)\&le;C_l\\ x\left(k\right),& if& x\left(k\right)>C_l\end{array}\right..$

5.根据权利要求1所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：第四步中，所述自相关序列的计算方法为：

当m≥0时，否则其中，m∈{-N+1,…,0,…,N-1}。

6.根据权利要求1所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：第五步中，所述死区非线性变换的阈值C₂＝aR_max，a为预设比例系数，a∈[0.6,0.9]，

7.根据权利要求1所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：第六步中，所述死区非线性变换的公式为：

${\stackrel{\&OverBar;}{R}}_{\stackrel{\&OverBar;}{xx}}\left(m\right)=\left\{\begin{array}{ccc}0,& if& R_{\stackrel{\&OverBar;}{xx}}\left(m\right)\&le;C_2\\ R_{\stackrel{\&OverBar;}{xx}}\left(m\right),& if& R_{\stackrel{\&OverBar;}{xx}}\left(m\right)>C_2\end{array}\right..$

8.根据权利要求1所述的一种基于心音信号自相关分析的心率计算方法，其特征在于：第七步中，所述帧长L_w＝[r_wN/f_s]，其中帧长系数r_w∈[3,8]。