working wavelet transform
Dependencies: CMSIS_DSP_5 include mbed
Fork of Nucleo-Heart-Rate by
main.cpp
- Committer:
- MockyBirdTwo
- Date:
- 2018-06-22
- Revision:
- 11:337c5ff16f27
- Parent:
- 10:f62efa525bb6
File content as of revision 11:337c5ff16f27:
#include "arm_common_tables.h" #include "arm_const_structs.h" #include "arm_math.h" #include "mbed.h" #include "vals.h" #include "../header/wavelib.h" DigitalOut myled(LED1); DigitalOut myled2(LED2); Serial usb_serial(SERIAL_TX, SERIAL_RX); // tx, rx const int baud_rate = 115200; // Baud rate. double* getPeaks(double* values, int window, double Tr_low) { static double result[4] = {0, 0, 0, 0}; // {Value min, index min, value max, index max} int wnd_max = 320;//Max window is 2 seconds = 2*fs = 320 int cont =0;// Variable to continue to max peak find int i,j; for(i = 0; i < wnd_max; i++) { if(values[i] < Tr_low) { //Find MIN peak below a threshold result[0] = values[i]; // Place where the threshold is met cont=1; //Continue to find MAX peak for(j=i; j<i+window; j++) {//Find MIN peak if(values[j] < result[0]) { result[0] = values[j]; //value MIN peak result[1] = (double) j; //index MIN peak } } break; } } if (cont==1) { for(i = result[1]; i < result[1]+window; i++) { result[2]=0; if (values[i]>result[2]) { result[2]=values[i]; result[3]=(double) i; } } } return result;//Will return result = {0,0,0,0} if nothing is found } int main() { usb_serial.baud(baud_rate); // Set serial USB connection baud rate (variable is declared in config part). wave_object wave_init_obj; wt_object wavelet; int N_samples = 2000; int wt_scale = 3; int i,j; static int fs =160;//Sampling freq. int wnd_offset = 0.25*fs; static int wnd=0.05*fs;//Window to find peaks in int Strtup =1;// Start with the starting sequence double Tr_low = -0.02; //Threshold MIN peak double Tr_temp = 0; //Temporary threshold MIN peak double Tr_scale = 0.45; //Starting threshold scale double *Peaks=(double*)calloc(4,sizeof(double)); //To save peak data double *PTD=(double*)calloc(3,sizeof(double)); //To save peak time difference int *Index_min=(int*)calloc(4,sizeof(int)); //To save MIN peak index int *Index_max=(int*)calloc(4,sizeof(int)); //To save MAX peak index double PTD_temp;//Temporary value of PTD double MTD2;//PTD average of last two PTD double MTD3;//PTD average of last 3 PTD double *inp = signal_in; //,*out,*diff; double *out =(double*)malloc(sizeof(double)* N_samples); double *wt_res=(double*)malloc(sizeof(double)* N_samples); double *t_wavelet=(double*)malloc(sizeof(double)* N_samples); //To save wavelet coefficients temporarily char *m_wavelet = "db4"; wave_init_obj = wave_init(m_wavelet);// Initialize the wavelet wavelet = wt_init(wave_init_obj, "modwt", N_samples, wt_scale);// Initialize the wavelet transform object modwt(wavelet, inp);// Perform MODW int wt_scale_start=2000;//Level 3 int wt_scale_end=4000; //Choose scales for (i=0,j=wt_scale_start; j<wt_scale_end; i++,j++) { t_wavelet[i]=wavelet->output[j]; } wavelet->output=t_wavelet; imodwt(wavelet,out);// Perform IMODWT for (i=0; i<N_samples; ++i) { wt_res[i]=wavelet->output[i]*fabs(wavelet->output[i]); } while (1) { if (Strtup == 1) { //Starting sequence //First heartbeat Peaks = getPeaks(wt_res, wnd, Tr_low);// Peaks={Value min, index min, value max, index max} Tr_temp = Tr_scale*Peaks[0]; //Set MIN peak threshold Index_min[0]= (int)Peaks[1]; //Save MIN peak index Index_max[0]=(int)Peaks[3]; //Save MAX peak index for (i=0; i<3; i++) { //Next 3 heartbeats Peaks = getPeaks(wt_res+(wnd_offset+Index_min[i]), wnd, Tr_temp);// Peaks={Value min, index min, value max, index max} Tr_temp = Tr_scale*Peaks[0];//Set MIN peak threshold Index_min[i+1]=wnd_offset+Index_min[i]+(int)Peaks[1];//Save MIN peak index + starting offset Index_max[i+1]=wnd_offset+Index_min[i]+(int)Peaks[3];//Save MAX peak index + starting offset PTD[i] = Index_max[i+1]-Index_max[i];//Determine PTD wnd_offset = 0.5*PTD[i]; //Set new window offset } //Set MTD2 & MTD3 MTD2 = (PTD[2]+PTD[1])/2; MTD3 = (PTD[2]+PTD[1]+PTD[0])/3; if(MTD2>(0.9*MTD3) && MTD2<(1.1*MTD3)) { //Signal stable, so heartbeat found Strtup = 0; Tr_low =Tr_temp; } else { //Unstable signal, restart startup Strtup = 1; Tr_low = -0.01; } } else if(Strtup==0) { //Normal sequence //First heartbeat Peaks = getPeaks(wt_res, wnd, Tr_low);// Peaks={Value min, index min, value max, index max} Tr_temp = Tr_scale*Peaks[0]; //Set MIN peak threshold Index_min[0]=(int)Peaks[1]; //Save MIN peak index Index_max[0]=(int)Peaks[3]; //Save MAX peak index usb_serial.printf("Normal sequence started \n"); for (i=0; i<2; i++) { Peaks = getPeaks(wt_res+(wnd_offset+Index_min[i]), wnd, Tr_temp);// Peaks={Value min, index min, value max, index max} Tr_temp = Tr_scale*Peaks[0]; //Set MIN peak threshold Index_min[i+1]=wnd_offset+Index_min[i]+(int)Peaks[1]; //Save MIN peak index Index_max[i+1]=wnd_offset+Index_min[i]+(int)Peaks[3]; //Save MAX peak index PTD_temp=Index_max[i+1]-Index_max[i]; if (PTD_temp>= (0.8*MTD2) && PTD_temp<= (1.2*MTD2)) {//Found heart rate Strtup = 0;//Continue normal sequence Tr_low =Tr_temp; //Set threshold //Pop off first PTD and add PTD_temp to the back of PTD PTD[1]=PTD[2]; PTD[0]=PTD[1]; PTD[2]=PTD_temp; //Set MTD2 & MTD3 MTD2 = (PTD[2]+PTD[1])/2; MTD3 = (PTD[2]+PTD[1]+PTD[0])/2; break; } Strtup = 1; } } else { //Error usb_serial.printf("Error in algorithm"); } break; } // while(1) { // for (i = 0; i < N_samples; ++i) { // usb_serial.printf("%d %e \n",i, wt_res[i]); // } // myled = !myled; // } // while(1) { // // for (i = 0; i < N_samples; ++i) { // usb_serial.printf("%d %e \n",i, wavelet->output[i]); // } // myled = !myled; // } //Algoritme }