% set up parameters 
N   = 64;

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

% set up signal and transform
t = ((0:(N-1))/N);
X =  cos(2*pi*1.5*t)- 0.5*cos(2*pi*10*t)+  0.5*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;


figure(1) 
for i = 1:N
   R = T(i,:);
   subplot(2,2,1)
   plot(1:N,R,'b',1:N,X,'r');
   axis([1,N,min([mn,mnX]),max([mx,mxX])]);
   title('Analysis waveform and signal')
   subplot(2,2,2)
   plot(1:N,Y,'b.',i,Y(i),'r+');
   axis([1,N,mnY-eps,mxY+eps]);
   title('Transform and DCT coefficient')
   
   
   C = Ti(:,i);
   Xr = Xr+Y(i)*C;
   subplot(2,2,3)
   plot(1:N,Ti(:,i));
   axis([1,N,mni-eps,mxi+eps]);
   title('Synthesis waveform')
   subplot(2,2,4)
   plot(1:N,Xr,'r',1:N,Y(i)*C,'b',1:N,X-Xr,'g');
   axis([1,N,mnX-eps,mxX+eps]);
   title('Partially reconstructed signal, current term and difference')
  
   drawnow;
   pause
   end;