% Example from 
% http://www.mathworks.com/access/helpdesk/help/toolbox/wavelet/ch04_ac2.html#134867
% Modified by DPO  02/2010 to facilitate the use of varying levels 
% and types of wavelets, to plot differently, and to compute 
% relative error and number of coefficients.

% Set the number of levels and the type of wavelets to use.
%
% Wavelet options include 'coif3' (Coiflets), 'db1' (Haar), 
%                         'db3' (Daubechies), 'bior2.4', 
%                     and 'dmey' (discrete Meyer).
% See 'doc wavedec' for more information.

nlev=4 
w = 'db1'

% Load original 1-D signal.

    load sumsin; s = sumsin;

% Perform the decomposition of s at level nlev. 

    [c,l] = wavedec(s,nlev,w);

% Reconstruct the approximation signals and detail signals at 
% levels 1 to nlev, using the wavelet decomposition structure [c,l].
for i = 1:nlev
    A(i,:) = wrcoef('a',c,l,w,i);
    D(i,:) = wrcoef('d',c,l,w,i);
end

% Plots. 
    t = 1:length(s);    % originally 100:900
    tmin = min(t);   
    tmax = max(t);   
    subplot(nlev+1,2,1); plot(t,s(t),'r'); 
    title(sprintf('Orig. signal and approx. 1 to %d.',nlev)); 
    axis([tmin tmax -3 3]) 
    subplot(nlev+1,2,2); 
    title(sprintf('Details 1 to %d.',nlev)); 
    for i = 1:nlev, 
        subplot(nlev+1,2,2*i+1); plot(t,A(i,t),'b'); 
        title(sprintf('A(%d,:) approximation',i)) 
        axis([tmin tmax -3 3]) 
        subplot(nlev+1,2,2*i+2); plot(t,D(i,t),'g'); 
        title(sprintf('D(%d,:) detail',i))     
    end


norms = norm(s);

for i=1:nlev,
  disp(sprintf('i = %d, Relative error = %f, Number of coefficients = %d',...
                     i, norm(s-A(i,:))/norms, l(nlev-i+2)))
end

disp(' ')
disp('Note that in the Matlab storage scheme, the most')
disp('detailed approximation (Level nlev) is in the 1st row of A')
disp('and the least detailed (Level 1) is in the last row.')
disp(' ')
disp('Also note that A(i,t) = A(i+1,t) + D(i+1,t).')