\setcounter{equation}{0}\section{Examples of Numerical Experiments}

An example not covered in this section is in
the program obtainable by {\tt ftp} as described in Section~5.
The  M-file {\tt ldemo.m} 
gives the data and explanation of 
a 10-fold compression of Lena.  With such a high compression
ratio, distortion is inevitable.  But the restored image 
contains fewer offensive visual artifacts such as blocking
effects.  This suggests that $W$-transforms may be an alternative
to DCT in achieving high ratio compression.

\subsection{{\tt compare}} 

Our first
experiment compares the performance of two $W$-transforms in 
compressing a sample signal {\tt x}.
The signal {\tt x} is transformed into
two vectors $ \yo $ and $ \yt $.  An approximate signal
is restore by retaining only some of the components
of $ \yt $.  The discrepancies between the original and
the approximate signals are then computed using the $ L^2 $
norm. To prepare for the experiment, one
chooses a signal {\tt x} and the $W$-transforms and then 
invokes {\tt compare.}
\begin{quote}
\verb+>> x = +
{\sl DEFINITION OF SIGNAL}
\LEFTLINE 
\verb+>> kw = [ 1 3 3 1 ]/4; kww = kwdau; compare+
\end{quote}

A sample output is as follows.
\def\ex{\vspace{-0.5\baselineskip}\endbrace
  \wd0=0.6\bxleng \fbox{\hskip0.7em\box0}}
\vspace*{0.5\baselineskip}

\LEFTLINE 
\bx
\begin{verbatim}
           W1OC        W1         W2

  # Coef             errors

     2   0.001031   0.001107   0.015413
     5   0.000928   0.001014   0.013358
     7   0.000818   0.000882   0.012025
    ..   ........   ........   ........
    35   0.000185   0.000207   0.001881
    37   0.000171   0.000191   0.001533
    40   0.000131   0.000148   0.000991
    ..   ........   ........   ........
\end{verbatim}
\ex
\vspace*{\baselineskip}
%     10   0.000729   0.000807   0.010355
%     32   0.000238   0.000266   0.002388
%     13   0.000653   0.000724   0.008730
%     27   0.000349   0.000388   0.003674
%     31   0.000259   0.000290   0.002729
%     16   0.000560   0.000624   0.007541
%     17   0.000536   0.000595   0.006775
%     20   0.000449   0.000491   0.005994
%     22   0.000422   0.000474   0.004916
%     25   0.000368   0.000411   0.004317
%     43   0.000087   0.000096   0.000736
%     45   0.000049   0.000055   0.000342
%     47   0.000018   0.000020   0.000118
%     50   0.000000   0.000000   0.000000

\def\ex{\vspace{-0.5\baselineskip}\endbrace
  \wd0=0.95\bxleng \begin{center}\fbox{\hskip0.7em\box0}\end{center}}

The first column gives the number of components of $ \yt $ retained
(these numbers can be controlled by changing the variable {\tt RR}).
The second column gives the errors when the signal
is compressed using the first $W$-transform with 
orthogonal compensation.
The third and fourth columns give the errors when the
signal is compressed using the two $W$-transforms without
orthogonal compensation.
The sample signal is the same signal {\tt x} used in Section~2,
and the $W$-transforms compared are the QS and $ D_4 $ transforms.
As seen in the above example, retaining just 2 components
of $ \yt $ in the QS transform is better than retaining 37
components in the $ D_4 $  transform.

The listing of the m-file follows.

\bx
\begin{verbatim}
      %%% compare.m

   1  if ~exist('RR'), RR = 0.05:0.05:0.5; end

   2  normx = norm(x);
   3  [y1 y2] = kwt(x);
   4  [d1 d2] = kwt(x,kww);

   5  disp(' ')
   6  disp('           W1OC        W1          W2')      
   7  disp(' ')
   8  disp('  # Coef             errors')
   9  disp(' ')

  10  for r=RR
  11    cy2 = largesta(r,y2); cx = ikwt(y1,cy2);     cE=cx-x;
  12    oy1 = oc(y1,y2-cy2);  ox = ikwt(oy1,cy2);    oE=ox-x;
  13    cd2 = largesta(r,d2); dx = ikwt(d1,cd2,kww); dE=dx-x;

  14    fprintf('   %3d   %f   %f   %f\n',...
  15    nnz(cy2),norm(oE)/normx,norm(cE)/normx,norm(dE)/normx)
  16  end
