function poly_smooth, data, width, DEGREE=degree, NLEFT=nl, NRIGHT=nr, $
DERIV_ORDER=order, COEFFICIENTS=filter_coef
;+
; NAME:
; POLY_SMOOTH
;
; PURPOSE:
; Apply a least-squares (Savitzky-Golay) polynomial smoothing filter
; EXPLANATION:
; 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
; coefficients.
; CALLING SEQUENCE:
;
; spectrum = poly_smooth( data, [ width, DEGREE = , NLEFT = , NRIGHT =
; DERIV_ORDER = ,COEFF = ]
;
; INPUTS:
; 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.
; RESULTS:
; Function returns the data convolved with polynomial filter coefs.
;
; EXAMPLE:
;
; Given a wavelength - flux spectrum (w,f), apply a 31 point quadratic
; smoothing filter and plot
;
; IDL> cgplot, w, poly_smooth(f,31)
; COMMON BLOCKS:
; common poly_smooth, degc, nlc, nrc, coefs, ordermax
;
; PROCEDURE:
; As described in Numerical Recipes, 2nd edition sec.14.8,
; Savitsky-Golay filter.
; 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.
;
; MODIFICATION HISTORY:
; 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
;-
compile_opt idl2
On_error,2
if N_params() LT 1 then begin
print,'Syntax - smoothdata = ' + $
'poly_smooth( data , width, [ DEGREE = , NLEFT = '
print,f='(35x,A)', 'NRIGHT = , DERIV_ORDER =, COEFFICIENT = ]'
return, -1
endif
common poly_smooth, degc, nlc, nrc, coefs, ordermax
if N_elements( degree ) NE 1 then degree = 2
if N_elements( order ) NE 1 then order = 0
order = ( order < (degree-1) ) > 0
if N_elements( width ) EQ 1 then begin
width = fix( width ) > 3
if (N_elements(nr) NE 1) AND (N_elements(nl) NE 1) then begin
nl = width/2
nr = width - nl -1
endif
endif
if N_elements( nr ) NE 1 then begin
if N_elements( nl ) EQ 1 then nr = nl else nr = 2
endif
if N_elements( nl ) NE 1 then begin
if N_elements( nr ) EQ 1 then nl = nr else nl = 2
endif
if N_elements( coefs ) LE 1 then begin
degc = 0
nlc = 0
nrc = 0
ordermax = 3
endif
if (degree NE degc) OR (nl NE nlc) OR (nr NE nrc) OR $
(order GT ordermax) then begin
degree = degree > 2
ordermax = ( ordermax < 3 ) > order
nj = degree+1
nl = nl > 0
nr = nr > 0
nrl = nr + nl + 1
if (nrl LE degree) then begin
message,"# of points in filter must be > degree",/INFO
return, data
endif
ATA = fltarr( nj, nj )
ATA[0,0] = 1
iaj = indgen( nj ) # replicate( 1, nj )
iaj = iaj + transpose( iaj )
m1_iaj = (-1)^iaj
for k = 1, nr>nl do begin
k_iaj = float( k )^iaj
CASE 1 OF
( k LE nr