function [PS,diff]=peakstats(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype,displayit) % function % PS=peakstats(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup, % smoothtype,displayit). Function to locate the positive peaks in a noisy % x-y time series data set and compute a table of summary statistics of the % peak intervals (the x-axis interval between adjacent detected peaks), % heights, widths, and areas, listing the maximum, minimum, average, and % percent standard deviation, median, and mode, and optionally displaying % the x,t data plot with numbered peaks and the histograms of the peak % intervals, heights, widths, and areas in figure window 2. Detects peaks % by looking for downward zero-crossings in the first derivative that % exceed SlopeThreshold. Returns peak statistics (PS) containing of the % maximumn, minimum, mean, median, mode, and percent standard deviation of % the peak intervals, heights, widths, in a matrix. Optionally returns the % peak intervals as a vector 'diff'. % Input 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 optional 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) % The optional last argument "displayit" =1 if the x,y data plot, numbered % peaks, and peak statistics table and the histogtams are to be displayed, % othersize not. % See http://terpconnect.umd.edu/~toh/spectrum/Smoothing.html and % http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm % T. C. O'Haver, 2013, 2016. Version 2. Last revised March, 2016 % % Examples: % x=[0:.1:100];y=5+5.*cos(x)+randn(size(x));peakstats(x,y,0,-1,11,19,3,1); % As above, but no display, and statistics table returned in matrix PS: % x=[0:.01:5]';PS=peakstats(x,x.*sin(x.^2).^2,0,-1,5,5); % x=[0:.1:100];y=5+5.*cos(x)+randn(size(x));PS=peakstats(x,y,0,-1,11,19,3); % x=[0:.01:5]';peakstats(x',x.*sin(x.^2).^2,0,-1,5,5,1,1) display % statistics table and histogram plots % % Related functions: % findvalleys.m, findpeaksL.m, findpeaksb.m, findpeaksplot.m, findpeaks.m, % findpeaksnr.m, findpeaksGSS.m, findpeaksLSS.m, findpeaksfit.m. if nargin==6;smoothtype=1;displayit=0;end % smoothtype=1, displayit=0 if not specified in 7th argument if nargin==7;displayit=0;end % displayit=0 if not specified in 8th 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),smoothwidth,smoothtype); 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, % 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, [Height, Position, Width]=gaussfit(xx,yy); PeakX=real(Position); % Compute peak position and height of fitted parabola PeakY=real(Height); MeasuredWidth=real(Width); % if the peak is too narrow for least-squares technique to work % well, just use the max value of y in the sub-group of points near peak. else PeakY=max(yy); pindex=val2ind(yy,PeakY); PeakX=xx(pindex(1)); MeasuredWidth=0; end % Construct matrix P. One row for each peak % detected, containing the peak number, peak % position (x-value) and peak height (y-value). % If peak measurements fails and results in NaN, skip this % peak if isnan(PeakX) || isnan(PeakY) || PeakYsa(2), sd=std(a); else sd=(std(a')); end; function isOctave = IsOctave() % Returns true if this code is being executed by Octave. % Returns false if this code is being executed by MATLAB, or any other MATLAB % variant. % % usage: isOctave = IsOctave() persistent octaveVersionIsBuiltIn; if (isempty(octaveVersionIsBuiltIn)) octaveVersionIsBuiltIn = (exist('OCTAVE_VERSION', 'builtin') == 5); % exist returns 5 to indicate a built-in function. end isOctave = octaveVersionIsBuiltIn; % If OCTAVE_VERSION is a built-in function, then we must be in Octave. % Since the result cannot change between function calls, it is cached in a % persistent variable. isOctave cannot be a persistent variable, because it % is the return value of the function, so instead the persistent result must % be cached in a separate variable.