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 <ecg/db.h>
 3  #include <ecg/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, nslope=0,
12          qtime, maxtime, t0, t1, t2, t3, t4, t5, t6, t7, t8, t9,
13          ms160, ms200, s2, scmax, scmin = 0, v[DB_MAXSIG];
14      DB_Anninfo a;
15      DB_Annotation annot;
16      static DB_Siginfo s[DB_MAXSIG];
17  
18      if (argc < 2) {
19          fprintf(stderr, "usage: %s record [threshold]\n", argv[0]);
20          exit(1);
21      }
22      a.name = argv[0]; a.stat = WRITE;
23      if (dbinit(argv[1], &a, 1, s, DB_MAXSIG) < 1) exit(2);
24      if (sampfreq(NULL) != 250.)
25          fprintf(stderr, "warning: %s is designed for 250 Hz input\n",
26                  argv[0]);
27      if (argc > 2) scmin = muvadu(0, atoi(argv[2]));
28      if (scmin < 1) scmin = muvadu(0, 1000);
29      slopecrit = scmax = 10 * scmin;
30      ms160 = strtim("0.16"); ms200 = strtim("0.2"); s2 = strtim("2");
31      annot.subtyp = annot.chan = annot.num = 0; annot.aux = NULL;
32      (void)getvec(v);
33      t9 = t8 = t7 = t6 = t5 = t4 = t3 = t2 = t1 = v[0];
34  
35      do {
36          filter = (t0 = v[0]) + 4*t1 + 6*t2 + 4*t3 + t4
37                  - t5         - 4*t6 - 6*t7 - 4*t8 - t9;
38          if (time % s2 == 0) {
39              if (nslope == 0) {
40                  slopecrit -= slopecrit >> 4;
41                  if (slopecrit < scmin) slopecrit = scmin;
42              }
43              else if (nslope >= 5) {
44                  slopecrit += slopecrit >> 4;
45                  if (slopecrit > scmax) slopecrit = scmax;
46              }
47          }
48          if (nslope == 0 && abs(filter) > slopecrit) {
49              nslope = 1; maxtime = ms160;
50              sign = (filter > 0) ? 1 : -1;
51              qtime = time;
52          }
53          if (nslope != 0) {
54              if (filter * sign < -slopecrit) {
55                  sign = -sign;
56                  maxtime = (++nslope > 4) ? ms200 : ms160;
57              }
58              else if (filter * sign > slopecrit &&
59                       abs(filter) > maxslope)
60                  maxslope = abs(filter);
61              if (maxtime-- < 0) {
62                  if (2 <= nslope && nslope <= 4) {
63                      slopecrit += ((maxslope>>2) - slopecrit) >> 3;
64                      if (slopecrit < scmin) slopecrit = scmin;
65                      else if (slopecrit > scmax) slopecrit = scmax;
66                      annot.time = strtim("i") - (time - qtime) - 4;
67                      annot.anntyp = NORMAL; (void)putann(0, &annot);
68                      time = 0;
69                  }
70                  else if (nslope >= 5) {
71                      annot.time = strtim("i") - (time - qtime) - 4;
72                      annot.anntyp = ARFCT; (void)putann(0, &annot);
73                  }
74                  nslope = 0;
75              }
76          }
77          t9 = t8; t8 = t7; t7 = t6; t6 = t5; t5 = t4;
78          t4 = t3; t3 = t2; t2 = t1; t1 = t0; time++;
79      } while (getvec(v) > 0);
80  
81      dbquit();
82      exit(0);
83  }

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 24--26:

Most of this program is independent of sampling frequency, but the filter (lines 36--37) 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). The program will work for MIT DB records sampled at 360 Hz, but better results can be obtained if the filter is optimized for the correct sampling frequency.

Lines 27--29:

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 32--33:

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

Lines 36--37:

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

Lines 38--47:

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

Lines 48--52:

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 53--76:

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 54--57), which must occur within a specified time. We record the maximum absolute value of the filter in maxslope (lines 58--60) for eventual use in updating the threshold (lines 63--65). Once a sufficient interval has elapsed following the last threshold crossing (line 61), if there were between 2 and 4 slopes, we have (apparently) found a QRS complex, and the program records a NORMAL annotation (lines 66--67). If there were 5 or more slopes, the program records an artifact annotation (lines 71--72). If only 1 slope was found, it is assumed to represent a baseline shift and no output is produced.

Lines 77--79:

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

Go to the first, previous, next, last section, table of contents.