% \def\bx{\begingroup\small}
% \def\ex{\endgroup}

\setcounter{equation}{0}\section{The MATLAB Programs}

In this section, we present the MATLAB codes of the 
commands given in Section~2.  The algorithms are explained in detail for
the benefit of those who need to adapt the codes to other
languages.

Each command is contained in a separate MATLAB M-file.  
We assume that
the readers are familiar with the rudiments of MATLAB programming.
The commands have not been written
in the most robust form; rigorous
checking for correct formats of input arguments is lacking in 
most cases.  Comments and feedback are appreciated.

When listing our M-files, we leave out the 
comments, except
those pertaining to the usage of the command.  
For reference, we have added line numbers.

\subsection{User-Friendliness}

In this subsection, we discuss our
goal and standard for writing user-friendly MATLAB commands.  
Our objective is to
help readers understand the design of the 
programs presented below.  Fortunately, the MATLAB language
is powerful enough to allow us to implement most of what we 
wanted.

The ability to 
allow some input arguments to be omitted and to assume default
values greatly simplifies the syntax of the commands when used
in common situations.  This feature is especially helpful
to beginners.  

Two techniques are employed  to assign default values.
The global variable {\tt kw} is used to designate the 
default $W$-transform used by {\tt kwt}, {\tt ikwt}, {\tt wma}, etc.
A user can interactively change the default by assigning
new values to {\tt kw}.  
Since many global variables can be dangerous, however,
we refrain from using more.  Instead, for other default
arguments, such as the number of levels used by {\tt wma} 
and the quantization level used by {\tt quant}, the default
values are built into the programs themselves.  A user cannot
conveniently change a default value without modifying the
M-file.  We believe that such
modification is not difficult to do, and we
have supplied easy-to-follow instructions in the M-files
using comments.

If possible, we make a command accept an argument in
any reasonable format.  For instance, {\tt kwt} does not fuss
if the input signal {\tt x} is a row or a column vector and
will produce a row or column output vector accordingly.

In a similar vein, we attempt to make a command 
more flexible in accommodating the order of the input
arguments.  For instance, {\tt wma(4,x)} and {\tt wma(x,4)} 
are both permitted and produce the same answer.
Another example is {\tt largesta(4,y)} and {\tt largesta(y,4)}.

We also try to pack more functionality into a single command.
The actual action and output of a command can be different
depending on the format of the
input or output arguments.  Some examples follow.

\begin{itemize}
\item Most toolboxes have different commands for 1D and 2D
transforms.  Our command {\tt kwt} will perform a 1D or a 2D transform 
according to whether the input signal is 1D or 2D.
Also {\tt kwt} will separate out the wavelet vectors{/}matrices
if there are two output arguments for 1D signals or
four output arguments for 2D signals.

\item The command {\tt ikwt(y1,y2)} or {\tt ikwt(y1,y2,y3,y4)} 
first assembles the input vector{/}matrices into a single vector{/}matrix
before performing the inverse $W$-transform.

\item The command {\tt largesta(r,y)} will extract the largest {\tt r} 
components of {\tt y} if {\tt r} is greater than 1 but the largest
fraction {\tt r} of the total number of components if {\tt r} is less than 1.

\item The command {\tt quant(y,10)} quantizes the components of {\tt y} to the nearest
10, 20, etc., while {\tt quant(y,'10')} quantizes the components of {\tt y} 
into approximately 10 levels.
\end{itemize}

Error-checking is 
part of being user-friendly.  
We repeat our earlier remark that we have
not spent sufficient effort on perfecting this aspect
of our commands (we will do our best when time permits).
We believe that the commands are sufficiently easy to use that
few users will really need extensive error-checking.

\newpage
\subsection{The Startup File}

\bx
\vspace*{1.2\baselineskip}
\begin{verbatim}
      %%% startup.m

   1  global kw
   2  kwqs = [ 1 3 3 1 ]/4;
   3  kw1 = [ 1 2 2 1 ];
   4  kw =kwqs;
   5  dr = sqrt(3);
   6  dau = [1 dr 2*dr-3 dr-2];
   7  ndau = dau/norm(dau);
   8  wdau = ndau(4:-1:1).*[1 -1 1 -1];
   9  kwdau = [ ndau ndau(3)/wdau(3) ndau(1)/wdau(1)];
  10  clear dr dau ndau wdau
\end{verbatim}
\ex

The commands in this file are to be appended to those in one's
startup file, which usually contains path setting and other
routine startup commands.

These commands declare a global variable {\tt kw}; define three
variables {\tt kwqs}, {\tt kwdau}, and {\tt kw1}; 
and set {\tt kw} initially to {\tt kwqs}.
The global variable {\tt kw} contains the parameter vector that determines
the default $W$-transform used by the commands {\tt kwt}, {\tt wma}, and their
families.  The variables {\tt kwqs}, {\tt kwdau}, and 
{\tt kw1} are parameter 
vectors for the QS, Daubechies, and a third sample
transform, respectively.   A general parameter vector has six components, as explained
in Section 3.3.  The fourth and fifth components, when omitted,
take the default values 1 and -1, respectively.
When one has developed her own favorite 
$W$-transforms, their corresponding parameter vectors can be added to the list,
simply following the format of line 2 or 3 to define more
variables.
Line 10 clears the temporary variables used only to construct
{\tt kwdau}.

