function [yds,ymean,annualsignal] = deseason(t,y)
% [yds,ymean,annualsignal] = 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
% annualsignal 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 annualsignal
% value of the 366th day of the year is calculated as a simple average of
% the annualsignal 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 annualsignal 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.
% Input check:
if length(t)~=length(y)
disp('Error: vectors t and y must be same length.')
return
end
t = floor(t); % ensures we're using simple daily values
tdoy = date2doy(t);
annualsignal = NaN(366,1);
yds = NaN(size(y)); % this'll be the deseasoned signal.
ymean = NaN(size(y)); % This'll be a time series of seasonal signal the length of the full record.
for d = 1:365
annualsignal(d) = mean(y(tdoy==d));
yds(tdoy==d) = y(tdoy==d)-annualsignal(d);
ymean(tdoy==d) = annualsignal(d);
end
% Linearly interpolate a mean value for Dec 31 of leap years:
annualsignal(366) = (annualsignal(365)+annualsignal(1))/2;
ymean(tdoy==366) = annualsignal(366);
yds(tdoy==366) = y(tdoy==366)-annualsignal(366);
% Restore the DC offset:
yds = yds + mean(y);
end