Skip to content

Commit

Permalink
Merge pull request #15 from aromring/HighHeartRate
Browse files Browse the repository at this point in the history 
Bug fix: incorrect heartbeat rates in some cases
  • Loading branch information
@aromring
aromring authored May 9, 2020
2 parents 42299b9 + f0e88cd commit e69b603
Show file tree
Hide file tree
Showing 2 changed files with 69 additions and 14 deletions.
73 changes: 66 additions & 7 deletions algorithm_by_RF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ void rf_heart_rate_and_oxygen_saturation(uint32_t *pun_ir_buffer, int32_t n_ir_b
*/
{
int32_t k;
static int32_t n_last_peak_interval=INIT_INTERVAL;
static int32_t n_last_peak_interval=LOWEST_PERIOD;
float f_ir_mean,f_red_mean,f_ir_sumsq,f_red_sumsq;
float f_y_ac, f_x_ac, xy_ratio;
float beta_ir, beta_red, x;
Expand Down Expand Up @@ -90,15 +90,24 @@ void rf_heart_rate_and_oxygen_saturation(uint32_t *pun_ir_buffer, int32_t n_ir_b

// Calculate Pearson correlation between red and IR
*correl=rf_Pcorrelation(an_x, an_y, n_ir_buffer_length)/sqrt(f_red_sumsq*f_ir_sumsq);

// Find signal periodicity
if(*correl>=min_pearson_correlation) {
// At the beginning of oximetry run the exact range of heart rate is unknown. This may lead to wrong rate if the next call does not find the _first_
// peak of the autocorrelation function. E.g., second peak would yield only 50% of the true rate.
if(LOWEST_PERIOD==n_last_peak_interval)
rf_initialize_periodicity_search(an_x, BUFFER_SIZE, &n_last_peak_interval, HIGHEST_PERIOD, min_autocorrelation_ratio, f_ir_sumsq);
// RF, If correlation os good, then find average periodicity of the IR signal. If aperiodic, return periodicity of 0
rf_signal_periodicity(an_x, BUFFER_SIZE, &n_last_peak_interval, LOWEST_PERIOD, HIGHEST_PERIOD, min_autocorrelation_ratio, f_ir_sumsq, ratio);
if(n_last_peak_interval!=0)
rf_signal_periodicity(an_x, BUFFER_SIZE, &n_last_peak_interval, LOWEST_PERIOD, HIGHEST_PERIOD, min_autocorrelation_ratio, f_ir_sumsq, ratio);
} else n_last_peak_interval=0;

// Calculate heart rate if periodicity detector was successful. Otherwise, reset peak interval to its initial value and report error.
if(n_last_peak_interval!=0) {
*pn_heart_rate = (int32_t)(FS60/n_last_peak_interval);
*pch_hr_valid = 1;
} else {
n_last_peak_interval=INIT_INTERVAL;
n_last_peak_interval=LOWEST_PERIOD;
*pn_heart_rate = -999; // unable to calculate because signal looks aperiodic
*pch_hr_valid = 0;
*pn_spo2 = -999 ; // do not use SPO2 from this corrupt signal
Expand Down Expand Up @@ -153,6 +162,56 @@ float rf_autocorrelation(float *pn_x, int32_t n_size, int32_t n_lag)
return sum/n_temp;
}

void rf_initialize_periodicity_search(float *pn_x, int32_t n_size, int32_t *p_last_periodicity, int32_t n_max_distance, float min_aut_ratio, float aut_lag0)
/**
* \brief Search the range of true signal periodicity
* \par Details
* Determine the range of current heart rate by locating neighborhood of
* the _first_ peak of the autocorrelation function. If at all lags until
* n_max_distance the autocorrelation is less than min_aut_ratio fraction
* of the autocorrelation at lag=0, then the input signal is insufficiently
* periodic and probably indicates motion artifacts.
* Robert Fraczkiewicz, 04/25/2020
* \retval Average distance between peaks
*/
{
int32_t n_lag;
float aut,aut_right;
// At this point, *p_last_periodicity = LOWEST_PERIOD. Start walking to the right,
// two steps at a time, until lag ratio fulfills quality criteria or HIGHEST_PERIOD
// is reached.
n_lag=*p_last_periodicity;
aut_right=aut=rf_autocorrelation(pn_x, n_size, n_lag);
// Check sanity
if(aut/aut_lag0 >= min_aut_ratio) {
// Either quality criterion, min_aut_ratio, is too low, or heart rate is too high.
// Are we on autocorrelation's downward slope? If yes, continue to a local minimum.
// If not, continue to the next block.
do {
aut=aut_right;
n_lag+=2;
aut_right=rf_autocorrelation(pn_x, n_size, n_lag);
} while(aut_right/aut_lag0 >= min_aut_ratio && aut_right<aut && n_lag<=n_max_distance);
if(n_lag>n_max_distance) {
// This should never happen, but if does return failure
*p_last_periodicity=0;
return;
}
aut=aut_right;
}
// Walk to the right.
do {
aut=aut_right;
n_lag+=2;
aut_right=rf_autocorrelation(pn_x, n_size, n_lag);
} while(aut_right/aut_lag0 < min_aut_ratio && n_lag<=n_max_distance);
if(n_lag>n_max_distance) {
// This should never happen, but if does return failure
*p_last_periodicity=0;
} else
*p_last_periodicity=n_lag;
}

void rf_signal_periodicity(float *pn_x, int32_t n_size, int32_t *p_last_periodicity, int32_t n_min_distance, int32_t n_max_distance, float min_aut_ratio, float aut_lag0, float *ratio)
/**
* \brief Signal periodicity
Expand All @@ -177,9 +236,9 @@ void rf_signal_periodicity(float *pn_x, int32_t n_size, int32_t *p_last_periodic
aut=aut_left;
n_lag--;
aut_left=rf_autocorrelation(pn_x, n_size, n_lag);
} while(aut_left>aut && n_lag>n_min_distance);
} while(aut_left>aut && n_lag>=n_min_distance);
// Restore lag of the highest aut
if(n_lag==n_min_distance) {
if(n_lag<n_min_distance) {
left_limit_reached=true;
n_lag=*p_last_periodicity;
aut=aut_save;
Expand All @@ -191,9 +250,9 @@ void rf_signal_periodicity(float *pn_x, int32_t n_size, int32_t *p_last_periodic
aut=aut_right;
n_lag++;
aut_right=rf_autocorrelation(pn_x, n_size, n_lag);
} while(aut_right>aut && n_lag<n_max_distance);
} while(aut_right>aut && n_lag<=n_max_distance);
// Restore lag of the highest aut
if(n_lag==n_max_distance) n_lag=0; // Indicates failure
if(n_lag>n_max_distance) n_lag=0; // Indicates failure
else n_lag--;
if(n_lag==*p_last_periodicity && left_limit_reached) n_lag=0; // Indicates failure
}
Expand Down
10 changes: 3 additions & 7 deletions algorithm_by_RF.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,13 +48,9 @@
// The sum is symmetrc, so you can evaluate it by multiplying its positive half by 2. It is precalcuated here for enhanced
// performance.
const float sum_X2 = 83325; // WARNING: you MUST recalculate this sum if you changed either ST or FS above!
#define MAX_HR 125 // Maximal heart rate. To eliminate erroneous signals, calculated HR should never be greater than this number.
// WARNING: The two parameters below are CRUCIAL! Proper HR evaluation depends on these.
#define MAX_HR 180 // Maximal heart rate. To eliminate erroneous signals, calculated HR should never be greater than this number.
#define MIN_HR 40 // Minimal heart rate. To eliminate erroneous signals, calculated HR should never be lower than this number.
// Typical heart rate. Set it to the upper value of the expected heart rate range in a given application. Obviously, it must be
// in between MIN_HR and MAX_HR. For example, if HR in an overnight measurement varies between 46 and 65, but 90% of the time
// stays between 50 and 60, then set it to 60.
// WARNING: This is a CRUCIAL parameter! Proper HR evaluation depends on it.
#define TYPICAL_HR 60
// Minimal ratio of two autocorrelation sequence elements: one at a considered lag to the one at lag 0.
// Good quality signals must have such ratio greater than this minimum.
const float min_autocorrelation_ratio = 0.5;
Expand All @@ -71,7 +67,6 @@ const int32_t BUFFER_SIZE = FS*ST; // Number of smaples in a single batch
const int32_t FS60 = FS*60; // Conversion factor for heart rate from bps to bpm
const int32_t LOWEST_PERIOD = FS60/MAX_HR; // Minimal distance between peaks
const int32_t HIGHEST_PERIOD = FS60/MIN_HR; // Maximal distance between peaks
const int32_t INIT_INTERVAL = FS60/TYPICAL_HR; // Seed value for heart rate determination routine
const float mean_X = (float)(BUFFER_SIZE-1)/2.0; // Mean value of the set of integers from 0 to BUFFER_SIZE-1. For ST=4 and FS=25 it's equal to 49.5.

void rf_heart_rate_and_oxygen_saturation(uint32_t *pun_ir_buffer, int32_t n_ir_buffer_length, uint32_t *pun_red_buffer, float *pn_spo2, int8_t *pch_spo2_valid, int32_t *pn_heart_rate,
Expand All @@ -80,6 +75,7 @@ float rf_linear_regression_beta(float *pn_x, float xmean, float sum_x2);
float rf_autocorrelation(float *pn_x, int32_t n_size, int32_t n_lag);
float rf_rms(float *pn_x, int32_t n_size, float *sumsq);
float rf_Pcorrelation(float *pn_x, float *pn_y, int32_t n_size);
void rf_initialize_periodicity_search(float *pn_x, int32_t n_size, int32_t *p_last_periodicity, int32_t n_max_distance, float min_aut_ratio, float aut_lag0);
void rf_signal_periodicity(float *pn_x, int32_t n_size, int32_t *p_last_periodicity, int32_t n_min_distance, int32_t n_max_distance, float min_aut_ratio, float aut_lag0, float *ratio);

#endif /* ALGORITHM_BY_RF_H_ */
Expand Down

0 comments on commit e69b603

Please sign in to comment.