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

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      long nsamp, t;
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(3);
36      }
37      if (isigopen(argv[1], s, nsig) != nsig)
38          exit(3);
39      if (isigsettime(strtim(argv[2])) < 0)
40          exit(4);
41      if ((nsamp = strtim(argv[3])) < 1) {
42          fprintf(stderr, "%s: inappropriate value for duration\n",
43                  argv[0]);
44          exit(5);
45      }
46      if (osigopen("16l", s, nsig) != nsig)
47          exit(6);
48  
49      (void)sample(0, 0L);
50      for (t = 0; t < nsamp && sample_valid(); t++) {
51          for (j = 0; j < nsig; j++) {
52              for (i = 0, vv = 0.; i < nc; i++)
53                  if (c[i] != 0.) vv += c[i]*sample(j, t+i);
54              v[j] = (int)vv;
55          }
56          if (putvec(v) < 0) break;
57      }
58  
59      (void)newheader("out");
60      wfdbquit();
61      exit(0);
62  }

(See http://www.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, or the start time is inappropriate, the program exits. See the previous example for details on how isigopen is used.

Lines 41--45:
The duration argument should be a time interval in HH:MM:SS format; strtim converts it to the appropriate number of samples.

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

Line 49:
Here, signal is invoked only for its side effect; assuming that any samples can be read from the specified record, sample(0, 0L) returns a valid sample, so that the value returned by sample_valid() is true (1).

Lines 50--57:
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 53, 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 52 before we begin, and is converted to an int in line 54 when we are finished. Once an entire sample vector has been filtered, it is written out in line 56. The entire process is repeated up to nsamp times, or until we run out of input samples.

Line 59:
The program creates a header file for record `out', using the signal specifications from record `16l' and the sampling frequency from the input record.

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

George B. Moody (george@mit.edu)