function linconv
% LINCONV Demo on linear convolution sum
%
%   Type h[n] and x[n] in the prompt, see their
%   convolution sum in the first figure.
%   In the second figure there is the convolution
%   step by step.
%
%   y[n] = SUM{k=-Inf..Inf}(h[k] x[n-k])
%
%   Maximum length of h or x is 30.
%
%   Example: h[n] = [1 1 1]
%            x[n} = [-2 3 3 1 2 4 5]
%   >> linconv
%    ... h[n]: [1 1 1]
%    ... x[n]: [-2 3 3 1 2 4 5]
%
% See also CONV, FLIPLR

% Jukka Parviainen, 18.10.1999, 2.2.2000

h = input('Type the impulse response h[n] : ');
h = check_input(h);
x = input('Type the input sequence x[n] : ');
x = check_input(x);

hlen = length(h);
xlen = length(x);
ylen = hlen+xlen-1;
y = conv(x,h);     % normal convolution is calculated in this way

% these vectors are for plotting
xv = [zeros(1,xlen+2) x zeros(1,ylen-xlen+2)];
hv = [zeros(1,xlen+2) h zeros(1,ylen-hlen+2)];
yv = [zeros(1,xlen+2) y 0 0];
nv = [-xlen-2:ylen+1];

% In the first figure h, x and y.
figure;
subplot(311),stem(nv,hv), title('Impulse response h[n]'); grid;
subplot(312),stem(nv,xv), title('Input sequence x[n]'); grid;
subplot(313),stem(nv,yv), title('Output y[n]=h[n]*x[n]'); grid;

fprintf(1, ['The convolution of input sequence x and impulse ' ...
        'response h\n']);
fprintf(1,'Press any key to see the convolution sum step by step...\n');
pause;

% In the second figure animation of constructing
% a convolution sum is shown

figure;
subplot(411), stem(nv,hv), title('Impulse response h[n]'); grid;

% x[n-k], x is inversed with fliplr
xf = fliplr(x);
ysv = zeros(1,2*xlen+hlen+3);

% animation steps. One extra step in the beginning
% and one extra in the end.
for i=-1:ylen
  subplot(412), 
  xv = [zeros(1,i+3) xf zeros(1,ylen+1-i)];  
  ynv = hv.*xv;
  yns = sum(ynv);
  ysv = ysv + [zeros(1,xlen+i+2) yns zeros(1, ylen+1-i)];
  
  subplot(412), stem(nv,xv); 
  title([ 'x[' num2str(i) '-k]' ]); grid;

  subplot(413), stem(nv,ynv); 
  axis([-xlen-3 ylen+2 min([min(y) 0]) max([max(y) 0])]); % ONKO JÄRKEVÄ
  title(['y[' num2str(i) '] = SUM(h[k] . x[' ...
         num2str(i) '-k]) = ' num2str(yns) ]); 
  grid;

  subplot(414), stem(nv,ysv); 
  axis([-xlen-3 ylen+2 min([min(y) 0]) max([max(y) 0])]); 
  title(['y[n], n = -Inf .. ' num2str(i)]); grid;

% A little red mark for the position mark
%  ydif = (max(ysv)-min(ysv))/20;
  ydif = 0;
%  if i>=0
    hl = line([i-.2 i+.2 i+.2 i-.2 i-.2], ...
              [ydif ydif -ydif -ydif ydif]);
    set(hl,'Color',[1 0 0]);
    set(hl,'Linewidth',4);
%  end

  fprintf(1,['A dot product of sequences h[k] and ' ...
             'x[' num2str(i) '-k] = ' num2str(yns) '.\n']);
  if i~=ylen
    fprintf(1,['Press any key to go to n=' num2str(i+1) '.\n\n']);
    pause;
  else
    subplot(414); title('y[n]');
  end;
end;

return;

% Subfunction for checking the arguments

function x = check_input(x)

len = length(size(x));
if len>2
  error('Too many dimensions in the argument');
end
mind = min(size(x));
if mind>1
  error('Argument must be a vector not a matrix');
end
mind = max(size(x));
if mind>30
  error('Maximum length of sequence is 30');
end
if length(x)>1
  if size(x,1) > 1
    x = x';
  end
end