


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 %%%%%%%%%% ARGUMENT CHECKING 0031 if (~isfield(spikes, 'threshT')) 0032 error('SS:threshT_Undefined', 'The SS object must define the sample index of threshold crossing.'); 0033 elseif (~isfield(spikes, 'threshV') || (length(spikes.threshV) ~= 2)) 0034 error('SS:threshV_Error', 'The SS object must define both low and high thresholds. Use +/- Inf if only one of the two was used.'); 0035 elseif ((spikes.threshV(2) - spikes.threshV(1)) < 0) 0036 error('SS:threshV_Illegal', 'The SS object high threshold must be greater than its low threshold.'); 0037 end 0038 0039 %%%%%%%%%% FIND PEAK REGIONS 0040 % This could be done with some straightforward code in a big, ugly 'for' 0041 % loop, but, hey, what's Matlab for if not replacing FOR loops with big, 0042 % ugly, vectorized code? At least its faster this way. 0043 0044 % Mark each time a waveform crosses threshold by taking the derivative 0045 % of the binary thresholded spike data. Zero-padding ensures that the 0046 % first and last sample points can be marked, with the result that 0047 % the # of times the signal crosses threshold (marked with +1) equals 0048 % the number of times it returns below it (marked with -1). 0049 lopeaks = diff(padmatrix((spikes.waveforms <= spikes.threshV(1)), [1 1 0 0]), 1, 2); 0050 hipeaks = diff(padmatrix((spikes.waveforms >= spikes.threshV(2)), [1 1 0 0]), 1, 2); 0051 0052 % The lo/hi markers can be combined; we'll remember which is which. 0053 isThreshLo = (spikes.waveforms(:, spikes.threshT) <= spikes.threshV(1)); 0054 allpeaks = zeros(size(lopeaks)); 0055 allpeaks( isThreshLo,:) = lopeaks( isThreshLo,:); 0056 allpeaks(~isThreshLo,:) = hipeaks(~isThreshLo,:); 0057 clear lopeaks hipeaks; % manual garbage collection . . . 0058 0059 % Read the marks into a list of peak start/stop indices. Note that the 0060 % +1's line up properly but the -1's correspond to the first sample _after_ 0061 % threshold recrossing; so we subtract 1 to get start/stop indices for the 0062 % contiguous block of supra-threshold samples. 0063 [start_r start_c] = find(allpeaks == 1); % pts that first exceed threshold 0064 [finis_r finis_c] = find(allpeaks == -1); % pts that just sank below 0065 start = sortrows([start_r start_c]); % sort both starts and ... 0066 finis = sortrows([finis_r finis_c]); % ... ends, so we can ... 0067 peak_inds = [start (finis(:,2)-1)]; % ... combine them (since they're 1:1). 0068 clear allpeaks; 0069 0070 % Next, select out the peaks that contain spikes.threshT. 0071 peaks_center = find((peak_inds(:,2) <= spikes.threshT) & (peak_inds(:,3) >= spikes.threshT)); 0072 if (length(peaks_center) ~= size(spikes.waveforms, 1)) 0073 error('SS:thresh_incorrect', ['The threshold information is invalid; some waveforms ' ... 0074 'do not have a central peak with these parameters.']); 0075 end 0076 0077 % Make the central peak mask; we need this to figure out the heights. 0078 mask = zeros(size(spikes.waveforms)); 0079 for k = 1:length(peaks_center) 0080 peakrow = peak_inds(peaks_center(k),:); 0081 mask(peakrow(1), peakrow(2):peakrow(3)) = 1; 0082 end 0083 0084 % OK, now heights are fairly easy. Just need to keep lo/high thresholds straight. 0085 peaks = mask .* spikes.waveforms; 0086 heights = zeros(size(spikes.waveforms, 1), 1); 0087 [loHeights, loWhere]= min(peaks, [], 2); 0088 [hiHeights, hiWhere] = max(peaks, [], 2); 0089 heights( isThreshLo) = loHeights( isThreshLo); 0090 heights(~isThreshLo) = hiHeights(~isThreshLo); 0091 0092 % And widths are even easier. 0093 widths = peak_inds(peaks_center, 3) - peak_inds(peaks_center, 2); 0094 0095 % Finally, we build the peak_locs by reorganizing information we already have. 0096 if (nargout > 2) 0097 peak_locs = zeros(size(spikes.waveforms, 1), 3); 0098 peak_locs(:, [1,3]) = peak_inds(peaks_center, [2,3]); % width info 0099 peak_locs( isThreshLo,2) = loWhere( isThreshLo); % height info 0100 peak_locs(~isThreshLo,2) = hiWhere(~isThreshLo); 0101 end 0102 0103 % Thats all, folks!