1、摘要
在上一篇文章中：具体介绍了如何用DMA和ADC采集心电数据并上传到上位机。在采集到心电数据后，紧接着微型处理器需要做的就是将数据转换成为可分析的生理指标——心率。顾名思义，心率指的就是一分钟内的心跳次数。计算心率的传统方法是计时一分钟，测量出一分钟内产生了多少次脉搏。这种做法的缺点十分明显，一是效率极低，每次更新心率测量的时间间隔非常长，二是容易解析失误，可能会忽略某些关键数据点。为了解决这种方式的弊端，本篇文章将介绍一种心率解析算法：动态阈值算法，用以从采样的心电信号中解析出实时心率，并且将实时心率大小上传到上位机进行观测。
2、算法核心思想与心率信号的有效特征点
算法的核心思想即为测量相邻的两次脉搏的时间间隔，再用心率的单位时间（一分钟）除以这个时间间隔得出心率。基于这种心率解析算法可以得到心率信号中的两个有效特征点
-
IBI(Inter-Beat Intervals)：心搏间隔，指的是两次有效脉搏之间的时间间隔。单位：ms。在心电图中的标记显示如下图（有多种选取标准，效果相同）：
-
BPM(Beat Per Minute)：每分钟节拍数，这里特指心率值，与IBI之间的转换关系为BPM = 60000 / IBI；
3、动态阈值算法分思路分析
无论是采用传统方法以及本篇文章中所提到的动态阈值算法，只有识别出一个脉搏，才能数出一分钟内脉搏数或者计算两个相邻的脉搏之间的时间间隔。我们立刻想到的方法可能会是从采集的电压波形数据中识别一个波峰/波谷，人为地在程序中设定一个阈值，当读取的电压数据超过这个阈值时系统便可认为检测到了一次脉搏。这种方法在应对输出的电压波形较为齐整时是可行的（波峰和波谷的值较为稳定且波形的频率分量较少），但是在心率的测算中以此种方式的测量会与实际的心率产生较大的误差，原因有以下两点：
- 心电波形的电压输出是不稳定的，振幅有大有小且极其容易发生变化；
- 心电波形中含有较为多的频率分量，采用常规方法进行计算容易识别到无效波形，导致结果与真实值间存在偏差。
既然固定阈值的方法不可取，那自然想到改变阈值——根据信号振幅调整阈值以适应不同信号的波峰检测。通过对一个周期内的信号多次采样，得出信号的最高与最低电压值，由此算出阈值，再用这个阈值对采集的电压值进行判断，考虑是否为波峰。也就是说电压信号的处理分两步，首先计算出参考阈值，然后用阈值对电压信号进行判定，识别一个波峰，实现思路如下图所示：
上面得出的是一段有效波形，而计算IBI只需要一个点。需要从一段有效信号上选取一个特征点，这个特征点的性质是：每一个有效波形仅有一个特征点存在，这个特征点代表着一个有效脉搏，只要能识别到这个特征点，就能在一个脉冲到来时触发任何动作（如计数、记录特征点时间）。通过记录两个相邻特征点的时间并求差值，计算IBI便水到渠成，在本节中我们将选取心电波形上升到振幅的一半作为特征点，我们可以捕获这个特征点并记录捕获时间进而计算IBI。
4、算法整体实现
根据上述分析，可得动态阈值算法的整体思路如下：
- 缓存一个波形周期内的多次采样值，求出最大最小值，计算出振幅中间值作为信号判定阈值；
- 通过把当前采样值和上一采样值与阈值作比较，寻找到信号上升到振幅中间位置的特征点，记录当前时间；
- 寻找下一个特征点并记录时间，算出两个点的时间差值，即相邻两次脉搏的时间间隔IBI；
- 由IBI计算心率值BPM；
算法代码实现如下：
/**
*@brief 寻找当前数组的最大值
*@attention 无
*@param data[] 16位无符号整型数组
*@retval 当前数组最大值
*/
static uint16_t MaxElement(uint16_t data[])
{
uint8_t i;
uint16_t max = data[0];
for(i = 1; i < Size; i++)
{
if(data[i] >= max)
{
max = data[i];
}
}
return max;
}
/**
*@brief 寻找当前数组的最小值
*@attention 无
*@param data[] 16位无符号整型数组
*@retval 当前数组最小值
*/
static uint16_t MinElement(uint16_t data[])
{
uint8_t i;
uint16_t min = data[0];
for(i = 1; i < Size; i++)
{
if(data[i] <= min)
{
min = data[i];
}
}
return min;
}
/**
*@brief 计算心率函数
*@attention 在主函数中以15Hz频率工作
*@param data 16位无符号整形数据
*@retval 无
*/
void Rate_Calculate(uint16_t data)
{
static uint8_t index; //当前数组下标
static uint8_t pulseflag; //特征点数目
static uint16_t Input_Data[Size]; //缓存数据数组
static uint16_t predata, redata; //一个波形周期内读取上一次值、当前值
uint16_t Maxdata, Mindata; //最大值，最小值
static uint16_t Middata; //平均值
static uint16_t range = 1024; //幅度
static uint8_t prepulse, pulse; //衡量指标值，用以标记特征点
static uint16_t time1, time2, time; //检测时间
uint16_t IBI, BIM; //两个特征点间隔时间，心率大小
predata = redata; //保存上一次的值
redata = data; //读取当前值
if(redata - predata < range) //滤除突变噪声信号干扰
{
Input_Data[index] = redata; //填充缓存数组
index++;
}
if(index >= Size)
{
index = 0; //覆盖数组
Maxdata = MaxElement(Input_Data);
Mindata = MinElement(Input_Data);
Middata = (Maxdata + Mindata) / 2; //更新参考阈值
range = (Maxdata - Mindata) / 2; //更新突变阈值
}
prepulse = pulse; //保存当前脉冲状态
pulse = (redata > Middata) ? 1 : 0; //采样值大于中间值为有效脉冲
if(prepulse == 0 && pulse == 1) //当检测到两个有效脉冲时为一个有效心率周期
{
pulseflag++;
pulseflag = pulseflag % 2;
if(pulseflag == 1) //检测到第一个有效脉冲
{
time1 = time; //标记第一个特征点检测时间
}
if(pulseflag == 0) //检测到第二个有效脉冲
{
time2 = time; //标记第二个特征点检测时间
time = 0; //清空计时
if(time1 < time2)
{
IBI = (time2 - time1) * Period; //计算两个有效脉冲的时间间隔（单位：ms）
BIM = 60000 / IBI; //1分钟有60000ms
/**限制BIM最高/最低值*/
if(BIM > 200)
{
BIM = 200;
}
if(BIM < 30)
{
BIM = 30;
}
printf("Currented Heart Rate:%d ", BIM);
}
}
}
time++;
}
注意事项：
- Size的大小合适选取，原则是尽可能使采样点覆盖一个完整的波形，但又不适宜过大，导致过多的数据点被缓存，降低准确率，同时要控制采样时间Period，不宜采样时间过快或过慢，这里给出两组合适的方案：Size = 66 Period = 15 or Size = 40 Period = 25。
- 传感器与手指之间的接触要保持稳定，频繁变动接触点及按压力度可能会对测量结果产生影响。
5、算法实现效果
按照正确的操作测量得到的心率数据在上位机中显示如下：
所测量的数据会自动保存在指定路径的txt文档中：
6、小结
本节介绍了采用动态阈值算法获取心率值的核心思想以及实现步骤，较为精确地测量出了实时心率数据。下一节将会介绍上位机的具体实现。
