function y = predict( x, p, n, flag )
% y = predict( x, p, n, flag)
% -----------------------
% Produce a series y whose elements are the predictions for
%  x, starting at sample n+1, based on a p-order model
%  trained on the first n samples of x.
%
%  x    - input time series, of length N
%  p    - order of predictor to be trained
%  n    - number of x samples to be used for training predictor
%  flag - if nonzero, predictions will be used in predicting
%         subsequent samples of y; otherwise, y will always be
%         estimated based on true samples (of x) alone
%
% (c) 2004 kornel laskowski cmu 15-782
%

% use lpc routine available in signal processing toolbox
% returns A=[ 1 A(2) ... A(N+1) ], where A are to be used
% in predictions as follows:
% Xp(n) = -A(2)*X(n-1) - A(3)*X(n-2) - ... - A(N+1)*X(n-N)

A = real(lpc((x(1:n)-mean(x(1:n))),p));

% set up initial values
Zi = zeros(1,p);
for i=1:p
	Zi(i) = x(n-p+i);
	for j=2:i
		Zi(i) = Zi(i) + A(j)*x(n-p+1+i-j);
	end
end

b = [ 0 -A(2:end)];

if flag == 0

	% since each prediction is based only on the true values of
	% x (ie. no previous predicted values), y is just the output
	% of an FIR filter applied to x.
	%
	% a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
	%                       - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

	b = [ 0 -A(2:end)];
	y = filter(b,1,x(n-p+1:end),Zi);

else

	% each prediction is based on past prediction values, except
	% for the initial conditions. so y is the output of an IIR
	% filter applied to a zero input. the initial energy in the
	% system is due to the initial conditions alone.
	%
	% a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
	%                       - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

	a = A;
	y = filter(0,a,zeros(1,length(x)-n+p),Zi);

end

y = y(p+1:end);