


SS_DEJITTER Aligns waveform peaks.
SPIKES = SS_DEJITTER(SPIKES) takes and returns a spike-sorting object
SPIKES and dejitters the spike waveforms using a 'center-of-mass'
method with a maximum shift of 3 (see below).
'Dejittering' refers to the registration of digitally sampled/thresholded
waveforms. The sampling/thresholding procedure introduces noise due to
(1) variation in relative timing between an analog waveform and the
sample clock and (2) variation in the time of threshold crossing due to
noise. The resulting misalignment even in otherwise identical analog
waveforms is known as jitter. Dejittering involves estimating the location
of a fiducial (e.g., the central peak) for each waveform and aligning these
markers across waveforms.
SPIKES = SS_DEJITTER(SPIKES, MAXSHIFT) also takes MAXSHIFT argument
(default: 3), specifying the maximum shift allowed during dejittering.
The number of samples per waveform will be altered by this function when
samples are made invalid by realignment (i.e., the amount of shift requires
extrapolating past the boundaries of the original waveform) and the returned
SPIKES object will have its waveforms clipped to exclude invalid regions.
The MAXSHIFT argument limits this by specifying the maximum allowed shift
(in samples). Waveforms requiring more than a MAXSHIFT samples realignment
are taken to be outliers and will not be shifted.
Dejittering requires the definition of threshold crossing time (threshT) and
threshold level vector (threshV) with both low and high thresholds (use +/-
Inf for either the low or the high threshold if it was not used in extracting
the waveforms). The 'center-of-mass' method considers the region following
threshT that remains below/above threshold and computes
(sum_t[t * (v_t - thresh)] / sum_t[(v_t - thresh)]).
to obtain the location of the peak. The interpolation to align the waveforms
is then performed using a cubic spline.
References:
Fee MS et al (1996). J. Neurosci Methods (69): 175-88
Sahani M (1999). PhD thesis, Pasadena, CA: Caltech.

