%  This script graphically shows the convolution operation between
%  two discrete sequences.  x is assumed the signal and h the system
%  impulse response to simulate the output of a linear system
%
%  Written by Kevin D. Donohue (donohue@engr.uky.edu) January 2009

x = [1 2 2 2 -2 -2 -2 1];  %  Input signal
h = [.5 1 0 1 .5];  % System response

%  Set range of convolution axes for full simulation
%  Assume index n=0 is is the first point in each vector 
tstart = -length(h);  %  Time axis begining
tend = length(x)+length(h);  %  Time axis end
cax = tstart:tend;  %  Convolution axis

% Embedd x and h in vector sychronized to convolution axis
xcax = zeros(size(cax));  %  Initialize signal vector on convolution axis with zeros
st = find(cax == 0);  %  Find index of 0 on axis
xcax(st:st+length(x)-1) = x;  %  embedd signal x on convolution axis
hcax = zeros(size(cax));   %  Initalize impulse response on convolution axis with zeros
hcax(1:length(h)) = fliplr(h);  %  time reverse h 
hshow = zeros(size(cax));   % 
hshow(st:st+length(h)-1) = h;  %  Not time reversed to plot h  
yout = zeros(size(cax));  % Initialize output vector
%  Plot signal 
figure(1)
stem(cax,xcax,'rx','Linewidth', 2)
title(['x(n)'])
xlabel('n')
set(gcf,'Position', [232   457   586   209])
%  Plot system response
figure(2)
stem(cax,hshow,'ko','Linewidth', 2)
title(['h(n)'])
xlabel('n')
set(gcf,'Position', [232   457   586   209])

%  Create figure window for display and fix axes to
%  cover the full simulation
figure(3)
g2 = subplot(2,1,2);  %  Dynamic system and input window
g1 = subplot(2,1,1);  %  Output window
set(gcf, 'Position', [123   214   720   381])

%  Loop throught each output sample showing 
%  the products between samples to be mulitplied and summed
for n=-1:length(x)+length(h)
    %  Plot function in convolution sum sychronized with convolution axis
  figure(3)
  subplot(2,1,2); stem(cax-.05,xcax,'rx', 'Linewidth', 2);
  hold on; stem(cax+.05,hcax,'ko', 'Linewidth', 2); hold off

  xlabel('m')
  title(['System response h(' int2str(n) '-m) in black X, Signal x(m) in red O'])
  %  Find index of current output point
  convind = find(cax == n);
  yout(convind) = sum(xcax.*hcax); %  Compute product and sum
  subplot(2,1,1); stem(cax(1:convind)-.05,yout(1:convind),'g>', 'Linewidth', 2);
  axis([tstart, tend, -4, 4])  %
  yl = get(g1,'Ylim');  %  get yaxis limits on plot to place text
  text(cax(1)+.5,yl(2)-.2*(yl(2)-yl(1)), ...
      ['$$ y(n)=  \sum_{m=-1}^{' int2str(length(x)+length(h)-1) '} x(m)h(n-m) $$ ' ],'Interpreter', 'Latex')
  title(['Convolution output up to n = ' int2str(n)])
  xlabel('n')
  pause  %  Wait for ket to be hit
  % slide h vector one sample to the right to denote shift from n increment
  % insert zeros in from the left
  hcax = [0 hcax(1:end-1)];
end