%% Peak Detection Tool
% with a selection of different peak detectors and optional sharpening and least-squares curve fitting (which requires peakfit.m in the path). To see the graphs on the right, right-click on the right panel and select "Disable synchronous scrolling". Tom O'Haver 5/21/2023

 file= uigetfile('*.csv;*.xlsx');% Click this button to load data
mydata=xlsread(file); %                       from disk in xlsx or csv format.

%

ListPeaks=true; % Print out list of detected peak parameters (position, height, width, area)
LabelPeaks=true; % Check this box to number all detected peaks on the plot% Set the x-axis scale expansion beginning and end points (0 - 100%)
startpc=20.8; % Percentage of data points to start data selection
endpc=50.9; % Percentage of data points to end data selection

% Double-click any of the sliders to change their ranges if necessary
SymmetrizeFactor=0; % First-derivative "de-tailing" for skewed peaks
SmoothWidth=5; % Smoothing to skip false noise peaks. 
NumPasses=4; % Number of sliding average passes (1-5)

SharpenPeaks=false; % Click to sharpen peaks by Fourier self-deconvolution
DecovolutionPeakShape="Gaussian";
PCGaussian=0; % Percent Gaussian (G/L blend)
% Increasing this increases the sharpening but also increases noise
fdw=0.68 ; % Deconvolution width, as a fraction of width of narrowest peak
DA=4.2; % Denominator addition, helps to reduce noise and ringing
ShowRawData=false; % Show raW data points (red) with processed data (blue)
PeakDetector="PeaksAboveThreshold" ;
threshold=0; % Amplitude threshold, % of maximun
PeakSlope=5; % Lower values detect wider peaks (findpeaksG only)
PeakGroup=9; % Number of data points that defines the peak top

% Controls for optional least-squares peak fitting
FitDetectedPeaks=true; % Perform curve fit only detected peaks BEFORE SMOOTHING
PeakShape="exponentially broadened Gaussian" ;
NumTrials=1; % Takes best of "NumTrials" curve fit trials
extra=12.2; % Shape parameter for variable-shape peaks (EMG, G/L blend, etc)
% -------------------------------------------------------------------------------
format short
NumPoints=length(mydata);
startpoint=round(startpc*NumPoints./100);
endpoint=round(endpc*NumPoints./100);
x=mydata(startpoint:endpoint,1);
y=mydata(startpoint:endpoint,2);


% Symmetrize by weighted first-derivative addition
ysym=y+SymmetrizeFactor.*derivxy(x,y);

% Recursive sliding average smoothing
sy=fastsmooth(ysym,SmoothWidth,NumPasses,1);


% Sharpen peaks
if SharpenPeaks==true
    MaxSignal=max(sy);
    P=peaksatG(x,sy,MaxSignal*threshold./100,PeakGroup);
    minw=min(P(:,4));
    dw=fdw.*minw;
    y=sy;
    if DecovolutionPeakShape=="Gaussian"
        df=gaussian(x,min(x),dw)+gaussian(x,max(x),dw);
        ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
    end

    if DecovolutionPeakShape=="Lorentzian"
        df=lorentzian(x,min(x),dw)+lorentzian(x,max(x),dw);
        ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
    end

    if DecovolutionPeakShape=="Logistic"
        df=logistic(x,min(x),dw)+logistic(x,max(x),dw);
        ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
    end

    if DecovolutionPeakShape=="GL blend"
        df=GL(x,min(x),dw,PCGaussian)+GL(x,max(x),dw,PCGaussian);
        ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
    end

    % Apply only denominator addition to reduce noise and ringing
    fftc=fft(df);
    % Deconvolution by fft/ifft with Denominator addition
    sy=ifft(fft(y)./(fftc+DA.*0.01.*max(fftc))).*sum(df);
end

if PeakDetector=="Allpeaks"
    P=allpeaksw(x,sy);
    NumPeaks=length(P);
    MaxSignal=max(sy);
    disp([file   '   Peak detector: allpeaksw.m'])
    format short
    disp('     peak    position    height    Width')
    disp(P);
