使用PPG信号计算脉率和血氧

Github代码地址：https://github.com/hzzhangqf0558/SPO2_HR

- PPG信号简介

脉搏波是心脏的搏动（振动）沿动脉血管和血流向外周传播而形成的。心脏是一个持续不断的振源，心室收缩时，血液快速射入主动脉致其基部压力骤增而膨胀；心室暂未射血时，主动脉基部压力下降，管壁弹性回缩，则恢复至原来位置。如此，主动脉管壁就因心室的舒缩而有节律地受迫振动。这种振动能沿弹性血管向末梢传播而形成脉搏波（横波）。在传播过程中，必然要受到血管壁弹性、血管周围阻力、血液黏度等因素的影响，从而使脉搏波形发生变化。

一个完整的脉搏波波形有 A，B，C，D 4 个重要特征点，其包含升支和降支。如图 1 所示，A 称作主波，B 称作潮波，C 称作重搏波峰，D 称作重搏波谷，OA 是主波上升支，OO‘是脉搏波周期，这个脉搏波曲线的变化体现了人体的生理病理变化。

OA 段波形简介
脉搏波形的上升支，是由于在动脉壁忽然扩展的时候，左心室射血，动脉血压快速升高，形成管壁扩张。O 点是心脏射血期的开始点，反映收缩期末血管内的压力和容积，A 点是主动脉压力最高点，反映动脉内压力与容积的最大值。心输出量、射血速度和阻力是影响上升支的幅度和斜率的主要因素，一般用脉搏波波形上升的斜率和幅度表示，心输出量较多、射血速度较快、主动脉弹性越小、阻力越小，则斜率较大，波幅较高；反之，则斜率较小，波幅较低。

AD 段波形简介
下降支的前段，是由于在心室射血的后面阶段，射血速度开始降低，造成主动脉流向周围的血量大于流进主动脉的血量，大动脉由扩张变成回缩，动脉血压逐步变低这个过程造成的。B 点是左心室射血停止点，是反射波的波峰点，也叫潮波波峰，反映动脉血管的张力、顺应性和外周阻力的大小。D 点是潮波波谷点，即心脏收缩与舒张的分界点，主要反映外周阻力的大小。

DO’段波形简介
下降支的后段，也叫重搏波，是由于心室扩张，动脉血压不断降低，主动脉内血液向心室方向反流形成的。反映主动脉的功能状况，血管弹性和血液流动状态。

PPG中的心率

(FS * 60 /n_peak_interval_sum)
计算出脉率-- beats per minutes 一分钟采集6000个样本，除以峰峰之间平均样本数
FS 频率；n_peak_interval_sum 一段PPG信号中的峰峰之间平均间隔点数；

PPG中的血氧

R = Red(ac)/Red(dc) / Ir(ac)/Ir(dc). Red 红光;Ir红外光; ac 交流分量;dc直流分量;
Spo2 = A*R*R + B * R + c (n阶多项式)

代码(已详细注释):

algorithm.h

