


THRESHOLDED_PEAKS Finds height/width of threshold crossing peaks.
[HEIGHTS, WIDTHS] = THRESHOLDED_PEAKS(SPIKES) takes a spike-sorting
SS object SPIKES with M spikes and returns two (M x 1) vectors
HEIGHTS and WIDTHS, containing measurements for each waveform in
the SPIKES object.
HEIGHTS and WIDTHS are determined from the central peak and
require both spikes.threshT and spikes.threshV(1:2) to be defined
(waveforms from low and high thresholds are handled automatically).
The central peak is taken as the contiguous region starting at the
spikes.threshT sample and ending when threshold (either low or
high) is recrossed in the opposite direction. Peak height
corresponds to the waveform value in this region that is farthest
from threshold and peak width is length of the region in samples.
[HEIGHTS, WIDTHS, PEAK_LOCS] = THRESHOLDED_PEAKS(SPIKES) returns
an (M x 3) matrix PEAK_LOCS in which the second column gives the
location (in samples) in the waveform that attained the HEIGHT,
while the first and third columns give the threshold cross and
recross locations (in samples) used to compute the WIDTHS.
[HEIGHTS, WIDTHS, PEAK_LOCS, MASK] = THRESHOLDED_PEAKS(SPIKES) further
returns an (M x N) matrix MASK, where (M x N) are the dimensions of
the spikes.waveforms matrix. MASK is a binary matrix with a
contiguous block of 1's in each row marking the samples contained in
the thresholded central peak for that row.

