%DCMDMALA Demonstrates Mallat algorithm in the Wavelet Toolbox. 

%       M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
%	Last Revision: 01-Oct-97.
%       Copyright (c) 1995-98 by The MathWorks, Inc.
% $Revision: 1.5 $

%DEMALLAT Short demo about basic steps of FWT 1D.
%Non documented function, demo script file.

% Demo initialization ====================
if ~exist('DRVShowGUIFlag'), figNumber = 0; end;

echo off
clc;
titre   = 'Mallat Algorithm (FWT)';

if wshowdrv('init',figNumber)
        if figNumber==0
		active_fig = wdfigutl('create',titre);
        else
                set(figNumber,'Name',titre);
                wshowdrv('setactfig',figNumber,figNumber);
                democmdm('startshow',figNumber,'mala');
                return;
        end
end
if exist('wait_mode')~=1 , wait_mode = 1; end
if figNumber , active_fig = wshowdrv('getactfig',figNumber); end

comment = str2mat( ...
                '% Construct elementary original 1-D signal.', ...
                ' ', ...
                '        s = 2 + kron(ones(1,8),[1 -1]) + ...' ,...
                '            ((1:16).^2)/32 + 0.2*randn(1,16);' ...
                );
wshowdrv('disp',figNumber,comment);

% Construct elementary original 1D signal.
%        s = 0.5*[1 -1 1 -1]; s = [s s s s] + (1:16)/16;
       s = 2 + kron(ones(1,8),[1 -1]) + ((1:16).^2)/32 + 0.2*randn(1,16);


figure(active_fig)
h1 = subplot(311); plot(s,'r');
title('Original signal s.');
set(h1,'Xlim',[1 16])

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% For a given orthogonal wavelet,', ...
                '% compute the two associate', ...
                '% decomposition filters.', ...
                ' ', ...
                '        [Lo_D,Hi_D] = wfilters(''db1'',''d'');' ...
                );
wshowdrv('disp',figNumber,comment);

% For a given orthogonal wavelet,
% compute the two associate
% decomposition filters.
        [Lo_D,Hi_D] = wfilters('db1','d');

x_fil = 1:length(Lo_D);
figure(active_fig)
h3 = subplot(323); wdstem(h3,x_fil,Lo_D,'y');
xlabel('Lo_D.');
figure(active_fig)
h4 = subplot(324); wdstem(h4,x_fil,Hi_D,'y');
xlabel('Hi_D.');

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Lo_D is the decomposition low-pass filter.', ...
                '% Check for sum.', ...
                ' ', ...
                '        sm = sum(Lo_D);' ...
                );
wshowdrv('disp',figNumber,comment);

% Lo_D is the decomposition low-pass filter.
% Check for sum.
        sm = sum(Lo_D);

figure(active_fig)
subplot(323); xlabel(['Lo_D : sum = ', num2str(sm)]);
 
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Lo_D is the decomposition low-pass filter.', ...
                '% Check for square norm.', ...
                ' ', ...
                '        nrm = norm(Lo_D);' ...
                );
wshowdrv('disp',figNumber,comment);

% Lo_D is the decomposition low-pass filter.
% Check for square norm.
        nrm = norm(Lo_D);

