[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