这里写自定义目录标题
Matlab v_findpeaks代码
写在前面
本函数主要用于寻找数据的上下极值点,可以用于求取包络
参考链接:http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/mdoc/v_mfiles/v_findpeaks.html#_top
代码
function [k,v]=v_findpeaks(y,m,w,x)
%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).
% Outputs are column vectors regardless of whether Y is row or column.
% If there is a plateau rather than a sharp peak, the routine will place the
% peak in the centre of the plateau. When the W input argument is specified,
% the routine will eliminate the lower of any pair of peaks whose separation
% is <=W; if the peaks have exactly the same height, the second one will be eliminated.
% Unless the 'f' or 'l' options are given, all peak locations satisfy 1<K<N.
%
% If no output arguments are specified, the results will be plotted.
%
% Copyright (C) Mike Brookes 2005
% Version: $Id: v_findpeaks.m 6564 2015-08-16 16:56:40Z dmb $
%
% VOICEBOX is a MATLAB toolbox for speech processing.
% Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You can obtain a copy of the GNU General Public License from
% http://www.gnu.org/copyleft/gpl.html or by writing to
% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2 || ~numel(m)
m=' ';
elseif nargin>3
x=x(:); % x must be a column vector
end
ny=length(y);
if any(m=='v')
y=-y(:); % invert y if searching for valleys
else
y=y(:); % force to be a column vector
end
dx=y(2:end)-y(1:end-1);
r=find(dx>0);
f=find(dx<0);
k=[]; % set defaults
v=[];
if ~isempty(r) && ~isempty(f) % we must have at least one rise and one fall
dr=r;
dr(2:end)=r(2:end)-r(1:end-1);
rc=ones(ny,1);
rc(r+1)=1-dr;
rc(1)=0;
rs=cumsum(rc); % = time since the last rise
df=f;
df(2:end)=f(2:end)-f(1:end-1);
fc=ones(ny,1);
fc(f+1)=1-df;
fc(1)=0;
fs=cumsum(fc); % = time since the last fall
rp=repmat(-1,ny,1);
rp([1; r+1])=[dr-1; ny-r(end)-1];
rq=cumsum(rp); % = time to the next rise
fp=repmat(-1,ny,1);
fp([1; f+1])=[df-1; ny-f(end)-1];
fq=cumsum(fp); % = time to the next fall
k=find((rs<fs) & (fq<rq) & (floor((fq-rs)/2)==0)); % the final term centres peaks within a plateau
v=y(k);
if any(m=='q') % do quadratic interpolation
if nargin>3
xm=x(k-1)-x(k);
xp=x(k+1)-x(k);
ym=y(k-1)-y(k);
yp=y(k+1)-y(k);
d=xm.*xp.*(xm-xp);
b=0.5*(yp.*xm.^2-ym.*xp.^2);
a=xm.*yp-xp.*ym;
j=(a>0); % j=0 on a plateau
v(j)=y(k(j))+b(j).^2./(a(j).*d(j));
k(j)=x(k(j))+b(j)./a(j); % x-axis position of peak
k(~j)=0.5*(x(k(~j)+fq(k(~j)))+x(k(~j)-rs(k(~j)))); % find the middle of the plateau
else
b=0.25*(y(k+1)-y(k-1));
a=y(k)-2*b-y(k-1);
j=(a>0); % j=0 on a plateau
v(j)=y(k(j))+b(j).^2./a(j);
k(j)=k(j)+b(j)./a(j);
k(~j)=k(~j)+(fq(k(~j))-rs(k(~j)))/2; % add 0.5 to k if plateau has an even width
end
elseif nargin>3 % convert to the x-axis using linear interpolation
k=x(k);
end
end
% add first and last samples if requested
if ny>1
if any(m=='f') && y(1)>y(2)
v=[y(1); v];
if nargin>3
k=[x(1); k];
else
k=[1; k];
end
end
if any(m=='l') && y(ny-1)<y(ny)
v=[v; y(ny)];
if nargin>3
k=[k; x(ny)];
else
k=[k; ny];
end
end
% now purge nearby peaks - note that the decision about which peaks to
% delete is not unique
if any(m=='m')
[v,iv]=max(v);
k=k(iv);
elseif nargin>2 && numel(w)==1 && w>0
j=find(k(2:end)-k(1:end-1)<=w);
while any(j)
j=j+(v(j)>=v(j+1));
k(j)=[];
v(j)=[];
j=find(k(2:end)-k(1:end-1)<=w);
end
end
elseif ny>0 && (any(m=='f') || any(m=='l'))
v=y;
if nargin>3
k=x;
else
k=1;
end
end
if any(m=='v')
v=-v; % invert peaks if searching for valleys
end
if ~nargout
if any(m=='v')
y=-y; % re-invert y if searching for valleys
ch='v';
else
ch='^';
end
plot(1:ny,y,'-',k,v,ch);
end