% disegna la risposta all'impulso del sistema LTI descritto
% dall'equazione alle differenze definita dai vettori di
% coefficienti a e b
%
function lti_demo1(a,b)

	% coefficienti dell'equazione alle
	% differenze che descrive il sistema
	if nargin < 1
% 		a = [1 -0.8];
% 		a = [1 -1.6 0.64];
% 		a = [1 1.6 0.64];
        a = [1 -1.5 1];
	end
	if nargin < 2
		b = [1 0];
% 		b = [0 1];
% 		b = [0.7 0.3];
	end

	% ricaviamo M ed N dalle dimensioni dei vettori
	M = numel(a) - 1;
	N = numel(b) - 1;

	% off_x serve a gestire gli indici
	% "negativi" con matlab (che ha indici
	% a partire da 1)
	off_x = 1+N;

	% sequenza di ingresso delta di Kronecker
	Nx = 30;
	x = zeros(1,Nx+N);
	x(off_x + 0) = 1; % x(0) = 1

	% alloca la sequenza di uscita
	Ny = Nx;
	y = zeros(1,Ny+M);
	off_y = 1+M;
	
	% condizioni iniziali nulle
	for n = -M:-1
		y(off_y+ n) = 0;
	end

	% loop
	for n = 0:Nx-1
		
		% y(n) = (b0 x(n) + b1 x(n-1) + ...
		%   + bN x(n-N) - a1 y(n-1) 
		%   - a2 y(n-2) - ... aM y(n-M))/a0
		tmp = 0;
		for k = 0:N
			tmp = tmp + b(1+ k) * x(off_x+ n-k);
		end
		for k = 1:M
			tmp = tmp - a(1+ k) * y(off_y+ n-k);
		end
		y(off_y+ n) = tmp/a(1+ 0);
		
	end
	
	% plot
	subplot(2,1,1)
	stem(-N:Nx-1,x)
	set(gca,'XLim',[-1,Nx+1])
	set(gca,'YLim',[min(x)-1,max(x)+1])
	xlabel('n')
	ylabel('x(n)')
	
	subplot(2,1,2)
	stem(-M:Nx-1,y)
	set(gca,'XLim',[-1,Nx+1])
	set(gca,'YLim',[min(y)-1,max(y)+1])
	xlabel('n')
	ylabel('y(n)')
end
