Home > voicebox > maxfilt.m

maxfilt

PURPOSE ^

MAXFILT find max of an exponentially weighted sliding window [Y,K,Y0]=(X,F,nn,D,X0)

SYNOPSIS ^

function [y,k,y0]=maxfilt(x,f,n,d,x0)

DESCRIPTION ^

MAXFILT find max of an exponentially weighted sliding window  [Y,K,Y0]=(X,F,nn,D,X0)

 Usage: (1) y=maxfilt(x)   % maximum filter along first non-singleton dimension
        (2) y=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=maxfilt([u v]);      [yu,ku,x0)=maxfilt(u);
                                        yv=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
             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 output

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [y,k,y0]=maxfilt(x,f,n,d,x0)
0002 %MAXFILT find max of an exponentially weighted sliding window  [Y,K,Y0]=(X,F,nn,D,X0)
0003 %
0004 % Usage: (1) y=maxfilt(x)   % maximum filter along first non-singleton dimension
0005 %        (2) y=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=maxfilt([u v]);      [yu,ku,x0)=maxfilt(u);
0008 %                                        yv=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 %             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]=maxfilt(x,1,2*w+1);p=find(((1:m)-k)==w);plot(1:m,x,'-',p-w,x(p-w),'+')
0028 
0029 %      Copyright (C) Mike Brookes 2003
0030 %      Version: $Id: maxfilt.m 4054 2014-01-12 19:11:46Z 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 s=size(x);
0052 if nargin<4 || isempty(d)
0053     d=find(s>1,1);                  % find first non-singleton dimension
0054     if isempty(d)
0055         d=1;
0056     end
0057 end
0058 if nargin>4 && numel(x0)>0          % initial values specified
0059     y=shiftdim(cat(d,x0,x),d-1);    % concatenate x0 and x along d
0060     nx0=size(x0,d);                 % number of values added onto front of data
0061 else                                % dimension specified, d
0062     y=shiftdim(x,d-1);
0063     nx0=0;
0064 end
0065 s=size(y);
0066 s1=s(1);
0067 if nargin<3 || isempty(n)
0068     n0=Inf;
0069 else
0070     n0=max(n,1);
0071 end
0072 if nargin<2 || isempty(f)
0073     f=1;
0074 end
0075 nn=n0;
0076 if nargout>2 % we need to output the tail for next time
0077     if n0<Inf
0078         ny0=min(s1,nn-1);
0079     else
0080         ny0=min(s1,1);
0081     end
0082     sy0=s;
0083     sy0(1)=ny0;
0084     if ny0<=0 || n0==Inf
0085         y0=zeros(sy0);
0086     else
0087         y0=reshape(y(1+s1-ny0:end,:),sy0);
0088         y0=shiftdim(y0,ndims(x)-d+1);
0089     end
0090 end
0091 nn=min(nn,s1);         % no point in having nn>s1
0092 k=repmat((1:s1)',[1 s(2:end)]);
0093 if nn>1
0094     j=1;
0095     j2=1;
0096     while j>0
0097         g=f^j;
0098         m=find(y(j+1:s1,:)<=g*y(1:s1-j,:));
0099         m=m+j*fix((m-1)/(s1-j));
0100         y(m+j)=g*y(m);
0101         k(m+j)=k(m);
0102         j2=j2+j;
0103         j=min(j2,nn-j2);                    % j approximately doubles each iteration
0104     end
0105 end
0106 if nargout==0
0107     if nargin<3
0108         x=shiftdim(x);
0109     else
0110         x=shiftdim(x,d-1);
0111     end
0112     ss=min(prod(s(2:end)),5);   % maximum of 5 plots
0113     plot(1:s1,reshape(y(nx0+1:end,1:ss),s1,ss),'-r',1:s1,reshape(x(:,1:ss),s1,ss),'-b');
0114 else
0115     if nargout>2 && n0==Inf && ny0==1 % if n0==Inf, we need to save the final output
0116         y0=reshape(y(end,:),sy0);
0117         y0=shiftdim(y0,ndims(x)-d+1);
0118     end
0119     if nx0>0                        % pre-data specified, x0
0120         s(1)=s(1)-nx0;
0121         y=shiftdim(reshape(y(nx0+1:end,:),s),ndims(x)-d+1);
0122         k=shiftdim(reshape(k(nx0+1:end,:),s),ndims(x)-d+1)-nx0;
0123     else                            % no pre-data
0124         y=shiftdim(y,ndims(x)-d+1);
0125         k=shiftdim(k,ndims(x)-d+1);
0126     end
0127 end

Generated on Fri 22-Sep-2017 19:37:38 by m2html © 2003