Home > voicebox > v_findpeaks.m

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:

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 General Public License as published by
0038 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0045 %
0046 %   You can obtain a copy of the GNU General Public License from
0047 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0048 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003