寻峰算法资料

function P=findpeaks（xySlopeThresholdAmpThresholdsmoothwidthpeakgroupsmoothtype）
% function P=findpeaks（xySlopeThresholdAmpThresholdsmoothwidthpeakgroupsmoothtype）
% Function to locate the positive peaks in a noisy x-y time series data
% set.  Detects peaks by looking for downward zero-crossings
% in the first derivative that exceed SlopeThreshold.
% Returns list （P） containing peak number and position 
% height width and area of each peak. Arguments “slopeThreshold“
% “ampThreshold“ and “smoothwidth“ control peak sensitivity.
% Higher values will neglect smaller features. “Smoothwidth“ is
% the width of the smooth applied before peak detection; larger
% values ignore narrow peaks. If smoothwidth=0 no smoothing
% is performed. “Peakgroup“ is the number points around the top 
% part of the peak that are taken for measurement. If Peakgroup=0
% the local maximum is takes as the peak height and position.
% The argument “smoothtype“ determines the smooth algorithm:
%   If smoothtype=1 rectangular （sliding-average or boxcar） 
%   If smoothtype=2 triangular （2 passes of sliding-average）
%   If smoothtype=3 pseudo-Gaussian （3 passes of sliding-average）
% See http://terpconnect.umd.edu/~toh/spectrum/Smoothing.html and 
% http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm
% T. C. O‘Haver 1995.  Version 5.1 Last revised December 2012
% Skip peaks if peak measurement results in NaN values
% Examples:
% findpeaks（0:.01:2humps（0:.01:2）0-155）
% x=[0:.01:50];findpeaks（xcos（x）0-155）
% x=[0:.01:5]‘;findpeaks（xx.*sin（x.^2）.^20-155）
if nargin~=7;smoothtype=1;end  % smoothtype=1 if not specified in argument
if smoothtype>3;smoothtype=3;end
if smoothtype<1;smoothtype=1;end 
smoothwidth=round（smoothwidth）;
peakgroup=round（peakgroup）;
if smoothwidth>1
    d=fastsmooth（deriv（y）smoothwidthsmoothtype）;
else
    d=y;
end
n=round（peakgroup/2+1）;
P=[0 0 0 0 0];
vectorlength=length（y）;
peak=1;
AmpTest=AmpThreshold;
for j=2*round（smoothwidth/2）-1:length（y）-smoothwidth
    if sign（d（j）） > sign （d（j+1）） % Detects zero-crossing
        if d（j）-d（j+1） > SlopeThreshold*y（j） % if slope of derivative is larger than SlopeThreshold
            if y（j） > AmpTest  % if height of peak is larger than AmpThreshold
                xx=zeros（size（peakgroup））;yy=zeros（size（peakgroup））;
                for k=1:peakgroup % Create sub-group of points near peak
                    groupindex=j+k-n+2;
                    if groupindex<1 groupindex=1;end
                    if groupindex>vectorlength groupindex=vectorlength;end
                    xx（k）=x（groupindex）;yy（k）=y（groupindex）;
                end
                if peakgroup>3
                    [coefSMU]=polyfit（xxlog（abs（yy））2）;  % Fit parabola to log10 of sub-group with centering and scaling
                    c1=coef（3）;c2=coef（2）;c3=coef（1）;
                    PeakX=-（（MU（2）.*c2/（2*c3））-MU（1））;   % Compute