# tests/testthat/test-filter.R In ravetools: Signal and Image Processing Toolbox for Analyzing Intracranial Electroencephalography Data

```test_that("C++ signal filter", {

initialize_filter <- function(b, a, x) {
# make sure a, b share the same order
na <- length(a)
nb <- length(b)

if( na > nb ) {
b <- c(b, rep(0, na - nb))
n <- na
} else {
a <- c(a, rep(0, nb - na))
n <- nb
}

# length of edge transients
nf <- max(1, 3 * (n - 1))

# compute the initial condition if n > 1
if( n > 1 ) {
z1 <- diag(1, n - 1) - cbind( -a[-1], rbind(diag(1, n - 2), 0) )
z2 <- b[-1] - b[1] * a[-1]
z <- solve(z1, z2)
} else {
z <- numeric(0)
}
list(
a = a,
b = b,
z = z,
nfilt = n,
nfact = nf
)
}
myFilter <- function(b, a, x, z) {
# make sure a, b share the same order
na <- length(a)
nb <- length(b)
if( na > nb ) {
b <- c(b, rep(0, na - nb))
n <- na
} else {
a <- c(a, rep(0, nb - na))
n <- nb
}
b <- b / a[1]
a <- a / a[1]
y <- rep(0, length(x))
if(missing(z)) {
z <- rep(0, n - 1)
}

for(m in seq_along(y)) {
xm <- x[m]
y[m] <- b[1] * xm + z[1]
ym <- y[m]
for( i in 2: (n-1)) {
z[ i-1 ] <- b[i] * xm + z[i] - a[i] * ym
}
z[n-1] <- b[n] * xm - a[n] * ym
}
list(y, z)
}

bf <- signal::butter(10, c(0.15, 0.3))
t <- seq(0, 1, by = 0.005)
x <- as.double(sin(2*pi*t*2.3)) + rnorm(length(t), mean = t)
b <- as.double(bf\$b)
a <- as.double(bf\$a)
nx <- length(x)
init <- initialize_filter(b, a, x)
nfact <- init\$nfact
z <- as.double(init\$z)

expected <- myFilter(b,a,x,z)
my_result <- ravetools:::cpp_filter(b,a,x,z)

testthat::expect_lt(
max(abs((expected[[1]] - my_result[[1]]) / (expected[[1]] + my_result[[1]] + 1))),
1e-3
)
testthat::expect_lt(
max(abs((expected[[2]] - my_result[[2]]) / (expected[[2]] + my_result[[2]] + 1))),
1e-3
)

})
```

## Try the ravetools package in your browser

Any scripts or data that you put into this service are public.

ravetools documentation built on June 22, 2024, 9:41 a.m.