mCoh.m

function [T,R,S] = mCoh(x,in,out,given,Nbt);
%Usage: [T,R,S] = mCoh(x,in,out,given,Nbt);
%   multiple coherence R of Output channel OUT with input
%   channels IN of multivariate time-series X as a function
%   of period T (in units of samples; multiply by dt
%   to get actual period); N is number of FCs averaged over/band
%  can also have "given" channels

%   following code allows logarithmically spaced bands;
%   hard coded at present, not an option
[nch,npts] = size(x);
if((min(in) <= 0) | (min(out) <=0) | ...
   (max(in) > nch)| (max(out) > nch))
   T = -1;
   return
end

logBands = 1;
alpha = .8;
minBW = 10;
[iband] = setBands(npts,Nbt,logBands,alpha,minBW);

temp =zeros(nch,1);
temp(in) = 1; temp(out) = 1;
temp(given) = 1;...
all=find(temp);
nAll = length(all);
nIn = length(in);
nOut = length(out);
nGiven = length(given);
Kin = zeros(nIn,1);
for j = 1:nIn
   for k = 1:nAll
      if(in(j) == all(k))
         Kin(j) = k;
      end
   end
end
Kout = zeros(nOut,1);
for j = 1:nOut
   for k = 1:nAll
      if(out(j) == all(k))
         Kout(j) = k;
      end
   end
end
Kgiven = [];
for j = 1:nGiven
   for k = 1:nAll
      if(given(j) == all(k))
         Kgiven(j) = k;
      end
   end
end

X = x(all,:);
%   get rid of NaN's, demean, taper ends of TS
Nend = fix(.1*npts);
for ic = 1:nAll
   xMean = mean(X(ic,isnan(X(ic,:))==0));
   X(ic,isnan(X(ic,:))) = xMean;
%   X(ic,:) = taperEnds(X(ic,:),Nend);
end

%   fft
X = fft(X,[],2);

%  compute sdm
[SDM] = xProd(X,iband);

%  compute multivariate squared correlations
[R,S] = mvcorr(SDM,Kout,Kin,Kgiven);

%  compute periods
T = npts*(1./mean(iband,1));