# CrossSpecMat

## PURPOSE

## SYNOPSIS

function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMat(data,win,params)

## DESCRIPTION

     Multi-taper cross-spectral matrix - another routine, this one allows for multiple trials and channels
     but does not do confidence intervals. Also this routine always averages over trials - continuous process

     Usage:

     [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMat(data,win,params)
     Input:
     Note units have to be consistent. See chronux.m for more information.
           data (in form samples x channels x trials)
           win  (duration of non-overlapping window)
           params: structure with fields tapers, pad, Fs, fpass
           - optional
               tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional.
                                                     If not specified, use [NW K]=[3 5]
                pad            (padding factor for the FFT) - optional. Defaults to 0.
                           e.g. For N = 500, if PAD = 0, we pad the FFT
                           to 512 points; if PAD = 2, we pad the FFT
                           to 2048 points, etc.
               Fs   (sampling frequency) - optional. Default 1.
               fpass    (frequency band to be used in the calculation in the form
                                       [fmin fmax])- optional.
                                       Default all frequencies between 0 and Fs/2
     Output:
           Sc (cross spectral matrix frequency x channels x channels)
           Cmat Coherence matrix frequency x channels x channels
           Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
           Cvec leading Eigenvector (frequency x channels)
           Cent A different measure of total coherence: GM/AM of SV^2s
           f (frequencies)

## CROSS-REFERENCE INFORMATION

This function calls:

- mtfftc Multi-taper fourier transform - continuous data
- dpsschk Helper function to calculate tapers and, if precalculated tapers are supplied,
- getfgrid Helper function that gets the frequency grid associated with a given fft based computation
- getparams Helper function to convert structure params to variables used by the

This function is called by:

## SOURCE CODE

    0001 function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMat(data,win,params)
    0002 %
    0003 %
    0004 % Multi-taper cross-spectral matrix - another routine, this one allows for multiple trials and channels
    0005 % but does not do confidence intervals. Also this routine always averages over trials - continuous process
    0006 %
    0007 % Usage:
    0008 %
    0009 % [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMat(data,win,params)
    0010 % Input:
    0011 % Note units have to be consistent. See chronux.m for more information.
    0012 %       data (in form samples x channels x trials)
    0013 %       win  (duration of non-overlapping window)
    0014 %       params: structure with fields tapers, pad, Fs, fpass
    0015 %       - optional
    0016 %           tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional.
    0017 %                                                 If not specified, use [NW K]=[3 5]
    0018 %            pad            (padding factor for the FFT) - optional. Defaults to 0.
    0019 %                       e.g. For N = 500, if PAD = 0, we pad the FFT
    0020 %                       to 512 points; if PAD = 2, we pad the FFT
    0021 %                       to 2048 points, etc.
    0022 %           Fs   (sampling frequency) - optional. Default 1.
    0023 %           fpass    (frequency band to be used in the calculation in the form
    0024 %                                   [fmin fmax])- optional.
    0025 %                                   Default all frequencies between 0 and Fs/2
    0026 % Output:
    0027 %       Sc (cross spectral matrix frequency x channels x channels)
    0028 %       Cmat Coherence matrix frequency x channels x channels
    0029 %       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
    0030 %       Cvec leading Eigenvector (frequency x channels)
    0031 %       Cent A different measure of total coherence: GM/AM of SV^2s
    0032 %       f (frequencies)
    0033 d=ndims(data);
    0034 if d<2, error('Need multidimensional array'); end
    0035 if d==2, [N,C]=size(data); end;
    0036 if d==3, [N,C,Ntr]=size(data); end;
    0037 if nargin < 3; params=[]; end;
    0038 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
    0039 clear err trialave params
    0040 nwin=round(win*Fs); nfft=2^(nextpow2(nwin)+pad);
    0041 [f,findx]=getfgrid(Fs,nfft,fpass);
    0042 tapers=dpsschk(tapers,nwin,Fs); % check tapers
    0043 Sc=zeros(length(findx),C,C);
    0044
    0045 Nwins=floor(N/nwin);
    0046
    0047 if d==3, % If there are multiple trials
    0048 for iwin=1:Nwins,
    0049     for i=1:Ntr,
    0050         data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:,i));
    0051         J1=mtfftc(detrend(data1),tapers,nfft,Fs);
    0052         J1=J1(findx,:,:);
    0053         for k=1:C,
    0054             for l=1:C,
    0055                 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
    0056             Sc(:,k,l)=Sc(:,k,l)+spec;
    0057             end
    0058         end
    0059     end
    0060 end
    0061 Sc=Sc/(Nwins*Ntr);
    0062 end
    0063
    0064 if d==2, % only one trial
    0065 for iwin=1:Nwins,
    0066         data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:));
    0067         J1=mtfftc(data1,tapers,nfft,Fs);
    0068         J1=J1(findx,:,:);
    0069         for k=1:C,
    0070             for l=1:C,
    0071             Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
    0072             end
    0073         end
    0074 end
    0075 Sc=Sc/Nwins;
    0076 end
    0077
    0078 Cmat=Sc;
    0079 Sdiag=zeros(length(findx),C);
    0080 for k=1:C,
    0081     Sdiag(:,k)=squeeze(Sc(:,k,k));
    0082 end
    0083
    0084 for k=1:C,
    0085     for l=1:C,
    0086         Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
    0087     end
    0088 end
    0089
    0090 Ctot=zeros(length(findx),1); Cent=Ctot;
    0091 Cvec=zeros(length(findx),C);
    0092 for i=1:length(findx),
    0093     [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
    0094     Ctot(i)=s(1).^2/sum(s.^2); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);
    0095     Cvec(i,:)=transpose(u(:,1));
    0096
    0097 end
    0098