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

Hans W Borchers Venherm.Borchers at t-online.de
Fri Feb 6 23:48:24 CET 2004

```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.

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)
}
#-----------------------------------------------------------------------

```