% set up parameters 
N   = 32;

%compute transform matrix and inverse
T   = dct(eye(N));
Ti = inv(T);

% set up signal and transform
t = ((0:(N-1))/N);
%X =  sin(2*pi*5.2*t);
X =  cos(pi*3.1*t);
%X = (1/N)*[1:(N/2),(N/2):(-1):1]-0.25; % linear tent
%X = rand(size(t));
X = X(:); % force a column vector
Xr = zeros(size(X)); % initialize a reconsruction vector
Y = T*X;

% 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,4,5]

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 DCT 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-eps,mxi+eps]);
   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(1)
      end
   end;