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).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