Home > chronux_1_50 > 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 (~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 numspikes = size(spikes.waveforms, 1);
0062 numsamples = size(spikes.waveforms, 2);
0063 
0064 %%%%%%%%%% FIND FIDUCIALS
0065 % first, get a mask over the peaks
0066 isThreshLo = (spikes.waveforms(:, spikes.threshT) < spikes.threshV(1));
0067 [h, w, pl, mask] = thresholded_peaks(spikes);
0068 
0069 % next, subtract appropriate thresholds, making all values in the peak
0070 % have positive sign, since we're about to treat them as 'mass'
0071 waves = spikes.waveforms;
0072 waves(~isThreshLo, :) = waves(~isThreshLo, :) - spikes.threshV(2);
0073 waves( isThreshLo, :) = spikes.threshV(1) - waves( isThreshLo, :);
0074 waves = waves .* mask;
0075 
0076 % now COM is straightforward
0077 fiducials = (waves * (1:numsamples)') ./ sum(waves, 2);
0078 
0079 clear mask waves;  % manual garbage collection
0080 
0081 %%%%%%%%%% REALIGN FIDUCIALS
0082 % We line up the fiducials around the sample that requires the least shifting.
0083 target = round(mean(fiducials));
0084 
0085 % Determine shifted indices for each waveform
0086 shifts = fiducials - target;
0087 shifts(abs(shifts) > maxshift) = 0;   % big shifts are outliers; don't alter these
0088 resample_inds = repmat((1:numsamples), [numspikes, 1]) + repmat(shifts, [1, numsamples]);
0089 
0090 % Which regions are invalid?
0091 left_valid = max(find(any(resample_inds < 1))) + 1;            if (isempty(left_valid)), left_valid = 1;  end;
0092 right_valid = min(find(any(resample_inds > numsamples))) - 1;  if (isempty(right_valid)), right_valid = numsamples; end;
0093 if (left_valid >= right_valid)
0094     error('SS:waveform_invalidated', 'Realignment invalidates all samples.  Try reducing MAXSHIFT.');
0095 end
0096 
0097 % Although 'spline' can find independent interpolants for each row of a matrix
0098 % (warning: the help for 'spline' is misleading on rows vs. cols), it does not
0099 % allow each row to be resampled on its own grid.  To get around this (since the
0100 % for-loop approach is ~4x slower), we obtain the piecewise polynomial coefficients
0101 % along each waveform, pad to handle endpoints, and then string them all together
0102 % into one long series.  This will let us compute the shifted waveforms in one
0103 % fell swoop . . . the only caveat is that requests to extrapolate beyond
0104 % [1:numsamples] become meaningless; but we're not interested in these anyway.
0105 pp = spline((1:numsamples), spikes.waveforms);
0106 
0107 % The coefficients from 'spline' come out ordered by column rather than by
0108 % row/waveform, so we reorder.
0109 pp.coefs = reshape(pp.coefs, numspikes, numsamples-1, []);
0110 pp.coefs = permute(pp.coefs, [2 1 3]);
0111 
0112 % DEBUGGING -- see the splie
0113 % example = 1;
0114 % figure;  cla;  plot(spikes.waveforms(example,:), 'm.-', 'LineWidth', 4);  hold on;
0115 % for t0 = 1:(numsamples-1)
0116 %     t = linspace(0,1,10);   v = polyval(squeeze(pp.coefs(t0,example,:)), t);  plot(t+t0,v,'k');
0117 % end
0118 % hold off;
0119 
0120 % 'ppeval' uses the right spline at endpoints; i.e., if we just strung together
0121 % the splines as is, the last sample of each waveform would be interpolated
0122 % using the first polynomial from the next waveform.  So we add a DC polynomial
0123 % to the end of each spike whose value equals to the last sample of the spike.
0124 % Now after stringing together, interpolation at ((r-1)*numsamples + [1:numsamples])
0125 % gives the rth spike (within roundoff error).
0126 padzeros = zeros(1, numspikes, 4);
0127 pp.coefs = cat(1, pp.coefs, padzeros);
0128 pp.coefs(numsamples,:,4) = spikes.waveforms(:,end)';
0129 
0130 % Make the 'strung-together' version.
0131 pp.coefs = reshape(pp.coefs, [], 4);
0132 
0133 % Rework the remaining information in the piecewise polynomial to be consistent.
0134 pp.pieces = numsamples * numspikes;   pp.dim = 1;   pp.breaks = (1:(pp.pieces+1));
0135 
0136 % Change resample indices to correspond to the strung-together piecewise polynomial.
0137 offset = ((1:numspikes) - 1)' * numsamples;
0138 resample_inds = resample_inds + repmat(offset, [1, numsamples]);
0139 
0140 % OK, finally.  Do the resample & clip the invalid regions
0141 resampled = ppval(pp, resample_inds);
0142 spikes.waveforms = resampled(:, left_valid:right_valid);
0143 
0144 % Also, the threshold crossing time has changed due to data clipping
0145 spikes.threshT = spikes.threshT - left_valid + 1;
0146 
0147 % Worse, the threshold time/value may have changed while shifting.  This
0148 % will make it hard to identify peaks later, so we need new values.
0149 meanlevel = mean(spikes.waveforms(:));   % estimate of baseline
0150 notallowed = [1:(spikes.threshT-maxshift-1),(spikes.threshT+maxshift+1):(size(spikes.waveforms, 2))];
0151 
0152 deflection = mean(abs(spikes.waveforms - meanlevel), 1);  % mean deviation from baseline
0153 deflection(notallowed) = 0;
0154 [junk, spikes.threshT] = max(deflection);  % max deflect near old threshT marks new threshold marker
0155 threshValues = spikes.waveforms(:, spikes.threshT);  % values at new marker
0156 spikes.threshV(1) = max([threshValues(threshValues < meanlevel); -Inf]); % value closest to baseline but below it
0157 spikes.threshV(2) = min([threshValues(threshValues > meanlevel); Inf]); % value closest to baseline but above it
0158 
0159 spikes.tictoc.dejitter = etime(clock, starttime);

Generated on Mon 09-Oct-2006 00:54:52 by m2html © 2003