% Mallat pyramid decomposition.

p = 3;
h = daub(2*p);
J = 5;
nlevels = 2;

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)')
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)
title(sprintf('Projection on to 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+1,1,1)
plot(x,PJ0f)
xlabel('x')
str = sprintf('P_%df(x)',J0);
ylabel(str)


i = 2;
total = PJ0f;
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+1,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;
  n = n*2;
  i = i+ 1;
end
subplot(nlevels+1,1,1)
title(strtot2)
pause

clf
subplot(211)
plot(x,total)
xlabel('x')
ylabel(strtot)

subplot(212)
plot(x,PJf)
xlabel('x')
ylabel(sprintf('P_%df(x)',J))

display(sprintf('Error is %e',norm(total-PJf)'))