# mtspectrumc_unequal_length_trials

## PURPOSE

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

## SYNOPSIS

function [ S, f, Serr ]= mtspectrumc_unequal_length_trials( data, movingwin, params, sMarkers )

## DESCRIPTION

     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

## CROSS-REFERENCE INFORMATION

This function calls:

- mtfftc Multi-taper fourier transform - continuous data
- 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
- specerr Function to compute lower and upper confidence intervals on the spectrum

This function is called by:

## SOURCE CODE

    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