


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