


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 form [NW K] e.g [3 5]) -- optional. If not
specified, use [NW K]=[3 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 form [NW K] e.g [3 5]) -- optional. If not 0026 % specified, use [NW K]=[3 5] 0027 % pad (padding factor for the FFT) - optional (can take values -1,0,1,2...). 0028 % -1 corresponds to no padding, 0 corresponds to padding 0029 % to the next highest power of 2 etc. 0030 % e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT 0031 % to 512 points, if pad=1, we pad to 1024 points etc. 0032 % Defaults to 0. 0033 % Fs (sampling frequency) - optional. Default 1. 0034 % fpass (frequency band to be used in the calculation in the form 0035 % [fmin fmax])- optional. 0036 % Default all frequencies between 0 and Fs/2 0037 % Output: 0038 % Sc (cross spectral matrix frequency x channels x channels) 0039 % Cmat Coherence matrix frequency x channels x channels 0040 % Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency) 0041 % Cvec leading Eigenvector (frequency x channels) 0042 % Cent A different measure of total coherence: GM/AM of SV^2s 0043 % f (frequencies) 0044 d=ndims(data); 0045 if size(d,1)==1; error('Need multiple channels; are you sure your format is channels x trials ?');end; 0046 [C,Ntr]=size(data); 0047 mintime=0; 0048 if nargin < 3; [mintime,maxtime]=minmaxsptimes(data);clear mintime; 0049 else maxtime=T; end; 0050 if nargin < 4; params=[]; end; 0051 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); 0052 clear err trialave params 0053 Nwin=round(Fs*win); % number of samples in window 0054 nfft=max(2^(nextpow2(Nwin)+pad),Nwin); 0055 [f,findx]=getfgrid(Fs,nfft,fpass); 0056 tapers=dpsschk(tapers,Nwin,Fs); % check tapers 0057 twin=linspace(0,win,Nwin); % times of occurrence of "samples" within window - times of evaluation of tapers 0058 0059 Sc=zeros(length(findx),C,C); 0060 tn=mintime:win:maxtime-win; 0061 Nwins=length(tn); 0062 if d==3, % If there are multiple trials 0063 for iwin=1:Nwins, 0064 t=[tn(iwin) tn(iwin)+T]; 0065 for i=1:Ntr, 0066 data1=data(:,i); 0067 data1=extractdatapt(data1,t,1); % extract spike times in window,reset times to be relative to beginning of window 0068 J1=mtfftpt(data1,tapers,nfft,twin,f,findx); 0069 for k=1:C, 0070 for l=1:C, 0071 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 0072 Sc(:,k,l)=Sc(:,k,l)+spec; 0073 end 0074 end 0075 end 0076 end 0077 Sc=Sc/(Nwins*Ntr); 0078 end 0079 0080 if d==2, % only one trial 0081 for iwin=1:Nwins, 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 Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 0088 end 0089 end 0090 end 0091 Sc=Sc/Nwins; 0092 end 0093 0094 Cmat=Sc; 0095 Sdiag=zeros(length(findx),C); 0096 for k=1:C, 0097 Sdiag(:,k)=squeeze(Sc(:,k,k)); 0098 end 0099 0100 for k=1:C, 0101 for l=1:C, 0102 Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l))); 0103 end 0104 end 0105 0106 Ctot=zeros(length(findx),1); Cent=Ctot; 0107 Cvec=zeros(length(findx),C); 0108 for i=1:length(findx), 0109 [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s); 0110 Ctot(i)=s(1).^2/sum(s.^2); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2); 0111 Cvec(i,:)=transpose(u(:,1)); 0112 0113 end 0114