pro sigclip, data_in, threshold=threshold, fwindow, order, degree, help=help, maxiter=maxiter, verbose=verbose, stddevin=stddevin, plot=plot ;NAME: ; pro sigclip.pro (version 2.1) ; ; PURPOSE: ; This program fits a smooth and continuous function ; to input data using a Savitsky-Golay filter and sigma ; clipping. Upon exit, the input data are overwritten with ; the smoothed data. ; ; CALLING SEQUENCE: ; sigclip, data_in [, threshold=threshold, /help] ; ; INPUTS: ; data_in - a 1-dimensional fltarr to be smoothed ; threshold - The sigma clipping threshold ; fwindow - Window size for SAVGOL filter ; order - Order for SAVGOL filter ; degree - Degree for SAVGOL filter ; maxiter - THIS DOES NOTHING. IT IS HERE FOR COMPATIBILITY WITH OLDER SOFTWARE ; verbose - THIS DOES NOTHING. IT IS HERE FOR COMPATIBILITY WITH OLDER SOFTWARE ; stddevin - THIS DOES NOTHING. IT IS HERE FOR COMPATIBILITY WITH OLDER SOFTWARE ; ; OUTPUTS ; ; KEYWORD PARAMETERS: ; help - Print a help message ; plot - Make a plot to show the fitting ; ; EXAMPLE ; ; ; ; REFERENCE: ; ; This program has been approved for release by ; B.J. Rauscher. It contains no proprietary information. ; ; ; MODIFICATION HISTORY: ; Written by: ; Modified: B.J. Rauscher, March 18, 2003 ; Initial release ; ; 4/8/03 -- B.J. Rauscher, Use interpolation to fill ; "gaps" of outliers. Use only one iterations ; since this seems to suffice. ; - ;; Handle request for help if (not keyword_set(help)) then begin help = 0 endif else begin stop, 'Usage: sigclip, data_in [, threshold=threshold, fwindow=fwindow, order=order, degree=degree, maxiter=maxiter, /verbose, /plot, /help]' endelse ;; Set defaults if (not(keyword_set(threshold))) then threshold = 5.0; Default sigma clipping threshold =3 sigma if (not(keyword_set(fwindow))) then fwindow = 64 if (not(keyword_set(order))) then order = 0 if (not(keyword_set(degree))) then degree = 3 if (not(keyword_set(plot))) then plot = 0; Default does not make a plot ;; Create the Savitsky-Golay kernel savgol_filter = savgol(fwindow, fwindow, order, degree) ;; Do the sigma clipping. It seems like only 1 iteration is ;; needed to get essentially everything. for niter = 0, 1 do begin ;; For the 0th iteration, use simple sigma clipping to ;; identify outliers. This works better in practice because ;; the first or last row in an SCA can occasionally be ;; anomalous. if (niter eq 0) then begin djs_iterstat, data_in, mean=mean, sigrej=3, mask=mask endif else begin ;; For subsequent iterations, apply a Savitsky-Golay filter and calculate residuals against this. ;; Smooth the data data_out = convol(reform(data_in), savgol_filter, /EDGE_TRUNCATE) ;; Calculate a residual array residual = sqrt((data_in - data_out)^2) ;; Identify statistical outliers. Use 3 sigma clipping djs_iterstat, residual, sigma=rms_error, sigrej=3, mask=mask endelse ;; The first and last pixels must be good. If they ;; are not, replace them with an average of neighboring pixels ;; Handle pixel 0 if (mask[0] eq 0) then begin djs_iterstat, data_in[0:fwindow/2-1], mean = meanpx, sigrej=3 data_in[0] = meanpx mask[0] = 1 endif ;; Handle last pixel if (mask[n_elements(data_in)-1] eq 0) then begin djs_iterstat, data_in[n_elements(data_in)-fwindow/2:n_elements(data_in)-1], mean = meanpx data_in[n_elements(data_in)-1] = meanpx mask[n_elements(data_in)-1] = 1 endif ;; Use linear interpolation to fill in outliers. ;; First look for statistical outliers for j = 0L, n_elements(data_in) - 2 do begin if (mask[j] eq 0) then begin ;; Find the next "good" pixel for k = j + 1, n_elements(data_in) - 1 do begin if (mask[k] eq 1) then begin slope = (data_in[k] - data_in[j-1])/(k - j + 1) for u = j, k - 1 do data_in[u] = data_in[j-1] + slope * (u - j + 1) j = k break endif endfor endif endfor endfor if (plot) then begin plot, data_in, yrange=[-10,10], psym=1 oplot, data_out, thick=3, color=2 endif ;; Overwrite the input data with the result data_in = data_out end