If these variable names conflict with ones normally used in
one's
own workspace, the reader can change his own notations or modify
all the M-files in this paper accordingly.

\newpage
\subsection{The $W$-transform \, {\tt kwt}} 

According to the general theory outlined in \cite{kt}, 
an arbitrary vector $ h=[h_1,h_2,h_3,h_4] $
and two constants $ c $ and $ d $ are first chosen.  The vector 
\begin{equation}    g=[g_1,g_2,g_3,g_4] = [h_1/c,h_2/c,h_3/d,h_4/d] \end{equation}
together with $ h $ then form a basic pair of vectors that are used
to construct the $W$-matrices
(of even and odd sizes, respectively):
\small 
\begin{equation}    \mbox{\normalsize  \bf W} =  \left( \begin{array}{cccccccccc}    g_1+g_2  & g_3  & g_4 \\   h_1+h_2  & h_3  & h_4 \\   &  g_1  & g_2  & g_3  & g_4 \\   &  h_1  & h_2  & h_3  & h_4 \\   &   &   &  g_1  & g_2  & g_3  & g_4 \\   &   &   &  h_1  & h_2  & h_3  & h_4  \\   &   &   &   &   &   &  \ddots  \\   &   &   &   &   & g_1  & g_2  & g_3+g_4 \\   &   &   &   &   & h_1  & h_2  & h_3+h_4 \end{array} \right )  \Label{w1}  \end{equation}
\normalsize 
and
\small 
\begin{equation}    \mbox{\normalsize  \bf W} =  \left( \begin{array}{cccccccccc}    g_1+g_2  & g_3  & g_4 \\   h_1+h_2  & h_3  & h_4 \\   &  g_1  & g_2  & g_3  & g_4 \\   &  h_1  & h_2  & h_3  & h_4 \\   &   &   &  g_1  & g_2  & g_3  & g_4 \\   &   &   &  h_1  & h_2  & h_3  & h_4  \\   &   &   &   &   &   &  \ddots  \\   &   &   &   &   & g_1  & g_2  & g_3 & g_4 \\   &   &   &   &   & h_1  & h_2  & h_3  & h_4  \\   &   &   &   &   &   &   & h_1  & h_2 +h_3 +h_4  \end{array} \right ).  \Label{w2}  \end{equation}
\normalsize 

The corresponding $W$-transform {\tt kwt(x)} 
of a column vector {\tt x} is computed by
choosing a {\bf W} of the appropriate size, forming
the matrix product {\bf W\tt x}, and then 
separating out the odd and even components of the resulting vector.
\begin{equation}    \mbox{\tt y1} = \mbox{odd \ components of } \mbox{\bf W\tt x} \end{equation}
\vspace*{-8mm}
\begin{equation}    \mbox{\tt y2} = \mbox{even components of } \mbox{\bf W\tt x} \end{equation}
We define {\tt kwt(x)} to be either the pair
{\tt y1} and {y2} or the vector obtained by appending
{\tt y2} to {\tt y1.} 

For a row vector {\tt x}, {\tt kwt(x)} is defined as
the transpose of the
transform of the transpose of {\tt x}.  For a matrix {\tt x},
its 2D $W$-transform is obtained by first applying the 1D $W$-transform to 
each column of {\tt x} and then applying the same $W$-transform to each
row of the resulting matrix.

Most researchers in the wavelet community are, however, more familiar
with another vector, which comprises the coefficients of the dilation
equation and is related to $ h $ by
\begin{equation}    \overline g = [ h_4, -h_3, h_2, -h_1 ] . \end{equation}
We program {\tt kwt} to use an input argument
{\tt k}, obtained by appending the constants $ c $ and $ d $ 
(if they are different from the default values of 1 and -1,
respectively) to $ \overline g $,
to designate the choice of the $W$-transform.
When {\tt k} is omitted, its value is
taken from the global variable {\tt kw}.
Initially, the default parameter vector is 
{\tt kw = [ 1 3 3 1 ]/4}, set by the {\tt startup.m} file.
Although any choice of {\tt k} will lead to a workable transform, 
only a carefully chosen one will result in a transform appropriate for 
compression.  

\bx
\begin{verbatim}
      %       y   = kwt(x)           default W-transform 
      %  [y1, y2] = kwt(x)           1D row/column x
      %  [y1, y2, y3, y4] = kwt(x)   2D matrix x
      %       y   = kwt(x,k)         general W-transform
      %                                determined by k
      
   1  function [y, yw, y3, y4] = kwt(x,k)

   2  if nargin == 1
   3    global kw
   4    if exist('kw') ~= 1, k = [1 3 3 1]/4; else k = kw; end
   5  end

   6  ss=size(x); 
   7  if ss(1) == 1, x = x(:); 
   8  elseif ss(2) > 1 
   9    if nargout == 4, [y, yw, y3, y4] = kwt2(x,k);
  10    else y = kwt2(x,k); end
  11    return
  12  end
