1. Preparation
Let’s load the package and generate a sample signal. The signal is a mixture of 2Hz, 10Hz, 60Hz, and noise
sample_rate <- 200
t <- seq(0, 10, by = 1 / sample_rate)
noise <- cumsum(rnorm(length(t), sd = 0.4)) / sqrt(seq_along(t))
x <- 0.2 * sin(t * 4 * pi) + 0.02 * sin(t * 20 * pi) +
0.005 * sin(t * 120 * pi) + noise
# Plot the signal
ravetools::diagnose_channel(
x,
srate = sample_rate,
name = "Sample signal",
window = 2 * sample_rate,
noverlap = sample_rate
)
2. FIR Filter
Construct a FIR filter with low pass at 3Hz and (around) 0.5Hz transition bandwidth:
# Low-pass filter
f1 <- ravetools::design_filter(
data_size = length(x),
sample_rate = sample_rate,
low_pass_freq = 3, low_pass_trans_freq = 0.5
)
ravetools::diagnose_filter(b = f1$b, a = f1$a, fs = sample_rate, n = 3000)
Apply the filter to the signal x using filtfilt method:
y1 <- ravetools::filtfilt(b = f1$b, a = f1$a, x = x)
par(mfrow = c(1, 2))
# compare the results
plot(t, x, type = "l", xlab = "Time", ylab = "",
main = "Mixture of 2, 10, and 60Hz", xlim = c(0, 1))
lines(t, y1, col = "red")
legend(
"topleft", c("Input", "Low: 3Hz"),
col = c(par("fg"), "red"), lty = 1
)
# plot welch-periodogram
ravetools::pwelch(x, fs = sample_rate, window = sample_rate * 2,
noverlap = sample_rate, plot = 1, ylim = c(-100, 10))
ravetools::pwelch(y1, fs = sample_rate, window = sample_rate * 2,
noverlap = sample_rate, plot = 2, col = "red")