diff --git a/app/src/app_service/ecg_service/algo/iflytop_simple_filter.c b/app/src/app_service/ecg_service/algo/iflytop_simple_filter.c index e653e86..fcf1fec 100644 --- a/app/src/app_service/ecg_service/algo/iflytop_simple_filter.c +++ b/app/src/app_service/ecg_service/algo/iflytop_simple_filter.c @@ -170,11 +170,6 @@ float median_filter_update(median_filter_t *filter, float val) { } } - // 4 --> 1,2 - // 3 --> 1 - // 2 --> 0,1 - // 1 --> 0 - if (filter->numVal % 2 == 0) { return (filter->temp[filter->numVal / 2 - 1] + filter->temp[filter->numVal / 2]) / 2; } else { diff --git a/app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.c b/app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.c new file mode 100644 index 0000000..5bfcd67 --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.c @@ -0,0 +1,128 @@ +#include "HC_Chen_detect.h" + +bool HC_Chen_detect(float signal) +{ + ecg_buff[ecg_buff_WR_idx++] = signal; + sample = ecg_buff_WR_idx; + ecg_buff_WR_idx %= (M+1); + + /* High pass filtering */ + if(number_iter < M){ + // first fill buffer with enough points for HP filter + hp_sum += ecg_buff[ecg_buff_RD_idx]; + hp_buff[hp_buff_WR_idx] = 0; + } + else{ + hp_sum += ecg_buff[ecg_buff_RD_idx]; + + int tmp = ecg_buff_RD_idx - M; + if(tmp < 0){ + tmp += M + 1; + } + + hp_sum -= ecg_buff[tmp]; + + float y1 = 0; + float y2 = 0; + + tmp = (ecg_buff_RD_idx - ((M+1)/2)); + if(tmp < 0){ + tmp += M + 1; + } + + y2 = ecg_buff[tmp]; + + y1 = HP_CONSTANT * hp_sum; + + hp_buff[hp_buff_WR_idx] = y2 - y1; + + } + + // done reading ECG buffer, increment position + ecg_buff_RD_idx++; + ecg_buff_RD_idx %= (M+1); + + // done writing to HP buffer, increment position + hp_buff_WR_idx++; + hp_buff_WR_idx %= (N+1); + + /* Low pass filtering */ + + // shift in new sample from high pass filter + lp_sum += hp_buff[hp_buff_RD_idx] * hp_buff[hp_buff_RD_idx]; + + if(number_iter < N){ + // first fill buffer with enough points for LP filter + next_eval_pt = 0; + } + else{ + // shift out oldest data point + int tmp = hp_buff_RD_idx - N; + if(tmp < 0){ + tmp += N+1; + } + + lp_sum -= hp_buff[tmp] * hp_buff[tmp]; + + next_eval_pt = lp_sum; + } + + // done reading HP buffer, increment position + hp_buff_RD_idx++; + hp_buff_RD_idx %= (N+1); + + /* Adapative thresholding beat detection */ + // set initial threshold + if(number_iter < window_size) { + if(next_eval_pt > treshold) { + treshold = next_eval_pt; + } + ++number_iter; + } + + // check if detection hold off period has passed + if(triggered){ + trig_time++; + + if(trig_time >= DELAY_TIME){ + triggered = false; + trig_time = 0; + } + } + + // find if we have a new max + if(next_eval_pt > win_max) win_max = next_eval_pt; + + // find if we are above adaptive threshold + if(next_eval_pt > treshold && !triggered) { + //result.push_back(true); + + last_qrs_point = sample; + + triggered = true; + return true; + } + else { + //result.push_back(false); + } + + // adjust adaptive threshold using max of signal found + // in previous window + if(win_idx++ >= window_size){ + // weighting factor for determining the contribution of + // the current peak value to the threshold adjustment + float gamma = (0.2f+0.15f)/2.0f; // 0.15~0.2 + + // forgetting factor - + // rate at which we forget old observations + float alpha = 0.01f + ( ((float) rand() / (float) RAND_MAX) * ((0.1f - 0.01f))); // 0~1 + //float alpha = 1.0f*exp(-0.00005f*(sample - last_qrs_point)); + + treshold = alpha * gamma * win_max + (1.0f - alpha) * treshold; + + // reset current window ind + win_idx = 0; + win_max = -10000000; + } + return false; +} diff --git a/app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.h b/app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.h new file mode 100644 index 0000000..cb9fa51 --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.h @@ -0,0 +1,53 @@ +#ifndef __HC_CHEN__ +#define __HC_CHEN__ + +#include +#include +#include +#include +#include + + +#include "QRS.h" + + +#define M 9 +#define N 54//SAMPLING_RATE * 0.15f +static const uint32_t window_size = SAMPLING_RATE; +static const float HP_CONSTANT = ((float) 1.0f / (float) M); +// circular buffer for input ecg signal +// we need to keep a history of M + 1 samples for HP filter +static float ecg_buff[M + 1] = {0}; +static int ecg_buff_WR_idx = 0; +static int ecg_buff_RD_idx = 0; + +// circular buffer for input ecg signal +// we need to keep a history of N+1 samples for LP filter +static float hp_buff[N + 1] = {0}; +static int hp_buff_WR_idx = 0; +static int hp_buff_RD_idx = 0; + +// LP filter outputs a single point for every input point +// This goes straight to adaptive filtering for eval +static float next_eval_pt = 0; + +// running sums for HP and LP filters, values shifted in FILO +static float hp_sum = 0; +static float lp_sum = 0; + +// parameters for adaptive thresholding +static float treshold = 0; +static bool triggered = false; +static int trig_time = 0; +static float win_max = 0; +static int win_idx = 0; +static int number_iter = 0; + +static int sample = 0; +static int last_qrs_point = 0; + +static const int DELAY_TIME = 180;//window_size * 0.5f; + +extern bool HC_Chen_detect(float); + +#endif diff --git a/app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.c b/app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.c new file mode 100644 index 0000000..35c1b9f --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.c @@ -0,0 +1,373 @@ +#include "Pan_Tompkins_detect.h" + +/* y(nT) = 1.875y(nT 每 T) 每 0.9219y(nT 每 2T) + x (nT) 每 x(nT 每 2T) */ +int TwoPoleRecursive(int data) +{ + static int xnt, xm1, xm2, ynt, ym1, ym2 = 0; + + xnt = data; + ynt = (ym1 + (ym1 >> 1) + (ym1 >> 2) + (ym1 >> 3)) + // 1.875 = 1 + 1/2 + 1/4 + 1/8 + (((ym2 >> 1) + (ym2 >> 2) + (ym2 >> 3) + (ym2 >> 5) + (ym2 >> 6)) + xnt - xm2); // 0.916 = 1 + 1/2 + 1/4 + 1/8 + 1/32 + 1/64 + xm2 = xm1; + xm1 = xnt; + xm2 = ym1; + ym2 = ym1; + ym1 = ynt; + return ynt; +} + +/* y(nT) = 2y(nT 每 T) 每 y(nT 每 2T) + x(nT) 每 2x(nT 每 6T) + x(nT 每 12T) */ +int LowPassFilter(int data) +{ + static int y1 = 0, y2 = 0, x[26], n = 12; + int y0; + + x[n] = x[n + 13] = data; + y0 = (y1 << 1) - y2 + x[n] - (x[n + 6] << 1) + x[n + 12]; + y2 = y1; + y1 = y0; + y0 >>= 5; + if(--n < 0){ + n = 12; + } + return y0; +} + +/* p(nT) = x(nT 每 16T) 每 32 [y(nT 每 T) + x(nT) 每 x(nT 每 32T)] */ +int HighPassFilter(int data) +{ + static int y1 = 0, x[66], n = 32; + int y0; + + x[n] = x[n + 33] = data; + y0 = y1 + x[n] - x[n + 32]; + y1 = y0; + if(--n < 0){ + n = 32; + } + return (x[n + 16] - (y0 >> 5)); +} + +/* y = 1/8 (2x( nT) + x( nT - T) - x( nT - 3T) - 2x( nT - 4T)) */ +int Derivative(int data) +{ + int y; + static int x_derv[4]; + + y = (data << 1) + x_derv[3] - x_derv[1] - ( x_derv[0] << 1); + y >>= 3; + for(int i = 0; i < 3; ++i){ + x_derv[i] = x_derv[i + 1]; + } + x_derv[3] = data; + return y; +} + +int Squar(int data) +{ + return (data * data); +} + +/* y(nT) = 1/N [x(nT 每 (N 每 1)T) + x(nT 每 (N 每 2)T) +...+ x(nT)] */ +int MovingWindowIntegral(int data) +{ + //static const int WINDOW_SIZE = SAMPLING_RATE * 0.2; + #define WINDOW_SIZE 72 + + static int x[WINDOW_SIZE], ptr = 0; + static long sum = 0; + long ly; + int y; + + if(++ptr == WINDOW_SIZE){ + ptr = 0; + } + sum -= x[ptr]; + sum += data; + x[ptr] = data; + ly = sum >> 5; + uint32_t MAX_INTEGRAL = 4096;//32400; + if(ly > MAX_INTEGRAL){ + y = MAX_INTEGRAL; + } + else{ + y = (int)ly; + } + return (y); +} + +SignalPoint ThresholdCalculate(int sample,float value,int bandpass,int square,int integral) +{ + //static const int QRS_TIME = SAMPLING_RATE * 0.1; + //static const int SEARCH_BACK_TIME = SAMPLING_RATE * 1.66f; + #define QRS_TIME 36 + #define SEARCH_BACK_TIME 598 + + static int bandpass_buffer[SEARCH_BACK_TIME],integral_buffer[SEARCH_BACK_TIME]; + static SignalPoint peak_buffer[SEARCH_BACK_TIME]; + static int square_buffer[QRS_TIME]; + static long unsigned last_qrs = 0, last_slope = 0, current_slope = 0; + static int peak_i = 0, peak_f = 0, threshold_i1 = 0, threshold_i2 = 0, threshold_f1 = 0, threshold_f2 = 0, spk_i = 0, spk_f = 0, npk_i = 0, npk_f = 0; + static bool qrs, regular = true, prev_regular; + static int rr1[8]={0}, rr2[8]={0}, rravg1, rravg2, rrlow = 0, rrhigh = 0, rrmiss = 0; + + SignalPoint result; + result.index = -1; + + peak_buffer[sample%SEARCH_BACK_TIME].index = sample; + peak_buffer[sample%SEARCH_BACK_TIME].value = value; + bandpass_buffer[sample%SEARCH_BACK_TIME] = bandpass; + integral_buffer[sample%SEARCH_BACK_TIME] = integral; + square_buffer[sample%QRS_TIME] = square; + + // If the current signal is above one of the thresholds (integral or filtered signal), it's a peak candidate. + if(integral >= threshold_i1 || bandpass >= threshold_f1){ + peak_i = integral; + peak_f = bandpass; + } + + // If both the integral and the signal are above their thresholds, they're probably signal peaks. + if((integral >= threshold_i1) && (bandpass >= threshold_f1)){ + // There's a 200ms latency. If the new peak respects this condition, we can keep testing. + if(sample > last_qrs + SAMPLING_RATE*0.2f){ + //if(sample > last_qrs + (SAMPLING_RATE*0.2f)){ + // If it respects the 200ms latency, but it doesn't respect the 360ms latency, we check the slope. + if(sample <= last_qrs + (long unsigned int)(0.36*SAMPLING_RATE)){ + // The squared slope is "M" shaped. So we have to check nearby samples to make sure we're really looking + // at its peak value, rather than a low one. + int current = sample; + current_slope = 0; + for(int j = current - QRS_TIME; j <= current; ++j){ + if(square_buffer[j%QRS_TIME] > current_slope){ + current_slope = square_buffer[j%QRS_TIME]; + } + } + //current_slope = square; + + if(current_slope <= (int)(last_slope/2)){ + qrs = false; + //return qrs; + } + + else{ + spk_i = 0.125*peak_i + 0.875*spk_i; + threshold_i1 = npk_i + 0.25*(spk_i - npk_i); + threshold_i2 = 0.5*threshold_i1; + + spk_f = 0.125*peak_f + 0.875*spk_f; + threshold_f1 = npk_f + 0.25*(spk_f - npk_f); + threshold_f2 = 0.5*threshold_f1; + + last_slope = current_slope; + qrs = true; + + result.value = value; + result.index = sample; + } + } + // If it was above both thresholds and respects both latency periods, it certainly is a R peak. + else{ + int current = sample; + current_slope = 0; + for(int j = current - QRS_TIME; j <= current; ++j){ + if(square_buffer[j%QRS_TIME] > current_slope){ + current_slope = square_buffer[j%QRS_TIME]; + } + } + //current_slope = square; + + spk_i = 0.125*peak_i + 0.875*spk_i; + threshold_i1 = npk_i + 0.25*(spk_i - npk_i); + threshold_i2 = 0.5*threshold_i1; + + spk_f = 0.125*peak_f + 0.875*spk_f; + threshold_f1 = npk_f + 0.25*(spk_f - npk_f); + threshold_f2 = 0.5*threshold_f1; + + last_slope = current_slope; + qrs = true; + + result.value = value; + result.index = sample; + } + } + // If the new peak doesn't respect the 200ms latency, it's noise. Update thresholds and move on to the next sample. + else{ + peak_i = integral; + npk_i = 0.125*peak_i + 0.875*npk_i; + threshold_i1 = npk_i + 0.25*(spk_i - npk_i); + threshold_i2 = 0.5*threshold_i1; + peak_f = bandpass; + npk_f = 0.125*peak_f + 0.875*npk_f; + threshold_f1 = npk_f + 0.25*(spk_f - npk_f); + threshold_f2 = 0.5*threshold_f1; + qrs = false; + /*outputSignal[current] = qrs; + if (sample > DELAY + BUFFSIZE) + output(outputSignal[0]); + continue;*/ + + //return qrs; + return result; + } + } + + // If a QRS complex was detected, the RR-averages must be updated. + if(qrs){ + // Add the newest RR-interval to the buffer and get the new average. + rravg1 = 0; + for (int i = 0; i < 7; ++i){ + rr1[i] = rr1[i+1]; + rravg1 += rr1[i]; + } + rr1[7] = sample - last_qrs; + last_qrs = sample; + rravg1 += rr1[7]; + rravg1 *= 0.125; + + // If the newly-discovered RR-average is normal, add it to the "normal" buffer and get the new "normal" average. + // Update the "normal" beat parameters. + if ( (rr1[7] >= rrlow) && (rr1[7] <= rrhigh) ){ + rravg2 = 0; + for (int i = 0; i < 7; ++i){ + rr2[i] = rr2[i+1]; + rravg2 += rr2[i]; + } + rr2[7] = rr1[7]; + rravg2 += rr2[7]; + rravg2 *= 0.125; + rrlow = 0.92*rravg2; + rrhigh = 1.16*rravg2; + rrmiss = 1.66*rravg2; + } + + prev_regular = regular; + if(rravg1 == rravg2){ + regular = true; + } + // If the beat had been normal but turned odd, change the thresholds. + else{ + regular = false; + if (prev_regular){ + threshold_i1 /= 2; + threshold_f1 /= 2; + } + } + } + // If no R-peak was detected, it's important to check how long it's been since the last detection. + else{ + int current = sample; + // If no R-peak was detected for too long, use the lighter thresholds and do a back search. + // However, the back search must respect the 200ms limit and the 360ms one (check the slope). + if((sample - last_qrs > (long unsigned int)rrmiss) && (sample > last_qrs + SAMPLING_RATE*0.2f)){ + //if((sample - last_qrs > (long unsigned int)rrmiss) && (sample > last_qrs + (SAMPLING_RATE*0.2f))){ + + // If over SEARCH_BACK_TIME of QRS complex + if((sample - last_qrs) > SEARCH_BACK_TIME){ + last_qrs = sample; + //return result; + } + + int qrs_last_index = 0; // Last point of QRS complex + + for(int i = current - (sample - last_qrs) + SAMPLING_RATE*0.2f; i < (long unsigned int)current; ++i){ + //for(int i = current - (sample - last_qrs) + (SAMPLING_RATE*0.2f); i < (long unsigned int)current; ++i){ + if((integral_buffer[i%SEARCH_BACK_TIME] > threshold_i2) && (bandpass_buffer[i%SEARCH_BACK_TIME] > threshold_f2)){ + current_slope = 0; + for(int j = current - QRS_TIME; j <= current; ++j){ + if(square_buffer[j%QRS_TIME] > current_slope){ + current_slope = square_buffer[j%QRS_TIME]; + } + } + //current_slope = square; + + if((current_slope < (int)(last_slope/2)) && (i + sample) < last_qrs + 0.36*last_qrs){ + qrs = false; + } + else if(i - last_qrs > 550){ + peak_i = integral_buffer[i%SEARCH_BACK_TIME]; + peak_f = bandpass_buffer[i%SEARCH_BACK_TIME]; + spk_i = 0.25*peak_i+ 0.75*spk_i; + spk_f = 0.25*peak_f + 0.75*spk_f; + threshold_i1 = npk_i + 0.25*(spk_i - npk_i); + threshold_i2 = 0.5*threshold_i1; + last_slope = current_slope; + threshold_f1 = npk_f + 0.25*(spk_f - npk_f); + threshold_f2 = 0.5*threshold_f1; + // If a signal peak was detected on the back search, the RR attributes must be updated. + // This is the same thing done when a peak is detected on the first try. + //RR Average 1 + rravg1 = 0; + for(int j = 0; j < 7; ++j){ + rr1[j] = rr1[j+1]; + rravg1 += rr1[j]; + } + rr1[7] = sample - (current - i) - last_qrs; + qrs = true; + qrs_last_index = i; + last_qrs = sample - (current - i); + rravg1 += rr1[7]; + rravg1 *= 0.125; + + //RR Average 2 + if((rr1[7] >= rrlow) && (rr1[7] <= rrhigh)){ + rravg2 = 0; + for (int i = 0; i < 7; ++i){ + rr2[i] = rr2[i+1]; + rravg2 += rr2[i]; + } + rr2[7] = rr1[7]; + rravg2 += rr2[7]; + rravg2 *= 0.125; + rrlow = 0.92*rravg2; + rrhigh = 1.16*rravg2; + rrmiss = 1.66*rravg2; + } + + prev_regular = regular; + if(rravg1 == rravg2){ + regular = true; + } + else{ + regular = false; + if(prev_regular){ + threshold_i1 /= 2; + threshold_f1 /= 2; + } + } + + break; + } + } + } + + if(qrs){ + //outputSignal[current] = false; + //outputSignal[i] = true; + //if (sample > DELAY + BUFFSIZE) + //output(outputSignal[0]); + //continue; + + //return qrs; + return peak_buffer[qrs_last_index%SEARCH_BACK_TIME]; + } + } + + // Definitely no signal peak was detected. + if(!qrs){ + // If some kind of peak had been detected, then it's certainly a noise peak. Thresholds must be updated accordinly. + if((integral >= threshold_i1) || (bandpass >= threshold_f1)){ + peak_i = integral; + npk_i = 0.125*peak_i + 0.875*npk_i; + threshold_i1 = npk_i + 0.25*(spk_i - npk_i); + threshold_i2 = 0.5*threshold_i1; + peak_f = bandpass; + npk_f = 0.125*peak_f + 0.875*npk_f; + threshold_f1 = npk_f + 0.25*(spk_f - npk_f); + threshold_f2 = 0.5*threshold_f1; + } + } + } + + return result; + } diff --git a/app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.h b/app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.h new file mode 100644 index 0000000..cbd65d0 --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.h @@ -0,0 +1,21 @@ +#ifndef __PAN_TOMPKINS__ +#define __PAN_TOMPKINS__ + +#include +#include +#include +#include +#include + +#include "QRS.h" +extern int TwoPoleRecursive(int); + +extern int LowPassFilter(int); +extern int HighPassFilter(int); + +extern int Derivative(int); +extern int Squar(int); +extern int MovingWindowIntegral(int); + +extern SignalPoint ThresholdCalculate(int,float,int,int,int); +#endif diff --git a/app/src/app_service/ecg_service/algo/qrs/QRS.h b/app/src/app_service/ecg_service/algo/qrs/QRS.h new file mode 100644 index 0000000..33aba1b --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/QRS.h @@ -0,0 +1,30 @@ +#pragma once +#include +//const uint32_t SAMPLING_RATE = 1000; +#define SAMPLING_RATE 360 + +typedef struct +{ + float value; + int32_t index; +}SignalPoint; + +enum +{ + NOTQRS, /* not-QRS (not a getann/putann code) */ + NORMAL, /* normal beat */ + LBBB, /* left bundle branch block beat */ + RBBB, /* right bundle branch block beat */ + ABERR, /* aberrated atrial premature beat */ + PVC, /* premature ventricular contraction */ + FUSION, /* fusion of ventricular and normal beat */ + NPC, /* nodal (junctional) premature beat */ + APC, /* atrial premature contraction */ + SVPB, /* premature or ectopic supraventricular beat */ + VESC, /* ventricular escape beat */ + NESC, /* nodal (junctional) escape beat */ + PACE, /* paced beat */ + UNKNOWN, /* unclassifiable beat */ + NOISE, /* signal quality change */ + ARFCT /* isolated QRS-like artifact */ +}; diff --git a/app/src/app_service/ecg_service/algo/qrs/So_Chen_detect.c b/app/src/app_service/ecg_service/algo/qrs/So_Chen_detect.c new file mode 100644 index 0000000..e329825 --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/So_Chen_detect.c @@ -0,0 +1,93 @@ +#include "So_Chen_detect.h" + +SignalPoint So_Chen_detect(SignalPoint signal,int initial_point,float threshold_parameter,float filter_parameter) +{ + /* init slop window pool, size = 5 */ + if(signal_window_count < signal_window_size){ + signal_window[signal_window_count%signal_window_size] = signal; + ++signal_window_count; + SignalPoint value; + value.index = -1; + return value; + } + else{ + signal_window[signal_window_count%signal_window_size] = signal; + ++signal_window_count; + SignalPoint value; + } + + /* calculate slop */ + uint32_t idx_for_slop = signal_window_count-2; + slop.value = ( (-2.0f * signal_window[(idx_for_slop-2)%signal_window_size].value) - signal_window[(idx_for_slop-1)%signal_window_size].value + signal_window[(idx_for_slop+1)%signal_window_size].value + (2.0f * signal_window[(idx_for_slop+2)%signal_window_size].value) ); + slop.index = signal_window[idx_for_slop%signal_window_size].index; + + /* init maxi */ + if(!so_chen_init_flag){ + if(!maxi_init){ + max.value = 0; + max.index = -1; + maxi = slop.value; + maxi_init = true; + } + ++init_count; + if(init_count > initial_point){ + so_chen_init_flag = true; + /* calculate slop threshold */ + slop_threshold = threshold_parameter / 16.0f * maxi; + } + if(slop.value > maxi){ + maxi = slop.value; + } + SignalPoint value; + value.index = -1; + return value; + } + + /* detect QRS complex on set */ + if(qrs_on_set_flag && (signal_window_count - last_point > enhanced_point)){ + if(!max_init){ + max = signal_window[(idx_for_slop)%signal_window_size]; + max_init = true; + } + if(signal_window[(idx_for_slop)%signal_window_size].value > max.value){ + max = signal_window[(idx_for_slop)%signal_window_size]; + max_slop = slop; + } + else if(signal_window[(idx_for_slop)%signal_window_size].value < max.value){ + last_point = signal_window_count; + qrs_on_set_flag = false; + max_init = false; + maxi = ((abs(max.value - qrs_onset_point.value) - maxi) / filter_parameter) + maxi; + slop_threshold = threshold_parameter / 16.0f * maxi; + last_maxi = maxi; + return max; + } + } + else{ + if(slop.value > slop_threshold){ + ++qrs_on_set_count; + } + else if(qrs_on_set_count){ + qrs_on_set_count = 0; + } + + if(qrs_on_set_count >= 2){ // is QRS complex on set + qrs_on_set_flag = true; + qrs_on_set_count = 0; + qrs_onset_idx = idx_for_slop; + qrs_onset_point = signal; + } + else if((signal_window_count - last_point > enhanced_point * 2) && (slop_threshold > 0)){ //decay threshold + + slop_threshold -= slop.value; + + if((signal_window_count - last_point > SAMPLING_RATE * 3)){ //threshold oscillating + slop_threshold -= ((signal_window_count - last_point) >> (int)threshold_parameter); + } + } + } + + SignalPoint value; + value.index = -1; + return value; +} diff --git a/app/src/app_service/ecg_service/algo/qrs/So_Chen_detect.h b/app/src/app_service/ecg_service/algo/qrs/So_Chen_detect.h new file mode 100644 index 0000000..a505775 --- /dev/null +++ b/app/src/app_service/ecg_service/algo/qrs/So_Chen_detect.h @@ -0,0 +1,42 @@ +#ifndef __SO_AND_CHEN__ +#define __SO_AND_CHEN__ + +#include +#include +#include +#include +#include + +#include "QRS.h" +static const uint32_t enhanced_point = SAMPLING_RATE * 0.35f; + +#define signal_window_size 5 +static SignalPoint signal_window[signal_window_size]; +static uint32_t signal_window_count = 0; + +static SignalPoint slop; + +static bool so_chen_init_flag = false; +static uint32_t init_count = 0; + +static bool maxi_init = false; +static float maxi; + +static float slop_threshold = 0; + +static SignalPoint qrs_onset_point; + +static int qrs_on_set_count = 0; +static int qrs_onset_idx = 0; +static bool qrs_on_set_flag = false; + +static SignalPoint max; +static SignalPoint max_slop; +static bool max_init = false; + +static float last_maxi = 0; +static uint32_t last_point = 0; + +SignalPoint So_Chen_detect(SignalPoint,int,float,float); + +#endif diff --git a/app/src/app_service/ecg_service/algo/zsimple_qrs.c b/app/src/app_service/ecg_service/algo/zsimple_qrs.c index b0d2507..b4b23fa 100644 --- a/app/src/app_service/ecg_service/algo/zsimple_qrs.c +++ b/app/src/app_service/ecg_service/algo/zsimple_qrs.c @@ -1,5 +1,31 @@ #include "zsimple_qrs.h" +void prv_zsimple_qrs_push_heartrate(zsimple_qrs_t* p_qrs, float heartrate) { + memmove(p_qrs->heartrate_cache, p_qrs->heartrate_cache + 1, sizeof(p_qrs->heartrate_cache[0]) * (HEART_RATE_CACHE_SIZE - 1)); + p_qrs->heartrate_cache[HEART_RATE_CACHE_SIZE - 1] = heartrate; +} + +void prv_zsimple_qrs_analys_heartrate(zsimple_qrs_t* p_qrs, float* diffavg) { + float avg = 0; + float sum = 0; + for (int i = 0; i < HEART_RATE_CACHE_SIZE; i++) { + sum += p_qrs->heartrate_cache[i]; + } + + avg = sum / HEART_RATE_CACHE_SIZE; + + float diffsum = 0; + for (int i = 0; i < HEART_RATE_CACHE_SIZE; i++) { + float diff = (p_qrs->heartrate_cache[i] - avg); + if (diff < 0) { + diff = -diff; + } + diffsum += diff; + } + + *diffavg = diffsum / HEART_RATE_CACHE_SIZE; +} + void zsimple_qrs_init(zsimple_qrs_t* p_qrs, float sample_rate_s) { p_qrs->last_peak_pos = 0; p_qrs->sample_rate_s = sample_rate_s; @@ -16,8 +42,9 @@ void zsimple_qrs_clear(zsimple_qrs_t* p_qrs) { p_qrs->index = 0; p_qrs->hasfindpeak = false; median_filter_init(&p_qrs->heart_rate_filter, 10); - + memset(p_qrs->heartrate_cache, 0, sizeof(p_qrs->heartrate_cache)); } + void zsimple_qrs_process_data(zsimple_qrs_t* p_qrs, int32_t indata, float min, float max, float refavg) { // p_qrs->index++; @@ -30,7 +57,17 @@ void zsimple_qrs_process_data(zsimple_qrs_t* p_qrs, int32_t indata, float min, f p_qrs->hasfindpeak = true; p_qrs->last_peak_pos = p_qrs->index; float heartrate = 60 / peak_interval_s; - p_qrs->heartrate = median_filter_update(&p_qrs->heart_rate_filter, heartrate); + + prv_zsimple_qrs_push_heartrate(p_qrs, heartrate); + float heartrate_af_median = median_filter_update(&p_qrs->heart_rate_filter, heartrate); + + float diffavg = 0; + prv_zsimple_qrs_analys_heartrate(p_qrs, &diffavg); + if (diffavg < 10) { + p_qrs->heartrate = heartrate_af_median; + }else{ + p_qrs->heartrate = 0; + } } } } else { diff --git a/app/src/app_service/ecg_service/algo/zsimple_qrs.h b/app/src/app_service/ecg_service/algo/zsimple_qrs.h index 1a9a8bb..47b970b 100644 --- a/app/src/app_service/ecg_service/algo/zsimple_qrs.h +++ b/app/src/app_service/ecg_service/algo/zsimple_qrs.h @@ -16,6 +16,8 @@ #define HEART_RATE_BUF_SIZE 10 +#define HEART_RATE_CACHE_SIZE 5 + typedef struct { int32_t last_peak_pos; float sample_rate_s; @@ -24,6 +26,9 @@ typedef struct { bool hasfindpeak; median_filter_t heart_rate_filter; + + float heartrate_cache[HEART_RATE_CACHE_SIZE]; + } zsimple_qrs_t; void zsimple_qrs_init(zsimple_qrs_t* p_qrs, float sample_rate_s); diff --git a/app/src/aproject_config/config.h b/app/src/aproject_config/config.h index dfc5ba1..3e7bea5 100644 --- a/app/src/aproject_config/config.h +++ b/app/src/aproject_config/config.h @@ -7,7 +7,7 @@ #define CATEGORY "M1003" // 等絳薊 #define MANUFACTURER_NAME "iflytop" -#define FIRMWARE_VERSION (505) +#define FIRMWARE_VERSION (506) #define BLESTACK_VERSION 1 #define BOOTLOADER_VERSION 1 #define HARDWARE_VERSION (2)