基于心音香农能量包络的实时心率检测的研究和实现*

范能胜， 张怀岺

(广东医科大学生物医学工程学院，东莞 523808)

1 引 言

心脏是人体的重要器官，其节律性跳动是由窦房结产生的电信号引起的，称之为心电信号。心脏的一些疾病会引起心电波形的改变，但也有些疾病只会引起心音的改变，如瓣膜类的疾病。心脏在工作时，由于瓣膜的开启和关闭，血液对心脏内壁、瓣膜和血管的冲击会发出节律性的声音，称之为“心音”。一般认为有四种心音[1]，分别是第一心音(S1)、第二心音(S2)、第三心音(S3)和第四心音(S4)。S1代表心室收缩期的开始，其振幅大，持续时间较长；S2代表心室舒张期的开始，时间较短，振幅低；S3和S4不易听到。在分析心音信号时，S1和S2是最重要的信号来源。

听诊器是常用的心音探测设备，其结构简单，携带方便，是医生的首选。但通过听诊诊断心脏疾病主观性较强，医生需要经过长期的实践，经验越丰富的医生诊断阳性率越高；随着科技的发展，用特定的传感器将心音拾取出来，通过放大和模数转换后变成数字信号，利用分析算法从信号中提取有用信息，可以避免医生主观因素的影响。一些典型病例的心音信号在数字化后，有利于保存和传播，可以用来培训医生；此外，在日常的诊疗活动中，有时需要知道患者的心率变化情况，虽然通过心电信号计算心率最为准确，但是需要连接心电图机，操作不便，医生一般通过计数患者颈部或腕部动脉在单位时间内的搏动次数来计算心率，方法简单，但耗时长。

本研究介绍了一种利用心音信号来计算心率的方法，通过心音传感器拾取心音信号并将其数字化，按照一定算法提取心音包络信号，根据包络信号计算心率，理论上只要测量两个以上的心动周期就可以将心率计算出来，最终在听诊的同时完成心率的测量。

2 材料和设备

2.1 心音传感器

正常人的心音信号频率一般在5～600 Hz[2]，测试系统选用了安徽合肥华科电子技术研究所生产的HKY-06B+型心音传感器作为心音拾取设备，内含微音传感元件，可以采集心脏和其它体表动脉搏动信号，经过内部集成化信号处理电路处理，输出模拟电信号，幅值为0.5～1.5 V，频率响应为10～600 Hz，使用单电源供电，电压范围3～5 V，传感器实物见图1。

2.2 实验硬件平台

由于在实现算法时需要做浮点数运算，同时还需要显示波形和其它信息，故本研究选用型号STM32F103VET6的单片机作为核心处理器，该单片机具有32位数据处理宽度，内部资源丰富，处理速度快，并自带静态存储控制器(flexible static memory controller,FSMC)接口功能，方便连接LCD液晶显示；采用2.8英寸并带触摸控制功能的液晶显示器进行操作控制和显示波形，核心实验电路见图2。

图1 心音传感器Fig.1 Heart sound sensor

由于单片机自带的12位A/D转换器，可以满足实验要求，不必再外扩A/D转换电路，心音传感器输出的模拟信号经单片机的PC1引脚输入，前辍为LCD的网络标号连接液晶显示器，前辍为TP和SPI的网络标号连接触摸屏。

2.3 实验数据来源

根据心音传感器频率响应和奈圭斯特(Nyquist)准则，采样频率必须是信号最高频率的两倍以上才可以还原出原始信号，因此设置采样频率为1 250 Hz。将心音传感器固定在体表合适位置，采样8 s，获得10 000个心音采样数据，每个采样数据占用两个字节，将采样数据保存在一个数组中，并在采样结束后导入matlab软件，作为算法分析验证的实验数据。

3 算法

