# v_findpeaks

## PURPOSE

V_FINDPEAKS finds peaks with optional quadratic interpolation [K,V]=(Y,M,W,X)

## SYNOPSIS

function [k,v]=v_findpeaks(y,m,w,x)

## DESCRIPTION

    V_FINDPEAKS finds peaks with optional quadratic interpolation [K,V]=(Y,M,W,X)

      Inputs:  Y(N,1)   is the input signal (does not work with UInt datatype)
               M        is mode:
                           'f' include the first sample if a downward initial slope
                           'l' include the last sample if an upward final slope
                           'm' return only the maximum peak
                           'q' performs quadratic interpolation
                           'v' finds valleys instead of peaks
               W        is the width tolerance; a peak will be eliminated if there is
                        a higher peak within +-W. Units are samples or X values
               X(N,1)   x-axis locations of Y values [default: 1:length(Y)]

     Outputs:  K(P,1)   are the positions in X of the peaks in Y (fractional if M='q')
               V(P,1)   are the peak amplitudes: if M='q' the amplitudes will be interpolated
                        whereas if M~='q' then V=Y(K).

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_fxpefac V_FXPEFAC PEFAC pitch tracker [FX,TT,PV,FV]=(S,FS,TINC,M,PP)
- v_fxrapt V_FXRAPT RAPT pitch tracker [FX,VUV]=(S,FS,M,Q)
- v_windinfo V_WINDINFO window information and figures of merit X=(W,FS)

