diff --git a/PULP/Makefile b/PULP/Makefile index 603266c..32cb8cd 100755 --- a/PULP/Makefile +++ b/PULP/Makefile @@ -1,22 +1,23 @@ PULP_APP = adaptive_rpeak_detection PULP_APP_FC_SRCS = main.c \ adaptive_Rpeak_detection/adaptiveRpeakDetection.c \ Morph_filt/morpho_filtering.c \ REWARD_R_peak_detection/relativeEnergy.c \ REWARD_R_peak_detection/peakDetection.c \ error_detection/error_detection.c \ profiling/profile.c \ profiling/profile_cl.c -PULP_APP_CL_OMP_SRCS = test_double_buffering.c +PULP_APP_CL_OMP_SRCS = kmeans_clustering/kmean/k_mean_functions.c \ + kmeans_clustering/search/search_functions.c CORES ?= 1 CTARGET ?= 1 stackSize ?= 2048 PULP_CFLAGS += -DTARGET=$(CTARGET) -DNUM_CORES=$(CORES) -O3 -g3 -w PULP_LDFLAGS = -lm #-lg #uncomment -lg for pulpissimo include $(PULP_SDK_HOME)/install/rules/pulp_rt.mk diff --git a/PULP/adaptive_Rpeak_detection/adaptiveRpeakDetection.c b/PULP/adaptive_Rpeak_detection/adaptiveRpeakDetection.c index cffa1ef..e4ed5ae 100755 --- a/PULP/adaptive_Rpeak_detection/adaptiveRpeakDetection.c +++ b/PULP/adaptive_Rpeak_detection/adaptiveRpeakDetection.c @@ -1,377 +1,407 @@ #include "rt/rt_api.h" #include "adaptiveRpeakDetection.h" #include "../defines.h" #include "../profiling/profile.h" #include "../profiling/profile_cl.h" #include "../profiling/defines.h" #include "../data/signal.h" #include "../Morph_filt/morpho_filtering.h" #include "../Morph_filt/defines_globals.h" #include "../error_detection/error_detection.h" -#include "../test_double_buffering.h" +#include "../kmeans_clustering/kmean/k_mean_functions.h" #define N_WINDOWS (int) (2*((ECG_VECTOR_SIZE-LONG_WINDOW)/DIM)+1)// Counting the worst case scenario when overlap is DIM RT_L2_DATA int16_t ecg_L2buff[(LONG_WINDOW+DIM)*(NLEADS+1)]; RT_L2_DATA rt_perf_t perf[NUM_CORES]; #ifdef MODULE_MF RT_L2_DATA int32_t *argMF[4]; RT_L2_DATA int32_t buffSize_windowMF; #endif #ifdef MODULE_RELEN RT_L2_DATA int32_t *argRelEn[4]; RT_L2_DATA int32_t start_RelEn = 1; RT_L2_DATA int32_t buffSize_windowRelEn; #endif #ifdef MODULE_RPEAK_REWARD RT_L2_DATA int32_t *argRW_Rpeak[3]; RT_L2_DATA int32_t indicesRpeaks[H_B+1]; #endif #ifdef MODULE_ERROR_DETECTION RT_L2_DATA int32_t *argErrDet[5]; RT_L2_DATA int32_t lastRpeak = 0; RT_L2_DATA int32_t lastRR = 0; RT_L2_DATA int32_t error_RWindow = 0; #endif #ifdef MODULE_CLUSTERING RT_L2_DATA int32_t rL2BufferIndex; RT_L2_DATA int32_t rL1BufferIndex = 0; RT_L1_DATA int16_t ecg_L1buff[DIM*(NLEADS+1)]; -RT_L1_DATA int32_t start_index_w = 0; -RT_L1_DATA int32_t end_main_loop; -RT_L2_DATA int32_t* argCL[3]; +RT_L2_DATA int32_t start_index_buff = 0; +RT_L2_DATA int32_t end_main_loop; +RT_L2_DATA int32_t flag_prev_error = 0; +RT_L2_DATA int32_t* argCL[7]; +RT_L1_DATA int32_t* indicesRpeaksCL[H_B+1]; +RT_L2_DATA int32_t rpeaks_counter_cl; RT_L2_DATA rt_event_sched_t * psched = 0; RT_L2_DATA int32_t done = 0; #endif RT_L2_DATA int32_t overlap; RT_L2_DATA int32_t rWindow; void clearRelEn() { clearAndResetRelEn(); resetPeakDetection(); } -// static void cluster_Rpeaks(int32_t *arg[]) -// { -// rt_team_fork(NUM_CORES, rpeaks, arg); -// } +#ifdef MODULE_CLUSTERING -static void cluster_test_doublebuff(int32_t *arg[]) +static void cluster_Rpeaks(int32_t *arg[]) { - rt_team_fork(NUM_CORES, testDoubleBuff, arg); + rt_team_fork(NUM_CORES, rpeaks, arg); } extern void end_of_call(void *arg) { done = 1; } static void fCore0_DmaTransfer_Window(void *arg) { rt_dma_copy_t dmaCp; // Copy data block window from L2 to shared L1 memory using the cluster DMA rt_dma_memcpy((unsigned int)&ecg_L2buff[rL2BufferIndex], (unsigned int)&ecg_L1buff[rL1BufferIndex], 2*DIM, RT_DMA_DIR_EXT2LOC, 0, &dmaCp); // Wait for dma to finish - rt_dma_wait(&dmaCp); - - // printf("rL1BufferIndex: %d endL1BufferIndex: %d\n",rL1BufferIndex,rL1BufferIndex+DIM ); - // for(int i=rL1BufferIndex; i= ECG_VECTOR_SIZE){ return; } if(rWindow > 0){ offset_window = LONG_WINDOW; } else{ offset_window = 0; } if (rWindow > 0) { for(int32_t i=0; i 0) start_RelEn = 0; argRelEn[0] = (int32_t*) &ecg_L2buff[offset_window + overlap]; argRelEn[1] = (int32_t*) &ecg_L2buff[offset_window + (LONG_WINDOW + DIM)*NLEADS+overlap]; buffSize_windowRelEn = LONG_WINDOW - offset_window + DIM-overlap; #endif relEn_w(argRelEn); #ifdef HWPERF_MODULES profile_stop(perf); #endif #ifdef PRINT_RELEN for(int32_t sample = offset_window + (LONG_WINDOW + DIM)*NLEADS; sample< (LONG_WINDOW + DIM)*(NLEADS+1); sample++) { printf("%d\n", ecg_L2buff[sample]); } #endif #endif #ifdef MODULE_RPEAK_REWARD #ifdef HWPERF_MODULE_RPEAK_REWARD || HWPERF_MODULES profile_start(perf); #endif getPeaks_w(argRW_Rpeak); rpeaks_counter = 0; while(indicesRpeaks[rpeaks_counter]!=0) { rpeaks_counter++; } #ifdef HWPERF_MODULE_RPEAK_REWARD || HWPERF_MODULES profile_stop(perf); #endif #ifdef PRINT_RPEAKS for(int32_t indR=0; indR 0 && error_RWindow == 1){ - #endif - // ----------------------------Copy current window ecg buffer from L2 to L1 memory if error was 1 ------------------------------ // + #else + if(rWindow == 0){ + rL1BufferIndex = 0; + start_index_buff = 0; + // -------------------------------- Copy first window ecg buffer from L2 to L1 memory ------------------------------------------ // // Initialize event event = rt_event_get_blocking(NULL); // Run function on Core 0 of the cluster rt_cluster_call(NULL, 0, fCore0_DmaTransfer_Window, NULL, NULL, STACK_SIZE, STACK_SIZE, 1, event); // Wait for event rt_event_wait(event); // ------------------------------------------------------------------------------------------------------------------------------// + }else{ + rL1BufferIndex = DIM; + start_index_buff = DIM; + #endif + // ----------------Copy current window ecg buffer from L2 to L1 memory if error was 1 or after first window--------------------- // + // Initialize event + event = rt_event_get_blocking(NULL); + // Run function on Core 0 of the cluster + rt_cluster_call(NULL, 0, fCore0_DmaTransfer_Window, NULL, NULL, STACK_SIZE, STACK_SIZE, 1, event); + + // Wait for event + rt_event_wait(event); + // ------------------------------------------------------------------------------------------------------------------------------// argCL[0] = (int32_t*) ecg_L1buff; //The start index is the one in argCL[1] - argCL[1] = &start_index_w; - argCL[2] = &end_main_loop; - rt_cluster_call(NULL, CID, cluster_test_doublebuff, argCL, NULL, 2048, 2048, NUM_CORES, rt_event_get(psched, end_of_call, (void *) CID)); + argCL[1] = &rWindow; + argCL[2] = &start_index_buff; + argCL[3] = &end_main_loop; + argCL[4] = &offset_ind; + argCL[5] = &flag_prev_error; + argCL[6] = indicesRpeaksCL; + rt_cluster_call(NULL, CID, cluster_Rpeaks, argCL, NULL, STACK_SIZE, STACK_SIZE, NUM_CORES, rt_event_get(psched, end_of_call, (void *) CID)); while(!done) rt_event_execute(psched, 1); done = 0; - #ifdef MODULE_ERROR_DETECTION + + // Move last sample to index DIM-1 of L1 buffer for the next window (this is necessary if the clustering module runs without error detection + // or if the previous error was 1 and the clustering must keep running) + // The clustering uses the approximated signal derivative (a diff function), so it needs the last sample of the previous window to not miss anything + ecg_L1buff[DIM-1] = ecg_L1buff[end_main_loop]; + + flag_prev_error = 1; + + #ifdef PRINT_RPEAKS_CL + rpeaks_counter_cl = 0; + + while(indicesRpeaksCL[rpeaks_counter_cl]!=0) { + rpeaks_counter_cl++; + } + + for(int ix_rr=0; ix_rr +#include "../search/search_functions.h" +#include "rt/rt_api.h" + +type_f genlogfunc(type_f x, type_f k, type_f b, type_f nu, type_f m){ + return (k/powf(1+expf(-b*(x-m)),(1/nu))); +} + + +type_f gaussian(type_f x, type_f mu, type_f sign){ + return expf(-powf((x - mu), 2)/(2*powf(sign, 2))); +} + +/* Function to sort an array using insertion sort*/ +void insertionSort(type_f arr[], int n) +{ + type_f key; + int i, j; + for (i = 1; i < n; i++) { + key = arr[i]; + j = i - 1; + + /* Move elements of arr[0..i-1], that are + greater than key, to one position ahead + of their current position */ + while (j >= 0 && arr[j] > key) { + arr[j + 1] = arr[j]; + j = j - 1; + } + arr[j + 1] = key; + } +} + +RT_L1_DATA type_f abs_s2_init[SIZE_PERCENTILE_ARR]; //it's the first 1000 samples of abs diff to sort and find the 99 percentile + +type_f percentile(int16_t *x, int end, int percToFind){ + + int indexPerc = (end-1)*percToFind/100; + for(int i = 0; i 0 && i==1){ + // hcentr = centroids_state[0]; + // lcentr = centroids_state[1]; + // h_el = (type_i) centroids_state[2]; + // l_el = (type_i) centroids_state[3]; + // } + //====== Uncomment if you are using end_last_peak for overlap ======// + + //====== Uncomment if you are using qrs_init for overlap ======// +// if(number_of_window > 0 && i==1 && last_peak > start_index_w){ +// last_peak = rpks[num_peaks-2]; +// mu = rpks[num_peaks-2] - rpks[num_peaks-3]; +// sd = std(rpks, num_peaks-2, num_peaks_total-2); +// #ifdef PRINT_DEBUG_CL_STD +// printf("last_peak: %d mu: %f sd: %f\n",last_peak, mu, sd); +// #endif +// if(sd < 10) +// sd = 10; +// if(sd > 50) +// sd = 50; + +// num_peaks = num_peaks-1; +// last_peak_w = last_peak; +// } + //====== Uncomment if you are using qrs_init for overlap ======// +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check) + printf("last_peak: %d input_value_gauss: %d mu: %f sd: %f\n",last_peak, start_index_w + (i-1) - last_peak, mu, sd); +#endif + bf = gaussian(start_index_w + (i-1) - last_peak, mu, sd); // i - last_peak_w +#ifdef PRINT_GAUSS + // DEBUG GAUSSIAN + printf("%f,", bf); + #endif +#ifdef PRINT_GAUSS_MU_SD + printf("%f %f %f \n", bf,mu,sd); +#endif + param_log = logf((hcentr - lcentr)/2)/(2*lcentr); + + param_log2 = log2f(hcentr/lcentr); + + bt = genlogfunc(x , hcentr, param_log, 1/param_log2, lcentr); + +#ifdef PRINT_GENLOGFUN + printf("%f,", bt); +#endif + + new_val = (bt * bf); + + if(x > new_val){ + partial_fit(&x, &hcentr, &lcentr, &h_el, &l_el, &label); +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check) + printf("index_all: %d i: %d x: %f hcentr: %f lcentr: %f label: %d\n", start_index_w + (i-1), i, x, hcentr, lcentr, label); +#endif +#ifdef PRINT_BAYFILT + // DEBUG BAYESIAN FILTERED SIGNAL + printf("%f,", x); +#endif + } + else{ + partial_fit(&new_val, &hcentr, &lcentr, &h_el, &l_el, &label); +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check) + printf("index_all: %d i: %d new_val: %f hcentr: %f lcentr: %f label: %d\n", start_index_w + (i-1), i, new_val, hcentr, lcentr, label); +#endif +#ifdef PRINT_BAYFILT + // DEBUG BAYESIAN FILTERED SIGNAL + printf("%f,", new_val); +#endif + } +#ifdef PRINT_CENTROIDS + // DEBUG CENTROIDS + printf("%f %f\n", lcentr,hcentr); +#endif + + if(in_qrs == 1){ + if(label == 0){ + zeroctr += 1; + }else{ + zeroctr = 0; + } + + s2[index_st] = (type_f) (ecg_buff[i] - ecg_buff[i-1]); + if(x > new_val){ + st[index_st] = x; + index_st++; + }else{ + st[index_st] = new_val; + index_st++; + } + +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check) + printf("qrs_init: %d zeroctr: %d qrs_init_total: %d\n", qrs_init, zeroctr, qrs_init_total); +#endif + // Once a QRS has been detected, we look for 120ms of 0 output (30 samples) or 200ms to consider that the QRS has finished. + if(zeroctr==30 || start_index_w + (i-1) - qrs_init_total > max_qrs_dur){ + +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check) + printf("start interval: %d stop interval: %d ecg_buff[start_interval]: %f i: %d zeroctr: %d label: %d index_st: %d stop_interval_st: %d \n",qrs_init, (i-1) - zeroctr + 1,ecg_buff[qrs_init],i, zeroctr, label,index_st-1,(index_st-1)-zeroctr+1 ); +#endif + + // Search for location of the new peak. Checking maximum and minimum slopes, and move towards the peak + end_qrs = (index_st-1)-zeroctr+1; + argMaxMinProductSignSearch(st, s2, 0, end_qrs,max_min_slope_productSign); + int max_slope = max_min_slope_productSign[0]; + int min_slope = max_min_slope_productSign[1]; + int max_slope_s2 = 0; + int min_slope_s2 = 0; + if (max_slope < min_slope){ + new_peak = qrs_init + max_slope + argMaxProductSignSearch(st, s2, max_slope, min_slope+1); + }else if (max_slope > min_slope){ + argMaxMinSearch(s2, 0, end_qrs, max_min_slope); + max_slope_s2 = max_min_slope[0]; + min_slope_s2 = max_min_slope[1]; + if(max_slope_s2 < min_slope_s2){ + // Q wave + new_peak = qrs_init + max_slope_s2 + argMaxSearch(s2, max_slope_s2, min_slope_s2+1); + }else{ + // S wave + new_peak = qrs_init + argMaxSearch(s2, 0, min_slope_s2+1); + } + }else{ + new_peak = qrs_init + argMaxSearch(st, 0, end_qrs); + } + +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check){ + printf("BEFORE CORRECTION new_peak: %d\n", new_peak); + if (max_slope > min_slope){ + printf("end_qrs: %d max_slope: %d min_slope: %d max_slope_s2: %d min_slope_s2: %d\n",end_qrs,max_slope,min_slope,max_slope_s2,min_slope_s2); + }else{ + printf("end_qrs: %d max_slope: %d min_slope: %d\n",end_qrs,max_slope,min_slope); + } + } +#endif + + if(ecg_buff[new_peak+1] > ecg_buff[new_peak]){ + while(ecg_buff[new_peak+1] > ecg_buff[new_peak]){ + new_peak += 1; + } + } + else{ + while(ecg_buff[new_peak] < ecg_buff[new_peak-1]){ + new_peak -= 1; + } + } + + +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check) + printf("AFTER CORRECTION new_peak: %d\n", new_peak); +#endif + + last_peak = new_peak + start_index_w; + if(last_peak > last_peak_w){ + if(last_peak!=rpks[num_peaks]){ + indicesRpeaks[ind_peak] = last_peak; + ind_peak++; +#ifdef PRINT_RPEAKS_DEBUG + printf("%d\n", last_peak); +#endif + } + if(num_peaks == MAX_PEAKS_IN_WIND){ + for(int ix=0; ix 1){ + mu = rpks[num_peaks]-rpks[num_peaks-1]; + sd = std(rpks, num_peaks, num_peaks_total); + #ifdef PRINT_DEBUG_CL_STD + if(number_of_window==num_wind_to_check) + printf("mu: %f sd: %f\n", mu, sd); + #endif + if(sd < 10) + sd = 10; + if(sd > 50) + sd = 50; + } + num_peaks++; + } + last_peak_w = last_peak; + end_last_peak = i-1; + centroids_state[0] = hcentr; + centroids_state[1] = lcentr; + centroids_state[2] = (type_f) h_el; + centroids_state[3] = (type_f) l_el; + zeroctr = 0; + in_qrs = 0; + index_st = 0; + } + + + }else{ + if(label==1 && start_index_w + (i-1) > last_peak + min_rr_dist){ + in_qrs = 1; + qrs_init = i-1; + qrs_init_total = start_index_w + (i-1); + s2[index_st] = (type_f) (ecg_buff[i] - ecg_buff[i-1]); + if(x > new_val){ + st[index_st] = x; + index_st++; + }else{ + st[index_st] = new_val; + index_st++; + } + } + } + } + +#ifdef PRINT_DEBUG_CL + if(number_of_window==num_wind_to_check-1) + printf("END OF WINDOW qrs_init: %d qrs_init_total: %d last_peak: %d mu: %f std: %f index_st: %d new_peak: %d\n",qrs_init,qrs_init_total,last_peak,mu,sd, index_st-1,new_peak); +#endif + +} + diff --git a/PULP/kmeans_clustering/kmean/k_mean_functions.h b/PULP/kmeans_clustering/kmean/k_mean_functions.h new file mode 100644 index 0000000..79c91fe --- /dev/null +++ b/PULP/kmeans_clustering/kmean/k_mean_functions.h @@ -0,0 +1,14 @@ +#include +#include +#include "../../defines.h" + +#define PEAKS_STD 5 +#define MAX_PEAKS_IN_WIND (N*15 + PEAKS_STD) + +type_f ganlogfunc(type_f x, type_f k, type_f b, type_f nu, type_f m); +type_f gaussian(type_f x, type_f mu, type_f sign); +type_f percentile(int16_t *x, int end, int percToFind); +void partial_fit(type_f *x, type_f *hcentr, type_f *lcentr, int *h_el, int *l_el, int *label); +type_f std(type_i *x, int num_peaks, int num_peaks_total); + +void rpeaks(int32_t* arg[]); diff --git a/PULP/kmeans_clustering/search/search_functions.c b/PULP/kmeans_clustering/search/search_functions.c new file mode 100644 index 0000000..bca7413 --- /dev/null +++ b/PULP/kmeans_clustering/search/search_functions.c @@ -0,0 +1,81 @@ +#include +#include "search_functions.h" + +type_i argMaxSearch(type_f* x, type_i ind_start, type_i ind_end){ + type_i argmax = 0; + type_f max_value = x[ind_start]; + for(type_i ix = ind_start+1; ix < ind_end; ix++){ + if(x[ix]>max_value){ + max_value = x[ix]; + argmax = ix; + } + } + return argmax; +} + +void argMaxMinSearch(type_f* x, type_i ind_start, type_i ind_end, type_i* argmax_argmin){ + argmax_argmin[0] = ind_start; + argmax_argmin[1] = ind_start; + type_f max_value = x[ind_start]; + type_f min_value = x[ind_start]; + for(type_i ix = ind_start+1; ix < ind_end; ix++){ + if(x[ix]>max_value){ + max_value = x[ix]; + argmax_argmin[0] = ix; + } + if(x[ix]0) return 1; + if(x<0) return -1; + + return 0; +} + +type_i argMaxProductSignSearch(type_f* x, type_f* x2, type_i ind_start, type_i ind_end){ + type_i argmax = 0; + type_f max_value = 0; + type_f product_arrays = 0; + + max_value = x[ind_start]*sign(x2[ind_start]); + + for(type_i ix = ind_start+1; ix < ind_end; ix++){ + product_arrays = x[ix]*sign(x2[ix]); + + if(product_arrays>max_value){ + max_value = product_arrays; + argmax = ix; + } + } + return argmax; +} + +void argMaxMinProductSignSearch(type_f* x, type_f* x2, type_i ind_start, type_i ind_end, type_i* argmax_argmin){ + argmax_argmin[0] = ind_start; + argmax_argmin[1] = ind_start; + type_f max_value = 0; + type_f min_value = 0; + type_f product_arrays = 0; + + max_value = x[ind_start]*sign(x2[ind_start]); + min_value = x[ind_start]*sign(x2[ind_start]); + + for(type_i ix = ind_start+1; ix < ind_end; ix++){ + product_arrays = x[ix]*sign(x2[ix]); + + if(product_arrays>max_value){ + max_value = product_arrays; + argmax_argmin[0] = ix; + } + + if(product_arrays - -void testDoubleBuff(int32_t* arg[]){ - int16_t* ecg_buff = (int16_t*) arg[0]; // pointing at ecg_buff from 0 and then using start_index_w and end_main_loop to use the correct part of the buffer - int32_t start_index_w = *(arg[1]); - int32_t end_main_loop = *(arg[2]); - - // printf("start_index_w: %d end_main_loop: %d\n",start_index_w,end_main_loop ); - // for(int i=start_index_w; i -#include "defines.h" - -void testDoubleBuff(int32_t* arg[]); \ No newline at end of file