Example 9: A Signal Averager

The following program is considerably more complex than the previous examples in this chapter. It reads an annotation file (for which the annotator name is specified in its first argument, and the record name in the second argument) and selects beats of a specified type to be averaged. The program selects segments of the signals that are within 50 milliseconds of the time of the specified beat annotations, subtracts a baseline estimate from each sample, and calculates an average waveform (by default, the average normal QRS complex).

 1  #include <stdio.h>
 2  #include <wfdb/wfdb.h>
 3  #include <wfdb/ecgmap.h>
 4  
 5  main(argc, argv)
 6  int argc;
 7  char *argv[];
 8  {
 9      int btype, i, j, nbeats = 0, nsig, hwindow, window;
10      long stoptime = 0L, **sum;
11      WFDB_Anninfo a;
12      WFDB_Annotation annot;
13      WFDB_Sample *v, *vb;
14      WFDB_Siginfo *s;
15      void *calloc();
16  
17      if (argc < 3) {
18          fprintf(stderr,
19                  "usage: %s annotator record [beat-type from to]\n",
20                  argv[0]);
21          exit(1);
22      }
23      a.name = argv[1]; a.stat = WFDB_READ;
24      if ((nsig = isigopen(argv[2], NULL, 0)) < 1) exit(2);
25      s = (WFDB_Siginfo *)malloc(nsig * sizeof(WFDB_Siginfo));
26      v = (WFDB_Sample *)malloc(nsig * sizeof(WFDB_Sample));
27      vb = (WFDB_Sample *)malloc(nsig * sizeof(WFDB_Sample));
28      sum = (long **)malloc(nsig * sizeof(long *));
29      if (s == NULL || v == NULL || vb == NULL || sum == NULL) {
30          fprintf(stderr, "%s: insufficient memory\n", argv[0]);
31          exit(2);
32      }
33      if (wfdbinit(argv[2], &a, 1, s, nsig) != nsig) exit(3);
34      hwindow = strtim(".05"); window = 2*hwindow + 1;
35      for (i = 0; i < nsig; i++)
36          if ((sum[i]=(long *)calloc(window,sizeof(long))) == NULL) {
37              fprintf(stderr, "%s: insufficient memory\n", argv[0]);
38              exit(2);
39          }
40      btype = (argc > 3) ? strann(argv[3]) : NORMAL;
41      if (argc > 4) iannsettime(strtim(argv[4]));
42      if (argc > 5) {
43          if ((stoptime = strtim(argv[5])) < 0L)
44              stoptime = -stoptime;
45          if (s[0].nsamp > 0L && stoptime > s[0].nsamp)
46              stoptime = s[0].nsamp;
47      }
48      else stoptime = s[0].nsamp;
49      if (stoptime > 0L) stoptime -= hwindow;
50      while (getann(0, &annot) == 0 && annot.time < hwindow)
51          ;
52      do {
53          if (annot.anntyp != btype) continue;
54          isigsettime(annot.time - hwindow - 1);
55          getvec(vb);
56          for (j = 0; j < window && getvec(v) > 0; j++)
57              for (i = 0; i < nsig; i++)
58                  sum[i][j] += v[i] - vb[i];
59          nbeats++;
60      } while (getann(0, &annot) == 0 &&
61               (stoptime == 0L || annot.time < stoptime));
62      if (nbeats < 1) {
63          fprintf(stderr, "%s: no `%s' beats found\n",
64                  argv[0], annstr(btype));
65          exit(4);
66      }
67      printf("Average of %d `%s' beats:\n", nbeats, annstr(btype));
68      for (j = 0; j < window; j++)
69          for (i = 0; i < nsig; i++)
70              printf("%g%c", (double)sum[i][j]/nbeats,
71                     (i == nsig-1) ? '\n' : '\t');
72      exit(0);
73  }

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

Notes:

Line 34:

The “half-window” is 50 milliseconds wide, and the “window” (the duration of a segment to be entered into the average) is one sample more than twice that amount (i.e., 50 milliseconds to either side of the fiducial point defined by the annotation).

Lines 35–39:

Here we allocate memory for the sum vectors that will be used to store the running totals. See K&R, page 167, for a description of calloc.

Line 40:

If a third argument is present on the command line, it is taken as an annotation code mnemonic for the desired beat type; otherwise, the program will average NORMAL QRS complexes.

Line 41:

If a fourth argument is present on the command line, it is taken as the start time; we arrange for the first annotation to be read by getann to be the first annotation that occurs after the chosen start time.

Lines 42–49:

This code similarly determines when the averaging should stop. Unless no stop time was specified on the command line and the signal length is not defined in the ‘hea’ file for the record, stoptime will have a positive value in line 49, which makes a tiny adjustment so that if a beat annotation occurs within 50 milliseconds of the end of the averaging period, the beat will not be included in the average.

Lines 50-51:

This code addresses the (admittedly unlikely) prospect that the first annotation(s) might occur within the first 50 milliseconds of the record; any such annotations will be excluded from the average.

Lines 52–61:

Here we read annotations (the first is already in annot when we enter the loop, and subsequent annotations are read in line 60); select the desired ones (line 53); skip to the correct spot in the signals (line 54; the sample selected there is the one just before the beginning of the window); read a sample from each signal (line 55) into the vb vector, which will be used as a crude baseline estimate; read window samples from each signal (line 56), subtracting the baseline from each and adding the result into the running totals; update a beat counter (line 59); and check for loop termination conditions (line 61).

Lines 62–71:

This is the output section. If no beats of type btype were found, obviously no average can be printed; note that the message goes to the standard error output, so the user will notice it even if the standard output has been redirected to a file. In the usual case, the averages are printed out as a table, with a column allocated to each signal. Note the cast in line 70 (necessary to preserve precision), and the trick used in line 71 to print a tab after each column but the last in each line.