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 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. 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 one of the following 0015 % forms: 0016 % (1) A numeric vector [TW K] where TW is the 0017 % time-bandwidth product and K is the number of 0018 % tapers to be used (less than or equal to 0019 % 2TW-1). 0020 % (2) A numeric vector [W T p] where W is the 0021 % bandwidth, T is the duration of the data and p 0022 % is an integer such that 2TW-p tapers are used. In 0023 % this form there is no default i.e. to specify 0024 % the bandwidth, you have to specify T and p as 0025 % well. Note that the units of W and T have to be 0026 % consistent: if W is in Hz, T must be in seconds 0027 % and vice versa. Note that these units must also 0028 % be consistent with the units of params.Fs: W can 0029 % be in Hz if and only if params.Fs is in Hz. 0030 % The default is to use form 1 with TW=3 and K=5 0031 % 0032 % pad (padding factor for the FFT) - optional. Defaults to 0. 0033 % e.g. For N = 500, if PAD = 0, we pad the FFT 0034 % to 512 points; if PAD = 2, we pad the FFT 0035 % to 2048 points, etc. 0036 % Fs (sampling frequency) - optional. Default 1. 0037 % fpass (frequency band to be used in the calculation in the form 0038 % [fmin fmax])- optional. 0039 % Default all frequencies between 0 and Fs/2 0040 % Output: 0041 % Sc (cross spectral matrix frequency x channels x channels) 0042 % Cmat Coherence matrix frequency x channels x channels 0043 % Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency) 0044 % Cvec leading Eigenvector (frequency x channels) 0045 % Cent A different measure of total coherence: GM/AM of SV^2s 0046 % f (frequencies) 0047 d=ndims(data); 0048 if d<2, error('Need multidimensional array'); end 0049 if d==2, [N,C]=size(data); end; 0050 if d==3, [N,C,Ntr]=size(data); end; 0051 if nargin < 3; params=[]; end; 0052 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); 0053 clear err trialave params 0054 nwin=round(win*Fs); nfft=max(2^(nextpow2(nwin)+pad),nwin); 0055 [f,findx]=getfgrid(Fs,nfft,fpass); 0056 tapers=dpsschk(tapers,nwin,Fs); % check tapers 0057 Sc=zeros(length(findx),C,C); 0058 0059 Nwins=floor(N/nwin); 0060 0061 if d==3, % If there are multiple trials 0062 for iwin=1:Nwins, 0063 for i=1:Ntr, 0064 data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:,i)); 0065 J1=mtfftc(detrend(data1),tapers,nfft,Fs); 0066 J1=J1(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=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:)); 0081 J1=mtfftc(data1,tapers,nfft,Fs); 0082 J1=J1(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)/sum(s); Cent(i)=exp(mean(log(s)))/mean(s); 0109 Cvec(i,:)=transpose(u(:,1)); 0110 end 0111