/** \file algorithm.h ******************************************************
*
*
* --------------------------------------------------------------------
*Default, this code runs in the x86/x64/linux system.
If this code runs into embeded system,
* this code should follow the following naming conventions:
*
*\n char              ch_pmod_value
*\n char (array)      s_pmod_s_string[16]
*\n float             f_pmod_value
*\n int32_t           n_pmod_value
*\n int32_t (array)   an_pmod_value[16]
*\n int16_t           w_pmod_value
*\n int16_t (array)   aw_pmod_value[16]
*\n uint16_t          uw_pmod_value
*\n uint16_t (array)  auw_pmod_value[16]
*\n uint8_t           uch_pmod_value
*\n uint8_t (array)   auch_pmod_buffer[16]
*\n uint32_t          un_pmod_value
*\n int32_t *         pn_pmod_value
*
* -------------------------------------------------------------------------
*/#ifndef ALGORITHM_H_
#define ALGORITHM_H_
//#include "mbed.h"#define true 1
#define false 0
#define FS 25//频率
#define BUFFER_SIZE  (FS* 12)//一段信号长度
#define HR_FIFO_SIZE 7//bufsize
#define MA4_SIZE  4// DO NOT CHANGE
#define HAMMING_SIZE  5// DO NOT CHANGE
#define min(x,y) ((x) < (y) ? (x) : (y))const double auw_hamm[HAMMING_SIZE] = { 40.9600000000000,276.480000000000,512,276.480000000000,40.9600000000000 }; //Hamm=  long16(512* hamming(5)');
//uch_spo2_table is computed as  -45.060*ratioAverage* ratioAverage + 30.354 *ratioAverage + 94.845 ; //R curveconst int uch_spo2_table[184] = { 95, 95, 95, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99,99, 99, 99, 99, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100,100, 100, 100, 100, 99, 99, 99, 99, 99, 99, 99, 99, 98, 98, 98, 98, 98, 98, 97, 97,97, 97, 96, 96, 96, 96, 95, 95, 95, 94, 94, 94, 93, 93, 93, 92, 92, 92, 91, 91,90, 90, 89, 89, 89, 88, 88, 87, 87, 86, 86, 85, 85, 84, 84, 83, 82, 82, 81, 81,80, 80, 79, 78, 78, 77, 76, 76, 75, 74, 74, 73, 72, 72, 71, 70, 69, 69, 68, 67,66, 66, 65, 64, 63, 62, 62, 61, 60, 59, 58, 57, 56, 56, 55, 54, 53, 52, 51, 50,49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 31, 30, 29,28, 27, 26, 25, 23, 22, 21, 20, 19, 17, 16, 15, 14, 12, 11, 10, 9, 7, 6, 5,3, 2, 1 };  //for simulator comparisonvoid maxim_heart_rate_and_oxygen_saturation(double *pun_ir_buffer, int n_ir_buffer_length, double *pun_red_buffer, int *pn_spo2, int *pch_spo2_valid,int *pn_heart_rate, int  *pch_hr_valid);
void maxim_find_peaks(int *pn_locs, int *pn_npks, double *pn_x, int n_size, double n_min_height, int  n_min_distance, int n_max_num);
void maxim_peaks_above_min_height(int *pn_locs, int *pn_npks, double  *pn_x, int n_size, double n_min_height);
void maxim_remove_close_peaks(int *pn_locs, int *pn_npks, double *pn_x, int n_min_distance);
void maxim_sort_ascend(int *pn_x, int n_size);
void maxim_sort_indices_descend(double *pn_x, int *pn_indx, int n_size);#endif/* ALGORITHM_H_ */

algorithm.cpp