## SOURCE CODE

    0001 function [k,v]=v_findpeaks(y,m,w,x)
    0002 %V_FINDPEAKS finds peaks with optional quadratic interpolation [K,V]=(Y,M,W,X)
    0003 %
    0004 %  Inputs:  Y(N,1)   is the input signal (does not work with UInt datatype)
    0005 %           M        is mode:
    0006 %                       'f' include the first sample if a downward initial slope
    0007 %                       'l' include the last sample if an upward final slope
    0008 %                       'm' return only the maximum peak
    0009 %                       'q' performs quadratic interpolation
    0010 %                       'v' finds valleys instead of peaks
    0011 %           W        is the width tolerance; a peak will be eliminated if there is
    0012 %                    a higher peak within +-W. Units are samples or X values
    0013 %           X(N,1)   x-axis locations of Y values [default: 1:length(Y)]
    0014 %
    0015 % Outputs:  K(P,1)   are the positions in X of the peaks in Y (fractional if M='q')
    0016 %           V(P,1)   are the peak amplitudes: if M='q' the amplitudes will be interpolated
    0017 %                    whereas if M~='q' then V=Y(K).
    0018
    0019 % Outputs are column vectors regardless of whether Y is row or column.
    0020 % If there is a plateau rather than a sharp peak, the routine will place the
    0021 % peak in the centre of the plateau. When the W input argument is specified,
    0022 % the routine will eliminate the lower of any pair of peaks whose separation
    0023 % is <=W; if the peaks have exactly the same height, the second one will be eliminated.
    0024 % Unless the 'f' or 'l' options are given, all peak locations satisfy 1<K<N.
    0025 %
    0026 % If no output arguments are specified, the results will be plotted.
    0027 %
    0028
    0029 %       Copyright (C) Mike Brookes 2005
    0030 %      Version: $Id: v_findpeaks.m 6564 2015-08-16 16:56:40Z dmb $
    0031 %
    0032 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0033 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0034 %
    0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0036 %   This program is free software; you can redistribute it and/or modify
    0037 %   it under the terms of the GNU Lesser General Public License as published by
    0038 %   the Free Software Foundation; either version 3 of the License, or
    0039 %   (at your option) any later version.
    0040 %
    0041 %   This program is distributed in the hope that it will be useful,
    0042 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0043 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0044 %   GNU Lesser General Public License for more details.
    0045 %
    0046 %   You can obtain a copy of the GNU Lesser General Public License from
    0047 %   https://www.gnu.org/licenses/ .
    0048 %
    0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0050
    0051 if nargin<2 || ~numel(m)
    0052     m=' ';
    0053 elseif nargin>3
    0054     x=x(:); % x must be a column vector
    0055 end
    0056 ny=length(y);
    0057 if any(m=='v')
    0058     y=-y(:);        % invert y if searching for valleys
    0059 else
    0060     y=y(:);        % force to be a column vector
    0061 end
    0062 dx=y(2:end)-y(1:end-1);
    0063 r=find(dx>0);
    0064 f=find(dx<0);
    0065 k=[]; % set defaults
    0066 v=[];
    0067 if ~isempty(r) && ~isempty(f)    % we must have at least one rise and one fall
    0068     dr=r;
    0069     dr(2:end)=r(2:end)-r(1:end-1);
    0070     rc=ones(ny,1);
    0071     rc(r+1)=1-dr;
    0072     rc(1)=0;
    0073     rs=cumsum(rc); % = time since the last rise
    0074     df=f;
    0075     df(2:end)=f(2:end)-f(1:end-1);
    0076     fc=ones(ny,1);
    0077     fc(f+1)=1-df;
    0078     fc(1)=0;
    0079     fs=cumsum(fc); % = time since the last fall
    0080     rp=repmat(-1,ny,1);
    0081     rp([1; r+1])=[dr-1; ny-r(end)-1];
    0082     rq=cumsum(rp);  % = time to the next rise
    0083     fp=repmat(-1,ny,1);
    0084     fp([1; f+1])=[df-1; ny-f(end)-1];
    0085     fq=cumsum(fp); % = time to the next fall
    0086     k=find((rs<fs) & (fq<rq) & (floor((fq-rs)/2)==0));   % the final term centres peaks within a plateau
    0087     v=y(k);
    0088     if any(m=='q')         % do quadratic interpolation
    0089         if nargin>3
    0090             xm=x(k-1)-x(k);
    0091             xp=x(k+1)-x(k);
    0092             ym=y(k-1)-y(k);
    0093             yp=y(k+1)-y(k);
    0094             d=xm.*xp.*(xm-xp);
    0095             b=0.5*(yp.*xm.^2-ym.*xp.^2);
    0096             a=xm.*yp-xp.*ym;
    0097             j=(a>0);            % j=0 on a plateau
    0098             v(j)=y(k(j))+b(j).^2./(a(j).*d(j));
    0099             k(j)=x(k(j))+b(j)./a(j); % x-axis position of peak
    0100             k(~j)=0.5*(x(k(~j)+fq(k(~j)))+x(k(~j)-rs(k(~j))));    % find the middle of the plateau
    0101         else
    0102             b=0.25*(y(k+1)-y(k-1));
    0103             a=y(k)-2*b-y(k-1);
    0104             j=(a>0);            % j=0 on a plateau
    0105             v(j)=y(k(j))+b(j).^2./a(j);
    0106             k(j)=k(j)+b(j)./a(j);
    0107             k(~j)=k(~j)+(fq(k(~j))-rs(k(~j)))/2;    % add 0.5 to k if plateau has an even width
    0108         end
    0109     elseif nargin>3 % convert to the x-axis using linear interpolation
    0110         k=x(k);
    0111     end
    0112 end
    0113 % add first and last samples if requested
    0114 if ny>1
    0115     if any(m=='f') && y(1)>y(2)
    0116         v=[y(1); v];
    0117         if nargin>3
    0118             k=[x(1); k];
    0119         else
    0120             k=[1; k];
    0121         end
    0122     end
    0123     if any(m=='l') && y(ny-1)<y(ny)
    0124         v=[v; y(ny)];
    0125         if nargin>3
    0126             k=[k; x(ny)];
    0127         else
    0128             k=[k; ny];
    0129         end
    0130     end
    0131
    0132     % now purge nearby peaks - note that the decision about which peaks to
    0133     % delete is not unique
    0134
    0135     if any(m=='m')
    0136         [v,iv]=max(v);
    0137         k=k(iv);
    0138     elseif nargin>2 && numel(w)==1 && w>0
    0139         j=find(k(2:end)-k(1:end-1)<=w);
    0140         while any(j)
    0141             j=j+(v(j)>=v(j+1));
    0142             k(j)=[];
    0143             v(j)=[];
    0144             j=find(k(2:end)-k(1:end-1)<=w);
    0145         end
    0146     end
    0147 elseif ny>0 && (any(m=='f') || any(m=='l'))
    0148     v=y;
    0149     if nargin>3
    0150         k=x;
    0151     else
    0152         k=1;
    0153     end
    0154 end
    0155 if any(m=='v')
    0156     v=-v;    % invert peaks if searching for valleys
    0157 end
    0158
    0159 if ~nargout
    0160     if any(m=='v')
    0161         y=-y;    % re-invert y if searching for valleys
    0162         ch='v';
    0163     else
    0164         ch='^';
    0165     end
    0166     plot(1:ny,y,'-',k,v,ch);
    0167 end