MAXREFDES117# Code Documentation  V01.00
Heart Rate / SpO2 Monitor
 All Files Functions Variables Macros Pages
algorithm.cpp
Go to the documentation of this file.
1 
27 /*******************************************************************************
28 * Copyright (C) 2016 Maxim Integrated Products, Inc., All Rights Reserved.
29 *
30 * Permission is hereby granted, free of charge, to any person obtaining a
31 * copy of this software and associated documentation files (the "Software"),
32 * to deal in the Software without restriction, including without limitation
33 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
34 * and/or sell copies of the Software, and to permit persons to whom the
35 * Software is furnished to do so, subject to the following conditions:
36 *
37 * The above copyright notice and this permission notice shall be included
38 * in all copies or substantial portions of the Software.
39 *
40 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
41 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
42 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
43 * IN NO EVENT SHALL MAXIM INTEGRATED BE LIABLE FOR ANY CLAIM, DAMAGES
44 * OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
45 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
46 * OTHER DEALINGS IN THE SOFTWARE.
47 *
48 * Except as contained in this notice, the name of Maxim Integrated
49 * Products, Inc. shall not be used except as stated in the Maxim Integrated
50 * Products, Inc. Branding Policy.
51 *
52 * The mere transfer of this software does not imply any licenses
53 * of trade secrets, proprietary technology, copyrights, patents,
54 * trademarks, maskwork rights, or any other form of intellectual
55 * property whatsoever. Maxim Integrated Products, Inc. retains all
56 * ownership rights.
57 *******************************************************************************
58 */
59 #include "algorithm.h"
60 #include "mbed.h"
61 
62 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,
63  int32_t *pn_heart_rate, int8_t *pch_hr_valid)
81 {
82  uint32_t un_ir_mean ,un_only_once ;
83  int32_t k ,n_i_ratio_count;
84  int32_t i, s, m, n_exact_ir_valley_locs_count ,n_middle_idx;
85  int32_t n_th1, n_npks,n_c_min;
86  int32_t an_ir_valley_locs[15] ;
87  int32_t an_exact_ir_valley_locs[15] ;
88  int32_t an_dx_peak_locs[15] ;
89  int32_t n_peak_interval_sum;
90 
91  int32_t n_y_ac, n_x_ac;
92  int32_t n_spo2_calc;
93  int32_t n_y_dc_max, n_x_dc_max;
94  int32_t n_y_dc_max_idx, n_x_dc_max_idx;
95  int32_t an_ratio[5],n_ratio_average;
96  int32_t n_nume, n_denom ;
97  // remove DC of ir signal
98  un_ir_mean =0;
99  for (k=0 ; k<n_ir_buffer_length ; k++ ) un_ir_mean += pun_ir_buffer[k] ;
100  un_ir_mean =un_ir_mean/n_ir_buffer_length ;
101  for (k=0 ; k<n_ir_buffer_length ; k++ ) an_x[k] = pun_ir_buffer[k] - un_ir_mean ;
102 
103  // 4 pt Moving Average
104  for(k=0; k< BUFFER_SIZE-MA4_SIZE; k++){
105  n_denom= ( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3]);
106  an_x[k]= n_denom/(int32_t)4;
107  }
108 
109  // get difference of smoothed IR signal
110 
111  for( k=0; k<BUFFER_SIZE-MA4_SIZE-1; k++)
112  an_dx[k]= (an_x[k+1]- an_x[k]);
113 
114  // 2-pt Moving Average to an_dx
115  for(k=0; k< BUFFER_SIZE-MA4_SIZE-2; k++){
116  an_dx[k] = ( an_dx[k]+an_dx[k+1])/2 ;
117  }
118 
119  // hamming window
120  // flip wave form so that we can detect valley with peak detector
121  for ( i=0 ; i<BUFFER_SIZE-HAMMING_SIZE-MA4_SIZE-2 ;i++){
122  s= 0;
123  for( k=i; k<i+ HAMMING_SIZE ;k++){
124  s -= an_dx[k] *auw_hamm[k-i] ;
125  }
126  an_dx[i]= s/ (int32_t)1146; // divide by sum of auw_hamm
127  }
128 
129 
130  n_th1=0; // threshold calculation
131  for ( k=0 ; k<BUFFER_SIZE-HAMMING_SIZE ;k++){
132  n_th1 += ((an_dx[k]>0)? an_dx[k] : ((int32_t)0-an_dx[k])) ;
133  }
134  n_th1= n_th1/ ( BUFFER_SIZE-HAMMING_SIZE);
135  // peak location is acutally index for sharpest location of raw signal since we flipped the signal
136  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
137 
138  n_peak_interval_sum =0;
139  if (n_npks>=2){
140  for (k=1; k<n_npks; k++)
141  n_peak_interval_sum += (an_dx_peak_locs[k]-an_dx_peak_locs[k -1]);
142  n_peak_interval_sum=n_peak_interval_sum/(n_npks-1);
143  *pn_heart_rate=(int32_t)(6000/n_peak_interval_sum);// beats per minutes
144  *pch_hr_valid = 1;
145  }
146  else {
147  *pn_heart_rate = -999;
148  *pch_hr_valid = 0;
149  }
150 
151  for ( k=0 ; k<n_npks ;k++)
152  an_ir_valley_locs[k]=an_dx_peak_locs[k]+HAMMING_SIZE/2;
153 
154 
155  // raw value : RED(=y) and IR(=X)
156  // we need to assess DC and AC value of ir and red PPG.
157  for (k=0 ; k<n_ir_buffer_length ; k++ ) {
158  an_x[k] = pun_ir_buffer[k] ;
159  an_y[k] = pun_red_buffer[k] ;
160  }
161 
162  // find precise min near an_ir_valley_locs
163  n_exact_ir_valley_locs_count =0;
164  for(k=0 ; k<n_npks ;k++){
165  un_only_once =1;
166  m=an_ir_valley_locs[k];
167  n_c_min= 16777216;//2^24;
168  if (m+5 < BUFFER_SIZE-HAMMING_SIZE && m-5 >0){
169  for(i= m-5;i<m+5; i++)
170  if (an_x[i]<n_c_min){
171  if (un_only_once >0){
172  un_only_once =0;
173  }
174  n_c_min= an_x[i] ;
175  an_exact_ir_valley_locs[k]=i;
176  }
177  if (un_only_once ==0)
178  n_exact_ir_valley_locs_count ++ ;
179  }
180  }
181  if (n_exact_ir_valley_locs_count <2 ){
182  *pn_spo2 = -999 ; // do not use SPO2 since signal ratio is out of range
183  *pch_spo2_valid = 0;
184  return;
185  }
186  // 4 pt MA
187  for(k=0; k< BUFFER_SIZE-MA4_SIZE; k++){
188  an_x[k]=( an_x[k]+an_x[k+1]+ an_x[k+2]+ an_x[k+3])/(int32_t)4;
189  an_y[k]=( an_y[k]+an_y[k+1]+ an_y[k+2]+ an_y[k+3])/(int32_t)4;
190  }
191 
192  //using an_exact_ir_valley_locs , find ir-red DC andir-red AC for SPO2 calibration ratio
193  //finding AC/DC maximum of raw ir * red between two valley locations
194  n_ratio_average =0;
195  n_i_ratio_count =0;
196 
197  for(k=0; k< 5; k++) an_ratio[k]=0;
198  for (k=0; k< n_exact_ir_valley_locs_count; k++){
199  if (an_exact_ir_valley_locs[k] > BUFFER_SIZE ){
200  *pn_spo2 = -999 ; // do not use SPO2 since valley loc is out of range
201  *pch_spo2_valid = 0;
202  return;
203  }
204  }
205  // find max between two valley locations
206  // and use ratio betwen AC compoent of Ir & Red and DC compoent of Ir & Red for SPO2
207 
208  for (k=0; k< n_exact_ir_valley_locs_count-1; k++){
209  n_y_dc_max= -16777216 ;
210  n_x_dc_max= - 16777216;
211  if (an_exact_ir_valley_locs[k+1]-an_exact_ir_valley_locs[k] >10){
212  for (i=an_exact_ir_valley_locs[k]; i< an_exact_ir_valley_locs[k+1]; i++){
213  if (an_x[i]> n_x_dc_max) {n_x_dc_max =an_x[i];n_x_dc_max_idx =i; }
214  if (an_y[i]> n_y_dc_max) {n_y_dc_max =an_y[i];n_y_dc_max_idx=i;}
215  }
216  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
217  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]) ;
218 
219 
220  n_y_ac= an_y[n_y_dc_max_idx] - n_y_ac; // subracting linear DC compoenents from raw
221  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
222  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]);
223  n_x_ac= an_x[n_y_dc_max_idx] - n_x_ac; // subracting linear DC compoenents from raw
224  n_nume=( n_y_ac *n_x_dc_max)>>7 ; //prepare X100 to preserve floating value
225  n_denom= ( n_x_ac *n_y_dc_max)>>7;
226  if (n_denom>0 && n_i_ratio_count <5 && n_nume != 0)
227  {
228  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) ;
229  n_i_ratio_count++;
230  }
231  }
232  }
233 
234  maxim_sort_ascend(an_ratio, n_i_ratio_count);
235  n_middle_idx= n_i_ratio_count/2;
236 
237  if (n_middle_idx >1)
238  n_ratio_average =( an_ratio[n_middle_idx-1] +an_ratio[n_middle_idx])/2; // use median
239  else
240  n_ratio_average = an_ratio[n_middle_idx ];
241 
242  if( n_ratio_average>2 && n_ratio_average <184){
243  n_spo2_calc= uch_spo2_table[n_ratio_average] ;
244  *pn_spo2 = n_spo2_calc ;
245  *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
246  }
247  else{
248  *pn_spo2 = -999 ; // do not use SPO2 since signal ratio is out of range
249  *pch_spo2_valid = 0;
250  }
251 }
252 
253 
254 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)
262 {
263  maxim_peaks_above_min_height( pn_locs, pn_npks, pn_x, n_size, n_min_height );
264  maxim_remove_close_peaks( pn_locs, pn_npks, pn_x, n_min_distance );
265  *pn_npks = min( *pn_npks, n_max_num );
266 }
267 
268 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)
276 {
277  int32_t i = 1, n_width;
278  *pn_npks = 0;
279 
280  while (i < n_size-1){
281  if (pn_x[i] > n_min_height && pn_x[i] > pn_x[i-1]){ // find left edge of potential peaks
282  n_width = 1;
283  while (i+n_width < n_size && pn_x[i] == pn_x[i+n_width]) // find flat peaks
284  n_width++;
285  if (pn_x[i] > pn_x[i+n_width] && (*pn_npks) < 15 ){ // find right edge of peaks
286  pn_locs[(*pn_npks)++] = i;
287  // for flat peaks, peak location is left edge
288  i += n_width+1;
289  }
290  else
291  i += n_width;
292  }
293  else
294  i++;
295  }
296 }
297 
298 
299 void maxim_remove_close_peaks(int32_t *pn_locs, int32_t *pn_npks, int32_t *pn_x, int32_t n_min_distance)
307 {
308 
309  int32_t i, j, n_old_npks, n_dist;
310 
311  /* Order peaks from large to small */
312  maxim_sort_indices_descend( pn_x, pn_locs, *pn_npks );
313 
314  for ( i = -1; i < *pn_npks; i++ ){
315  n_old_npks = *pn_npks;
316  *pn_npks = i+1;
317  for ( j = i+1; j < n_old_npks; j++ ){
318  n_dist = pn_locs[j] - ( i == -1 ? -1 : pn_locs[i] ); // lag-zero peak of autocorr is at index -1
319  if ( n_dist > n_min_distance || n_dist < -n_min_distance )
320  pn_locs[(*pn_npks)++] = pn_locs[j];
321  }
322  }
323 
324  // Resort indices longo ascending order
325  maxim_sort_ascend( pn_locs, *pn_npks );
326 }
327 
328 void maxim_sort_ascend(int32_t *pn_x,int32_t n_size)
329 
336 {
337  int32_t i, j, n_temp;
338  for (i = 1; i < n_size; i++) {
339  n_temp = pn_x[i];
340  for (j = i; j > 0 && n_temp < pn_x[j-1]; j--)
341  pn_x[j] = pn_x[j-1];
342  pn_x[j] = n_temp;
343  }
344 }
345 
346 void maxim_sort_indices_descend(int32_t *pn_x, int32_t *pn_indx, int32_t n_size)
354 {
355  int32_t i, j, n_temp;
356  for (i = 1; i < n_size; i++) {
357  n_temp = pn_indx[i];
358  for (j = i; j > 0 && pn_x[n_temp] > pn_x[pn_indx[j-1]]; j--)
359  pn_indx[j] = pn_indx[j-1];
360  pn_indx[j] = n_temp;
361  }
362 }
363