% set up parameters 
N   = 128;
lev = 3;
dt = 0.1;

%set up analysis filters
%[la,ha] = wfilters('db4')
la = [-1,2,6,2,-1]; ha = [-1,2,-1]; 
%la = [1,1]/2; ha = [1 -1]/2; 


%compute transform matrix and inverse
T = fbanlcol(eye(N),la,ha,lev); 
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)*[1:(N/2),(N/2):(-1):1]-0.25; % linear tent
X = 127*(8*t.*(1-t)-1);  % quadratic tent
%X = rand(size(t));
X = X(:); % force a column vector
Xr = zeros(size(X)); % initialize a reconstruction vector
Y = T*X;
% eps = .05, Y = eps*round(Y/eps); % quantization

%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
waitforkey = [1,2,3]

figure(1) 
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;