# CrossSpecMatpt

## PURPOSE

## SYNOPSIS

function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,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 - point process as times

     Usage:

     [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
     Input:
     Note units have to be consistent. See chronux.m for more information.
           data (as a struct array with dimensions channels x trials) - note
           that times of measurement have to be consistent, we assume all
           times are specified relative to the start time of the trials which
           are taken to be zero.
           win  (duration of non-overlapping window)
           trialduration (since it is not possible to infer trial duration
                   from spike times, this is an optional argument. If not specified
                   the routine uses the minimum and maximum spike time (across all
                   channels and trials) as the window of calculation.) -
                   optional
           params: structure with fields tapers, pad, Fs, fpass
           - optional
               tapers : precalculated tapers from dpss or in the one of the following
                        forms:
                        (1) A numeric vector [TW K] where TW is the
                            time-bandwidth product and K is the number of
                            tapers to be used (less than or equal to
                            2TW-1).
                        (2) A numeric vector [W T p] where W is the
                            bandwidth, T is the duration of the data and p
                            is an integer such that 2TW-p tapers are used. In
                            this form there is no default i.e. to specify
                            the bandwidth, you have to specify T and p as
                            well. Note that the units of W and T have to be
                            consistent: if W is in Hz, T must be in seconds
                            and vice versa. Note that these units must also
                            be consistent with the units of params.Fs: W can
                            be in Hz if and only if params.Fs is in Hz.
                            The default is to use form 1 with TW=3 and K=5

                pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
                        -1 corresponds to no padding, 0 corresponds to padding
                        to the next highest power of 2 etc.
                           e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
                           to 512 points, if pad=1, we pad to 1024 points etc.
                           Defaults to 0.
               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:

- 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
- extractdatapt Extract segements of spike times between t(1) and t(2)
- minmaxsptimes Find the minimum and maximum of the spike times in each channel
- mtfftpt Multi-taper fourier transform for point process given as times

This function is called by:

## SOURCE CODE

    0001 function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,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
    0006 % over trials - point process as times
    0007 %
    0008 % Usage:
    0009 %
    0010 % [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatpt(data,win,T,params)
    0011 % Input:
    0012 % Note units have to be consistent. See chronux.m for more information.
    0013 %       data (as a struct array with dimensions channels x trials) - note
    0014 %       that times of measurement have to be consistent, we assume all
    0015 %       times are specified relative to the start time of the trials which
    0016 %       are taken to be zero.
    0017 %       win  (duration of non-overlapping window)
    0018 %       trialduration (since it is not possible to infer trial duration
    0019 %               from spike times, this is an optional argument. If not specified
    0020 %               the routine uses the minimum and maximum spike time (across all
    0021 %               channels and trials) as the window of calculation.) -
    0022 %               optional
    0023 %       params: structure with fields tapers, pad, Fs, fpass
    0024 %       - optional
    0025 %           tapers : precalculated tapers from dpss or in the one of the following
    0026 %                    forms:
    0027 %                    (1) A numeric vector [TW K] where TW is the
    0028 %                        time-bandwidth product and K is the number of
    0029 %                        tapers to be used (less than or equal to
    0030 %                        2TW-1).
    0031 %                    (2) A numeric vector [W T p] where W is the
    0032 %                        bandwidth, T is the duration of the data and p
    0033 %                        is an integer such that 2TW-p tapers are used. In
    0034 %                        this form there is no default i.e. to specify
    0035 %                        the bandwidth, you have to specify T and p as
    0036 %                        well. Note that the units of W and T have to be
    0037 %                        consistent: if W is in Hz, T must be in seconds
    0038 %                        and vice versa. Note that these units must also
    0039 %                        be consistent with the units of params.Fs: W can
    0040 %                        be in Hz if and only if params.Fs is in Hz.
    0041 %                        The default is to use form 1 with TW=3 and K=5
    0042 %
    0043 %            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
    0044 %                    -1 corresponds to no padding, 0 corresponds to padding
    0045 %                    to the next highest power of 2 etc.
    0046 %                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
    0047 %                       to 512 points, if pad=1, we pad to 1024 points etc.
    0048 %                       Defaults to 0.
    0049 %           Fs   (sampling frequency) - optional. Default 1.
    0050 %           fpass    (frequency band to be used in the calculation in the form
    0051 %                                   [fmin fmax])- optional.
    0052 %                                   Default all frequencies between 0 and Fs/2
    0053 % Output:
    0054 %       Sc (cross spectral matrix frequency x channels x channels)
    0055 %       Cmat Coherence matrix frequency x channels x channels
    0056 %       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
    0057 %       Cvec leading Eigenvector (frequency x channels)
    0058 %       Cent A different measure of total coherence: GM/AM of SV^2s
    0059 %       f (frequencies)
    0060 d=ndims(data);
    0061 if size(d,1)==1; error('Need multiple channels; are you sure your format is channels x trials ?');end;
    0062 [C,Ntr]=size(data);
    0063 mintime=0;
    0064 if nargin < 3; [mintime,maxtime]=minmaxsptimes(data);clear mintime;
    0065 else maxtime=T; end;
    0066 if nargin < 4; params=[]; end;
    0067 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
    0068 clear err trialave params
    0069 Nwin=round(Fs*win); % number of samples in window
    0070 nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
    0071 [f,findx]=getfgrid(Fs,nfft,fpass);
    0072 tapers=dpsschk(tapers,Nwin,Fs); % check tapers
    0073 twin=linspace(0,win,Nwin); % times of occurrence of "samples" within window - times of evaluation of tapers
    0074
    0075 Sc=zeros(length(findx),C,C);
    0076 tn=mintime:win:maxtime-win;
    0077 Nwins=length(tn);
    0078 if d==3, % If there are multiple trials
    0079 for iwin=1:Nwins,
    0080     t=[tn(iwin) tn(iwin)+T];
    0081     for i=1:Ntr,
    0082         data1=data(:,i);
    0083         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
    0084         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
    0085         for k=1:C,
    0086             for l=1:C,
    0087                 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
    0088             Sc(:,k,l)=Sc(:,k,l)+spec;
    0089             end
    0090         end
    0091     end
    0092 end
    0093 Sc=Sc/(Nwins*Ntr);
    0094 end
    0095
    0096 if d==2, % only one trial
    0097 for iwin=1:Nwins,
    0098         data1=data(:,i);
    0099         data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window
    0100         J1=mtfftpt(data1,tapers,nfft,twin,f,findx);
    0101         for k=1:C,
    0102             for l=1:C,
    0103             Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
    0104             end
    0105         end
    0106 end
    0107 Sc=Sc/Nwins;
    0108 end
    0109
    0110 Cmat=Sc;
    0111 Sdiag=zeros(length(findx),C);
    0112 for k=1:C,
    0113     Sdiag(:,k)=squeeze(Sc(:,k,k));
    0114 end
    0115
    0116 for k=1:C,
    0117     for l=1:C,
    0118         Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
    0119     end
    0120 end
    0121
    0122 Ctot=zeros(length(findx),1); Cent=Ctot;
    0123 Cvec=zeros(length(findx),C);
    0124 for i=1:length(findx),
    0125     [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
    0126     Ctot(i)=s(1)/sum(s); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);
    0127     Cvec(i,:)=transpose(u(:,1));
    0128
    0129 end
    0130