SAVGOL

The SAVGOL function returns the coefficients of a Savitzky-Golay smoothing filter, which can then be applied using the CONVOL function. The Savitzky-Golay smoothing filter, also known as least squares or DISPO (digital smoothing polynomial), can be used to smooth a noisy signal.

The filter is defined as a weighted moving average with weighting given as a polynomial of a certain degree. The returned coefficients, when applied to a signal, perform a polynomial least-squares fit within the filter window. This polynomial is designed to preserve higher moments within the data and reduce the bias introduced by the filter. The filter can use any number of points for this weighted average.

This filter works especially well when the typical peaks of the signal are narrow. The heights and widths of the curves are generally preserved.

This routine is written in the IDL language. Its source code can be found in the file savgol.pro in the lib subdirectory of the IDL distribution.

SAVGOL is based on the Savitzky-Golay Smoothing Filters described in section 14.8 of Numerical Recipes in C: The Art of Scientific Computing (Second Edition), published by Cambridge University Press, and is used by permission.

Syntax

Return Value

This function returns an array of floating-point numbers that are the coefficients of the smoothing filter.

Arguments

Nleft

An integer specifying the number of data points to the left of each point to include in the filter.

Nright

An integer specifying the number of data points to the right of each point to include in the filter.

Order

An integer specifying the order of the derivative desired. For smoothing, use order 0. To find the smoothed first derivative of the signal, use order 1, for the second derivative, use order 2, etc.

Degree

An integer specifying the degree of smoothing polynomial. Typical values are 2 to 4. Lower values for Degree will produce smoother results but may introduce bias, higher values for Degree will reduce the filter bias, but may "over fit" the data and give a noisier result.

Keywords

DOUBLE

Set this keyword to force the computation to be done in double-precision arithmetic.

Examples

The following example creates a noisy 400-point vector with 4 Gaussian peaks of decreasing width. It then plots the original vector, the vector smoothed with a 33-point Boxcar smoother (the SMOOTH function), and the vector smoothed with 33-point wide Savitzky-Golay filter of degree 4. The bottom plot shows the first derivative of the noisy signal and the first derivative using the Savitzky-Golay filter of degree 4:

n = 401 ; number of points 
np = 4  ; number of peaks 
; Form the baseline: 
y = REPLICATE(0.5, n) 
; Index the array: 
x = FINDGEN(n) 
; Add each Gaussian peak: 
FOR i=0, np-1 DO BEGIN 
   c = (i + 0.5) * FLOAT(n)/np ; Center of peak 
   peak = -(3 * (x-c) / (75. / 1.5 ^ i))^2 
   ; Add Gaussian. Cutoff of -50 avoids underflow errors for 
   ; tiny exponentials: 
   y = y + EXP(peak>(-50)) 
ENDFOR 
; Add noise: 
y1 = y + 0.10 * RANDOMN(-121147, n) 
 
 
!P.MULTI=[0,1,3] 
 
; Boxcar smoothing width 33: 
PLOT, x, y1, TITLE='Signal+Noise; Smooth (width33)' 
OPLOT, SMOOTH(y1, 33, /EDGE_TRUNCATE), THICK=3 
 
; Savitzky-Golay with 33, 4th degree polynomial: 
savgolFilter = SAVGOL(16, 16, 0, 4) 
PLOT, x, y1, TITLE='Savitzky-Golay (width 33, 4th degree)' 
OPLOT, x, CONVOL(y1, savgolFilter, /EDGE_TRUNCATE), THICK=3 
 
; Savitzky-Golay width 33, 4th degree, 1st derivative: 
savgolFilter = SAVGOL(16, 16, 1, 4) 
PLOT, x, DERIV(y1), YRANGE=[-0.2, 0.2], TITLE=$ 
 'First Derivative: Savitzky-Golay(width 33, 4th degree, order 1)' 
OPLOT, x, CONVOL(y1, savgolFilter, /EDGE_TRUNCATE), THICK=3

The following is the resulting plot. Notice how the Savitzky-Golay filter preserves the high peaks but does not do as much smoothing on the flatter regions. Note also that the Savitzky-Golay filter is able to construct a good approximation of the first derivative.