% prediction error filter - laguardia airport
clear all;

% load laguardia airport temperature data
T = load('lgameanT.txt');
N = length(T);
dt = 1.0;
tmin = 0.0;
t = tmin + dt*[1:N]';
tmax = tmin + dt*(N-1);

% total power in signal
P = T'*T;

M=20;
G = zeros(N+1,M);
d = zeros(N+1,1);

% matrix formulation of temperature
for p = [1:M]
    G(p:N,p) = T(1:N-p+1);
    d(p)=0;
end

% single constrait with large epsilon forcing f(1) to be (-1)
G(N+1,1)=1e6;
d(N+1)=-1e6;

% generalized least squares inversion
f = inv(G'*G)*G'*d;

% compute predcion error filter * temperature data
% (cheat by doing convolution)
y = conv(f,T);

% plot prediction error filter
figure(1);
clf;
plot( t(1:M), f );
hold on;
xlabel( 'time, days' );
ylabel( 'filter' );

% plot temperature and error
figure(2);
clf;
Np = N;
plot( t(1:Np), T(1:Np), 'b' );
hold on;
plot( t(1:Np), y(1:Np), 'r' );
axis( [tmin, dt*Np, -50, 150]' );
xlabel( 'time, days' );
ylabel( 'temp, degF' );
E = y'*y;
title( sprintf('root mean square error %f percent', 100.0*sqrt(E/P) ));



