add peak frequency function

Author: Bart

NAMESPACE

@@ -14,6 +14,7 @@ export(is_acc)
 export(is_uniform)
 export(n_axis)
 export(n_samples)
+export(peak_frequency)
 export(plot_time)
 import(vctrs)
 importFrom(stats,na.omit)

R/peak_frequency.R

@@ -0,0 +1,48 @@
+
+#' Calculate the peak frequency per axis for accelration bursts
+#'
+#' @param x An `acc` vector
+#' @param resolution A scalar with the units Hertz
+#'
+#' @returns returns a list with the same length as `x` with the peak frequency per axis
+#'
+#' @details
+#'
+#' To increase the resolution of the result zero padding can be used. This can be controlled uing the resolution argument.
+#' Note that increasing resolution without increasing the number of samples in a acceleration burst has a limited ability to get closer to the true frequency.
+#'
+#' @export
+#'
+#' @examples
+#'   a<-acc(list(cbind(z=cos(1:200/(80/(pi*2))),
+#'   x=sin(1:200/(5/(pi*2))))), units::set_units(400,'Hz'))
+#'   peak_frequency(a)
+#'   peak_frequency(a, units::set_units(.25, "Hz"))
+#'   # Increasing resolution more
+#'   peak_frequency(a, units::set_units(.005, "Hz"))
+#'   a<-acc(list(cbind(z=cos(80+1:200/(80/(pi*2))),
+#'   x=sin((1:200)/(5/(pi*2))))), units::set_units(400,'Hz'))
+#'   peak_frequency(a, units::set_units(.005, "Hz"))
+
+
+
+peak_frequency<-function(x, resolution=NA){
+  assertthat::assert_that(inherits(x,'acc'))
+  s<-!is.na(x)
+  res<-rep(list(NULL),length(x))
+  if(all(!s)){
+    return(res)
+  }
+  b<-field(x[s],"bursts")
+  b<-    purrr::map2(b,purrr::map(b, inherits,'units'), ~if(.y){units::drop_units(.x)}else{.x})
+  b_centered<-purrr::map2(b,purrr::map(b, colMeans), ~ t((.x))-.y)
+  if(!is.na(resolution)){
+    to_pad<-units::drop_units(field(x[s],"frequency")/resolution)-n_samples(x[s])
+    b_centered<-purrr::map2(b_centered, to_pad, ~cbind(.x,matrix(0, ncol=.y, nrow=nrow(.x))))
+  }
+  b_mod<-purrr::map(b_centered, ~ do.call(rbind,lapply(apply(.x,1, fft, simplify = F),  Mod))[,1:ceiling(ncol(.x)/2), drop=F])
+  peak<- purrr::map(b_mod, ~ apply(.x,1, which.max ))
+  peak_freq<-purrr::pmap(list(peak,field(x[s],"frequency"),purrr::map(b_mod, ncol)), ~ (..1-1)*(..2/..3/2))
+  res[s]<-peak_freq
+  res
+}

man/peak_frequency.Rd

