cohgrampt

 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)

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