mtspectrumc_unequal_length_trials

 This routine computes the multi-taper spectrum 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 spectrum is evaluated for each window, and then the
 window spectrum 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:

       S       frequency x channels
       f       frequencies x 1
       Serr    (error bars) only for err(1)>=1

0001 function [ S, f, Serr ]= mtspectrumc_unequal_length_trials( data, movingwin, params, sMarkers )
0002 
0003 % This routine computes the multi-taper spectrum 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 spectrum is evaluated for each window, and then the
0008 % window spectrum 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 %
0026 %       S       frequency x channels
0027 %       f       frequencies x 1
0028 %       Serr    (error bars) only for err(1)>=1
0029 %
0030 %
0031 
0032 iwAvg = 1; % 0=no weighted average, 1=weighted average
0033 debug = 0; % will display intermediate calcs.
0034 
0035 if nargin < 2; error('Unequal length trials:: Need data and window parameters'); end;
0036 if nargin < 3; params=[]; end;
0037 if isempty( sMarkers ), error( 'Unequal length trials:: Need Markers...' ); end
0038 [ tapers, pad, Fs, fpass, err, trialave, params ] = getparams( params );
0039 if nargout > 2 && err(1)==0; 
0040 %   Cannot compute error bars with err(1)=0. change params and run again.
0041     error('Unequal length trials:: When Serr is desired, err(1) has to be non-zero.');
0042 end;
0043 
0044 % Set moving window parameters to no-overlapping
0045 if abs(movingwin(2) - movingwin(1)) >= 1e-6, disp( 'avgSpectrum:: Warming: Window parameters for averaging should be non-overlapping. Set movingwin(2) = movingwin(1).' ); end
0046 
0047 wLength = round( Fs * movingwin(1) ); % number of samples in window
0048 wStep = round( movingwin(2) * Fs ); % number of samples to step through
0049 
0050 % Check whether window lengths satify segment length > NW/2
0051 if ( wLength < 2*tapers(1) ), error( 'avgSpectrum:: movingwin(1) > 2*tapers(1)' ); end
0052 
0053 % Left align segment markers for easier coding
0054 sM = ones( size( sMarkers, 1 ), 2 ); 
0055 sM( :, 2 ) = sMarkers( :, 2 ) - sMarkers( :, 1 ) + 1;
0056 
0057 % min-max segments
0058 Nmax = max( sM(:,2) ); Nmin = min( sM(:,2) );
0059 if ( Nmin < 2*tapers(1) ), error( 'avgSpectrum:: Smallest segment length > 2*tapers(1). Change taper settings' ); end
0060 
0061 % max time-sample length will be the window length.
0062 nfft = 2^( nextpow2( wLength ) + pad );
0063 [ f, findx ] = getfgrid( Fs, nfft, fpass); 
0064 
0065 % Precompute all the tapers
0066 sTapers = tapers;
0067 sTapers = dpsschk( sTapers, wLength, Fs ); % compute tapers for window length
0068 
0069 nChannels = size( data, 2 ); 
0070 nSegments = size( sMarkers, 1 );
0071 
0072 if debug
0073     disp( ['Window Length = ' num2str(wLength)] );
0074     disp( ['Window Step = ' num2str(wStep)] );
0075     disp( ' ' );
0076 end
0077 
0078 s = zeros( length(f), nChannels );
0079 serr = zeros( 2, length(f), nChannels );
0080 S = zeros( length(f), nChannels );
0081 Serr = zeros( 2, length(f), nChannels );
0082 nWins = 0;
0083 for sg = 1 : nSegments
0084     % Window lengths & steps fixed above
0085     % For the given segment, compute the positions & number of windows
0086     N = sM(sg,2); 
0087     wStartPos = 1 : wStep : ( N - wLength + 1 );
0088     nWindows = length( wStartPos );
0089     if nWindows
0090         nWins = nWins + nWindows; % for averaging purposes
0091 
0092         w=zeros(nWindows,2);
0093         for n = 1 : nWindows
0094             w(n,:) = [ wStartPos(n), (wStartPos(n) + wLength - 1) ]; % nWindows x 2. just like segment end points
0095         end
0096 
0097         % Shift window limits back to original sample-stamps
0098         w(:, 1) = w(:,1) + (sMarkers( sg, 1 ) - 1);
0099         w(:, 2) = w(:,2) + (sMarkers( sg, 1 ) - 1);
0100 
0101         if debug
0102             disp( ['Segment Start/Stop = ' num2str( w(1,1) ) ' ' num2str( w(end,2) ) ] );
0103             disp( ['Min / Max Window Positions = ' num2str( min(w(:,1)) ) ' ' num2str( max(w(:,1)) ) ] );
0104             disp( ['Total Number of Windows = ' num2str(nWindows) ]);
0105             disp( ' ' );
0106         end
0107 
0108         % Pile up window segments similar to segment pileup
0109         wData = zeros( wLength, nChannels, nWindows ); %initialize to avoid fragmentation
0110         for n = 1:nWindows
0111             %wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ), 'constant' );
0112             wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ) );
0113         end
0114 
0115         % J1 = frequency x taper x nWindows
0116         % J2 = frequency x taper x nWindows x nChannels
0117         J2 = zeros( length(f), tapers(2), nWindows, nChannels ); J2 = complex( J2, J2 );
0118         for c = 1 : nChannels
0119             J1 = mtfftc( squeeze(wData( :, c, : )), sTapers, nfft, Fs ); % FFT for the tapered data
0120             J2( :, :, :, c ) = J1(findx,:,:);
0121         end
0122         % J2 = frequency x taper x nWindows x nChannels
0123         % Inner mean = Average over tapers => frequency x nWindows x nChannels
0124         % Outer mean = Average over windows => frequency x nChannels
0125         dim1 = [length(f), nWindows, nChannels];
0126         dim2 = [length(f), nChannels];
0127         % s = frequency x nChannels
0128         s = reshape( squeeze( mean( reshape( squeeze( mean( conj(J2).*J2, 2 ) ), dim1), 2 ) ), dim2 );
0129 
0130         % Now treat the various "windowed data" as "trials"
0131         % serr = 2 x frequency x channels. Output from specerr = 2 x frequency x 1
0132         for c = 1 : nChannels
0133             serr( :, :, c ) = specerr( squeeze( s(:, c ) ), squeeze( J2(:,:,:, c ) ), err, 1 );
0134         end
0135         
0136         if iwAvg
0137             % Segment Weighted error estimates.
0138             S = S + nWindows*s;
0139             Serr = Serr + nWindows*serr;
0140         else
0141             S = S + s;
0142             Serr = Serr + serr;
0143         end
0144 
0145     else
0146         if debug, disp(['avgSpectrum:: Zero windows for segment: ' num2str(sg) ]); end
0147     end
0148 end
0149 
0150 % Segment Weighted error estimates.
0151 % Only over those that had non-zero windows
0152 if nWins && iwAvg
0153     S=S/nWins; Serr=Serr/nWins;
0154 end
0155 if ~nWins
0156     if debug, disp(['avgCoherence:: No segment long enough with movingwin parameters found. Reduce movingwin.' ]); end
0157 end
0158 
0159 
0160 
0161