Chanel's edits
Dependencies: max32630fthr USBDevice
panTompkins.h
- Committer:
- saleiferis
- Date:
- 2020-03-11
- Revision:
- 16:f6bfa6b66e96
- Parent:
- 15:b15b4b6c6da8
File content as of revision 16:f6bfa6b66e96:
/** * ------------------------------------------------------------------------------* * File: panTompkins.c * * ANSI-C implementation of Pan-Tompkins real-time QRS detection algorithm * * Author: Rafael de Moura Moreira <rafaelmmoreira@gmail.com> * * License: MIT License * * ------------------------------------------------------------------------------* * ---------------------------------- HISTORY ---------------------------------- * * date | author | description * * ----------| -------------| ---------------------------------------------------* * 2019/04/11| Rafael M. M. | - Fixed moving-window integral. * * | | - Fixed how to find the correct sample with the * * | | last QRS. * * | | - Replaced constant value in code by its #define. * * | | - Added some casting on comparisons to get rid of * * | | compiler warnings. * * 2019/04/15| Rafael M. M. | - Removed delay added to the output by the filters.* * | | - Fixed multiple detection of the same peak. * * 2019/04/16| Rafael M. M. | - Added output buffer to correctly output a peak * * | | found by back searching using the 2nd thresholds. * * 2019/04/23| Rafael M. M. | - Improved comparison of slopes. * * | | - Fixed formula to obtain the correct sample from * * | | the buffer on the back search. * * ------------------------------------------------------------------------------* * MIT License * * * * Copyright (c) 2018 Rafael de Moura Moreira * * * * Permission is hereby granted, free of charge, to any person obtaining a copy * * of this software and associated documentation files (the "Software"), to deal * * in the Software without restriction, including without limitation the rights * * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * * copies of the Software, and to permit persons to whom the Software is * * furnished to do so, subject to the following conditions: * * * * The above copyright notice and this permission notice shall be included in all* * copies or substantial portions of the Software. * * * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * * SOFTWARE. * *-------------------------------------------------------------------------------* * Description * * * * The main goal of this implementation is to be easy to port to different opera-* * ting systems, as well as different processors and microcontrollers, including * * embedded systems. It can work both online or offline, depending on whether all* * the samples are available or not - it can be adjusted on the input function. * * * * The main function, panTompkings(), calls input() to get the next sample and * * store it in a buffer. Then it runs through a chain of filters: DC block, low * * pass @ 15 Hz and high pass @ 5Hz. The filtered signal goes both through a de- * * rivative filter, which output is then squared, and through a windowed-integra-* * tor. * * * * For a signal peak to be recognized as a fiducial point, its correspondent va- * * lue on both the filtered signal and the integrator must be above a certain * * threshold. Additionally, there are time-restraints to prevent a T-wave from * * being mistakenly identified as an R-peak: a hard 200ms restrain (a new peak * * 200ms from the previous one is, necessarily, a T-wave) and a soft 360ms res- * * train (the peak's squared slope must also be very high to be considered as a * * real peak). * * * * When a peak candidate is discarded, its value is used to update the noise * * thresholds - which are also used to estimate the peak thresholds. * * * * Two buffers keep 8 RR-intervals to calculate RR-averages: one of them keeps * * the last 8 RR-intervals, while the other keeps only the RR-intervals that res-* * pect certain restrictions. If both averages are equal, the heart pace is con- * * sidered normal. If the heart rate isn't normal, the thresholds change to make * * it easier to detect possible weaker peaks. If no peak is detected for a long * * period of time, the thresholds also change and the last discarded peak candi- * * date is reconsidered. * *-------------------------------------------------------------------------------* * Instructions * * * * Here's what you should change to adjust the code to your needs: * * * * On panTompkins.h: * * - typedef int dataType; * * Change it from 'int' to whatever format your data is (float, unsigned int etc)* * * * On both panTompkins.h and panTompkins.c: * * - void init(char file_in[], char file_out[]); * * This function is meant to do any kind of initial setup, such as starting a * * serial connection with an ECG sensor. Change its parameters to whatever info * * you need and its content. The test version included here loads 2 plain text * * files: an input file, with the signal as a list of integer numbers in ASCII * * format and an output file where either 0 or 1 will be written for each sample,* * whether a peak was detected or not. * * * * On panTompkins.c: * * - #define WINDOWSIZE * * Defines the size of the integration window. The original authors suggest on * * their 1985 paper a 150ms window. * * * * - #define FS * * Defines the sampling frequency. * * * * - #define NOSAMPLE * * A value to indicate you don't have any more samples to read. Choose a value * * which a sample couldn't possibly have (e.g.: a negative value if your A/D con-* * verter only works with positive integers). * * * * - #define BUFFSIZE * * The size of the signal buffers. It should fit at least 1.66 times an RR-inter-* * val. Heart beats should be between 60 and 80 BPS for humans. So, considering * * 1.66 times 1 second should be safe. * * * * - #define DELAY 22 * * The delay introduced to the output signal. The first DELAY samples will be ig-* * nored, as the filters add a delay to the output signal, causing a mismatch * * between the input and output signals. It's easier to compare them this way. * * If you need them both to have the same amount of samples, set this to 0. If * * you're working with different filters and/or sampling rates, you might need to* * adjust this value. * * * * - #include <stdio.h> * * The file, as it is, both gets its inputs and sends its outputs to files. It * * works on both Windows and Linux. If your source isn't a file, and/or your sys-* * tem doesn't have the <stdio.h> header, remove it. * * Include any other headers you might need to make your implementation work, * * such as hardware libraries provided by your microcontroller manufacturer. * * * * - The input() function * * Change it to get the next sample from your source (a file, a serial device etc* * previously set up in your init() function. Return the sample value or NOSAMPLE* * if there are no more available samples. * * * * - The output() function * * Change it to output whatever you see fit, however you see fit: an RR-interval * * (which can be sent as a parameter to your function using the RR arrays), the * * index of sample or timestamp which caused a R peak, whether a sample was a R * * peak or not etc, and it can be written on a file, displayed on screen, blink a* * LED etc. * * * * - The panTompkins() function * * This function is almost entirely ANSI C, which means it should work as is on * * most platforms. The only lines you really have to change are the fclose() ones* * at the very end, which are only here to allow testing of the code on systems * * such as Windows and Linux for PC. You may wish to create extra variables or * * increase the buffers' size as you see fit to extract different kinds of infor-* * mation to output, or fine tune the detection as you see fit - such as adding * * different conditions for verification, initializing self-updating variables * * with empirically-obtained values or changing the filters. * *-------------------------------------------------------------------------------* */ #ifndef PAN_TOMPKINS #define PAN_TOMPKINS //#define WINDOWSIZE 20 // Integrator window size, in samples. The article recommends 150ms. So, FS*0.15. // However, you should check empirically if the waveform looks ok. //#define NOSAMPLE -32000 // An indicator that there are no more samples to read. Use an impossible value for a sample. #define FS 200 // Sampling frequency. #define BUFFSIZE 136 // The size of the buffers (in samples). Must fit more than 1.66 times an RR interval, which // typically could be around 1 second. //#define DELAY 22 // Delay introduced by the filters. Filter only output samples after this one. // Set to 0 if you want to keep the delay. Fixing the delay results in DELAY less samples // in the final end result. #include "panTompkins.h" //#include <stdio.h> // Remove if not using the standard file functions. #define BUFF_SIZE 136 //FILE *fin, *fout; // Remove them if not using files and <stdio.h>. extern float signal[BUFF_SIZE]; //store signal segment to be filtered extern float bp_signal[BUFF_SIZE]; extern float derivative[BUFF_SIZE]; extern float squared[BUFF_SIZE]; extern float integral[BUFF_SIZE]; extern bool outputSignal[BUFF_SIZE]; // extern int rr1[8], rr2[8], rravg1, rravg2, rrlow , rrhigh , rrmiss ; extern long unsigned int i, j, sample , current , lastQRS , lastSlope , currentSlope ; extern float peak_i , peak_f , threshold_i1 , threshold_i2 , threshold_f1 , threshold_f2 , spk_i , spk_f , npk_i , npk_f ; extern bool qrs, regular, prevRegular; // Variables for preprocessing extern float y ; //filtered sample to be displayed extern float prev_y ; // keep track of previous output to calculate derivative extern float sq_y; extern float movmean; // result of moving window integration //extern Serial pc; /* Use this function for any kind of setup you need before getting samples. This is a good place to open a file, initialize your hardware and/or open a serial connection. Remember to update its parameters on the panTompkins.h file as well. */ void panTompkinsInit() { //fin = fopen(file_in, "r"); //fout = fopen(file_out, "w+"); // Initializing the RR averages for (i = 0; i < 8; i++) { rr1[i] = 0; rr2[i] = 0; } //pc.printf("panTompkinsInit() ..."); } /* Use this function to read and return the next sample (from file, serial, A/D converter etc) and put it in a suitable, numeric format. Return the sample, or NOSAMPLE if there are no more samples. */ /* float input() { int num = NOSAMPLE; if (!feof(fin)) fscanf(fin, "%d", &num); return num; }*/ /* Use this function to output the information you see fit (last RR-interval, sample index which triggered a peak detection, whether each sample was a R peak (1) or not (0) etc), in whatever way you see fit (write on screen, write on file, blink a LED, call other functions to do other kinds of processing, such as feature extraction etc). Change its parameters to receive the necessary information to output. */ /* This is the actual QRS-detecting function. It's a loop that constantly calls the input and output functions and updates the thresholds and averages until there are no more samples. More details both above and in shorter comments below. */ void panTompkinsIter() { //pc.printf("PanTompkinsIter"); // The main loop where everything proposed in the paper happens. Ends when there are no more signal samples. qrs = false; // If the current signal is above one of the thresholds (integral or filtered signal), it's a peak candidate. if (integral[current] >= threshold_i1 || bp_signal[current] >= threshold_f1) { peak_i = integral[current]; peak_f = bp_signal[current]; } // If both the integral and the signal are above their thresholds, they're probably signal peaks. if ((integral[current] >= threshold_i1) && (bp_signal[current] >= threshold_f1)) { // There's a 200ms latency. If the new peak respects this condition, we can keep testing. if (sample > lastQRS + FS/5) { // If it respects the 200ms latency, but it doesn't respect the 360ms latency, we check the slope. if (sample <= lastQRS + (long unsigned int)(0.36*FS)) { // 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. currentSlope = 0; for (j = current - 10; j <= current; j++) if (squared[j] > currentSlope) currentSlope = squared[j]; if (currentSlope <= (float)(lastSlope/2)) { qrs = false; } 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; lastSlope = currentSlope; qrs = true; } } // If it was above both thresholds and respects both latency periods, it certainly is a R peak. else { currentSlope = 0; for (j = current - 10; j <= current; j++) if (squared[j] > currentSlope) currentSlope = squared[j]; 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; lastSlope = currentSlope; qrs = true; } } // 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[current]; 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 = bp_signal[current]; 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]);*/ return; } } // If a R-peak 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 (i = 0; i < 7; i++) { rr1[i] = rr1[i+1]; rravg1 += rr1[i]; } rr1[7] = sample - lastQRS; lastQRS = 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 (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; } prevRegular = regular; if (rravg1 == rravg2) { regular = true; } // If the beat had been normal but turned odd, change the thresholds. else { regular = false; if (prevRegular) { 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 { // 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 - lastQRS > (long unsigned int)rrmiss) && (sample > lastQRS + FS/5)) { for (i = current - (sample - lastQRS) + FS/5; i < (long unsigned int)current; i++) { if ( (integral[i] > threshold_i2) && (bp_signal[i] > threshold_f2)) { currentSlope = 0; for (j = i - 10; j <= i; j++) if (squared[j] > currentSlope) currentSlope = squared[j]; if ((currentSlope < (float)(lastSlope/2)) && (i + sample) < lastQRS + 0.36*lastQRS) { qrs = false; } else { peak_i = integral[i]; peak_f = bp_signal[i]; 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; lastSlope = currentSlope; 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 (j = 0; j < 7; j++) { rr1[j] = rr1[j+1]; rravg1 += rr1[j]; } rr1[7] = sample - (current - i) - lastQRS; qrs = true; lastQRS = sample - (current - i); rravg1 += rr1[7]; rravg1 *= 0.125; //RR Average 2 if ( (rr1[7] >= rrlow) && (rr1[7] <= rrhigh) ) { rravg2 = 0; for (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; } prevRegular = regular; if (rravg1 == rravg2) { regular = true; } else { regular = false; if (prevRegular) { threshold_i1 /= 2; threshold_f1 /= 2; } } break; } } } if (qrs) { outputSignal[current] = false; outputSignal[i] = true; /*if (sample > DELAY + BUFFSIZE) output(outputSignal[0]);*/ return; } } // 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[current] >= threshold_i1) || (bp_signal[current] >= threshold_f1)) { peak_i = integral[current]; 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 = bp_signal[current]; 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; } } } // The current implementation outputs '0' for every sample where no peak was detected, // and '1' for every sample where a peak was detected. It should be changed to fit // the desired application. // The 'if' accounts for the delay introduced by the filters: we only start outputting after the delay. // However, it updates a few samples back from the buffer. The reason is that if we update the detection // for the current sample, we might miss a peak that could've been found later by backsearching using // lighter thresholds. The final waveform output does match the original signal, though. outputSignal[current] = qrs; } /* if (sample > DELAY + BUFFSIZE) output(outputSignal[0]); */ //} while(1);//while (signal[current] != NOSAMPLE); // Output the last remaining samples on the buffer //for (i = 1; i < BUFFSIZE; i++) // output(outputSignal[i]); // These last two lines must be deleted if you are not working with files. //fclose(fin); //fclose(fout); #endif