


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 = 1; % will display intermediate calcs. 0034 0035 if nargin < 2; error('avgSpectrum:: Need data and window parameters'); end; 0036 if nargin < 3; params=[]; end; 0037 if isempty( sMarkers ), error( 'avgSpectrum:: 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('avgSpectrum:: 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