Home > chronux > spectral_analysis > continuous > coherencyc_unequal_length_trials.m

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: 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

Generated on Fri 28-Sep-2012 12:34:30 by m2html © 2005