优化心率检测算法

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 {

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;
+}

app/src/app_service/ecg_service/algo/qrs/HC_Chen_detect.h

@@ -0,0 +1,53 @@
+#ifndef __HC_CHEN__
+#define __HC_CHEN__
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <math.h>
+#include <stdint.h>
+
+
+#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

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 ¨C T) ¨C 0.9219y(nT ¨C 2T) + x (nT) ¨C x(nT ¨C 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 ¨C T) ¨C y(nT ¨C 2T) + x(nT) ¨C 2x(nT ¨C 6T) + x(nT ¨C 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 ¨C 16T) ¨C 32 [y(nT ¨C T) + x(nT) ¨C x(nT ¨C 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 ¨C (N ¨C 1)T) + x(nT ¨C (N ¨C 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;
+	 }

app/src/app_service/ecg_service/algo/qrs/Pan_Tompkins_detect.h

@@ -0,0 +1,21 @@
+#ifndef __PAN_TOMPKINS__
+#define __PAN_TOMPKINS__
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdbool.h>
+#include <math.h>
+
+#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

app/src/app_service/ecg_service/algo/qrs/QRS.h

@@ -0,0 +1,30 @@
+#pragma once
+#include <stdint.h>
+//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 */   
+};

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;
+}

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 <stdlib.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <math.h>
+
+#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

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 {

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);

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)