;
;----------------------------------------------------------------
; PV-WAVE: Signal Processing Toolkit                             |
; Copyright (C) 1994, Visual Numerics, Inc.                      |
; All rights reserved.  Unauthorized reproduction prohibited.    |
;----------------------------------------------------------------
;
FUNCTION WAVELET, param1, param2, param3, Backwards=backwards
;
; $Id: wavelet.pro,v 1.6 2005/08/26 23:55:44 allanw Exp $
;
;+
; Name:
;	WAVELET
; Purpose:
;	Performs a wavelet transform of a signal
; Category:
;	Signal Processing
; Calling Sequence:
;	result = WAVELET(h, x, n)
;	result = WAVELET(h, waveletstruct, /Backwards)
; Inputs:
;	Forwards Wavelet Transform (default):
;	h = a digital filter structure designed by QMFDESIGN
;	x = a digital signal to be filtered
;	n = the number of QMF stages to go through
;
;	Backwards Wavelet Transform (with /Backwards):
;	h = same as above
;	waveletstruct = the output of the forward transform
; Keyword parameters:
;	Backward = if this is set, the original signal is reconstructed
;		from the wavelet data structure (see below) and the QMF
; Outputs:
;	In forward mode, the function returns an unnamed structure which 
;	contains data in a special format.
;
;	When /Backwards is used, that structure is passed in as input
;	and the function returns a reconstruction of the original input
;	signal.
; Common blocks:
;	sigpro_common
; Side effects:
;	None
; Restrictions:
; Procedure:
;	For each n, add another QMF structure on the chain.  Pass the
;	signal through those nodes, saving the 'G' output along the way
;	and passing the 'H' output into the next QMF node.  See the
;	reference for the detailed information on the QMF tree that this
;	creates.
;-
;
@sigpro_common
   IF NOT KEYWORD_SET(backwards) THEN BEGIN
      ;
      ; The forward wavelet transform.
      ;
      IF N_PARAMS() LT 3 THEN MESSAGE, 'Forward transform needs 3 arguments.'
      h = param1
      x = param2
      n = param3
      ;
      ; First off, is there even enough data for the value of 'n'?
      ;
      nx = N_ELEMENTS(x)
      np = N_ELEMENTS(h.b)
      ;
      IF (nx + n*np/2)  LT 2l^n THEN MESSAGE, 'Input data is too short.'
      ;
      ; Make a structure to hold the x vectors at each stage.
      ;
      length = NINT((nx + np)/2.0 - 0.1, /Long) 
      tmpstring = '{'
      FOR i = 1, n DO BEGIN
         tmpstring = tmpstring + ',hminus' + STRING(2l^i - 2) + $
           ':DBLARR(' + STRING(length) + ')'
         IF i LT n THEN BEGIN 
            length = NINT((length + np)/2.0 - 0.1, /Long)
         ENDIF
      ENDFOR
      ;
      ; Now, add the H part of the last QMF structure
      ;
      tmpstring = tmpstring + ',hminus' + STRING(2l^n - 1) + ':DBLARR(' + $
                                                                     STRING(length) + ')' 
      ;
      ; And add a vector of LONGs.  This vector says how long the input
      ; data is at each stage.  When reconstructing the input signal, it
      ; is mathematically impossible to know exactly how long the input
      ; to each stage of the QMF was.  So, we need this additional in-
      ; formation to reconstruct it.
      ;
      tmpstring = tmpstring + ', xlength: LONARR(' + STRING(n) + ')}'
      tmpstring = STRCOMPRESS(tmpstring, /Remove_All)
      status = EXECUTE('result = ' + tmpstring)
      ;
      ; Now that the structure is set up, it's time to fill it in.
      ; For each level of QMF, we take the lowpass output of the
      ; previous filter, and pass it back into QMF, storing the term-
      ; inated outputs in the structure
      ;
      nexth = h
      nextx = x
      FOR i = 1, n DO BEGIN
         ;
         ; Put the length of the input into the xinlength vector
         ;
         result.xlength(i-1) = N_ELEMENTS(nextx)
         ;
         ; Do the QMF
         ;
         QMF, nexth, nextx, result.xlength(i-1), xh, xg
         ;
         ; Take the G part and stash it away in the structure.
         ;
         tmpstring = 'result.hminus' + STRING(2l^i - 2) + '= xg'
         status = EXECUTE(STRCOMPRESS(tmpstring, /Remove_All))
         IF i LT n THEN nextx = xh
      ENDFOR
      tmpstring = 'result.hminus' + STRING(2l^n - 1) + '= xh'
      status = EXECUTE(STRCOMPRESS(tmpstring, /Remove_All))
      ;
   ENDIF ELSE BEGIN
      ;
      ; This is the /Backwards case
      ;
      IF N_PARAMS() LT 2 THEN MESSAGE, 'Backward transform needs 2 arguments.'
      h = param1
      wstruct = param2
      ;
      ; First, let's try to figure out what 'n' was.  Take the number
      ; of tags of the structure and subtract one for the array of 'nx'
      ; and one for the final lowpass output.
      ;
      n = N_TAGS(wstruct) - 2
      status = EXECUTE(STRCOMPRESS('xh = wstruct.hminus' + $
                                   STRING(2l^n - 1), /Remove_All))
      FOR i = n, 1, -1 DO BEGIN
         status = EXECUTE(STRCOMPRESS('xg = wstruct.hminus' + $
                                      STRING(2l^i - 2), /Remove_All))
         ;
         ; Do the backwards QMF
         ;
         QMF, h, newx, wstruct.xlength(i-1), xh, xg, /Backwards
         xh = newx
      ENDFOR
      ;
      ; Now, all iterations are over.  The 'lastxh' is actually the
      ; original input signal (to within some machine noise).
      ;
      result = xh
   ENDELSE
;
   RETURN, result
END