\setcounter{equation}{0}\section{Overview of Commands}

We assume that the readers are familiar with MATLAB, both as
an interactive package and as a programming language, and that
the M-files of the commands discussed in this paper have been properly
installed (see Section~5).  The best way to proceed is
to invoke a MATLAB session and try the
following hands-on tutorial.

First let us generate a sample signal to experiment with.

\begin{quote}
\verb+>> t = 0:0.01:1;+ 
\LEFTLINE 
\verb+>> x = t .* sin(20*t);+
\end{quote}

The first line generates a linear vector {\tt t}, useful as
a shorthand in the next line to construct {\tt x} (or additional
sample signals in the future).  Another interesting sample
signal is {\tt  x = t .* cos(10*t) .* sin(20*t)}.

The command to invoke a $W$-transform is {\tt kwt}.  In its simplest form
\begin{quote}
\verb+>> y = kwt(x);+
\end{quote}
transforms {\tt x} into a new vector {\tt y} using the default QS transform.
To see what {\tt kwt} does to {\tt x}, we graph both the input and
output vectors using
\begin{quote}
\verb+>> plot(x), plot(y)+
\end{quote}

\begin{center}
{\mbox{\psfig{file=plotx.eps,height=2.0in,width=5.3in}}}
\end{center}
\begin{center}
Figure 1. Plot of original signal {\tt x = t sin(20t)}
\end{center}
\begin{center}
\hspace*{2mm} {\mbox{\psfig{file=ploty.eps,height=2.0in,width=5.18in}}}
\end{center}
\begin{center}
Figure 2. Plot of the QS transform of {\tt x} 
\end{center}

\vspace*{\baselineskip}

One sees clearly that {\tt y} is made up of two parts.  Its first
half resembles the original signal {\tt x}, while the second half
consists of relatively small components.  Indeed, we can ask
{\tt kwt} to split the output signal into two parts using the following
form of the command:
\begin{quote}
\verb+>> [y1, y2] = kwt(x);+
\end{quote}

The commands
\begin{quote}
\verb+>> length(x), length(y), length(y1), length(y2)+
\end{quote}
give the lengths of  {\tt x}, {\tt y}, {\tt y1}, and {\tt y2;} they are
101, 101, 51, and 50, respectively.  It is a general
rule that the input and output vectors are equal in length.
The split output vectors are each half the length of
the input vector if it is even; otherwise, {\tt y1} has one more component
than {\tt y2.} 

The versatility of the command {\tt kwt} lies in its ability to handle
correctly input signals of different formats: row or column vectors
as well as matrices.  In the last case, a two-dimensional $W$-transform (once in the
horizontal direction and once in the vertical direction) is
performed to produce four submatrices.  If one has a matrix $ L $ representing
a grayscale image, such as the popular Lena (intensity
range: 0--255), one should try
\begin{quote}
\verb+>> image(kwt(L)/16), axis('image'), colormap(gray(256))+
\end{quote}

The inverse $W$-transform is one of the following:
\begin{quote}
\verb+>> ikwt(y)+
\LEFTLINE 
\verb+>> ikwt(y1,y2)+
\end{quote}
One should verify that each of these is practically (within
the accuracy of numerical error) identical to the original 
vector {\tt x}.

Lossy compression of a signal is achieved by discarding small
components of {\tt y2}, or by quantizing.  The command
\begin{quote}
\verb+>> cy2 = largesta(10, y2);+
\end{quote}
produces a new vector that retains the largest (in absolute
value) ten components of {\tt y2} and truncates the
rest to 0.  The actual number of components retained may be different from 10
if more than one component ties for the tenth place.
The vector {\tt cy = [y1, cy2]} is now a (lossily) compressed
representation of the original signal {\tt x}.  It has only 61 nonzero
components.  It takes, however, more than 61 data units to record {\tt cy},
since both the location and magnitude of the nonzero elements
of {\tt y2} need to be stored.  
Furthermore, the magnitude of the components of $ \yo $ is no longer
within the range assigned to the original signal
and thus may require more bits to represent.
With more advanced techniques, the result can usually be stored by using
just a few extra units in addition to the 61 units.
The compression ratio achieved is thus approximately $ 100/65 $.

Alternatively, the command
\begin{quote}
\verb+>> cy2 = largesta(0.2,y2);+
\end{quote}
retains the largest (in absolute value) two-tenths of the components of $ \yt $.
The compression ratio can be raised by lowering the $ 0.2 $ threshold, the
first input argument to {\tt largesta}.

Another tool for compression is the command
\begin{quote}
\verb+>> quant(y2,5)+
\end{quote}
which quantizes the components of {\tt y2} by using a quantization interval of
length~5.

A much higher compression ratio is achievable by repeating the 
compression process
on {\tt y1}, once or several more times, leading naturally to
the method of multiresolution analysis.  In theory, one can continue the process to as many levels
as possible until the most recently obtained {\tt y1} is of length 2.
In practice, however, after a few levels, further compression
produces either intolerable errors or very little 
additional benefit.
To obtain a three-level multiresolution analysis on {\tt x} is simple:
\begin{quote}
\verb+>> z = wma(x,3);+
\end{quote}
When {\tt x} is a row vector, the structure of {\tt z} is
{\tt [y31 y32 y22 y12]},
where {\tt y31} denotes the third-level {\tt y1} and so on.  

