树莓派pico+MAX30102实现心率脉搏检测——数据处理

前言：

        经过上文我们已经能够通过I2C配置传感器模块，并且将光反射强度的原始数据得到。那么这篇文章便是基于前文，对原始数据进行处理，得到我们想要的脉搏和心率。上一篇文章：链接同时为了方便大家理解我还对之前的文章进行一些补充，比如之前只大概讲了一下脉搏检测模块的原理，以及具体的实现。起主要作用的便是配置寄存器，所以在这里给大家推荐一篇文章个人认为该博主对寄存器的讲解很全面。max30102寄存器

正文：

        对于对原始数据的处理都是用得一套算法。主要便是通过计算公式得出结果，但是再利用公式之前我们会对数据进行处理，比如平滑信号，筛除不合理数据等等一系列操作。其实个人认为也可以根据计算公式直接将原始数据进行计算得到结果不过这样误差可能会很大。最后对于代码绝大部分我都有注释，大家可以结合我的注释分析。

# -*-coding:utf-8
# 数据采集频率(25HZ,0.04秒一个数据)
SAMPLE_FREQ = 25
# 四个数据一组，平滑信号
MA_SIZE = 4
# 红外数据与红光数据的数据长度最大值
BUFFER_SIZE = 100
def answer(ir_data, red_data):
    # 求平均值
    num = len(ir_data)
    ir_mean = sum(ir_data) // num

    # 所有的ir_data数据减去平均值在取反
    x = [0]*num
    for i in range(num):
        x[i] = -1*(ir_data[i]-ir_mean)

    # 4个数据一组，对每一组求平均值，信号波平滑处理
    for i in range(len(x) - MA_SIZE):
        x[i] = sum(x[i:i+MA_SIZE]) // MA_SIZE

    # 阈值合理化
    n_th = sum(x)//len(x)
    n_th = 30 if n_th < 30 else n_th  # min allowed
    n_th = 60 if n_th > 60 else n_th  # max allowed

    # 得到合适的波峰与波谷数值数组所对应的下标，n_peaks一定是偶数（当len( ir_valley_locs)偶数则n_peaks与其相等，否则减1）
    ir_valley_locs, n_peaks = find_peaks(x, BUFFER_SIZE, n_th, 4, 15)
    peak_interval_sum = 0
    if n_peaks >= 2:
        for i in range(1, n_peaks):
            # 求波峰与波谷之间的差值的累加和，即心脏收与缩时光反射强度差的所有和
            peak_interval_sum += ir_valley_locs[i] - ir_valley_locs[i-1]
        peak_interval_sum = int(peak_interval_sum / (n_peaks - 1))   # 求波峰的平均时间间隔(平均滤波)
        hr = int(SAMPLE_FREQ * 60 / peak_interval_sum)  # 心率计算公式
        hr_valid = True
    else:
        # 标记错误的结果
        hr = -999 
        hr_valid = False

    # find precise min near ir_valley_locs (???)
    exact_ir_valley_locs_count = n_peaks

    # FIXME: needed??
    for i in range(exact_ir_valley_locs_count):
        # 下标位置不可能大于100
        if ir_valley_locs[i] > BUFFER_SIZE:
            spo2 = -999 
            spo2_valid = False
            return hr, hr_valid, spo2, spo2_valid

    i_ratio_count = 0
    ratio = []


    red_dc_max_index = -1
    ir_dc_max_index = -1
    for k in range(exact_ir_valley_locs_count-1):
        red_dc_max = -16777216
        ir_dc_max = -16777216
        if ir_valley_locs[k+1] - ir_valley_locs[k] > 3:
            # 找出小标是ir_valley_locs[k], ir_valley_locs[k+1]之间的最大red,ir.
            for i in range(ir_valley_locs[k], ir_valley_locs[k+1]):
                if ir_data[i] > ir_dc_max:
                    ir_dc_max = ir_data[i]
                    ir_dc_max_index = i
                if red_data[i] > red_dc_max:
                    red_dc_max = red_data[i]
                    red_dc_max_index = i
            # 求心率公式为:R=((ir_max+ir_min)(red_max-red_min))/((red_max+red_min)(ir_max-ir_min));
            red_ac = int((red_data[ir_valley_locs[k+1]] - red_data[ir_valley_locs[k]]) * (red_dc_max_index - ir_valley_locs[k]))
            red_ac = red_data[ir_valley_locs[k]] + int(red_ac / (ir_valley_locs[k+1] - ir_valley_locs[k]))
            red_ac = red_data[red_dc_max_index] - red_ac 

            ir_ac = int((ir_data[ir_valley_locs[k+1]] - ir_data[ir_valley_locs[k]]) * (ir_dc_max_index - ir_valley_locs[k]))
            ir_ac = ir_data[ir_valley_locs[k]] + int(ir_ac / (ir_valley_locs[k+1] - ir_valley_locs[k]))
            ir_ac = ir_data[ir_dc_max_index] - ir_ac  

            nume = red_ac * ir_dc_max
            denom = ir_ac * red_dc_max
            if (denom > 0 and i_ratio_count < 5) and nume != 0:
                ratio.append(int(((nume * 100) & 0xffffffff) / denom))
                i_ratio_count += 1


    ratio = sorted(ratio) 
    mid_index = int(i_ratio_count / 2)

    ratio_ave = 0
    if mid_index > 1:
        ratio_ave = int((ratio[mid_index-1] + ratio[mid_index])/2)
    else:
        if len(ratio) != 0:
            ratio_ave = ratio[mid_index]

    # why 184?
    # print("ratio average: ", ratio_ave)
    if 2 < ratio_ave < 184:
        # -45.060 * ratioAverage * ratioAverage / 10000 + 30.354 * ratioAverage / 100 + 94.845
        # 心率计算公式:SpO2 = ((-45.060)RR + 30.354 * R + 94.845).
        spo2 = -45.060 * (ratio_ave**2) / 10000.0 + 30.054 * ratio_ave / 100.0 + 94.845
        spo2_valid = True
    else:
        spo2 = -999
        spo2_valid = False

    return hr, hr_valid, spo2, spo2_valid