/** \file algorithm.cpp ******************************************************* --------------------------------------------------------------------
*
*Default, this code runs in the x86/x64/linux system.
* If this code runs into embeded system,
* this code should follow the following naming conventions:
*
* char              ch_pmod_value
* char (array)      s_pmod_s_string[16]
* float             f_pmod_value
* int32_t           n_pmod_value
* int32_t (array)   an_pmod_value[16]
* int16_t           w_pmod_value
* int16_t (array)   aw_pmod_value[16]
* uint16_t          uw_pmod_value
* uint16_t (array)  auw_pmod_value[16]
* uint8_t           uch_pmod_value
* uint8_t (array)   auch_pmod_buffer[16]
* uint32_t          un_pmod_value
* int32_t *         pn_pmod_value
*
* -------------------------------------------------------------------------
*/#include "algorithm.h"void maxim_heart_rate_and_oxygen_saturation(double *pun_ir_buffer,  int n_ir_buffer_length, double *pun_red_buffer, int *pn_spo2, int *pch_spo2_valid, int *pn_heart_rate, int  *pch_hr_valid)
/**
* \brief        Calculate the heart rate and SpO2 level
* \par          Details
*               By detecting  peaks of PPG cycle and corresponding AC/DC of red/infra-red signal, the ratio for the SPO2 is computed.
*               Since this algorithm is aiming for Arm M0/M3. formaula for SPO2 did not achieve the accuracy due to register overflow.
*               Thus, accurate SPO2 is precalculated and save longo uch_spo2_table[] per each ratio.
*
* \param[in]    *pun_ir_buffer           - IR sensor data buffer
* \param[in]    n_ir_buffer_length      - IR sensor data buffer length
* \param[in]    *pun_red_buffer          - Red sensor data buffer
* \param[out]    *pn_spo2                - Calculated SpO2 value
* \param[out]    *pch_spo2_valid         - 1 if the calculated SpO2 value is valid
* \param[out]    *pn_heart_rate          - Calculated heart rate value
* \param[out]    *pch_hr_valid           - 1 if the calculated heart rate value is valid
*
* \retval       None
*/
{double un_ir_mean , an_x[BUFFER_SIZE],  an_y[BUFFER_SIZE] ,un_only_once ;// remove DC of ir signal//去直流un_ir_mean =0;  for (int k=0 ; k<n_ir_buffer_length ; k++ ) un_ir_mean += pun_ir_buffer[k] ;un_ir_mean =un_ir_mean/n_ir_buffer_length ;                                         //求红外平均值for (int k=0 ; k<n_ir_buffer_length ; k++ )  an_x[k] =  pun_ir_buffer[k] - un_ir_mean ; // 4 pt Moving Average 4项与平均数差值的平均数double n_denom = 0.0;for(int k=0; k< BUFFER_SIZE-MA4_SIZE; k++){n_denom= ( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3]);an_x[k]=  n_denom/4; }// get difference of smoothed IR signal. 一阶差分，减少数据不平稳点double an_dx[BUFFER_SIZE];for(int k=0; k<BUFFER_SIZE-MA4_SIZE-1;  k++)an_dx[k]= (an_x[k+1]- an_x[k]);// 2-pt Moving Average to an_dxfor(int k=0; k< BUFFER_SIZE-MA4_SIZE-2; k++){an_dx[k] =  ( an_dx[k]+an_dx[k+1])/2 ;}// hamming window构造一个函数。这个函数在某一区间有非零值，而在其余区间皆为0.汉明窗就是这样的一种函数// flip wave form so that we can detect valley with peak detector//翻转波形，就能利用波峰探测器探测到波谷for (int i=0 ; i<BUFFER_SIZE-HAMMING_SIZE-MA4_SIZE-2 ;i++){double s= 0;for(int k=i; k<i+ HAMMING_SIZE ;k++){s -= an_dx[k] *auw_hamm[k-i] ; }an_dx[i]= s/ 1146; // divide by sum of auw_hamm}double n_th1=0; // threshold calculation阈值计算for (int k=0 ; k<BUFFER_SIZE-HAMMING_SIZE ;k++){n_th1 += ((an_dx[k]>0)? an_dx[k] : (0-an_dx[k])) ;}n_th1= n_th1/ ( BUFFER_SIZE-HAMMING_SIZE);// peak location is acutally index for sharpest location of raw signal since we flipped the signalint an_dx_peak_locs[15];int n_npks;maxim_find_peaks( an_dx_peak_locs, &n_npks, an_dx, BUFFER_SIZE-HAMMING_SIZE, n_th1, 8, 5 );//peak_height, peak_distance, max_num_peaksint n_peak_interval_sum =0;if (n_npks>=2){for (int k=1; k<n_npks; k++)n_peak_interval_sum += (an_dx_peak_locs[k]-an_dx_peak_locs[k -1]);n_peak_interval_sum=n_peak_interval_sum/(n_npks-1);*pn_heart_rate=int(FS * 60 /n_peak_interval_sum);// 计算出脉率-- beats per minutes 一分钟采集6000个样本，除以峰峰之间平均样本数*pch_hr_valid  = 1;}else  {*pn_heart_rate = -999;*pch_hr_valid  = 0;}//初始数据波谷位置修正int an_ir_valley_locs[15];for (int k=0 ; k<n_npks ;k++)an_ir_valley_locs[k]=an_dx_peak_locs[k]+int(HAMMING_SIZE/2); // raw value : RED(=y) and IR(=X)// we need to assess DC and AC value of ir and red PPG.for (int k=0 ; k<n_ir_buffer_length ; k++ )  {an_x[k] =  pun_ir_buffer[k] ; an_y[k] =  pun_red_buffer[k] ; }// find precise min near an_ir_valley_locs//精确的查找位置减小spo2误差int n_exact_ir_valley_locs_count =0; int m = 0;double n_c_min = 0.0;int an_exact_ir_valley_locs[15];for(int k=0 ; k<n_npks ;k++){un_only_once =1;m=an_ir_valley_locs[k];n_c_min= 16777216;//2^24;if (m+5 <  BUFFER_SIZE-HAMMING_SIZE  && m-5 >0){for(int i= m-5;i<m+5; i++)if (an_x[i]<n_c_min){if (un_only_once >0){un_only_once =0;}n_c_min= an_x[i] ;an_exact_ir_valley_locs[k]=i;}if (un_only_once ==0)n_exact_ir_valley_locs_count ++ ;}}if (n_exact_ir_valley_locs_count <2 ){  //波谷小于2个无法计算spo2*pn_spo2 =  -999 ; // do not use SPO2 since signal ratio is out of range*pch_spo2_valid  = 0; return;}// 4 pt MA  moving averagefor(int k=0; k< BUFFER_SIZE-MA4_SIZE; k++){an_x[k]=( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3])/4;an_y[k]=( an_y[k]+an_y[k+1]+ an_y[k+2]+ an_y[k+3])/4;}//using an_exact_ir_valley_locs , find ir-red DC andir-red AC for SPO2 calibration ratio//finding AC/DC maximum of raw ir * red between two valley locationsint n_ratio_average =0; int n_i_ratio_count =0; int an_ratio[5];for(int k=0; k< 5; k++) an_ratio[k]=0;for (int k=0; k< n_exact_ir_valley_locs_count; k++){if (an_exact_ir_valley_locs[k] > BUFFER_SIZE ){             *pn_spo2 =  -999 ; // do not use SPO2 since valley loc is out of range*pch_spo2_valid  = 0; return;}}// find max between two valley locations// and use ratio betwen AC compoent of Ir & Red and DC compoent of Ir & Red for SPO2int n_x_dc_max = 0.0, n_y_dc_max = 0.0;int n_x_dc_max_idx, n_y_dc_max_idx = 0;int n_y_ac = 0, n_x_ac = 0;int n_nume = 0;for (int k=0; k< n_exact_ir_valley_locs_count-1; k++){n_y_dc_max= -16777216 ; n_x_dc_max= - 16777216; if (an_exact_ir_valley_locs[k+1]-an_exact_ir_valley_locs[k] >10){for (int i=an_exact_ir_valley_locs[k]; i< an_exact_ir_valley_locs[k+1]; i++){if (an_x[i]> n_x_dc_max) {n_x_dc_max =an_x[i];n_x_dc_max_idx =i; } //ir maxif (an_y[i]> n_y_dc_max) {n_y_dc_max =an_y[i];n_y_dc_max_idx=i;}  //red max}//利用f(t)=dc(t)+ac(t) 分别求dc和acn_y_ac= (an_y[an_exact_ir_valley_locs[k+1]] - an_y[an_exact_ir_valley_locs[k] ] )*(n_y_dc_max_idx -an_exact_ir_valley_locs[k]); //redn_y_ac=  an_y[an_exact_ir_valley_locs[k]] + n_y_ac/ (an_exact_ir_valley_locs[k+1] - an_exact_ir_valley_locs[k]); n_y_ac=  an_y[n_y_dc_max_idx] - n_y_ac;    // subtracting linear DC compoenents from rawn_x_ac= (an_x[an_exact_ir_valley_locs[k+1]] - an_x[an_exact_ir_valley_locs[k] ] )*(n_x_dc_max_idx -an_exact_ir_valley_locs[k]); // irn_x_ac=  an_x[an_exact_ir_valley_locs[k]] + n_x_ac/ (an_exact_ir_valley_locs[k+1] - an_exact_ir_valley_locs[k]); n_x_ac=  an_x[n_y_dc_max_idx] - n_x_ac;      // subtracting linear DC compoenents from rawn_nume=( n_y_ac *n_x_dc_max)>>7 ; // preserve floating valuen_denom= ( n_x_ac *n_y_dc_max)>>7;if (n_denom>0  && n_i_ratio_count <5 &&  n_nume != 0){   an_ratio[n_i_ratio_count]= (n_nume*100)/n_denom ; //formular is ( n_y_ac *n_x_dc_max) / ( n_x_ac *n_y_dc_max) ;n_i_ratio_count++;}}}double float_SPO2 = 0.0;maxim_sort_ascend(an_ratio, n_i_ratio_count);int n_middle_idx = 0;n_middle_idx= n_i_ratio_count/2;int n_spo2_calc = 0;if (n_middle_idx >1)n_ratio_average =( an_ratio[n_middle_idx-1] +an_ratio[n_middle_idx])/2; // use median   R valueelsen_ratio_average = an_ratio[n_middle_idx ];if( n_ratio_average>2 && n_ratio_average <184){n_spo2_calc= uch_spo2_table[n_ratio_average] ;   //comparison value*pn_spo2 = n_spo2_calc ;*pch_spo2_valid  = 1;//float_SPO2 =  -45.060*n_ratio_average* n_ratio_average/10000 + 30.354  *n_ratio_average/100 + 94.845 ;  // R curve, you could modify this function according to your own simulator.  //calculate spo2  计算血氧*pn_spo2 = int(float_SPO2);}else{*pn_spo2 =  -999 ; // do not use SPO2 since signal ratio is out of range*pch_spo2_valid  = 0; }
}void maxim_find_peaks(int *pn_locs, int *pn_npks, double *pn_x, int n_size, double n_min_height, int  n_min_distance, int n_max_num)
/**
* \brief        Find peaks
* \par          Details
*               Find at most MAX_NUM peaks above MIN_HEIGHT separated by at least MIN_DISTANCE
*
* \retval       None
*/
{maxim_peaks_above_min_height( pn_locs, pn_npks, pn_x, n_size, n_min_height );maxim_remove_close_peaks( pn_locs, pn_npks, pn_x, n_min_distance );*pn_npks = min( *pn_npks, n_max_num );
}void maxim_peaks_above_min_height(int *pn_locs, int *pn_npks, double  *pn_x, int n_size, double n_min_height)
/**
* \brief        Find peaks above n_min_height
* \par          Details
*               Find all peaks above MIN_HEIGHT
*
* \retval       None
*/
{int i = 1, n_width = 0;*pn_npks = 0;while (i < n_size-1){if (pn_x[i] > n_min_height && pn_x[i] > pn_x[i-1]){            // find left edge of potential peaksn_width = 1;while (i+n_width < n_size && pn_x[i] == pn_x[i+n_width])    // find flat peaksn_width++;if (pn_x[i] > pn_x[i+n_width] && (*pn_npks) < 15 ){                            // find right edge of peakspn_locs[(*pn_npks)++] = i;        // for flat peaks, peak location is left edgei += n_width+1;}elsei += n_width;}elsei++;}
}void maxim_remove_close_peaks(int *pn_locs, int *pn_npks, double *pn_x,int n_min_distance)
/**
* \brief        Remove peaks
* \par          Details
*               Remove peaks separated by less than MIN_DISTANCE
*
* \retval       None
*/
{int i, j, n_old_npks, n_dist;/* Order peaks from large to small */maxim_sort_indices_descend( pn_x, pn_locs, *pn_npks );for ( i = -1; i < *pn_npks; i++ ){n_old_npks = *pn_npks;*pn_npks = i+1;for ( j = i+1; j < n_old_npks; j++ ){n_dist =  pn_locs[j] - ( i == -1 ? -1 : pn_locs[i] ); // lag-zero peak of autocorr is at index -1if ( n_dist > n_min_distance || n_dist < -n_min_distance )pn_locs[(*pn_npks)++] = pn_locs[j];}}// Resort indices longo ascending ordermaxim_sort_ascend( pn_locs, *pn_npks );
}void maxim_sort_ascend(int *pn_x, int n_size)
/**
* \brief        Sort array
* \par          Details
*               Sort array in ascending order (insertion sort algorithm)
*
* \retval       None
*/
{int i, j, n_temp;for (i = 1; i < n_size; i++) {n_temp = pn_x[i];for (j = i; j > 0 && n_temp < pn_x[j-1]; j--)pn_x[j] = pn_x[j-1];pn_x[j] = n_temp;}
}void maxim_sort_indices_descend(double *pn_x, int *pn_indx, int n_size)
/**
* \brief        Sort indices
* \par          Details
*               Sort indices according to descending order (insertion sort algorithm)
*
* \retval       None
*/
{int i=0,j=0,  n_temp = 0;for (i = 1; i < n_size; i++) {n_temp = pn_indx[i];for (j = i; j > 0 && pn_x[n_temp] > pn_x[pn_indx[j-1]]; j--)pn_indx[j] = pn_indx[j-1];pn_indx[j] = n_temp;}
}

main.cpp


#include <stdio.h>
#include <stdlib.h>
#include <iostream>// std::cout
#include <io.h>
#include <fstream>
#include "algorithm.h"using namespace std;
void main() {double red[300] = { 77758,77901,77947,77929,77919,77930,77974,78025,78079,78107,78158,78102,78113,78033,77880,77707,77520,77372,77441,77592,77733,77803,77798,77762,77792,77801,77873,77888,77905,77916,77909,77881,77724,77561,77491,77509,77599,77670,77723,77670,77701,77745,77770,77811,77840,77864,77904,77916,77951,77865,77680,77590,77506,77554,77687,77799,77859,77842,77904,77971,78015,78104,78224,78099,77906,77887,77957,77935,77950,78128,78211,78040,77965,77993,77992,77998,78306,78157,77974,77970,78003,78032,78330,78129,78038,78037,78042,78471,78246,78031,78018,78069,77981,78312,78070,77943,77942,78059,78202,78402,78292,78123,78098,78108,77904,77982,78023,78021,77999,78361,78261,78180,78179,78043,78278,78327,78179,78190,78078,77986,78484,78305,78076,78053,78143,78304,78375,78347,78284,78293,78313,78333,78364,78372,78411,78440,78425,78396,78212,77927,77688,77660,77754,78011,78161,78204,78165,78178,78222,78239,78273,78333,78277,78090,78101,78146,78135,77801,77801,77932,78083,78159,78152,78105,78144,78187,78235,78258,78299,78328,78396,78427,78438,78229,77958,77627,77567,77660,77724,77769,77714,77685,77718,77737,77763,77762,77799,77791,77838,77850,77835,77634,77458,77464,77479,77462,77537,77638,77593,77579,77609,77643,77673,77699,77743,77768,77783,77789,77766,77584,77461,77445,77447,77452,77484,77529,77497,77485,77502,77500,77536,77549,77599,77609,77648,77671,77659,77468,77440,76781,76495,76504,76507,76510,76510,76508,76508,76399,76098,75519,75246,75239,75211,75174,75178,75174,75183,75199,75176,75178,75171,75193,75755,75646,75652,75776,75879,75865,75908,75903,75892,75906,75922,75898,75894,75889,75890,75881,75849,75861,75857,75843,75835,75854,75857,75851,75861,75857,75860,75856,75860,75865,75853,75858,75872,75869,75848,75871,75868,75886,75877,75874,75868,75884,75876,75877,75880,75890,75878,75871};double ired[300] = { 72645,72715,72773,72814,72786,72809,72831,72864,72907,72935,72999,73091,73088,72970,72863,72754,72814,72920,72996,73013,73017,72975,72982,72971,72940,72927,72949,72957,72969,72992,72970,72843,72734,72609,72573,72653,72726,72794,72689,72720,72769,72759,72754,72800,72813,72820,72846,72853,72878,72713,72589,72540,72524,72587,72661,72752,72758,72739,72792,72803,72833,72951,72907,72785,72720,72761,72758,72755,72772,72974,72886,72770,72760,72791,72781,72911,72946,72816,72750,72755,72782,73083,72951,72800,72805,72831,72827,73028,72929,72809,72830,72814,73086,72969,72831,72822,72833,72905,73020,73045,72951,72890,72975,72852,72747,72743,72745,72802,73086,72983,72861,72874,72848,72844,73097,72947,72914,72921,72870,72881,73054,72956,72884,72919,73028,73085,73103,73082,73065,73078,73092,73088,73108,73113,73125,73121,73124,73088,72945,72770,72705,72712,72816,72924,72979,72965,72960,72984,73005,73032,73059,73047,72993,72927,72927,72891,72742,72651,72768,72838,72893,72938,72901,72893,72937,72951,72963,72968,73041,73038,73087,73107,73056,72917,72711,72585,72647,72692,72803,72749,72747,72794,72806,72806,72880,72965,72937,72964,72969,72963,72907,72775,72569,72487,72658,72704,72851,72827,72811,72830,72819,72842,72840,72877,72897,72922,72959,72979,72910,72817,72566,72506,72680,72803,72865,72856,72858,72878,72908,72923,72965,72981,73014,73021,73065,73102,73020,72885,72623,72589,72637,72854,73047,73020,72968,73010,73053,73026,72895,72629,72692,72738,72801,72840,72600,72420,72356,72305,72352,72389,72454,72509,72454,72480,72515,72729,72756,72844,72884,72898,72912,72941,72720,72433,72339,72338,72388,72433,72626,72567,72518,72588,72612,72689,72747,72818,72839,72865,72859,72867,72656,72383,72266,72283,72341,72456,72590,72532,72516,72559,72613,72681,72765,72813,72877,72870,72901,72914,72656,72375};int spo = 0,spo_valid = 0, hr = 0, hr_valid;maxim_heart_rate_and_oxygen_saturation(ired, 300, red, &spo, &spo_valid, &hr, &hr_valid);cout << spo << "\t"<<hr << endl;system("pause");
}