\end{verbatim}
\ex
\bx
\begin{verbatim}
  13  if length(k) == 4, c=1; d=-1;
  14  elseif length(k) == 5, c=k(5); d=-1;
  15  elseif length(k) == 6, c=k(5); d = k(6); 
  16  else   k = [1 3 3 1]/4; c=1; d=-1;
  17  end

  18  C=-k(4)/c; D=k(3)/c; A=-k(2)/d; B=k(1)/d; 

  19  s=size(x);
  20  n=floor(s(1)/2);
  21  x1=x(1:2:s(1),:);
  22  x2=x(2:2:s(1),:);

  23  if s(1) == 2*n
  24  y=[ (C+D)*x1(1,:)+A*x2(1,:)+B*x1(2,:)
  25       C*x2(1:n-2,:)+D*x1(2:n-1,:)+A*x2(2:n-1,:)+B*x1(3:n,:)
  26       C*x2(n-1,:)+D*x1(n,:)+(A+B)*x2(n,:) ];
  27  yw=[ c*(C+D)*x1(1,:)+d*(A*x2(1,:)+B*x1(2,:))
  28       c*(C*x2(1:n-2,:)+D*x1(2:n-1,:))+d*(A*x2(2:n-1,:) ...
  29       +B*x1(3:n,:))
  30       c*(C*x2(n-1,:)+D*x1(n,:))+d*(A+B)*x2(n,:) ];

  31  else
  32  y=[ (C+D)*x1(1,:)+A*x2(1,:)+B*x1(2,:)
  33       C*x2(1:n-1,:)+D*x1(2:n,:)+A*x2(2:n,:)+B*x1(3:n+1,:)
  34       C*x2(n,:)+(A+B+D)*x1(n+1,:) ];
  35  yw=[ c*(C+D)*x1(1,:)+d*(A*x2(1,:)+B*x1(2,:))
  36       c*(C*x2(1:n-1,:)+D*x1(2:n,:))+d*(A*x2(2:n,:) ...
  37       +B*x1(3:n+1,:)) ];
  38  end

  39  nargout <= 1
  40    if ss(1) == 1, y = [y' yw'];
  41    else y = [y; yw]; end
  42  elseif ss(1) == 1,
  43    y = y'; yw = yw';
  44  end
\end{verbatim}
\ex

Line 1 declares {\tt kwt} to be a function that has at most two
input arguments and at most four output arguments.  The
case of having four output arguments occurs when {\tt x} is a matrix 
and the user wants the 2D transform to be separated out into four
submatrices.

Lines 2--5 check whether the argument {\tt k} is present.  If not,
its value is then taken from the global variable {\tt kw}.  To
take care of the case in which the global variable {\tt kw} may have 
been inadvertently erased, line 4 sets {\tt kw} to the default QS
transform.

Lines 6--12 examine the format of {\tt x} and choose the
various paths.  Line 7 changes {\tt x} to a column vector if it
is input as a row vector.  (Lines 42 and 45 convert
the output back into a row vector to match the input format.)
Line~8 checks whether {\tt x} is a 2D matrix.  If it is, the
work is delegated to the M-file {\tt kwt2.m}.  A user
has no need to use {\tt kwt2} directly; it is provided merely to
make the coding easier to understand.  We omit the listing of 
{\tt kwt2.m}.

Lines 13--17 check the format of the parameter vector {\tt k,}
supplying the default values for $ c $ and $ d $ if necessary.

Line 18 computes the components of $ \overline g $ {\tt = [ C, D, A, B ]}.
The names used to denote the components have
been retained from previous versions of the program.

Lines 21 and 22 extract the odd and even components of 
{\tt x}.  These are used to simplify the computation of the
$W$-transform.

Line 23 checks whether the length of {\tt x} is even.  If so,
lines 24--30 compute the transform.  Note that in the implementation
of the $W$-transform, we do not use matrix multiplication.  Lines 32--37
compute the transform for odd-length~{\tt x}.

Finally, lines 39--44 put the output in the correct format
according to the number of output arguments requested and
the format of the input signal.

\newpage
\subsection{The Inverse $W$-transform \, {\tt ikwt}} 

\bx
\vspace*{1.2\baselineskip}
\begin{verbatim}
      %  y = ikwt(x)              default inverse W-transform 
      %  y = ikwt(x1,x2)          1D x1, x2
      %  y = ikwt(x1,x2,x3,x4)    2D matrices
      %  y = ikwt(x,k)            general inverse W-transform 

   1  function y=ikwt(x1,x2,x3,x4,k)

   2  ss=size(x1); 

   3  if nargin == 1
   4    if ss(1) == 1, x = x1(:); 
   5    elseif ss(2) > 1 & nargin == 1, y = ikwt2(x1); return
   6    else x = x1;
   7    end
   8  elseif nargin == 2
   9    if min(ss) > 1
  10      y = ikwt2(x1,x2); return
  11    elseif size(x2,1) == 1 & size(x2,2) < 7 & size(x2,2) > 3
  12      x = x1(:); k = x2;
  13    else x = [x1(:); x2(:)];
  14    end
  15  elseif nargin == 3
  16    x = [x1(:); x2(:)]; k = x3;
  17  elseif nargin == 4
  18    y = ikwt2([x1 x2; x3 x4]); return
  19  elseif nargin == 5
  20    y = ikwt2([x1 x2; x3 x4],k); return
  21  end

  22  if exist('k') ~=1
  23    global kw
  24    if exist('kw') ~= 1, k = [1 3 3 1]/4; else k = kw; end
  25  end

  26  if length(k) == 4, c=1; d=-1;
  27  elseif length(k) == 5, c=k(5); d=-1;
  28  elseif length(k) == 6, c=k(5); d = k(6); 
  29  else   k = [1 3 3 1]/4; c=1; d=-1;
  30  end
