function atrous, spec, filter = filter, n_scales = n_scales, check = check
;+
; NAME:
; ATROUS
; PURPOSE:
; Perform a 1-D "a trous" wavelet decomposition on an image.
;
; CALLING SEQUENCE:
; decomposition = ATROUS( spectrum [, filter = filter, $
; n_scales = n_scales, /check]
;
; INPUTS:
; SPECTRUM -- A 1-D Image to Filter
;
; KEYWORD PARAMETERS:
; FILTER -- A 1D (!!) array representing the filter to be used.
; Defaults to a B_3 spline filter (see Stark & Murtaugh
; "Astronomical Image and Data Analysis", Spinger-Verlag,
; 2002, Appendix A)
; N_SCALES -- Set to name of variable to receive the number of
; scales performed by the decomposition.
; CHECK -- Check number of scales to be performed and return.
; OUTPUTS:
; DECOMPOSITION -- A 2-D array with scale running along the 2nd axis
; (large scales -> small scales). The first plane
; of the array is the smoothed image. To recover
; the image at any scale, just total the array
; along the 2nd dimension up to the scale desired.
;
;
; RESTRICTIONS:
; Uses FFT convolutions which edge-wrap instead of mirroring the
; edges as suggested by Stark & Mutaugh. Wait for it.
;
; MODIFICATION HISTORY:
;
; Mon Oct 6 11:49:50 2003, Erik Rosolowsky
; Written.
; 24 April 2005, JohnJohn: Rewrote as a function, converted to work
; with 1D spectra rather than 2D images.
;
;-
on_error, 2
if n_elements(filter) eq 0 then filter = 1./[16, 4, 8/3., 4, 16]
nel = n_elements(spec)
if 1-keyword_set(n_scales) then $
n_scales = floor((alog(nel)/n_elements(filter))/alog(2))
decomp = fltarr(nel, n_scales+1)
sp = spec
for k = 0, n_scales-1 do begin
; Smooth image with a convolution by a filter
smooth = convol(sp, filter, /edge_wrap)
; num_conv, sp, filter, smooth
decomp[0, n_scales-k] = sp-smooth
sp = smooth
; Generate new filter
newfilter = fltarr(2*n_elements(filter)-1)
newfilter[2*findgen(n_elements(filter))] = filter
; note filter is padded with zeros between the images
filter = newfilter
endfor
; Stick last smoothed image into end of array
decomp[*, 0] = smooth
return, decomp
end