%		     The University of Queensland
%     School of Information Technology & Electrical Engineering
%	     COMS3100/7100 Introduction to Communications

% Generate graphs which demonstrate the Gibbs effect for the sinc
% function

% Generate the time-domain signal for which we are taking successive
% approximations of the DTFT

N = 31;
n0 = -N:N;
x = sinc(n0 / 2);

% Define the array of radian frequencies at which we will sample
% the successive approximations to the DTFT

omega = -pi:pi/1024:pi;

% The DTFT of the sinc function is a square wave in the frequency
% domain --- we define this now

X = zeros(1, 2049);
X(513:1536) = 2;

% Now for the main loop in which each of the graphs are drawn

i = 1;
while 1

  % Get the value of M for this graph and generate n, the array of time
  % indices

  M = 2*i - 1;
  n = -M:M;

  % Generate the subset of time-domain samples from which we will create
  % a partial DTFT

  xM = sinc(n / 2);

  % Generate the matrix whose elements are n omega, n increasing down
  % the columns, omega across the rows

  nomega = n' * omega;

  % To generate X_M, we need only multiply a row vector of the
  % the time-domain samples of x[n] by the matrix whose rows are
  % complex exponentials (the complex exponents of the elements
  % of nomega, defined above)

  XM = xM * exp(-j * nomega);

  % Now create two subplots.  In the first, plot the time-domain samples
  % we are using.

  subplot(2, 1, 1);

  % Firstly plot the `complete' time-domain signal

  stem(n0, x, 'k:');
  hold on;

  % Then plot the subset of samples we are using for the partial DTFT

  stem(n, xM, 'r', 'filled');
  hold off;

  % Set up labels and ensure uniform scaling on each graph

  xlabel('\fontname{times}\itn');
  ylabel('\fontname{times}\itx\rm[\itn\rm]');
  xlim([-N, N]);
  ylim([-0.25, 1.05]);

  % The second subplot contains the frequency domain representation

  subplot(2, 1, 2);

  % First plot the Mth order approximation to the DTFT, X_M(e^jw)

  plot(omega / pi, real(XM), 'b');
  hold on;

  % Then also plot the function X(e^jw) to which the approximations
  % converge

  plot(omega / pi, X, 'k:');
  hold off;

  % Set up labels and ensure uniform scaling on each graph

  xlabel('\fontname{times}\it\pi\omega');
  ylabel(['\fontname{times}\itX\rm_{', num2str(M), '}(\ite^{j \omega}\rm)']);
  xlim([-1, 1]);
  ylim([-0.5, 2.5]);

  pause;

  i = i + 1;

end;