V_EWGRPDEL calculates the energy weighted group delay waveform Y=(X,W,M) For each sample, x(n), this routine calculates the energy-weighted average group delay over frequency using a window centred on x(n). This is equal to the delay from the window centre to the centre of gravity of energy in the window. Inputs: x is the input signal w is the window (or just the length of a Hamming window) m is the sample of w to use as the centre [default=(length(w)+1)/2] mm the actual value of m used. Output point y(i) is based on x(i+m-w:i+m-1). If w is odd and m has its default value, then an impulse at x(i) will result in a negative-going zero crossing at y(i). More generally, if the window is symmetrical it will result in a negative-going zero crossing at y(i+m-(w+1)/2).
0001 function [y,mm]=v_ewgrpdel(x,w,m) 0002 %V_EWGRPDEL calculates the energy weighted group delay waveform Y=(X,W,M) 0003 % For each sample, x(n), this routine calculates the energy-weighted average 0004 % group delay over frequency using a window centred on x(n). This is equal to 0005 % the delay from the window centre to the centre of gravity of energy in the window. 0006 % 0007 % Inputs: x is the input signal 0008 % w is the window (or just the length of a Hamming window) 0009 % m is the sample of w to use as the centre [default=(length(w)+1)/2] 0010 % 0011 % mm the actual value of m used. Output point y(i) is based on x(i+m-w:i+m-1). 0012 % 0013 % If w is odd and m has its default value, then an impulse at x(i) will 0014 % result in a negative-going zero crossing at y(i). More generally, if the 0015 % window is symmetrical it will result in a negative-going zero crossing at 0016 % y(i+m-(w+1)/2). 0017 0018 % Example: x=zeros(1000,1); x(100:100:900)=1; v_ewgrpdel(x,99); 0019 0020 % Copyright (C) Mike Brookes 2003 0021 % Version: $Id: v_ewgrpdel.m 10865 2018-09-21 17:22:45Z dmb $ 0022 % 0023 % VOICEBOX is a MATLAB toolbox for speech processing. 0024 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0025 % 0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0027 % This program is free software; you can redistribute it and/or modify 0028 % it under the terms of the GNU General Public License as published by 0029 % the Free Software Foundation; either version 2 of the License, or 0030 % (at your option) any later version. 0031 % 0032 % This program is distributed in the hope that it will be useful, 0033 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0035 % GNU General Public License for more details. 0036 % 0037 % You can obtain a copy of the GNU General Public License from 0038 % http://www.gnu.org/copyleft/gpl.html or by writing to 0039 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0041 0042 if nargin < 2 0043 w=hamming(length(x)); 0044 elseif length(w)==1 0045 w=hamming(w); 0046 end 0047 w=w.^2; 0048 lw=length(w); 0049 if nargin < 3 0050 m=(1+lw)/2; 0051 end 0052 m=max(round(m),1); 0053 mm=m; 0054 wn=w(:).*(m-(1:lw))'; 0055 x2=[x(:); zeros(m-1,1)].^2; 0056 yn=filter(wn,1,x2); 0057 yd=filter(w,1,x2); 0058 yd(yd<eps)=1; 0059 y=yn(m:end)./yd(m:end); 0060 if nargout==0 0061 plot(y); 0062 hold on; 0063 plot(x,'r'); 0064 hold off; 0065 end