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)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;