# 寻找波峰与波谷
def find_peaks(x, size, min_height, min_dist, max_num):
    """
    Find at most MAX_NUM peaks above MIN_HEIGHT separated by at least MIN_DISTANCE
    """
    ir_valley_locs, n_peaks = find_peaks_above_min_height(x, size, min_height, max_num)
    ir_valley_locs, n_peaks = remove_close_peaks(n_peaks, ir_valley_locs, x, min_dist)

    n_peaks = min(n_peaks, max_num)

    return ir_valley_locs, n_peaks

# 寻找波峰与波谷
def find_peaks_above_min_height(x, size, min_height, max_num):
    """
    Find all peaks above MIN_HEIGHT
    """
    i = 0
    n_peaks = 0  # 平坦峰的个数，即波峰和者波谷的总个数
    ir_valley_locs = []  # 波峰或者波谷的下标（位置）
    while i < size - 1:
        if x[i] > min_height and x[i] > x[i-1]:  # find the left edge of potential peaks
            n_width = 1
            # i到n_width之间值相等，平坦峰，没有波峰与波谷
            while i + n_width < size - 1 and x[i] == x[i+n_width]:
                n_width += 1
            if x[i] > x[i+n_width] and n_peaks < max_num:  # find the right edge of peaks
                # ir_valley_locs[n_peaks] = i
                ir_valley_locs.append(i)
                n_peaks += 1  # original uses post increment
                i += n_width + 1
            else:
                i += n_width
        else:
            i += 1
    return ir_valley_locs, n_peaks

# 删除不满足条件的波峰与波谷（波峰与波谷的差值不满足条件）
#  即是删除不准确的测量数据
def remove_close_peaks(n_peaks, ir_valley_locs, x, min_dist):
    # 按照列表x中的下标的位置排序(从大到小)
    sorted_indices = sorted(ir_valley_locs, key=lambda i: x[i],reverse=True)
    i = -1
    while i < n_peaks:
        old_n_peaks = n_peaks
        n_peaks = i + 1
        j = i + 1
        # 快慢双指针找到满足合适范围的波峰波谷，筛选提高准确度
        while j < old_n_peaks:
            if i != -1:
                n_dist = sorted_indices[j] - sorted_indices[i]
            else:
                n_dist = sorted_indices[j] + 1
            if n_dist > min_dist or n_dist < -1 * min_dist:
                sorted_indices[n_peaks] = sorted_indices[j]
                n_peaks += 1
            j += 1
        i += 1

    sorted_indices[:n_peaks] = sorted(sorted_indices[:n_peaks])
    return sorted_indices, n_peaks

由于是五一假期期间创作的，所以介绍的也比较粗糙，对于心率计算公式或者其它细节如何实现的网上也很详细。祝大家五一快乐！最后希望各位使用了我代码的友友们能够三连走一波。