Alexandru Radu
gg
Fork of
RD117_MBED
by 
Maxim Integrated
Diff: algorithm/algorithm.cpp
- Revision:
- 4:5273ab1085ab
- Parent:
- 3:7c0fb55eb3ff
diff -r 7c0fb55eb3ff -r 5273ab1085ab algorithm/algorithm.cpp
--- a/algorithm/algorithm.cpp Thu Apr 21 19:38:17 2016 +0000
+++ b/algorithm/algorithm.cpp Thu May 19 22:04:06 2016 +0000
@@ -4,31 +4,28 @@
* Filename: algorithm.cpp
* Description: This module calculates the heart rate/SpO2 level
*
-* Revision History:
-*\n 1-18-2016 Rev 01.00 SK Initial release.
-*\n
*
* --------------------------------------------------------------------
*
* This code follows 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
+* 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
*
* ------------------------------------------------------------------------- */
/*******************************************************************************
-* Copyright (C) 2015 Maxim Integrated Products, Inc., All Rights Reserved.
+* Copyright (C) 2016 Maxim Integrated Products, Inc., All Rights Reserved.
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
@@ -62,8 +59,8 @@
#include "algorithm.h"
#include "mbed.h"
-void maxim_heart_rate_and_oxygen_saturation(uint32_t *pun_ir_buffer , int32_t n_ir_buffer_length, uint32_t *pun_red_buffer , int32_t *pn_spo2, int8_t *pch_spo2_valid ,
- int32_t *pn_heart_rate , int8_t *pch_hr_valid)
+void maxim_heart_rate_and_oxygen_saturation(uint32_t *pun_ir_buffer, int32_t n_ir_buffer_length, uint32_t *pun_red_buffer, int32_t *pn_spo2, int8_t *pch_spo2_valid,
+ int32_t *pn_heart_rate, int8_t *pch_hr_valid)
/**
* \brief Calculate the heart rate and SpO2 level
* \par Details
@@ -82,191 +79,179 @@
* \retval None
*/
{
-
+ uint32_t un_ir_mean ,un_only_once ;
+ int32_t k ,n_i_ratio_count;
+ int32_t i,s ,m, n_exact_ir_valley_locs_count ,n_middle_idx;
+ int32_t n_th1, n_npks,n_c_min;
+ int32_t an_ir_valley_locs[15] ;
+ int32_t an_exact_ir_valley_locs[15] ;
+ int32_t an_dx_peak_locs[15] ;
+ int32_t n_peak_interval_sum;
+
+ int32_t n_y_ac, n_x_ac;
+ int32_t n_spo2_calc;
+ int32_t n_y_dc_max, n_x_dc_max;
+ int32_t n_y_dc_max_idx, n_x_dc_max_idx;
+ int32_t an_ratio[5],n_ratio_average;
+ int32_t n_nume, n_denom ;
+ // remove DC of ir signal
+ un_ir_mean =0;
+ for (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 (k=0 ; k<n_ir_buffer_length ; k++ ) an_x[k] = pun_ir_buffer[k] - un_ir_mean ;
+
+ // 4 pt Moving Average
+ for(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/(int32_t)4;
+ }
- uint32_t un_ir_mean ,un_only_once ;
- int32_t k ,n_i_ratio_count;
- int32_t i,s ,m, n_exact_ir_valley_locs_count ,n_middle_idx;
- int32_t n_th1, n_npks,n_c_min;
- int32_t an_ir_valley_locs[15] ;
- int32_t an_exact_ir_valley_locs[15] ;
- int32_t an_dx_peak_locs[15] ;
- int32_t n_peak_interval_sum;
-
- int32_t n_y_ac, n_x_ac;
- int32_t n_spo2_calc;
- int32_t n_y_dc_max, n_x_dc_max;
- int32_t n_y_dc_max_idx, n_x_dc_max_idx;
- int32_t an_ratio[5],n_ratio_average;
- int32_t n_nume, n_denom ;
- // remove DC of ir signal
- un_ir_mean =0;
- for (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 (k=0 ; k<n_ir_buffer_length ; k++ ) an_x[k] = pun_ir_buffer[k] - un_ir_mean ;
-
- // 4 pt Moving Average
- for(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/(int32_t)4;
- }
-
- // get difference of smoothed IR signal
-
- for( k=0; k<BUFFER_SIZE-MA4_SIZE-1; k++)
- an_dx[k]= (an_x[k+1]- an_x[k]);
+ // get difference of smoothed IR signal
+
+ for( 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(k=0; k< BUFFER_SIZE-MA4_SIZE-2; k++){
- an_dx[k] = ( an_dx[k]+an_dx[k+1])/2 ;
- }
-
- // hamming window
- // flip wave form so that we can detect valley with peak detector
- for ( i=0 ; i<BUFFER_SIZE-HAMMING_SIZE-MA4_SIZE-2 ;i++){
- s= 0;
- for( k=i; k<i+ HAMMING_SIZE ;k++){
- s -= an_dx[k] *auw_hamm[k-i] ;
- }
- an_dx[i]= s/ (int32_t)1146; // divide by sum of auw_hamm
- }
+ // 2-pt Moving Average to an_dx
+ for(k=0; k< BUFFER_SIZE-MA4_SIZE-2; k++){
+ an_dx[k] = ( an_dx[k]+an_dx[k+1])/2 ;
+ }
-
- n_th1=0; // threshold calculation
- for ( k=0 ; k<BUFFER_SIZE-HAMMING_SIZE ;k++){
- n_th1 += ((an_dx[k]>0)? an_dx[k] : ((int32_t)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
- 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
+ // hamming window
+ // flip wave form so that we can detect valley with peak detector
+ for ( i=0 ; i<BUFFER_SIZE-HAMMING_SIZE-MA4_SIZE-2 ;i++){
+ s= 0;
+ for( k=i; k<i+ HAMMING_SIZE ;k++){
+ s -= an_dx[k] *auw_hamm[k-i] ;
+ }
+ an_dx[i]= s/ (int32_t)1146; // divide by sum of auw_hamm
+ }
- n_peak_interval_sum =0;
- if (n_npks>=2){
- for (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 =(int32_t)( 6000/ n_peak_interval_sum );// beats per minutes
- //prlongf(">>> *pn_heart_rate= %d \n", *pn_heart_rate) ;
- *pch_hr_valid = 1;
- }
- else {
- *pn_heart_rate = -999;
- *pch_hr_valid = 0;
- }
-
- for ( k=0 ; k<n_npks ;k++)
- an_ir_valley_locs[k]= an_dx_peak_locs[k] +HAMMING_SIZE /2;
+
+ n_th1=0; // threshold calculation
+ for ( k=0 ; k<BUFFER_SIZE-HAMMING_SIZE ;k++){
+ n_th1 += ((an_dx[k]>0)? an_dx[k] : ((int32_t)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
+ 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
+
+ n_peak_interval_sum =0;
+ if (n_npks>=2){
+ for (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=(int32_t)(6000/n_peak_interval_sum);// beats per minutes
+ *pch_hr_valid = 1;
+ }
+ else {
+ *pn_heart_rate = -999;
+ *pch_hr_valid = 0;
+ }
+
+ for ( k=0 ; k<n_npks ;k++)
+ an_ir_valley_locs[k]=an_dx_peak_locs[k]+HAMMING_SIZE/2;
- // raw value : RED(=y) and IR(=X)
- // we need to assess DC and AC value of ir and red PPG.
- for (k=0 ; k<n_ir_buffer_length ; k++ ) {
- an_x[k] = pun_ir_buffer[k] ;
- an_y[k] = pun_red_buffer[k] ;
- }
+ // raw value : RED(=y) and IR(=X)
+ // we need to assess DC and AC value of ir and red PPG.
+ for (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
+ n_exact_ir_valley_locs_count =0;
+ for(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(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 ){
+ *pn_spo2 = -999 ; // do not use SPO2 since signal ratio is out of range
+ *pch_spo2_valid = 0;
+ return;
+ }
+ // 4 pt MA
+ for(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])/(int32_t)4;
+ an_y[k]=( an_y[k]+an_y[k+1]+ an_y[k+2]+ an_y[k+3])/(int32_t)4;
+ }
- // find precise min near an_ir_valley_locs
- n_exact_ir_valley_locs_count =0;
- for ( 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(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 ++ ;
+ //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
+ n_ratio_average =0;
+ n_i_ratio_count =0;
+
+ for(k=0; k< 5; k++) an_ratio[k]=0;
+ for (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
+
+ for (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 (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; }
+ if (an_y[i]> n_y_dc_max) {n_y_dc_max =an_y[i];n_y_dc_max_idx=i;}
+ }
+ 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; // subracting 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; // subracting linear DC compoenents from raw
+ n_nume=( n_y_ac *n_x_dc_max)>>7 ; //prepare X100 to 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++;
}
}
- if (n_exact_ir_valley_locs_count <2 ){
- *pn_spo2 = -999 ; // do not use SPO2 since signal ratio is out of range
- *pch_spo2_valid = 0;
- return;
- }
- // 4 pt MA
- for(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])/(int32_t)4;
- an_y[k]=( an_y[k]+an_y[k+1]+ an_y[k+2]+ an_y[k+3])/(int32_t)4;
-
- }
+ }
+
+ maxim_sort_ascend(an_ratio, n_i_ratio_count);
+ n_middle_idx= n_i_ratio_count/2;
+
+ if (n_middle_idx >1)
+ n_ratio_average =( an_ratio[n_middle_idx-1] +an_ratio[n_middle_idx])/2; // use median
+ else
+ n_ratio_average = an_ratio[n_middle_idx ];
- //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
- n_ratio_average =0;
- n_i_ratio_count = 0;
-
- for(k=0; k< 5; k++) an_ratio[k]=0;
- for (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
-
- for (k=0; k< n_exact_ir_valley_locs_count-1; k++){
- n_y_dc_max= -16777216 ;
- n_x_dc_max= - 16777216;
- // printf("range=%d: %d\n ", an_exact_ir_valley_locs[k], an_exact_ir_valley_locs[k+1]);
- if (an_exact_ir_valley_locs[k+1]-an_exact_ir_valley_locs[k] >10){
- for (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; }
- if (an_y[i]> n_y_dc_max) {n_y_dc_max =an_y[i];n_y_dc_max_idx=i;}
- }
- 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; // subracting 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; // subracting linear DC compoenents from raw
- n_nume=( n_y_ac *n_x_dc_max)>>7 ; //prepare X100 to 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++;
- }
- }
-
- // prlongf("an_ratio[%d]= %d an_exact_ir_valley_locs[k] =%d , an_exact_ir_valley_locs[%d] =%d \n",k, an_ratio[k] ,an_exact_ir_valley_locs[k] ,k+1, an_exact_ir_valley_locs[k+1] ) ;
- // prlongf("n_nume= %d ,n_denom= %d n_y_ac = %d, n_x_dc_max = %d, n_x_ac= %d, n_y_dc_max = %d\n",n_nume, n_denom, n_y_ac ,n_x_dc_max ,n_x_ac ,n_y_dc_max );
-
- }
-
- maxim_sort_ascend(an_ratio, n_i_ratio_count);
- n_middle_idx= n_i_ratio_count/2;
-
- if (n_middle_idx >1)
- n_ratio_average =( an_ratio[n_middle_idx-1] +an_ratio[n_middle_idx])/2; // use median
- 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] ;
- *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 ; // for comparison with table
- }
- else{
- *pn_spo2 = -999 ; // do not use SPO2 since signal ratio is out of range
- *pch_spo2_valid = 0;
- }
-
-
+ if( n_ratio_average>2 && n_ratio_average <184){
+ n_spo2_calc= uch_spo2_table[n_ratio_average] ;
+ *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 ; // for comparison with table
+ }
+ else{
+ *pn_spo2 = -999 ; // do not use SPO2 since signal ratio is out of range
+ *pch_spo2_valid = 0;
+ }
}
-void maxim_find_peaks( int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x, int32_t n_size, int32_t n_min_height, int32_t n_min_distance, int32_t n_max_num )
+void maxim_find_peaks(int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x, int32_t n_size, int32_t n_min_height, int32_t n_min_distance, int32_t n_max_num)
/**
* \brief Find peaks
* \par Details
@@ -280,7 +265,7 @@
*pn_npks = min( *pn_npks, n_max_num );
}
-void maxim_peaks_above_min_height( int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x, int32_t n_size, int32_t n_min_height )
+void maxim_peaks_above_min_height(int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x, int32_t n_size, int32_t n_min_height)
/**
* \brief Find peaks above n_min_height
* \par Details
@@ -299,7 +284,7 @@
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
+ // for flat peaks, peak location is left edge
i += n_width+1;
}
else
@@ -311,7 +296,7 @@
}
-void maxim_remove_close_peaks( int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x, int32_t n_min_distance )
+void maxim_remove_close_peaks(int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x,int32_t n_min_distance)
/**
* \brief Remove peaks
* \par Details
@@ -340,7 +325,7 @@
maxim_sort_ascend( pn_locs, *pn_npks );
}
-void maxim_sort_ascend(int32_t *pn_x, int32_t n_size)
+void maxim_sort_ascend(int32_t *pn_x,int32_t n_size)
/**
* \brief Sort array
* \par Details
@@ -358,7 +343,7 @@
}
}
-void maxim_sort_indices_descend( int32_t *pn_x, int32_t *pn_indx, int32_t n_size)
+void maxim_sort_indices_descend(int32_t *pn_x, int32_t *pn_indx, int32_t n_size)
/**
* \brief Sort indices
* \par Details