algorithm.cpp

00001 /** \file algorithm.cpp ******************************************************
00002 *
00003 * Project: MAXREFDES117#
00004 * Filename: algorithm.cpp
00005 * Description: This module calculates the heart rate/SpO2 level
00006 *
00007 *
00008 * --------------------------------------------------------------------
00009 *
00010 * This code follows the following naming conventions:
00011 *
00012 * char              ch_pmod_value
00013 * char (array)      s_pmod_s_string[16]
00014 * float             f_pmod_value
00015 * int32_t           n_pmod_value
00016 * int32_t (array)   an_pmod_value[16]
00017 * int16_t           w_pmod_value
00018 * int16_t (array)   aw_pmod_value[16]
00019 * uint16_t          uw_pmod_value
00020 * uint16_t (array)  auw_pmod_value[16]
00021 * uint8_t           uch_pmod_value
00022 * uint8_t (array)   auch_pmod_buffer[16]
00023 * uint32_t          un_pmod_value
00024 * int32_t *         pn_pmod_value
00025 *
00026 * ------------------------------------------------------------------------- */
00027 /*******************************************************************************
00028 * Copyright (C) 2016 Maxim Integrated Products, Inc., All Rights Reserved.
00029 *
00030 * Permission is hereby granted, free of charge, to any person obtaining a
00031 * copy of this software and associated documentation files (the "Software"),
00032 * to deal in the Software without restriction, including without limitation
00033 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
00034 * and/or sell copies of the Software, and to permit persons to whom the
00035 * Software is furnished to do so, subject to the following conditions:
00036 *
00037 * The above copyright notice and this permission notice shall be included
00038 * in all copies or substantial portions of the Software.
00039 *
00040 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
00041 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
00042 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
00043 * IN NO EVENT SHALL MAXIM INTEGRATED BE LIABLE FOR ANY CLAIM, DAMAGES
00044 * OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
00045 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
00046 * OTHER DEALINGS IN THE SOFTWARE.
00047 *
00048 * Except as contained in this notice, the name of Maxim Integrated
00049 * Products, Inc. shall not be used except as stated in the Maxim Integrated
00050 * Products, Inc. Branding Policy.
00051 *
00052 * The mere transfer of this software does not imply any licenses
00053 * of trade secrets, proprietary technology, copyrights, patents,
00054 * trademarks, maskwork rights, or any other form of intellectual
00055 * property whatsoever. Maxim Integrated Products, Inc. retains all
00056 * ownership rights.
00057 *******************************************************************************
00058 */
00059 #include "algorithm.h"
00060 #include "mbed.h"
00061 
00062 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, 
00063                               int32_t *pn_heart_rate, int8_t  *pch_hr_valid)
00064 /**
00065 * \brief        Calculate the heart rate and SpO2 level
00066 * \par          Details
00067 *               By detecting  peaks of PPG cycle and corresponding AC/DC of red/infra-red signal, the ratio for the SPO2 is computed.
00068 *               Since this algorithm is aiming for Arm M0/M3. formaula for SPO2 did not achieve the accuracy due to register overflow.
00069 *               Thus, accurate SPO2 is precalculated and save longo uch_spo2_table[] per each ratio.
00070 *
00071 * \param[in]    *pun_ir_buffer           - IR sensor data buffer
00072 * \param[in]    n_ir_buffer_length      - IR sensor data buffer length
00073 * \param[in]    *pun_red_buffer          - Red sensor data buffer
00074 * \param[out]    *pn_spo2                - Calculated SpO2 value
00075 * \param[out]    *pch_spo2_valid         - 1 if the calculated SpO2 value is valid
00076 * \param[out]    *pn_heart_rate          - Calculated heart rate value
00077 * \param[out]    *pch_hr_valid           - 1 if the calculated heart rate value is valid
00078 *
00079 * \retval       None
00080 */
00081 {
00082     uint32_t un_ir_mean ,un_only_once ;
00083     int32_t k ,n_i_ratio_count;
00084     int32_t i,s ,m, n_exact_ir_valley_locs_count ,n_middle_idx;
00085     int32_t n_th1, n_npks,n_c_min;      
00086     int32_t an_ir_valley_locs[15] ;
00087     int32_t an_exact_ir_valley_locs[15] ;
00088     int32_t an_dx_peak_locs[15] ;
00089     int32_t n_peak_interval_sum;
00090     
00091     int32_t n_y_ac, n_x_ac;
00092     int32_t n_spo2_calc; 
00093     int32_t n_y_dc_max, n_x_dc_max; 
00094     int32_t n_y_dc_max_idx, n_x_dc_max_idx; 
00095     int32_t an_ratio[5],n_ratio_average; 
00096     int32_t n_nume,  n_denom ;
00097     // remove DC of ir signal    
00098     un_ir_mean =0; 
00099     for (k=0 ; k<n_ir_buffer_length ; k++ ) un_ir_mean += pun_ir_buffer[k] ;
00100     un_ir_mean =un_ir_mean/n_ir_buffer_length ;
00101     for (k=0 ; k<n_ir_buffer_length ; k++ )  an_x[k] =  pun_ir_buffer[k] - un_ir_mean ; 
00102     
00103     // 4 pt Moving Average
00104     for(k=0; k< BUFFER_SIZE-MA4_SIZE; k++){
00105         n_denom= ( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3]);
00106         an_x[k]=  n_denom/(int32_t)4; 
00107     }
00108 
00109     // get difference of smoothed IR signal
00110     
00111     for( k=0; k<BUFFER_SIZE-MA4_SIZE-1;  k++)
00112         an_dx[k]= (an_x[k+1]- an_x[k]);
00113 
00114     // 2-pt Moving Average to an_dx
00115     for(k=0; k< BUFFER_SIZE-MA4_SIZE-2; k++){
00116         an_dx[k] =  ( an_dx[k]+an_dx[k+1])/2 ;
00117     }
00118     
00119     // hamming window
00120     // flip wave form so that we can detect valley with peak detector
00121     for ( i=0 ; i<BUFFER_SIZE-HAMMING_SIZE-MA4_SIZE-2 ;i++){
00122         s= 0;
00123         for( k=i; k<i+ HAMMING_SIZE ;k++){
00124             s -= an_dx[k] *auw_hamm[k-i] ; 
00125                      }
00126         an_dx[i]= s/ (int32_t)1146; // divide by sum of auw_hamm 
00127     }
00128 
00129  
00130     n_th1=0; // threshold calculation
00131     for ( k=0 ; k<BUFFER_SIZE-HAMMING_SIZE ;k++){
00132         n_th1 += ((an_dx[k]>0)? an_dx[k] : ((int32_t)0-an_dx[k])) ;
00133     }
00134     n_th1= n_th1/ ( BUFFER_SIZE-HAMMING_SIZE);
00135     // peak location is acutally index for sharpest location of raw signal since we flipped the signal         
00136     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 
00137 
00138     n_peak_interval_sum =0;
00139     if (n_npks>=2){
00140         for (k=1; k<n_npks; k++)
00141             n_peak_interval_sum += (an_dx_peak_locs[k]-an_dx_peak_locs[k -1]);
00142         n_peak_interval_sum=n_peak_interval_sum/(n_npks-1);
00143         *pn_heart_rate=(int32_t)(6000/n_peak_interval_sum);// beats per minutes
00144         *pch_hr_valid  = 1;
00145     }
00146     else  {
00147         *pn_heart_rate = -999;
00148         *pch_hr_valid  = 0;
00149     }
00150             
00151     for ( k=0 ; k<n_npks ;k++)
00152         an_ir_valley_locs[k]=an_dx_peak_locs[k]+HAMMING_SIZE/2; 
00153 
00154 
00155     // raw value : RED(=y) and IR(=X)
00156     // we need to assess DC and AC value of ir and red PPG. 
00157     for (k=0 ; k<n_ir_buffer_length ; k++ )  {
00158         an_x[k] =  pun_ir_buffer[k] ; 
00159         an_y[k] =  pun_red_buffer[k] ; 
00160     }
00161 
00162     // find precise min near an_ir_valley_locs
00163     n_exact_ir_valley_locs_count =0; 
00164     for(k=0 ; k<n_npks ;k++){
00165         un_only_once =1;
00166         m=an_ir_valley_locs[k];
00167         n_c_min= 16777216;//2^24;
00168         if (m+5 <  BUFFER_SIZE-HAMMING_SIZE  && m-5 >0){
00169             for(i= m-5;i<m+5; i++)
00170                 if (an_x[i]<n_c_min){
00171                     if (un_only_once >0){
00172                        un_only_once =0;
00173                    } 
00174                    n_c_min= an_x[i] ;
00175                    an_exact_ir_valley_locs[k]=i;
00176                 }
00177             if (un_only_once ==0)
00178                 n_exact_ir_valley_locs_count ++ ;
00179         }
00180     }
00181     if (n_exact_ir_valley_locs_count <2 ){
00182        *pn_spo2 =  -999 ; // do not use SPO2 since signal ratio is out of range
00183        *pch_spo2_valid  = 0; 
00184        return;
00185     }
00186     // 4 pt MA
00187     for(k=0; k< BUFFER_SIZE-MA4_SIZE; k++){
00188         an_x[k]=( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3])/(int32_t)4;
00189         an_y[k]=( an_y[k]+an_y[k+1]+ an_y[k+2]+ an_y[k+3])/(int32_t)4;
00190     }
00191 
00192     //using an_exact_ir_valley_locs , find ir-red DC andir-red AC for SPO2 calibration ratio
00193     //finding AC/DC maximum of raw ir * red between two valley locations
00194     n_ratio_average =0; 
00195     n_i_ratio_count =0; 
00196     
00197     for(k=0; k< 5; k++) an_ratio[k]=0;
00198     for (k=0; k< n_exact_ir_valley_locs_count; k++){
00199         if (an_exact_ir_valley_locs[k] > BUFFER_SIZE ){             
00200             *pn_spo2 =  -999 ; // do not use SPO2 since valley loc is out of range
00201             *pch_spo2_valid  = 0; 
00202             return;
00203         }
00204     }
00205     // find max between two valley locations 
00206     // and use ratio betwen AC compoent of Ir & Red and DC compoent of Ir & Red for SPO2 
00207 
00208     for (k=0; k< n_exact_ir_valley_locs_count-1; k++){
00209         n_y_dc_max= -16777216 ; 
00210         n_x_dc_max= - 16777216; 
00211         if (an_exact_ir_valley_locs[k+1]-an_exact_ir_valley_locs[k] >10){
00212             for (i=an_exact_ir_valley_locs[k]; i< an_exact_ir_valley_locs[k+1]; i++){
00213                 if (an_x[i]> n_x_dc_max) {n_x_dc_max =an_x[i];n_x_dc_max_idx =i; }
00214                 if (an_y[i]> n_y_dc_max) {n_y_dc_max =an_y[i];n_y_dc_max_idx=i;}
00215             }
00216             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
00217             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])  ; 
00218         
00219         
00220             n_y_ac=  an_y[n_y_dc_max_idx] - n_y_ac;    // subracting linear DC compoenents from raw 
00221             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
00222             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]); 
00223             n_x_ac=  an_x[n_y_dc_max_idx] - n_x_ac;      // subracting linear DC compoenents from raw 
00224             n_nume=( n_y_ac *n_x_dc_max)>>7 ; //prepare X100 to preserve floating value
00225             n_denom= ( n_x_ac *n_y_dc_max)>>7;
00226             if (n_denom>0  && n_i_ratio_count <5 &&  n_nume != 0)
00227             {   
00228                 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) ;
00229                 n_i_ratio_count++;
00230             }
00231         }
00232     }
00233 
00234     maxim_sort_ascend(an_ratio, n_i_ratio_count);
00235     n_middle_idx= n_i_ratio_count/2;
00236 
00237     if (n_middle_idx >1)
00238         n_ratio_average =( an_ratio[n_middle_idx-1] +an_ratio[n_middle_idx])/2; // use median
00239     else
00240         n_ratio_average = an_ratio[n_middle_idx ];
00241 
00242     if( n_ratio_average>2 && n_ratio_average <184){
00243         n_spo2_calc= uch_spo2_table[n_ratio_average] ;
00244         *pn_spo2 = n_spo2_calc ;
00245         *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
00246     }
00247     else{
00248         *pn_spo2 =  -999 ; // do not use SPO2 since signal ratio is out of range
00249         *pch_spo2_valid  = 0; 
00250     }
00251 }
00252 
00253 
00254 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)
00255 /**
00256 * \brief        Find peaks
00257 * \par          Details
00258 *               Find at most MAX_NUM peaks above MIN_HEIGHT separated by at least MIN_DISTANCE
00259 *
00260 * \retval       None
00261 */
00262 {
00263     maxim_peaks_above_min_height( pn_locs, pn_npks, pn_x, n_size, n_min_height );
00264     maxim_remove_close_peaks( pn_locs, pn_npks, pn_x, n_min_distance );
00265     *pn_npks = min( *pn_npks, n_max_num );
00266 }
00267 
00268 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)
00269 /**
00270 * \brief        Find peaks above n_min_height
00271 * \par          Details
00272 *               Find all peaks above MIN_HEIGHT
00273 *
00274 * \retval       None
00275 */
00276 {
00277     int32_t i = 1, n_width;
00278     *pn_npks = 0;
00279     
00280     while (i < n_size-1){
00281         if (pn_x[i] > n_min_height && pn_x[i] > pn_x[i-1]){            // find left edge of potential peaks
00282             n_width = 1;
00283             while (i+n_width < n_size && pn_x[i] == pn_x[i+n_width])    // find flat peaks
00284                 n_width++;
00285             if (pn_x[i] > pn_x[i+n_width] && (*pn_npks) < 15 ){                            // find right edge of peaks
00286                 pn_locs[(*pn_npks)++] = i;        
00287                 // for flat peaks, peak location is left edge
00288                 i += n_width+1;
00289             }
00290             else
00291                 i += n_width;
00292         }
00293         else
00294             i++;
00295     }
00296 }
00297 
00298 
00299 void maxim_remove_close_peaks(int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x,int32_t n_min_distance)
00300 /**
00301 * \brief        Remove peaks
00302 * \par          Details
00303 *               Remove peaks separated by less than MIN_DISTANCE
00304 *
00305 * \retval       None
00306 */
00307 {
00308     
00309     int32_t i, j, n_old_npks, n_dist;
00310     
00311     /* Order peaks from large to small */
00312     maxim_sort_indices_descend( pn_x, pn_locs, *pn_npks );
00313 
00314     for ( i = -1; i < *pn_npks; i++ ){
00315         n_old_npks = *pn_npks;
00316         *pn_npks = i+1;
00317         for ( j = i+1; j < n_old_npks; j++ ){
00318             n_dist =  pn_locs[j] - ( i == -1 ? -1 : pn_locs[i] ); // lag-zero peak of autocorr is at index -1
00319             if ( n_dist > n_min_distance || n_dist < -n_min_distance )
00320                 pn_locs[(*pn_npks)++] = pn_locs[j];
00321         }
00322     }
00323 
00324     // Resort indices longo ascending order
00325     maxim_sort_ascend( pn_locs, *pn_npks );
00326 }
00327 
00328 void maxim_sort_ascend(int32_t *pn_x,int32_t n_size) 
00329 /**
00330 * \brief        Sort array
00331 * \par          Details
00332 *               Sort array in ascending order (insertion sort algorithm)
00333 *
00334 * \retval       None
00335 */
00336 {
00337     int32_t i, j, n_temp;
00338     for (i = 1; i < n_size; i++) {
00339         n_temp = pn_x[i];
00340         for (j = i; j > 0 && n_temp < pn_x[j-1]; j--)
00341             pn_x[j] = pn_x[j-1];
00342         pn_x[j] = n_temp;
00343     }
00344 }
00345 
00346 void maxim_sort_indices_descend(int32_t *pn_x, int32_t *pn_indx, int32_t n_size)
00347 /**
00348 * \brief        Sort indices
00349 * \par          Details
00350 *               Sort indices according to descending order (insertion sort algorithm)
00351 *
00352 * \retval       None
00353 */ 
00354 {
00355     int32_t i, j, n_temp;
00356     for (i = 1; i < n_size; i++) {
00357         n_temp = pn_indx[i];
00358         for (j = i; j > 0 && pn_x[n_temp] > pn_x[pn_indx[j-1]]; j--)
00359             pn_indx[j] = pn_indx[j-1];
00360         pn_indx[j] = n_temp;
00361     }
00362 }
00363