include physical units treatment

DESCRIPTION

@@ -25,7 +25,6 @@ Depends:
     R (>= 4.1)
 Suggests:
     ggplot2 (>= 3.3.0),
-    hms (>= 1.1.0),
     viridis (>= 0.6.0),
     rlang (>= 1.0.0),
     testthat (>= 3.0.0)
@@ -38,3 +37,5 @@ Config/testthat/edition: 3
 URL: https://lschneiderbauer.github.io/fCWTr/, https://github.com/lschneiderbauer/fCWTr
 Roxygen: list(markdown = TRUE)
 BugReports: https://github.com/lschneiderbauer/fCWTr/issues
+Imports: 
+    units (>= 0.8)

NAMESPACE

@@ -4,6 +4,7 @@ S3method(as.data.frame,fcwtr_scalogram)
 S3method(plot,fcwtr_scalogram)
 export(fcwt)
 export(fcwt_batch)
+export(u)
 importFrom(graphics,plot)
 importFrom(utils,head)
 importFrom(utils,setTxtProgressBar)

NEWS.md

@@ -1,5 +1,7 @@
 # fcwtr 0.2.9999
 
+-   Include improved physical unit treatment with the 'units' package. Frequency and time parameters can now be passed a bestowed unit.
+
 -   Overhaul `fcwt_batch()`: in some cases, time slices were accidentally discarded, this should now be fixed.
 
 -   Fix a bug where `fcwt(..., n_freqs = 1, ...)` leads to a faulty data representation.

R/fcwt.R

@@ -22,7 +22,9 @@
 #' @param sample_freq
 #'  Sampling rate of input time series. This number primarily establishes
 #'  a connection to physical units which is used in other frequency definitions
-#'  as well as the units of the output data.
+#'  as well as the units of the output data. Expects either a value with frequency
+#'  units, generated with `u()`, or a pure number, in which case it is
+#'  interpreted in units of 'Hertz'.
 #'
 #' @param n_freqs
 #'  Number of frequency bins generated by the CWT. The frequencies
@@ -33,7 +35,8 @@
 #'  Optionally specifies the frequency range `[freq_end, freq_begin]`. If not
 #'  specified the maximal meaningful frequency range, depending on the input signal,
 #'  is taken.
-#'  The range and `sample_freq` need to be specified in the same units.
+#'  A frequency-valued number, generated with `u()`, or a pure number, that is
+#'  interpreted in units of 'Hertz'.
 #'
 #' @param freq_scale
 #' Should the frequency scale be linear or logarithmic?
@@ -73,9 +76,9 @@
 #' res <-
 #'   fcwt(
 #'     ts_sin_440,
-#'     sample_freq = 44100,
-#'     freq_begin = 50,
-#'     freq_end = 1000,
+#'     sample_freq = u(44.1, "kHz"),
+#'     freq_begin = u(50, "Hz"),
+#'     freq_end = u(1000, "Hz"),
 #'     n_freqs = 10,
 #'     sigma = 5
 #'   )
