Home > chronux_2_00 > spikesort > ss_dejitter.m

ss_dejitter

PURPOSE ^

SS_DEJITTER Aligns waveform peaks.

SYNOPSIS ^

function spikes = ss_dejitter(spikes, maxshift)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 (spikes.threshT > size(spikes.waveforms, 2))
0054     error('SS:threshT_invalid', 'The threshold time ThreshT does not fall within the supplied data waveform.');
0055 elseif (~isfield(spikes, 'threshV') || (length(spikes.threshV) ~= 2))
0056     error('SS:threshV_invalid', 'The SS object must define both low and high thresholds.  Use +/- Inf if only one of the two was used.'); 
0057 elseif ((spikes.threshV(2) - spikes.threshV(1)) < 0)
0058     error('SS:threshV_illegal', 'The SS object high threshold must be greater than its low threshold.');
0059 end
0060 
0061 
0062 %%%%%%%%%% CONSTANTS
0063 % spline_chunk = 5000;   % # of splines to do at a time; keeps memory use down
0064 numspikes = size(spikes.waveforms, 1);
0065 numsamples = size(spikes.waveforms, 2);
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);

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003