技术领域
[0001] 本
发明涉及医疗器械与医学
信号处理领域，更具体地说，涉及一种基于心音自相关函数的心率计算算法。
背景技术
[0002] 心率是反映身体健康状况的重要指标之一，心率监测对于部分心脏病患者用药、家庭护理等均具有重要意义。
[0003] 随着
远程医疗的普及，
电子听诊器将进入更多社区医院和家庭，心音的远程听诊也将日益成为包括心脏病患者家庭护理在内的众多远程医疗项目的重要手段之一。相对于基于
脉搏测量心率等常用方法，利用心音信号计算心率准确度更高，使用方便且成本低廉。因此，基于心音计算心率势必成为今后心率监测的主要方式之一。
[0004] 基于心音信号计算心率通常可通过确
定心音包络峰值点间距离或心音自相关函数峰值点间距离来实现。但这些方法在实际应用中存在以下不利因素：
[0005] 1)包络和自相关函数均会因
呼吸音、摩擦音等因素的干扰而出现异常的波峰，导致心率计算错误；例如，即便是很浅的呼吸都可能因胸脯的轻微起伏使听诊器移位而出现摩擦音；而呼吸较急促时，这种干扰将更显著；
[0006] 2)有的心音同时存在第一和第二心音，由于第一和第二心音
波形相似，若直接应用传统的相关分析方法计算心率，容易得到虚假的周期值；虽然通过识别并分割第一和第二心音可克服这种因素的影响，但是相当部分心音第一和第二心音的分界并不明显，因此，不能确保通过去掉第二心音所致包络峰值或自相关函数峰值的形式来确定心率；
[0007] 3)
采样时间不可能过长；一般来说，听诊持续时间通常不超过半分钟，有时甚至更短，所以基于心音信号监测心率仅能根据数个心跳周期的心音信号来进行心率估算，此时现有的自相关函数的计算量极大，当采样
频率较高时，为在短时间内完成自相关函数计算，往往需配备计算能
力较强的
硬件设备，导致硬件成本偏高。
[0008] 故此，采用现有的方案和手段来实现基于心音估算心率，还存在很多困难。
发明内容
[0009] 本发明的目的在于克服
现有技术中的缺点与不足，提供一种基于心音自相关函数计算、可克服呼吸音等因素的干扰、只利用心音序列即可准确计算出心率值的心率计算算法。
[0010] 为了达到上述目的，本发明通过下述技术方案予以实现：一种基于心音自相关函数的心率计算算法，其特征在于：包括如下步骤：
[0011] 第一步，以
采样频率fs获取持续时间为t、长度为N的心音序列x(k)，其中， k＝1，...，N，N为正整数，N＝tfs；
[0012] 第二步，对心音序列x(k)作滤波处理得到序列
[0013] 第三步，以频率 对序列 欠采样得到序列s(i)， 其中，round()为四舍五入取整函数；
[0014] 以LT为半窗长求序列s(i)的包络En(m)；其中，
[0015] 第四步，设定移位序号为n；在移位序号的搜索范围内将包络直接进行自相关函数R
(n)计算；或者是将包络进行处理后，再在移位序号
的搜索范围内将处理后的包络进行自相关函
数R(n)计算；其中，c1为上限系数，c2为下限系数；
[0016] 第五步，初始化移位序号n、峰值点序号nP和峰值RP：令移位序号峰值点序号 峰值RP＝0；
[0017] 第六步，将移位序号n自加1；判断当前移位序号n和对应的自相关函数R(n) 数值大小：
[0018] 若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为差值下限；之后跳至第七步；
[0019] 否则直接跳至第七步；
[0020] 第七步，判断当前移位序号n的大小：若 则返回至第六步；否则跳转至第八步；
[0021] 第八步，计算心率值： 单位为次/分钟。
[0022] 本发明计算算法的原理是：根据心音信号的特点，通过限定自相关函数峰值搜索范围、增加滤波等方式抑制呼吸音、摩擦音等因素引起的干扰。在搜索范围内计算自相关函数，通过搜索自相关函数的最大峰值点，并采用出现最大峰值点时的移位距离来确定心率周期。该方式的好处是：第一心音和第二心音产生机制相近，两者波形亦相似，导致心音自相关函数在单个心音周期内会出现因第一、第二心音之间的互相关而形成的多个峰值。但由于第一心音和第二心音波形存在一定差别，其因互相关所形成的峰值幅值较低。故此，本发明第六步中寻找R(n)的最大峰值点，并在第八步通过该峰值点的序号计算心率，而不是简单地通过自相关函数相邻峰值点间距离计算心率，以避免因第一、第二心音波形相似造成的心率计算错误。由于周期信号的自相关函数在基频的倍频处会出现峰值；若心音信号
质量不佳，极可能出现自相关函数在倍频处的峰值高于基频处的峰值，而导致心率计算出现错误。故此，本发明第六步中采用 |n-ainP|＞nc，ai＝2，3，..作为峰值点选择的必要条件之一，可避免该错误的发生。此外本发明计算算法限定自相关函数峰值搜索范围，既可准确计算心率值，又可降低计算算法计算量，快速得出结果，有利于降低
硬件实现要求和生产成本。
[0023] 当采样频率较高时，直接计算实测听诊信号自相关函数的计算量极大；本发明在第三步通过欠采样和包络计算来降低所需计算自相关函数的序列长度。以采样频率为22050Hz并取 且LT＝10为例，采用欠采样和包络计算，则处理前后所需计算自相关函数的序列的长度比高达220∶1；因此本发明计算算法计算量大大降低，从而降低对硬件实现的要求。
[0024] 该计算算法既可以在小型可穿戴设备和小型电子听诊器中实现，又可以在大型电子听诊设备中实现。
[0025] 优选地，所述第一步中，持续时间为t的取值范围为：2.5s＜t＜4s。
[0026] 正常心率值范围为45＜HR＜180次/分钟，而呼吸频率一般不会超过40次/ 分钟，呼吸频率超过24次/分钟已属呼吸急促。故此，本发明第一步中选取持续时间落于2.5秒到4秒之间的心音序列来计算心率，可保证所选取的信号段至少包含两次以上心跳，满足自相关函数计算心率所需的心跳数；同时，也保证了在所选信号段内不会出现两次或两次以上呼吸，有效避免呼吸音对于心音周期计算的影响。
[0027] 优选地，所述第二步中，对心音序列x(k)作滤波处理得到序列 是指：对心音序列x(k)作截止频率为fc的低通滤波得到序列 以抑制呼吸音和摩擦音对心音信号的干扰。
[0028] 优选地，所述截止频率fc的取值范围为：fc∈[150，250]Hz；低通滤波是指二阶以上的低通滤波。相对于心音信号，呼吸音信号和摩擦音信号在高于250Hz 的频段中成份占比较高。而决定心音包络的成份则主要集中于心音信号的低频段[0，250]Hz，尤其是[0，150]Hz频段内。即仅根据心音信号在[0，150]Hz频段内的成份即可有效确定心率。故此，本发明第二步中对心音序列作截止频率为 fc∈[150，250]Hz的低通滤波，有助于降低呼吸音和摩擦音的幅值，抑制其对于心音包络计算以及对于心率计算的干扰。
[0029] 优选地，所述第三步中，频率 的取值范围为：
[0030] 半窗长LT的取值范围为：
[0031] 包络En(m)的计算方法为：
[0032]
[0033] 或者是
[0034] 优选地，所述第四步中，上限系数c1和下限系数c2的取值范围分别为： c1∈[0.25，0.27]，c2∈[1.2，1.5]。
[0035] 考虑到心率值的范围为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。
[0036] 其中一种方案是：所述第四步中，在移位序号的搜索范围内将包络直接进行自相关函数R
(n)计算是指：计算包络En(m)的自相关函数R(n)：
[0037]
[0038] 其中，
[0039] 另一种方案是：所述第四步中，将包络进行处理后，再在移位序号的搜索范围内将处理后的包络进行自相关函
数R(n)计算，是指：
[0040] 对包络En(m)作归一化处理： 其中，max(En)为En(m) 的最大值；之后在移位序号 的搜索范围内计算包
络 的自相关函数R(n)。
[0041] 优选地，所述在移位序号 的搜索范围内计算包络 的自相关函数R(n)，是指：
[0042]
[0043] 其中，
[0044] 与现有技术相比，本发明具有如下优点与有益效果：
[0045] 1、本发明计算算法可克服呼吸音、摩擦音以及第一、第二心音波形相似等因素的干扰，只利用心音信号即可实现心率计算，准确程度高；
[0046] 2、本发明计算算法计算量小，对硬件实现的要求低，非常便于在小型可穿戴设备和小型电子听诊器中实现，同时本发明计算算法也特别适用于在医院和家庭用电子听诊辅助诊疗系统中实现心率监测。
附图说明
[0047] 图1是
实施例一基于心音自相关函数的心率计算算法的
流程图；
[0048] 图2(a)～图2(e)是本发明实施例中计算算法各步骤所得结果；
[0049] 图3(a)～图3(h)是本发明实施例中滤波前后呼吸音、心音信号曲线及相应
频谱；
[0050] 图4是实施例二基于心音自相关函数的心率计算算法的流程图。
具体实施方式
[0051] 下面结合附图与具体实施方式对本发明作进一步详细的描述。
[0052] 实施例一
[0053] 本实施例基于心音自相关函数的心率计算算法，其工作流程如图1所示，包括如下步骤：
[0054] 第一步，以采样频率fs获取持续时间为t、长度为N的心音序列x(k)，其中 k＝1，...，N，N为正整数，N＝tfs；持续时间为t的取值范围优选为： 2.5s＜t＜4s；
[0055] 第二步，对心音序列x(k)作滤波处理得到序列 优选是对心音序列x(k) 作截止频率为fc的低通滤波得到序列 截止频率fc的取值范围优选为：fc∈[150，250]Hz；低通滤波可采用现有低通滤波方法，优选采用二阶以上低通滤波；采用其它可抑制呼吸音和摩擦音干扰的滤波处理也可；
[0056] 第三步，以频率 对序列 欠采样得到序列s(i)， 其中，以LT为半窗长求序列s(i)的包络 En
(m)；其中， round()为四舍五入取整函数：
[0057] 频率 的取值范围优选为：
[0058] 半窗长LT的取值范围优选为：
[0059] 包络En(m)的计算方法优选为：
[0060] 或者是
[0061] 第四步，设定移位序号为n ；将包络进行处理后 ，再在移位序号的搜索范围内将处理后的包络进行自相关函
数R(n)计算；其中，c1为上限系数，c2为下限系数；上限系数c1和下限系数c2的取值范围分别优选为：c1∈[0.25，0.27]，c2∈[1.2，1.5]；
[0062] 将包络进行处理是指：对包络En(m)作归一化处理： 其中，max(En)为En(m)的最大值；之后在移位序号
的搜索范围内计算包络 的自相关函数R(n)：
[0063]
[0064] 实际应用中自相关函数也可以采用其它计算公式。
[0065] 第五步，初始化移位序号n、峰值点序号nP和峰值RP：令移位序号峰值点序号 峰值RP＝0；
[0066] 第六步，将移位序号n自加1；判断当前移位序号n和对应的自相关函数R(n) 数值大小：
[0067] 若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]；之后跳至第七步；
[0068] 否则直接跳至第七步；
[0069] 第七步，判断移位序号n的大小：若 则返回至第六步，否则跳转至第八步；
[0070] 第八步，计算心率值： 单位为次/分钟。
[0071] 图2(a)是实测听诊信号，图2(b)是低通滤波后信号，图2(c)是欠采样所得信号，图2(d)是欠采样信号的包络，图2(e)包络的自相关函数。
[0072] 图3(a)是兼有心音和呼吸音的信号，图3(b)是兼有心音和呼吸音的信号频谱，图3(c)是兼有心音和呼吸音的信号低通滤波后曲线，图3(d)是兼有心音和呼吸音的信号低通滤波后频谱，图3(e)是滤波前心音信号，图3(f)是滤波前心音频谱，图3(g)是低通滤波后心音信号，图3(h)是低通滤波后心音频谱。
[0073] 本发明计算算法的原理是：
[0074] 如图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)含呼吸音和摩擦音的
片段在低通滤波前后曲线的幅值)，抑制其对于心音包络计算以及对于心率计算的干扰。
[0075] 正常心率值范围为45＜HR＜180次/分钟，而呼吸频率一般不会超过40次/ 分钟，呼吸频率超过24次/分钟已属呼吸急促。故此，本发明第一步中选取持续时间落于2.5秒到4秒之间的心音序列来计算心率，可保证所选取的信号段至少包含两次以上心跳，满足自相关函数计算心率所需的心跳数；同时，也保证了在所选信号段内不会出现两次或两次以上呼吸，有效避免呼吸音对于心音周期计算的影响。
[0076] 考虑到心率值的范围为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)中所示峰值点搜索范围)。该上下限的选择有利于排除呼吸音的影响，避免将呼吸频率误作为心率。
[0077] 周期信号的自相关函数在基频的倍频处会出现峰值。若心音信号质量不佳，极可能出现自相关函数在倍频处的峰值高于基频处的峰值(见图2(e)中所示心音包络自相关函数倍频处波峰)，而导致心率计算出现错误。故此，本发明第六步中采用|n-ainP|＞nc，ai＝2，3，..作为峰值点选择的必要条件之一，可避免该错误的发生。
[0078] 第一心音和第二心音产生机制相近，两者波形亦相似，导致心音自相关函数在单个心音周期内会出现因第一、第二心音之间的互相关而形成的多个峰值(见图2(e)中所示因第一、第二心音波形相似所致心音包络自相关函数峰值点)。但由于第一心音和第二心音波形存在一定差别，其因互相关所形成的峰值幅值较低。故此，本发明第六步中寻找R(n)的最大峰值点，并在第八步通过该峰值点的序号计算心率，而不是简单地通过R(n)相邻峰值点间距离计算心率，以避免因第一、第二心音波形相似造成的心率计算错误。
[0079] 当采样频率较高时，直接计算实测听诊信号自相关函数的计算量极大；本发明第三步通过包络计算或者是欠采样和包络计算来降低所需计算自相关函数的序列长度；其中，通过欠采样和包络计算可进一步降低所需计算自相关函数的序列长度。以采样频率为22050Hz并取 且LT＝10为例，采用S1 步到S3步中的欠采样和包络计算，则处理前后所需计算自相关函数的序列的长度比高达220：1；因此本发明计算算法计算量大大降低，从而降低对硬件实现的要求。
[0080] 下面以具体示例对本发明计算算法进行说明：
[0081] 一种基于心音自相关函数的心率计算算法，包括如下步骤：
[0082] 第一步：取长度为N＝77175点、采样频率为fs＝22050Hz的心音序列x(k)，其中k＝1，...，N，持续时间t＝N/fs，显然所取持续时间t满足条件2.5s＜t＜4s，即所分析的心音序列持续时间落于2.5秒到4秒之间。
[0083] 应说明的是，心率计算对心音的保真程度要求不高，采样频率可行的取值范围可宽至fs∈[2000，48000]Hz。
[0084] 第二步：对x(k)作截止频率为fc＝200Hz、二阶低通滤波得滤波后序列[0085]
[0086] 其中，低通滤波参数采用巴特沃斯低通
滤波器设计计算算法，由截止频率fc＝200Hz和采样频率fs＝22050Hz共同决定。
[0087] 第三步，取频率 则 以 对 欠采样得 其中，i＝1，...，3500。
[0088] 取LT＝10，求s(i)的包络 其中，
[0089] 第四步，对En(m)作归一化处理：
[0090] 分别取c1＝0.25，c2＝1.5，则计算 的自相关函数：
[0091]
[0092] 其中，
[0093] 第五步，分别取 峰值点序号峰值RP＝0。
[0094] 第六步，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；之后跳至第七步；否则直接跳至第七步。
[0095] 第七步，若 则返回到第六步，否则跳转至第八步。
[0096] 第八步，计算心率得 次/分钟。
[0097] 实施例二
[0098] 本实施例基于心音自相关函数的心率计算算法，其工作流程如图4所示，与实施例一的区别在于：本实施例中，第四步，在移位序号的搜索范围内将包络En(m)直接进行自相关
函数R(n)计算。
[0099] 具体地说，是指：计算包络En(m)的自相关函数R(n)：
[0100]
[0101] 其中，
[0102] 本实施例的其余步骤与实施例一相同。
[0103] 上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
