% This program is used to compute the multitaper coherence between two real
% signals (sig1 and sig2).
% format [squa,phase]=mtcoherence(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, Zhao M, Chao BF, Interannual polar motion with relation to 
%      the North Atlantic Oscillation. Global and Planet Change, 18: 79-84, 1998
%
% Percival D. B. and A. T. Walden, Spectral Analysis for Physical Applications: 
%     Multitaper and Conventional Univariate Techniques. Cambridge, England: Cambridge 
%     University Press£¬1993.
% .............................................................................
%
% Yonghong Zhou, July 2007, at Shanghai Astronomical Observatory 
%
function [squa,phase]=mtcoherence(sig1,sig2)
coh=mtcoh(sig1,sig2);
number=length(coh);
interval=input('sampling interval = ');
fnequist=1/(2*interval);
freq=fnequist*(0:(number-1))'/number;
squa=[freq,abs(coh).^2];
phase=[freq,angle(coh)*180/pi];
%
%-------------
figure(9),subplot(211),plot(squa(:,1),squa(:,2));
ylabel('Squared Coherence');
figure(9),subplot(212),plot(phase(:,1),phase(:,2));
ylabel('Coherence Phase(degree)');
xlabel('frequency');
% For example, if the unit of sampling interval is "year",
% then, the unit of frequency is "cycles per year". 
%--------------







% This program computes the 4pi-multitapered coherence spectrum 
% of real-valued time series datax and datay of length n.
% The coherence spectrum is complex-valued, so:
% abs(coh).^2 gives the squared-coherence, and angle(coh) gives the
% phase angle.
%
% format: function coh=mtcoh(datax,datay)
%
function coh=mtcoh(datax,datay)
n=length(datax);
m=7;
nw=(m+1)/2;   % m=2*nw-1
[v,lambda]=dpss(n,nw);
disp('fft coherence...')
sxx=zeros(n,1);
syy=sxx;
ssxy=sxx+i*sxx;
for j=1:m
   sx=fft(v(:,j).*datax);
   sy=fft(v(:,j).*datay);
   sxy=sx.*conj(sy);
   sx=sx.*conj(sx);
   sy=sy.*conj(sy);
   sxx=sxx+sx;
   syy=syy+sy;
   ssxy=ssxy+sxy;
end
n1=ceil(n/2);
sxx=sxx(1:n1);
syy=syy(1:n1);
ssxy=ssxy(1:n1);
coh=ssxy./sqrt((sxx.*syy));