% Problem 9: Wiener predictor.
% Rich Kozick, Spring 1997

% Signal to predict
y = [6.47 -1.28 3.75 1.96 -.77 4.15 -4.05 1.66 -2.06 -3.33 .09 -6.07 -.41 -3.62 -1.53 .48 -1.94 2.36 .15 2.22 2.26 1.07 2.27 .69 .76 -.50 -.64 -1.99 -1.60 -3.23 -2.22 -2.42]';

Ly = length(y);

% Signal delayed by D
D = 1;
x = [zeros(D,1); y(1:Ly-D)];

N = 1;		% Length of FIR Wiener prediction filter

% Compute Wiener filter coefficients

P = zeros(N,1);
R = zeros(N,N);

for k=N:Ly
  xk = x(k:-1:(k-N+1));
  P = P + y(k)*xk;
  R = R + xk*xk';
end

w = R \ P;

yhat = filter(w, 1, x);

ii = 1:Ly;
plot(ii, y(ii), '-', ii, y(ii), 'o', ii, yhat(ii), 'x')
legend('Original', 'Original', 'Prediction')