\end{verbatim}
\ex

\newpage
\subsection{Decompression \, {\tt decp}} 

The decompression function first
compresses the $W$-transformed output matrices {\tt Y1}, 
{\tt Y2}, {\tt Y3}, and {\tt Y4} (obtained from an
earlier $W$-transform), either by
discarding small components or by quantizing the 
detail submatrices (according to whether the input
argument {\tt r} is less than or greater than 1), and then
reconstructs an approximate image.
The optional output arguments {\tt cr} and {\tt nz} 
give the {\sl compression ratio} and {\sl number of nonzero elements
retained}.  These terms do not have the precise meaning
as commonly used and are provided for reference only.

\bx
\begin{verbatim}
     1	function [X, cr, nz]= decp(Y1,Y2,Y3,Y4,r)

     2	s = size(Y1);

     3	if nargin == 1, Y2 = zeros(s); Y3 = Y2; Y4 = Y2; 
     4	elseif nargin ==5 & r <= 1 , 
     5	                    Y = [ 0*Y1 Y2; Y3 Y4 ];
     6	                    Y = largesta(r,Y);
     7	                    Y(s) = Y1;
     8	elseif nargin == 5, Y2 = quant(Y2,r);
     9	                    Y3 = quant(Y3,r);
    10	                    Y4 = quant(Y4,4*r);
    11	  nz = prod(size(Y1))+nnz(Y2)+nnz(Y3)+nnz(Y4);
    12	  Y = [Y1 Y2;Y3 Y4];
    13	end

    14	X = ikwt2(Y);
    15	cr = 4*s(1)*s(2)/nz;
\end{verbatim}
\ex

\newpage
\subsection{Image Compression \,  {\tt imcomp}} 

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

   1  if exist('IMCOMPR') ~= 1
   2  if ~exist('r3'), r3 = [ .01 .05 .1 ]; end
   3  fprintf('current r3 = ')
   4  fprintf(' %f ',r3)
   5  tr3 = input('Enter new r3: ','s');
   6  if length(tr3) > 0
   7    eval([ 'tr3 = [' tr3 '];']);
   8    if length(tr3) == 2
   9      eval([ 'r3(round(' num2str(tr3(1)) ')) = '...
  10             num2str(tr3(2)) ';'])
  11    else r3 = tr3;
  12    end
  13  end
  14  else
  15    r3 = IMCOMPR;
  16  end

  17  kwold = kw;
  18  kw = kwqs;

  19  LEVELS = length(r3);
  20  if LEVELS == 4
  21  [cWL31 cr4 nz4] = decp(WL41,WL42,WL43,WL44,r3(4));
  22  nz4 = nz4 - prod(size(WL41)); 
  23  disp([ 'Level 4 decompression done   nz = ' int2str(nz4) ])
  24  else
  25  cWL31 = WL31;
  26  end

  27  [cWL21 cr3 nz3] = decp(cWL31,WL32,WL33,WL34,r3(3));
  28  N3 = prod(size(WL31)); nz3 = nz3 - N3; 
  29  disp([ 'Level 3 decompression done   nz = ' int2str(nz3) ])

  30  [cWL11 cr2 nz2] = decp(cWL21,WL22,WL23,WL24,r3(2));
  31  nz2 = nz2 - prod(size(cWL21)); 
  32  disp([ 'Level 2 decompression done   nz = ' int2str(nz2) ])
