# v_momfilt

## PURPOSE

V_MOMFILT calculates moments of a signal using a sliding window Y=(X,R,W,M)

## SYNOPSIS

function [y,mm]=v_momfilt(x,r,w,m)

## DESCRIPTION

    V_MOMFILT calculates moments of a signal using a sliding window Y=(X,R,W,M)

     Inputs: x    is the input signal
             r    is a list of moments to calculate
                  (+ve = relative to mean, -ve = relative to zero)
             w    is the window (or just the length of a Hamming window)
                  Note: If the window is asymmetric, you should be aware that it gets
                  flipped in the convolution process
             m    is the sample of w to use as the centre [default=ceil(length(w)/2+0.5)]

             mm   the actual value of m used. Output point y(i) is based on x(i+m-w:i+m-1).

     Example:
      To calculate a running kurtosis using a Hamming window of length 30:
                 y=v_momfilt(x,[2 4],30); k=y(:,2)./y(:,1).^2

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function [y,mm]=v_momfilt(x,r,w,m)
    0002 %V_MOMFILT calculates moments of a signal using a sliding window Y=(X,R,W,M)
    0003 %
    0004 % Inputs: x    is the input signal
    0005 %         r    is a list of moments to calculate
    0006 %              (+ve = relative to mean, -ve = relative to zero)
    0007 %         w    is the window (or just the length of a Hamming window)
    0008 %              Note: If the window is asymmetric, you should be aware that it gets
    0009 %              flipped in the convolution process
    0010 %         m    is the sample of w to use as the centre [default=ceil(length(w)/2+0.5)]
    0011 %
    0012 %         mm   the actual value of m used. Output point y(i) is based on x(i+m-w:i+m-1).
    0013 %
    0014 % Example:
    0015 %  To calculate a running kurtosis using a Hamming window of length 30:
    0016 %             y=v_momfilt(x,[2 4],30); k=y(:,2)./y(:,1).^2
    0017
    0018 %       Copyright (C) Mike Brookes 2007
    0019 %      Version: $Id: v_momfilt.m 10865 2018-09-21 17:22:45Z dmb $
    0020 %
    0021 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0022 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0023 %
    0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0025 %   This program is free software; you can redistribute it and/or modify
    0026 %   it under the terms of the GNU Lesser General Public License as published by
    0027 %   the Free Software Foundation; either version 3 of the License, or
    0028 %   (at your option) any later version.
    0029 %
    0030 %   This program is distributed in the hope that it will be useful,
    0031 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0032 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0033 %   GNU Lesser General Public License for more details.
    0034 %
    0035 %   You can obtain a copy of the GNU Lesser General Public License from
    0036 %   https://www.gnu.org/licenses/ .
    0037 %
    0038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0039 if nargin < 3
    0040    w=hamming(length(x));
    0041 elseif  length(w)==1
    0042    w=hamming(w);
    0043 end
    0044 lw=length(w);
    0045 w=w(:);                 % force to a column vector
    0046 if nargin < 4
    0047    m=(1+lw)/2;
    0048 end
    0049 m=max(round(m),1);
    0050 mm=m;
    0051 r=round(r(:))';             % force integer row vector of moments
    0052
    0053 lx=prod(size(x));
    0054 lxx=lx+m-1;
    0055 xx=zeros(lxx,1);
    0056 xx(1:lx)=x(:);              % extend with zeros so filter() works correctly
    0057
    0058 cw=cumsum(w);
    0059 sw=cw(end);
    0060 y0=repmat(sw,lxx,1);
    0061 lxw=min(lxx,lw);
    0062 y0(1:lxw)=cw(1:lxw);
    0063 y0(lx+1:lx+m-1)=y0(lx+1:lx+m-1)-cw(1:m-1);      % equivalent to y0=filter(w,1,xx^0);
    0064 yd=y0(m:end);
    0065 yd(abs(yd)<eps)=1;
    0066
    0067 nr=length(r);
    0068 wlx=ones(lx,1);
    0069 wlxx=ones(lxx,1);
    0070 y=zeros(lx,nr);
    0071 mr=max(abs(r));                         % max moment to calculate
    0072 mk=zeros(1,mr);                       % list of moments
    0073 mk(-r(r<0))=1;                         % choose the moments we need to calculate
    0074 maxr=max(r);
    0075 if maxr>1
    0076     mk(1:maxr)=1;
    0077 end
    0078 ml=find(mk>0);
    0079 lml=length(ml);
    0080 if lml
    0081     mlx=mk;
    0082     mlx(ml)=1:lml;                        % mapping from moment into ml
    0083     wlml=ones(1,lml);
    0084     xm=filter(w,1,xx(:,wlml).^ml(wlxx,:));    % calculate all the moments
    0085     xm=xm(m:end,:)./yd(:,wlml);                             % remove the useless start values and normalize
    0086 end
    0087 fr=find(r<0);
    0088 if length(fr)
    0089     y(:,fr)=xm(:,mlx(-r(fr)));    % zero-centred moments
    0090 end
    0091 fr=find(r==0);                              % 0'th moment
    0092 if length(fr)
    0093     y(:,fr)=1;
    0094 end
    0095 fr=find(r==1);                              % 1'st moment about mean
    0096 if length(fr)
    0097     y(:,fr)=0;
    0098 end
    0099 fr=find(r==2);                              % 2'nd moment about mean
    0100 if length(fr)
    0101     yfr=xm(:,2)-xm(:,1).^2;
    0102     y(:,fr)=yfr(:,ones(1,length(fr)));
    0103 end
    0104 if maxr>2
    0105     mon=[1 -1];
    0106     bc=[1 -2 1];
    0107     am=zeros(lx,maxr);
    0108     am(:,1)=xm(:,1);                            % copy the mean across
    0109     ml=2:maxr;
    0110     wlml=ones(1,length(ml));
    0111     am(:,2:end)=xm(:,wlml).^ml(ones(lx,1),:);              % calculate powers of the mean
    0112     for i=3:maxr
    0113         bc=conv(bc,mon);                        % calculate binomial coefficients
    0114         fr=find(r==i);
    0115         if length(fr)
    0116             yfr=xm(:,i)+sum(xm(:,i-1:-1:2).*am(:,1:i-2).*bc(wlx,2:i-1),2)+am(:,i)*(bc(i)+bc(i+1));
    0117             y(:,fr)=yfr(:,ones(1,length(fr)));
    0118         end
    0119     end
    0120 end