%% analsynDWT.m
% DWT analysis and synthesis waveforms for a level 3 filter bank
% show analysis and synthesis waveforms for the DWT 
% show stepwise reconstruction
% uses wavelet toolbox

% S. Allen Broughton  6 Nov 2010

%% set up parameters, filters and transform matrix and inverse
N = 128; % this should be even
lev = 3;
dt = 0.1;

%% set up analysis filters
% uncomment the desired filter bank pair, not all work
% la = [1,1]/2; ha = [1,-1]/2;  % scaled Haar 
% la = [1,1]/sqrt(2); ha = [1,-1]/sqrt(2); % orthogonal Haar = db1
% la = [1+sqrt(3),3+sqrt(3),3-sqrt(3),1-sqrt(3)]/(4*sqrt(2)); ...
% ha =  [1-sqrt(3),-3+sqrt(3),3+sqrt(3),-1-sqrt(3)]/(4*sqrt(2)); % db2
[la ha] = wfilters('db3');
% la = [-1,2,6,2,-1 0]/8; ha = [0 0 -1,2,-1 0]/2; % Legall
% la = [-1,2,6,2,-1 0]/sqrt(1+4+36+4+1); ha = [0 0 -1,2,-1 0]/sqrt(1+4+1); % unit Legall
% la = [1,2,1 0]/4; ha = [1 -2 1 0]/4; % worthless
la
ha

%% compute transform matrix and inverse
T = dwtmatrix(N,lev,la,ha); 
Ti = inv(T); % works for non-orthogonal case

%% set up signal and transform  
t = ((0:(N-1))/N);
% X =  sin(2*pi*5.2*t)+0.3*rand(size(t));
% X =  [1:(N/2),(N/2):(-1):1]; X = X/max(abs(X))-0.5;% linear tent
% X =  t.*(1-t);  X = X/max(abs(X))-0.5; % quadratic tent
% X = rand(size(t));  X = X/max(abs(X))-0.5;  % random noise
 X = t.*(t-0.5).*(t-1); X =X/max(abs(X)); %cubic
X = X(:); % force a column vector
Xr = zeros(size(X)); % initialize a reconstruction vector
Y = T*X;

%% plot DFT of la and  ha
lax = [la, zeros(1,N-length(la))];
hax = [ha, zeros(1,N-length(ha))];
figure(1)
subplot(2,1,1)
FD = (-N/2):(N/2-1);
plot(FD,fftshift(abs(fft(lax))))
title('|DFT| of low pass filter');
subplot(2,1,2)
axis tight
plot(FD,fftshift(abs(fft(hax))))
title('|DFT| of high pass filter');
axis tight

%% graph window parameters
mx = max(T(:));
mn = min(T(:));
mxi = max(Ti(:));
mni = min(Ti(:));
mxX = max(X);
mnX = min(X);
mxY = max(Y);
mnY = min(Y);
eps = 0.1;

%%  iteration parameters and iterate
waitforkey = [1,2,3,4,5]

figure(2) 
for i = 1:N
   D = T;
   D(i,:)= mn*ones(1,N);
   R = T(i,:);
   subplot(2,3,1)
   imshow(D,[]); 
   title('Analysis matrix')
   subplot(2,3,2)
   plot(1:N,R,1:N,X);
   axis([1,N,min([mn,mnX]),max([mx,mxX])]);
   title('Analysis waveform and signal')
   subplot(2,3,3)
   plot(1:N,Y,'b.',i,Y(i),'r+');
   axis([1,N,mnY-eps,mxY+eps]);
   title('Transform and DWT coefficient')
   
   
   D=Ti;  
   D(:,i)= mni*ones(N,1);
   C = Ti(:,i);
   Xr = Xr+Y(i)*C;
   subplot(2,3,4)
   imshow(D,[]); 
   title('Synthesis matrix')
   subplot(2,3,5)
   plot(1:N,Ti(:,i));
   axis([1,N,mni,mxi]);
   title('Synthesis waveform')
   subplot(2,3,6)
   plot(1:N,Xr,'b',1:N,Y(i)*C,'r');
   axis([1,N,mnX-eps,mxX+eps]);
   title('Partially reconstructed signal and current term')
   drawnow;
   pauseforkey = sum(find(waitforkey==i)); 
   if pauseforkey > 0
      disp('Press a key to proceed'), pause,
      else pause(dt)
      end
   end;