\end{verbatim}
\ex
\bx
\begin{verbatim}
  33  [cL cr1 nz1] = decp(cWL11,WL12,WL13,WL14,r3(1));
  34  NN = prod(size(cWL11)); nz1 = nz1 - NN; 
  35  disp([ 'Level 1 decompression done   nz = ' int2str(nz1) ])

  36  if r3(1) < 1
  37    if LEVELS == 4
  38      nz = nz4+nz3+nz2+nz1+prod(size(WL41));
  39    else
  40      nz = nz3+nz2+nz1+prod(size(WL31));
  41    end

  42  cr = prod(size(L))/nz;
  43  fprintf('\nCompression Ratio: %d/%d = %f\n',prod(size(L)),nz,cr)

  44  TITLE = [ 'Compressed Image, CR = ' num2str(cr) '(' ...
  45    num2str(r3(1)) ', ' num2str(r3(2)) ', ' num2str(r3(3)) ];
  46  XLABEL = [ int2str(nz1) '   ' int2str(nz2) '   ' ...
  47            int2str(nz3) ];

  48  if LEVELS == 3
  49    TITLE = [ TITLE ')' ];
  50  else
  51    TITLE = [ TITLE ', ' num2str(r3(4)) ')' ];
  52    XLABEL = [ XLABEL '   ' int2str(nz4) ];
  53  end

  54  else
  55    cr = 8*NN/(N3*10 + nz3*(22-log(r3(3))/log(2)) +...
  56         nz2*(18-log(r3(2))/log(2)) + nz1*(14-log(r3(1))/log(2)));
  57    fprintf('\ncr = %f\n',cr)

  58  TITLE = [ 'Quantized Image, cr = ' num2str(cr) '  (' ...
  59    int2str(r3(1)) ', ' int2str(r3(2)) ', ' int2str(r3(3)) ')' ];
  60  XLABEL = [ int2str(nz1) '   ' int2str(nz2) '   ' int2str(nz3) ];
  61  end

  62  image(cL), colormap(gray(256))
  63  axis image
  64  title(TITLE)
  65  xlabel(XLABEL)
\end{verbatim}
\ex

Our third  experiment uses a 256-level grayscale image
stored in a matrix~{\tt L}.  The command {\tt wmai} should
be used first to produce a three-
or four-level multiresolution analysis of {\tt L}.
One then invokes {\tt incomp}.
One will be prompted to input {\tt r,} comprising
three or four values to be used
as the fractions of components of the detail matrices to be
retained in each of the 
levels.  Experience shows that the
first-level details are less important, however, hence, only a
relatively small fraction of the components need be kept 
to give a satisfactory approximate image.  The
suggested default is {\tt r = 0.01, 0.05, 0.1} using three
levels.

The compression ratio displayed by the program
is actually based on the ratio of the size
of the original image to the number of nonzero
elements retained in the compressed representation,
adjusted by the qunatization level.
The QS transform, for one, leads to transformed
matrices with elements much larger than the original
matrix.  Hence, a precise bit rate count of the compressed data
has to be performed to give an accurate compression
ratio.  On the other hand, we have not used
entropy compression, such
as Huffman coding or arithmetic coding, that can 
further increase the compression ratio.

\newpage
\subsection{Finding the Optimal $W$-transform \, {\tt opt}} 

Our fourth experiment compresses a signal {\tt x} by taking the
$W$-transform specified by {\tt kww = [ 1 r r 1 ]} 
and then discarding all except
the largest (in absolute value) 4 components of the resulting
{\tt y2}.  The value of {\tt r} is varied (as specified by
the vector {\tt RRopt}) and the corresponding errors
printed.  
\bx
\begin{verbatim}
   1  normx = norm(x);
   2  if ~exist('RRopt'), RRopt = 2.5:0.05:3.5; end
   3    disp(' ')
   4    disp(' k vector used is [ 1 r r 1 ]')
   5    disp(' ')
   6    disp('    r         WOC         W')
   7    disp(' ')
   8  for r=RRopt
   9    kww = [1 r r 1];
  10    [y1 y2] = kwt(x,kww);
  11    cy2 = largesta(4,y2);
  12    cx = ikwt(y1,cy2,kww);
  13    cE=cx-x;
  14    oy1 = oc(y1,y2-cy2,kww);
  15    ox = ikwt(oy1,cy2,kww);
  16    oE=ox-x;
  17    fprintf('%f   %f   %f\n',r,norm(oE)/normx,norm(cE)/normx)
  18  end
\end{verbatim}
\ex

The following output shows that the optimal {\tt r} is about $ 2.95 $.
\def\ex{\vspace{-0.5\baselineskip}\endbrace
  \wd0=0.5\bxleng \fbox{\hskip0.7em\box0}}
\vspace*{0.48\baselineskip}

\LEFTLINE 
\bx
\begin{verbatim}
    r        WOC          W
........   ........   ........
2.850000   0.002690   0.002926
2.900000   0.001450   0.001571
2.950000   0.000277   0.000294
3.000000   0.000984   0.001060
3.050000   0.002113   0.002270
........   ........   ........
\end{verbatim}
\ex

To determine a more accurate result, change {\tt RRopt} to
{\tt 2.90:0.01:3.0.} and invoke {\tt opt} again.
