function [yds,ymean,annsig] = deseason(t,y)
% [yds,ymean,annsig] = deseason(t,y) returns a time series yds of the
% same length as the multi-year daily time series y.  ymean is a vector the
% same length as t and y, containing only the seasonal signal.  The vector
% annsig is of length 366, the annual signal of daily values for day
% of year 1 to 366. Note that this will only work for daily records with t
% in datenum format.  
% 
% This function works by finding the mean value of the signal corresponding to 
% for each day of the year. That mean value is then subtracted from the raw signal 
% corresponding to that day of every year. For leap years, the annsig
% value of the 366th day of the year is calculated as a simple average of
% the annsig for Dec. 31 and Jan. 1. 
% 
% Please note that use of this function may be considered somewhat sloppy
% science. It's a clunky method of removing seasonal signals, but if your
% time series is long enough this deseason function is a quick and easy way
% to compare the magnitude of seasonal variations versus other variations
% in your signal. Or simply plotting annsig is a fair qualitative
% measure of the overall shape of a seasonal signal.  
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
% REQUIREMENTS: 
% 
% This function requires Anthony Kendall's date2doy function available on the 
% Mathworks FEX site: http://www.mathworks.com/matlabcentral/fileexchange/18316
% 
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
% EXAMPLE: Deseason some daily time series of temperature or something: 
% 
% % Declare some input data with seasonal signal in it: 
% t = (datenum('Jan 5, 2000'):datenum('May 15, 2016'))'; 
% temp = 25*cos(2*pi*date2doy(t)/365.25) + 30*rand(size(t))-15; 
% 
% % Deseason the data:
% [yds,ymean] = deseason(t,temp);
% 
% % Plot the raw and deseasoned data:
% figure
% set(gcf,'position',[100 100 600 250]); 
% plot(t,temp)
% datetick
% ylabel('temperature or something');
% axis tight
% box off
% hold on;
% plot(t,yds,'k')
% plot(t,ymean,'r')
% legend('raw daily data','deseasoned data','seasonal signal',...
%     'location','eastoutside'); 
% legend boxoff
% 
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
% Created by Chad A. Greene, August 2013. 
% version 2: now restores DC offset to yds, January 2014. 
%
% Nathan Steiger, LDEO, Columbia University, March 2019
% Updated code to smooth the annual cycle and make it faster

% Input check:
if length(t)~=length(y)
    disp('Error: vectors t and y must be same length.')
    return
end
y=y(:);
t = floor(t); % ensures we're using simple daily values 
tdoy = date2doy(t); 
annsig = NaN(366,1); 
for d = 1:365
    annsig(d) = mean(y(tdoy==d)); 
end
% Apply smooth annual cycle (it's likely not smooth because of sampling)
annsig=smooth(annsig,14,'lowess');
% Estimate leap day because of few samples
annsig(366) = (annsig(365)+annsig(1))/2;

% Mean annual cycle repeated
ymean=annsig(tdoy);
%yds=y-ymean+mean(y);
yds=y-ymean;

end
