[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.
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)
}
#-----------------------------------------------------------------------
More information about the R-help
mailing list