


Multi-taper cross-spectral matrix - another routine, allows for multiple trials and channels
Does not do confidence intervals. Also this routine always averages over trials - continuous process
Usage:
[Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatc(data,win,params)
Input:
Note units have to be consistent. See chronux.m for more information.
data (in form samples x channels x trials)
win (duration of non-overlapping window)
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]=CrossSpecMatc(data,win,params) 0002 % Multi-taper cross-spectral matrix - another routine, allows for multiple trials and channels 0003 % Does not do confidence intervals. Also this routine always averages over trials - continuous process 0004 % 0005 % Usage: 0006 % 0007 % [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMatc(data,win,params) 0008 % Input: 0009 % Note units have to be consistent. See chronux.m for more information. 0010 % data (in form samples x channels x trials) 0011 % win (duration of non-overlapping window) 0012 % params: structure with fields tapers, pad, Fs, fpass 0013 % - optional 0014 % tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional. If not 0015 % specified, use [NW K]=[3 5] 0016 % pad (padding factor for the FFT) - optional. Defaults to 0. 0017 % e.g. For N = 500, if PAD = 0, we pad the FFT 0018 % to 512 points; if PAD = 2, we pad the FFT 0019 % to 2048 points, etc. 0020 % Fs (sampling frequency) - optional. Default 1. 0021 % fpass (frequency band to be used in the calculation in the form 0022 % [fmin fmax])- optional. 0023 % Default all frequencies between 0 and Fs/2 0024 % Output: 0025 % Sc (cross spectral matrix frequency x channels x channels) 0026 % Cmat Coherence matrix frequency x channels x channels 0027 % Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency) 0028 % Cvec leading Eigenvector (frequency x channels) 0029 % Cent A different measure of total coherence: GM/AM of SV^2s 0030 % f (frequencies) 0031 d=ndims(data); 0032 if d<2, error('Need multidimensional array'); end 0033 if d==2, [N,C]=size(data); end; 0034 if d==3, [N,C,Ntr]=size(data); end; 0035 if nargin < 3; params=[]; end; 0036 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); 0037 clear err trialave params 0038 nwin=round(win*Fs); nfft=max(2^(nextpow2(nwin)+pad),nwin); 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]=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 end 0095