tests/testthat/test_Filter_Analysis_Functions.R

In gjmvanboxtel/gsignal: Signal Processing

# gsignal filter analysis functions
library(gsignal)
library(testthat)

tol <- 1e-6

# -----------------------------------------------------------------------
# freqs()

test_that("parameters to freqs() are correct", {
  expect_error(freqs())
  expect_error(freqs(1))
  expect_error(freqs(1, 2))
  expect_error(freqs(1, 2, 3, plot = 'invalid'))
})

test_that("freqs() tests are correct", {
  h <- freqs(1, 1, 1, plot = FALSE)
  expect_equal(h, 1 + 0i)

  h <- freqs(c(1, 2), c(1, 1), 1:4, plot = FALSE)
  expect_equal(h, c(1.5-0.5i, 1.2-0.4i, 1.1-0.3i, 1.058824-0.235294i), tolerance = tol)
  
})

# -----------------------------------------------------------------------
# fwhm()

test_that("parameters to fwhm() are correct", {
  expect_error(fwhm())
  expect_error(fwhm(1))
  expect_error(fwhm(1, 2, 3, 4, 5))
  expect_error(fwhm(array(1)))
  expect_error(fwhm(1, c(1,2)))
  expect_error(fwhm(1, array(1:12, dim = c(2, 3, 2))))
  expect_error(fwhm(1, 2, ref = 'invalid'))
  expect_error(fwhm(1, 2, level = 'invalid'))
})

test_that("fwhm() tests are correct", {

  x <- seq(-pi, pi, 0.001)
  y <- cos(x)
  expect_equal(fwhm(x, y), 2 * pi / 3, tolerance = tol)
  
  expect_equal(fwhm(y = -10:10), 0L)
  expect_equal(fwhm(y = rep(1L, 50)), 0L)

  x <- seq(-20, 20, 1)
  y1 <- -4 + rep(0L, length(x)); y1[4:10] <- 8
  y2 <- -2 + rep(0L, length(x)); y2[4:11] <- 2
  y3 <-  2 + rep(0L, length(x)); y3[5:13] <- 10
  expect_equal(fwhm(x, cbind(y1, y2, y3)), c(20 / 3, 7.5, 9.25),
               tolerance = tol)

  x <- 1:3
  y <- c(-1, 3, -1)
  expect_equal(fwhm(x, y), 0.75)
  expect_equal(fwhm(x, y, 'max'), 0.75)
  expect_equal(fwhm(x, y, 'zero'), 0.75)
  expect_equal(fwhm(x, y, 'middle'), 1L)
  expect_equal(fwhm(x, y, 'min'), 1L)

  x <- 1:3
  y <- c(-1, 3, -1)
  expect_equal(fwhm(x, y, level = 0.1), 1.35)
  expect_equal(fwhm(x, y, ref = 'max', level = 0.1), 1.35)
  expect_equal(fwhm(x, y, ref = 'min', level = 0.1), 1.40)
  expect_equal(fwhm(x, y, ref = 'abs', level = 2.5), 0.25)
  expect_equal(fwhm(x, y, ref = 'abs', level = -0.5), 1.75)
  
  x <- -5:5
  y <- 18 - x * x
  expect_equal(fwhm(y = y), 6)
  expect_equal(fwhm(x, y), 6)
  expect_equal(fwhm(x, y, 'min'), 7)
  
})

# -----------------------------------------------------------------------
# freqz()

test_that("parameters to freqz() are correct", {
  expect_error(freqz())
  expect_error(freqz('invalid'))
  expect_error(freqz(NA, 1))
  expect_error(freqz(1, NA))
  expect_error(freqz(1, 1, -1))
  expect_error(freqz(1, 1, 1, FALSE, 0))
})

