Apply a least-squares (Savitzky-Golay) polynomial smoothing filter
Reduce noise in 1-D data (e.g. time-series, spectrum) but retain
dynamic range of variations in the data by applying a least squares
smoothing polynomial filter,
Also called the Savitzky-Golay smoothing filter, cf. Numerical
Recipes (Press et al. 1992, Sec.14.8)
The low-pass filter coefficients are computed by effectively
least-squares fitting a polynomial in moving window,
centered on each data point, so the new value will be the
zero-th coefficient of the polynomial. Approximate first derivates
of the data can be computed by using first degree coefficient of
each polynomial, and so on. The filter coefficients for a specified
polynomial degree and window width are computed independent of any
data, and stored in a common block. The filter is then convolved
with the data array to result in smoothed data with reduced noise,
but retaining higher order variations (better than SMOOTH).
This procedure became partially obsolete in IDL V5.4 with the
introduction of the SAVGOL function, which computes the smoothing
spectrum = poly_smooth( data, [ width, DEGREE = , NLEFT = , NRIGHT =
DERIV_ORDER = ,COEFF = ]
data = 1-D array, such as a spectrum or time-series.
width = total number of data points to use in filter convolution,
(default = 5, using 2 past and 2 future data points),
must be larger than DEGREE of polynomials, and a guideline is to
make WIDTH between 1 and 2 times the FWHM of desired features.
Optional Input Keywords
DEGREE = degree of polynomials to use in designing the filter
via least squares fits, (default DEGREE = 2)
The higher degrees will preserve sharper features.
NLEFT = # of past data points to use in filter convolution,
excluding current point, overrides width parameter,
so that width = NLEFT + NRIGHT + 1. (default = NRIGHT)
NRIGHT = # of future data points to use (default = NLEFT).
DERIV_ORDER = order of derivative desired (default = 0, no derivative).
Optional Output Keyword
COEFFICIENTS = optional output of the filter coefficients applied,
but they are all stored in common block for reuse, anyway.
Function returns the data convolved with polynomial filter coefs.
Given a wavelength - flux spectrum (w,f), apply a 31 point quadratic
smoothing filter and plot
IDL> cgplot, w, poly_smooth(f,31)
common poly_smooth, degc, nlc, nrc, coefs, ordermax
As described in Numerical Recipies, 2nd edition sec.14.8,
Matrix of normal eqs. is formed by starting with small terms
and then adding progressively larger terms (powers).
The filter coefficients of up to derivative ordermax are stored
in common, until the specifications change, then recompute coefficients.
Coefficients are stored in convolution order, zero lag in the middle.
Written, Frank Varosi NASA/GSFC 1993.
Converted to IDL V5.0 W. Landsman September 1997
Use /EDGE_TRUNCATE keyword to CONVOL W. Landsman March 2006