% prediction error filter - neuse river hydrograph

clear all;
% load Neuse River discharge data
D = load('neuse.txt');
T = D(:,2);
P = T'*T; % overall power in data
N = length(T);
dt = 1.0;
tmin = 0.0;
t = tmin + dt*[1:N]';
tmax = tmin + dt*(N-1);

% least squares matrices
M=100;
G = zeros(N+1,M);
d = zeros(N+1,1);
for p = [1:M]
    G(p:N,p) = T(1:N-p+1);
    d(p)=0;
end

% impose f(1) = (-1) with a single constraint with large epsilon
G(N+1,1)=1e6;
d(N+1)=-1e6;

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

% compute prediction error filter * discharge
% cheat - use convolution
y = conv(f,T);

% total prediction error
E = y'*y;

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



figure(2);
clf;
Np = 365;
plot( t(1:Np), T(1:Np), 'b' );
hold on;
plot( t(1:Np), y(1:Np), 'r' );
axis( [tmin, dt*Np, -5000, 25000]' );
xlabel( 'time, days' );
ylabel( 'discharge, cfs' );
title( sprintf('root mean square error %f percent', 100.0*sqrt(E/P) ));

figure(3);
clf;
Np = 365;
plot( t(1:Np), T(1:Np), 'b' );
hold on;
plot( t(1:Np), y(1:Np), 'r' );
axis( [tmin, dt*Np, -5000, 25000]' );
xlabel( 'time, days' );
ylabel( 'discharge, cfs' );
title( sprintf('root mean square error %f percent', 100.0*sqrt(E/P) ));



E