% This program is used to make the cross-correlation analysis, i.e., 
% to compute the cross-correlation between two signals (sig1 and sig2). 
% format corr=correlation(sig1,sig2)
% Note: Each signal (sig1 or sig2) should be given as a column vector.
%
% .......................................................................
% References: 
%
% Zhou YH, SalsteinDA, Chen JL, Revised atmospheric excitation function series 
%     related to Earth's variable rotation under consideration of surface topography, 
%     Journal of Geophysical Research, 111, D12108, doi:10.1029/2005JD006608, 2006.
%
% Zhou YH, Zheng DW, Liao XH, Wavelet analysis of interannual LOD, AAM 
%     and ENSO: 1997~98 El Nino and 1998~99 La Nino signals, Journal of Geodesy, 
%     Vol.75, Page 164-168, 2001.
% .......................................................................
%
% Yonghong Zhou, July 2007, at Shanghai Astronomical Observatory 
%
function corr=correlation(sig1,sig2)
sig1=detrend(sig1,0);
sig2=detrend(sig2,0);
corr=xcorr(sig1,sig2,'coeff');
maximum=max(corr);
fprintf('\nThe maximum cross-coeffient =                            %g',maximum);
[mm,nn]=size(corr);
corr=[(-(mm-1)/2:(mm-1)/2)',corr];
corr1=corr(:,1);
corr2=corr(:,2);
%xx=find(corr2>maximum-10^(-10));
xx=find(corr2==maximum);
if length(xx)>1 
   fprintf('\n please check!\n');
end
lag=corr1(xx);
%The POSITIVE lag means that sig1 leads sig2, and vice versa.
fprintf('\nThe LAG value that corresponds to the maximum cross-coeffient = %g',lag);
fprintf('\nNow, plot cross-coeffient')
h=input('input figure handle ----');
figure(h),plot(corr(:,1),corr(:,2));
xlabel('Time Lag');
ylabel('Cross-Correlation');