在对心音信号进行采集时，还可能引入其它的噪声信号，如呼吸音信号、传感器和胸壁摩擦产生的噪声信号等。处理算法的目的是衰减各种噪声信号的同时，将心音信号的包络提取出来作为分析处理使用。由于本研究的目的是利用心音信号来计算心率，所以需要从心音信号包络中准确识别心音峰值。根据相邻峰值之间的距离以及采样参数，计算瞬时心率值。有关心音信号处理的流程和算法，前人做了许多研究，包括心音信号的预处理、求取香农能量包络、寻找心音峰值等。结合这些研究成果和实际要求，我们采用如下方法对原始心音信号进行处理。

3.1 信号归一化预处理

为便于后续处理数据，需对原始输入信号进行归一化处理[4]，按式(1)进行:

(1)

图2 实验电路Fig.2 Main circuit

式中，xk为原始采样数据，max (|xi|)为输入序列数据绝对值的最大值，xnorm为归一化计算后的输出信号。

由于传感器的输出电压在0.5～1.5 V，所以输入信号振幅的峰-峰值为1.0 V，为了实时处理的需要，将式(1)中的max(|xi|)固定设为峰-峰值的一半，即0.5 V。在实时处理时，先对数字化的信号减去一个固定数值(基线下移)，处理过程见图3，然后再按照式(1)进行归一化运算。

3.2 计算香农能量包络

使用式(1)对原始信号进行归一化处理后，可以采用不同的方法来计算心音信号的包络，这些方法[4-8]包括：

绝对值：E=|x|

平方能量：E=x2

香农熵：E=-|x|log |x|

香农能量：E=-x2log (x2)

图3 原始信号处理过程(a).原始信号;(b).基线下移Fig.3 Process of original signal (a).original signal;(b).baseline down

图4显示了不同计算方法的情况下归一化原始信号幅值和包络输出值之间的关系曲线，由于对称性，负值部分被忽略了。

图4 不同包络侦测方法的比较Fig.4 The comparision of different envelope detection methods

由图4可知，绝对值的方法对输入信号的所有组成部分采用了相同的权重，难以将需要的信号凸显出来；平方能量对高强度信号部分给予了较大的权重，不利于低强度信号的识别；香农熵方法对低强度信号给予了较大的权重，衰减了高强度信号；香农能量的方法对中等强度信号给予了更多权重，对低强度信号比高强度信号衰减得更多，这样增强了中等偏高强度的信号。

对于心率测量只需关注第一心音S1即可，其它信号可以尽可能的衰减，根据文献[9]的研究,采用式(2)求取香农能量的平均值Eavg：

(2)

式中，N为平均滤波的数据点数，x(i)为输入信号，n为幂次。幂次n越高，计算所得的包络波形越平滑。在本研究中，经过验证，当n=4时可以取得较好的抑噪效果。

3.3 香农能量平均滤波的实现

用式(2)进行平均滤波时，N取值越大，得到的波形越平滑[10]，图5是N=201时的平均滤波效果，由图中可以看出，大部分的小幅波动被有效抑制，突出了第一心音信号。

图5 信号处理过程(a).归一化处理后的信号;(b).E=-x4log (x4);(c).Fig.5 Signal processing(a).normalization;(b).E=-x4log (x4);(c).

在实时处理时，可以设置一个环型数组用来保存图4(b)中计算出来的香农能量数据，数组的长度为平均滤波的点数，每次有新的数据进入时，让数组下标加一再保存。在此以7点平均滤波为例说明移动平均滤波的实现过程。

定义一个长度为7个元素的浮点数组Array，输入的数据按先后顺序依次保存在数组Array[0]～Array[6]，当第8个数据到来时覆盖最早保存的数据，将其存在Array[0]单元，以后的数据保存方法以此类推，环型数组见图6。

采用迭代递归的方式进行运算，环型数组每次存入的数值由式(3)给出：

(3)

