% 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)'))