BeginPackage["Wavelet1Dv1`"]
(* This package performs both a forward and inverse 1D
wavelet transform with perfect reconstruction and the
ability to decompose a signal of length 2^M into the
maximum of M levels of decomposition. *)
(* This program was written in May 1993 and revised
April 1994 by:
James F. Scholl
Rockwell Science Center
1049 Camino Dos Rios
Thousand Oaks, CA 91360
email: jfs@risc.rockwell.com *)
(* This is the Mathematica v.2.x implementation of a
version of the FWT developed by N. Getz in Dec. 1992
(Memo No. UCB/ERL M92/138, Electronics Research Lab. U. of
California, Berkeley CA 94720)... This is an independent
version of his code, not a Mathematica translation of
his MATLAB code in the appendix of his report! Keep
in mind this may not be the best looking or most
efficient Mathematica code in the world so helpful
comments from you users is always welcome. *)
(* If interested further in Wavelets, consult the two
books both published by the Society of Industrial and
Applied Mathematics (SIAM), 3600 University City Science
Center, Philadelphia, PA 19104-2688, USA:
Wavelets: Algorithms & Applications by Yves Meyer
Ten Lectures on Wavelets by Ingrid Daubechies
and references therein. *)
(* This version of the FWT (called the Fast Discrete
Periodic Wavelet Transform) implicitly pads the finite
signal periodically, but due to Getz's ingenious
implementation, there is no data padding at all.
Furthermore, the wavelet filters corresponding to any
particular orthonormal basis systems adaptively change
their length and the coefficients when the length of the
signal gets smaller than that of the signal. For further
details see the above reference. *)
(* Some Caveats: Only even length filters can be used
here if one wishes to do a full-level
analysis.
Also: block coding with this algorithm
may introduce spurious high frequency
coefficients - there is still perfect
reconstruction. One must be a little
careful when splitting up a large chunk
of data!. This is not a major problem,
just a minor inconvenience. *)
(* I appreciate any comments, suggestions on making this
code run better and faster (future versions of Mathematica
will enable this code to be compiled and thus run much
faster). Please contact me at the above mail or electronic
address with these suggestions or comments. *)
DPWT1D::usage = "DPWT1D[dat,lfil,lev] performs a 1D discrete
wavelet transform of the input signal dat w.r.t the wavelet
basis system represented by the input low pass filter lfil.
The high pass filter hfil is automatically calculated from
lfil. The transform is performed by doing lev levels of
decomposition using a recursive pyramid filtering /
downsampling scheme. It is preferred that the
input data length be a power of 2, but even-length data can
be used, but the number of levels of decomposition may be
limited. For J levels of decomposition and signal length 2^L,
the output is in the form {AJ,{DJ,D(J-1),...,D2,D1}},
where AJ is the approximation subsignal of length 2^(L-J)
and Di is the ith detail subsignal of length 2^(L-i). "
IDPWT1D::usage = "IDPWT1D[wdat,lfil] performs the inverse 1D
discrete wavelet transform on the input coefficients wdat.
The input wdat must in the form {AJ,{DJ,D(J-1),...D2,D1}},
where AJ and Di are the approximation subsignal and the ith
detail subsignal for a J-level forward decomposition. The
number of levels in the recursive pyramid upsampling / filter
procedure is automatically determined from the length of wdat.
The wavelet basis system is represented by the low pass filter
coefficients lfil; the high pass filter is calculated within
the routine."
Begin["`Private`"]
(* Here are some examples of filters representing the scaling
functions of some systems: Haar, Daubechies' D4, and
the 2nd and 4th order Coifman systems C2, and C4.
These are low-pass filters. The high-pass filters
corresponding to the primary wavelets are automatically
calculated from the low-pass filter coeffcients in the
routines:
Haar = {0.5, 0.5},
D4 =0.5 * {N[(1 + Sqrt[3])/4], N[(3 + Sqrt[3])/4],
N[(3 - Sqrt[3])/4], N[(1 - Sqrt[3])/4]},
Coif2 = {-0.051429728471,0.238929728471,
0.602859456942, 0.272140543058, -0.051429972847,
-0.011070271529},
Coif4 = {0.011587596739, -0.02932013798, -0.04763959031,
0.273021046535, 0.574682393857, 0.294867193696,
-0.054085607092, -0.042026480461, 0.016744410163,
0.003967883613, -0.001289203356, -0.000509505399}. *)
LHtilde[low_List,high_List,ls_Integer]:=
Module[{ltil,htil,Nit2,kk,ii},
(* This routine adjusts the filter coefficients for
data sets that are smaller than the filter length --
for data of length two the filter / basis system
becomes the Haar. *)
Nit2 = Length[low];
If[Nit2 <= ls, Return[{low, high}],
ltil = Table[Sum[low[[1 + Mod[ii,ls] + kk*ls]],
{kk,0,Floor[(Nit2 - 1 - Mod[ii,ls])/ls]}],
{ii,0,Min[Nit2,ls]-1}];
htil = Table[Sum[high[[1 + Mod[ii,ls] + kk*ls]],
{kk,0,Floor[(Nit2 - 1 - Mod[ii,ls])/ls]}],
{ii,0,Min[Nit2,ls]-1}]]; {ltil,htil}]
indexer[arr_List, ind_List]:= Module[{la=Length[ind],
arr2,iii},
(* This routine does filter / signal bookkeeping for
both the forward and inverse filter convolutions. *)
arr2 = arr[[ Table[ind[[iii]]+1 , {iii, 1, la}] ]];
arr2]
ForwardConv[low_List, high_List, dat_List] :=
Module[{ltl,htl,lh=Length[low],lnf=Length[dat],
lnf2,ii,jj,mung,arfd,atmp,dtmp},
(* This routine performs one level of analysis filtering-or
decomposition-or forward filtering-whatever your
terminology) and downsampling with a factor of 2.
The covolutions of the data with the low and high pass
filters are done using dot products and
some nifty bookkeeping with the "indexer" routine. *)
lnf2 = lnf/2;mung=Min[lnf,lh]-1;
{ltl,htl} = LHtilde[low,high,lnf];
arfd = Table[ Table[(2*ii)+jj, {jj, 0, mung}],
{ii, 0, lnf2-1}]; arfd = Mod[arfd,lnf];
atmp = Table[ltl . indexer[dat, arfd[[ii]] ],
{ii, 1, lnf2}];
dtmp = Table[htl . indexer[dat, arfd[[ii]] ],
{ii, 1, lnf2}]; {atmp, dtmp} /.(0->(N[10^-20])) ]
InverseConv[low_List, high_List, atmp_List, dtmp_List] :=
Module[{ltl,htl,lh=Length[low],lh2,lh22,lnf=Length[atmp],
lnf2,ii,jj,dat,even,odd,arfd,arff,mung},
(* This routine performs one level of synthesis filtering-or
reconstruction-or inverse filtering-whatever your
terminology) and upsampling by a factor of 2.
The convolutions of the data with the low and high pass
filters are done using dot products and
some nifty bookkeeping with the "indexer" routine. *)
lnf2 = 2*lnf; lh2=lh/2; mung=Min[lh2,lnf];
{ltl,htl} = LHtilde[low,high,lnf2];lh22=Length[ltl]/2;
even = Table[2*ii, {ii, 0, lh22-1}];
odd = Table[(2*ii)+1, {ii, 0, lh22-1}];
arff = Flatten[Table[{even,odd},{lnf}],1];
arfd = Table[ Table[Ceiling[(ii+1)/2]-jj, {jj,1,mung}],
{ii,0,lnf2-1}]; arfd = Mod[arfd,lnf];
dat = Table[
((indexer[ltl,arff[[ii]]] . indexer[atmp,arfd[[ii]]]) +
(indexer[htl,arff[[ii]]] . indexer[dtmp,arfd[[ii]]]) ),
{ii, 1, lnf2}]; 2.0*dat]
DPWT1D[dat_List,lfil_List,lev_Integer] :=
Module[{lh=Length[lfil],ls=Length[dat],hfil,
app={},detl={},kk,ii,tad,tadt,tdd},
(* This routine transform codes a signal with a wavelet
basis system represented by the input array "lfil"
which is the low pass filter corresponding to the scaling
function...the high pass filter "hfil" corresponding
to the primary wavelet is automatically calculated from
lfil. The number of decomposition steps in the
transform is input as "lev." The maximum number that
lev has is L, where the length of the data is 2^L. *)
(* The output data is given in the form of a list of lists.
For example, if we do a four level decomposition on a
signal of length 2^M, we get the subsignals A4, D4, D3,
D2, D1, listed in increasing resolution. The output is
given in the format {A3, {D4, D3, D2, D1}} so that one
can easily pick out and examine the subsignals of
interest. *)
(* A check for length and levels: *)
If[lev > Log[2,Length[dat]], Return["lev too large"]];
hfil = Table[((-1)^(kk+1))*lfil[[1+(lh-1)-kk]],
{kk,0,lh-1}];
For[ii=1;{tad,tdd}=ForwardConv[lfil,hfil,dat];
tadt=tad;app=AppendTo[app,tad];
detl=AppendTo[detl,tdd], ii < lev, ii++,
{tad,tdd}=ForwardConv[lfil,hfil,tadt];tadt=tad;
app=AppendTo[app,tad];detl=AppendTo[detl,tdd]];
{Last[app],Reverse[detl]}]
IDPWT1D[wdat_List,lfil_List] :=
Module[{lh=Length[lfil],ls=Length[dat],hfil,
temp,lev,kk,ii},
(* This routine inverse transform codes a signal with a
wavelet basis system represented by the input array "lfil"
which is the low pass filter corresponding to the scaling
function...the high pass filter "hfil" corresponding
to the primary wavelet is automatically calculated from
lfil. The number of decomposition steps in the inverse
transform is calculated by the length of the detail
function array. *)
lev = Length[ wdat[[2]] ];
hfil = Table[((-1)^(kk+1))*lfil[[1+(lh-1)-kk]],
{kk,0,lh-1}];
For[ii=2;temp = InverseConv[lfil,hfil,wdat[[1]],wdat[[2,1]]],
ii <= lev, ii++,
temp = InverseConv[lfil,hfil,temp,wdat[[2,ii]] ]];
temp]
End[]
EndPackage[]