Since a multiresolution analysis produces a variable number of  wavelets, it is
not feasible to ask {\tt wma} to output the separate wavelets
automatically.  We provide another function to carve out
specified wavelets from the output of {\tt wma}.  For 1D
signals,
\begin{quote}
\verb+>> w = maw(z,2);+
\end{quote}
gives the second-level multiresolution analysis wavelet from {\tt z},
obtained above by using {\tt wma.} Similarly,
\begin{quote}
\verb+>> v = maw(z,3,1);+
\end{quote}
gives the {\tt y1} vector at the top level of a three-level multiresolution analysis.

For 2D signals, the third input argument specifies which of the
four submatrices is to be carved out.
One can also call the interactive m-file {\tt wmai} 
(which is not a function) to do the multiresolution analysis.
MATLAB then prompts for the signal and the number of 
levels.
\begin{quote}
\verb+Enter signal to be analyzed: +{\sl L}
\LEFTLINE 
\verb+Enter number of levels     : +{\sl 3}
\end{quote}
The results of the three-level analysis are stored in the signals
{\tt WL11}, {\tt WL12}, {\tt WL21}, etc.

\newpage
For visualization, 
\begin{quote}
\verb+>> wplot(z,3)+
\end{quote}
plots the components {\tt y31,} {\tt y32,} {\tt y22,} {\tt y12} of {\tt z} to give the
familiar third-level wavelet coefficient diagram.  Since the wavelet coefficients are
usually too small to show up significantly on the graph, we
use the command
\begin{quote}
\verb+>> wplot(z,3,0.01)+
\end{quote}
to magnify the second-level wavelet coefficients $ 100=1/0.01 $ times,
the third-level coefficients $ 100^2 $ times, and so on.

\begin{center}
\fbox{
\parbox{5in}{
\vspace*{\baselineskip}
\psfig{file=wplot.eps,height=2.5in,width=5.0in}
\vspace*{\baselineskip}
}
}
\end{center}
\begin{center}
Figure 3. Plot of multiresolution analysis wavelet components
\end{center}

Inverse multiresolution analysis is as expected:
\begin{quote}
\verb+>> X = iwma(z,3);+
\end{quote}
The multilevel wavelet representation of the signal must be assembled into
a single vector or matrix, such as {\tt z}, before {\tt iwma} is called.

The default $W$-transform used by {\tt kwt} and {\tt wma} is controlled by a global
variable {\tt kw} set by the startup file {\tt startup.m}.  Typing
\verb+ >> kw +
reveals the default value {\tt [1 3 3 1]/4}, 
comprising the coefficients in the 2-scale dilation equation for the
wavelet associated with the QS transform.  
In \cite{kt}, the QS transform is defined using {\tt [1 3 3 1]},
but we discover that, in practice, it is more convenient 
to use one-fourth of the vector.  This only affects the
transform by a scaling.
By setting {\tt kw} to other values,
such as 
\begin{quote}
\verb+>> kw = [1, 2, 2, 1]+
\end{quote}
one can choose a different
$W$-transform.  The particular choice {\tt kw = kwdau} ({\tt kwdau} is another vector
automatically set by the startup file {\tt startup.m}) gives the Daubechies
transform of order 4.  Our implementation of the
Daubechies transform differs from the conventional one in the
way it treats the endpoint components; see \cite{kt}.  
Unlike {\tt kwqs} and {\tt kw1}, the vector {\tt kwdau} has six components.

\begin{quote}
\verb+kwdau = [0.4830 0.8365 0.2241 -0.1294 0.2679 -3.7321]+
\end{quote}

A second method to call a specific $W$-transform, without altering the default
value of {\tt kw,} is to supply the parameter vector as an additional
argument to the function invoked.  The commands
\begin{quote}
\verb+>> y = kwt(x,kwdu);+
\LEFTLINE 
\verb+>> y = wma(x,kwdu,3);+
\end{quote}
use the Daubechies transform but retain the QS transform as the default for
future calls.

The decomposition effected by
a $W$-transform is in general not orthogonal.
In \cite{kt}, it
is shown that for the QS transform this is not a real drawback
in practice.  For most reasonable signals,
the error incurred by discarding small components of $ \yt $
is on the same order of magnitude as the optimal
compression obtainable by using orthogonal decomposition.
For the perfectionist, we provide
an {\sl orthogonal compensation} procedure
that can further improve the compression performance.  
Let {\tt e} be
the part of {\tt y2} to be discarded, for instance, 
\begin{quote}
\verb+>> e = y2 - cy1;+
\end{quote}
where {\tt cy1 = largesta(0.1, y2)} as defined before.  Then
\begin{quote}
\verb+>> cy1 = oc(y1,e);+
\end{quote}
produces a vector {\tt cy1} to be used in place of {\tt y1,} in such a way
that the pair {\tt (cy1, cy2)} gives a better approximation to 
{\tt x} than the pair {\tt (y1, cy2)}.

The commands mentioned above provide an adequate set of
building blocks for writing
more complicated M-files for experimentation.
Examples given in Section~4 include M-files to compare two
$W$-transforms, to compress 2D images, and to search for an optimal
$W$-transform.  Other commands in the toolbox include {\tt sf},
which constructs
the scaling function; {\tt wf},
which  constructs the wavelet function, and
{\tt sdil}, which solves a 2-scale dilation equation.