0001 function [heights, widths, peak_locs, mask] = thresholded_peaks(spikes) 0002 0003 % THRESHOLDED_PEAKS Finds height/width of threshold crossing peaks. 0004 % [HEIGHTS, WIDTHS] = THRESHOLDED_PEAKS(SPIKES) takes a spike-sorting 0005 % SS object SPIKES with M spikes and returns two (M x 1) vectors 0006 % HEIGHTS and WIDTHS, containing measurements for each waveform in 0007 % the SPIKES object. 0008 % 0009 % HEIGHTS and WIDTHS are determined from the central peak and 0010 % require both spikes.threshT and spikes.threshV(1:2) to be defined 0011 % (waveforms from low and high thresholds are handled automatically). 0012 % The central peak is taken as the contiguous region starting at the 0013 % spikes.threshT sample and ending when threshold (either low or 0014 % high) is recrossed in the opposite direction. Peak height 0015 % corresponds to the waveform value in this region that is farthest 0016 % from threshold and peak width is length of the region in samples. 0017 % 0018 % [HEIGHTS, WIDTHS, PEAK_LOCS] = THRESHOLDED_PEAKS(SPIKES) returns 0019 % an (M x 3) matrix PEAK_LOCS in which the second column gives the 0020 % location (in samples) in the waveform that attained the HEIGHT, 0021 % while the first and third columns give the threshold cross and 0022 % recross locations (in samples) used to compute the WIDTHS. 0023 % 0024 % [HEIGHTS, WIDTHS, PEAK_LOCS, MASK] = THRESHOLDED_PEAKS(SPIKES) further 0025 % returns an (M x N) matrix MASK, where (M x N) are the dimensions of 0026 % the spikes.waveforms matrix. MASK is a binary matrix with a 0027 % contiguous block of 1's in each row marking the samples contained in 0028 % the thresholded central peak for that row. 0029 0030 % Last Modified By: sbm on Thu Oct 6 20:30:16 2005 0031 0032 %%%%%%%%%% ARGUMENT CHECKING 0033 if (~isfield(spikes, 'threshT')) 0034 error('SS:threshT_Undefined', 'The SS object must define the sample index of threshold crossing.'); 0035 elseif (~isfield(spikes, 'threshV') | (length(spikes.threshV) ~= 2)) 0036 error('SS:threshV_Error', 'The SS object must define both low and high thresholds. Use +/- Inf if only one of the two was used.'); 0037 elseif ((spikes.threshV(2) - spikes.threshV(1)) < 0) 0038 error('SS:threshV_Illegal', 'The SS object high threshold must be greater than its low threshold.'); 0039 end 0040 0041 %%%%%%%%%% FIND PEAK REGIONS 0042 % This could be done with some straightforward code in a big, ugly 'for' 0043 % loop, but, hey, what's Matlab for if not replacing FOR loops with big, 0044 % ugly, vectorized code? At least its faster this way. 0045 0046 % Mark each time a waveform crosses threshold by taking the derivative 0047 % of the binary thresholded spike data. Zero-padding ensures that the 0048 % first and last sample points can be marked, with the result that 0049 % the # of times the signal crosses threshold (marked with +1) equals 0050 % the number of times it returns below it (marked with -1). 0051 lopeaks = diff(padmatrix([spikes.waveforms <= spikes.threshV(1)], [1 1 0 0]), 1, 2); 0052 hipeaks = diff(padmatrix([spikes.waveforms >= spikes.threshV(2)], [1 1 0 0]), 1, 2); 0053 0054 % The lo/hi markers can be combined; we'll remember which is which. 0055 isThreshLo = (spikes.waveforms(:, spikes.threshT) <= spikes.threshV(1)); 0056 allpeaks = zeros(size(lopeaks)); 0057 allpeaks( isThreshLo,:) = lopeaks( isThreshLo,:); 0058 allpeaks(~isThreshLo,:) = hipeaks(~isThreshLo,:); 0059 clear lopeaks hipeaks; % manual garbage collection . . . 0060 0061 % Read the marks into a list of peak start/stop indices. Note that the 0062 % +1's line up properly but the -1's correspond to the first sample _after_ 0063 % threshold recrossing; so we subtract 1 to get start/stop indices for the 0064 % contiguous block of supra-threshold samples. 0065 [start_r start_c] = find(allpeaks == 1); % pts that first exceed threshold 0066 [finis_r finis_c] = find(allpeaks == -1); % pts that just sank below 0067 start = sortrows([start_r start_c]); % sort both starts and ... 0068 finis = sortrows([finis_r finis_c]); % ... ends, so we can ... 0069 peak_inds = [start (finis(:,2)-1)]; % ... combine them (since they're 1:1). 0070 clear allpeaks; 0071 0072 % Next, select out the peaks that contain spikes.threshT. 0073 peaks_center = find((peak_inds(:,2) <= spikes.threshT) & (peak_inds(:,3) >= spikes.threshT)); 0074 if (length(peaks_center) ~= size(spikes.waveforms, 1)) 0075 error('SS:thresh_incorrect', ['The threshold information is invalid; some waveforms ' ... 0076 'do not have a central peak with these parameters.']); 0077 end 0078 0079 % Make the central peak mask; we need this to figure out the heights. 0080 mask = zeros(size(spikes.waveforms)); 0081 for k = 1:length(peaks_center) 0082 peakrow = peak_inds(peaks_center(k),:); 0083 mask(peakrow(1), peakrow(2):peakrow(3)) = 1; 0084 end 0085 0086 % OK, now heights are fairly easy. Just need to keep lo/high thresholds straight. 0087 peaks = mask .* spikes.waveforms; 0088 heights = zeros(size(spikes.waveforms, 1), 1); 0089 [loHeights, loWhere]= min(peaks, [], 2); 0090 [hiHeights, hiWhere] = max(peaks, [], 2); 0091 heights( isThreshLo) = loHeights( isThreshLo); 0092 heights(~isThreshLo) = hiHeights(~isThreshLo); 0093 0094 % And widths are even easier. 0095 widths = peak_inds(peaks_center, 3) - peak_inds(peaks_center, 2); 0096 0097 % Finally, we build the peak_locs by reorganizing information we already have. 0098 if (nargout > 2) 0099 peak_locs = zeros(size(spikes.waveforms, 1), 3); 0100 peak_locs(:, [1,3]) = peak_inds(peaks_center, [2,3]); % width info 0101 peak_locs( isThreshLo,2) = loWhere( isThreshLo); % height info 0102 peak_locs(~isThreshLo,2) = hiWhere(~isThreshLo); 0103 end 0104 0105 % Thats all, folks!