;IDL Program: wtestwwb.pro
;PURPOSE: This IDL routine reads in Sky's data, performs WWB functions 
;(DWT etc), and writes the results out into a file. After running this on 
;several platforms, the output can then be compared to check that the WWB 
;functions are working the same for each platform. Timing is performed on
;each of the functions as well to test speed on the different platforms. 
;The program: waveset.pro must be run first to set the COMMON BLOCK variables for WWB.
;INPUT: None (but run waveset.pro first)
;OUTPUT: A Text file with the filename of your choice containing
;results that look like this:
;-----------------------------------------------------------------------
;Amara Graps		2-4-96		Updated: 
;-----------------------------------------------------------------------

;Set the Wavelet to Haar1
WaveType = 'Haar'
ParVal = 1

;Set the file up for output
GET_LUN, wunit
outfile = '                                                     '
READ, 'Enter name with path for output file: ', outfile
outfile = STRTRIM(outfile, 2)
OPENW, wunit, outfile

;Print header information
PRINTF, wunit, ' '
PRINTF, wunit, 'Results of IDL program  wtestwwb.pro'
PRINTF, wunit, 'Written by Amara Graps  February 1996'
PRINTF, wunit, systime()		;current system date and time
format = '$("Platform = ", A5, 2X, "IDL Version = ", A5)'
PRINTF, wunit, format, !version.os, !version.release ;Platform and Version
PRINTF, wunit, ' '


;Read in the dataset
st = 'skywave'
WREADDAT, 2, st, len, wd, sig
pt = 'Data Set: ' + st
WPLOTIT, pt, 1, sig

;Set the QMF,
WINTWAVE, sig, WaveType, ParVal, QMF, LD

;-----------------------------------------
;Perform the Discrete Wavelet Transform
;-----------------------------------------
PRINTF, wunit, 'Discrete Wavelet Transform.'
tt = SYSTIME(1)		;Start timing
wc = WFWTPO(sig,LD,QMF)
delta_tt =  SYSTIME(1) - tt	;End timing

;Write out the transformed coefficients
FOR i = 0, len-1 DO PRINTF, wunit, wc(i)
format = '$("Elapsed Time = ", 1X, F9.6, " seconds. ")'
PRINTF, wunit, format, delta_tt
PRINTF, wunit, ' '
PRINTF, wunit, ' '


;-----------------------------------------
;Perform the Scalegram
;-----------------------------------------
PRINTF, wunit, 'Scalegram.'
tt = SYSTIME(1)		;Start timing
sg = WCALCSCG(sig, QMF)
delta_tt =  SYSTIME(1) - tt	;End timing

;Write out the scalegram values
lensg = n_elements(sg)
FOR i = 0, lensg-1 DO PRINTF, wunit, sg(i)
format = '$("Elapsed Time = ", 1X, F9.6, " seconds. ")'
PRINTF, wunit, format, delta_tt
PRINTF, wunit, ' '
PRINTF, wunit, ' '


;---------------------------------------------
;Perform the Scalogram
;---------------------------------------------
;Note- the current implementation is just
;the discrete wavelet transform coeffs plotted
;in a surface plot,so there's nothing new
;here. A future implementation will use the
;continuous wavelet transform instead.


;-----------------------------------------
;Perform the Multi-Resolution Analysis
;-----------------------------------------
;Note: More than half of this analysis is a
;dyad-split PLOT of the forward/inverse signal
;convolved with the original dataset for EACH
;dyad level. The following routine is the plotted 
;sections taken out and for one dyad only.

PRINTF, wunit, 'Multi-Resolution Analysis.'
tt = SYSTIME(1)		;Start timing
mra = WCALCMRA(sig, LD, QMF)
delta_tt =  SYSTIME(1) - tt	;End timing

;Write out the multi-resolution analysis values
lenmra = n_elements(mra)
FOR i = 0, lenmra-1 DO PRINTF, wunit, mra(i)
format = '$("Elapsed Time = ", 1X, F9.6, " seconds. ")'
PRINTF, wunit, format, delta_tt
PRINTF, wunit, ' '
PRINTF, wunit, ' '

;-----------------------------------------
;Perform the Wavelet Packet Analysis
;-----------------------------------------
;Note: More than half of this analysis are
;dyad-split PLOTS and TREES of the forward/inverse 
;signals. The following routine is just the
;analysis part, not the trees or entropy calculation.

PRINTF, wunit, 'Wavelet Packet Analysis.'
Dsplit = 6	;fineness of the frequency splitting
tt = SYSTIME(1)		;Start timing
wpa = WPANALYS(sig,Dsplit,QMF)
delta_tt =  SYSTIME(1) - tt	;End timing

;Write out the multi-resolution analysis values
lenwpa = n_elements(wpa)
FOR i = 0, lenwpa-1 DO PRINTF, wunit, wpa(i)
format = '$("Elapsed Time = ", 1X, F9.6, " seconds. ")'
PRINTF, wunit, format, delta_tt
PRINTF, wunit, ' '
PRINTF, wunit, ' '


;Close up output file
CLOSE, wunit
FREE_LUN, wunit


END
;************************************************************************