test_that("freqz() tests are correct", {
  
  # test correct values and fft-polyval consistency
  # butterworth filter, order 2, cutoff pi/2 radians
  b <- c(0.292893218813452, 0.585786437626905, 0.292893218813452)
  a <- c(1, 0, 0.171572875253810)
  hw <- freqz(b, a, 32)
  expect_equal(Re(hw$h[1]), 1)
  expect_equal(abs(hw$h[17])^2, 0.5)
  expect_equal(hw$h, freqz(b, a, hw$w)$h,
               tolerance = tol)  # fft should be consistent with polyval
  
  # test whole-half consistency
  b <- c(1, 1, 1)/3  # 3-sample average
  hw <- freqz(b, 1, 32, whole = TRUE)
  expect_equal(hw$h[2:16], Conj(hw$h[32:18]), tolerance = tol)
  hw2 <- freqz(b, 1, 16, whole = FALSE)
  expect_equal(hw$h[1:16], hw2$h, tolerance = tol)
  expect_equal(hw$w[1:16], hw2$w, tolerance = tol)
  
  # test sampling frequency properly interpreted
  b <- c(1, 1, 1) / 3; a <- c(1, 0.2)
  hw <- freqz(b, a, 16, fs = 320)
  expect_equal(hw$w, (0:15) * 10)
  hw2 <- freqz(b, a, (0:15) * 10, fs = 320)
  expect_equal(hw2$w, (0:15) * 10)
  expect_equal(hw$h, hw2$h, tolerance = tol)
  hw3 <- freqz(b, a, 32, whole = TRUE, fs = 320)
  expect_equal(hw3$w, (0:31) * 10)
  
  # Github Issue #9
  sys <- Arma(b = c(0.000731569700125585, 0, -0.00292627880050234, 0, 
                    0.00438941820075351, 0, -0.00292627880050234,
                    0, 0.000731569700125585),
              a = c(1, -7.04203950456817, 21.7283807572297,
                    -38.3801847495478, 42.4586614063362, -30.1289923970756,
                    13.3939745471932, -3.41070179210893, 0.380901732541468)
              )
  hw <- freqz(sys, 5L)
  expect_equal(hw$h, c(3.255208e-06+0.000000e+00i, -9.759891e-04+1.062529e-01i,
                       3.328587e-03+2.658062e-03i,  3.101129e-04+1.143374e-04i,
                       1.305432e-05+2.075128e-06i), tol = 1e-3)
  
  hw <- freqz(sys, c(0, 1))
  expect_equal(hw$w, c(0, 1))
  expect_equal(hw$h, c(3.25520833e-06+0i, 8.25367425e-03+0.01047173i), tol = 1e-5)
  
})


# -----------------------------------------------------------------------
# grpdelay()

test_that("parameters to grpdelay() are correct", {
  expect_error(grpdelay())
  expect_error(grpdelay('invalid'))
})

test_that("grpdelay() tests are correct", {
  
  gd1 <- grpdelay(c(0, 1))
  gd2 <- grpdelay(c(0, 1), 1)
  expect_equal(gd1$gd, gd2$gd)

  gd <- grpdelay(c(0, 1), 1, 4)
  expect_equal(gd$gd, rep(1L, 4))
  expect_equal(gd$w, pi/4 * 0:3, tolerance = tol)

  gd <- grpdelay(c(0, 1), 1, 4, whole = TRUE)
  expect_equal(gd$gd, rep(1L, 4))
  expect_equal(gd$w, pi/2 * 0:3, tolerance = tol)

  gd <- grpdelay(c(0, 1), 1, 4, fs = 0.5)
  expect_equal(gd$gd, rep(1L, 4))
  expect_equal(gd$w, 1/16 * 0:3, tolerance = tol)

  gd <- grpdelay(c(0, 1), 1, 4, TRUE, 1)
  expect_equal(gd$gd, rep(1L, 4))
  expect_equal(gd$w, 1/4 * 0:3)

  gd <- grpdelay(c(1, -0.9i), 1, 4, TRUE, 1)
  gd0 <- 0.447513812154696; gdm1 <- 0.473684210526316
  expect_equal(gd$gd, c(gd0, -9, gd0, gdm1), tolerance = tol)
  expect_equal(gd$w, 1/4 * 0:3)
  
  gd <- grpdelay(1, c(1, .9), n = 2 * pi * c(0, 0.125, 0.25, 0.375))
  expect_equal(gd$gd, c(-0.47368, -0.46918, -0.44751, -0.32316), tolerance = 1e-5)
  
  gd <- grpdelay(1, c(1, .9), c(0, 0.125, 0.25, 0.375), fs = 1)
  expect_equal(gd$gd, c(-0.47368, -0.46918, -0.44751, -0.32316), tolerance = 1e-5)
  
  gd <- grpdelay(c(1, 2), c(1, 0.5, .9), 4)
  expect_equal(gd$gd, c(-0.29167, -0.24218, 0.53077, 0.40658), tolerance = 1e-5)

  b1 <- c(1, 2); a1f <- c(0.25, 0.5, 1); a1 <- rev(a1f)
  gd1 <- grpdelay(b1, a1, 4)$gd
  gd <- grpdelay(conv(b1, a1f), 1, 4)$gd - 2
  expect_equal(gd, gd1, tolerance = 1e-5)
  expect_equal(gd, c(0.095238, 0.239175, 0.953846, 1.759360), tolerance = 1e-5)
  
  a <- c(1, 0, 0.9)
  b <- c(0.9, 0, 1)
  dh <- grpdelay(b, a, 512, 'whole')$gd
  da <- grpdelay(1, a, 512, 'whole')$gd
  db <- grpdelay(b, 1, 512, 'whole')$gd
  expect_equal(dh, db + da, ttolerance = 1e-5)

})

# -----------------------------------------------------------------------
# impz()

test_that("parameters to impz() are correct", {
  expect_error(impz())
  expect_error(impz('invalid'))
})

test_that("impz() tests are correct", {
  
  xt <- impz(1, c(1, -1, 0.9), 100)
  expect_equal(length(xt$t), 100L)
  expect_equal(xt$t, 0:99)
  
  xt <- impz(1, c(1, -1, 0.9), 0:101)
  expect_equal(length(xt$t), 102L)
  expect_equal(xt$t, 0:101)
  
})