R/smooth.R

In MsCoreUtils: Core Utils for Mass Spectrometry Data

#' @title Smoothing
#'
#' @description
#' This function smoothes a numeric vector.
#'
#' @param x `numeric`, i.e. m/z values.
#' @param cf `matrix`, a coefficient matrix generated by `coefMA`, `coefWMA` or
#' `coefSG`.
#'
#' @return `smooth`: A `numeric` of the same length as `x`.
#'
#' @author Sebastian Gibb, Sigurdur Smarason (weighted moving average)
#' @aliases smooth
#' @family noise estimation and smoothing functions
#' @importFrom stats filter
#' @export
#' @examples
#' x <- c(1:10, 9:1)
#' plot(x, type = "b", pch = 20)
#' cf <- list(MovingAverage = coefMA(2),
#'            WeightedMovingAverage = coefWMA(2),
#'            SavitzkyGolay = coefSG(2))
#' for (i in seq_along(cf)) {
#'     lines(smooth(x, cf[[i]]), col = i + 1, pch = 20, type = "b")
#' }
#' legend("bottom", legend = c("x", names(cf)), pch = 20,
#'        col = seq_len(length(cf) + 1))
smooth <- function(x, cf) {
    d <- dim(cf)
    if (!is.matrix(cf) || d[1L] != d[2L] || d[1L] < 3)
        stop("'cf' has to be matrix with equal number of rows and colums.")
    n <- length(x)
    w <- dim(cf)[1L]
    .validateWindow(w, n)
    hws <- trunc(w / 2L)
    y <- filter(x = x, filter = cf[hws + 1L, ], sides = 2L)
    attributes(y) <- NULL

    ## Implemention of left/right extrema based on:
    ## sgolay in signal 0.7-3/R/sgolay.R by Paul Kienzle <pkienzle@users.sf.net>
    ## modified by Sebastian Gibb <mail@sebastiangibb.de>
    ## fix left/right extrema
    lhws <- seq_len(hws)
    y[lhws] <- cf[lhws, , drop = FALSE] %*% x[seq_len(w)]
    y[seq.int(to = n, length.out = hws)] <-
        cf[seq.int(to = w, length.out = hws), , drop = FALSE] %*%
        x[seq.int(to = n, length.out = w)]
    y
}

#' @describeIn smooth Simple Moving Average
#'
#' This function calculates the coefficients for a simple moving average.
#'
#' @note
#' The `hws` depends on the used method ((weighted) moving
#' average/Savitzky-Golay).
#'
#' @param hws `integer(1)`, half window size, the resulting window reaches from
#' `(i - hws):(i + hws)`.
#'
#' @return `coefMA`: A `matrix` with coefficients for a simple moving average.
#' @export
coefMA <- function(hws) {
    w <- 2L * hws + 1L
    matrix(1L / w, nrow = w, ncol = w)
}

#' @describeIn smooth Weighted Moving Average
#'
#' This function calculates the coefficients for a weighted moving average with
#' weights depending on the distance from the center calculated as
#' `1/2^abs(-hws:hws)` with the sum of all weigths normalized to 1.
#'
#' @return `coefWMA`: A `matrix` with coefficients for a weighted moving average.
#' @export
coefWMA <- function(hws) {
    k <- 1 / 2^abs(-hws:hws)
    w <- 2L * hws + 1L
    matrix(k / sum(k), nrow = w, ncol = w, byrow = TRUE)
}

#' @describeIn smooth Savitzky-Golay-Filter
#'
#' This function calculates the Savitzky-Golay-Coefficients. The additional
#' argument `k` controls the order of the used polynomial. If `k` is set to zero
#' it yield a simple moving average.
#'
#' @param k `integer(1)`, set the order of the polynomial used to calculate the
#' coefficients.
#'
#' @details
#' For the Savitzky-Golay-Filter the `hws` should be smaller than
#' *FWHM* of the peaks (full width at half maximum; please find details in
#' Bromba and Ziegler 1981).
#'
#' In general the `hws` for the (weighted) moving average (`coefMA`/`coefWMA`)
#' has to bemuch smaller than for the Savitzky-Golay-Filter to conserve the
#' peak shape.
#'
#' @return `coefSG`: A `matrix` with *Savitzky-Golay-Filter* coefficients.
#'
#' @references
#' A. Savitzky and M. J. Golay. 1964.
#' Smoothing and differentiation of data by simplified least squares procedures.
#' Analytical chemistry, 36(8), 1627-1639.
#'
#' M. U. Bromba and H. Ziegler. 1981.
#' Application hints for Savitzky-Golay digital smoothing filters.
#' Analytical Chemistry, 53(11), 1583-1586.
#'
#' Implementation based on:
#' Steinier, J., Termonia, Y., & Deltour, J. (1972). Comments on Smoothing and
#' differentiation of data by simplified least square procedure.
#' Analytical Chemistry, 44(11), 1906-1909.
#'
#' @export
coefSG <- function(hws, k = 3L) {
    if (length(k) != 1L || !is.integer(k) || k < 0L)
        stop("'k' has to be an integer of length 1 and larger 0.")

    nk <- k + 1L
    k <- seq_len(nk) - 1L
    w <- 2L * hws + 1L

    if (w < nk)
        stop("The window size has to be larger than the polynomial order.")

    K <- matrix(k, nrow = w, ncol = nk, byrow = TRUE)

    ## filter is applied to -hws:hws around current data point
    ## to avoid removing (NA) of left/right extrema
    ## lhs: 0:(2 * hws)
    ## rhs: (n - 2 * hws):n

    ## filter matrix contains 2 * hws + 1 rows
    ## row 1:hws == lhs coef
    ## row hws + 1 == typical sg coef
    ## row (n - hws - 1):n == rhs coef
    F <- matrix(NA_real_, nrow = w, ncol = w)
    for (i in seq_len(hws + 1L)) {
        M <- matrix(seq_len(w) - i, nrow = w, ncol = nk, byrow = FALSE)
        X <- M^K
        F[i, ] <- (solve(t(X) %*% X) %*% t(X))[1L, ]
    }
    ## rhs (row (n-m):n) are equal to reversed lhs
    F[seq.int(to = w, length.out = hws), ] <- rev(F[seq_len(hws), ])
    F
}