Example 7: A General-Purpose FIR Filter

This program illustrates a useful technique for obtaining something close to 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 <ecg/db.h>
 3  #define BUFLN 512
 4  int sample_ok = 1;
 5
 6  DB_Sample sample(s, t)
 7  DB_Signal s;
 8  DB_Time t;
 9  {
10      static DB_Sample sbuf[BUFLN][DB_MAXSIG];
11      static DB_Time tt = -1L;
12
13      if (t <= tt - BUFLN)
14          fprintf(stderr, "sample: buffer too short\n");
15      while (t > tt)
16          if (getvec(sbuf[(++tt)&(BUFLN-1)]) < 0) sample_ok = 0;
17      return (sbuf[t&(BUFLN-1)][s]);
18  }
19
20  main(argc, argv)
21  int argc;
22  char *argv[];
23  {
24      double *c, one = 1.0, vv, atof();
25      int i, j, nc = argc - 4, nsig, v[DB_MAXSIG];
26      long nsamp, t;
27      static DB_Siginfo s[DB_MAXSIG];
28
29      if (argc < 4) {
30          fprintf(stderr,
31            "usage: %s record start duration [ coefficients ... ]\n",
32                  argv[0]);
33          exit(1);
34      }
35      if (nc < 1) {
36          nc = 1; c = &one;
37      }
38      else if (nc >= BUFLN ||
39               (c = (double *)calloc(nc, sizeof(double))) == NULL) {
40          fprintf(stderr, "%s: too many coefficients\n", argv[0]);
41          exit(2);
42      }
43      for (i = 0; i < nc; i++)
44          c[i] = atof(argv[i+4]);
45      if ((nsig = isigopen(argv[1], s, DB_MAXSIG)) < 1)
46          exit(3);
47      if (isigsettime(strtim(argv[2])) < 0)
48          exit(4);
49      if ((nsamp = strtim(argv[3])) < 1) {
50          fprintf(stderr, "%s: inappropriate value for duration\n",
51                  argv[0]);
52          exit(5);
53      }
54      if (osigopen("16l", s, nsig) != nsig)
55          exit(6);
56
57      for (t = 0; t < nsamp && sample_ok; t++) {
58          for (j = 0; j < nsig; j++) {
59              for (i = 0, vv = 0.; i < nc; i++)
60                  if (c[i] != 0.) vv += c[i]*sample(j, t+i);
61              v[j] = (int)vv;
62          }
63          if (putvec(v) < 0) break;
64      }
65
66      (void)newheader("out");
67      dbquit();
68      exit(0);
69  }

Notes:

Line 4:

BUFLN must be a power of 2 (why? see lines 16 and 17), and it should be larger than the length of the filter (i.e., the caller should not look back further than BUFLN-1 samples into the past, relative to the most recent sample that has been read).

Lines 6--18:

This function supplies input samples to the main routine as needed, and frees the main routine of the need to read them in strict time order. The sample function returns the sample from signal s with adjusted sample number t (i.e., relative to the beginning of the segment to be processed), either by retrieving it from a circular buffer of samples recently read, or by reading it using getvec. The ugly-looking getvec argument is simply the next slot in the circular buffer; note that tt is an internal "clock" for sample, which records the adjusted sample number of the most recently read sample. In this program, the test in line 13 is redundant (why?) and might be removed for efficiency's sake.

Lines 35--37:

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

Lines 38--44:

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 43 and 44, the ASCII strings that represent the coefficients are converted to double format and stored in the coefficient vector.

Lines 45--48:

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.

Lines 49--53:

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

Lines 54--55:

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 57--64:

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 60, we retrieve an input sample, multiply it by a filter coefficient, and add it to a running sum. The sum is initialized to zero in line 59 before we begin, and is converted to an int in line 61 when we are finished. Once an entire sample vector has been filtered, it is written out in line 63. The entire process is repeated up to nsamp times, or until we run out of input samples.

Line 66:

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

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