% Examples of LSQ for empirical model Fourier Series

addpath ../MiscUtil/
clearf

dt=15; % 15 days
T=360*24; % 24 years
t=[1:dt:T]';
w=2*pi/T;

% intraseasonal component
n1=T/(30*3); % 3 months
var_a1=1; a1=sqrt(2*var_a1);

% seasonal component
n2=T/(360);  % 12 months
var_a2=2; a2=sqrt(2*var_a2);

% interannual component
n3=T/(360*2); % 2 years
var_a3=0.8; a3=sqrt(2*var_a3);

% longer-term variations
n4=T/(360*5); % 5 years
var_a4=0.5; a4=sqrt(2*var_a4);

% very long unresolved frequency
n5=T/(360*50); % 50 years
var_a5=0.9; a5=sqrt(2*var_a5);

% very short unresolved frequency
n6=T/(10); % 10 days
var_a6=0.9; a6=sqrt(2*var_a6);

% CASE 1 : no noise, only wave components.
y= a1*sin(n1*w*t) + a2*cos(n2*w*t) + a3*sin(n3*w*t) + a4*cos(n4*w*t) ;

% CASE 2 : add noise
n=randn(length(t),1); n=n/var(n)*10;
y= a1*sin(n1*w*t) + a2*cos(n2*w*t) + a3*sin(n3*w*t) + a4*cos(n4*w*t) + n;

% CASE 3 : unresolved frequencies
n=randn(length(t),1); n=n/var(n)*0;
y= 0*[a1*sin(n1*w*t) + a2*cos(n2*w*t) + a3*sin(n3*w*t) + a4*cos(n4*w*t)] + ...
   [a5*sin(n5*w*t) + a6*cos(n6*w*t)] + n;

% CASE 4 : El Nino
load nino
dt=30;
t=[1:702]*30;
y=n.nino34 - mean(n.nino34);
load PacClimateIND.mat
y=ind.enso.ba.cpi;
%y=ind.pdv.pdo;
y=detrend(y);
y=nn(y-mean(y));
t=[1:length(y)]*30;
icase=4;



% Fit empirical model
[an, bn, period, yhat, yhat_filter, E, amp] = utl_sincos(y,dt);

% Filter the data
%bandperiod(1) = 360*0;
%bandperiod(2) = 360*1;
%yhat= utl_sincos_filter(y,dt, E, amp, bandperiod);



clf; 
subplot(3,1,1); 
plot(t/360, y); xlabel('years'); ylabel('SST'); hold on
plot(t/360, yhat,'r'); xlabel('years'); ylabel('SST'); 
legend 'y' 'yhat'

subplot(3,1,2)
plot( [an.^2 + bn.^2]/2); xlabel('coefficient number'); ylabel('variance');

subplot(3,1,3)
plot(period/360, [an.^2 + bn.^2]/2); xlabel('period years'); ylabel('variance');
hold on; plot(period/360, [an.^2 + bn.^2]/2,'*');
set(gca, 'xlim', [0 8]);
if icase==4, set(gca, 'xlim', [0 20]); end


return
mfig;
[sigh2, freq , period, sigh2_AR1, sigh2_AR1_90, phase ] = myspec(y,dt,1);
plot_myspec (sigh2, sigh2_AR1, sigh2_AR1_90, freq)