%function  [s, c, ph, ci, phi] = cmtm(x,y,dt,NW,qbias,confn,qplot);
%
%Multi-taper method coherence using adaptive weighting and correcting 
%for the bias inherent to coherence estimates.  The 95% coherence
%confidence level is computed by cohconf.m.  In addition, a built-in
%Monte Carlo estimation procedure is available to estimate phase 95%
%confidence limits. 
%
% Inputs:
%         x     - Input data vector 1.  
%         y     - Input data vector 2.  
%         dt    - Sampling interval (default 1) 
%         NW    - Number of windows to use (default 8) 
%         qbias - Correct coherence estimate for bias (yes, 1)  (no, 0, default).
%         confn - Number of iterations to use in estimating phase uncertainty using a Monte Carlo
%                 method. (default 0)
%         qplot - Plot the results, (yes, 1), (No, 0, default).  The upper tickmarks indicate the
%                 bandwidth of the coherence and phase estimates.  
%
% Outputs:
%         s       - frequency
%         c       - coherence
%         ph      - phase
%         ci      - 95% coherence confidence level
%         phi     - 95% phase confidence interval, bias corrected
%                   (add and subtract phi from ph).
%
%
%required files: cohconf.m, cohbias.m, cohbias.mat, Matlab signal processing toolbox.
%
%Peter Huybers
%MIT, 2003
%phuyber@mit.edu

function [s, c, ph, ci, phi] = cmtm(x,y,dt,NW,qbias,confn,qplot);

%check input
if nargin<2, help cmtm; return; end;
if nargin<7, qplot=0;  end;
if nargin<6, confn=0;  end;
if nargin<5, qbias=0;  end;
if nargin<4, NW=8;     end;
if length(NW)==0, NW=8;end;
if nargin<3, dt=1;     end;
if length(dt)==0, dt=1;end;


if NW<1.5, disp('Warning: NW must be greater or equal to 1.5'); return; end;
if nargin>4, 
	  disp('-------------------------');
disp(['Number of windows: ',num2str(NW)]);
if qbias==1, disp('Bias correction:   On');
else,        disp('Bias correction:   Off'); end;
disp(['Confidence Itera.: ',num2str(confn)]);
if qplot==1, disp('Plotting:          On');
else,        disp('Plotting:          Off'); end;
disp('-------------------------');
end;

x=x(:)-mean(x); 
y=y(:)-mean(y);
if length(x)~=length(y), disp('Warning: the lengths of x and y must be equal.'); return; end;

%define some parameters
N   = length(x);
k   = min(round(2*NW),N); 
k   = max(k-1,1);
s   = (0:1/(N*dt):1/dt-1/(N*dt))';
pls=2:(N+1)/2+1;
v   = (2*NW-1); %approximate degrees of freedom

if rem(length(y),2)==1; pls=pls(1:end-1); end;

%Compute the discrete prolate spheroidal sequences, requires the spectral analysis toolbox.
[E,V]=dpss(N,NW,k);

%Compute the windowed DFTs.
fkx=fft(E(:,1:k).*x(:,ones(1,k)),N);
fky=fft(E(:,1:k).*y(:,ones(1,k)),N);

Pkx=abs(fkx).^2;
Pky=abs(fky).^2;

