V_MAXFILT find max of an exponentially weighted sliding window [Y,K,Y0]=(X,F,nn,D,X0)
Usage: (1) y=v_maxfilt(x) % maximum filter along first non-singleton dimension
(2) y=v_maxfilt(x,0.95) % use a forgetting factor of 0.95 (= time const of -1/log(0.95)=19.5 samples)
(3) Two equivalent methods (i.e. you can process x in chunks):
y=v_maxfilt([u v]); [yu,ku,x0)=v_maxfilt(u);
yv=v_maxfilt(v,[],[],[],x0);
y=[yu yv];
Inputs: X Vector or matrix of input data
F exponential forgetting factor in the range 0 (very forgetful) to 1 (no forgetting)
F=exp(-1/T) gives a time constant of T samples [default = 1]
n Length of sliding window [default = Inf (equivalent to [])]
D Dimension for work along [default = first non-singleton dimension]
X0 Initial values placed in front of the X data
Outputs: Y Output matrix - same size as X
K Index array: Y=X(K). (Note that these value may be <=0 if input X0 is present)
Y0 Last nn-1 values (used to initialize a subsequent call to
v_maxfilt()) (or last output if n=Inf)
This routine calaulates y(p)=max(f^r*x(p-r), r=0:n-1) where x(r)=-inf for r<1
y=x(k) on output0001 function [y,k,y0]=v_maxfilt(x,f,n,d,x0) 0002 %V_MAXFILT find max of an exponentially weighted sliding window [Y,K,Y0]=(X,F,nn,D,X0) 0003 % 0004 % Usage: (1) y=v_maxfilt(x) % maximum filter along first non-singleton dimension 0005 % (2) y=v_maxfilt(x,0.95) % use a forgetting factor of 0.95 (= time const of -1/log(0.95)=19.5 samples) 0006 % (3) Two equivalent methods (i.e. you can process x in chunks): 0007 % y=v_maxfilt([u v]); [yu,ku,x0)=v_maxfilt(u); 0008 % yv=v_maxfilt(v,[],[],[],x0); 0009 % y=[yu yv]; 0010 % 0011 % Inputs: X Vector or matrix of input data 0012 % F exponential forgetting factor in the range 0 (very forgetful) to 1 (no forgetting) 0013 % F=exp(-1/T) gives a time constant of T samples [default = 1] 0014 % n Length of sliding window [default = Inf (equivalent to [])] 0015 % D Dimension for work along [default = first non-singleton dimension] 0016 % X0 Initial values placed in front of the X data 0017 % 0018 % Outputs: Y Output matrix - same size as X 0019 % K Index array: Y=X(K). (Note that these value may be <=0 if input X0 is present) 0020 % Y0 Last nn-1 values (used to initialize a subsequent call to 0021 % v_maxfilt()) (or last output if n=Inf) 0022 % 0023 % This routine calaulates y(p)=max(f^r*x(p-r), r=0:n-1) where x(r)=-inf for r<1 0024 % y=x(k) on output 0025 0026 % Example: find all peaks in x that are not exceeded within +-w samples 0027 % w=4;m=100;x=rand(m,1);[y,k]=v_maxfilt(x,1,2*w+1);p=find(((1:m)-k)==w);plot(1:m,x,'-',p-w,x(p-w),'+') 0028 0029 % Bugs/suggestion 0030 % (1) if n<0 then centre filter support of the current sample (or make it +-n) 0031 0032 % Copyright (C) Mike Brookes 2003-2014 0033 % Version: $Id: v_maxfilt.m 10865 2018-09-21 17:22:45Z dmb $ 0034 % 0035 % VOICEBOX is a MATLAB toolbox for speech processing. 0036 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0037 % 0038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0039 % This program is free software; you can redistribute it and/or modify 0040 % it under the terms of the GNU General Public License as published by 0041 % the Free Software Foundation; either version 2 of the License, or 0042 % (at your option) any later version. 0043 % 0044 % This program is distributed in the hope that it will be useful, 0045 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0046 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0047 % GNU General Public License for more details. 0048 % 0049 % You can obtain a copy of the GNU General Public License from 0050 % http://www.gnu.org/copyleft/gpl.html or by writing to 0051 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0053 0054 s=size(x); 0055 if nargin<4 || isempty(d) % if d is unspecified 0056 d=find(s>1,1); % find first non-singleton dimension 0057 if isempty(d) 0058 d=1; % else defailt to dimension 1 0059 end 0060 end 0061 if nargin>4 && numel(x0)>0 % initial values specified 0062 y=shiftdim(cat(d,x0,x),d-1); % concatenate x0 and x along d 0063 nx0=size(x0,d); % number of values added onto front of data 0064 else 0065 y=shiftdim(x,d-1); % dimension specified, d 0066 nx0=0; % no values added onto front of data 0067 end 0068 s=size(y); 0069 s1=s(1); % size of active dimension 0070 if nargin<3 || isempty(n) 0071 n0=Inf; % default filter support is Inf 0072 else 0073 n0=max(n,1); 0074 end 0075 if nargin<2 || isempty(f) 0076 f=1; 0077 end 0078 nn=n0; 0079 if nargout>2 % we need to output the tail for next time 0080 if n0<Inf 0081 ny0=min(s1,nn-1); 0082 else 0083 ny0=min(s1,1); 0084 end 0085 sy0=s; 0086 sy0(1)=ny0; 0087 if ny0<=0 || n0==Inf 0088 y0=zeros(sy0); 0089 else 0090 y0=reshape(y(1+s1-ny0:end,:),sy0); 0091 y0=shiftdim(y0,ndims(x)-d+1); 0092 end 0093 end 0094 nn=min(nn,s1); % no point in having nn>s1 0095 k=repmat((1:s1)',[1 s(2:end)]); 0096 if nn>1 0097 j=1; 0098 j2=1; 0099 while j>0 0100 g=f^j; 0101 m=find(y(j+1:s1,:)<=g*y(1:s1-j,:)); 0102 m=m+j*fix((m-1)/(s1-j)); 0103 y(m+j)=g*y(m); 0104 k(m+j)=k(m); 0105 j2=j2+j; 0106 j=min(j2,nn-j2); % j approximately doubles each iteration 0107 end 0108 end 0109 if nargout==0 0110 if nargin<3 0111 x=shiftdim(x); 0112 else 0113 x=shiftdim(x,d-1); 0114 end 0115 ss=min(prod(s(2:end)),5); % maximum of 5 plots 0116 plot(1:s1,reshape(y(nx0+1:end,1:ss),s1,ss),'-r',1:s1,reshape(x(:,1:ss),s1,ss),':b'); 0117 else 0118 if nargout>2 && n0==Inf && ny0==1 % if n0==Inf, we need to save the final output 0119 y0=reshape(y(end,:),sy0); 0120 y0=shiftdim(y0,ndims(x)-d+1); 0121 end 0122 if nx0>0 % pre-data specified, x0 0123 s(1)=s(1)-nx0; 0124 y=shiftdim(reshape(y(nx0+1:end,:),s),ndims(x)-d+1); 0125 k=shiftdim(reshape(k(nx0+1:end,:),s),ndims(x)-d+1)-nx0; 0126 else % no pre-data 0127 y=shiftdim(y,ndims(x)-d+1); 0128 k=shiftdim(k,ndims(x)-d+1); 0129 end 0130 end