;
; Time-stamp: <2000-11-01 16:48:54 baron>
;
pro wsmooth,name,thresfac
;This example demonstrates the use of IDL's discrete wavelet transform
;to remove noise from an observed spectrum.
;
;This is likely the best smoothing algorithm available
; 
;It has been improved by the addition of a chi-square test
;to guarantee that no important information is lost
;
; Begin by choosing the number of wavelet coefficients to use and a
; threshold value:
; Open the data file, read an spectrum
; and close the file:
;
filename1=name
filelength,filename1,npts1
;
readit,filename1,npts1,wl1,flux1
;
;
; Expand the image to the nearest power of two using cubic
; convolution, and transform the image into its wavelet
; representation using the WTN function:
;

base=floor(alog(npts1)/alog(2))+1
pwr = 2^base
print,'pwr = ',pwr
flux1_e = CONGRID(flux1, pwr, /CUBIC)
wl1_e= congrid(wl1,pwr,/cubic)
coeffs1 = 12 & thres1 = thresfac*max(flux1_e)
wtn_flux1 = WTN(flux1_e, coeffs1)
;
; now we do the compression
;
wtn_flux_1 = fltarr(N_ELEMENTS(flux1_e))
wtn_noise=wtn_flux_1
index_good = where(abs(wtn_flux1) ge thres1) 
index_bad = where(abs(wtn_flux1) lt thres1) 
wtn_flux_1[index_good] = wtn_flux1[index_good]
wtn_flux_1[index_bad] = 0.
wtn_noise[index_good] = 0.
wtn_noise[index_bad] = wtn_flux1[index_bad]
print,n_elements(flux1_e),n_elements(wtn_flux_1)
;
PRINT, 'Percentage of elements under threshold: ',$
   100.*N_ELEMENTS(WHERE(ABS(wtn_flux1) LT thres1, $
   count1)) / N_ELEMENTS(flux1_e)
openw,1,'original.dat'
for i=0,pwr-1 do printf,1,wl1_e(i),flux1_e(i)
close,1
openw,1,'wavelets.dat'
for i=0,pwr-1 do printf,1,wl1_e(i),wtn_flux_1(i)
close,1
;
;
;
; Apply the inverse wavelet transform to the image:
;
flux_1_c = WTN(wtn_flux_1, COEFFS1, /INVERSE)
noise_flux =  WTN(wtn_noise, COEFFS1, /INVERSE)
; Calculate and print the amount of data used in reconstruction of
; the image:
PRINT, 'The spectrum is reconstructed from:', $
   100.0 - (100.* count1/N_ELEMENTS(flux1_e)),$
   '% of original image data.'
;
; check the chi^2 of the original and recontructed spectra
;
rel=0.1
fluxchi2,flux1_e,flux_1_c,wl1_e,rel,chisq,prob
print,"chi_sqr = ",chisq," prob = ",prob
;
;
; Finally, write out the reconstructed spectrum
;
;
outputfile='original.swv.dat'
openw,1,outputfile
for i=0,pwr-1 do printf,1,wl1_e(i),flux_1_c(i)
close,1

outputfile='noise.dat'
openw,1,outputfile
for i=0,pwr-1 do printf,1,wl1_e(i),noise_flux(i)
close,1
;
end