BeginPackage["Waves`"]
(* Author: Brani Vidakovic, ISDS, Duke University ;
Functions Dec and Comp are based on M. V. Wickerhauser's
mathematica program; December 1994*)

Mirror::usage = "Mirror[_filter_] gives the mirror  \
filter for the input _filter_. This is an adjoint  \
operator H* of the operator H corresponding to  _filter_." 

WT::usage = "WT[_vector_,  _filter_]  performs the direct \
 wavelet transformation of the data vector  _vector_. \
The  wavelet  base is chosen by _filter_. The length \
of the vector  _vector_  has to be a degree of 2."

WR::usage = "WR[_vector_,  _filter_] gives the wavelet \
reconstruction algorithm. From the set of  wavelet \
coefficients _vector_ the data set is reconstructed. \
The wavelet base is chosen by _filter_."

Dec::usage = "An auxiliary function needed for the \
direct wavelet transformation. See WT."

Comp::usage = "An auxiliary function needed for the \
inverse wavelet transformation (wavelet reconstruction  \
algorithm). See WR."

 Begin["`Private`"]


Mirror[ filter_List]:= Module[{fl=Length[filter]},
Table[ -(-1)^i filter[[fl+1-i]], {i, 1, fl}]];


Dec[ vector_List, filter_List]:=  Module[
 {vl= Length[vector], fl=Length[filter]},
 Table[ 
      Sum[ filter[[m]] vector[[Mod[2 k+m - 3, vl]+1 ]],
         {m,1,fl}],
      {k,1,vl/2}] 
 ];


Comp[ vector_List, filter_List]:= Module[
 { temp=Table[0,{i,1,2 Length[vector]}],
 vl=Length[vector], fl=Length[filter]},
 Do[ temp[[ Mod[2 j + i -3, 2 vl]+1]] +=
        vector[[j]] filter[[i]],
      {j, 1, vl}, {i, 1, fl}];
  temp];


WT[ vector_List, filter_List]:=
Module[ { wav={}, c,d, ve=vector, H=filter, 
G=Mirror[filter]},
While[ Length[ve] > 1,
     lev=Log[2,Length[ve]]-1;
     c = Dec[ve, H];
     d = Dec[ve, G];
     wav= Join[ wav, d ];
     ve = c]; Join[wav, c] ];

WR[ vector_List, filter_List]:=
Module[ {i=1, vl=Length[vector], c=Take[vector,-1],
        d=Take[RotateRight[vector,1],-1],
        mirrorf=Mirror[filter], cn, dn, k=1},
While[ i <= vl/2 ,  
       k += i; 
       i= 2 i;
       cn=Comp[c, filter]+Comp[d, mirrorf];
       dn=Take[RotateRight[vector, k], -i ];
       c=cn;
       d=dn;
     ];
c  ];

End[ ]

EndPackage[ ]