此处的out(i)是本次运算的输出结果，out(i-1)为前一次运算的输出结果，m为平均滤波长度，in(i+p)表示在(i+p)刻的输入值，可理解为最新的输入值，in(i-q)表示在(i-q)时刻的输入值，in(i+p)-in(i-q)在环型数组中表示(i+p)刻的输入值覆盖掉(i-q)刻的输入值，当数组长度m=7时，由式(3)计算得到：

图6 环型数组Fig.6 Circular array

(4)

式(5)-式(8)按顺序演示了环型数组保存时的迭代过程，如果括弧内的值为负，表示该项的值为0。

out(0)=out(-1)+in(3)-in(-4)=in(3)

(5)

out(1)=out(0)+in(4)-in(-3)=in(3)+in(4)

(6)

out(2)=out(1)+in(5)-in(-2)=in(3)+in(4)+in(5)

(7)

…

out(7)=out(6)+in(10)-in(0)

(8)

in(10)表示新输入的第8个数据，in(10)-in(0)表示第8个输入的数据覆盖掉下标为0的数组元素，平均滤波的结果为out(i)/m。由于整个运算过程都是简单的四则运算，所以运算速度很快，完全满足信号实时处理的要求。

3.4 寻找峰值

要计算实时心率，必须寻找心音包络信号的峰值并确定它在数列中的坐标[11-12]，由于采样率已知，则只需知道相邻两个峰之间距离(数据点数)就可计算出瞬时心率，计算公式如下：

(9)

式中，HR表示心率，N为相邻两个心音峰值之间的数据点数，等于两个峰值的坐标之差，SR为采样率。

峰值指单位时间内的最大值，此“单位时间”称之为时间窗口。窗口的宽度需要根据采样率和心率来确定，本研究中将窗口宽度设置为0.32 S。同时考虑噪声的影响，还需设置一个阈值，在窗口范围内，只有大于这个阈值的最大值才能认为是心音峰值，见图7。

图7 寻找峰值(方框为时间窗口)Fig.7 Finding peak(red rectangle present time window)

在程序实现时，将新的输入值与前一次的值进行比较，前者大于后者，则进行最大值替换，保存数据坐标，并清零窗口计数器，反之，则不替换，并启动窗口计数器；窗口计数器达到窗口计数宽度，则认为当前的最大值为心音峰值，寻找峰值的算法如下：

//WinCnt: 窗口计数器；WINWIDTH: 窗口宽度；cnt：输入值计数器(数值坐标)

//CurrVal: 当前值; MaxVal: 最大值；MaxPos: 最大值坐标；ThresholdVal: 阈值If(WinCnt

{

if(CurrVal >MaxVal AND CurrVal>ThresholdVal)

{

MaxVal = CurrVal;

MaxPos = cnt;

WinCnt = 0;

}

else

{

WinCnt++;

}

}

else

{

return MaxPos;//返回峰值的坐标

}

根据以上算法，每找到一个心音峰值，返回其坐标，同时复位相关变量，然后根据式(9)即可计算出心率。实际应用时，可以每三次瞬时心率取平均值作为心率，既可满足实时显示的要求，又可以平滑计算引入的误差。

4 结果和讨论

心音信号因受到测量位置、呼吸音、环境等因素的影响，信号噪声较大，容易受干扰，理论上通过心音信号计算出来的心率比心电信号计算出来的心率误差大。为了验证算法的准确性，使用一台心电图机(型号ECG-101，深圳邦健公司)作为对照设备，实验设备和对照设备见图8。

图8 对照设备和实验设备(a).对照设备;(b)实验设备;(c).结果对照Fig.8 Reference device and experimental device(a).reference device;(b)experimental device;(c).results contrast

选择了六位测量对象，使用两种设备同时对每个人的心率进行测量，实验时测量的结果见表1。

由表1可以看出，与心电图机测量的心率对照，实验设备测量值的准确率达到97%以上。由于条件的限制，本次测试的对象均为健康成年人，未对有心脏疾病的患者进行测试，假如对这类患者进行测试，需要进一步完善算法，以提高测量的准确性。