load Metdata2003rodan
d=tair;
t=day-day(1);
plot(t,d)

% smooth with a 1-d quadratic loess filter
help loess1d

dsm=loess1d(t,d,t,0.5);
plot(t,d,t,dsm);
han=plot(t,d,'r.',t,dsm);
set(han(2),'linew',2)
legend('data','loess 0.5 day timescale')
title('hourly data and quadratic loess (0.05 day) smoothed')
figure(gcf)
disp('Paused...');pause

% now decimate the data and use loess to fills the gaps
N = length(d);
rm = sort(round(rand([1 0.8*N])*N));
dups = find(diff(rm)==0);
rm(dups) = [];
rm(find(rm==0))=[];
disp(['removing ' int2str(100*length(rm)/length(d)) '% of data']);

d2 = d;
t2 = t;
d2(rm) = [];
t2(rm) = [];

han=plot(t,d,'.',t2,d2,'.');
set(han(1),'color',0.7*[1 1 1],'markersize',5)
set(han(2),'color','b','markersize',5)
legend('deleted data','remaining data')
title('decimated data')
figure(gcf)
disp('Paused...');pause

% rerun loess with 0.5 day half-width
dsm2 = loess1d(t2,d2,t,0.5,0);

han=plot(t,dsm,t,dsm2,t,d,'r.',t2,d2,'r.');
set(han([1 2]),'linew',2)
set(han(3),'color',0.7*[1 1 1],'markersize',5)
set(han(4),'color','r','markersize',5)
legend('loess all 0.5 day','loess some 0.5 day',...
   'deleted data','remaining data')
title('loess smoothed fits to full and decimated data sets')
figure(gcf)
disp('Paused...');pause


% but what time scales are appropriate to use for smoothing?

% look at autocorrelation
cf1 = xcorr(d);
cf1(1:length(d)-1)=[];
plot(t,cf1/cf1(1))
legend('all data')
xlabel('lag (days)')
title('normalized autocorelation of data')
figure(gcf)
disp('Paused...');pause

% seasonal scale shows
% remove it
fann = 2*pi/365.25;
E = [ones(size(t)) cos(fann*t) sin(fann*t)];
a = E\d;
dfit = a(1) + a(2)*cos(fann*t) + a(3)*sin(fann*t);

% or more easily:
dfit = E*a;

% look at autocorrelation
cf2 = xcorr(d-dfit);
cf2(1:length(d)-1)=[];
plot(t,cf1/cf1(1),t,cf2/cf2(1))
legend('all data','annual cycle removed')
set(gca,'xlim',[0 10])
xlabel('lag (days)')
title('autocorelation of data')
figure(gcf)
disp('Paused...');pause

% could also remove diurnal cycle
fann = 2*pi/365.25;
f24h = 2*pi/1;
E = [ones(size(t)) cos(fann*t) sin(fann*t) cos(f24h*t) sin(f24h*t)];
a = E\d;

dfit = E*a;

% look at autocorrelation
cf3 = xcorr(d-dfit);
cf3(1:length(d)-1)=[];
plot(t,cf1/cf1(1),t,cf2/cf2(1),t,cf3/cf3(1))
legend('all data','annual cycle removed','annual and diurnal removed')
set(gca,'xlim',[0 10])
xlabel('lag (days)')
title('autocorelation of data')
grid on
figure(gcf)
disp('done...');

% try removing all long time scale (10-day) variability with loess, and
% recompute autocorrelation



dsm3=loess1d(t,d,t,30,0);

E=[ones(size(t)) cos(2*pi*t) sin(2*pi*t)];
a=E\(d-dsm3);
dfit=E*a;

cf4=xcorr(d-dsm3-dfit);
cf4(1:length(d)-1)=[];
plot(t,cf1/cf1(1),t,cf2/cf2(1),t,cf3/cf3(1),t,cf4/cf4(1))
legend('all data','annual cycle removed','annual and diurnal removed',...
   'diurnal and 30-day smoothed removed')
set(gca,'xlim',[0 10])
xlabel('lag (days)')
title('autocorelation of data')
figure(gcf)