# cohgrampt

## PURPOSE

Multi-taper time-frequency coherence - two point processes given as times

## SYNOPSIS

function [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr)

## DESCRIPTION

     Multi-taper time-frequency coherence - two point processes given as times

     Usage:

     [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr)
     Input:
     Note units have to be consistent. Thus, if movingwin is in seconds, Fs
     has to be in Hz. see chronux.m for more information.

           data1  (structure array of spike times with dimension trials; also accepts 1d array of spike times) -- required
           data2  (structure array of spike times with dimension trials; also accepts 1d array of spike times) -- required
           movingwin (in the form [window winstep] -- required
           params: structure with fields tapers, pad, Fs, fpass, err, trialave
           - 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
                        Note that T has to be equal to movingwin(1).

                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
               err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
                                       [0 p] or 0 - no error bars) - optional. Default 0.
               trialave (average over trials when 1, don't average when 0) - optional. Default 0
           fscorr   (finite size corrections, 0 (don't use finite size corrections) or
                     1 (use finite size corrections) - optional
                    (available only for spikes). Defaults 0.
     Output:
           C (magnitude of coherency time x frequencies x trials for trialave=0;
                  time x frequency for trialave=1)
           phi (phase of coherency time x frequencies x trials for no trial averaging;
                  time x frequency for trialave=1)
           S12 (cross spectrum - time x frequencies x trials for no trial averaging;
                  time x frequency for trialave=1)
           S1 (spectrum 1 - time x frequencies x trials for no trial averaging;
                  time x frequency for trialave=1)
           S2 (spectrum 2 - time x frequencies x trials for no trial averaging;
                  time x frequency for trialave=1)
           t (time)
           f (frequencies)
           zerosp (1 for windows and trials where spikes were absent (in either channel),zero otherwise)
           confC (confidence level for C at 1-p %) - only for err(1)>=1
           phistd - theoretical/jackknife (depending on err(1)=1/err(1)=2) standard deviation for phi
                    Note that phi + 2 phistd and phi - 2 phistd will give 95% confidence
                    bands for phi - only for err(1)>=1
           Cerr  (Jackknife error bars for C - use only for Jackknife - err(1)=2)

## CROSS-REFERENCE INFORMATION

This function calls:

- check_consistency Helper routine to check consistency of data dimensions
- 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
- coherencypt Multi-taper coherency - point process times
- 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

This function is called by:

- testscript function testscript(pname,direction,movingwin,segave,params,fscorr)

## SOURCE CODE

    0001 function [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr)
    0002 % Multi-taper time-frequency coherence - two point processes given as times
    0003 %
    0004 % Usage:
    0005 %
    0006 % [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr)
    0007 % Input:
    0008 % Note units have to be consistent. Thus, if movingwin is in seconds, Fs
    0009 % has to be in Hz. see chronux.m for more information.
    0010 %
    0011 %       data1  (structure array of spike times with dimension trials; also accepts 1d array of spike times) -- required
    0012 %       data2  (structure array of spike times with dimension trials; also accepts 1d array of spike times) -- required
    0013 %       movingwin (in the form [window winstep] -- required
    0014 %       params: structure with fields tapers, pad, Fs, fpass, err, trialave
    0015 %       - optional
    0016 %           tapers : precalculated tapers from dpss or in the one of the following
    0017 %                    forms:
    0018 %                   (1) A numeric vector [TW K] where TW is the
    0019 %                       time-bandwidth product and K is the number of
    0020 %                       tapers to be used (less than or equal to
    0021 %                       2TW-1).
    0022 %                   (2) A numeric vector [W T p] where W is the
    0023 %                       bandwidth, T is the duration of the data and p
    0024 %                       is an integer such that 2TW-p tapers are used. In
    0025 %                       this form there is no default i.e. to specify
    0026 %                       the bandwidth, you have to specify T and p as
    0027 %                       well. Note that the units of W and T have to be
    0028 %                       consistent: if W is in Hz, T must be in seconds
    0029 %                       and vice versa. Note that these units must also
    0030 %                       be consistent with the units of params.Fs: W can
    0031 %                       be in Hz if and only if params.Fs is in Hz.
    0032 %                       The default is to use form 1 with TW=3 and K=5
    0033 %                    Note that T has to be equal to movingwin(1).
    0034 %
    0035 %            pad            (padding factor for the FFT) - optional (can take values -1,0,1,2...).
    0036 %                    -1 corresponds to no padding, 0 corresponds to padding
    0037 %                    to the next highest power of 2 etc.
    0038 %                       e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
    0039 %                       to 512 points, if pad=1, we pad to 1024 points etc.
    0040 %                       Defaults to 0.
    0041 %           Fs   (sampling frequency) - optional. Default 1.
    0042 %           fpass    (frequency band to be used in the calculation in the form
    0043 %                                   [fmin fmax])- optional.
    0044 %                                   Default all frequencies between 0 and Fs/2
    0045 %           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
    0046 %                                   [0 p] or 0 - no error bars) - optional. Default 0.
    0047 %           trialave (average over trials when 1, don't average when 0) - optional. Default 0
    0048 %       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
    0049 %                 1 (use finite size corrections) - optional
    0050 %                (available only for spikes). Defaults 0.
    0051 % Output:
    0052 %       C (magnitude of coherency time x frequencies x trials for trialave=0;
    0053 %              time x frequency for trialave=1)
    0054 %       phi (phase of coherency time x frequencies x trials for no trial averaging;
    0055 %              time x frequency for trialave=1)
    0056 %       S12 (cross spectrum - time x frequencies x trials for no trial averaging;
    0057 %              time x frequency for trialave=1)
    0058 %       S1 (spectrum 1 - time x frequencies x trials for no trial averaging;
    0059 %              time x frequency for trialave=1)
    0060 %       S2 (spectrum 2 - time x frequencies x trials for no trial averaging;
    0061 %              time x frequency for trialave=1)
    0062 %       t (time)
    0063 %       f (frequencies)
    0064 %       zerosp (1 for windows and trials where spikes were absent (in either channel),zero otherwise)
    0065 %       confC (confidence level for C at 1-p %) - only for err(1)>=1
    0066 %       phistd - theoretical/jackknife (depending on err(1)=1/err(1)=2) standard deviation for phi
    0067 %                Note that phi + 2 phistd and phi - 2 phistd will give 95% confidence
    0068 %                bands for phi - only for err(1)>=1
    0069 %       Cerr  (Jackknife error bars for C - use only for Jackknife - err(1)=2)
    0070
    0071 if nargin < 3; error('Need data1 and data2 and window parameters'); end;
    0072 if nargin < 4; params=[]; end;
    0073
    0074 if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    0075     error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
    0076 end
    0077
    0078 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
    0079 if nargin < 5 || isempty(fscorr); fscorr=0; end;
    0080
    0081 if nargout > 10 && err(1)~=2;
    0082     error('Cerr computed only for Jackknife. Correct inputs and run again');
    0083 end;
    0084 if nargout > 8 && err(1)==0;
    0085     error('Errors computed only if err(1) is not equal to zero');
    0086 end;
    0087
    0088 [N,Ch]=check_consistency(data1,data2);
    0089 [mintime1,maxtime1]=minmaxsptimes(data1);
    0090 [mintime2,maxtime2]=minmaxsptimes(data2);
    0091 mintime=min(mintime1,mintime2);
    0092 maxtime=max(maxtime1,maxtime2);
    0093
    0094 tn=mintime+movingwin(1)/2:movingwin(2):maxtime-movingwin(1)/2;
    0095 Nwin=round(Fs*movingwin(1)); % number of samples in window
    0096 % Nstep=round(movingwin(2)*Fs); % number of samples to step through
    0097 nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
    0098 f=getfgrid(Fs,nfft,fpass); Nf=length(f);
    0099 params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers
    0100 nw=length(tn);
    0101 if trialave;
    0102    C=zeros(nw,Nf);
    0103    S12=zeros(nw,Nf);
    0104    S1=zeros(nw,Nf);
    0105    S2=zeros(nw,Nf);
    0106    phi=zeros(nw,Nf);
    0107    Cerr=zeros(2,nw,Nf);
    0108 %    phierr=zeros(2,nw,Nf);
    0109    phistd=zeros(nw,Nf);
    0110 else
    0111    C=zeros(nw,Nf,Ch);
    0112    S12=zeros(nw,Nf,Ch);
    0113    S1=zeros(nw,Nf,Ch);
    0114    S2=zeros(nw,Nf,Ch);
    0115    phi=zeros(nw,Nf,Ch);
    0116    Cerr=zeros(2,nw,Nf,Ch);
    0117 %    phierr=zeros(2,nw,Nf,Ch);
    0118    phistd=zeros(nw,Nf,Ch);
    0119 end;
    0120 zerosp=zeros(nw,Ch);
    0121
    0122 for n=1:nw;
    0123    t=linspace(tn(n)-movingwin(1)/2,tn(n)+movingwin(1)/2,Nwin);
    0124    datawin1=extractdatapt(data1,[t(1) t(end)]);datawin2=extractdatapt(data2,[t(1) t(end)]);
    0125    if nargout==11;
    0126      [c,ph,s12,s1,s2,f,zsp,confc,phie,cerr]=coherencypt(datawin1,datawin2,params,fscorr,t);
    0127 %      phierr(1,n,:,:)=squeeze(phie(1,:,:));
    0128 %      phierr(2,n,:,:)=squeeze(phie(2,:,:));
    0129      phistd(n,:,:)=phie;
    0130      Cerr(1,n,:,:)=squeeze(cerr(1,:,:));
    0131      Cerr(2,n,:,:)=squeeze(cerr(2,:,:));
    0132    elseif nargout==10;
    0133      [c,ph,s12,s1,s2,f,zsp,confc,phie]=coherencypt(datawin1,datawin2,params,fscorr,t);
    0134 %      phierr(1,n,:,:)=squeeze(phie(1,:,:));
    0135 %      phierr(2,n,:,:)=squeeze(phie(2,:,:));
    0136      phistd(n,:,:)=phie;
    0137    else
    0138      [c,ph,s12,s1,s2,f,zsp]=coherencypt(datawin1,datawin2,params,fscorr,t);
    0139    end;
    0140    C(n,:,:)=c;
    0141    phi(n,:,:)=ph;
    0142    S12(n,:,:)=s12;
    0143    S1(n,:,:)=s1;
    0144    S2(n,:,:)=s2;
    0145    zerosp(n,:)=zsp;
    0146 end;
    0147 t=tn;
    0148 C=squeeze(C); phi=squeeze(phi);S12=squeeze(S12); S1=squeeze(S1); S2=squeeze(S2);zerosp=squeeze(zerosp);
    0149 if nargout > 9; confC=confc; end;
    0150 if nargout==11;Cerr=squeeze(Cerr);end;
    0151 % if nargout==10; phierr=squeeze(phierr);end
    0152 if nargout==10; phistd=squeeze(phistd);end