Finds the peak locations in a signal
peakIdxs = find_peaks(inSig,params)
Accepts an input signal and identifies the peak index locations.
A threshold parameter (param.thresh) should be passed in to
indicate the threshold under which peaks should be ignored.
The threshold is useful if working with a signal that is
noisy. The threshold should be chosen sufficiently above the
noise floor. The signal may optionally be low pass filtered
before performing peak detection. A window parameter
(params.peak_height_window) may also be passed in, to indicate the size
of a sliding window of samples within which to only accept the highest
peak. For instance, for a peak window of 1 second, and given a signal
with many peaks that occur within a few tenths of a second of each other,
only the highest peak within a given 1 second window will be retained,
and all other peaks within that window will be discarded (equivalent to
taking the highest peak within +/- 0.5 s of any given peak).
INPUTS
params.perform: (optional) A cell array of strings of processes to perform.
'lowPass' - perform a low pass filtering before peak detection.
'peak_heights' - returns peak height values in second output argument
'peak_widths' - returns peak widths in third output argument
'peak_areas' - returns peak areas in fourth output argument
params.thresh: (optional) Value between 0 and 1. Fraction of maximum
amplitude to use as a threshold below which peak
detecting will not be performed. This may be useful if
your signal has many small peaks with low amplitudes.
params.Fs: Sampling rate of signal. Required if performing low pass
filtering, or if specifying a peak window.
params.lpFreq: Low pass cutoff frequency. Only needed if performing low pass
filtering. Requires params.Fs.
params.peak_height_thresh: (optional) peaks with individual heights less
than peak_height_thresh will not be retained. This
differs from params.thresh, in that params.thresh is a
global amplitude threshold, and peak_height_thresh is a
local peak height threshold.
params.peak_height_window: (optional) sliding window in which only the
highest peak is retained. for each peak found, if there
are other peaks within +/- peak_height_window seconds of
that peak, all but the highest peak will be discarded.
params.peak_heights, params.peak_widths, and params.peak_areas
are parameter substructures that are passed to peak_heights,
peak_widths, and peak_areas. If they are not specified, empty
parameter structures are used.
OUTPUTS
existance of particular outputs necessarily depend upon params.perform,
but the indices of outputs when they exist are such:
varargout{1} - peak indexes
varargout{2} - peak heights
varargout{3} - base indexes
varargout{4} - peak widths
varargout{5} - peak areas
Copyright (c) 2007-2012 The Regents of the University of California
All Rights Reserved
Author:
Stefan Tomic 5/07
06/17/09 - Stefan Tomic:
- added handling of xVals in find_peaks
- modified input argument formatting to accept
param/value pairs with compatibility to old function calls
- now returns peak heights, peak widths, and peak
areas as additional output arguments
06/18/09 - Fred Barrett
- added peak_height_window and peak_height_thresh
06/21/09 - Fred Barrett
- documnted outputs
- added "base indexes" to varargout, thus shifting peak widths
and peak areas (formerly varargout{3:4}) to varargout{4:5}
10/28/09 - Stefan Tomic:
- fixed handling of flat peaks (peaks that have consecutive
indices with the same Y values). find_peaks now returns
the middle index instead of the left-most index. In the
case of an even number of indices, it returns the left
of the two middle indices.
- Fixed a bug that was causing find_peaks to return a peak
index at the end of the signal when the end of the signal
was flat (e.g. silence at the end of an audio signal).0001 function varargout = find_peaks(varargin) 0002 % Finds the peak locations in a signal 0003 % 0004 % peakIdxs = find_peaks(inSig,params) 0005 % 0006 % Accepts an input signal and identifies the peak index locations. 0007 % A threshold parameter (param.thresh) should be passed in to 0008 % indicate the threshold under which peaks should be ignored. 0009 % The threshold is useful if working with a signal that is 0010 % noisy. The threshold should be chosen sufficiently above the 0011 % noise floor. The signal may optionally be low pass filtered 0012 % before performing peak detection. A window parameter 0013 % (params.peak_height_window) may also be passed in, to indicate the size 0014 % of a sliding window of samples within which to only accept the highest 0015 % peak. For instance, for a peak window of 1 second, and given a signal 0016 % with many peaks that occur within a few tenths of a second of each other, 0017 % only the highest peak within a given 1 second window will be retained, 0018 % and all other peaks within that window will be discarded (equivalent to 0019 % taking the highest peak within +/- 0.5 s of any given peak). 0020 % 0021 % INPUTS 0022 % params.perform: (optional) A cell array of strings of processes to perform. 0023 % 'lowPass' - perform a low pass filtering before peak detection. 0024 % 'peak_heights' - returns peak height values in second output argument 0025 % 'peak_widths' - returns peak widths in third output argument 0026 % 'peak_areas' - returns peak areas in fourth output argument 0027 % params.thresh: (optional) Value between 0 and 1. Fraction of maximum 0028 % amplitude to use as a threshold below which peak 0029 % detecting will not be performed. This may be useful if 0030 % your signal has many small peaks with low amplitudes. 0031 % 0032 % params.Fs: Sampling rate of signal. Required if performing low pass 0033 % filtering, or if specifying a peak window. 0034 % params.lpFreq: Low pass cutoff frequency. Only needed if performing low pass 0035 % filtering. Requires params.Fs. 0036 % 0037 % params.peak_height_thresh: (optional) peaks with individual heights less 0038 % than peak_height_thresh will not be retained. This 0039 % differs from params.thresh, in that params.thresh is a 0040 % global amplitude threshold, and peak_height_thresh is a 0041 % local peak height threshold. 0042 % 0043 % params.peak_height_window: (optional) sliding window in which only the 0044 % highest peak is retained. for each peak found, if there 0045 % are other peaks within +/- peak_height_window seconds of 0046 % that peak, all but the highest peak will be discarded. 0047 % 0048 % params.peak_heights, params.peak_widths, and params.peak_areas 0049 % are parameter substructures that are passed to peak_heights, 0050 % peak_widths, and peak_areas. If they are not specified, empty 0051 % parameter structures are used. 0052 % 0053 % OUTPUTS 0054 % existance of particular outputs necessarily depend upon params.perform, 0055 % but the indices of outputs when they exist are such: 0056 % varargout{1} - peak indexes 0057 % varargout{2} - peak heights 0058 % varargout{3} - base indexes 0059 % varargout{4} - peak widths 0060 % varargout{5} - peak areas 0061 % 0062 % Copyright (c) 2007-2012 The Regents of the University of California 0063 % All Rights Reserved 0064 % 0065 % Author: 0066 % Stefan Tomic 5/07 0067 % 06/17/09 - Stefan Tomic: 0068 % - added handling of xVals in find_peaks 0069 % - modified input argument formatting to accept 0070 % param/value pairs with compatibility to old function calls 0071 % - now returns peak heights, peak widths, and peak 0072 % areas as additional output arguments 0073 % 06/18/09 - Fred Barrett 0074 % - added peak_height_window and peak_height_thresh 0075 % 06/21/09 - Fred Barrett 0076 % - documnted outputs 0077 % - added "base indexes" to varargout, thus shifting peak widths 0078 % and peak areas (formerly varargout{3:4}) to varargout{4:5} 0079 % 10/28/09 - Stefan Tomic: 0080 % - fixed handling of flat peaks (peaks that have consecutive 0081 % indices with the same Y values). find_peaks now returns 0082 % the middle index instead of the left-most index. In the 0083 % case of an even number of indices, it returns the left 0084 % of the two middle indices. 0085 % - Fixed a bug that was causing find_peaks to return a peak 0086 % index at the end of the signal when the end of the signal 0087 % was flat (e.g. silence at the end of an audio signal). 0088 0089 %this section handles compatibility with the old way of calling 0090 %find_peaks, which just had two input arguments (inData,params) 0091 if(~ischar(varargin{1})) 0092 inSig = varargin{1}; 0093 if(nargin > 1) 0094 params = varargin{2}; 0095 end 0096 0097 %set varargin to new input argument formatting, so that we can 0098 %forward the arguments to peak_heights, peak_widths, and peak_areas 0099 varargin = {'signal',inSig,'params',params}; 0100 0101 else 0102 0103 %this is the new way of calling find_peaks, with tag/value pairs. 0104 %the main reason for going this direction is that peak_widths, 0105 %peak_heights and peak_areas use tag value pairs and it is easy to just forward 0106 %the arguments. 0107 for iarg = 1:2:nargin 0108 switch varargin{iarg} 0109 case 'signal' 0110 inSig = varargin{iarg+1}; 0111 case 'params' 0112 params = varargin{iarg+1}; 0113 case 'xVals' 0114 xVals = varargin{iarg+1}; 0115 case 'perform' 0116 perform = varargin{iarg+1}; 0117 end 0118 end 0119 end 0120 0121 try 0122 perform; 0123 catch 0124 perform = {'peak_heights','peak_widths','peak_areas'}; 0125 end 0126 0127 try 0128 xVals; 0129 catch 0130 %find_peaks didn't receive xVals, so peak indices will be returned 0131 xVals = [1:length(inSig)]'; 0132 end 0133 0134 %make sure that we are dealing with a column vector 0135 if(size(inSig,1) > 1 && size(inSig,2) > 1) 0136 error('Input signal must be a vector'); 0137 elseif(size(inSig,1) > 1 && size(inSig,2) == 1) 0138 inSig = inSig; 0139 else 0140 inSig = inSig'; 0141 end 0142 0143 if(~isfield(params,'thresh')) 0144 params.thresh = 0; 0145 end 0146 0147 if(isfield(params,'perform') && ismember('lowPass',params.perform)) 0148 Fs = params.Fs; 0149 lpFreq = params.lpFreq; 0150 W = lpFreq/Fs*2; 0151 [B,A] = butter(12,W,'low'); 0152 inSig = filtfilt(B,A,inSig); 0153 end 0154 0155 %the 0 at the beginning takes care of the diff offset and 0156 %solves for possible peaks at the very beginning of the signal. The 0157 %-1 value at the end takes care of possible peaks at the end of the 0158 %signal. 0159 diffSig = diff(inSig,1); 0160 sigIsFlat = (diffSig == 0); 0161 sigNotFlatIdxs = [1 ; find(~sigIsFlat) + 1]; 0162 sigIsFlatIdxs = find(sigIsFlat) + 1; 0163 sigNotFlat = inSig(sigNotFlatIdxs); 0164 diffSigNotFlat = diff(sigNotFlat,1); 0165 0166 peakRefMask = [0 ; (diffSigNotFlat(1:end-1) > 0 & diffSigNotFlat(2:end) < 0)]; 0167 0168 peakRefs = find(peakRefMask); 0169 0170 leftPeakIdxs = sigNotFlatIdxs(peakRefs); 0171 0172 %deal with edge conditions 0173 if(diffSig(1) < 0) 0174 leftPeakIdxs = [1 ; leftPeakIdxs]; 0175 end 0176 0177 if(diffSig(end) > 0) 0178 leftPeakIdxs = [leftPeakIdxs ; length(inSig)]; 0179 end 0180 0181 possiblePeakMask = zeros(size(inSig)); 0182 %check if any detected possible peaks are on flat locations 0183 %detected points will be on left side of flat regions, so 0184 %we need to check flatness of right side 0185 nLeftPeaks = length(leftPeakIdxs); 0186 for iLeftPeak = 1:nLeftPeaks 0187 thisPeakIdx = leftPeakIdxs(iLeftPeak); 0188 %if the next not-flat index is not adjacent to detected 0189 %then we are on a flat location 0190 nextIdx = sigNotFlatIdxs(find(sigNotFlatIdxs > thisPeakIdx,1,'first')); 0191 if(isempty(nextIdx)) 0192 thisPeakAdjIdxs = thisPeakIdx; 0193 else 0194 thisPeakAdjIdxs = [thisPeakIdx:nextIdx-1]; 0195 end 0196 0197 minAdjIdx = min(thisPeakAdjIdxs); 0198 maxAdjIdx = max(thisPeakAdjIdxs); 0199 0200 middleIdx = floor((maxAdjIdx - minAdjIdx)/2) + minAdjIdx; 0201 0202 possiblePeakMask(middleIdx) = 1; 0203 end 0204 0205 %only return portions of the signal that are above the threshold. 0206 threshVal = min(inSig) + (max(inSig) - min(inSig))*params.thresh; 0207 peakIdxs = find((inSig >= threshVal) & possiblePeakMask); 0208 0209 try params.peak_heights; catch params.peak_heights = []; end 0210 [pheights,bidxs] = peak_heights('signal',inSig,'params',... 0211 params.peak_heights,'peakIdxs',peakIdxs,'xVals',xVals); 0212 0213 if isfield(params,'peak_height_thresh') 0214 rmidxs = pheights < params.peak_height_thresh; 0215 peakIdxs(rmidxs) = []; 0216 pheights(rmidxs) = []; 0217 bidxs(rmidxs) = []; 0218 end 0219 0220 % remove peaks within +/- peak_window seconds of each peak 0221 if isfield(params,'peak_height_window') 0222 if isfield(params,'Fs') 0223 for ip=length(peakIdxs):-1:1 0224 0225 if ip > length(peakIdxs), continue, end 0226 % get peak window 0227 win_start=peakIdxs(ip)-round(params.Fs*params.peak_height_window); 0228 win_end =peakIdxs(ip)+round(params.Fs*params.peak_height_window); 0229 if win_start < 0, win_start = 0; end 0230 if win_end > length(inSig), win_end = length(inSig); end 0231 window = (win_start:win_end); 0232 0233 % get peaks within window 0234 widxs = find(ismember(peakIdxs,window)); 0235 wheights = pheights(widxs); 0236 0237 % don't keep the max height widxs 0238 maxidx = find(wheights==max(wheights)); 0239 if length(maxidx) > 1, maxidx = maxidx(1); end 0240 widxs(maxidx) = []; 0241 0242 % remove widxs from peakIdxs 0243 if ~isempty(widxs) && ~isempty(peakIdxs) 0244 peakIdxs(widxs) = []; 0245 pheights(widxs) = []; 0246 bidxs(widxs) = []; 0247 end 0248 end 0249 else 0250 warning('peak_height_window not applied: no sampling rate specified'); 0251 end 0252 end 0253 0254 varargout{1} = xVals(peakIdxs); 0255 0256 if(ismember('peak_heights',perform)) 0257 varargout{2} = pheights; 0258 varargout{3} = bidxs; 0259 end 0260 0261 if(ismember('peak_widths',perform)) 0262 try params.peak_widths; catch params.peak_widths = []; end; 0263 pwidths = peak_widths('signal',inSig,'params',params.peak_widths,'peakIdxs',peakIdxs,'xVals',xVals); 0264 %peak widths returns the borders, so calculate actual widths here 0265 varargout{4} = pwidths(2:2:end) - pwidths(1:2:end); 0266 end 0267 0268 if(ismember('peak_areas',perform)) 0269 try params.peak_areas; catch params.peak_areas = []; end; 0270 varargout{5} = peak_areas('signal',inSig,'params',params.peak_areas,'peakIdxs',peakIdxs,'xVals',xVals); 0271 end 0272 0273 return 0274 0275