\end{verbatim}
\ex
\bx
\begin{verbatim}
  31  C=-k(4)/c; D=k(3)/c; A=-k(2)/d; B=k(1)/d; 
  32  DD = (A*D-B*C); 

  33  s=size(x); n=floor(s(1)/2);

  34  if s(1) == n*2
  35    x1=x(1:n,:);
  36    x2=x(n+1:s(1),:);
  37    y1=[ (d*x1(1,:)-x2(1,:))/(C+D)
  38         (C*(c*x1(1:n-1,:)-x2(1:n-1,:))...
  39         +A*(d*x1(2:n,:)-x2(2:n,:)))/DD ]/(d-c);
  40    y2=[ (D*(-c*x1(1:n-1,:)+x2(1:n-1,:))...
  41         +B*(-d*x1(2:n,:)+x2(2:n,:)))/DD
  42         (-c*x1(n,:)+x2(n,:))/(A+B) ]/(d-c);
  43    y=zeros(s);
  44    y(1:2:2*n,:)=y1; y(2:2:s(1),:)=y2;

  45  else
  46    x1=x(1:n+1,:);
  47    x2=x(n+2:s(1),:);
  48    y1=[ (d*x1(1,:)-x2(1,:))/(C+D)
  49         (C*(c*x1(1:n-1,:)-x2(1:n-1,:))...
  50         +A*(d*x1(2:n,:)-x2(2:n,:)))/DD ]/(d-c);
  51    DDD = A*D+A^2+A*B-B*C;
  52    y1=[y1; (C*(c*x1(n,:)-x2(n,:))/(d-c)+A*x1(n+1,:))/DDD];
  53    y2=(D*(-c*x1(1:n-1,:)+x2(1:n-1,:))...
  54         +B*(-d*x1(2:n,:)+x2(2:n,:)))/DD/(d-c);
  55    y2=[y2; (-(D+A+B)*(c*x1(n,:)-x2(n,:))/(d-c)...
  56         -B*x1(n+1,:))/DDD];
  57    y=zeros(s);
  58    y(1:2:s(1),:)=y1; y(2:2:s(1),:)=y2;
  59  end

  60  if ss(1) == 1, y = y'; end
\end{verbatim}
\ex

A remarkable and useful property of a $W$-matrix is that its inverse
is also generated by a pair of vectors of length 4.  Furthermore,
all $W$-matrices of even, or odd, sizes have the same structure; they
differ only in the number of pairs of basic vectors used.

Using MAPLE, we are able to obtain the formula for the
inverse matrix of a general $W$-matrix.  The algorithm for computing inverse $W$-transform 
consists of merging {\tt y1} and {\tt y2} by interlacing their components
and then multiply by the inverse $W$-matrix.

To make {\tt ikwt} more user-friendly, we allow the input of 1D
signals either as two separate vectors {\tt y1} and {\tt y2} or as one
single vector obtained by appending {\tt y2} to {\tt y1,} and the input
of 2D images either as four separates matrices or as one single
matrix.  The parameter
vector	{\tt k} may be supplied as an argument or omitted ({\tt kw} 
is used as default).  Ambiguity arises and {\sl CAUTION} 
must be exercised when there are two
input arguments each of length between 4 and 6, since the
second argument can be interpreted either as {\tt y2} or as {\tt k}.
In such a case, the latter interpretation prevails.


Lines 3--21 examine the number of input arguments and
choose the appropriate path, in particular, delegating the work
to {\tt ikwt2} if there is indication that the input signal
is a matrix.  The listing of {\tt ikwt2} is omitted.

Lines 22--32 check the existence and format of {\tt k} 
and extract from it the coefficients used in the inverse transform.

Lines 34--44 perform the inverse transform on even-lengthed
signals and lines 46--58 on odd-lengthed signals.
The formulas used are slightly more complicated than
the corresponding ones for the forward transform.  This is because
the inverse matrix has a factor the reciprocal of the
determinant of the $W$-matrix and components in
both the first two rows and 
the last two rows have rather involved expressions.
Other than that, the formulas are straightforward.
Line 60 puts the output in a row vector if the input
signal is a row vector.

\newpage
\subsection{The $ W $-Multiresolution Analysis \, {\tt wma}} 

The command {\tt wma} carries out a multiresolution analysis on
a signal {\tt x} up to {\tt n} levels using the parameter vector {\tt k}.
As usual, the default {\tt k} is taken from {\tt kw.}  The default
{\tt n} is 3. The arguments {\tt n} and {\tt k} can be input in either
order.

\bx
\begin{verbatim}
      %  y = wma(x)         default W-multiresol. anal., n = 3
      %  y = wma(x,n)       n-level analysis
      %  y = wma(x,n,k)     general W-multiresol. anal.
      %
      %  [y, ind] = wma(x)  ind can be used to separate out the
      %                       wavelet vectors

   1  function [y, ind] = wma(x,n,k)

   2  if nargin == 1, n = 3; 
   3  elseif nargin ==2, if length(n) > 1, k = n; n = 3; end
   4  else  if length(k)==1, tmp = n; n = k; k = tmp; end
   5  end

   6  if exist('k') ~=1
   7    global kw, 
   8    if exist('kw') ~= 1, k = [1 3 3 1]/4; else k = kw; end
   9  end

  10  ss = size(x); ind = ss;

  11  for ii = 1:n
  12    y(1:ss(1),1:ss(2)) = kwt(x,k);
  13    ss = ceil(ss/2); ind = [ss; ind];
  14    x = y(1:ss(1),1:ss(2));
  15  end
