diff --git a/DESCRIPTION b/DESCRIPTION
index 1e92777..b102dfd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: PaleoSpec
 Title: Spectral tools for the ECUS group
-Version: 0.31
+Version: 0.32
 Authors@R: c(
   person("Thomas", "Laepple", email = "tlaepple@awi.de", role = c("aut", "cre")),
   person("Thomas", "Muench", email = "tmuench@awi.de", role = c("aut")),
@@ -14,7 +14,7 @@ Depends:
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Imports: 
     multitaper,
     zoo,
diff --git a/NEWS.md b/NEWS.md
index 45945d6..f4c4208 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,10 @@
+# PaleoSpec 0.32
+
+* SpecACF can now use slepian tapers
+* Small fix to rfreq calculation in SpecACF
+* RemoveFirst and RemoveLast in gg_spec work on a per spectrum basis for list
+  and spec_df objects
+
 # PaleoSpec 0.31
 
 * Improve linking between function documentation using @family
diff --git a/R/SpecACF.R b/R/SpecACF.R
index d9dcf7f..dace8c0 100644
--- a/R/SpecACF.R
+++ b/R/SpecACF.R
@@ -1,7 +1,9 @@
-#' Estimate Power Spectra via the Autocovariance Function
+#' Estimate Power Spectra via the Autocovariance Function With Optional Slepian
+#' Tapers
 #'
 #' @param x a vector or matrix of binned values, possibly with gaps
-#' @param bin.width the width of the bins, effectively delta_t
+#' @param deltat,bin.width the time-step of the timeseries, equivalently the
+#' width of the bins in a binned timeseries, set only one
 #' @param pos.f.only return only positive frequencies, defaults to TRUE If TRUE,
 #'   freq == 0, and frequencies higher than 1/(2*bin.width) which correspond to
 #'   the negative frequencies are removed
@@ -20,6 +22,8 @@
 #'   indicates that the estimates are unbiased such that smoothing across
 #'   frequencies to remove negative estimates results in an unbiased power
 #'   spectrum.
+#' @inheritParams SpecMTM
+#' @importFrom multitaper spec.mtm
 #' @author Torben Kunz and Andrew Dolman <andrew.dolman@awi.de>
 #' @return a spec object (list)
 #' @family functions to estimate power spectra
@@ -27,92 +31,221 @@
 #'
 #' @examples
 #' set.seed(20230312)
-#' x <- cumsum(rgamma(200, shape = 1.5, rate = 1.5/10))
-#' y <- SimProxySeries(a = 0.1, b = 1, t.smpl = x, nt = 2000,
-#'  smth.lab = list(type = "rect", tau = 1))
-#' y_binned <- BinTimeseries(x, y, bin.width = 15)
-#' sp1 <- SpecACF(y_binned$mean.value, bin.width = 15)
-#' sp2 <- LogSmooth(sp1)
-#' LPlot(sp1)
-#' LLines(sp2, col = "red")
-SpecACF <- function(x, bin.width,
-                    demean = TRUE, detrend = TRUE,
-                    TrimNA = TRUE,
-                    pos.f.only = TRUE,
-                    return.working = FALSE){
-
-   # Convert ts to vector
-
-  if (is.ts(x)){
+#'
+#' # Comparison with SpecMTM
+#'
+#' tsM <- replicate(2, SimPLS(1e03, 1, 0.1))
+#' spMk3 <- SpecACF(tsM, bin.width = 1, k = 3, nw = 2)
+#' spMk1 <- SpecACF(tsM, bin.width = 1, k = 1, nw = 0)
+#'
+#' spMTMa <- SpecMTM(tsM[,1], deltat = 1)
+#' spMTMb <- SpecMTM(tsM[,2], deltat = 1)
+#' spMTM <- spMTMa
+#' spMTM$spec <- (spMTMa$spec + spMTMb$spec)/2
+#'
+#' gg_spec(list(
+#'   `ACF k=1` = spMk1,
+#'   `ACF k=3` = spMk3,
+#'   `MTM k=3` = spMTM
+#' ), alpha.line = 0.75) +
+#'   ggplot2::facet_wrap(~spec_id)
+#'
+#' ## No gaps
+#'
+#' ts1 <- SimPLS(1000, 1, 0.1)
+#'
+#' sp_ACF1 <- SpecACF(ts1, 1, k = 1)
+#' sp_MTM7 <- SpecMTM(ts1, nw = 4, k = 7, deltat = 1)
+#' sp_ACF7 <- SpecACF(ts1, 1, k = 7, nw = 4)
+#'
+#' gg_spec(list(
+#'   `ACF k=1` = sp_ACF1, `ACF k=7` = sp_ACF7, `MTM k=7` = sp_MTM7
+#' ))
+#'
+#' # With Gaps
+#'
+#' gaps <- (arima.sim(list(ar = 0.5), n = length(ts1))) > 1
+#' table(gaps)
+#' ts1_g <- ts1
+#' ts1_g[gaps] <- NA
+#'
+#' sp_ACF1_g <- SpecACF(ts1_g, 1)
+#' sp_ACFMTM1_g <- SpecACF(ts1_g, bin.width = 1, nw = 4, k = 7)
+#'
+#' gg_spec(list(
+#'   ACF_g = sp_ACF1_g,
+#'   ACF_g_smoothed = FilterSpecLog(sp_ACF1_g),
+#'   ACF_g_tapered = sp_ACFMTM1_g
+#' ), conf = FALSE) +
+#'   ggplot2::geom_abline(intercept = log10(0.1), slope = -1, lty = 2)
+#'
+#'
+#'
+#' ## AR4
+#' arc_spring <- c(2.7607, -3.8106, 2.6535, -0.9238)
+#'
+#' tsAR4 <- arima.sim(list(ar = arc_spring), n = 1e03) + rnorm(1e03, 0, 10)
+#' plot(tsAR4)
+#' spAR4_ACF <- SpecACF(tsAR4, 1)
+#' spAR4_MTACF <- SpecACF(as.numeric(tsAR4), 1, k = 15, nw = 8)
+#'
+#' gg_spec(list(#'
+#'   `ACF k=1` = spAR4_ACF,
+#'   `ACF k=15` = spAR4_MTACF)
+#' )
+#'
+#' ## Add gaps to timeseries
+#'
+#' gaps <- (arima.sim(list(ar = 0.5), n = length(tsAR4))) > 2
+#' table(gaps)
+#' tsAR4_g <- tsAR4
+#' tsAR4_g[gaps] <- NA
+#'
+#' plot(tsAR4, col = "green")
+#' lines(tsAR4_g, col = "blue")
+#'
+#' table(tsAR4_g > 0, useNA = "always")
+#'
+#' spAR4_ACF_g <- SpecACF(as.numeric(tsAR4_g), 1)
+#' spAR4_MTACF_g <- SpecACF(as.numeric(tsAR4_g), 1, nw = 8, k = 15)
+#'
+#' table(spAR4_ACF_g$spec < 0)
+#' table(spAR4_MTACF_g$spec < 0)
+#'
+#' gg_spec(list(
+#'   `ACF gaps k=1` = spAR4_ACF_g,
+#'   `ACF gaps k = 15` = spAR4_MTACF_g,
+#'   `ACF full k = 15` = spAR4_MTACF
+#' )
+#' )
+SpecACF <- function(x,
+                      deltat = NULL, bin.width = NULL,
+                      k = 3, nw = 2,
+                      demean = TRUE, detrend = TRUE,
+                      TrimNA = TRUE,
+                      pos.f.only = TRUE,
+                      return.working = FALSE) {
+  if (is.null(deltat) & is.null(bin.width) & is.ts(x) == FALSE) {
+    stop("One of deltat or bin.width must be set")
+  }
+
+  # Convert ts to vector
+  if (is.ts(x)) {
     d <- dim(x)
+    dt_ts <- deltat(x)
+
     x <- as.vector(x)
 
-    if (is.null(d) == FALSE){
+    if (is.null(deltat) == FALSE) {
+      if (dt_ts != deltat) stop("timeseries deltat does not match argument deltat")
+    }
+
+    if (is.null(bin.width) == FALSE) {
+      if (dt_ts != bin.width) stop("timeseries deltat does not match argument bin.width")
+    }
+
+    if (is.null(deltat)) deltat <- dt_ts
+    if (is.null(bin.width)) bin.width <- dt_ts
+
+    if (is.null(d) == FALSE) {
       dim(x) <- d
     }
+  }
 
+  if (is.null(bin.width) & is.null(deltat) == FALSE) {
+    bin.width <- deltat
   }
 
+  if (is.null(bin.width) == FALSE & is.null(deltat)) {
+    deltat <- bin.width
+  }
 
   # Convert vector to matrix
-  if (is.vector(x)){
+  if (is.vector(x)) {
     x <- matrix(x, ncol = 1)
   }
 
-
-
-  if (is.data.frame(x) == TRUE){
+  if (is.data.frame(x) == TRUE) {
     x <- as.matrix(x)
   }
 
-  if (TrimNA){
+  if (TrimNA) {
     x <- TrimNA(x)
   }
 
-  if (detrend){
-    i <- seq_along(x[,1])
+  if (detrend) {
+    i <- seq_along(x[, 1])
     x <- apply(x, 2, function(y) {
       stats::residuals(stats::lm(y ~ i, na.action = "na.exclude"))
     })
   }
 
   # remove mean from each record
-  if (demean){
+  if (demean) {
     x <- x - colMeans(x, na.rm = TRUE)
   }
 
-  lag = (0:(nrow(x)-1))
+  lag <- (0:(nrow(x) - 1))
 
-  if (return.working == TRUE){
-    working <- mean.ACF(x, return.working = TRUE)
+  ncolx <- ncol(x)
+
+  ## Tapering
+  if (k > 1) {
+    n <- nrow(x)
+
+    dpssIN <- multitaper:::dpss(n,
+                                k = k, nw = nw,
+                                returnEigenvalues = TRUE
+    )
+    dw <- dpssIN$v #* sqrt(bin.width)
+    ev <- dpssIN$eigen
+
+    x <- matrix(unlist(lapply(1:ncol(x), function(i) {
+      lapply(1:ncol(dw), function(j) {
+        x[, i] * dw[, j]
+      })
+    })), nrow = nrow(x), byrow = FALSE)
+  }
+
+  if (return.working == TRUE) {
+    working <- PaleoSpec:::mean.ACF(x, return.working = TRUE)
     acf <- working[["acf"]]
-  } else{
-    acf <- mean.ACF(x)
+  } else {
+    acf <- PaleoSpec:::mean.ACF(x)
+  }
+
+  freq <- (lag / (nrow(x))) / bin.width
+
+  if (k > 1){
+    spec <- Re(stats::fft(acf)) * bin.width * nrow(x)
+  } else {
+    spec <- Re(stats::fft(acf)) * bin.width
   }
 
-  freq = (lag / (nrow(x))) / bin.width
-  spec = Re(stats::fft(acf)) * bin.width
 
-  rfreq <- (utils::tail(freq, 1) + freq[2])/2
 
-  if (pos.f.only){
-    spec <- spec[freq > 0 & freq <= 1/(2*bin.width)]
-    freq <- freq[freq > 0 & freq <= 1/(2*bin.width)]
+
+  rfreq <- max(freq)
+
+  if (pos.f.only) {
+    spec <- spec[freq > 0 & freq <= 1 / (2 * bin.width)]
+    freq <- freq[freq > 0 & freq <= 1 / (2 * bin.width)]
   }
 
-  out <- list(bin.width = bin.width,
-              rfreq = rfreq,
-              nrec = ncol(x),
-              lag = lag, acf = acf,
-              freq = freq, spec = spec, f.length = 1,
-              dof = rep(2*ncol(x), length(freq)))
+  out <- list(
+    bin.width = bin.width,
+    rfreq = rfreq,
+    nrec = ncolx,
+    lag = lag, acf = acf,
+    freq = freq, spec = spec, f.length = 1,
+    dof = rep(2 * ncolx * k, length(freq))
+  )
 
-  class(out) <-  c("SpecACF", "spec")
+  class(out) <- c("SpecACF", "spec")
 
-  if (return.working){
+  if (return.working) {
     out <- list(working = working, spec = out)
   }
-
   return(out)
 }
 
diff --git a/R/gg_spec.R b/R/gg_spec.R
index 357805d..e674410 100644
--- a/R/gg_spec.R
+++ b/R/gg_spec.R
@@ -86,24 +86,42 @@ gg_spec <- function(x, gg = NULL, conf = TRUE,
   }
 
 
-  if (removeFirst > 0) {
-    df <- df[rank(df$freq) > removeFirst, ]
-  }
-
-  if (removeLast > 0) {
-    df <- df[rank(-df$freq) > removeLast,]
-  }
+  ## Remove first and or last freqs on a per spec_id basis
 
   if (exists("spec_id", df) == FALSE){
     df$spec_id <- 1
-    }
-
+  }
 
   if (is.numeric(df$spec_id)){
     df$spec_id <- as.character(df$spec_id)
   }
 
 
+  if (removeFirst > 0) {
+    #df <- df[rank(df$freq) > removeFirst, ]
+
+    # Convert 'df' to a data.table
+    data.table::setDT(df)
+
+    # Group by 'spec_id' and filter rows where 'rank(freq)' is greater than 'removeFirst'
+    df <- df[, .SD[rank(freq) > removeFirst], by = spec_id]
+
+    # return to data.frame
+    df <- as.data.frame(df)
+  }
+
+  if (removeLast > 0) {
+    #df <- df[rank(-df$freq) > removeLast,]
+    # Convert 'df' to a data.table
+    data.table::setDT(df)
+
+    # Group by 'spec_id' and filter rows where 'rank(-freq)' is greater than 'removeLast'
+    df <- df[, .SD[rank(-freq) > removeLast], by = spec_id]
+
+    # return to data.frame
+    df <- as.data.frame(df)
+  }
+
 
   if (is.null(gg)) {
     p <- ggplot2::ggplot(data = df)
diff --git a/README.Rmd b/README.Rmd
index 9b30e35..b26b8a2 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -84,7 +84,7 @@ where: $S(f) = \alpha f^{-\beta}$
 set.seed(20221109)
 
 # length of the time series
-N <- 1e04
+N <- 1e03
 
 # parameters of the powerlaw spectrum
 alpha <- 0.1
diff --git a/README.md b/README.md
index 69775bd..a712c0e 100644
--- a/README.md
+++ b/README.md
@@ -88,7 +88,7 @@ powerlaw like properties, where: $S(f) = \alpha f^{-\beta}$
 set.seed(20221109)
 
 # length of the time series
-N <- 1e04
+N <- 1e03
 
 # parameters of the powerlaw spectrum
 alpha <- 0.1
diff --git a/man/LogSmooth.Rd b/man/LogSmooth.Rd
index 3ef596b..c51539b 100644
--- a/man/LogSmooth.Rd
+++ b/man/LogSmooth.Rd
@@ -51,8 +51,8 @@ LLines(LogSmooth(spec,df.log=0.1, removeFirst = 3, removeLast = 20),lwd=2,col='r
 }
 \seealso{
 Other functions to filter / smooth spectra: 
-\code{\link{FilterSpecLog}()},
-\code{\link{FilterSpec}()}
+\code{\link{FilterSpec}()},
+\code{\link{FilterSpecLog}()}
 }
 \author{
 Thomas Laepple
diff --git a/man/SpecACF.Rd b/man/SpecACF.Rd
index 845ea60..a26a491 100644
--- a/man/SpecACF.Rd
+++ b/man/SpecACF.Rd
@@ -2,11 +2,15 @@
 % Please edit documentation in R/SpecACF.R
 \name{SpecACF}
 \alias{SpecACF}
-\title{Estimate Power Spectra via the Autocovariance Function}
+\title{Estimate Power Spectra via the Autocovariance Function With Optional Slepian
+Tapers}
 \usage{
 SpecACF(
   x,
-  bin.width,
+  deltat = NULL,
+  bin.width = NULL,
+  k = 3,
+  nw = 2,
   demean = TRUE,
   detrend = TRUE,
   TrimNA = TRUE,
@@ -17,7 +21,13 @@ SpecACF(
 \arguments{
 \item{x}{a vector or matrix of binned values, possibly with gaps}
 
-\item{bin.width}{the width of the bins, effectively delta_t}
+\item{deltat, bin.width}{the time-step of the timeseries, equivalently the
+width of the bins in a binned timeseries, set only one}
+
+\item{k}{a positive integer, the number of tapers, often 2*nw.}
+
+\item{nw}{a positive double precision number, the time-bandwidth
+parameter.}
 
 \item{demean}{remove the mean from each record (column) in x, defaults to
 TRUE. If detrend is TRUE, mean will be removed during detrending regardless
@@ -47,14 +57,93 @@ Estimates the power spectrum from a single time series, or the
 }
 \examples{
 set.seed(20230312)
-x <- cumsum(rgamma(200, shape = 1.5, rate = 1.5/10))
-y <- SimProxySeries(a = 0.1, b = 1, t.smpl = x, nt = 2000,
- smth.lab = list(type = "rect", tau = 1))
-y_binned <- BinTimeseries(x, y, bin.width = 15)
-sp1 <- SpecACF(y_binned$mean.value, bin.width = 15)
-sp2 <- LogSmooth(sp1)
-LPlot(sp1)
-LLines(sp2, col = "red")
+
+# Comparison with SpecMTM
+
+tsM <- replicate(2, SimPLS(1e03, 1, 0.1))
+spMk3 <- SpecACF(tsM, bin.width = 1, k = 3, nw = 2)
+spMk1 <- SpecACF(tsM, bin.width = 1, k = 1, nw = 0)
+
+spMTMa <- SpecMTM(tsM[,1], deltat = 1)
+spMTMb <- SpecMTM(tsM[,2], deltat = 1)
+spMTM <- spMTMa
+spMTM$spec <- (spMTMa$spec + spMTMb$spec)/2
+
+gg_spec(list(
+  `ACF k=1` = spMk1,
+  `ACF k=3` = spMk3,
+  `MTM k=3` = spMTM
+), alpha.line = 0.75) +
+  ggplot2::facet_wrap(~spec_id)
+
+## No gaps
+
+ts1 <- SimPLS(1000, 1, 0.1)
+
+sp_ACF1 <- SpecACF(ts1, 1, k = 1)
+sp_MTM7 <- SpecMTM(ts1, nw = 4, k = 7, deltat = 1)
+sp_ACF7 <- SpecACF(ts1, 1, k = 7, nw = 4)
+
+gg_spec(list(
+  `ACF k=1` = sp_ACF1, `ACF k=7` = sp_ACF7, `MTM k=7` = sp_MTM7
+))
+
+# With Gaps
+
+gaps <- (arima.sim(list(ar = 0.5), n = length(ts1))) > 1
+table(gaps)
+ts1_g <- ts1
+ts1_g[gaps] <- NA
+
+sp_ACF1_g <- SpecACF(ts1_g, 1)
+sp_ACFMTM1_g <- SpecACF(ts1_g, bin.width = 1, nw = 4, k = 7)
+
+gg_spec(list(
+  ACF_g = sp_ACF1_g,
+  ACF_g_smoothed = FilterSpecLog(sp_ACF1_g),
+  ACF_g_tapered = sp_ACFMTM1_g
+), conf = FALSE) +
+  ggplot2::geom_abline(intercept = log10(0.1), slope = -1, lty = 2)
+
+
+
+## AR4
+arc_spring <- c(2.7607, -3.8106, 2.6535, -0.9238)
+
+tsAR4 <- arima.sim(list(ar = arc_spring), n = 1e03) + rnorm(1e03, 0, 10)
+plot(tsAR4)
+spAR4_ACF <- SpecACF(tsAR4, 1)
+spAR4_MTACF <- SpecACF(as.numeric(tsAR4), 1, k = 15, nw = 8)
+
+gg_spec(list(#'
+  `ACF k=1` = spAR4_ACF,
+  `ACF k=15` = spAR4_MTACF)
+)
+
+## Add gaps to timeseries
+
+gaps <- (arima.sim(list(ar = 0.5), n = length(tsAR4))) > 2
+table(gaps)
+tsAR4_g <- tsAR4
+tsAR4_g[gaps] <- NA
+
+plot(tsAR4, col = "green")
+lines(tsAR4_g, col = "blue")
+
+table(tsAR4_g > 0, useNA = "always")
+
+spAR4_ACF_g <- SpecACF(as.numeric(tsAR4_g), 1)
+spAR4_MTACF_g <- SpecACF(as.numeric(tsAR4_g), 1, nw = 8, k = 15)
+
+table(spAR4_ACF_g$spec < 0)
+table(spAR4_MTACF_g$spec < 0)
+
+gg_spec(list(
+  `ACF gaps k=1` = spAR4_ACF_g,
+  `ACF gaps k = 15` = spAR4_MTACF_g,
+  `ACF full k = 15` = spAR4_MTACF
+)
+)
 }
 \seealso{
 Other functions to estimate power spectra: 
diff --git a/man/figures/README-unnamed-chunk-10-1.png b/man/figures/README-unnamed-chunk-10-1.png
index 75cd0a9..4561da1 100644
Binary files a/man/figures/README-unnamed-chunk-10-1.png and b/man/figures/README-unnamed-chunk-10-1.png differ
diff --git a/man/figures/README-unnamed-chunk-3-1.png b/man/figures/README-unnamed-chunk-3-1.png
index f53dcf5..3e8a96a 100644
Binary files a/man/figures/README-unnamed-chunk-3-1.png and b/man/figures/README-unnamed-chunk-3-1.png differ
diff --git a/man/figures/README-unnamed-chunk-4-1.png b/man/figures/README-unnamed-chunk-4-1.png
index 46c9e7a..0d098e2 100644
Binary files a/man/figures/README-unnamed-chunk-4-1.png and b/man/figures/README-unnamed-chunk-4-1.png differ
diff --git a/man/figures/README-unnamed-chunk-5-1.png b/man/figures/README-unnamed-chunk-5-1.png
index 050c59d..bd6a23e 100644
Binary files a/man/figures/README-unnamed-chunk-5-1.png and b/man/figures/README-unnamed-chunk-5-1.png differ
diff --git a/man/figures/README-unnamed-chunk-6-1.png b/man/figures/README-unnamed-chunk-6-1.png
index 5eaca85..79f9046 100644
Binary files a/man/figures/README-unnamed-chunk-6-1.png and b/man/figures/README-unnamed-chunk-6-1.png differ
diff --git a/man/figures/README-unnamed-chunk-7-1.png b/man/figures/README-unnamed-chunk-7-1.png
index 0cec6ea..5cbaad8 100644
Binary files a/man/figures/README-unnamed-chunk-7-1.png and b/man/figures/README-unnamed-chunk-7-1.png differ
diff --git a/man/figures/README-unnamed-chunk-8-1.png b/man/figures/README-unnamed-chunk-8-1.png
index c3d603f..331d386 100644
Binary files a/man/figures/README-unnamed-chunk-8-1.png and b/man/figures/README-unnamed-chunk-8-1.png differ
diff --git a/man/figures/README-unnamed-chunk-9-1.png b/man/figures/README-unnamed-chunk-9-1.png
index 1540779..70409b2 100644
Binary files a/man/figures/README-unnamed-chunk-9-1.png and b/man/figures/README-unnamed-chunk-9-1.png differ
diff --git a/tests/testthat/test-FilterSpec.R b/tests/testthat/test-FilterSpec.R
index 3bef6a2..8d910dd 100644
--- a/tests/testthat/test-FilterSpec.R
+++ b/tests/testthat/test-FilterSpec.R
@@ -148,7 +148,7 @@ test_that("FilterSpecLog handles df.log parameter correctly", {
 test_that("FilterSpecLog handles edge cases appropriately", {
   empty_spec <- list(freq = numeric(0), spec = numeric(0), dof = numeric(0), shape = numeric(0))
 
-  expect_error(FilterSpecLog(empty_spec))
+  expect_warning(expect_error(FilterSpecLog(empty_spec)))
   expect_error(FilterSpecLog("not a spec"))
 })
 
diff --git a/tests/testthat/test-SpecACF.R b/tests/testthat/test-SpecACF.R
index 86dc681..40bdcf7 100644
--- a/tests/testthat/test-SpecACF.R
+++ b/tests/testthat/test-SpecACF.R
@@ -13,9 +13,20 @@ test_that("SpecACF", {
   expect_length(sp1$spec, n/2)
   expect_length(sp1$dof, n/2)
 
+  # test deltat argument
+  ts2 <- ts(rnorm(n), deltat = 3)
+  sp1_bw <- SpecACF(ts2, bin.width = 3)
+  sp1_dt <- SpecACF(ts2, deltat = 3)
+  expect_s3_class(sp1_dt, "spec")
+  expect_equal(sp1_bw$freq, sp1_dt$freq)
+
+  # test using tapers
+  sp1_tap <- SpecACF(ts1, deltat = 1, k = 3, nw = 2)
+  expect_s3_class(sp1_tap, "spec")
+  expect_equal(sp1_tap$dof, rep(3*2, length(sp1_tap$freq)))
 
   # test for neg freq
-  sp2 <- SpecACF(ts1, bin.width = 1, pos.f.only = FALSE)
+  sp2 <- SpecACF(ts1, bin.width = 1, k = 1, pos.f.only = FALSE)
 
   expect_s3_class(sp2, "spec")
 
@@ -28,7 +39,7 @@ test_that("SpecACF", {
 
   m <- matrix(rnorm(3*100), ncol = 3)
 
-  spm <- SpecACF(m, bin.width = 1)
+  spm <- SpecACF(m, bin.width = 1, k = 1)
 
   #LPlot(spm)
 
@@ -45,7 +56,7 @@ test_that("SpecACF", {
 
 
   # do not demean
-  sp1ndm <- SpecACF(rnorm(100, mean = 100), bin.width = 1, demean = FALSE,
+  sp1ndm <- SpecACF(rnorm(100, mean = 100), bin.width = 1, k = 1, demean = FALSE,
                     detrend = FALSE, return.working = TRUE, pos.f.only = FALSE)
   expect_true(mean(sp1ndm$working$x) > 90)
   expect_true(sp1ndm$spec[1] > sp2$spec[1])
diff --git a/tests/testthat/test-gg_spec.R b/tests/testthat/test-gg_spec.R
index 3c35d3a..efe38d8 100644
--- a/tests/testthat/test-gg_spec.R
+++ b/tests/testthat/test-gg_spec.R
@@ -49,3 +49,18 @@ test_that("plots have correct no colours", {
   expect_equal(nrow(unique(gg1_df2_b$data[[1]]["colour"])), 2)
 
 })
+
+
+test_that("removeFirst and removeLast operate on a per spec_id basis", {
+  sp1 <- SpecMTM(ts(rnorm(100)))
+  sp2 <- SpecMTM(ts(rnorm(1000)))
+
+  gf1 <- gg_spec(list(sp1, sp2), removeFirst = 1)
+
+  g0 <- gg_spec(list(sp1, sp2), removeFirst = 0)
+
+  gl3 <- gg_spec(list(sp1, sp2), removeLast = 3)
+
+  expect_that(nrow(g0$data) - nrow(gl3$data), equals(6))
+
+})