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).^20001 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 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