% eda11_04
% fit of a sinusoid of unknown (but approximately annual)
% frequency to Black Rock Forest temperature data using the
% iterative least squares method
clear all;

% load Black Rock Forest temperature data
D=load('brf_temp.txt');
t=D(:,1);
dobs=D(:,2);
Ns=size(D);

% clean the data by deleting outliers and zero (= no data) values
k = find( (dobs>-35) & (dobs<40) & (dobs~=0.0) );
t = D(k,1);
dobs = D(k,2);
N = length(dobs);

% some statistics
dmean = mean(dobs);
amp = 0.5*(max(dobs)-min(dobs));

% Annual period, frequency, angular frequency
Ta = 365;
fa = 1/Ta;
wa = 2*pi*fa;

% nonlinear theory and its derivative
% note model parameters are scaled to be order unity
% d = m1*dmean + amp*m2*sin(m4 wa t) + amp*m3*cos(m4 wa t)
% dd/dm1 =  dmean
% dd/dm2 =  amp*sin(m4 wa t)
% dd/dm3 =  amp*cos(m4 wa t)
% dd/dm4 =  amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)

% inital guesses
M=4;
m = zeros(M,1);
m(1) =  1;          % average level is mean
m(2) =  0.01;       % phase with max in summer
m(3) = -0.99;       % so only -cos() term
m(4) =  1;         % annual frequency 

s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;

figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
plot(t,dpre,'-','LineWidth',2,'Color',[0.7,0.7,0.7]);
xlabel('time (days)');
ylabel('temperature (deg C)');

figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
xlabel('time (days)');
ylabel('temperature (deg C)');

s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;
dd = dobs - dpre;
E = dd'*dd;
for i=[1:10]
    Eprev = E;
    G = [ dmean*ones(N,1), s, c, m(2)*wa*t.*c-m(3)*wa*t.*s ];
    dm = (G'*G)\(G'*dd);
    m = m + dm;
    s = amp*sin(m(4)*wa*t);
    c = amp*cos(m(4)*wa*t);
    dpre = m(1)*dmean + m(2)*s + m(3)*c;
    dd = dobs - dpre;
    E = dd'*dd;
    if( (abs(E-Eprev)/E) < 1.0e-6 )
        break
    end
end

figure(1);
plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);
figure(2);
plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);

ffinal = (wa*m(4))/(2*pi);
Tfinal = (2*pi)/(wa*m(4));
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
axis([320,420,0,1]);
plot( [Tfinal, Tfinal]', [0, 0.75]', 'k-', 'LineWidth', 3 );
xlabel('period (days)');
ylabel('amplitude spectral density');

% covariance
% f = f0 + df   and   T = T0 + dT    and   T0 = 1 / f0
% T = 1 / (f0 + df ) = (1/f0) (1+df/f0)^-1
%   = (1/f0) (1-df/f0) = (1/f0) - df/(f0^2)
% dT = (-1/(f0^2)) df
% var(dT) = (-1/(f0^2))^2 var(df)

G = [ dmean*ones(N,1), s, c, m(2)*wa*t.*c-m(3)*wa*t.*s ];
sigma2d = E/(N-M);
GTGI = inv(G'*G);
sigma2m4 = sigma2d*GTGI(4,4);
sigma2f = (wa^2)*sigma2m4/((2*pi)^2);
sigma2T = sigma2f / (ffinal^4);
sigmaT = sqrt(sigma2T);
DT95 = 2*sigmaT;
plot( [Tfinal-DT95, Tfinal+DT95]', [0.7, 0.7]', 'k-', 'LineWidth', 3 );

fprintf('posterior sqrt(sigma2d) %f\n', sqrt(sigma2d));
fprintf('Estimated period %.3f +/- %.3f days\n', Tfinal, DT95 );

% compare with result of FFT

% standard time & frequency definitions
Dt = t(2)-t(1);
tstart = t(1);
tend = t(N);
N2 = floor((tend-tstart)/Dt)+1;
N2 = 2*floor(N2/2); % insure even
t2 = tstart + Dt*[0:N2-1]';
fny = 1/(2*Dt);
Nf = N2/2+1;
Df = fny / (N2/2);
f = Df*[0:Nf]';

% interpolate data to even sampling
% but note that data gaps filled in with stright lines
d2 = interp1(t,dobs,t2,'linear','extrap');

% demean, take FFT and compute amplitude spectral density
d2 = d2-mean(d2);
d2tilda = Dt*fft(d2);
d2tilda = d2tilda(1:Nf);
asd = abs(d2tilda);
maxasd = max(asd);

plot( 1./f(2:Nf), 0.9*asd(2:Nf)/maxasd, 'k-', 'LineWidth', 2 );
plot( 1./f(2:Nf), 0.9*asd(2:Nf)/maxasd, 'ko', 'LineWidth', 2 );




