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