[R] Savitzky-Golay smoothing -- an R implementation

[Hans W Borchers]

As the request for the Savitzky-Golay Algorithm in R has come up several
times, I here include my implementation based on code written for Matlab.
Savitzky-Golay uses the pseudo-inverse pinv() of a matrix. There is an
'generalized inverse' ginv() in the MASS package, but I use a simpler form
because I didn't want to 'require' MASS any time I apply Savitzky-Golay.

Savitzky-Golay is not only a good method for chemical engineering, it can
successfully be applied to smooth process data. One approach is to determine
the noise level in a time series (ACF, winGamma, ...) and then choose the
parameter fl such that the difference between the time series and its
Savitzky-Golay approximation reflects the noise level.

I would be glad to hear about comments and improvements.

Hans W. Borchers
ABB Corporate Research

P. S.:   Example:

    t  <- sin(2*pi*(1:1000)/200)
    t1 <- t + rnorm(1000)/10
    t2 <- sav.gol(t1, 51)
    plot(1:1000, t1)
    lines(1:1000, t,  col="blue")
    lines(1:1000, t2, col="red")

# ----------------------------------------------------------------------
#   Savitzky-Golay Algorithm
# ----------------------------------------------------------------------
# T2 <- sav.gol(T, fl, forder=4, dorder=0);
#
# Polynomial filtering method of Savitzky and Golay
# See Numerical Recipes, 1992, Chapter 14.8, for details.
#
# T      = vector of signals to be filtered
#          (the derivative is calculated for each ROW)
# fl     = filter length (for instance fl = 51..151)
# forder = filter order (2 = quadratic filter, 4= quartic)
# dorder = derivative order (0 = smoothing, 1 = first derivative, etc.)
#
sav.gol <- function(T, fl, forder=4, dorder=0)
{
    m <- length(T)
    dorder <- dorder + 1

    # -- calculate filter coefficients --
    fc <- (fl-1)/2                          # index: window left and right
    X  <- outer(-fc:fc, 0:forder, FUN="^")  # polynomial terms and 
coefficients
    Y  <- pinv(X);                          # pseudoinverse

    # -- filter via convolution and take care of the end points --
    T2 <- convolve(T, rev(Y[dorder,]), type="o")    # convolve(...)
    T2 <- T2[(fc+1):(length(T2)-fc)]
}
#-----------------------------------------------------------------------
#   *** PseudoInvers of a Matrix ***
#   using singular value decomposition
#
pinv <- function (A)
{
    s <- svd(A)
    # D <- diag(s$d); Dinv <- diag(1/s$d)
    # U <- s$u; V <- s$v
    # A = U D V'
    # X = V Dinv U'
    s$v %*% diag(1/s$d) %*% t(s$u)
}
#-----------------------------------------------------------------------