\end{verbatim}
\ex

Line 1 indicates that the function has an optional output {\tt ind},
which gives the indices marking the ends of the various 
segments of the transformed signal at all the levels.
For instance, if {\tt ind} is the matrix
\begin{quote}
\verb+[ 1 13; 1 26; 1 51; 1 101 ]+,
\end{quote}
then the first-level {\tt y2} is made up of the 52-nd to
the 101-st components, the second-level {\tt y2} is of the 27-th
to the 51-st components, the third-level {\tt y2} of the 14-th
to the 26-th components, and finally the third-level {\tt y1} 
of the 1-st to the 13-th components.
That the first column of {\tt ind} are all 1 indicates that
the signal is one dimensional.  If {\tt x} is a matrix, the first column
of {\tt ind} will give the corresponding indices in the vertical
direction.

Lines 2--5 check the number and format of the input arguments,
figuring out which argument is {\tt n} and which is {\tt k} and
supplying the default value for {\tt n} if needed.  Lines 6--9
supply the default value for {\tt k} if needed.

Lines 11--16 form the main loop that calls {\tt kwt} to perform the
transform and then find the size and
signal to be transformed at the next level.

\newpage
\subsection{The Inverse $ W $-Multiresolution Analysis \, {\tt iwma}} 

The inverse $ W $-multiresolution analysis is similar to the forward analysis except
that {\tt ikwt} is invoked instead of {\tt kwt}, and the highest
level is dealt with first.
The arguments {\tt n} and {\tt k} can be input in either order.

The various wavelet vectors{/}submatrices must be first 
assembled into one single signal {\tt x} before input into
the function.

\bx
\begin{verbatim}
        %  y = iwma(x)       default inverse W-multiresol. transf.
        %                      n = 3
        %  y = iwma(x,n)     n-level default inverse W-ma
        %  y = iwma(x,n,k)   general inverse W-ma

     1	function y = iwma(x,n,k)

     2	if nargin == 1, n = 3; 
     3	elseif nargin == 2, if length(n) > 1, k = n; n = 3; end
     4	else  if length(k)==1, tmp = n; n = k; k = tmp; end
     5	end

     6	if exist('k') ~=1
     7	  global kw, 
     8	  if exist('kw') ~= 1, k = [1 3 3 1]/4; else k = kw; end
     9	end

    10	ss(n,:) = size(x);
    11	for ii = n-1:-1:1
    12	  ss(ii,:) = ceil(ss(ii+1,:)/2);
    13	end

    14	for ii = 1:n
    15	  x(1:ss(ii,1),1:ss(ii,2)) = ...
    16      ikwt(x(1:ss(ii,1),1:ss(ii,2)),k);
    17	end
    18	y = x;
\end{verbatim}
\ex

\newpage
\subsection{Extracting Multiresolution Analysis Wavelets \, {\tt maw}} 

\bx
\vspace*{1.2\baselineskip}
\begin{verbatim}
      %  z = maw(y,j)     j-th level wavelet vector for 1D y
      %  z = maw(y,j,1)   (top) j-th level y1 vector for 1D y
      %
      %  z = maw(y,j,h)   h-th submatrix (out of 4) at j-th level
      %                     for 2D y

   1  function z = maw(y,j,h)

   2  ind = size(y);

   3  if min(ind) == 1
   4    ind = length(y);
   5    for i = 1:j-1, ind = ceil(ind/2); end
   6    if nargin == 2, z = y(ceil(ind/2)+1:ind); 
   7    else z = y(1:ceil(ind/2)); end
   8  else

   9    for i = 1:j-1, ind = ceil(ind/2); end
  10    ind0 = ceil(ind/2);
  11    if ~exist('h'), h = 1; end

  12    if h == 1,     z = y(1:ind0(1),1:ind0(2));
  13    elseif h == 2, z = y(1:ind0(1),ind0(2)+1:ind(2));
  14    elseif h == 3, z = y(ind0(1)+1:ind(1),1:ind0(2));
  15    else z = y(ind0(1)+1:ind(1),ind0(2)+1:ind(2));
  16    end
  17  end
\end{verbatim}
\ex

For one-dimensional signals, we need only one index to specify
the level at which to extract the wavelet. Thus {\tt maw(y,3)} 
gives the level 3 wavelet.  Supply a second index if one
requires the top-level {\tt y1} vector.

For two-dimensional signals, we need one index to specify the
level and another index to specify which of the four
submatrices is needed.

Lines 3--7 handle the 1D case and lines 8--17 the 2D case.

\newpage
\subsection{Interactive Multiresolution Analysis \, wmai}

