# v_maxfilt

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

    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 output

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_activlev V_ACTIVLEV Measure active speech level as in ITU-T P.56 [LEV,AF,FSO]=(sp,FS,MODE)
- v_dypsa V_DYPSA Derive glottal closure instances from speech [gci,goi] = (s,fs)

## SOURCE CODE

    0001 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 Lesser General Public License as published by
    0041 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0048 %
    0049 %   You can obtain a copy of the GNU Lesser General Public License from
    0050 %   https://www.gnu.org/licenses/ .
    0051 %
    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