[R] Designing and incorporating a digital filter

Jonathan Williams jonathan.williams at pharmacology.oxford.ac.uk
Mon Aug 11 18:42:19 CEST 2003


I have a time series of data from an electroencephalogram (EEG).
I wish to filter the data to get rid of 50Hz mains 'hum'. I have
'designed' a combination bandpass and notch filter using a web-
site. The site returns the filter in "ANSI C" source code. It is:-

/* Digital filter designed by mkfilter/mkshape/gencode   A.J. Fisher
   Command line: /www/usr/fisher/helpers/mkfilter -Bu -Bp -o 5 -a
1.0000000000e-02 1.0000000000e-01 -Z 2.0000000000e-01 -l */

#define NZEROS 12
#define NPOLES 12
#define GAIN   1.548068051e+03

static float xv[NZEROS+1], yv[NPOLES+1];

static void filterloop()
  { for (;;)
      { xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; xv[4] =
xv[5]; xv[5] = xv[6]; xv[6] = xv[7]; xv[7] = xv[8]; xv[8] = xv[9]; xv[9] =
xv[10]; xv[10] = xv[11]; xv[11] = xv[12];
        xv[12] = `next input value' / GAIN;
        yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; yv[4] =
yv[5]; yv[5] = yv[6]; yv[6] = yv[7]; yv[7] = yv[8]; yv[8] = yv[9]; yv[9] =
yv[10]; yv[10] = yv[11]; yv[11] = yv[12];
        yv[12] =   (xv[12] - xv[0]) +   0.6180339888 * (xv[1] - xv[11]) + 4
* (xv[2] - xv[10])
                     +   3.0901699437 * (xv[9] - xv[3]) + 5 * (xv[8] -
xv[4]) +   6.1803398875 * (xv[5] - xv[7])
                     - 1.77636e-15 * xv[6]
                     + ( -0.0000000000 * yv[0]) + ( -0.0000000000 * yv[1])
                     + ( -0.1555326685 * yv[2]) + (  1.8017745109 * yv[3])
                     + ( -9.4758632417 * yv[4]) + ( 29.7961492230 * yv[5])
                     + (-62.0371858570 * yv[6]) + ( 89.3642283590 * yv[7])
                     + (-90.1867413610 * yv[8]) + ( 62.9470473210 * yv[9])
                     + (-29.0651076170 * yv[10]) + (  8.0112312881 *
yv[11]);
        `next output value' = yv[12];
      }
  }

I'd like to know how to incorporate this into my R code (or, better, I'd
like to know how to design bandpass filters directly in R - the help file
for the filter() function does not detail how to do this)

Thank you
Jonathan Williams

Jonathan Williams
OPTIMA
Radcliffe Infirmary
Woodstock Road
OXFORD OX2 6HE
Tel +1865 (2)24356




More information about the R-help mailing list