# cohgrampb

## PURPOSE

Multi-taper time-frequency coherence,cross-spectrum and individual spectra - two binned point processes

## SYNOPSIS

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

## DESCRIPTION

     Multi-taper time-frequency coherence,cross-spectrum and individual spectra - two binned point processes

     Usage:

     [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampb(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 (binned point process data in form samples x trials) -- required
           data2 (binned point process data in form samples x trials) -- 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 - jackknife/theoretical 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
- coherencycpb Multi-taper coherency,cross-spectrum and individual spectra - continuous and binned point process data
- coherencypb Multi-taper coherency,cross-spectrum and individual spectra - binned point process

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]=cohgrampb(data1,data2,movingwin,params,fscorr)
    0002 % Multi-taper time-frequency coherence,cross-spectrum and individual spectra - two binned point processes
    0003 %
    0004 % Usage:
    0005 %
    0006 % [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampb(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 (binned point process data in form samples x trials) -- required
    0012 %       data2 (binned point process data in form samples x trials) -- 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) -
    0048 %           optional. Default 0
    0049 %       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
    0050 %                 1 (use finite size corrections) - optional
    0051 %                (available only for spikes). Defaults 0.
    0052 % Output:
    0053 %       C (magnitude of coherency time x frequencies x trials for trialave=0;
    0054 %              time x frequency for trialave=1)
    0055 %       phi (phase of coherency time x frequencies x trials for no trial averaging;
    0056 %              time x frequency for trialave=1)
    0057 %       S12 (cross spectrum - time x frequencies x trials for no trial averaging;
    0058 %              time x frequency for trialave=1)
    0059 %       S1 (spectrum 1 - time x frequencies x trials for no trial averaging;
    0060 %              time x frequency for trialave=1)
    0061 %       S2 (spectrum 2 - time x frequencies x trials for no trial averaging;
    0062 %              time x frequency for trialave=1)
    0063 %       t (time)
    0064 %       f (frequencies)
    0065 %       zerosp (1 for windows and trials where spikes were absent (in either channel),zero otherwise)
    0066 %       confC (confidence level for C at 1-p %) - only for err(1)>=1
    0067 %       phistd - jackknife/theoretical standard deviation for phi - Note that
    0068 %                phi + 2 phistd and phi -2 phistd will give 95% confidence bands for phi
    0069 %                 - only for err(1)>=1
    0070 %       Cerr  (Jackknife error bars for C - use only for Jackknife - err(1)=2)
    0071
    0072 if nargin < 3; error('Need data1 and data2 and window parameters'); end;
    0073 if nargin < 4; params=[]; end;
    0074
    0075 if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    0076     error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
    0077 end
    0078
    0079 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
    0080 if nargin < 5 || isempty(fscorr); fscorr=0; end;
    0081
    0082 if nargout > 8 && err(1)==0;
    0083     error('When errors are desired, err(1) has to be non-zero.');
    0084 end;
    0085 if nargout > 10 && err(1)~=2;
    0086     error('Cerr computed only for Jackknife. Correct inputs and run again');
    0087 end;
    0088 [N,Ch]=check_consistency(data1,data2);
    0089
    0090 Nwin=round(Fs*movingwin(1)); % number of samples in window
    0091 Nstep=round(movingwin(2)*Fs); % number of samples to step through
    0092 nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
    0093 f=getfgrid(Fs,nfft,fpass);
    0094 Nf=length(f);
    0095 params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers
    0096
    0097 winstart=1:Nstep:N-Nwin+1;
    0098 nw=length(winstart);
    0099 if trialave;
    0100    C=zeros(nw,Nf);
    0101    S12=zeros(nw,Nf);
    0102    S1=zeros(nw,Nf);
    0103    S2=zeros(nw,Nf);
    0104    phi=zeros(nw,Nf);
    0105    Cerr=zeros(2,nw,Nf);
    0106 %   phierr=zeros(2,nw,Nf);
    0107    phistd=zeros(nw,Nf);
    0108 else
    0109    C=zeros(nw,Nf,Ch);
    0110    S12=zeros(nw,Nf,Ch);
    0111    S1=zeros(nw,Nf,Ch);
    0112    S2=zeros(nw,Nf,Ch);
    0113    phi=zeros(nw,Nf,Ch);
    0114    Cerr=zeros(2,nw,Nf,Ch);
    0115 %   phierr=zeros(2,nw,Nf,Ch);
    0116    phistd=zeros(nw,Nf,Ch);
    0117 end;
    0118 zerosp=zeros(nw,Ch);
    0119
    0120 for n=1:nw;
    0121    indx=winstart(n):winstart(n)+Nwin-1;
    0122    datawin1=data1(indx,:);datawin2=data2(indx,:);
    0123    if nargout==11;
    0124      [c,ph,s12,s1,s2,f,zsp,confc,phie,cerr]=coherencypb(datawin1,datawin2,params,fscorr);
    0125 %      phierr(1,n,:,:)=squeeze(phie(1,:,:));
    0126 %      phierr(2,n,:,:)=squeeze(phie(2,:,:));
    0127      phistd(n,:,:)=phie;
    0128      Cerr(1,n,:,:)=squeeze(cerr(1,:,:));
    0129      Cerr(2,n,:,:)=squeeze(cerr(2,:,:));
    0130    elseif nargout==10;
    0131      [c,ph,s12,s1,s2,f,zsp,confc,phie]=coherencypb(datawin1,datawin2,params,fscorr);
    0132 %      phierr(1,n,:,:)=squeeze(phie(1,:,:));
    0133 %      phierr(2,n,:,:)=squeeze(phie(2,:,:));
    0134      phistd(n,:,:)=phie;
    0135    else
    0136      [c,ph,s12,s1,s2,f,zsp]=coherencycpb(datawin1,datawin2,params,fscorr);
    0137    end;
    0138    C(n,:,:)=c;
    0139    phi(n,:,:)=ph;
    0140    S12(n,:,:)=s12;
    0141    S1(n,:,:)=s1;
    0142    S2(n,:,:)=s2;
    0143    zerosp(n,:)=zsp;
    0144 end;
    0145 C=squeeze(C); phi=squeeze(phi);S12=squeeze(S12); S1=squeeze(S1); S2=squeeze(S2);zerosp=squeeze(zerosp);
    0146 if nargout > 9; confC=confc; end;
    0147 if nargout==11;Cerr=squeeze(Cerr);end;
    0148 % if nargout==10; phierr=squeeze(phierr);end
    0149 if nargout==10; phistd=squeeze(phistd);end
    0150 winmid=winstart+round(Nwin/2);
    0151 t=winmid/Fs;