


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 - continuous process
Usage:
[Sc,Cmat,Ctot,f]=CrossSpecMat(data,win,params)
Input:
Note units have to be consistent. See chronux.m for more information.
data (in form samples x channels x trials)
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]=CrossSpecMat(data,win,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 over trials - continuous process 0006 % 0007 % Usage: 0008 % 0009 % [Sc,Cmat,Ctot,f]=CrossSpecMat(data,win,params) 0010 % Input: 0011 % Note units have to be consistent. See chronux.m for more information. 0012 % data (in form samples x channels x trials) 0013 % params: structure with fields tapers, pad, Fs, fpass 0014 % - optional 0015 % tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional. If not 0016 % specified, use [NW K]=[3 5] 0017 % pad (padding factor for the FFT) - optional. Defaults to 0. 0018 % e.g. For N = 500, if PAD = 0, we pad the FFT 0019 % to 512 points; if PAD = 2, we pad the FFT 0020 % to 2048 points, etc. 0021 % Fs (sampling frequency) - optional. Default 1. 0022 % fpass (frequency band to be used in the calculation in the form 0023 % [fmin fmax])- optional. 0024 % Default all frequencies between 0 and Fs/2 0025 % Output: 0026 % Sc (cross spectral matrix frequency x channels x channels) 0027 % Cmat Coherence matrix frequency x channels x channels 0028 % Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency) 0029 % Cvec leading Eigenvector (frequency x channels) 0030 % Cent A different measure of total coherence: GM/AM of SV^2s 0031 % f (frequencies) 0032 d=ndims(data); 0033 if d<2, error('Need multidimensional array'); end 0034 if d==2, [N,C]=size(data); end; 0035 if d==3, [N,C,Ntr]=size(data); end; 0036 if nargin < 3; params=[]; end; 0037 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); 0038 nwin=round(win*Fs); nfft=2^(nextpow2(nwin)+pad); 0039 [f,findx]=getfgrid(Fs,nfft,fpass); 0040 tapers=dpsschk(tapers,nwin,Fs); % check tapers 0041 Sc=zeros(length(findx),C,C); 0042 0043 Nwins=floor(N/nwin); 0044 0045 if d==3, % If there are multiple trials 0046 for iwin=1:Nwins, 0047 for i=1:Ntr, 0048 data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:,i)); 0049 J1=mtfftc(detrend(data1),tapers,nfft,Fs); 0050 J1=J1(findx,:,:); 0051 for k=1:C, 0052 for l=1:C, 0053 spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 0054 Sc(:,k,l)=Sc(:,k,l)+spec; 0055 end 0056 end 0057 end 0058 end 0059 Sc=Sc/(Nwins*Ntr); 0060 end 0061 0062 if d==2, % only one trial 0063 for iwin=1:Nwins, 0064 data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:)); 0065 J1=mtfftc(data1,tapers,nfft,Fs); 0066 J1=J1(findx,:,:); 0067 for k=1:C, 0068 for l=1:C, 0069 Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2)); 0070 end 0071 end 0072 end 0073 Sc=Sc/Nwins; 0074 end 0075 0076 Cmat=Sc; 0077 Sdiag=zeros(length(findx),C); 0078 for k=1:C, 0079 Sdiag(:,k)=squeeze(Sc(:,k,k)); 0080 end 0081 0082 for k=1:C, 0083 for l=1:C, 0084 Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l))); 0085 end 0086 end 0087 0088 Ctot=zeros(length(findx),1); Cent=Ctot; 0089 Cvec=zeros(length(findx),C); 0090 for i=1:length(findx), 0091 [u s v]=svd(squeeze(Sc(i,:,:)));s=diag(s); 0092 Ctot(i)=s(1).^2/sum(s.^2); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2); 0093 Cvec(i,:)=transpose(u(:,1)); 0094 0095 end 0096