function linerem,lamda,spec,sig=sig,cell=cell,nsigma=nsigma,$ bkgval=bkgval,bkgerr=bkgerr,quiet=quiet,posve=posve,negve=negve,$ _extra=e ;+ ;function linerem ; remove lines from input spectrum and return the "cleaned" spectrum ; ;syntax ; cleanspec=linerem(lamda,spec,sig=sig,cell=cell,nsigma=nsigma,$ ; bkgval=bkgval,bkgerr=bkgerr,/quiet,/posve,/negve) ; ;parameters ; lamda [INPUT; required] wavelengths at which spectrum ; is defined. ; spec [INPUT; optional] the spectrum. ; NOTE: if not given, LAMDA is taken to be SPEC and ; the array indices are taken to be LAMDA ; ;keywords ; sig [INPUT] error at each point; if not given, the ; errors are taken to be 1+sqrt(abs(SPEC)+0.75). ; nsigma [INPUT; default: 4] multiple of SIG to consider ; as a threshold for detection of features ; cell [INPUT] 1D filter to use in computing the background ; * default is [1,1,0,0,0,1,1] ; * if scalar, then [IC+1,ICC,IC+1], where ; IC=intarr(2*CELL), ICC=intarr(2*CELL+1) ; bkgval [OUTPUT] final background values at each bin ; NOTE: This is essentially a smoothed version of ; the cleaned spectrum! ; bkgerr [OUTPUT] error estimates on BKGVAL ; quiet [INPUT] if set, doesn't show, doesn't tell ; posve [INPUT] if set, removes only +ve deviations ; negve [INPUT] if set, removes only -ve deviations ; _extra [JUNK] here only to prevent crashes ; ;description ; 1. make sure that the spectrum is defined on a regular grid. ; (if not, rebin [NOT IMPLEMENTED!]) ; 2. convolve the spectrum with background cell to determine ; local background. ; 3. also propagate errors (assume gaussian; if poisson, use ; gaussian approximation) ; 4. compare local value with local background ; 5. flag those bins which are significantly different from local ; background (use +-NSIGMA*SIG as a threshold value) ; 6. reset the flagged bin values to local background values ; 7. repeat 2-6 until no new bins are flagged ; ;these are all the smoothing tools in PINTofALE ; ALALOESS() makes a loess curve ; CLTSMOOTH() removes -ves from background subtracted spectra or light curves ; CONV_RMF convolves with specified RMF ; HAARTRAN() smoothes by filtering on Haar wavelet coefficients ; LINEREM() removes lines from spectrum by iterative S/N filtering ; NOISMOOTH() does boxcar accumulation a la Ebeling's asmooth ; REGROUP() accumulates from one end ; SCRMF does fast convolution using optimized RMF ; SMOOTHIE does peak-descent accumulation to build up S/N ; SPLAC computes a segmented piecewise linear approximations ; UNKINK() removes sharp kinks from a curve ; VARSMOOTH() does local smoothing at specified scales ; VOORSMOOTH() does fast local smoothing at specified scales ; ;history ; vinay kashyap (Oct96) ; changed keyword SIGMA to SIG (VK; Jun98) ; CELL was being altered if input as scalar -- corrected; altered ; behavior of output when no deviations are found (VK; Oct98) ; now if QUIET>1, won't ask to quit if many bins are found (VK; 99Aug) ; changed POS and NEG to POSVE and NEGVE (VK; MMJul) ;- np=n_params(0) if np eq 0 then begin print, 'Usage: clean_spectrum=linerem(lamda,spectrum,sig=sig,nsigma=nsigma,$' print, ' cell=cell,bkgval=bkgval,bkgerr=bkgerr,/quiet,/posve,/negve)' print, ' removes lines from spectrum and returns continuum' return,-1L endif ; relabel x=lamda & nx=n_elements(x) & ny=n_elements(spec) if ny ne nx then begin y=x & x=lindgen(nx) endif else y=spec ; check keywords if not keyword_set(sig) then sig=sqrt(abs(y)+0.75)+1. if not keyword_set(nsigma) then nsigma=4. if keyword_set(cell) then begin csize=cell nc=n_elements(cell) if nc eq 1 then begin cc=abs(cell(0)) ic=intarr((2*cc)>1) & icc=intarr(2*cc+1) & csize=[ic+1,icc,ic+1] nc=n_elements(csize) endif if nc gt nx then begin message,'Background cell too large',/info & return,y endif endif else csize=[1,1,0,0,0,1,1] if not keyword_set(quiet) then quiet=0 if keyword_set(posve) and keyword_set(negve) then begin posve=0 & negve=0 endif ; initialize norm=total(csize) if norm le 0. then begin c1='background cell has zero area! ignoring normalization' message,c1,/info endif ; check to see that the spectrum is defined on a uniform grid dx=x(1:*)-x & ddx=dx(uniq(dx,sort(dx))) & nddx=n_elements(ddx) odx=where(abs(dx-median(dx))/median(dx) gt 1e-3,modx) if modx gt 1 then begin c1='spectrum not on uniform grid... convolutions may result in nonsense' message,c1,/info if not quiet then print, 'there are '+strtrim(nddx,2)+' bin sizes:',ddx if quiet lt 2 then begin c1='hit any key to continue, q to return, x to stop' & print,c1 c1=get_kbrd(1) endif else c1='' if strlowcase(c1) eq 'q' then return,y if strlowcase(c1) eq 'x' then begin print,'there are '+strtrim(nx,2)+' bins and '+strtrim(nddx,2)+$ ' unique bin widths' help,x,y,dx,ddx stop,'type RETURN,Y to return sans change' endif endif ; convolve spectrum with cell to get local background bkg=convol(y,csize,/edge_truncate) & if norm ne 0 then bkg=bkg/norm ; propagate errors bge=convol(sig^2,abs(csize),/edge_truncate) & if norm ne 0 then bge=bge/norm bge=sqrt(bge) ; find significant deviations dely=(y-bkg) & nok=where(abs(dely) gt nsigma*abs(bge)) if keyword_set(posve) then nok=where(dely gt nsigma*abs(bge)) if keyword_set(negve) then nok=where(dely lt -nsigma*abs(bge)) if nok(0) eq -1 then begin message,'no significant deviations.. spectrum is the continuum?',/info bkgval=bkg & bkgerr=bge return,y ;hey. endif y(nok)=bkg(nok) inew=lonarr(nx) & inew(nok)=1 ; show if not quiet then begin print,strtrim(long(total(inew)),2)+' bins reset' print,'plotting lambda v/s spectrum, line-removed spectrum, bkg w. errors' plot,x,spec,psym=10,/xs & oplot,x,y,col=100 & oplot,x,bkg,col=150 oplot,x,bkg+nsigma*bge,col=150,line=2 oplot,x,bkg-nsigma*bge,col=150,line=2 endif ; iterate y1=y & bkg1=bkg & sigg=sig & sig(nok)=bge(nok) while total(inew) gt 0 do begin ; all as above... bkg1=convol(y1,csize,/edge_truncate) & if norm ne 0 then bkg1=bkg1/norm bge1=convol(sigg^2,abs(csize),/edge_truncate) if norm ne 0 then bge1=bge1/norm bge1=sqrt(bge1) dely=(y1-bkg1) & nok=where(abs(dely) gt nsigma*abs(bge1)) if keyword_set(posve) then nok=where(dely gt nsigma*abs(bge)) if keyword_set(negve) then nok=where(dely lt -nsigma*abs(bge)) inew=lonarr(nx) if nok(0) ne -1 then begin inew(nok)=1 & y1(nok)=bkg1(nok) & sigg(nok)=bge(nok) endif if not quiet then begin print,strtrim(long(total(inew)),2)+' bins reset in this iteration' plot,x,spec,psym=10,/xs & oplot,x,y1,col=100 & oplot,x,bkg1,col=150 oplot,x,bkg1+nsigma*bge1,col=150,line=2 oplot,x,bkg1-nsigma*bge1,col=150,line=2 endif endwhile bkgval=bkg1 & bkgerr=sigg return,y1 end