0001 function spikes = ss_dejitter(spikes, maxshift) 0002 0003 % SS_DEJITTER Aligns waveform peaks. 0004 % SPIKES = SS_DEJITTER(SPIKES) takes and returns a spike-sorting object 0005 % SPIKES and dejitters the spike waveforms using a 'center-of-mass' 0006 % method with a maximum shift of 3 (see below). 0007 % 0008 % 'Dejittering' refers to the registration of digitally sampled/thresholded 0009 % waveforms. The sampling/thresholding procedure introduces noise due to 0010 % (1) variation in relative timing between an analog waveform and the 0011 % sample clock and (2) variation in the time of threshold crossing due to 0012 % noise. The resulting misalignment even in otherwise identical analog 0013 % waveforms is known as jitter. Dejittering involves estimating the location 0014 % of a fiducial (e.g., the central peak) for each waveform and aligning these 0015 % markers across waveforms. 0016 % 0017 % SPIKES = SS_DEJITTER(SPIKES, MAXSHIFT) also takes MAXSHIFT argument 0018 % (default: 3), specifying the maximum shift allowed during dejittering. 0019 % 0020 % The number of samples per waveform will be altered by this function when 0021 % samples are made invalid by realignment (i.e., the amount of shift requires 0022 % extrapolating past the boundaries of the original waveform) and the returned 0023 % SPIKES object will have its waveforms clipped to exclude invalid regions. 0024 % The MAXSHIFT argument limits this by specifying the maximum allowed shift 0025 % (in samples). Waveforms requiring more than a MAXSHIFT samples realignment 0026 % are taken to be outliers and will not be shifted. 0027 % 0028 % Dejittering requires the definition of threshold crossing time (threshT) and 0029 % threshold level vector (threshV) with both low and high thresholds (use +/- 0030 % Inf for either the low or the high threshold if it was not used in extracting 0031 % the waveforms). The 'center-of-mass' method considers the region following 0032 % threshT that remains below/above threshold and computes 0033 % (sum_t[t * (v_t - thresh)] / sum_t[(v_t - thresh)]). 0034 % to obtain the location of the peak. The interpolation to align the waveforms 0035 % is then performed using a cubic spline. 0036 % 0037 % References: 0038 % Fee MS et al (1996). J. Neurosci Methods (69): 175-88 0039 % Sahani M (1999). PhD thesis, Pasadena, CA: Caltech. 0040 0041 % Last Modified By: sbm on Thu Oct 6 20:30:26 2005 0042 0043 starttime = clock; 0044 0045 %%%%%%%%%% DEFAULTS 0046 if (nargin < 2), maxshift = 3; end; % default max shift 0047 0048 %%%%%%%%%% ARGUMENT CHECKING 0049 if (~isfield(spikes, 'waveforms') | (size(spikes.waveforms, 1) < 1)) 0050 error('SS:waveforms_undefined', 'The SS object does not contain any waveforms!'); 0051 elseif (~isfield(spikes, 'threshT')) 0052 error('SS:threshT_undefined', 'The SS object must define the sample index of threshold crossing.'); 0053 elseif (~isfield(spikes, 'threshV') | (length(spikes.threshV) ~= 2)) 0054 error('SS:threshV_error', 'The SS object must define both low and high thresholds. Use +/- Inf if only one of the two was used.'); 0055 elseif ((spikes.threshV(2) - spikes.threshV(1)) < 0) 0056 error('SS:threshV_illegal', 'The SS object high threshold must be greater than its low threshold.'); 0057 end 0058 0059 %%%%%%%%%% CONSTANTS 0060 spline_chunk = 5000; % # of splines to do at a time; keeps memory use down 0061 upsample = 10; % amount of upsampling to find extrema using 'peak' method 0062 numspikes = size(spikes.waveforms, 1); 0063 numsamples = size(spikes.waveforms, 2); 0064 0065 maxchunk = floor(((numspikes - 1) / spline_chunk)); 0066 0067 %%%%%%%%%% FIND FIDUCIALS 0068 % first, get a mask over the peaks 0069 isThreshLo = (spikes.waveforms(:, spikes.threshT) < spikes.threshV(1)); 0070 [h, w, pl, mask] = thresholded_peaks(spikes); 0071 0072 % next, subtract appropriate thresholds, making all values in the peak 0073 % have positive sign, since we're about to treat them as 'mass' 0074 waves = spikes.waveforms; 0075 waves(~isThreshLo, :) = waves(~isThreshLo, :) - spikes.threshV(2); 0076 waves( isThreshLo, :) = spikes.threshV(1) - waves( isThreshLo, :); 0077 waves = waves .* mask; 0078 0079 % now COM is straightforward 0080 fiducials = (waves * [1:numsamples]') ./ sum(waves, 2); 0081 0082 clear mask waves; % manual garbage collection 0083 0084 %%%%%%%%%% REALIGN FIDUCIALS 0085 % We line up the fiducials around the sample that requires the least shifting. 0086 target = round(mean(fiducials)); 0087 0088 % Determine shifted indices for each waveform 0089 shifts = fiducials - target; 0090 shifts(abs(shifts) > maxshift) = 0; % big shifts are outliers; don't alter these 0091 resample_inds = repmat([1:numsamples], [numspikes, 1]) + repmat(shifts, [1, numsamples]); 0092 0093 % Which regions are invalid? 0094 left_valid = max(find(any(resample_inds < 1))) + 1; if (isempty(left_valid)), left_valid = 1; end; 0095 right_valid = min(find(any(resample_inds > numsamples))) - 1; if (isempty(right_valid)), right_valid = numsamples; end; 0096 if (left_valid >= right_valid) 0097 error('SS:waveform_invalidated', 'Realignment invalidates all samples. Try reducing MAXSHIFT.'); 0098 end 0099 0100 % Although 'spline' can find independent interpolants for each row of a matrix 0101 % (warning: the help for 'spline' is misleading on rows vs. cols), it does not 0102 % allow each row to be resampled on its own grid. To get around this (since the 0103 % for-loop approach is ~4x slower), we obtain the piecewise polynomial coefficients 0104 % along each waveform, pad to handle endpoints, and then string them all together 0105 % into one long series. This will let us compute the shifted waveforms in one 0106 % fell swoop . . . the only caveat is that requests to extrapolate beyond 0107 % [1:numsamples] become meaningless; but we're not interested in these anyway. 0108 pp = spline([1:numsamples], spikes.waveforms); 0109 0110 % The coefficients from 'spline' come out ordered by column rather than by 0111 % row/waveform, so we reorder. 0112 pp.coefs = reshape(pp.coefs, numspikes, numsamples-1, []); 0113 pp.coefs = permute(pp.coefs, [2 1 3]); 0114 0115 % DEBUGGING -- see the splie 0116 % example = 1; 0117 % figure; cla; plot(spikes.waveforms(example,:), 'm.-', 'LineWidth', 4); hold on; 0118 % for t0 = 1:(numsamples-1) 0119 % t = linspace(0,1,10); v = polyval(squeeze(pp.coefs(t0,example,:)), t); plot(t+t0,v,'k'); 0120 % end 0121 % hold off; 0122 0123 % 'ppeval' uses the right spline at endpoints; i.e., if we just strung together 0124 % the splines as is, the last sample of each waveform would be interpolated 0125 % using the first polynomial from the next waveform. So we add a DC polynomial 0126 % to the end of each spike whose value equals to the last sample of the spike. 0127 % Now after stringing together, interpolation at ((r-1)*numsamples + [1:numsamples]) 0128 % gives the rth spike (within roundoff error). 0129 padzeros = zeros(1, numspikes, 4); 0130 pp.coefs = cat(1, pp.coefs, padzeros); 0131 pp.coefs(numsamples,:,4) = spikes.waveforms(:,end)'; 0132 0133 % Make the 'strung-together' version. 0134 pp.coefs = reshape(pp.coefs, [], 4); 0135 0136 % Rework the remaining information in the piecewise polynomial to be consistent. 0137 pp.pieces = numsamples * numspikes; pp.dim = 1; pp.breaks = [1:(pp.pieces+1)]; 0138 0139 % Change resample indices to correspond to the strung-together piecewise polynomial. 0140 offset = ([1:numspikes] - 1)' * numsamples; 0141 resample_inds = resample_inds + repmat(offset, [1, numsamples]); 0142 0143 % OK, finally. Do the resample & clip the invalid regions 0144 resampled = ppval(pp, resample_inds); 0145 spikes.waveforms = resampled(:, left_valid:right_valid); 0146 0147 % Also, the threshold crossing time has changed due to data clipping 0148 spikes.threshT = spikes.threshT - left_valid + 1; 0149 0150 % Worse, the threshold time/value may have changed while shifting. This 0151 % will make it hard to identify peaks later, so we need new values. 0152 meanlevel = mean(spikes.waveforms(:)); % estimate of baseline 0153 notallowed = [1:(spikes.threshT-maxshift-1),(spikes.threshT+maxshift+1):(size(spikes.waveforms, 2))]; 0154 0155 deflection = mean(abs(spikes.waveforms - meanlevel), 1); % mean deviation from baseline 0156 deflection(notallowed) = 0; 0157 [junk, spikes.threshT] = max(deflection); % max deflect near old threshT marks new threshold marker 0158 threshValues = spikes.waveforms(:, spikes.threshT); % values at new marker 0159 spikes.threshV(1) = max([threshValues(threshValues < meanlevel); -Inf]); % value closest to baseline but below it 0160 spikes.threshV(2) = min([threshValues(threshValues > meanlevel); Inf]); % value closest to baseline but above it 0161 0162 spikes.tictoc.dejitter = etime(clock, starttime);