end
if PeakDetector=="PeaksAboveThreshold"
    MaxSignal=max(sy);
    P=peaksatG(x,sy,MaxSignal*threshold./100,PeakGroup);
    sizeP=size(P);
    NumPeaks=sizeP(1);
    disp([file   '   Peak detector: peaksatG.m'])
    if ListPeaks==true
        format short
        disp('    peak #   Position   Height     Width     Area')
        for n=1:NumPeaks
            disp(['    ' sprintf('%0.4g',n) '         ' sprintf('%0.4g',P(n,2)) '    '  sprintf('%0.4g',P(n,3))  '    '  sprintf('%0.4g',P(n,4)) '     '  sprintf('%0.4g',P(n,5))]);
        end
    end
end
if PeakDetector=="AutoPeaks"
    MaxSignal=max(sy);
    P=autopeaks(x,sy,PeakGroup);
    sizeP=size(P);
    NumPeaks=sizeP(1);
    disp([file   '   Peak detector: AutoPeaks.m'])
    if ListPeaks==true
        format short
        disp('Peak      Position    PeakMax   Peak-valley  Perp drop  Tan skim');
        for n=1:NumPeaks
            disp(['    ' sprintf('%0.4g',n) '         ' sprintf('%0.4g',P(n,2)) '    '  sprintf('%0.4g',P(n,3))  '    '  sprintf('%0.4g',P(n,4)) '      '  sprintf('%0.4g',P(n,5)) '      '  sprintf('%0.4g',P(n,6))]);
        end
    end
end

if PeakDetector=="findpeaksG"
    MaxSignal=max(sy);
    WidthPoints=length(x)/PeakSlope;
    SlopeThreshold=0.7*WidthPoints^-2;
    P=findpeaksG(x,sy,SlopeThreshold,MaxSignal*threshold./100,1,PeakGroup);
    sizeP=size(P);
    NumPeaks=sizeP(1);
    disp([file   '   Peak detector: findpeaksG.m'])
    if ListPeaks==true
        format short
        disp('    peak #   Position   Height     Width     Area')
        for n=1:NumPeaks
            disp(['    ' sprintf('%0.4g',n) '         ' sprintf('%0.4g',P(n,2)) '    '  sprintf('%0.4g',P(n,3))  '    '  sprintf('%0.4g',P(n,4)) '     '  sprintf('%0.4g',P(n,5))]);
        end
    end
end

figure(1)
clf
if LabelPeaks==true
    if ShowRawData==true
        plot(x,sy,x,y,'.')
        legend('processed data','Raw data')
    else
        plot(x,sy);
    end
    if SharpenPeaks==true
        title('Smoothed and Sharpened signal');
    else
        title('Selected data region');
    end
    grid
    text(P(:,2),P(:,3),num2str(P(:,1))) % Numbers peaks
else
    if ShowRawData==true
        plot(x,sy,x,y,'.')
        legend('processed data','Raw data')
    else
        plot(x,sy);
    end
    if SharpenPeaks==true
        title('Smoothed and Sharpened signal');
    else
        title('Selected data region');
    end
    grid
end
xlabel(['Number of peaks detected: ' num2str(NumPeaks) ])

