Home > voicebox > momfilt.m

momfilt

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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=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]=momfilt(x,r,w,m)
0002 %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=momfilt(x,[2 4],30); k=y(:,2)./y(:,1).^2
0017 
0018 %       Copyright (C) Mike Brookes 2007
0019 %      Version: $Id: momfilt.m 713 2011-10-16 14:45:43Z 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 General Public License as published by
0027 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0034 %
0035 %   You can obtain a copy of the GNU General Public License from
0036 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0037 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003