Yu anning
/
LPC_MAX30102
the MAX30102 driver for LPC1768 board
Embed:
(wiki syntax)
Show/hide line numbers
algorithm.cpp
Go to the documentation of this file.
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
Generated on Mon Aug 1 2022 16:46:04 by 1.7.2