if FitDetectedPeaks==true
    P=peaksatG(x,sy,MaxSignal*threshold./100,PeakGroup);
    % Calculate first guess vector from peak detector results
    FirstGuess=[];
    for peak=1:NumPeaks
        FirstGuess=[FirstGuess P(peak,2) P(peak,4)];
    end
    % disp(FirstGuess)
    switch PeakShape
        case 'Gaussian'
            ShapeNumber=1 ;
        case 'Lorentzian'
            ShapeNumber=2 ;
        case 'Logistic distribution'
            ShapeNumber=3 ;
        case 'Pearson'
            ShapeNumber=4 ;
        case 'exponentially broadened Gaussian'
            ShapeNumber=5 ;
        case 'equal-width Gaussians'
            ShapeNumber=6 ;
        case 'equal-width Lorentzians'
            ShapeNumber=7 ;
        case 'exponentially broadened equal-width Gaussians'
            ShapeNumber=8 ;
        case 'exponential pulse'
            ShapeNumber=9 ;
        case 'Gaussian/Lorentzian blend'
            ShapeNumber=13 ;
    end
    figure(2)
    if SharpenPeaks==true
        [P,GOF]=peakfit([x,sy],0,0,NumPeaks,ShapeNumber,extra,NumTrials,FirstGuess,0,0,1);
    else
        if SymmetrizeFactor==0
            [P,GOF]=peakfit([x,y],0,0,NumPeaks,ShapeNumber,extra,NumTrials,FirstGuess,0,0,1);
        else
            [P,GOF]=peakfit([x,ysym],0,0,NumPeaks,ShapeNumber,extra,NumTrials,FirstGuess,0,0,1);
        end
    end

    if ListPeaks==true
        disp(['Fitting error=' num2str(GOF(1)) '%    R2= ' num2str(GOF(2))])
        disp('   peak #    Position  Height    Width     Area')
        for n=1:NumPeaks
            disp(['    ' sprintf('%0.4g',n) '         ' sprintf('%0.4g',P(n,2)) '    '  sprintf('%0.4g',P(n,3))  '    '  sprintf('%0.4g',P(n,4)) '     '  sprintf('%0.4g',P(n,5))]);
        end
    end
end

%

function SmoothY=fastsmooth(Y,w,type,ends)
% fastsmooth(Y,w,type,ends) smooths vector Y with smooth
%  of width w. Version 3.0, October 2016.
% The argument "type" determines the smooth type:
%   If type=1, rectangular (sliding-average or boxcar)
%   If type=2, triangular (2 passes of sliding-average)
%   If type=3, pseudo-Gaussian (3 passes of sliding-average)
%   If type=4, pseudo-Gaussian (4 passes of same sliding-average)
%   If type=5, multiple-width (4 passes of different sliding-average)
% The argument "ends" controls how the "ends" of the signal
% (the first w/2 points and the last w/2 points) are handled.
%   If ends=0, the ends are zero.  (In this mode the elapsed
%     time is independent of the smooth width). The fastest.
%   If ends=1, the ends are smoothed with progressively
%     smaller smooths the closer to the end. (In this mode the
%     elapsed time increases with increasing smooth widths).
% fastsmooth(Y,w,type) smooths with ends=0.
% fastsmooth(Y,w) smooths with type=1 and ends=0.
%
% Examples:
% fastsmooth([1 1 1 10 10 10 1 1 1 1],3)= [0 1 4 7 10 7 4 1 1 0]
% fastsmooth([1 1 1 10 10 10 1 1 1 1],3,1,1)= [1 1 4 7 10 7 4 1 1 1]
%
% x=1:100;
% y=randn(size(x));
% plot(x,y,x,fastsmooth(y,5,3,1),'r')
% xlabel('Blue: white noise.    Red: smoothed white noise.')
%
% Demonstration scripts and functions: SmoothExperiment.m,
% SmoothWidthTest.m, smoothdemo.m, SmoothOptimization.m,
% DemoSegmentedSmooth.m, DeltaTest.m, SmoothVsCurvefit.m
%
% Related functions: SegmentedSmooth.m, medianfilter.m, killspikes.m
% Download from
% https://terpconnect.umd.edu/~toh/spectrum/functions.html

% Copyright (c) 2012,2018 Thomas C. O'Haver
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
if nargin==2, ends=0; type=1; end
if nargin==3, ends=0; end
switch type
    case 1
        SmoothY=sa(Y,w,ends);
    case 2
        SmoothY=sa(sa(Y,w,ends),w,ends);
    case 3
        SmoothY=sa(sa(sa(Y,w,ends),w,ends),w,ends);
    case 4
        SmoothY=sa(sa(sa(sa(Y,w,ends),w,ends),w,ends),w,ends);
    case 5
        SmoothY=sa(sa(sa(sa(Y,round(1.6*w),ends),round(1.4*w),ends),round(1.2*w),ends),w,ends);
end
end

function SmoothY=sa(Y,smoothwidth,ends)
w=round(smoothwidth);
SumPoints=sum(Y(1:w));
s=zeros(size(Y));
halfw=round(w/2);
L=length(Y);
for k=1:L-w
    s(k+halfw-1)=SumPoints;
    SumPoints=SumPoints-Y(k);
    SumPoints=SumPoints+Y(k+w);
