[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

Example 10: A QRS Detector

This program reads a single ECG signal, attempts to detect QRS complexes, and records their locations in an annotation file. The detector algorithm is based on a Pascal program written by W.A.H. Engelse and C. Zeelenberg, “A single scan algorithm for QRS-detection and feature extraction”, Computers in Cardiology 6:37-42 (1979).

 
 1  #include <stdio.h>
 2  #include <wfdb/wfdb.h>
 3  #include <wfdb/ecgcodes.h>
 4  
 5  #define abs(A)  ((A) >= 0 ? (A) : -(A))
 6  
 7  main(argc, argv)
 8  int argc;
 9  char *argv[];
10  {
11      int filter, time=0, slopecrit, sign, maxslope=0, nsig, nslope=0,
12          qtime, maxtime, t0, t1, t2, t3, t4, t5, t6, t7, t8, t9,
13          ms160, ms200, s2, scmax, scmin = 0;
14      WFDB_Anninfo a;
15      WFDB_Annotation annot;
16      WFDB_Sample *v;
17      WFDB_Siginfo *s;
18  
19      if (argc < 2) {
20          fprintf(stderr, "usage: %s record [threshold]\n", argv[0]);
21          exit(1);
22      }
23      a.name = "qrs"; a.stat = WFDB_WRITE;
24      if ((nsig = isigopen(argv[1], NULL, 0)) < 1) exit(2);
25      s = (WFDB_Siginfo *)malloc(nsig * sizeof(WFDB_Siginfo));
26      v = (WFDB_Sample *)malloc(nsig * sizeof(WFDB_Sample));
27      if (s == NULL || v == NULL) {
28          fprintf(stderr, "%s: insufficient memory\n", argv[0]);
29          exit(2);
30      }
31      if (isigopen(argv[1], s, nsig) != nsig) exit(2);
32      if (sampfreq((char *)NULL) < 240. ||
33          sampfreq((char *)NULL) > 260.)
34          setifreq(250.);
35      if (annopen(argv[1], &a, 1) < 0) exit(2);
36      if (argc > 2) scmin = muvadu(0, atoi(argv[2]));
37      if (scmin < 1) scmin = muvadu(0, 1000);
38      slopecrit = scmax = 10 * scmin;
39      ms160 = strtim("0.16"); ms200 = strtim("0.2"); s2 = strtim("2");
40      annot.subtyp = annot.chan = annot.num = 0; annot.aux = NULL;
41      (void)getvec(v);
42      t9 = t8 = t7 = t6 = t5 = t4 = t3 = t2 = t1 = v[0];
43  
44      do {
45          filter = (t0 = v[0]) + 4*t1 + 6*t2 + 4*t3 + t4
46                  - t5         - 4*t6 - 6*t7 - 4*t8 - t9;
47          if (time % s2 == 0) {
48              if (nslope == 0) {
49                  slopecrit -= slopecrit >> 4;
50                  if (slopecrit < scmin) slopecrit = scmin;
51              }
52              else if (nslope >= 5) {
53                  slopecrit += slopecrit >> 4;
54                  if (slopecrit > scmax) slopecrit = scmax;
55              }
56          }
57          if (nslope == 0 && abs(filter) > slopecrit) {
58              nslope = 1; maxtime = ms160;
59              sign = (filter > 0) ? 1 : -1;
60              qtime = time;
61          }
62          if (nslope != 0) {
63              if (filter * sign < -slopecrit) {
64                  sign = -sign;
65                  maxtime = (++nslope > 4) ? ms200 : ms160;
66              }
67              else if (filter * sign > slopecrit &&
68                       abs(filter) > maxslope)
69                  maxslope = abs(filter);
70              if (maxtime-- < 0) {
71                  if (2 <= nslope && nslope <= 4) {
72                      slopecrit += ((maxslope>>2) - slopecrit) >> 3;
73                      if (slopecrit < scmin) slopecrit = scmin;
74                      else if (slopecrit > scmax) slopecrit = scmax;
75                      annot.time = strtim("i") - (time - qtime) - 4;
76                      annot.anntyp = NORMAL; (void)putann(0, &annot);
77                      time = 0;
78                  }
79                  else if (nslope >= 5) {
80                      annot.time = strtim("i") - (time - qtime) - 4;
81                      annot.anntyp = ARFCT; (void)putann(0, &annot);
82                  }
83                  nslope = 0;
84              }
85          }
86          t9 = t8; t8 = t7; t7 = t6; t6 = t5; t5 = t4;
87          t4 = t3; t3 = t2; t2 = t1; t1 = t0; time++;
88      } while (getvec(v) > 0);
89  
90      wfdbquit();
91      exit(0);
92  }

(See http://physionet.org/physiotools/wfdb/examples/example10.c for a copy of this program.)

Notes:

Line 5:

A macro that evaluates to the absolute value of its argument.

Lines 11–12:

The names of these variables match those in the original Pascal program.

Lines 32–35:

Most of this program is independent of sampling frequency, but the filter (lines 45–46) and the threshold are as specified by the authors of the original program for human ECGs sampled at 250 Hz (e.g., the AHA DB). If the sampling frequency of the input record is significantly different, we use setifreq to specify that we want getvec to give us data resampled at 250 Hz. The output annotation file is created in line 35 only after invoking setifreq if necessary.

Lines 36–38:

The threshold is actually a slope criterion (with units of amplitude/time); these lines normalize the threshold with respect to the signal gain. The default value is used unless the user supplies an acceptable alternative. The variables scmin and scmax are lower and upper bounds for the adaptive threshold slopecrit.

Lines 41–42:

Here we read the first sample and copy it into the variables that will be used to store the ten most recent samples.

Lines 45–46:

This FIR filter differentiates and low-pass filters the input signal.

Lines 47–56:

Here we adjust the threshold if more than two seconds have elapsed since a QRS was detected. In line 49, slopecrit is set to 15/16 of its previous value if no slopes have been found; in line 53, it is set to 17/16 of its previous value if 5 or more slopes were found (suggesting the presence of noise).

Lines 57–61:

If the condition in line 48 is satisfied, we may have found the beginning of a QRS complex. We record that a slope has been found, set the timer maxtime to 160 msec, and save the sign of the slope and the current time relative to the previous beat.

Lines 62–85:

This code is executed once we have found a slope. Each time the filter output crosses the threshold, we record another slope and begin looking for a threshold crossing of the opposite sign (lines 63–66), which must occur within a specified time. We record the maximum absolute value of the filter in maxslope (lines 67–69) for eventual use in updating the threshold (lines 72–74). Once a sufficient interval has elapsed following the last threshold crossing (line 70), if there were between 2 and 4 slopes, we have (apparently) found a QRS complex, and the program records a NORMAL annotation (lines 75–76). If there were 5 or more slopes, the program records an artifact annotation (lines 80–81). If only 1 slope was found, it is assumed to represent a baseline shift and no output is produced.

Lines 86–88:

At the end of the loop, the samples are shifted through the tn variables and another sample is read.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]

George B. Moody (george@mit.edu)