function R=f_autocorr(SDATA, SF)
 
%** function to compute autocorrelation following Davis, 1986
%
%   SDATA   (detrended) time series
%   SF      sample frequency
%   tau     lag, use 1/4 time series length
%   R       correlation 
% 

 
N=length(SDATA);
tau=floor(N/4);
 
R=zeros(tau,2);
R(:,1)=[1:tau]*SF; 

for t=1:tau
   R(t,2)=sum(SDATA(1+t:N).*SDATA(1:N-t))-sum(SDATA(1+t:N))*sum(SDATA(1:N-t))/ ...
       sqrt( (sum(SDATA(1+t:N).^2)- sum(SDATA(1+t:N))^2) * (sum(SDATA(1:N-t).^2)- sum(SDATA(1:N-t))^2) );
end