eegkit (version 1.0-5)

eegfilter: Filters EEG Data

Description

Low-pass, high-pass, or band-pass filter EEG data using either a Butterworth filter (default) or a finite impulse response (FIR) filter.

Usage

eegfilter(x, Fs, lower, upper, method = "butter",
          order = 3L, forwardreverse = TRUE, 
          scale = FALSE, plot = FALSE)

Value

Filtered version of input data.

Arguments

x

Vector or matrix (time by channel) of EEG data with n time points.

Fs

Sampling rate of x in Hz.

lower

Lower band in Hz. Smallest frequency to keep.

upper

Upper band in Hz. Largest frequency to keep.

method

Filtering method. Either "butter" for a Butterworth filter or "fir1" for a FIR filter.

order

Order of the filter. See corresponding argument of butter or fir1.

forwardreverse

If TRUE (default), the data are forward and reverse filtered via filtfilt. Otherwise the data are (forward) filtered via filter.

scale

If FALSE (default), the filter is not normalized. Otherwise the magnitude of the center of the first passband is normalized to 1.

plot

If TRUE, the filter is plotted via freqz_plot.

Author

Nathaniel E. Helwig <[e-mail]>

Details

For a low-pass filter, only enter the upper frequency to keep. For a high-pass filter, only enter the lower frequency to keep. For a band-pass filter, enter both the lower and upper frequency bounds.

References

http://en.wikipedia.org/wiki/Butterworth_filter

http://en.wikipedia.org/wiki/Fir_filter

See Also

filter, filtfilt, butter, fir1

Examples

##########   EXAMPLE   ##########

# create data generating signals
n <- 1000                             # 1000 Hz signal
s <- 2                                # 2 seconds of data
t <- seq(0, s, length.out = s * n)    # time vector
s1 <- sin(2*pi*t)                     # 1 Hz sinusoid
s5 <- sin(2*pi*t*5)                   # 5 Hz sinusoid
s10 <- sin(2*pi*t*10)                 # 10 Hz sinusoid
s20 <- sin(2*pi*t*20)                 # 20 Hz sinusoid

# create data
set.seed(1)                           # set random seed
e <- rnorm(s * n, sd = 0.25)          # Gaussian error
mu <- s1 + s5 + s10 + s20             # 1 + 5 + 10 + 20 Hz mean
y <- mu + e                           # data = mean + error

# 4-th order Butterworth filter (2 to 15 Hz band-pass)
yf.but <- eegfilter(y, Fs = n, lower = 2, upper = 15, method = "butter", order = 4)

# 350-th order FIR filter (2 to 15 Hz band-pass)
yf.fir <- eegfilter(y, Fs = n, lower = 2, upper = 15, method = "fir1", order = 350)

# check quality of results
yftrue <- s5 + s10                    # true (filtered) mean signal
mean((yf.but - yftrue)^2)             # mse between yf.but and yftrue
mean((yf.fir - yftrue)^2)             # mse between yf.fir and yftrue

# plot true and estimated filtered signals
plot(t, yftrue, type = "l", lty = 1, lwd = 2, ylim = c(-3, 3))
lines(t, yf.but, col = "blue", lty = 2, lwd = 2)
lines(t, yf.fir, col = "red", lty = 3, lwd = 2)
legend("topright", legend = c("Truth", "Butterworth", "FIR"), 
       lty = 1:3, lwd = 2, col = c("black", "blue", "red"), bty = "n")

# power spectral density before and after filtering (dB)
par(mfrow=c(1,3), mar = c(5, 4.5, 4, 2) + 0.1)
eegpsd(y, Fs = n, upper = 50, t = "b",
       main = "Before Filtering", lwd = 2)
rect(2, -63, 15, 1, col = rgb(0.5,0.5,0.5,1/4))
legend("topright", legend = "2-15 Hz Filter", 
       fill = rgb(0.5,0.5,0.5,1/4), bty = "n")
eegpsd(yf.but, Fs = n, upper = 50, t = "b",
       main = "After Butterworth Filter", lwd = 2)
eegpsd(yf.fir, Fs = n, upper = 50, t = "b",
       main = "After FIR Filter", lwd = 2)
       
# power spectral density before and after filtering (mv^2)
par(mfrow=c(1,3), mar = c(5, 4.5, 4, 2) + 0.1)
eegpsd(y, Fs = n, upper = 50, unit = "mV^2", t = "b",
       main = "Before Filtering", lwd = 2)
rect(2, 0, 15, 1.05, col = rgb(0.5,0.5,0.5,1/4))
legend("topright", legend = "2-15 Hz Filter", 
       fill = rgb(0.5,0.5,0.5,1/4), bty = "n")
eegpsd(yf.but, Fs = n, upper = 50, unit = "mV^2", t = "b",
       main = "After Butterworth Filter", lwd = 2)
eegpsd(yf.fir, Fs = n, upper = 50, unit = "mV^2", t = "b",
       main = "After FIR Filter", lwd = 2)