end
s(k+halfw)=sum(Y(L-w+1:L));
SmoothY=s./w;
% Taper the ends of the signal if ends=1.
if ends==1
    startpoint=(smoothwidth + 1)/2;
    SmoothY(1)=(Y(1)+Y(2))./2;
    for k=2:startpoint
        SmoothY(k)=mean(Y(1:(2*k-1)));
        SmoothY(L-k+1)=mean(Y(L-2*k+2:L));
    end
    SmoothY(L)=(Y(L)+Y(L-1))./2;
end
end

function P=peaksatG(x,y,AmpThreshold,peakgroup)
% function P=peaksatG(x,y,AmpThreshold,peakgroup)
% peaksatG (peaks Above Threshold/Gaussian) lists every y value that (a)
% has lower y values on both sides and (b) is above the specified
% threshold. Returns a 2 x n matrix P with the x and y values of each peak,
% where n is the number of peaks detected. (Noisy data must be smoothed
% first). "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.
% See http://terpconnect.umd.edu/~toh/spectrum/Smoothing.html and
% http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm
%
% Examples:
% peaksatG(0:.01:2,humps(0:.01:2),0,5)
% x=[0:.01:50];y=(1+cos(x)).^2;P=peaksatG(x,y,0,5)
% x=[0:.01:5];y=x.*sin(x.^2).^2;P=peaksatG(x,y,-1,5)
% x=[-10:.1:10];y=exp(-(x).^2);P=peaksatG(x,y,0.3,3)
%
% Find, measure, and plot peaks with variable unknown positions
%   x=-50:.2:50;
%   y=exp(-(x).^2)+exp(-(x+50*rand()).^2);
%   plot(x,y,'m.')
%   P=peaksatG(x,y,0.2,5);
%   text(P(:,2),P(:,3),num2str(P(:,1)))
%   disp('           peak #    Position    Height        Width         Area')
%   disp(P)
%
% Related functions:
% findpeaksG.m, findpeaksSG.m, findvalleys.m, findpeaksL.m, findpeaksb.m,
% findpeaksb3.m, findpeaksplot.m, peakstats.m, findpeaksnr.m,
% findpeaksGSS.m, ipeak.m findpeaksLSS.m, findpeaksfit.m, findsteps.m,
% findsquarepulse.m, idpeaks.m

