% This demonstration takes a pulse p_1(t) and its Fourier transform
% sinc(w/(2*pi)) and attempts to show how Fourier reconstruction happens

fh = figure;

% Original time domain signal
tv = -5:0.01:5;
xv = zeros(size(tv));  xv(find(abs(tv)<1/2)) = 1;
subplot(4,1,1);
plot(tv,xv);
axis([-5 5 -0.5 1.5]);
xlabel('t');  ylabel('x(t)');

% Frequency domain representation
wv = -30:0.01:30;
for i=1:length(wv)
  w = wv(i);
  if w~=0
    Xv(i) = sin(w/2)./(w/2);  
  else
    Xv(i) = 1;
  end
end
subplot(4,1,2);
plot(wv,Xv);
axis([-30 30 -0.5 1.5]);
xlabel('\omega');  ylabel('X(\omega)');

ws = 1;
xrv = zeros(size(xv));
for n = 0:floor(30/ws)
  w = n*ws;
  lw = w - ws/2;  uw = w + ws/2;  
  if w~=0
    Xm = sin(w/2)./(w/2);
  else
    Xm = 1;
  end
  ph1 = patch([lw lw uw uw], [-0.5 1.5 1.5 -0.5], 'r');
  ph2 = patch([-lw -lw -uw -uw], [-0.5 1.5 1.5 -0.5], 'r');
  
  % Reconstructed sample for frequency band
  xbv = ws/pi*Xm*cos(n*ws*tv);
  subplot(4,1,3);
  plot(tv,xbv);
  axis([-5 5 -1 1]);
  
  % Accumulated
  xrv = xrv + xbv;
  subplot(4,1,4);
  plot(tv,xrv);
  axis([-5 5 -0.5 1.5]);
  
  pause;
  subplot(4,1,2);
  delete(ph1);  delete(ph2);
end