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