This routine computes the average multi-taper coherence for a given set of unequal length segments. It is
based on modifications to the Chronux routines. The segments are continuously structured in the
data matrix, with the segment boundaries given by markers. Below,
movingwin is used in a non-overlaping way to partition each segment into
various windows. Th coherence is evaluated for each window, and then the
window coherence estimates averaged. Further averaging is conducted by
repeating the process for each segment.
Inputs:
data = data( samples, channels )- here segments must be stacked
as explained in the email
movingwin = [window winstep] i.e length of moving
window and step size. Note that units here have
to be consistent with units of Fs. If Fs=1 (ie normalized)
then [window winstep]should be in samples, or else if Fs is
unnormalized then they should be in time (secs).
sMarkers = N x 2 array of segment start & stop marks. sMarkers(n, 1) = start
sample index; sMarkers(n,2) = stop sample index for the nth segment
params = see Chronux help on mtspecgramc
Output:
Cmn magnitude of coherency - frequencies x iChPairs
Phimn phase of coherency - frequencies x iChPairs
Smn cross spectrum - frequencies x iChPairs
Smm spectrum m - frequencies x channels
f frequencies x 1
ConfC 1 x iChPairs; confidence level for Cmn at 1-p % - only for err(1)>=1
PhiStd frequency x iChPairs; error bars for phimn - only for err(1)>=1
Cerr 2 x frequency x iChPairs; Jackknife error bars for Cmn - use only for Jackknife - err(1)=2
Here iChPairs = indices corresponding to the off-diagonal terms of the
lower half matrix. iChPairs = 1 : nChannels*(nChannels-1)/2. So,
iChPairs=1,2,3,4,...correspond to C(2,1), C(3,1), C(3,2), C(4,1), etc.
The mapping can be obtained as follows:
C(i,j) = Cmn(:,k) where k = j + (1/2)*(i-1)*(i-2)
The above also applies to phimn, Smn
Note:
segment length >= NW/2 where NW = half bandwidth parameter (see dpss). So the power spectrum will
be computed only for those segments whose length > NW/2. For that reason, the routine returns the
indices for segments for which the spectra is computed. This check is
done here since pSpecgramAvg calls it.0001 function [Cmn,Phimn,Smn,Smm,f,ConfC,PhiStd,Cerr] = coherencyc_unequal_length_trials( data, movingwin, params, sMarkers ) 0002 0003 % This routine computes the average multi-taper coherence for a given set of unequal length segments. It is 0004 % based on modifications to the Chronux routines. The segments are continuously structured in the 0005 % data matrix, with the segment boundaries given by markers. Below, 0006 % movingwin is used in a non-overlaping way to partition each segment into 0007 % various windows. Th coherence is evaluated for each window, and then the 0008 % window coherence estimates averaged. Further averaging is conducted by 0009 % repeating the process for each segment. 0010 % 0011 % Inputs: 0012 % 0013 % data = data( samples, channels )- here segments must be stacked 0014 % as explained in the email 0015 % movingwin = [window winstep] i.e length of moving 0016 % window and step size. Note that units here have 0017 % to be consistent with units of Fs. If Fs=1 (ie normalized) 0018 % then [window winstep]should be in samples, or else if Fs is 0019 % unnormalized then they should be in time (secs). 0020 % sMarkers = N x 2 array of segment start & stop marks. sMarkers(n, 1) = start 0021 % sample index; sMarkers(n,2) = stop sample index for the nth segment 0022 % params = see Chronux help on mtspecgramc 0023 % 0024 % Output: 0025 % Cmn magnitude of coherency - frequencies x iChPairs 0026 % Phimn phase of coherency - frequencies x iChPairs 0027 % Smn cross spectrum - frequencies x iChPairs 0028 % Smm spectrum m - frequencies x channels 0029 % f frequencies x 1 0030 % ConfC 1 x iChPairs; confidence level for Cmn at 1-p % - only for err(1)>=1 0031 % PhiStd frequency x iChPairs; error bars for phimn - only for err(1)>=1 0032 % Cerr 2 x frequency x iChPairs; Jackknife error bars for Cmn - use only for Jackknife - err(1)=2 0033 % 0034 % Here iChPairs = indices corresponding to the off-diagonal terms of the 0035 % lower half matrix. iChPairs = 1 : nChannels*(nChannels-1)/2. So, 0036 % iChPairs=1,2,3,4,...correspond to C(2,1), C(3,1), C(3,2), C(4,1), etc. 0037 % The mapping can be obtained as follows: 0038 % 0039 % C(i,j) = Cmn(:,k) where k = j + (1/2)*(i-1)*(i-2) 0040 % 0041 % The above also applies to phimn, Smn 0042 % 0043 % Note: 0044 % segment length >= NW/2 where NW = half bandwidth parameter (see dpss). So the power spectrum will 0045 % be computed only for those segments whose length > NW/2. For that reason, the routine returns the 0046 % indices for segments for which the spectra is computed. This check is 0047 % done here since pSpecgramAvg calls it. 0048 0049 iwAvg = 1; % 0=no weighted average, 1=weighted average 0050 debug = 1; % will display intermediate calcs. 0051 0052 if nargin < 2; error('avgCoherence:: Need data and window parameters'); end; 0053 if nargin < 3; params=[]; end; 0054 [ tapers, pad, Fs, fpass, err, trialave, params ] = getparams( params ); 0055 if isempty( sMarkers ), error( 'avgCoherence:: Need Markers...' ); end 0056 % Not designed for "trialave" so set to 0 0057 params.trialave = 0; 0058 [ tapers, pad, Fs, fpass, err, trialave, params ] = getparams( params ); 0059 if nargout > 7 && err(1)~=2; 0060 error('avgCoherence:: Cerr computed only for Jackknife. Correct inputs and run again'); 0061 end; 0062 if nargout > 5 && err(1)==0; 0063 % Errors computed only if err(1) is nonzero. Need to change params and run again. 0064 error('avgCoherence:: When errors are desired, err(1) has to be non-zero.'); 0065 end; 0066 if size(data,2)==1, error('avgCoherence:: Need more than 1 channel to compute coherence'); end 0067 0068 % Set moving window parameters to no-overlapping 0069 if abs(movingwin(2) - movingwin(1)) >= 1e-6, disp( 'avgCoherence:: Warming: Window parameters for averaging should be non-overlapping. Set movingwin(2) = movingwin(1).' ); end 0070 0071 wLength = round( Fs * movingwin(1) ); % number of samples in window 0072 wStep = round( movingwin(2) * Fs ); % number of samples to step through 0073 0074 % Check whether window lengths satify segment length > NW/2 0075 if ( wLength < 2*tapers(1) ), error( 'avgCoherence:: movingwin(1) > 2*tapers(1)' ); end 0076 0077 % Left align segment markers for easier coding 0078 sM = ones( size( sMarkers, 1 ), 2 ); 0079 sM( :, 2 ) = sMarkers( :, 2 ) - sMarkers( :, 1 ) + 1; 0080 0081 % min-max segments 0082 Nmax = max( sM(:,2) ); Nmin = min( sM(:,2) ); 0083 if ( Nmin < 2*tapers(1) ), error( 'avgCoherence:: Smallest segment length > 2*tapers(1). Change taper settings' ); end 0084 0085 % max time-sample length will be the window length. 0086 nfft = 2^( nextpow2( wLength ) + pad ); 0087 [ f, findx ] = getfgrid( Fs, nfft, fpass); 0088 0089 % Precompute all the tapers 0090 sTapers = tapers; 0091 sTapers = dpsschk( sTapers, wLength, Fs ); % compute tapers for window length 0092 0093 nChannels = size( data, 2 ); 0094 nSegments = size( sMarkers, 1 ); 0095 iChPairs = ceil( nChannels*(nChannels-1)/2 ); 0096 0097 if debug 0098 disp( ['Window Length = ' num2str(wLength)] ); 0099 disp( ['Window Step = ' num2str(wStep)] ); 0100 disp( ' ' ); 0101 end 0102 0103 % 0104 % coherr outputs such that: 0105 % confc is has dimensions [1 size(cmn,2)] => confc = 1 x iChPairs 0106 % phistd has dimensions [f size(cmn,2)] => phistd = frequencies x iChPairs = size( cmn ) 0107 % cerr has dimensions [2 size(cmn)] => cerr = 2 x frequencies x iChPairs 0108 % 0109 cerr = zeros( 2, length(f), iChPairs ); confc = zeros(1,iChPairs); phistd=zeros( length(f), iChPairs ); 0110 Cerr = zeros( 2, length(f), iChPairs ); ConfC = zeros(1,iChPairs); PhiStd=zeros( length(f), iChPairs ); 0111 0112 %serr = zeros( 2, length(f), nChannels ); 0113 smm = zeros( length(f), nChannels ); 0114 smn = zeros( length(f), iChPairs ); cmn=smn; phimn=smn; smn = complex( smn, smn ); 0115 Smm = smm; Smn = smn; Cmn = cmn; Phimn = phimn; 0116 nWins = 0; 0117 for sg = 1 : nSegments 0118 % Window lengths & steps fixed above 0119 % For the given segment, compute the positions & number of windows 0120 N = sM(sg,2); 0121 wStartPos = 1 : wStep : ( N - wLength + 1 ); 0122 nWindows = length( wStartPos ); 0123 if nWindows 0124 nWins = nWins + nWindows; % for averaging purposes 0125 0126 w=zeros(nWindows,2); 0127 for n = 1 : nWindows 0128 w(n,:) = [ wStartPos(n), (wStartPos(n) + wLength - 1) ]; % nWindows x 2. just like segment end points 0129 end 0130 0131 % Shift window limits back to original sample-stamps 0132 w(:, 1) = w(:,1) + (sMarkers( sg, 1 ) - 1); 0133 w(:, 2) = w(:,2) + (sMarkers( sg, 1 ) - 1); 0134 0135 if debug 0136 disp( ['Segment Start/Stop = ' num2str( w(1,1) ) ' ' num2str( w(end,2) ) ] ); 0137 disp( ['Min / Max Window Positions = ' num2str( min(w(:,1)) ) ' ' num2str( max(w(:,1)) ) ] ); 0138 disp( ['Total Number of Windows = ' num2str(nWindows) ]); 0139 disp( ' ' ); 0140 end 0141 0142 % Pile up window segments similar to segment pileup 0143 wData = zeros( wLength, nChannels, nWindows ); %initialize to avoid fragmentation 0144 for n = 1:nWindows 0145 %wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ), 'constant' ); 0146 wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ) ); 0147 end 0148 0149 % J1 = frequency x taper x nWindows 0150 % J2 = frequency x taper x nWindows x nChannels 0151 J2 = zeros( length(f), tapers(2), nWindows, nChannels ); J2 = complex( J2, J2 ); 0152 for c = 1 : nChannels 0153 J1 = mtfftc( squeeze(wData( :, c, : )), sTapers, nfft, Fs ); % FFT for the tapered data 0154 J2( :, :, :, c ) = J1(findx,:,:); 0155 end 0156 % J2 = frequency x taper x nWindows x nChannels 0157 % Inner mean = Average over tapers => frequency x nWindows x nChannels 0158 % Outer mean = Average over windows => frequency x nChannels 0159 % smm = diagonal terms, ie power spectrum 0160 % 0161 dim1 = [length(f), nWindows, nChannels]; 0162 dim2 = [length(f), nChannels]; 0163 % s = frequency x nChannels 0164 smm = reshape( squeeze( mean( reshape( squeeze( mean( conj(J2).*J2, 2 ) ), dim1), 2 ) ), dim2 ); 0165 0166 % 0167 % Compute only the lower off-diagonal terms 0168 % smn = Cross Spectrum terms = complex 0169 % cmn = abs( coherence ); phimn = phase( coherence ) 0170 % 0171 0172 % 0173 % coherr outputs such that: 0174 % confc is has dimensions [1 size(cmn,2)] => confc = 1 x iChPairs 0175 % phistd has dimensions [f size(cmn,2)] => phistd = frequencies x iChPairs = size( cmn ) 0176 % cerr has dimensions [2 size(cmn)] => cerr = 2 x frequencies x iChPairs 0177 % 0178 cerr = zeros( 2, length(f), iChPairs ); confc = zeros(1,iChPairs); phistd=zeros( length(f), iChPairs ); 0179 0180 dim = [length(f), tapers(2), nWindows]; 0181 id = 1; 0182 for m=2:nChannels 0183 Jm = reshape( squeeze( J2(:,:,:,m) ), dim ); % frequency x taper x nWindows 0184 for n=1:m-1 % since we want off-diagonal terms only 0185 Jn = reshape( squeeze( J2(:,:,:,n) ), dim ); % frequency x taper x nWindows 0186 % 0187 % Average the Cross-Spectrum, Smn, over the windows 0188 % smn = complex 0189 % First average over tapers, then over windows 0190 % 0191 smn(:,id) = squeeze( mean( squeeze( mean( conj(Jm).*Jn, 2 ) ), 2 ) ); % frequency x iChPairs 0192 % 0193 % Coh = Coherence = complex = size( smn ) = frequency x iChPairs 0194 % 0195 Coh = smn(:,id) ./ sqrt( smm(:,m) .* smm(:,n) ); 0196 cmn(:,id) = abs(Coh); % frequencies x iChPairs 0197 phimn(:,id) = angle(Coh); % frequencies x iChPairs 0198 0199 % Since we've averaged over segments, set trialave = 1 0200 % 0201 % coherr outputs: 0202 % confc is has dimensions [1 size(Cmn(:,1),2)] => confC = 1 x iChPairs 0203 % phierr has dimensions [f size(Cmn(:,1),2)] => phistd = frequencies x iChPairs = size( Cmn ) 0204 % cerr has dimensions [2 size(Cmn(:,1))] => Cerr = 2 x frequencies x iChPairs 0205 % 0206 % Now treat the various "windowed data" as "trials" 0207 [ cconfc, cphistd, ccerr ] = coherr( cmn(:,id), Jm, Jn, err, 1 ); 0208 cerr(:,:,id ) = ccerr; 0209 confc(id) = cconfc; 0210 %size(cphistd), size(phistd) 0211 0212 phistd(:,id) = cphistd; % frequencies x iChPairs 0213 0214 id = id + 1; 0215 end 0216 end 0217 0218 if iwAvg 0219 % Segment Weighted error estimates. 0220 Smm = Smm + nWindows*smm; 0221 Smn = Smn + nWindows*smn; 0222 Cmn = Cmn + nWindows*cmn; 0223 Phimn = Phimn + nWindows*phimn; 0224 PhiStd = PhiStd + nWindows*phistd; 0225 ConfC = ConfC + nWindows*confc; 0226 Cerr = Cerr + nWindows*cerr; 0227 else 0228 Smm = Smm + smm; 0229 Smn = Smn + smn; 0230 Cmn = Cmn + cmn; 0231 Phimn = Phimn + phimn; 0232 PhiStd = PhiStd + phistd; 0233 ConfC = ConfC + confc; 0234 Cerr = Cerr + cerr; 0235 end 0236 0237 else 0238 if debug, disp(['avgCoherence:: Zero windows for segment: ' num2str(sg) ]); end 0239 end 0240 end 0241 0242 % Segment Weighted error estimates. 0243 % Only over those that had non-zero windows 0244 if nWins && iwAvg 0245 Smn=Smn/nWins; Smm=Smm/nWins; Cmn=Cmn/nWins; Phimn=Phimn/nWins; PhiStd=PhiStd/nWins; ConfC=ConfC/nWins; Cerr=Cerr/nWins; 0246 end 0247 if ~nWins 0248 if debug, disp(['avgCoherence:: No segment long enough with movingwin parameters found. Reduce movingwin.' ]); end 0249 end 0250 0251 0252 0253