function [dof,Tscale,xc,lags] = degreesfreedom(D);
% [dof,Tscale,xc,lags] = degreesfreedom(D);
% Estimate degrees of freedom and autocorrelation scale for the ensemble 
% of time series in rows of data matrix D.
%
% The method assumes the rows of D are a set of time series of 
% length N = size(D,2) at M = size(D,1) different locations.
%
% John Wilkin, April 2004

% Compute the covariance of each time series
%   If rows of D are all zero there will be a warning from xcov
%   and the output will be NaN
warning off MATLAB:divideByZero
clear xc
for m=1:size(D,1)
   % 'coeff' gives normalized function xc(0)=1
   [xc(m,:),lags] = xcov(D(m,:),'coeff');
end
warning on MATLAB:divideByZero

% remove the rows with NaNs
rm = find(any(isnan(xc),2)==1);
xc(rm,:) = [];

% Take mean over all M time series to form the ensemble autocorelation
% function
xc = mean(xc,1);

% Retain half of cross-correlation sequence that starts from lag=0 (because
% it is symmetric)
xc = xc(find(lags==0):end);
lags = lags(find(lags==0):end);

% Integrate to first zero crossing, at t=t*, and find autocorrelation scale 
% Tscale such that: Tscale = integral_0^t* xc(t') dt' 
% (This assumes time series are recorded at uniform sampling intervals.)
% Tscale will be in the normalized units of the sampling interval.
tstar = min(find(xc<=0))-1;
Tscale = mean(xc(1:tstar))*tstar;

% Degrees of Freedom is estimated as the number of correlation time scales 
% in length of the data record
dof = round(size(D,2)/Tscale);