% Mallat pyramid decomposition.

p = 1;
p = input('Order of wavelet, p (defaults to 1) = ');
if isempty(p)
  p = 1;
end 
h = daub(2*p);
%J = 7;
J = input('Finest resolution, J (defaults to 7) = ');
if isempty(J)
  J = 7;
end 
%nlevels = 4;
nlevels = input('Number of levels, nlevels (defaults to 4) = ');
if isempty(nlevels)
  nlevels = 4;
end 

clf
J1 = J+2;;
nx = 2^J1;
dx = 1/nx;
x = (0:nx-dx/2)'/nx;
f = exp(-50*(x-0.5+dx/2).^2);
%f = ones(nx,1);
%f = x;
%f = sin(3*pi*x);
plot(x,f);
title('Original function, f(x)')
axis([0 1 0 1.1])
xlabel('x')
ylabel('f(x)')
pause

% Projection on to V_J.
cJ = scalecoeffs(f,1,h,J,J1);
PJf = expand(cJ,1,h,J,J1);
plot(x,PJf)
axis([0 1 0 1.1])
title(sprintf('Projection on to finest scale, V_%d',J))
xlabel('x')
ylabel(sprintf('P_%df(x)',J))
pause

% Create the Mallat pyramid.
cJdec = fwt(cJ,nlevels,h,0);

% Extract the wavelet coefficients and compute the
% projections on to V_J0, W_J0, W_J0+1, ..., W_J
J0 = J - nlevels;
n = 2^J0;
cJ0 = cJdec(1:n);
PJ0f = expand(cJ0,1,h,J0,J1);
subplot(nlevels+2,1,1)
plot(x,PJ0f)
xlabel('x')
str = sprintf('P_%df(x)',J0);
ylabel(str)


i = 2;
total = PJ0f;
subplot(nlevels+2,1,nlevels+2);
plot(x,total)
ylabel('Total')
xlabel('x')
pause
strtot = str;
strtot2 = sprintf('Projections on to V_%d',J0);
for j = J0:J-1
  dj = cJdec(n+1:2*n);
  Qjf = wexpand(dj,1,h,j,J1);
  subplot(nlevels+2,1,i)
  plot(x,Qjf)
  xlabel('x')
  str = sprintf('Q_%df(x)',j);
  ylabel(str)
  strtot = strcat(strtot, ' + ', str);
  strtot2 = strcat(strtot2, sprintf(', W_%d',j));
  total = total + Qjf;
  subplot(nlevels+2,1,nlevels+2);
  plot(x,total)
  ylabel('Total')
  xlabel('x')
  if (j < J-1)
    pause
  end
  n = n*2;
  i = i+ 1;
end
subplot(nlevels+2,1,1)
title(strtot2)
pause

clf
subplot(211)
plot(x,total,'-',x,f,':')
xlabel('x')
ylabel(strtot)
title('Multiple scales combined are equivalent to single scale projection')
legend('Projection','Original function')

subplot(212)
plot(x,PJf,'-',x,f,':')
xlabel('x')
ylabel(sprintf('P_%df(x)',J))
title('Single scale projection')
legend('Projection','Original function')

display(sprintf('Approximation error in either case is %e',norm(total-f)'))