panTompkins.h

00001 /**
00002  * ------------------------------------------------------------------------------*
00003  * File: panTompkins.c                                                           *
00004  *       ANSI-C implementation of Pan-Tompkins real-time QRS detection algorithm *
00005  * Author: Rafael de Moura Moreira <rafaelmmoreira@gmail.com>                    *
00006  * License: MIT License                                                          *
00007  * ------------------------------------------------------------------------------*
00008  * ---------------------------------- HISTORY ---------------------------------- *
00009  *    date   |    author    |                     description                    *
00010  * ----------| -------------| ---------------------------------------------------*
00011  * 2019/04/11| Rafael M. M. | - Fixed moving-window integral.                    *
00012  *           |              | - Fixed how to find the correct sample with the    *
00013  *           |              | last QRS.                                          *
00014  *           |              | - Replaced constant value in code by its #define.  *
00015  *           |              | - Added some casting on comparisons to get rid of  *
00016  *           |              | compiler warnings.                                 *
00017  * 2019/04/15| Rafael M. M. | - Removed delay added to the output by the filters.*
00018  *           |              | - Fixed multiple detection of the same peak.       *
00019  * 2019/04/16| Rafael M. M. | - Added output buffer to correctly output a peak   *
00020  *           |              | found by back searching using the 2nd thresholds.  *
00021  * 2019/04/23| Rafael M. M. | - Improved comparison of slopes.                   *
00022  *           |              | - Fixed formula to obtain the correct sample from  *
00023  *           |              | the buffer on the back search.                     *
00024  * ------------------------------------------------------------------------------*
00025  * MIT License                                                                   *
00026  *                                                                               *
00027  * Copyright (c) 2018 Rafael de Moura Moreira                                    *
00028  *                                                                               *
00029  * Permission is hereby granted, free of charge, to any person obtaining a copy  *
00030  * of this software and associated documentation files (the "Software"), to deal *
00031  * in the Software without restriction, including without limitation the rights  *
00032  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell     *
00033  * copies of the Software, and to permit persons to whom the Software is         *
00034  * furnished to do so, subject to the following conditions:                      *
00035  *                                                                               *
00036  * The above copyright notice and this permission notice shall be included in all*
00037  * copies or substantial portions of the Software.                               *
00038  *                                                                               *
00039  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR    *
00040  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,      *
00041  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE   *
00042  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER        *
00043  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, *
00044  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE *
00045  * SOFTWARE.                                                                     *
00046  *-------------------------------------------------------------------------------*
00047  * Description                                                                   *
00048  *                                                                               *
00049  * The main goal of this implementation is to be easy to port to different opera-*
00050  * ting systems, as well as different processors and microcontrollers, including *
00051  * embedded systems. It can work both online or offline, depending on whether all*
00052  * the samples are available or not - it can be adjusted on the input function.  *
00053  *                                                                               *
00054  * The main function, panTompkings(), calls input() to get the next sample and   *
00055  * store it in a buffer. Then it runs through a chain of filters: DC block, low  *
00056  * pass @ 15 Hz and high pass @ 5Hz. The filtered signal goes both through a de- *
00057  * rivative filter, which output is then squared, and through a windowed-integra-*
00058  * tor.                                                                          *
00059  *                                                                               *
00060  * For a signal peak to be recognized as a fiducial point, its correspondent va- *
00061  * lue on both the filtered signal and the integrator must be above a certain    *
00062  * threshold. Additionally, there are time-restraints to prevent a T-wave from   *
00063  * being mistakenly identified as an R-peak: a hard 200ms restrain (a new peak   *
00064  * 200ms from the previous one is, necessarily, a T-wave) and a soft 360ms res-  *
00065  * train (the peak's squared slope must also be very high to be considered as a  *
00066  * real peak).                                                                   *
00067  *                                                                               *
00068  * When a peak candidate is discarded, its value is used to update the noise     *
00069  * thresholds - which are also used to estimate the peak thresholds.             *
00070  *                                                                               *
00071  * Two buffers keep 8 RR-intervals to calculate RR-averages: one of them keeps   *
00072  * the last 8 RR-intervals, while the other keeps only the RR-intervals that res-*
00073  * pect certain restrictions. If both averages are equal, the heart pace is con- *
00074  * sidered normal. If the heart rate isn't normal, the thresholds change to make *
00075  * it easier to detect possible weaker peaks. If no peak is detected for a long  *
00076  * period of time, the thresholds also change and the last discarded peak candi- *
00077  * date is reconsidered.                                                         *
00078  *-------------------------------------------------------------------------------*
00079  * Instructions                                                                  *
00080  *                                                                               *
00081  * Here's what you should change to adjust the code to your needs:               *
00082  *                                                                               *
00083  * On panTompkins.h:                                                             *
00084  * - typedef int dataType;                                                       *
00085  * Change it from 'int' to whatever format your data is (float, unsigned int etc)*
00086  *                                                                               *
00087  * On both panTompkins.h and panTompkins.c:                                      *
00088  * - void init(char file_in[], char file_out[]);                                 *
00089  * This function is meant to do any kind of initial setup, such as starting a    *
00090  * serial connection with an ECG sensor. Change its parameters to whatever info  *
00091  * you need and its content. The test version included here loads 2 plain text   *
00092  * files: an input file, with the signal as a list of integer numbers in ASCII   *
00093  * format and an output file where either 0 or 1 will be written for each sample,*
00094  * whether a peak was detected or not.                                           *
00095  *                                                                               *
00096  * On panTompkins.c:                                                             *
00097  * - #define WINDOWSIZE                                                          *
00098  * Defines the size of the integration window. The original authors suggest on   *
00099  * their 1985 paper a 150ms window.                                              *
00100  *                                                                               *
00101  * - #define FS                                                                  *
00102  * Defines the sampling frequency.                                               *
00103  *                                                                               *
00104  * - #define NOSAMPLE                                                            *
00105  * A value to indicate you don't have any more samples to read. Choose a value   *
00106  * which a sample couldn't possibly have (e.g.: a negative value if your A/D con-*
00107  * verter only works with positive integers).                                    *
00108  *                                                                               *
00109  * - #define BUFFSIZE                                                            *
00110  * The size of the signal buffers. It should fit at least 1.66 times an RR-inter-*
00111  * val. Heart beats should be between 60 and 80 BPS for humans. So, considering  *
00112  * 1.66 times 1 second should be safe.                                           *
00113  *                                                                               *
00114  * - #define DELAY 22                                                            *
00115  * The delay introduced to the output signal. The first DELAY samples will be ig-*
00116  * nored, as the filters add a delay to the output signal, causing a mismatch    *
00117  * between the input and output signals. It's easier to compare them this way.   *
00118  * If you need them both to have the same amount of samples, set this to 0. If   *
00119  * you're working with different filters and/or sampling rates, you might need to*
00120  * adjust this value.                                                            *
00121  *                                                                               *
00122  * - #include <stdio.h>                                                          *
00123  * The file, as it is, both gets its inputs and sends its outputs to files. It   *
00124  * works on both Windows and Linux. If your source isn't a file, and/or your sys-*
00125  * tem doesn't have the <stdio.h> header, remove it.                             *
00126  * Include any other headers you might need to make your implementation work,    *
00127  * such as hardware libraries provided by your microcontroller manufacturer.     *
00128  *                                                                               *
00129  * - The input() function                                                        *
00130  * Change it to get the next sample from your source (a file, a serial device etc*
00131  * previously set up in your init() function. Return the sample value or NOSAMPLE*
00132  * if there are no more available samples.                                       *
00133  *                                                                               *
00134  * - The output() function                                                       *
00135  * Change it to output whatever you see fit, however you see fit: an RR-interval *
00136  * (which can be sent as a parameter to your function using the RR arrays), the  *
00137  * index of sample or timestamp which caused a R peak, whether a sample was a R  *
00138  * peak or not etc, and it can be written on a file, displayed on screen, blink a*
00139  * LED etc.                                                                      *
00140  *                                                                               *
00141  * - The panTompkins() function                                                  *
00142  * This function is almost entirely ANSI C, which means it should work as is on  *
00143  * most platforms. The only lines you really have to change are the fclose() ones*
00144  * at the very end, which are only here to allow testing of the code on systems  *
00145  * such as Windows and Linux for PC. You may wish to create extra variables or   *
00146  * increase the buffers' size as you see fit to extract different kinds of infor-*
00147  * mation to output, or fine tune the detection as you see fit - such as adding  *
00148  * different conditions for verification, initializing self-updating variables   *
00149  * with empirically-obtained values or changing the filters.                     *
00150  *-------------------------------------------------------------------------------*
00151  */
00152  
00153 #ifndef PAN_TOMPKINS
00154 #define PAN_TOMPKINS
00155 
00156 //#define WINDOWSIZE 20   // Integrator window size, in samples. The article recommends 150ms. So, FS*0.15.
00157                         // However, you should check empirically if the waveform looks ok.
00158 //#define NOSAMPLE -32000 // An indicator that there are no more samples to read. Use an impossible value for a sample.
00159 #define FS 200          // Sampling frequency.
00160 #define BUFFSIZE 136    // The size of the buffers (in samples). Must fit more than 1.66 times an RR interval, which
00161                         // typically could be around 1 second.
00162 
00163 //#define DELAY 22        // Delay introduced by the filters. Filter only output samples after this one.
00164                         // Set to 0 if you want to keep the delay. Fixing the delay results in DELAY less samples
00165                         // in the final end result.
00166 
00167 #include "panTompkins.h"
00168 
00169 //#include <stdio.h>      // Remove if not using the standard file functions.
00170 #define BUFF_SIZE 136 
00171 
00172 //FILE *fin, *fout;       // Remove them if not using files and <stdio.h>.
00173 extern float signal[BUFF_SIZE]; //store signal segment to be filtered
00174 extern float bp_signal[BUFF_SIZE];
00175 extern float derivative[BUFF_SIZE];
00176 extern float squared[BUFF_SIZE];
00177 extern float integral[BUFF_SIZE];
00178 extern bool outputSignal[BUFF_SIZE];
00179 
00180 //
00181 extern int rr1[8], rr2[8], rravg1, rravg2, rrlow , rrhigh , rrmiss ;
00182 extern long unsigned int i, j, sample , current , lastQRS , lastSlope , currentSlope ;
00183 extern float peak_i , peak_f , threshold_i1 , threshold_i2 , threshold_f1 , threshold_f2 , spk_i , spk_f , npk_i , npk_f ;
00184 extern bool qrs, regular, prevRegular;
00185 
00186 // Variables for preprocessing
00187 extern float y ; //filtered sample to be displayed
00188 extern float prev_y ; // keep track of previous output to calculate derivative
00189 extern float sq_y;
00190 extern float movmean; // result of moving window integration
00191 //extern Serial pc;
00192 
00193 /*
00194     Use this function for any kind of setup you need before getting samples.
00195     This is a good place to open a file, initialize your hardware and/or open
00196     a serial connection.
00197     Remember to update its parameters on the panTompkins.h file as well.
00198 */
00199 void panTompkinsInit()
00200 {
00201     //fin = fopen(file_in, "r");
00202     //fout = fopen(file_out, "w+");
00203       // Initializing the RR averages
00204     for (i = 0; i < 8; i++)
00205     {
00206         rr1[i] = 0;
00207         rr2[i] = 0;
00208     }
00209     //pc.printf("panTompkinsInit() ...");
00210 }
00211 
00212 /*
00213     Use this function to read and return the next sample (from file, serial,
00214     A/D converter etc) and put it in a suitable, numeric format. Return the
00215     sample, or NOSAMPLE if there are no more samples.
00216 */
00217 /*
00218 float input()
00219 {
00220     int num = NOSAMPLE;
00221     if (!feof(fin))
00222         fscanf(fin, "%d", &num);
00223 
00224     return num;
00225 }*/
00226 
00227 /*
00228     Use this function to output the information you see fit (last RR-interval,
00229     sample index which triggered a peak detection, whether each sample was a R
00230     peak (1) or not (0) etc), in whatever way you see fit (write on screen, write
00231     on file, blink a LED, call other functions to do other kinds of processing,
00232     such as feature extraction etc). Change its parameters to receive the necessary
00233     information to output.
00234 */
00235 
00236 /*
00237     This is the actual QRS-detecting function. It's a loop that constantly calls the input and output functions
00238     and updates the thresholds and averages until there are no more samples. More details both above and in
00239     shorter comments below.
00240 */
00241 
00242 
00243     
00244 
00245   
00246 void panTompkinsIter()
00247 {
00248     //pc.printf("PanTompkinsIter");
00249     // The main loop where everything proposed in the paper happens. Ends when there are no more signal samples.
00250 
00251         qrs = false;
00252 
00253         // If the current signal is above one of the thresholds (integral or filtered signal), it's a peak candidate.
00254         if (integral[current] >= threshold_i1 || bp_signal[current] >= threshold_f1)
00255         {
00256             peak_i = integral[current];
00257             peak_f = bp_signal[current];
00258         }
00259 
00260         // If both the integral and the signal are above their thresholds, they're probably signal peaks.
00261         if ((integral[current] >= threshold_i1) && (bp_signal[current] >= threshold_f1))
00262         {
00263             // There's a 200ms latency. If the new peak respects this condition, we can keep testing.
00264             if (sample > lastQRS + FS/5)
00265             {
00266                 // If it respects the 200ms latency, but it doesn't respect the 360ms latency, we check the slope.
00267                 if (sample <= lastQRS + (long unsigned int)(0.36*FS))
00268                 {
00269                     // The squared slope is "M" shaped. So we have to check nearby samples to make sure we're really looking
00270                     // at its peak value, rather than a low one.
00271                     currentSlope = 0;
00272                     for (j = current - 10; j <= current; j++)
00273                         if (squared[j] > currentSlope)
00274                             currentSlope = squared[j];
00275 
00276                     if (currentSlope <= (float)(lastSlope/2))
00277                     {
00278                         qrs = false;
00279                     }
00280 
00281                     else
00282                     {
00283                         spk_i = 0.125*peak_i + 0.875*spk_i;
00284                         threshold_i1 = npk_i + 0.25*(spk_i - npk_i);
00285                         threshold_i2 = 0.5*threshold_i1;
00286 
00287                         spk_f = 0.125*peak_f + 0.875*spk_f;
00288                         threshold_f1 = npk_f + 0.25*(spk_f - npk_f);
00289                         threshold_f2 = 0.5*threshold_f1;
00290 
00291                         lastSlope = currentSlope;
00292                         qrs = true;
00293                     }
00294                 }
00295                 // If it was above both thresholds and respects both latency periods, it certainly is a R peak.
00296                 else
00297                 {
00298                     currentSlope = 0;
00299                     for (j = current - 10; j <= current; j++)
00300                         if (squared[j] > currentSlope)
00301                             currentSlope = squared[j];
00302 
00303                     spk_i = 0.125*peak_i + 0.875*spk_i;
00304                     threshold_i1 = npk_i + 0.25*(spk_i - npk_i);
00305                     threshold_i2 = 0.5*threshold_i1;
00306 
00307                     spk_f = 0.125*peak_f + 0.875*spk_f;
00308                     threshold_f1 = npk_f + 0.25*(spk_f - npk_f);
00309                     threshold_f2 = 0.5*threshold_f1;
00310 
00311                     lastSlope = currentSlope;
00312                     qrs = true;
00313                 }
00314             }
00315             // If the new peak doesn't respect the 200ms latency, it's noise. Update thresholds and move on to the next sample.
00316             else
00317             {
00318                 peak_i = integral[current];
00319                 npk_i = 0.125*peak_i + 0.875*npk_i;
00320                 threshold_i1 = npk_i + 0.25*(spk_i - npk_i);
00321                 threshold_i2 = 0.5*threshold_i1;
00322                 peak_f = bp_signal[current];
00323                 npk_f = 0.125*peak_f + 0.875*npk_f;
00324                 threshold_f1 = npk_f + 0.25*(spk_f - npk_f);
00325                 threshold_f2 = 0.5*threshold_f1;
00326                 qrs = false;
00327                 outputSignal[current] = qrs;
00328                 /*if (sample > DELAY + BUFFSIZE)
00329                     output(outputSignal[0]);*/
00330                 return;
00331             }
00332 
00333         }
00334 
00335         // If a R-peak was detected, the RR-averages must be updated.
00336         if (qrs)
00337         {
00338             // Add the newest RR-interval to the buffer and get the new average.
00339             rravg1 = 0;
00340             for (i = 0; i < 7; i++)
00341             {
00342                 rr1[i] = rr1[i+1];
00343                 rravg1 += rr1[i];
00344             }
00345             rr1[7] = sample - lastQRS;
00346             lastQRS = sample;
00347             rravg1 += rr1[7];
00348             rravg1 *= 0.125;
00349 
00350             // If the newly-discovered RR-average is normal, add it to the "normal" buffer and get the new "normal" average.
00351             // Update the "normal" beat parameters.
00352             if ( (rr1[7] >= rrlow) && (rr1[7] <= rrhigh) )
00353             {
00354                 rravg2 = 0;
00355                 for (i = 0; i < 7; i++)
00356                 {
00357                     rr2[i] = rr2[i+1];
00358                     rravg2 += rr2[i];
00359                 }
00360                 rr2[7] = rr1[7];
00361                 rravg2 += rr2[7];
00362                 rravg2 *= 0.125;
00363                 rrlow = 0.92*rravg2;
00364                 rrhigh = 1.16*rravg2;
00365                 rrmiss = 1.66*rravg2;
00366             }
00367 
00368             prevRegular = regular;
00369             if (rravg1 == rravg2)
00370             {
00371                 regular = true;
00372             }
00373             // If the beat had been normal but turned odd, change the thresholds.
00374             else
00375             {
00376                 regular = false;
00377                 if (prevRegular)
00378                 {
00379                     threshold_i1 /= 2;
00380                     threshold_f1 /= 2;
00381                 }
00382             }
00383         }
00384         // If no R-peak was detected, it's important to check how long it's been since the last detection.
00385         else
00386         {
00387             // If no R-peak was detected for too long, use the lighter thresholds and do a back search.
00388             // However, the back search must respect the 200ms limit and the 360ms one (check the slope).
00389             if ((sample - lastQRS > (long unsigned int)rrmiss) && (sample > lastQRS + FS/5))
00390             {
00391                 for (i = current - (sample - lastQRS) + FS/5; i < (long unsigned int)current; i++)
00392                 {
00393                     if ( (integral[i] > threshold_i2) && (bp_signal[i] > threshold_f2))
00394                     {
00395                         currentSlope = 0;
00396                         for (j = i - 10; j <= i; j++)
00397                             if (squared[j] > currentSlope)
00398                                 currentSlope = squared[j];
00399 
00400                         if ((currentSlope < (float)(lastSlope/2)) && (i + sample) < lastQRS + 0.36*lastQRS)
00401                         {
00402                             qrs = false;
00403                         }
00404                         else
00405                         {
00406                             peak_i = integral[i];
00407                             peak_f = bp_signal[i];
00408                             spk_i = 0.25*peak_i+ 0.75*spk_i;
00409                             spk_f = 0.25*peak_f + 0.75*spk_f;
00410                             threshold_i1 = npk_i + 0.25*(spk_i - npk_i);
00411                             threshold_i2 = 0.5*threshold_i1;
00412                             lastSlope = currentSlope;
00413                             threshold_f1 = npk_f + 0.25*(spk_f - npk_f);
00414                             threshold_f2 = 0.5*threshold_f1;
00415                             // If a signal peak was detected on the back search, the RR attributes must be updated.
00416                             // This is the same thing done when a peak is detected on the first try.
00417                             //RR Average 1
00418                             rravg1 = 0;
00419                             for (j = 0; j < 7; j++)
00420                             {
00421                                 rr1[j] = rr1[j+1];
00422                                 rravg1 += rr1[j];
00423                             }
00424                             rr1[7] = sample - (current - i) - lastQRS;
00425                             qrs = true;
00426                             lastQRS = sample - (current - i);
00427                             rravg1 += rr1[7];
00428                             rravg1 *= 0.125;
00429 
00430                             //RR Average 2
00431                             if ( (rr1[7] >= rrlow) && (rr1[7] <= rrhigh) )
00432                             {
00433                                 rravg2 = 0;
00434                                 for (i = 0; i < 7; i++)
00435                                 {
00436                                     rr2[i] = rr2[i+1];
00437                                     rravg2 += rr2[i];
00438                                 }
00439                                 rr2[7] = rr1[7];
00440                                 rravg2 += rr2[7];
00441                                 rravg2 *= 0.125;
00442                                 rrlow = 0.92*rravg2;
00443                                 rrhigh = 1.16*rravg2;
00444                                 rrmiss = 1.66*rravg2;
00445                             }
00446 
00447                             prevRegular = regular;
00448                             if (rravg1 == rravg2)
00449                             {
00450                                 regular = true;
00451                             }
00452                             else
00453                             {
00454                                 regular = false;
00455                                 if (prevRegular)
00456                                 {
00457                                     threshold_i1 /= 2;
00458                                     threshold_f1 /= 2;
00459                                 }
00460                             }
00461 
00462                             break;
00463                         }
00464                     }
00465                 }
00466 
00467                 if (qrs)
00468                 {
00469                     outputSignal[current] = false;
00470                     outputSignal[i] = true;
00471                     /*if (sample > DELAY + BUFFSIZE)
00472                         output(outputSignal[0]);*/
00473                     return;
00474                 }
00475             }
00476 
00477             // Definitely no signal peak was detected.
00478             if (!qrs)
00479             {
00480                 // If some kind of peak had been detected, then it's certainly a noise peak. Thresholds must be updated accordinly.
00481                 if ((integral[current] >= threshold_i1) || (bp_signal[current] >= threshold_f1))
00482                 {
00483                     peak_i = integral[current];
00484                     npk_i = 0.125*peak_i + 0.875*npk_i;
00485                     threshold_i1 = npk_i + 0.25*(spk_i - npk_i);
00486                     threshold_i2 = 0.5*threshold_i1;
00487                     peak_f = bp_signal[current];
00488                     npk_f = 0.125*peak_f + 0.875*npk_f;
00489                     threshold_f1 = npk_f + 0.25*(spk_f - npk_f);
00490                     threshold_f2 = 0.5*threshold_f1;
00491                 }
00492             }
00493         }
00494         // The current implementation outputs '0' for every sample where no peak was detected,
00495         // and '1' for every sample where a peak was detected. It should be changed to fit
00496         // the desired application.
00497         // The 'if' accounts for the delay introduced by the filters: we only start outputting after the delay.
00498         // However, it updates a few samples back from the buffer. The reason is that if we update the detection
00499         // for the current sample, we might miss a peak that could've been found later by backsearching using
00500         // lighter thresholds. The final waveform output does match the original signal, though.
00501         outputSignal[current] = qrs;
00502     
00503     
00504 }
00505         /*
00506         if (sample > DELAY + BUFFSIZE)
00507             output(outputSignal[0]);
00508 
00509             */
00510 //} while(1);//while (signal[current] != NOSAMPLE);
00511 
00512     // Output the last remaining samples on the buffer
00513     //for (i = 1; i < BUFFSIZE; i++)
00514     //    output(outputSignal[i]);
00515 
00516     // These last two lines must be deleted if you are not working with files.
00517     //fclose(fin);
00518     //fclose(fout);
00519 #endif