figure(active_fig)
subplot(323); xlabel(['Lo_D : norm = ', num2str(nrm)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Compute approx. coef. of s using db1 : first step.', ...
                '% Convolve s and Lo_D the', ...
                '% decomposition low-pass filter.', ...
                ' ', ...
                '        tempo = conv(s,Lo_D);' ...
                );
wshowdrv('disp',figNumber,comment);

% Compute approx. coef. of s
% using db1 : first step.
% Convolve s and Lo_D the
% decomposition low-pass filter.

        tempo = conv(s,Lo_D);

figure(active_fig)
set(h1,'Xticklabel',[]);
h2 = subplot(312); plot(tempo);
title('s and Lo_D convolution.');
set(h2,'Xlim',[1 length(tempo)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Compute approx. coef. of s using db1 : second step.', ...
                '% Decimate convolution result', ...
                '% and keep central part.', ...
                ' ', ...
                '        ca1 = dyaddown(tempo);' ...
                );
wshowdrv('disp',figNumber,comment);

% Compute approx. coef. of s
% using db1 : second step.
% Decimate convolution result
% and keep central part.
        ca1 = dyaddown(tempo);

figure(active_fig)
h5 = subplot(325); plot(ca1);
xlabel('Approx. coef. : ca1.');
set(h5,'Xlim',[1 length(ca1)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

figure(active_fig)
h3 = subplot(323); wdstem(h3,x_fil,Lo_D,'y');
% xlabel('Lo_D.');
figure(active_fig)
h4 = subplot(324); wdstem(h4,x_fil,Hi_D,'y');
xlabel('Hi_D.');

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Hi_D is the decomposition high-pass filter.', ...
                '% Check for sum.', ...
                ' ', ...
                '        sm = sum(Hi_D);' ...
                );
wshowdrv('disp',figNumber,comment);

% Hi_D is the decomposition high-pass filter.
% Check for sum.
        sm = sum(Hi_D);

figure(active_fig)
subplot(324); xlabel(['Hi_D : sum = ', num2str(sm)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Hi_D is the decomposition high-pass filter.', ...
                '% Check for square norm.', ...
                ' ', ...
                '        nrm = norm(Hi_D);' ...
                );
wshowdrv('disp',figNumber,comment);

% Hi_D is the decomposition high-pass filter.
% Check for square norm.
        nrm = norm(Hi_D);

figure(active_fig)
subplot(324); xlabel(['Hi_D : norm = ', num2str(nrm)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Compute detail coef. of s using db1 : first step.', ...
                '% Convolve s and Lo_D the', ...
                '% decomposition high-pass filter.', ...
                ' ', ...
                '        tempo = conv(s,Hi_D);' ...
                );
wshowdrv('disp',figNumber,comment);

% Compute detail coef. of s
% using db1 : first step.
% Convolve s and Hi_D the
% decomposition low-pass filter.
        tempo = conv(s,Hi_D);

figure(active_fig)
h2 = subplot(312); plot(tempo);
title('s and Hi_D convolution.');
set(h2,'Xlim',[1 length(tempo)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Compute detail coef. of s using db1 : second step.', ...
                '% Decimate convolution result', ...
                '% and keep central part.', ...
                ' ', ...
                '        cd1 = dyaddown(tempo);' ...
                );
wshowdrv('disp',figNumber,comment);

% Compute detail coef. of s
% using db1 : second step.
% Decimate convolution result
% and wkeep central part.
        cd1 = dyaddown(tempo);

figure(active_fig)
h6 = subplot(326); plot(cd1);
xlabel('Detail coef. : cd1.');
set(h6,'Xlim',[1 length(cd1)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% This sequence is the one step', ...
                '% discrete wavelet transform of s.', ...
                '% M-file dwt performs it directly.', ...
                ' ', ...
                '        [ca1,cd1] = dwt(s,''db1'');' ...
                );
wshowdrv('disp',figNumber,comment);

% This sequence is the one step
% discrete wavelet transform of s.
% M-file dwt performs it directly.
        [ca1,cd1] = dwt(s,'db1');

figure(active_fig)
h5 = subplot(325); plot(ca1,'g');
xlabel('Approx. coef. : ca1.');
set(h5,'Xlim',[1 length(ca1)]);
figure(active_fig)
h6 = subplot(326); plot(cd1,'g');
xlabel('Detail coef. : cd1.');
set(h6,'Xlim',[1 length(cd1)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Now how to reconstruct ?', ...
                '% For a given orthogonal wavelet,', ...
                '% compute the two associate', ...
                '% reconstruction filters.', ...
                ' ', ...
                '        [Lo_R,Hi_R] = wfilters(''db1'',''r'');' ...
                );
wshowdrv('disp',figNumber,comment);

% Now how to reconstruct ?
% For a given orthogonal wavelet,
% compute the two associate
% reconstruction filters.
        [Lo_R,Hi_R] = wfilters('db1','r');

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

figure(active_fig)
h3 = subplot(323); wdstem(h3,x_fil,Lo_R,'y'); % plot(Lo_R);
xlabel('Lo_R.');
figure(active_fig)
h4 = subplot(324); wdstem(h4,x_fil,Hi_R,'y'); % plot(Hi_R);
xlabel('Hi_R.');
figure(active_fig)
h1 = subplot(321); plot(ca1,'g');
set(h1,'Xlim',[1 length(ca1)]);
title('Approx. coef. : ca1.');
figure(active_fig)
h2 = subplot(322); plot(cd1,'g');
title('Detail coef. : cd1.');
set(h2,'Xlim',[1 length(cd1)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
              '% Reconstruct detail : first step.', ...
              '% Upsample cd1 inserting zeros.', ...
              ' ', ...
              '        tempo = dyadup(cd1);' ...
                );
wshowdrv('disp',figNumber,comment);

% Reconstruct detail : first step.
% Upsample cd1 inserting zeros.
        tempo = dyadup(cd1);

figure(active_fig)
subplot(323); xlabel('');
subplot(324); xlabel('');
figure(active_fig)
h3 = subplot(313); plot(tempo);
xlabel('Upsampled det. coeff. cd1.');
set(h3,'Xlim',[1 length(tempo)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
              '% Reconstruct detail : second step.', ...
              '% Convolves upsampled detail coeff.', ...
              '% and Hi_R reconstruction high-pass filter.', ...
              ' ', ...
              '        tempo = conv(tempo,Hi_R);' ...
                );
wshowdrv('disp',figNumber,comment);

% Reconstruct detail : second step.
% Convolves upsampled detail coeff.
% and Hi_Rreconstruction high-pass filter.
        tempo = conv(tempo,Hi_R);

figure(active_fig)
subplot(313); plot(tempo);
xlabel('Upsampled detail coeff. and Hi_R convolution.');
set(h3,'Xlim',[1 length(tempo)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
              '% Reconstruct detail : last step.', ...
              '% Take central part of length 16', ...
              '% of upsampled detail coeff. and Hi_R', ...
              '% convolution.', ...
              ' ', ...
              '        d1 = wkeep(tempo,16);' ...
                );
wshowdrv('disp',figNumber,comment);

% Reconstruct detail : last step.
% Take central part of length 16
% of upsampled detail coeff. and Hi_R
% convolution.
        d1 = wkeep(tempo,16);


figure(active_fig)
subplot(313); plot(d1);
xlabel('Reconstructed detail d1.');
set(h3,'Xlim',[1 length(d1)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
            '% Using the same line reconstruct', ...
            '% approximation. First step. Upsample', ...
            '% ca1 inserting zeros.', ...
            ' ', ...
            '        tempo = dyadup(ca1);' ...
                );
wshowdrv('disp',figNumber,comment);

% Using the same line reconstruct
% approximation. First step. Upsample
% ca1 inserting zeros.
        tempo = dyadup(ca1);


figure(active_fig)
subplot(313); plot(tempo);
xlabel('Upsampled approx. coeff. ca1');
set(h3,'Xlim',[1 length(tempo)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
            '% Reconstruct approx. : second step.', ...
            '% Convolves upsampled approx coeff.', ...
            '% and Lo_R reconstruction low-pass filter.', ...
            ' ', ...
            '        tempo = conv(tempo,Lo_R);' ...
                );
wshowdrv('disp',figNumber,comment);

% Reconstruct approx. : second.
% Convolves upsampled approx coeff.
% and Lo_R reconstruction low-pass filter.
        tempo = conv(tempo,Lo_R);

figure(active_fig)
subplot(313); plot(tempo);
xlabel('Upsampled approx coeff. and Lo_R convolution.');
set(h3,'Xlim',[1 length(tempo)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
              '% Reconstruct approx : last step.', ...
              '% Take central part of length 16', ...
              '% of upsampled approx coeff. and Lo_R', ...
              '% convolution.', ...
              ' ', ...
              '        a1 = wkeep(tempo,16);' ...
              );
wshowdrv('disp',figNumber,comment);

% Reconstruct approx : last step
% Take central part of length 16
% of upsampled approx coeff. and Lo_R
% convolution.
        a1 = wkeep(tempo,16);

subplot(313); plot(a1);
xlabel('Reconstructed approx. a1.');
set(h3,'Xlim',[1 length(a1)]);

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
            '% Now add a1 and d1 in order to', ...
            '% recover original signal.', ...
            ' ', ...
            '        a0 = a1 + d1;' ...
            );
wshowdrv('disp',figNumber,comment);

% Now add a1 and d1 in order to
% recover original signal.
        a0 = a1 + d1;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

figure(active_fig)
h1 = subplot(311); plot(s,'r');
title('Original signal s.');
set(h1,'Xlim',[1 16]);
pause(1)
set(h1,'Xticklabel',[]);
figure(active_fig)
h2 = subplot(312); plot(1:16,a1,1:16,d1);
title('a1 and d1.');
set(h2,'Xlim',[1 16]);
pause(1)
set(h2,'Xticklabel',[]);
figure(active_fig)
h3 = subplot(313); plot(a0);
title('Reconstructed signal : a0 = a1 + d1.');
set(h3,'Xlim',[1 16]);
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
             '% This sequence is the one step', ...
             '% inverse discrete wavelet transform of s.', ...
             '% M-file idwt performs it directly.', ...
             ' ', ...
             '        a0 = idwt(ca1,cd1,''db1'',16);' ...
            );
wshowdrv('disp',figNumber,comment);

% This sequence is the one step
% inverse discrete wavelet transform
% of s. M-file idwt performs it
% directly.
        a0 = idwt(ca1,cd1,'db1',16);

figure(active_fig)
h3 = subplot(313); plot(a0,'g');
title('Reconstructed signal a0.');
set(h3,'Xlim',[1 16]);

if figNumber==0
        wshowdrv('pause',figNumber,wait_mode);
        delete(active_fig);
end
