Example 7: A General-Purpose FIR Filter

This program illustrates the use of sample to obtain random access to signals, a technique that is particularly useful for implementing digital filters. The first argument is the record name, the second and third arguments are the start time and the duration of the segment to be filtered, and the rest of the arguments are finite-impulse-response (FIR) filter coefficients. For example, if this program were compiled into an executable program called ‘filter’, it might be used by

filter 100 5:0 20 .2 .2 .2 .2 .2

which would apply a five-point moving average (rectangular window) filter to 20 seconds of record ‘100’, beginning 5 minutes into the record. The output of the program is readable as record ‘out’, for which a header file is created in the current directory.

 1  #include <stdio.h>
 2  #include <wfdb/wfdb.h>
 3  
 4  main(argc, argv)
 5  int argc;
 6  char *argv[];
 7  {
 8      double *c, one = 1.0, vv, atof();
 9      int i, j, nc = argc - 4, nsig;
10      WFDB_Time nsamp, t, t0, t1;
11      static WFDB_Sample *v;
12      static WFDB_Siginfo *s;
13  
14      if (argc < 4) {
15          fprintf(stderr,
16            "usage: %s record start duration [ coefficients ... ]\n",
17                  argv[0]);
18          exit(1);
19      }
20      if (nc < 1) {
21          nc = 1; c = &one;
22      }
23      else if ((c = (double *)calloc(nc, sizeof(double))) == NULL) {
24          fprintf(stderr, "%s: too many coefficients\n", argv[0]);
25          exit(2);
26      }
27      for (i = 0; i < nc; i++)
28          c[i] = atof(argv[i+4]);
29      if ((nsig = isigopen(argv[1], NULL, 0)) < 1)
30          exit(3);
31      s = (WFDB_Siginfo *)malloc(nsig * sizeof(WFDB_Siginfo));
32      v = (WFDB_Sample *)malloc(nsig * sizeof(WFDB_Sample));
33      if (s == NULL || v == NULL) {
34          fprintf(stderr, "insufficient memory\n");
35          exit(4);
36      }
37      if (isigopen(argv[1], s, nsig) != nsig)
38          exit(5);
39      t0 = strtim(argv[2]);
40      if (t0 < (WFDB_Time)0) t0 = -t0;
41      (void)sample(0, t0);
42      if (!sample_valid()) {
43          fprintf(stderr, "%s: inappropriate value for start time\n",
44                 argv[0]);
45          exit(6);
46      }
47      if ((nsamp = strtim(argv[3])) < 1) {
48          fprintf(stderr, "%s: inappropriate value for duration\n",
49                  argv[0]);
50          exit(7);
51      }
52      t1 = t0 + nsamp;
53      if (osigopen("16l", s, nsig) != nsig)
54          exit(8);
55  
56      for (t = t0; t < t1 && sample_valid(); t++) {
57          for (j = 0; j < nsig; j++) {
58              for (i = 0, vv = 0.; i < nc; i++)
59                  if (c[i] != 0.) vv += c[i]*sample(j, t+i);
60              v[j] = (WFDB_Sample)vv;
61          }
62          if (putvec(v) < 0) break;
63      }
64  
65      (void)newheader("out");
66      wfdbquit();
67      exit(0);
68  }

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

Notes:

Lines 20–22:

If no coefficients are provided on the command line, the program will simply copy the selected segment of the input signals.

Lines 23–28:

If there are more coefficients than there are samples in the circular buffer, or if memory cannot be allocated for the coefficient vector, the program cannot work properly, so it exits with an error message. In lines 27 and 28, the ASCII strings that represent the coefficients are converted to double format and stored in the coefficient vector.

Lines 29–40:

The record name is argv[1], and the start time is argv[2]; if the record can’t be opened, the program exits. See the previous example for details on how isigopen is used. If the user provides an absolute start time (see section timstr and strtim), the negative value returned by strtim is converted to a sample number in line 40.

Lines 41–46:

Here, sample is invoked only for its side effect; if any samples can be read from the specified record beginning at sample number t0, then sample(0, 0L) returns a valid sample, so that the value returned by sample_valid is true (1). If not, the program exits.

Lines 47–52:

The duration argument should be a time interval in HH:MM:SS format; strtim converts it to the appropriate number of samples, and t1 is set to the calculated end time in line 52.

Lines 53–54:

The output signals will be written to files in the current directory according to the specifications for record ‘16l’ (see section Piped and Local Records). If we can’t write as many output signals as there are input signals, the program exits.

Lines 56–63:

Here’s where the work is done. The outer loop is executed once per sample vector, the middle loop once per signal, and the inner loop once per coefficient. In line 59, we retrieve an input sample, multiply it by a filter coefficient, and add it to a running sum. The sum (vv) is initialized to zero in line 58 before we begin, and is converted to a WFDB_Sample in line 60 when we are finished. Once an entire output sample vector ia ready, it is written in line 62. The entire process is repeated until we reach input sample number t1, or we run out of input samples.

Line 65:

The program creates a header file for record ‘out’, using the signal specifications from record ‘16l’ and the sampling frequency from the input record.