Adds FilterSpec and FilterLogSpec

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: PaleoSpec
 Title: Spectral tools for the ECUS group
-Version: 0.2.91
+Version: 0.3
 Authors@R: c(
   person("Thomas", "Laepple", email = "[email protected]", role = c("aut", "cre")),
   person("Thomas", "Muench", email = "[email protected]", role = c("aut")),

NAMESPACE

@@ -11,6 +11,8 @@ export(BinTimeseries)
 export(ClosestElement)
 export(ColTransparent)
 export(DF2Spec)
+export(FilterSpec)
+export(FilterSpecLog)
 export(FirstElement)
 export(GetTransferFunction)
 export(GetVarFromSpectra)
@@ -40,6 +42,7 @@ export(SpecACF)
 export(SpecInterpolate)
 export(SpecMTM)
 export(SubsampleTimeseriesBlock)
+export(TrimNA)
 export(as.spec)
 export(as_spec_df)
 export(gg_spec)

NEWS.md

@@ -1,3 +1,12 @@
+# PaleoSpec 0.3
+
+* Add FilterSpec and FilterSpecLog as replacements / alternatives to LogSmooth
+
+ FilterSpec and FilterSpecLog uses ApplyFilter to (optionally) avoid loosing the
+ highest and lowest frequencies when smoothing a spectrum. Default behaviour is 
+ method 3: min roughness (reflect and invert ends). Avoids the artefact seen at 
+ high freqs with LogSmooth.
+
 # PaleoSpec 0.2.91
 
 * Return vector of dof from SpecMTM even when all == 2

R/FilterSpec.R

@@ -0,0 +1,250 @@
+#' Filter a Power Spectrum Object
+#'
+#' @param spec A spec object
+#' @param keep_low_f Keep filtered (smoothed) low frequencies or replace with unfiltered
+#' @inheritParams stats::spec.pgram
+#' @inheritParams ApplyFilter
+#' @return A spec object (list)
+#' @seealso [FilterSpecLog] for filtering with filter widths equal in log-space
+#' @export
+#'
+#' @examples
+#' ## Comparison of the four methods - for power spectra, methods 0, 2 or 3 make the most sense
+#' library(PaleoSpec)
+#'
+#' a <- 100
+#' b <- 1
+#' N <- 1e03
+#' set.seed(20230625)
+#' ts1 <- SimPLS(N, beta = b, alpha = a)
+#' sp1 <- SpecMTM(ts(ts1), bin.width = 1)
+#' LPlot(sp1)
+#' abline(log10(a), -b, col = "green")
+
+#' fl <- seq(3, 9, by = 2)
+#' sp1_f3_0 <- FilterSpec(sp1, spans = fl, method = 0)
+#' sp1_f3_1 <- FilterSpec(sp1, spans = fl, method = 1)
+#' sp1_f3_2 <- FilterSpec(sp1, spans = fl, method = 2)
+#' sp1_f3_3 <- FilterSpec(sp1, spans = fl, method = 3)
+#' sp1_f3_4 <- FilterSpec(sp1, spans = fl, method = 4)
+#'
+#' LPlot(sp1)
+#' LLines(sp1_f3_0, col = "blue")
+#' LLines(sp1_f3_1, col = "green", lty = 2)
+#' LLines(sp1_f3_2, col = "red", lty = 3)
+#' LLines(sp1_f3_3, col = "orange", lty = 4)
+#' LLines(sp1_f3_4, col = "gold", lty = 5)
+#'
+#' ## Comparison of keeping the filtered values in the reflected end portions or not
+#' sp1_f3_0T <- FilterSpec(sp1, spans = fl, method = 0, keep_low_f = TRUE)
+#' sp1_f3_0F <- FilterSpec(sp1, spans = fl, method = 0, keep_low_f = FALSE)
+
+#' LPlot(sp1_f3_0F)
+#' LLines(sp1_f3_0T, col = "red")
+
+#' sp1_f3_2T <- FilterSpec(sp1, spans = fl, method = 2, keep_low_f = TRUE)
+#' sp1_f3_2F <- FilterSpec(sp1, spans = fl, method = 2, keep_low_f = FALSE)
+
+#' LPlot(sp1_f3_2F)
+#' LLines(sp1_f3_2T, col = "red")
+FilterSpec <- function(spec, spans, method = 3, keep_low_f = TRUE) {
+  if (length(spec$dof) == 1) {
+    spec$dof <- rep(spec$dof, length(spec$freq))
+  }
+
+  dof0 <- spec$dof
+
+  kernel <- stats::kernel("modified.daniell", spans %/% 2)
+  filter <- kernel[-kernel$m:kernel$m]
+
+  spec_filt <- ApplyFilter(spec$spec, filter = filter, method = method)
+
+  if (keep_low_f == FALSE) {
+    # replace filtered spec with original in area where freqs have been reflected
+    i <- 1:ceiling(length(filter) / 2)
+    spec_filt[i] <- spec$spec[i]
+
+    iend <- length(spec$freq) - (i-1)
+
+    spec_filt[iend] <- spec$spec[iend]
+
+  }
+
+  spec$spec <- as.numeric(spec_filt)
+
+  # degrees of freedom of the kernel
+  df.kern <- stats::df.kernel(kernel)
+
+  spec$dof <- df.kern * spec$dof / 2
+
+  if (keep_low_f == FALSE) {
+
+    i <- 1:ceiling(length(filter) / 2)
+    spec$dof[i] <- dof0[i]
+
+    iend <- length(spec$freq) - (i-1)
+    spec$dof[iend] <- dof0[iend]
+
+  }
+
+  # Adjust DOF in reflected filter region
+  if (keep_low_f == TRUE){
+
+    fl <- length(filter)
+    i <- 1:ceiling(fl / 2)
+    iend <- length(spec$freq) - (i-1)
+
+
+    if (method %in% c(2,3)){
+      scl <- 2 * (fl - (i - 1)) / fl
+      spec$dof[i] <- spec$dof[i] / scl
+      spec$dof[iend] <- spec$dof[iend] / scl
+
+    }
+
+    if (method == 0){
+      # remove NA portion
+      spec$freq <- spec$freq[is.na(spec$spec) == FALSE]
+      spec$dof <- spec$dof[is.na(spec$spec) == FALSE]
+      spec$shape <- spec$shape[is.na(spec$spec) == FALSE]
+      spec$spec <- spec$spec[is.na(spec$spec) == FALSE]
+    }
+
+  }
+
+  spec$shape <- spec$dof / 2
+
+  spec <- AddConfInterval(spec)
+
+
+  return(spec)
+}
+
+
+
+#' Smooth a Spectrum with Evenly Spaced Bins in Logspace
+#'
+#' @param spec A spec object
+#' @inheritParams LogSmooth
+#' @inheritParams ApplyFilter
+#'
+#' @return A spec object (list)
+#' @seealso [LogSmooth()] for an alternative implementation of log spaced filtering
+#' @export
+#' @examples
+#' library(PaleoSpec)
+#'
+#' # simulate a timeseries with powerlaw power spectrum
+#' a <- 100
+#' b <- 1
+#' N <- 1e03
+#'
+#' set.seed(20230625)
+#' ts1 <- SimPLS(N, beta = b, alpha = a)
+#' sp1 <- SpecMTM(ts(ts1), bin.width = 1)
+#' LPlot(sp1)
+#' abline(log10(a), -b, col = "green")
+#' #
+#' sp1_f3_0 <- FilterSpecLog(sp1, method = 0)
+#' sp1_f3_2 <- FilterSpecLog(sp1, method = 2)
+#'
+#' LPlot(sp1)
+#' LLines(sp1_f3_0, col = "blue")
+#' LLines(sp1_f3_2, col = "green", lty = 3)
+#'
+#' sp1_df0.05 <- FilterSpecLog(sp1)
+#' sp1_df0.1 <- FilterSpecLog(sp1, df.log = 0.1)
+#'
+#' LPlot(sp1)
+#' LLines(sp1_df0.05, col = "blue")
+#' LLines(sp1_df0.1, col = "red")
+#'
+#' ## A combination of FilterSpec and FilterSpecLog
+#'
+#' sp1_FSL <- FilterSpecLog(sp1)
+#' sp1_FSL_FS <- FilterSpec(FilterSpecLog(sp1), spans = c(3, 5))
+#' sp1_FS_FSL <- FilterSpecLog(FilterSpec(sp1, spans = c(3, 5)))
+#' LPlot(sp1)
+#' LLines(sp1_FSL, col = "blue")
+#' LLines(sp1_FSL_FS, col = "red")
+#' LLines(sp1_FS_FSL, col = "green")
+FilterSpecLog <- function(spec,
+                          df.log = 0.05,
+                          spans = NULL,
+                          method = 3, f.res = 10){
+
+  GetFW <- function(spec, df.log) {
+    ((exp(df.log) - 1) * max(spec$freq)) / min(spec$freq)
+  }
+
+  if (length(spec$dof) == 1){
+    spec$dof <- rep(spec$dof, length(spec$freq))
+  }
+
+  if (is.null(spans)){
+    spans <- GetFW(spec, df.log = df.log)
+  }
+
+  # interpolate spectrum onto equal in log space freq axis
+  delta_f <- min(spec$freq)
+  logfreq <- log(spec$freq)
+
+  freq_logspace <- (seq(min(logfreq), max(logfreq)+delta_f, length.out = f.res*length(spec$freq)))
+  spec_loginterp <- stats::approx(logfreq, spec$spec, xout = freq_logspace, rule = 2)$y
+
+  spans_adj <- spans * f.res
+
+  # DOF of filter
+  kernel <- stats::kernel("daniell", spans_adj %/% 2)
+  filter <- kernel[-kernel$m:kernel$m]
+  df.kern <- stats::df.kernel(kernel)
+
+  # DOF of a boxcar filter the same width
+  kernal.flat <- stats::kernel("daniell", length(filter) %/% 2)
+  df.kern.flat <- stats::df.kernel(kernal.flat)
+
+  # modify for non boxcar filters
+  df.mod <- df.kern / df.kern.flat
+
+  # smooth/filter in log space
+  spec_filt <- ApplyFilter(spec_loginterp, filter = filter, method = method)
+
+  # re-interpolate back to original freq axis
+  spec3 <- stats::approx(freq_logspace, spec_filt, xout = logfreq)$y
+
+  # overwrite spec with filtered spec
+  spec$spec <- spec3
+
+  # keep old DOF
+  dof0 <- spec$dof
+
+  # Gets the difference in delta_f for the log and standard freq axis
+  NpF <- function(freq, fw, df){
+
+    posdiff <- (exp(log(freq) + df) - freq)
+    negdiff <- (freq - exp(log(freq) - df))
+
+    fdiff <- rowMeans(cbind(negdiff, posdiff))
+
+    2 * fw * (fdiff/df) * 1/(2*max(freq))
+  }
+  df.logkern <- NpF(spec$freq, length(filter), df = diff(freq_logspace[1:2]))
+
+  spec$dof <- spec$dof + df.mod * df.logkern * spec$dof/2
+  spec$shape <- spec$dof/2
+  spec$spans <- paste(spans, collapse = ",")
+
+  if (method == 0){
+    # remove NA portion
+    spec$freq <- spec$freq[is.na(spec$spec) == FALSE]
+    spec$dof <- spec$dof[is.na(spec$spec) == FALSE]
+    spec$shape <- spec$shape[is.na(spec$spec) == FALSE]
+    spec$spec <- spec$spec[is.na(spec$spec) == FALSE]
+  }
+
+  spec <- AddConfInterval(spec)
+
+
+
+  return(spec)
+}

R/SpecACF.R

@@ -217,24 +217,43 @@ mvacf.by.fft <- function(x){
 
 #' Remove leading and trailing rows of all NA
 #'
-#' @param m a numeric matrix
-#'
-#' @return a numeric matrix
-#' @keywords internal
+#' @param m a numeric matrix, data.frame or vector
+#' @param trim trim leading and trailing rows of "all" NA or containing "any" NA values
+#' @return a numeric matrix, data.frame or vector
+#' @export
 #' @examples
-#' m <- matrix(c(NA, NA, NA, 1:9, NA,NA,NA, 10:12, NA,NA,NA), ncol = 3, byrow = TRUE)
+#' m <- matrix(c(NA, NA, NA, 1, NA, NA, NA, 1, 1, NA, NA, NA, 1:9, NA,NA,NA, 10:12, NA, 1, NA, NA,NA,NA), ncol = 3, byrow = TRUE)
 #' m
-#' PaleoSpec:::TrimNA(m)
-TrimNA <- function(m){
+#' TrimNA(m)
+#' TrimNA(m, trim = "any")
+TrimNA <- function(m, trim = c("all", "any")) {
+  trim <- match.arg(trim)
+
+  # make it work on vectors
+  class_m <- class(m)
+  if (class_m[1] == "numeric") {
+    m <- cbind(m)
+  }
 
-  empty.row <- is.nan(rowMeans(m, na.rm = TRUE))
-  rank.good <- (empty.row == FALSE) * 1:length(empty.row)
+  if (trim == "all") {
+    empty.row <- is.nan(rowMeans(m, na.rm = TRUE))
+    rank.good <- (empty.row == FALSE) * 1:length(empty.row)
+  } else if (trim == "any") {
+    empty.row <- is.na(rowMeans(m))
+    rank.good <- (empty.row == FALSE) * 1:length(empty.row)
+  }
 
   first.good <- which.min(empty.row * 1:length(empty.row))
   last.good <- which.max(rank.good)
 
-  m[first.good:last.good, , drop = FALSE]
+  m <- m[first.good:last.good, , drop = FALSE]
+
+  # return to a vector
+  if (class_m[1] == "numeric") {
+    m <- as.numeric(m)
+  }
 
+  return(m)
 }
 
 

README.Rmd

@@ -27,8 +27,8 @@ PaleoSpec is an R package to assist the analysis of variance and power spectra.
 You can install the development version of PaleoSpec from [GitHub](https://github.com/) with:
 
 ``` r
-# install.packages("devtools")
-devtools::install_github("EarthSystemDiagnostics/paleospec")
+# install.packages("remotes")
+remotes::install_github("EarthSystemDiagnostics/paleospec")
 ```
 
 ## Usage

README.md

@@ -18,8 +18,8 @@ You can install the development version of PaleoSpec from
 [GitHub](https://github.com/) with:
 
 ``` r
-# install.packages("devtools")
-devtools::install_github("EarthSystemDiagnostics/paleospec")
+# install.packages("remotes")
+remotes::install_github("EarthSystemDiagnostics/paleospec")
 ```
 
 ## Usage