Home > chronux_2_00 > 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 %%%%%%%%%% 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!

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