@@ -0,0 +1,72 @@
+test_that("Correct frequency found", {
+  a<-new_acc(list(matrix(sin(1:200/(50/(pi*2))))), units::set_units(200,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(4,'Hz')))
+})
+
+test_that("Multiple axis peak freq and changing freq", {
+  a<-new_acc(list(cbind(z=cos(1:200/(100/(pi*2))),
+                        x=sin(1:200/(5/(pi*2))))), units::set_units(200,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=2,x=40),'Hz')))
+
+  a<-new_acc(list(cbind(z=cos(1:200/(100/(pi*2))),
+                        x=sin(1:200/(5/(pi*2))))), units::set_units(100,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=1,x=20),'Hz')))
+  a<-new_acc(list(cbind(z=cos(1:200/(100/(pi*2))),
+                        x=sin(1:200/(5/(pi*2))))), units::set_units(400,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=4,x=80),'Hz')))
+})
+
+test_that("length does not influnce result", {
+  a<-new_acc(list(cbind(z=cos(1:199/(100/(pi*2))),
+                        x=sin(1:199/(5/(pi*2))))), units::set_units(200,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=2,x=40),'Hz')))
+  a<-new_acc(list(cbind(z=cos(1:199/(100/(pi*2))),
+                        x=sin(1:199/(5/(pi*2))))), units::set_units(100,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=1,x=20),'Hz')))
+  a<-new_acc(list(cbind(z=cos(1:199/(100/(pi*2))),
+                        x=sin(1:199/(5/(pi*2))))), units::set_units(400,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=4,x=80),'Hz')))
+})
+
+test_that("Multiple axis peak freq intercept does not matter", {
+  a<-new_acc(list(cbind(z=-3+.1*cos(1:200/(100/(pi*2))),x=3*(2+sin(1:200/(50/(pi*2)))))), units::set_units(200,'Hz'))
+  expect_equal(peak_frequency(a), list(units::set_units(c(z=2,x=4),'Hz')))
+})
+
+test_that("NA returns empty",{
+  expect_equal(peak_frequency(new_acc(list(NULL), frequency = NA)), list(NULL))
+  expect_equal(peak_frequency(new_acc(list(NULL, NULL), frequency = c(NA, NA))), list(NULL, NULL))
+  expect_equal(peak_frequency(new_acc(list(NULL,matrix(sin(1:200/(50/(pi*2)))),  NULL),
+                                      frequency = units::set_units(c(NA,200, NA),'Hz'))),
+               list(NULL,units::set_units(4,"Hz"), NULL))
+
+})
+
+test_that("Resolution alows to identify partial frequencies", {
+  a<-new_acc(list(cbind(z=cos(1:200/(80/(pi*2))),
+                        x=sin(1:200/(5/(pi*2))))), units::set_units(200,'Hz'))
+  expect_equal(peak_frequency(a),
+               list(units::set_units(c(z=3,x=40),'Hz')))
+  expect_equal(peak_frequency(a, resolution = units::set_units(.5,'Hz')),
+               list(units::set_units(c(z=2.5,x=40),'Hz')))
+  expect_equal(peak_frequency(a, resolution = units::set_units(.25,'Hz')),
+               list(units::set_units(c(z=2.5,x=40),'Hz')))
+})
+
+
+test_that("Resolution alows to identify partial frequencies", {
+  a<-new_acc(list(matrix(runif(100), ncol=10)), units::set_units(23,'Hz'))
+  expect_equal((((unlist(peak_frequency(a, resolution = units::set_units(.005,'Hz')))/.005)+.5)%%1)-.5,
+               rep(0,10) )
+  expect_equal((((unlist(peak_frequency(a, resolution = units::set_units(.025,'Hz')))/.025)+.5)%%1)-.5,
+               rep(0,10) )
+#  expect_equal((((unlist(peak_frequency(a, resolution = units::set_units(.035,'Hz')))/.035)+.5)%%1)-.5,
+#               rep(0,10) )
+  a<-new_acc(list(matrix(runif(1000), ncol=10)), units::set_units(23,'Hz'))
+  expect_equal((((unlist(peak_frequency(a, resolution = units::set_units(.005,'Hz')))/.005)+.5)%%1)-.5,
+               rep(0,10) )
+  expect_equal((((unlist(peak_frequency(a, resolution = units::set_units(.025,'Hz')))/.025)+.5)%%1)-.5,
+               rep(0,10) )
+
+  })
+

tests/testthat/test-peak_frequency.R

@@ -83,14 +83,14 @@ The same can be done with ornitella data.
 ```{r ornitella}
 d <- as.Date("2021-3-3")
 
-lbbg_data <- move2::movebank_download_study("LBBG_ZEEBRUGGE - Lesser black-backed gulls ",
+lbbg_data <- movebank_download_study("LBBG_ZEEBRUGGE - Lesser black-backed gulls ",
   sensor_type_id = c("gps", "acceleration"),
   timestamp_start = as.POSIXct(d),
   timestamp_end = as.POSIXct(d + 1)
 )
 lbbg_data <- lbbg_data %>%
-  mutate(acceleration = as_acc(.))  %>%
-  filter(!is.na(acceleration) | sensor_type_id==653)%>%
+  mutate(acceleration = as_acc(.)) %>%
+  filter(!is.na(acceleration) | sensor_type_id==653) %>%
   select(acceleration)
 lbbg_data
 ```
@@ -113,6 +113,19 @@ Alternatively you can plot one burst directly, in this case showing the clear wi
 ```{r plot_one, fig.width=6}
 plot_time(lbbg_data$acceleration[340], mt_time(lbbg_data)[340])
 ```
+
+Using the `peak_frequency` function we can identify the wing beat frequency in this acceleration burst.
+```{r pf}
+peak_frequency(lbbg_data$acceleration[340])
+```
+
+Due to the nature of the fast fourier transform there is a limit to the frequencies that can be identified. To resolve this zero padding can be used, this is controlled by the resolution argument. By setting this a more detailed estimate of the wing beat frequency can be obtained.
+
+```{r pf_pad}
+peak_frequency(lbbg_data$acceleration[340], 
+               resolution = units::set_units(.01,"Hz"))
+```
+
 We can also visualize one burst in three dimensions showing the cyclic patterns in the acceleration.
 
 ```{r plot3d, fig.width=6, fig.height=6}