\bx
\vspace*{1.2\baselineskip}
\begin{verbatim}
      %%% wmai.m

   1  SIGNAL = input('Enter signal to be analyzed: ','s');
   2  LEVELS = input('Enter the number of levels : ');
   3  if length(LEVELS) == 0, LEVELS = 3; end
   4  SIZES = 0; eval([ 'SIZES = size(' SIGNAL ');' ])
   5  if min(SIZES) == 1, 
   6    eval([ '[W' SIGNAL '11, W' SIGNAL '12] = ...
   7         kwt(' SIGNAL ');' ] )
   8  else
   9    eval([ '[W' SIGNAL '11, W' SIGNAL '12, W' SIGNAL '13,  ...
  10         W' SIGNAL '14] = kwt(' SIGNAL ');' ] )
  11  end
  12  for ii = 2:LEVELS
  13    II = int2str(ii);
  14    if min(SIZES) == 1, 
  15      eval([ '[W' SIGNAL II '1, W' SIGNAL II '2] = ...
  16           kwt(W' SIGNAL int2str(ii-1) '1);' ] )
  17    else
  18      eval([ '[W' SIGNAL II '1, W' SIGNAL II '2, ...
  19           W' SIGNAL II '3, W' SIGNAL II '4] = ...
  20           kwt(W' SIGNAL int2str(ii-1) '1);' ] )
  21    end
  22  end
\end{verbatim}
\ex

The MATLAB command {\tt eval} is used to concatenate strings
to construct variable names for the wavelet vector{/}matrix.

A disadvantage of using a script, such as {\tt wmai},
instead of a function is that a few variables, such as {\tt SIGNAL} 
and {\tt LEVELS}, will be left dangling in the workspace.

{\sl Beware of conflicting variable names.}
\newpage
\subsection{Visualization \, {\tt wplot}\, and \, {\tt lbar}} 

Using {\tt wplot}, one can visualize
the top-level {\tt y1} and the {\tt y2} of all levels 
in the traditional plot (Figure 3) used by 
many wavelet researchers.  Built into {\tt wplot} is the same
procedure used by {\tt maw} to carve out the various
segments from the input signal (lines 13--16).  Each
segment is then plotted in a subplot using the subroutine
{\tt lbar}.

\bx
\begin{verbatim}
      %  wplot(x)        plot 3-level multiresol. anal.
      %  wplot(x,n)      plot n-level multiresol. anal.
      %  wplot(x,n,r)    use r < 1 as magnification factor
      %  wplot(x,n,r,c)  use color c

   1  function wplot(x,n,r,c)

   2  if exist('n') ~= 1, n = 3; r = 1; 
   3  elseif exist('r') ~=1, if n > 1, n = ceil(n); r = 1;
   4                         else r = n; n = 3; end
   5  else if n > 1, n = ceil(n); 
   6       else tmp = r; r = n; n = max(2,tmp); end
   7  end
   8  if exist('c') ~= 1, c = 'y'; end

   9  ind = length(x); 
  10  mx = min(min(x),0); Mx = max(max(x),0);
  11  figure(gcf)

  12  for ii = n+1:-1:1
  13    if ii > 2, xx = x(ceil(ind/2)+1:ind)/r^(ii-2);
  14    elseif ii > 1, xx = x(ceil(ind/2)+1:ind);
  15    else xx = x(1:ind); end
  16    ind = ceil(ind/2);

  17    h = subplot(n+1,1,ii);
  18    lbar(xx,c)
  19    axis([ 1 length(xx) mx Mx ])
  20    axis off
  21    rect = get(h,'pos');
  22    set(h,'pos',[ rect(1) rect(2)-(n+1-ii)*rect(4)/10 ...
  23        rect(3) rect(4) ])
  24  end
\end{verbatim}
\ex

Lines 2--8 check the input arguments.  The argument {\tt n} 
specifies the number of levels (defaults to 3) and argument
{\tt r} is a magnification factor used for lower level wavelets;
{\tt n} and {\tt r} can be input in either order.  The argument {\tt c} 
specifies the color (defaults to yellow) of the plot and
can be input only after {\tt n} and {\tt r}.

Line 10 determines the lower and upper limits on the
vertical axis used by each subplot. Line 11 raises the
current figure to the foreground.

Lines 17--20 do the plotting and adjust the axis limits.
Lines 21--23 move the subplots closer together for
aesthetic reasons.

By convention, a wavelet plot is constructed by drawing a 
vertical line at each integral point on the horizontal
with height proportional to the appropriate component of the
wavelet vector.  In essence, it is a bar graph in which each bar
is a vertical line.  MATLAB does not seem to have a
command to draw such a plot, and we therefore write our own.
Expecting it to be of interest in other contexts, we
make the command {\tt lbar} slightly more general.  It can
take two optional arguments, one to specify the range to be
labeled on the horizontal axis, and another
to specify the color of the graph.

\bx
\begin{verbatim}
      %  lbar(y)        plot y as line-bar graph
      %  lbar(x,y)      plot y VS x as line-bar graph
      %  lbar(x,y,c)    use color c

   1  function lbar(x,y,c)

   2  n = length(x);

   3  if nargin == 1, y = x(:); x = (1:n)'; c = 'y';
   4  elseif nargin == 2, 
   5    if length(y) == 1, c = y; y = x(:); x = (1:n)';
   6    else x = x(:); y = y(:); c = 'y'; end
   7  end

   8  n = 3*n;
   9  xx = zeros(n,1); yy = xx;
  10  xx(1:3:n) = x; xx(2:3:n) = x; xx(3:3:n) = x;
  11  yy(2:3:n) = y;
  12  plot(xx,yy,c)
\end{verbatim}
\ex

\newpage
\subsection{Compression Utilities \, {\tt largesta} and \, {\tt quant}} 

\bx
\vspace*{1.2\baselineskip}
\begin{verbatim}
      %   y     = largesta(n,x)   n > 1, largest (absolute value)
      %                                  n components of x
      %  [y, N] = largesta(r,x)   r < 1, largest fraction r
      %                           N = number of components retained

   1  function [ y, N ] = largesta(n,x)

   2  if length(n) > 1, tmp = n; n = x(1); x = tmp; end
   3  ax = abs(x); 
   4  M = max(max(ax)); 
   5  a = 0; b = M; n1 = prod(size(x(:))); n2 = nnz(ax == M);

   6  if  n == 1, y = x; n = n1; return
   7  elseif n <= 0, y = 0*x; n = 0; return
   8  elseif n < 1, n = min(n1-1,round(n1*n)); 
   9  end

  10  while n1 ~= n & n2 ~= n & n1 ~= n2 & b-a > 0.000001*M
  11    c = (a+b)/2;
  12    nc = nnz(ax >= c);
  13    if nc > n, a = c; n1 = nc;
  14    else       b = c; n2 = nc;
  15    end
  16  end
  17  if n1 == n
  18    small = find(ax < b);
  19    % y = x.*(ax > b);
  20  else
  21    small = find(ax < c);
  22    % y = x.*(ax > c);
  23  end

  24  y = x;
  25  y(small) = zeros(size(small));
  26  N = nnz(y);
\end{verbatim}
\ex

The command {\tt largesta} retains the largest (in absolute values)
{\tt n} elements of the matrix {\tt x} and truncates the rest to 0.
If {\tt n} is larger than 1, it is used as the actual count of
elements to be retained.  If {\tt n} is less than or equal to 1,
it is interpreted as the fraction of the total number of 
elements to be retained. 
The number of elements actually retained may differ
from the requested number.  The algorithm uses a bisection
procedure.  A sister command {\tt largest} retains the largest
(in algebraic value) elements.

The command {\tt quant(x,n)} quantizes the components to take
values in the set $ \left\{0,\pm n,\pm 2n,\pm 3n,\cdots \right\} $.
If {\tt n} is a string containing a number, quant will
produce approximately {\tt n} quantization levels.
If {\tt n} is not specified, approximately ten quantization 
levels will be used.

\bx
\begin{verbatim}
      %    y    = quant(x)      quantize x (approximately 10 levels)
      %    y    = quant(x,q)    use q as quantization interval
      %  [y, Q] = quant(x,'n')  use approximately n levels
      %                         Q = actual quantization interval

   1  function [y, q] = quant(x,q)
      
   2  if nargin == 1
   3    q = '10';
   4  elseif length(x) == 1 | isstr(x)
   5    tmp = x; x = q; q = tmp;
   6  end

   7  if isstr(q)
   8    eval([ 'q = ' q ';' ])
   9    r = (max(max(x)) - min(min(x)))/q*1.5;
  10    p = floor(log(r)/log(10));
  11    q = floor(r/10^p)*10^p;
  12  end

  13  y = round(x/q)*q;
\end{verbatim}
\ex

In principle, neither {\tt largesta} nor {\tt quant} is the ultimate
utility for automating compression.  The ideal utility can automatically
determine the cutoff threshold, such as the article {\tt n} or {\tt r} to
{\tt largesta} and the article {\tt q}  to {\tt quant}, based upon a specified
compression ratio or compression performance.

\newpage
\subsection{Orthogonal Compensation \, {\tt oc}} 

\def\yo{\mbox{\tt y1}}
\def\yt{\mbox{\tt y2}}
\def\bg{\mbox{\bf v}}
\def\bh{\mbox{\bf w}}

The $ W $-transform associated with a $W$-matrix
$ W $ can be summarized in the equations
\begin{equation}    \mbox{\tt y} = \mbox{\bf W} \mbox{\tt x} , \quad  \mbox{\tt x} = \mbox{\bf W}^{-1} \mbox{\tt y} .  \Label{qs}  \end{equation}
The odd and even components of {\tt y} form the pair
of vectors
\begin{equation}    \yo = [ y_{11},\, y_{12},\, \cdots  ]' , \quad  \yt = [ y_{21},\, y_{22},\, \cdots  ]' . \end{equation}
Let us denote the columns of the matrix $ \mbox{\bf W}^{-1} $ as
\begin{equation}     \left[ \begin{array}{cccccccccc}    \bg_1 &  \bh_1 &  \bg_2 &  \bh_2 &  \cdots   \end{array} \right ] . \end{equation}
Then the second equation in (\ref{qs}) has the equivalent form
\begin{equation}    \mbox{\tt x} = ( y_{11}\bg_1 + y_{12}\bg_2 +\cdots   ) +  ( y_{21}\bh_1 + y_{22}\bh_2 +\cdots   ) .  \Label{qs2}  \end{equation}
This equation suggests that the $ W $-transform can be interpreted as
the decomposition of {\tt x} along the subspaces $ G $ and
$ H $, spanned by $ \bg $ and $ \bh $, respectively.  

In the analogous interpretation of the
Haar and Daubechies $ D_4 $  
transforms, the linear subspaces $ G $ and $ H $ are orthogonal
to each other.  In addition, the one-dimensional subspaces
generated by all the $ \bg_i $ and $ \bh_i $ are mutually
orthogonal.  When some of the components in $ \yt $ are 
discarded, the compressed vector is then the unique
signal, in the space spanned by the remaining base vectors,
that best approximates the original signal.

For a general $ W $-transform, $ G $ and $ H $ are not necessarily
orthogonal.  This fact seems to argue against the use
of general $ W $-transforms.  In practice, a reasonable signal (one that
is not wildly oscillating or badly degraded by noise)
usually has such small coefficients in the $ H $ subspace decomposition that 
even if we do not take additional steps to optimize the
approximation, the error incurred in simply discarding
them is negligible.  This is true, in particular, for the QS transform.

We give below the method of orthogonal compensation to enhance
the approximation when discarding some of the components of
$ \yt $.  Let {\bf d} be the vector to be discarded.  It is
likely to be a partial sum of the expression in the
second pair of parentheses in (\ref{qs2}).  We decompose {\bf d}
into a linear combination of the vectors $ \bg_i $ and
an error vector that is orthogonal to $ G $:

\begin{equation}    \mbox{\bf d} = ( a_1\bg_1 + a_2\bg_2 +\cdots  ) + \mbox{\bf e}.  \Label{oc}  \end{equation}

After determining $ a_i $, they are added
to the corresponding $ y_{1i} $, so that the actual part that
is discarded is {\bf e}, which is orthogonal to~$ G $.
To this end, we take inner products of {\bf d} with each
of $ \bg_i $.  One can easily verify that $ a_i $ is the
solution to the tridiagonal system of linear equations
\begin{equation}     \left( \begin{array}{cccccccccc}    \alpha _{11} & \alpha _{12} \\ \alpha _{21} & \alpha _{22} & \alpha _{23} \\   & \alpha _{32} & \alpha _{33} & \alpha _{34} \\   &   &   & \ddots  \\   &   &   &   & \ddots  \end{array} \right )  \left( \begin{array}{cccccccccc}    a_1 \\ a_2 \\ a_3 \\ \vdots \\ \vdots  \end{array} \right ) =  \left( \begin{array}{cccccccccc}    <\mbox{\bf d},\bg_1>  \\  <\mbox{\bf d},\bg_2>  \\ <\mbox{\bf d},\bg_3>  \\  \vdots \\  \vdots  \end{array} \right ) , \end{equation}
where $ \alpha _{ij}=<\bg_i,\bg_j> $.
The matrix in the equation is square, with width equal to the length
of $ \yo $.  

\bx
\begin{verbatim}
      %  y = oc(y1,ye)      orthogonal compensation, default W-tranf.
      %  y = oc(y1,ye,k)    general W-transform

   1  function y = oc(y1,ye,k)

   2  if nargin == 2
   3    global kw
   4    if exist('kw') ~= 1, k = [1 3 3 1]/4; else k = kw; end
   5  end

   6  if length(k) == 4, c=1; d=-1;
   7  elseif length(k) == 5, c=k(5); d=-1;
   8  elseif length(k) == 6, c=k(5); d = k(6); 
   9  else   k = [1 3 3 1]/4; c=1; d=-1;
  10  end

  11  k1=k(1); k2=k(2); k3=k(3); k4=k(4);
  12  h1=-k1/d; h2=-k2/d; h3=-k3/c; h4=-k4/c;
\end{verbatim}
\ex
\bx
\begin{verbatim}
  13  s1 = size(y1); se = size(ye);
  14  y1 = y1(:); ye = ye(:);
  15  m = length(y1); n = length(ye);

  16  A = k1*k3+k2*k4;
  17  B = k1^2+k2^2+k3^2+k4^2;
  18  C = k3*k1+k4*k2;
  19  B1 = (k2*k3-k1*k4)^2/(k2-k1)^2 +k3^2+k4^2;
  20  Bn = k1^2+k2^2+(k2*k3-k1*k4)^2/(k4-k3)^2;

  21  D = k1*h3+k2*h4;
  22  E = k1*h1+k2*h2+k3*h3+k4*h4; 
  23  F = k3*h1+k4*h2; 
  24  E1 = -((k2*k3-k1*k4)/(k2-k1)/c)^2*d+k3*h3+k4*h4; 

  25  if m == n
  26    En = k1*h1+k2*h2-((k2*k3-k1*k4)/(k4-k3)/d)^2*c; 
  27    z = ([ 0;D*ye(1:n-1) ] ...
  28         + [ E1*ye(1);E*ye(2:n-1);En*ye(n) ] ...
  29         + [ F*ye(2:n);0 ]);
  30    GH = toeplitz([ B A zeros(1,n-2) ],...
  31         [ B C zeros(1,n-2) ]);
  32    GH(1,1) = B1;  GH(n,n) = Bn;

  33  else
  34    r = -h2*h3+h1*h4;
  35    r = (r/(r-h1*h2+h2^2))^2;
  36    Dn = (c-d)*r*( (h1-h2+h3)*h1 - h2*k4 );
  37    En = k1*h1+k2*h2-r*c*( (h1-h2+h3)^2 + (k4)^2 );
  38    z = ([ 0;D*ye(1:n-1);Dn*ye(n) ] ...
  39         + [ E1*ye(1);E*ye(2:n-1);En*ye(n);0 ] ...
  40         + [ F*ye(2:n);0;0 ]);
  41    GH = toeplitz([ B A zeros(1,n-1) ],...
  42         [ B C zeros(1,n-1) ]);
  43    GH(1,1) = B1;  
  44    GH(n,n) = k1^2+k2^2+r*c^2*( (h1-h2+h3)^2 + (k4)^2 );
  45    GH(n+1,n+1) = r*(d-c)^2*(h1^2+h2^2);
  46    GH(n,n+1) = -c*Dn;
  47    GH(n+1,n) = GH(n,n+1);
  48  end

  49  y = y1+GH\z;
  50  if s1(1) ==1, y = y'; end
\end{verbatim}
\ex