%Itteration to determine adaptive weights:    
for i1=1:2,
  if i1==1,   vari=x'*x/N; Pk=Pkx; end;
  if i1==2,   vari=y'*y/N; Pk=Pky; end;
  P    = (Pk(:,1)+Pk(:,2))/2;   % initial spectrum estimate
  Ptemp= zeros(N,1);
  P1   = zeros(N,1);
  tol  = .0005*vari/N;          % usually within 'tol'erance in about three iterations, see equations from [2] (P&W pp 368-370).   
  a    = vari*(1-V);
  while sum(abs(P-P1)/N)>tol            
    b=(P*ones(1,k))./(P*V'+ones(N,1)*a'); % weights
    wk=(b.^2).*(ones(N,1)*V');            % new spectral estimate
    P1=(sum(wk'.*Pk')./ sum(wk'))';
    Ptemp=P1; P1=P; P=Ptemp;              % swap P and P1
  end
  if i1==1, 
    fkx=sqrt(k)*sqrt(wk).*fkx./repmat(sum(sqrt(wk'))',1,k);
    Fx=P;  %Power density spectral estimate of x
  end;
  if i1==2, 
    fky=sqrt(k)*sqrt(wk).*fky./repmat(sum(sqrt(wk'))',1,k);
    Fy=P;  %Power density spectral estimate of y
  end;
end;
%As a check, the quantity sum(abs(fkx(pls,:))'.^2) is the same as Fx and
%the spectral estimate from pmtmPH.

%Compute coherence
Cxy= sum([fkx.*conj(fky)]');
ph = angle(Cxy)*180/pi;
c  = abs(Cxy)./sqrt(sum(abs(fkx').^2).*sum(abs(fky').^2));

%correct for the bias of the estimate
if qbias==1,
  c=cohbias(v,c)'; 
end;

%Phase uncertainty estimates via Monte Carlo analysis. 
if confn>1,
  cb=cohbias(v,c)';
for iter=1:confn;
  if rem(iter,10)==0, disp(['phase confidence iteration: ',num2str(iter)]); end;
  fx=fft(randn(size(x))+1);
  fx=fx/sum(abs(fx));
  fy=fft(randn(size(y))+1); 
  fy=fy/sum(abs(fy));
  ys =real(ifft(fy.*sqrt(1-cb'.^2)));    
  ys =ys+real(ifft(fx.*cb')); 
  xs =real(ifft(fx));
  [si, ciph(iter,:), phi(iter,:)]=cmtm(xs,ys,dt,NW);
end;
  pl=round(.975*iter);     %sorting and averaging to determine confidence levels.
  phi=sort(phi);  
  phi=[phi(pl,:); -phi(iter-pl+1,:)];
  phi=mean(phi);
  phi=conv(phi(1:end),[1 1 1]/3);  
  phi=phi(2:end-1);
else,
  phi=zeros(size(pls)); 
end;
  
%Cut to one-sided funtions
c = c(pls);
s = s(pls)';
ph=ph(pls);
phl=ph-phi;
phu=ph+phi;

%Coherence confidence level
ci=cohconf(v,.95);  %not corrected for bias, this is conservative.
ci=ci*ones(size(c));

%plotting
if qplot==1,
  %coherence
  figure(gcf); clf;
  subplot(211); hold on;
  plot(s,c);
  h=ylabel('coherence');
  h=xlabel('frequency');
  plot(s,ci,'k--');
  pl=find(c>ci(1));
  title(['mean is ',num2str(mean(c),2),'   ',num2str(100*length(pl)/length(c),2),'% of estimates above 95% confidence level']);
  axis tight; h=axis; axis([h(1:2) 0 1.025]);
  w  = NW/(dt*N);   %half-bandwidth of the dpss
  plot([s(1) h(2)],[1.02 1.02],'k');
  for ds=min(s):2*w:max(s);
    plot([ds ds],[.98 1.02],'k');
  end;
  
  %phase
  subplot(212); hold on;
  plot(s,ph);
  if confn>0,
    col=[.9 .9 .9];
    h=fill([s(1) s(1:end) fliplr([s(1:end) s(end)])],[phu(1) phl fliplr([phu phl(end)])],col);
    set(h,'edgecolor',col);

    pl=find(phu<=180); phu(pl)=-180;
    pl=find(phu> 180); phu(pl)=phu(pl)-360;
    phlt=-180*ones(size(phl));
    h=fill([s(1) s(1:end) fliplr([s(1:end) s(end)])],[phu(1) phlt fliplr([phu phlt(end)])],col);
    set(h,'edgecolor',col);

    pl=find(phl>=-180); phl(pl)=180;
    pl=find(phl< -180); phl(pl)=phl(pl)+360;
    phut=180*ones(size(phl));
    h=fill([s(1) s(1:end) fliplr([s(1:end) s(end)])],[phut(1) phl fliplr([phut phl(end)])],col);    
    set(h,'edgecolor',col);        
  end;
  h=plot(s,ph); 
  plot(s,zeros(size(s)),'k--');
  axis tight; h=axis; axis([h(1:2) -180 180]);
  h=xlabel('frequency'); 
  h=ylabel('phase'); 
end;