% Compute discrete wavelet transform of the test signals created by routine
% signals_comp.m. We use Daubechies 2 scaling function and mother wavelet.
% 
% Here the decomposition is done only at one level, so no multiresolution.
%
% Boundaries are taken into account simply by zero extension, leading to 
% artefacts near the endpoints of the interval.
%
% Samuli Siltanen April 2004

% Convolution kernels for the Daubechies 2 wavelet
h = [1+sqrt(3) 3+sqrt(3) 3-sqrt(3) 1-sqrt(3)]/(4*sqrt(2));
g = [h(4)      -h(3)     h(2)      -h(1)];

% Loop over signals
for sss = 1:5
    
    % Load precomputed signal from file
    loadcommand = ['load s', num2str(sss), ' t s'];
    eval(loadcommand)
    
    % Compute smooth part of the signal and downsample
    sm = conv2(s,fliplr(h),'same');
    sm = sm(1:2:end);
    
    % Compute details of the signal and downsample
    de = conv2(s,fliplr(g),'same');
    de = de(1:2:end);
    
    
    
    
    % Reconstruct the signal by upsampling and convolution
    sm2 = [zeros(size(sm));sm];
    sm2 = sm2(:).';
    sm2 = conv2(sm2,h,'same');
    de2 = [zeros(size(de));de];
    de2 = de2(:).';
    de2 = conv2(de2,g,'same');
    recon = sm2 + de2;
    
    % Plot the result
    figure(sss)
    clf
    hold on
    plot(t,s,'b')
    plot(t,recon,'r')
    
end


