Home > chronux_1_1 > spikesort > helper > thresholded_peaks.m

thresholded_peaks

PURPOSE ^

THRESHOLDED_PEAKS Finds height/width of threshold crossing peaks.

SYNOPSIS ^

function [heights, widths, peak_locs, mask] = thresholded_peaks(spikes)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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!

Generated on Sun 13-Aug-2006 11:49:44 by m2html © 2003