# coherencyc_unequal_length_trials

## PURPOSE

This routine computes the average multi-taper coherence for a given set of unequal length segments. It is

## SYNOPSIS

function [Cmn,Phimn,Smn,Smm,f,ConfC,PhiStd,Cerr] = coherencyc_unequal_length_trials( data, movingwin, params, sMarkers )

## DESCRIPTION

     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.

## CROSS-REFERENCE INFORMATION

This function calls:

- mtfftc Multi-taper fourier transform - continuous data
- coherr Function to compute lower and upper confidence intervals on the coherency
- dpsschk Helper function to calculate tapers and, if precalculated tapers are supplied,
- getfgrid Helper function that gets the frequency grid associated with a given fft based computation
- getparams Helper function to convert structure params to variables used by the

This function is called by:

## SOURCE CODE

    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