spo2_algorithm.cpp

00001 #include "spo2_algorithm.h"
00002 
00003 
00004 
00005 const uint16_t auw_hamm[31]={41,276,512,276,41 }; //Hamm=  long16(512* hamming(5)');
00006 //SPO2table is computed as  -45.060*ratioAverage* ratioAverage + 30.354 *ratioAverage + 94.845 ;
00007 const uint8_t uch_spo2_table[184]={ 95, 95, 95, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 
00008                             99, 99, 99, 99, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
00009                             100, 100, 100, 100, 99, 99, 99, 99, 99, 99, 99, 99, 98, 98, 98, 98, 98, 98, 97, 97, 
00010                             97, 97, 96, 96, 96, 96, 95, 95, 95, 94, 94, 94, 93, 93, 93, 92, 92, 92, 91, 91, 
00011                             90, 90, 89, 89, 89, 88, 88, 87, 87, 86, 86, 85, 85, 84, 84, 83, 82, 82, 81, 81, 
00012                             80, 80, 79, 78, 78, 77, 76, 76, 75, 74, 74, 73, 72, 72, 71, 70, 69, 69, 68, 67, 
00013                             66, 66, 65, 64, 63, 62, 62, 61, 60, 59, 58, 57, 56, 56, 55, 54, 53, 52, 51, 50, 
00014                             49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 31, 30, 29, 
00015                             28, 27, 26, 25, 23, 22, 21, 20, 19, 17, 16, 15, 14, 12, 11, 10, 9, 7, 6, 5, 
00016                             3, 2, 1 } ;
00017 static  int32_t an_dx[ BUFFER_SIZE-MA4_SIZE]; // delta
00018 static  int32_t an_x[ BUFFER_SIZE]; //ir
00019 static  int32_t an_y[ BUFFER_SIZE]; //red
00020 
00021 void spo2_algorithm::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, 
00022                               int32_t *pn_heart_rate, int8_t  *pch_hr_valid)
00023 /**
00024 * \brief        Calculate the heart rate and SpO2 level
00025 * \par          Details
00026 *               By detecting  peaks of PPG cycle and corresponding AC/DC of red/infra-red signal, the ratio for the SPO2 is computed.
00027 *               Since this algorithm is aiming for Arm M0/M3. formaula for SPO2 did not achieve the accuracy due to register overflow.
00028 *               Thus, accurate SPO2 is precalculated and save longo uch_spo2_table[] per each ratio.
00029 *
00030 * \param[in]    *pun_ir_buffer           - IR sensor data buffer
00031 * \param[in]    n_ir_buffer_length      - IR sensor data buffer length
00032 * \param[in]    *pun_red_buffer          - Red sensor data buffer
00033 * \param[out]    *pn_spo2                - Calculated SpO2 value
00034 * \param[out]    *pch_spo2_valid         - 1 if the calculated SpO2 value is valid
00035 * \param[out]    *pn_heart_rate          - Calculated heart rate value
00036 * \param[out]    *pch_hr_valid           - 1 if the calculated heart rate value is valid
00037 *
00038 * \retval       None
00039 */
00040 {
00041     uint32_t un_ir_mean ,un_only_once ;
00042     int32_t k ,n_i_ratio_count;
00043     int32_t i,s ,m, n_exact_ir_valley_locs_count ,n_middle_idx;
00044     int32_t n_th1, n_npks,n_c_min;      
00045     int32_t an_ir_valley_locs[15] ;
00046     int32_t an_exact_ir_valley_locs[15] ;
00047     int32_t an_dx_peak_locs[15] ;
00048     int32_t n_peak_interval_sum;
00049     
00050     int32_t n_y_ac, n_x_ac;
00051     int32_t n_spo2_calc; 
00052     int32_t n_y_dc_max, n_x_dc_max; 
00053     int32_t n_y_dc_max_idx, n_x_dc_max_idx; 
00054     int32_t an_ratio[5],n_ratio_average; 
00055     int32_t n_nume,  n_denom ;
00056     // remove DC of ir signal    
00057     un_ir_mean =0; 
00058     for (k=0 ; k<n_ir_buffer_length ; k++ ) un_ir_mean += pun_ir_buffer[k] ;
00059     un_ir_mean =un_ir_mean/n_ir_buffer_length ;
00060     for (k=0 ; k<n_ir_buffer_length ; k++ )  an_x[k] =  pun_ir_buffer[k] - un_ir_mean ; 
00061     
00062     // 4 pt Moving Average
00063     for(k=0; k< BUFFER_SIZE-MA4_SIZE; k++){
00064         n_denom= ( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3]);
00065         an_x[k]=  n_denom/(int32_t)4; 
00066     }
00067 
00068     // get difference of smoothed IR signal
00069     
00070     for( k=0; k<BUFFER_SIZE-MA4_SIZE-1;  k++)
00071         an_dx[k]= (an_x[k+1]- an_x[k]);
00072 
00073     // 2-pt Moving Average to an_dx
00074     for(k=0; k< BUFFER_SIZE-MA4_SIZE-2; k++){
00075         an_dx[k] =  ( an_dx[k]+an_dx[k+1])/2 ;
00076     }
00077     
00078     // hamming window
00079     // flip wave form so that we can detect valley with peak detector
00080     for ( i=0 ; i<BUFFER_SIZE-HAMMING_SIZE-MA4_SIZE-2 ;i++){
00081         s= 0;
00082         for( k=i; k<i+ HAMMING_SIZE ;k++){
00083             s -= an_dx[k] *auw_hamm[k-i] ; 
00084                      }
00085         an_dx[i]= s/ (int32_t)1146; // divide by sum of auw_hamm 
00086     }
00087 
00088  
00089     n_th1=0; // threshold calculation
00090     for ( k=0 ; k<BUFFER_SIZE-HAMMING_SIZE ;k++){
00091         n_th1 += ((an_dx[k]>0)? an_dx[k] : ((int32_t)0-an_dx[k])) ;
00092     }
00093     n_th1= n_th1/ ( BUFFER_SIZE-HAMMING_SIZE);
00094     // peak location is acutally index for sharpest location of raw signal since we flipped the signal         
00095     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 
00096 
00097     n_peak_interval_sum =0;
00098     if (n_npks>=2){
00099         for (k=1; k<n_npks; k++)
00100             n_peak_interval_sum += (an_dx_peak_locs[k]-an_dx_peak_locs[k -1]);
00101         n_peak_interval_sum=n_peak_interval_sum/(n_npks-1);
00102         *pn_heart_rate=(int32_t)(6000/n_peak_interval_sum);// beats per minutes
00103         *pch_hr_valid  = 1;
00104     }
00105     else  {
00106         *pn_heart_rate = -999;
00107         *pch_hr_valid  = 0;
00108     }
00109             
00110     for ( k=0 ; k<n_npks ;k++)
00111         an_ir_valley_locs[k]=an_dx_peak_locs[k]+HAMMING_SIZE/2; 
00112 
00113 
00114     // raw value : RED(=y) and IR(=X)
00115     // we need to assess DC and AC value of ir and red PPG. 
00116     for (k=0 ; k<n_ir_buffer_length ; k++ )  {
00117         an_x[k] =  pun_ir_buffer[k] ; 
00118         an_y[k] =  pun_red_buffer[k] ; 
00119     }
00120 
00121     // find precise min near an_ir_valley_locs
00122     n_exact_ir_valley_locs_count =0; 
00123     for(k=0 ; k<n_npks ;k++){
00124         un_only_once =1;
00125         m=an_ir_valley_locs[k];
00126         n_c_min= 16777216;//2^24;
00127         if (m+5 <  BUFFER_SIZE-HAMMING_SIZE  && m-5 >0){
00128             for(i= m-5;i<m+5; i++)
00129                 if (an_x[i]<n_c_min){
00130                     if (un_only_once >0){
00131                        un_only_once =0;
00132                    } 
00133                    n_c_min= an_x[i] ;
00134                    an_exact_ir_valley_locs[k]=i;
00135                 }
00136             if (un_only_once ==0)
00137                 n_exact_ir_valley_locs_count ++ ;
00138         }
00139     }
00140     if (n_exact_ir_valley_locs_count <2 ){
00141        *pn_spo2 =  -999 ; // do not use SPO2 since signal ratio is out of range
00142        *pch_spo2_valid  = 0; 
00143        return;
00144     }
00145     // 4 pt MA
00146     for(k=0; k< BUFFER_SIZE-MA4_SIZE; k++){
00147         an_x[k]=( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3])/(int32_t)4;
00148         an_y[k]=( an_y[k]+an_y[k+1]+ an_y[k+2]+ an_y[k+3])/(int32_t)4;
00149     }
00150 
00151     //using an_exact_ir_valley_locs , find ir-red DC andir-red AC for SPO2 calibration ratio
00152     //finding AC/DC maximum of raw ir * red between two valley locations
00153     n_ratio_average =0; 
00154     n_i_ratio_count =0; 
00155     
00156     for(k=0; k< 5; k++) an_ratio[k]=0;
00157     for (k=0; k< n_exact_ir_valley_locs_count; k++){
00158         if (an_exact_ir_valley_locs[k] > BUFFER_SIZE ){             
00159             *pn_spo2 =  -999 ; // do not use SPO2 since valley loc is out of range
00160             *pch_spo2_valid  = 0; 
00161             return;
00162         }
00163     }
00164     // find max between two valley locations 
00165     // and use ratio betwen AC compoent of Ir & Red and DC compoent of Ir & Red for SPO2 
00166 
00167     for (k=0; k< n_exact_ir_valley_locs_count-1; k++){
00168         n_y_dc_max= -16777216 ; 
00169         n_x_dc_max= - 16777216; 
00170         if (an_exact_ir_valley_locs[k+1]-an_exact_ir_valley_locs[k] >10){
00171             for (i=an_exact_ir_valley_locs[k]; i< an_exact_ir_valley_locs[k+1]; i++){
00172                 if (an_x[i]> n_x_dc_max) {n_x_dc_max =an_x[i];n_x_dc_max_idx =i; }
00173                 if (an_y[i]> n_y_dc_max) {n_y_dc_max =an_y[i];n_y_dc_max_idx=i;}
00174             }
00175             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
00176             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])  ; 
00177         
00178         
00179             n_y_ac=  an_y[n_y_dc_max_idx] - n_y_ac;    // subracting linear DC compoenents from raw 
00180             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
00181             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]); 
00182             n_x_ac=  an_x[n_y_dc_max_idx] - n_x_ac;      // subracting linear DC compoenents from raw 
00183             n_nume=( n_y_ac *n_x_dc_max)>>7 ; //prepare X100 to preserve floating value
00184             n_denom= ( n_x_ac *n_y_dc_max)>>7;
00185             if (n_denom>0  && n_i_ratio_count <5 &&  n_nume != 0)
00186             {   
00187                 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) ;
00188                 n_i_ratio_count++;
00189             }
00190         }
00191     }
00192 
00193     maxim_sort_ascend(an_ratio, n_i_ratio_count);
00194     n_middle_idx= n_i_ratio_count/2;
00195 
00196     if (n_middle_idx >1)
00197         n_ratio_average =( an_ratio[n_middle_idx-1] +an_ratio[n_middle_idx])/2; // use median
00198     else
00199         n_ratio_average = an_ratio[n_middle_idx ];
00200 
00201     if( n_ratio_average>2 && n_ratio_average <184){
00202         n_spo2_calc= uch_spo2_table[n_ratio_average] ;
00203         *pn_spo2 = n_spo2_calc ;
00204         *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
00205     }
00206     else{
00207         *pn_spo2 =  -999 ; // do not use SPO2 since signal ratio is out of range
00208         *pch_spo2_valid  = 0; 
00209     }
00210 }
00211 
00212 
00213 void spo2_algorithm::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)
00214 /**
00215 * \brief        Find peaks
00216 * \par          Details
00217 *               Find at most MAX_NUM peaks above MIN_HEIGHT separated by at least MIN_DISTANCE
00218 *
00219 * \retval       None
00220 */
00221 {
00222     maxim_peaks_above_min_height( pn_locs, pn_npks, pn_x, n_size, n_min_height );
00223     maxim_remove_close_peaks( pn_locs, pn_npks, pn_x, n_min_distance );
00224     *pn_npks = min( *pn_npks, n_max_num );
00225 }
00226 
00227 void spo2_algorithm::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)
00228 /**
00229 * \brief        Find peaks above n_min_height
00230 * \par          Details
00231 *               Find all peaks above MIN_HEIGHT
00232 *
00233 * \retval       None
00234 */
00235 {
00236     int32_t i = 1, n_width;
00237     *pn_npks = 0;
00238     
00239     while (i < n_size-1){
00240         if (pn_x[i] > n_min_height && pn_x[i] > pn_x[i-1]){            // find left edge of potential peaks
00241             n_width = 1;
00242             while (i+n_width < n_size && pn_x[i] == pn_x[i+n_width])    // find flat peaks
00243                 n_width++;
00244             if (pn_x[i] > pn_x[i+n_width] && (*pn_npks) < 15 ){                            // find right edge of peaks
00245                 pn_locs[(*pn_npks)++] = i;        
00246                 // for flat peaks, peak location is left edge
00247                 i += n_width+1;
00248             }
00249             else
00250                 i += n_width;
00251         }
00252         else
00253             i++;
00254     }
00255 }
00256 
00257 
00258 void spo2_algorithm::maxim_remove_close_peaks(int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x,int32_t n_min_distance)
00259 /**
00260 * \brief        Remove peaks
00261 * \par          Details
00262 *               Remove peaks separated by less than MIN_DISTANCE
00263 *
00264 * \retval       None
00265 */
00266 {
00267     
00268     int32_t i, j, n_old_npks, n_dist;
00269     
00270     /* Order peaks from large to small */
00271     maxim_sort_indices_descend( pn_x, pn_locs, *pn_npks );
00272 
00273     for ( i = -1; i < *pn_npks; i++ ){
00274         n_old_npks = *pn_npks;
00275         *pn_npks = i+1;
00276         for ( j = i+1; j < n_old_npks; j++ ){
00277             n_dist =  pn_locs[j] - ( i == -1 ? -1 : pn_locs[i] ); // lag-zero peak of autocorr is at index -1
00278             if ( n_dist > n_min_distance || n_dist < -n_min_distance )
00279                 pn_locs[(*pn_npks)++] = pn_locs[j];
00280         }
00281     }
00282 
00283     // Resort indices longo ascending order
00284     maxim_sort_ascend( pn_locs, *pn_npks );
00285 }
00286 
00287 void spo2_algorithm::maxim_sort_ascend(int32_t *pn_x,int32_t n_size) 
00288 /**
00289 * \brief        Sort array
00290 * \par          Details
00291 *               Sort array in ascending order (insertion sort algorithm)
00292 *
00293 * \retval       None
00294 */
00295 {
00296     int32_t i, j, n_temp;
00297     for (i = 1; i < n_size; i++) {
00298         n_temp = pn_x[i];
00299         for (j = i; j > 0 && n_temp < pn_x[j-1]; j--)
00300             pn_x[j] = pn_x[j-1];
00301         pn_x[j] = n_temp;
00302     }
00303 }
00304 
00305 void spo2_algorithm::maxim_sort_indices_descend(int32_t *pn_x, int32_t *pn_indx, int32_t n_size)
00306 /**
00307 * \brief        Sort indices
00308 * \par          Details
00309 *               Sort indices according to descending order (insertion sort algorithm)
00310 *
00311 * \retval       None
00312 */ 
00313 {
00314     int32_t i, j, n_temp;
00315     for (i = 1; i < n_size; i++) {
00316         n_temp = pn_indx[i];
00317         for (j = i; j > 0 && pn_x[n_temp] > pn_x[pn_indx[j-1]]; j--)
00318             pn_indx[j] = pn_indx[j-1];
00319         pn_indx[j] = n_temp;
00320     }
00321 }