% Copyright (c), 2022 Thomas C. O'Haver
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%
peakgroup=round(peakgroup);
n=round(peakgroup/2+1);
P=[0 0 0 0 0];
vectorlength=length(y);
peak=1;
for j=2:length(y)-1
    if y(j) > y(j-1)
        if y(j) > y(j+1)
            if y(j) > AmpThreshold  % 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 % for k=1:peakgroup,...
                if peakgroup>2
                    [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 % if peakgroup>2,...
                % Construct matrix P. One row for each peak detected,
                % containing the peak number, peak position (x-value) and
                % peak height (y-value). If peak measurement fails and
                % results in NaN, or if the measured peak height is less
                % than AmpThreshold, skip this peak
                if isnan(PeakX) || isnan(PeakY) || PeakY<AmpThreshold
                    % Skip this peak
                else % Otherwise count this as a valid peak
                    P(peak,:) = [round(peak) PeakX PeakY MeasuredWidth  1.0646.*PeakY*MeasuredWidth];
                    peak=peak+1; % Move on to next peak
                end % if isnan...
            end % if y(j) > AmpThreshold,...
        end
    end
end % for j=2*round(smoothwidth/2)-1:...
end

function [index,closestval]=val2ind(x,val)
% Returns the index and the value of the element of vector x that is closest to val
% If more than one element is equally close, returns vectors of indicies and values
% Tom O'Haver (toh@umd.edu) October 2006
% Examples: If x=[1 2 4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5 9]
% [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3]
dif=abs(x-val);
index=find((dif-min(dif))==0);
closestval=x(index);
end

function [Height, Position, Width]=gaussfit(x,y)
% Converts y-axis to a log scale, fits a parabola
% (quadratic) to the (x,ln(y)) data, then calculates
% the position, width, and height of the
% Gaussian from the three coefficients of the
% quadratic fit.  This is accurate only if the data have
% no baseline offset (that is, trends to zero far off the
% peak) and if there are no zeros or negative values in y.
%
% Example 1: Simplest Gaussian data set
% [Height, Position, Width]=gaussfit([1 2 3],[1 2 1])
%    returns Height = 2, Position = 2, Width = 2
%
% Example 2: best fit to synthetic noisy Gaussian
% x=50:150;y=100.*gaussian(x,100,100)+10.*randn(size(x));
% [Height,Position,Width]=gaussfit(x,y)
%   returns [Height,Position,Width] clustered around 100,100,100.
%
% Example 3: plots data set as points and best-fit Gaussian as line
% x=[1 2 3 4 5];y=[1 2 2.5 2 1];
% [Height,Position,Width]=gaussfit(x,y)
% plot(x,y,'o',linspace(0,8),Height.*gaussian(linspace(0,8),Position,Width))
% Copyright (c) 2012, Thomas C. O'Haver

% To prevent problems from taking the log of zero or negative values,
% make the lowest value of y equal to 1% of the maximum value.
maxy=max(y);
for p=1:length(y)
    if y(p)<(maxy/100),y(p)=maxy/100;end
end % end of for p=1:length(y),
logabsy=log(abs(y));
sizex=size(x);
sizey=size(logabsy);
if (sizex(1)==sizey(1))
    [coef,~,MU]=polyfit(x',logabsy',2);
else
    [coef,~,MU]=polyfit(x,logabsy',2);
end
c1=coef(3);c2=coef(2);c3=coef(1);
% Compute peak position and height or fitted parabola
Position=-((MU(2).*c2/(2*c3))-MU(1));
Height=exp(c1-c3*(c2/(2*c3))^2);
Width=norm(MU(2).*2.35703/(sqrt(2)*sqrt(-1*c3)));
end

function [M,A]=autopeaks(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype)
% Peak detection and height and area measurement for peaks of arbitrary
% shape in x,y time series data. Returns the peak table in the matrix P and
% the peak detection parameters in the vector A. The syntax autopeaks(x,y)
% works well in some cases, but if not try autopeaks(x,y,n), using
% different values of n (roughly the number of peaks that would fit into
% the signal record) until it detects the peaks that you want to measure.
% For the most precise control over peak detection, you can specify all the
% peak detection parameters by typing
% M=autopeaks(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup)
%
% The syntax is similar syntax to autofindpeaks.m, in that that the peak
% detection parameters SlopeThreshold, AmpThreshold, smoothwidth,
% peakgroup, smoothtype can be omitted and the function will calculate
% initial values based on the number of data points and/or the optional
% peak density n (the last of 3 arguments after x and y). If you do specify
% all the peak detection parameters, they can be scalars or vectors with
% one entry for each segment (like findpeaskSG.m).
%
% Locates and measures the positive peaks in a noisy x-y time series data.
% Detects peaks by looking for downward zero-crossings in the first
% derivative whose upward slopes exceed SlopeThreshold. Returns matrix (M)
% containing peak number, the absolute peak height, peak-valley difference,
% perpendicular drop area, and the tangent skim area of each peak.
% Arguments "SlopeThreshold", "AmpThreshold" and "smoothwidth" control peak
% sensitivity of each segment. Higher values will neglect smaller features.
% "smoothwidth" is a vector of the widths of the smooths applied before
% peak detection; larger values ignore narrow peaks. If smoothwidth=0, no
% smoothing is performed. "peakgroup" is a vector of the number of points
% around the top part of the peak that are taken for initial estimate of
% the peak center and width (minimum 3). Note: this function uses
% smoothing only for peak detection; it performs measurements on the raw
% unsmoothed y data. If the data are noisy, it may be beneficial to smooth
% the y data yourself before calling autopeaks.m, using any smooth function
% of your choice.
%
% See https://terpconnect.umd.edu/~toh/spectrum/Integration.html
% http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm
% (c) T.C. O'Haver, 2017.  Version 1, January, 2017
%
% The script testautopeaks.m will run all of the following examples with a
% 1-second pause between each, plotting the data and the peaks (by itself
% autopeaks.m does not plot, use autopeaksplot.m for that):
%
% Example 1: sin(x).^2 has theoretical peaks at x=0.5pi, 1.5pi, 2.5pi,
% 3.5pi..., with peak heights of 1.0 and peak areas of pi/2 = 1.5708.
%  disp('          Peak      Position    PeakMax    Peak-valley    Perp drop   Tan skim')
%  x=[0:.01:20]';y=sin(x).^2;autopeaks(x,y,0,0,5,5,1)
%
% Example 2: The built-in "humps" function has two peaks, no noise.
%  autopeaks(0:.01:2,humps(0:.01:2),0,0,1,1,1)
%
% Example 3: Series of peaks that get progressively taller and wider.
%  x=[0:.01:5]';y=x.*sin(x.^2).^2;autopeaks(x,y,0,0,5,5,1)
%
% Example 4: Like example 3, with random white noise added.
%  x=[0:.01:5]';y=.1.*randn(size(x))+x.*sin(-x.^2).^2;
%  autopeaks(x,y,.001,.5,15,15,1)
%
% Example 5; Like example 3, with added rising baseline.
%  x=[0:.01:5]';y=x+x.*sin(x.^2).^2;autopeaks(x,y,0,0,5,5,1)
%
% Example 6: Gaussian on linear baseline, theoretical area = 1.7725
%  x=[0:.1:10];y=2+x/10+exp(-(x-5).^2)+.01.*randn(size(x));
%  disp('          Peak      Position    PeakMax    Peak-valley    Perp drop   Tan skim')
%  autopeaks(x,y,0.0001,0.3,3,5,1)
%  Peak-valley and Tan skim are most accurate height and area measures.
%
% Example 7: Two overlapping Gaussians, zero baseline, theoretical area = 1.7725
%  x=[0:.05:10];y=exp(-(x-6).^2)+exp(-(x-3.5).^2)+.01.*randn(size(x));
%  disp('          Peak      Position    PeakMax    Peak-valley    Perp drop   Tan skim')
%  autopeaks(x,y,0.0001,0.5,8,6,1)
%  PeakMax and Perp drop are most accurate height and area measures.
%
% Example 8: Narrow Gaussian peak on sloping linear baseline, followed by
% much broader peak, theoretical heights 1 and 1, areas 1.6 and 31.9.
%  x=1:.2:150;
%  y=(150-x)./200+gaussian(x,30,1.5)+gaussian(x,90,30);
%  P=autopeaks(x,y,[0.001 .00005],[.6 .6],[5 120],[8 150],1)
%  PeakAreaError1=100.*((1.6-P(1,6))./1.6)
%  PeakAreaError2=100.*((31.9-P(2,6))./31.9)
%
% Similar functions: measurepeaks.m, autopeaksplot.m, autofindpeaks.m

% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%
switch nargin
    % 'nargin' is the number of input arguments
    case 1  % One argument only
        % Assumne that the argument must be a matrix of data.
        % If DataMatrix is in the wrong transposition, fix it.
        datasize=size(x);
        if datasize(1)<datasize(2),x=x';end
        datasize=size(x);
        if datasize(2)==1, %  Must be ipeak(y-vector)
            y=x;
            x=[1:length(x)]'; % Create an independent variable vector
        else
            % Must be ipeak(DataMatrix)
            y=x(:,2);
            x=x(:,1); % Split matrix argument
        end
        % Calculate default values of peak detection parameters
        PeakDensity=20;
        % Estimate approx number of points in a peak half-width
        WidthPoints=length(y)/PeakDensity;
        SlopeThreshold=WidthPoints^-2;
        AmpThreshold=abs(min(y)+0.5*(max(y)-min(y)));
        smoothwidth=10;
        peakgroup=10;
        if peakgroup>100,peakgroup=100;end   % Keep peakgroup below 100
        if smoothwidth>100,smoothwidth=100;end   % Keep smoothwidth below 100
        smoothtype=3;
        A=[SlopeThreshold AmpThreshold smoothwidth peakgroup];
    case 2
        PeakDensity=20;
        % Estimate approx number of points in a peak half-width
        WidthPoints=length(y)/PeakDensity;
        SlopeThreshold=WidthPoints^-2;
        AmpThreshold=abs(min(y)+0.1*(max(y)-min(y)));
        smoothwidth=round(WidthPoints/3);
        peakgroup=round(WidthPoints/3);
        if peakgroup>100,peakgroup=100;end   % Keep peakgroup below 100
        if smoothwidth>100,smoothwidth=100;end   % Keep smoothwidth below 100
        smoothtype=3;
        A=[SlopeThreshold AmpThreshold smoothwidth peakgroup];
    case 3
        % Must be separate x and y data vectors plus a peak density
        % estimate.
        % Calculate values of peak detection parameters
        % arguments based on the peak density, PeakD
        PeakDensity=SlopeThreshold;
        % Estimate approx number of points in a peak half-width
        WidthPoints=length(y)/PeakDensity;
        SlopeThreshold=WidthPoints^-2;
        AmpThreshold=abs(min(y)+0.02*(max(y)-min(y)));
        smoothwidth=round(WidthPoints/3);
        peakgroup=round(WidthPoints/3);
        if peakgroup>100,peakgroup=100;end   % Keep peakgroup below 100
        if smoothwidth>100,smoothwidth=100;end   % Keep smoothwidth below 100
        smoothtype=3;
        A=[SlopeThreshold AmpThreshold smoothwidth peakgroup];
end

if peakgroup<3,peakgroup=3;end
% Find the peaks in the data
P=findpeaksSG(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype);
sizeP=size(P);
NumPeaks=sizeP(1);
clear PeakHeight
clear M

% M=zeros(NumPeaks,6);
goodpeak=1;
for PeakNumber=1:NumPeaks
    PeakPos=P(PeakNumber,2);
    PeakWidth=P(PeakNumber,4);
    PeakIndex=val2ind(x,PeakPos);
    x1=val2ind(x,PeakPos-1.3.*PeakWidth); % Left-most end of peak
    x2=val2ind(x,PeakPos+1.3.*PeakWidth);  % Right-most end of peak
    Halfway=val2ind(x,PeakPos); % x-value of center of peak
    miny=min(y(x1:x2)); % Minimum y value of that region
    yy(x1:x2)=y(x1:x2)-miny;  % Subtract minimum from that region
    % Calculate height, peak-valley, and perpendicular drop area
    PeakMax=y(PeakIndex);
    LeftSideMinimum=min(yy(x1:Halfway));
    RightSideMinimum=min(yy(Halfway:x2));
    LeftValleyIndex=x1+val2ind(yy(x1:Halfway),LeftSideMinimum);
    RightValleyIndex=PeakIndex+val2ind(yy(Halfway:x2),RightSideMinimum);
    if RightValleyIndex(1)>length(x),RightValleyIndex(1)=length(x);end
    PDarea=sum(y(LeftValleyIndex(1):RightValleyIndex(1))).*(x(2)-x(1)); % Perpendicular Drop area
    PVHeight=PeakMax-(interp1([x(LeftValleyIndex(1)) x(RightValleyIndex(1))],[y(LeftValleyIndex(1)) y(RightValleyIndex(1))],PeakPos));
    interpy=interp1([x(LeftValleyIndex(1)) x(RightValleyIndex(1))],[y(LeftValleyIndex(1)) y(RightValleyIndex(1))],x(LeftValleyIndex(1):RightValleyIndex(1)));
    TSarea=sum((y(LeftValleyIndex(1):RightValleyIndex(1))-interpy).*(x(x1+1)-x(x1)));
    % [PeakNumber PeakPos PeakMax PVHeight]
    if PVHeight>AmpThreshold,
        M(goodpeak,:)=[goodpeak PeakPos  PeakMax PVHeight PDarea TSarea];
        goodpeak=goodpeak+1;
    end
end % for PeakNumber
x=1:.2:100;
y=gaussian(x,40,10)+gaussian(x,50,10)+.01.*randn(size(x));
A=[SlopeThreshold AmpThreshold smoothwidth peakgroup];
end