@@ -91,22 +94,27 @@ fcwt <- function(signal,
                  remove_coi = TRUE,
                  n_threads = 2L) {
   stopifnot(is.numeric(signal))
-  stopifnot(is.numeric(sample_freq), sample_freq > 0)
-  stopifnot(is.numeric(freq_begin), freq_begin > 0)
   stopifnot(is.numeric(freq_end), freq_end > freq_begin)
   stopifnot(is.numeric(n_freqs), n_freqs > 0)
   stopifnot(is.numeric(sigma), sigma > 0)
   stopifnot(is.numeric(n_threads))
-  # stopifnot(is.logical(abs))
+
+  sample_freq <- as_freq(sample_freq)
+  freq_begin <- as_freq(freq_begin)
+  freq_end <- as_freq(freq_end)
+  stopifnot(sample_freq > u(0, "Hz"))
+  stopifnot(freq_begin > u(0, "Hz"))
 
   freq_scale <- match.arg(freq_scale)
 
   freq_scale_lgl <- (freq_scale == "linear")
   stopifnot(isTRUE(freq_scale_lgl) || isFALSE(freq_scale_lgl))
 
+  # we expect the numbers in Hz here
   output <-
     fcwt_raw(
-      as.numeric(signal), as.integer(sample_freq), freq_begin, freq_end,
+      as.numeric(signal), as.integer(hz(sample_freq)),
+      hz(freq_begin), hz(freq_end),
       as.integer(n_freqs), sigma, as.integer(n_threads),
       freq_scale_lgl,
       FALSE,

R/fcwt_batch.R

@@ -4,9 +4,9 @@
 #' processing junks of the input signal and keeping only low-resolution output
 #' data to preserve memory.
 #' This is only useful for very long signals whose output does not fit into
-#' memory as a whole.
-#' It should not be used on short signals since boundary artefacts are discarded
-#' (and those potentially dominate for short sequences).
+#' the available memory as a whole.
+#' It should not be used on short signals since boundary artefacts are
+#' automatically discarded (and those potentially dominate for short signals).
 #'
 #' @details
 #' In case of input sequences that exceed the a certain size, the output
@@ -27,9 +27,9 @@
 #' This function splits up the input sequence into batches, processes each
 #' batch separately, reduces the time resolution, and adds the outputs together.
 #'
-#' Attention: Boundary artefacts are removed, so some high frequency information
-#' at the beginning and the end of the sequence is lost. (The amount depends on
-#' the minimal frequency captured `min_freq`.)
+#' Attention: In contrast to `fcwt()` boundary artefacts are automatically removed,
+#' so some high frequency information at the beginning and the end of the
+#' sequence is lost. (The amount depends on the minimal frequency captured `min_freq`.)
 #'
 #' @param max_batch_size
 #'  The maximal batch size that is used for splitting up the input sequence.
@@ -38,9 +38,10 @@
 #'  processing it.
 #'  The actual batch size depends on the requested `time_resolution`.
 #' @param time_resolution
-#'  The time resolution in inverse units of `sample_freq` of the result.
+#'  The time resolution in physical units, generated by `u()`, or a pure number,
+#'  which is interpreted in unit of seconds.
 #'  Memory consumption is directly related to that.
-#'  Can not be higher than the time resolution of the input signal.
+#'  Can not be smaller than the time resolution of the input signal.
 #' @param progress_bar
 #'  Monitoring progress can sometimes be useful when performing time consuming
 #'  operations. Setting `progress_bar = TRUE` enables printing a progress
@@ -83,6 +84,16 @@ fcwt_batch <- function(signal,
                        n_threads = 2L,
                        progress_bar = FALSE) {
 
+  time_resolution <- as_sec(time_resolution)
+  sample_freq <- as_freq(sample_freq)
+  freq_begin <- as_freq(freq_begin)
+  freq_end <- as_freq(freq_end)
+
+  stopifnot(
+    time_resolution > u(0, "s"),
+    time_resolution >= 1 / sample_freq
+  )
+
   # aggregation window
   w <- wnd_from_resolution(time_resolution, sample_freq)
 

R/fcwtr_scalogram.R

@@ -1,6 +1,5 @@
 new_fcwtr_scalogram <- function(matrix, coi_mask, sample_freq, freq_begin, freq_end,
                                 freq_scale, sigma) {
-
   obj <-
     structure(
       matrix,
@@ -13,26 +12,16 @@ new_fcwtr_scalogram <- function(matrix, coi_mask, sample_freq, freq_begin, freq_
       sigma = sigma
     )
 
-  dimnames(obj) <-
-    list(
-      seq_along(matrix[, 1]) - 1,
-      seq2(
-        freq_end, freq_begin, length.out = dim(matrix)[[2]],
-        scale = freq_scale
-      )
-    )
-
   obj
 }
 
 fcwtr_scalogram <- function(matrix, sample_freq, freq_begin, freq_end,
                             freq_scale, sigma) {
-
   stopifnot(is.matrix(matrix))
   stopifnot(freq_scale %in% c("linear", "log"))
-  stopifnot(is.numeric(sample_freq))
-  stopifnot(is.numeric(freq_begin))
-  stopifnot(is.numeric(freq_end))
+  stopifnot(inherits(sample_freq, "units"))
+  stopifnot(inherits(freq_begin, "units"))
+  stopifnot(inherits(freq_end, "units"))
   stopifnot(is.numeric(sigma))
 
   coi_mask <-
@@ -42,11 +31,14 @@ fcwtr_scalogram <- function(matrix, sample_freq, freq_begin, freq_end,
       sample_freq = sample_freq,
       freq_begin = freq_begin,
       freq_end = freq_end,
+      freq_scale = freq_scale,
       sigma = sigma
     )
 
-  new_fcwtr_scalogram(matrix, coi_mask, sample_freq, freq_begin, freq_end,
-                      freq_scale, sigma)
+  new_fcwtr_scalogram(
+    matrix, coi_mask, sample_freq, freq_begin, freq_end,
+    freq_scale, sigma
+  )
 }
 
 sc_set_coi_na <- function(x) {
@@ -66,19 +58,20 @@ sc_coi_mask <- function(x) {
 #' @return A boolean matrix of the same dimensions as `x`. `TRUE` values
 #'         indicate values inside the boundary "cone of influence".
 #' @noRd
-coi_mask <- function(dim_t, dim_f, sample_freq, freq_begin, freq_end, sigma) {
+coi_mask <- function(dim_t, dim_f, sample_freq, freq_begin, freq_end,
+                     freq_scale, sigma) {
   # The standard deviation Σ of a the Gauß like wave packet at frequency f
   # and sampling frequency f_s with given σ is given by
   # Σ = σ / sqrt(2) f_s / f
   # we choose 4Σ to define the support of a wave packet
   # (and so boundary effects are expected to occur until 2Σ)
-  coi_pred <- \(f, t) t * f < sqrt(2) * sigma
+  coi_pred <- \(f, t) ddu(t * f) < sqrt(2) * sigma
 
   # express in dimensionless quantities
   t <- rep(1:dim_t, times = dim_f)
   f <-
     rep(
-      seq(freq_end, freq_begin, length.out = dim_f) / sample_freq,
+      seq2(freq_end, freq_begin, length.out = dim_f, scale = freq_scale) / sample_freq,
       each = dim_t
     )
 
@@ -108,7 +101,7 @@ sc_dim_time <- function(x) {
 sc_coi_time_interval <- function(x) {
   stopifnot(inherits(x, "fcwtr_scalogram"))
 
-  #unique(which(is.na(x), arr.ind = TRUE)[, 1])
+  # unique(which(is.na(x), arr.ind = TRUE)[, 1])
 
   full_info_rows <- which(rowSums(sc_coi_mask(x)) == 0)
 
@@ -227,7 +220,7 @@ tbind <- function(..., deparse.level = 1) {
 #' @return
 #'  A [data.frame()] object representing the scalogram data with four columns:
 #' \describe{
-#'   \item{time_ind}{An integer index uniquely identifying time slices.}
+#'   \item{time_index}{An integer index uniquely identifying time slices.}
 #'   \item{time}{The time difference to the first time slice in physical units.
 #'               The time unit is the inverse of the frequency unit chosen by the user
 #'               for the `sample_freq` argument of [fcwt()].}
@@ -243,19 +236,25 @@ tbind <- function(..., deparse.level = 1) {
 #'   sample_freq = 44100,
 #'   n_freqs = 10
 #' ) |>
-#' as.data.frame() |>
-#' head()
+#'   as.data.frame() |>
+#'   head()
 #'
 #' @export
 as.data.frame.fcwtr_scalogram <- function(x, ...) {
-  df <- as.data.frame(as.table(x), stringsAsFactors = FALSE)
-  names(df) <- c("time_ind", "freq", "value")
-  df[["time_ind"]] <- as.integer(df[["time_ind"]])
-  df[["freq"]] <- as.numeric(df[["freq"]])
-
-  df[["time"]] <- df[["time_ind"]] / attr(x, "sample_freq")
+  dim_t <- seq_along(x[, 1])
+  dim_f <-
+    seq2(
+      attr(x, "freq_end"), attr(x, "freq_begin"),
+      length.out = dim(x)[[2]],
+      scale = attr(x, "freq_scale")
+    )
 
-  df[, c("time_ind", "time", "freq", "value")]
+  data.frame(
+    time_index = dim_t[row(x)],
+    time = dim_t[row(x)] / attr(x, "sample_freq"),
+    freq = dim_f[col(x)],
+    value = c(x)
+  )
 }
 
 #' Scalogram plotting
@@ -313,20 +312,22 @@ autoplot.fcwtr_scalogram <- function(object, n = 1000, ...) {
 
   freq_scale <- attr(object, "freq_scale")
 
-  scale_y <-
+  transform <-
     if (freq_scale == "log") {
-      ggplot2::scale_y_log10
+      "log10"
     } else {
-      ggplot2::scale_y_continuous
+      "identity"
     }
 
+  df <- as.data.frame(sc_agg(object, wnd_from_target_size(n, object)))
+
   # first aggregate the time series,
   # since we cannot really see too much resolution anyways
-  as.data.frame(sc_agg(object, wnd_from_target_size(n, object))) |>
+  df |>
     ggplot(aes(x = .data$time, y = .data$freq, fill = .data$value)) +
     geom_raster() +
     viridis::scale_fill_viridis(discrete = FALSE) +
-    scale_y(name = "Frequency") +
-    ggplot2::scale_x_time(name = "Time") +
+    units::scale_x_units(name = "Time", unit = "s") +
+    units::scale_y_units(name = "Frequency", transform = transform) +
     ggplot2::theme_minimal()
 }