原网址(侵删):使用PPG信号计算脉率和血氧_2-pt moving average to an_dx-CSDN博客
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 curve
const 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 comparison
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);
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_dx
for(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 signal
int 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_peaks
int 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 average
for(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 locations
int 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 SPO2
int 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 max
if (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和ac
n_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]); //red
n_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 raw
n_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]); // ir
n_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 raw
n_nume=( n_y_ac *n_x_dc_max)>>7 ; // preserve floating value
n_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 value
else
n_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 peaks
n_width = 1;
while (i+n_width < n_size && pn_x[i] == pn_x[i+n_width]) // find flat peaks
n_width++;
if (pn_x[i] > pn_x[i+n_width] && (*pn_npks) < 15 ){ // find right edge of peaks
pn_locs[(*pn_npks)++] = i;
// for flat peaks, peak location is left edge
i += n_width+1;
}
else
i += n_width;
}
else
i++;
}
}
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 -1
if ( n_dist > n_min_distance || n_dist < -n_min_distance )
pn_locs[(*pn_npks)++] = pn_locs[j];
}
}
// Resort indices longo ascending order
maxim_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");
}