Linear predictive extrapolation

extrapolation_demo.m

% This code answers the question at http://dsp.stackexchange.com/a/110/64

N = 150;    % Order of LPC auto-regressive model
P = 500;    % Number of samples in the extrapolated time series
M = 150;    % Point at which to start predicting

t = 1:P;

x = 5*sin(t/3.7+.3)+3*sin(t/1.3+.1)+2*sin(t/34.7+.7); %This is the measured signal

a = arburg(x, N);

y = zeros(1, P);

% fill in the known part of the time series
y(1:M) = x(1:M);

% Run the initial timeseries through the filter to get the filter state 
[~, zf] = filter(-[0 a(2:end)], 1, x(1:M));     

% Now use the filter as an IIR to extrapolate
y((M+1):P) = filter([0 0], -a, zeros(1, P-M), zf); 

% The above filter command can be replaced with this for-loop:
%
% for ii=(M+1):P      
%   y(ii) = -sum(a(2:end) .* y((ii-1):-1:(ii-N)));
% end

plot(t, x,      ...    % signal
     t, y,      ...    % extrapolation
     t, x - y);        % difference
   
l = line(M*[1 1], get(gca, 'ylim'));
set(l, 'color', [0,0,0]);
legend('actual signal', 'extrapolated signal', 'error','start of extrapolation');