octave_packages/signal-1.1.3/grpdelay.m

   1 ## Copyright (C) 2000 Paul Kienzle  <pkienzle@users.sf.net>
   2 ## Copyright (C) 2004 Julius O. Smith III <jos@ccrma.stanford.edu>
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## Compute the group delay of a filter.
  18 ##
  19 ## [g, w] = grpdelay(b)
  20 ##   returns the group delay g of the FIR filter with coefficients b.
  21 ##   The response is evaluated at 512 angular frequencies between 0 and
  22 ##   pi. w is a vector containing the 512 frequencies.
  23 ##   The group delay is in units of samples.  It can be converted
  24 ##   to seconds by multiplying by the sampling period (or dividing by
  25 ##   the sampling rate fs).
  26 ##
  27 ## [g, w] = grpdelay(b,a)
  28 ##   returns the group delay of the rational IIR filter whose numerator
  29 ##   has coefficients b and denominator coefficients a.
  30 ##
  31 ## [g, w] = grpdelay(b,a,n)
  32 ##   returns the group delay evaluated at n angular frequencies.  For fastest
  33 ##   computation n should factor into a small number of small primes.
  34 ##
  35 ## [g, w] = grpdelay(b,a,n,'whole')
  36 ##   evaluates the group delay at n frequencies between 0 and 2*pi.
  37 ##
  38 ## [g, f] = grpdelay(b,a,n,Fs)
  39 ##   evaluates the group delay at n frequencies between 0 and Fs/2.
  40 ##
  41 ## [g, f] = grpdelay(b,a,n,'whole',Fs)
  42 ##   evaluates the group delay at n frequencies between 0 and Fs.
  43 ##
  44 ## [g, w] = grpdelay(b,a,w)
  45 ##   evaluates the group delay at frequencies w (radians per sample).
  46 ##
  47 ## [g, f] = grpdelay(b,a,f,Fs)
  48 ##   evaluates the group delay at frequencies f (in Hz).
  49 ##
  50 ## grpdelay(...)
  51 ##   plots the group delay vs. frequency.
  52 ##
  53 ## If the denominator of the computation becomes too small, the group delay
  54 ## is set to zero.  (The group delay approaches infinity when
  55 ## there are poles or zeros very close to the unit circle in the z plane.)
  56 ##
  57 ## Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}],  is the rate of change of
  58 ## phase with respect to frequency.  It can be computed as:
  59 ##
  60 ##               d/dw H(e^-jw)
  61 ##        g(w) = -------------
  62 ##                 H(e^-jw)
  63 ##
  64 ## where
  65 ##         H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).
  66 ##
  67 ## By the quotient rule,
  68 ##                    A(z) d/dw B(z) - B(z) d/dw A(z)
  69 ##        d/dw H(z) = -------------------------------
  70 ##                               A(z) A(z)
  71 ## Substituting into the expression above yields:
  72 ##                A dB - B dA
  73 ##        g(w) =  ----------- = dB/B - dA/A
  74 ##                    A B
  75 ##
  76 ## Note that,
  77 ##        d/dw B(e^-jw) = sum(k b_k e^-jwk)
  78 ##        d/dw A(e^-jw) = sum(k a_k e^-jwk)
  79 ## which is just the FFT of the coefficients multiplied by a ramp.
  80 ##
  81 ## As a further optimization when nfft>>length(a), the IIR filter (b,a)
  82 ## is converted to the FIR filter conv(b,fliplr(conj(a))).
  83 ## For further details, see
  84 ## http://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html
  85
  86 function [gd,w] = grpdelay(b,a=1,nfft=512,whole,Fs)
  87
  88   if (nargin<1 || nargin>5)
  89     print_usage;
  90   end
  91   HzFlag=0;
  92   if length(nfft)>1
  93     if nargin>4
  94       print_usage();
  95     elseif nargin>3  % grpdelay(B,A,F,Fs)
  96       Fs = whole;
  97       HzFlag=1;
  98     else % grpdelay(B,A,W)
  99       Fs = 1;
 100     end
 101     w = 2*pi*nfft/Fs;
 102     nfft = length(w)*2;
 103     whole = '';
 104   else
 105     if nargin<5
 106       Fs=1; % return w in radians per sample
 107       if nargin<4, whole='';
 108       elseif ~ischar(whole)
 109         Fs = whole;
 110         HzFlag=1;
 111         whole = '';
 112       end
 113       if nargin<3, nfft=512; end
 114       if nargin<2, a=1; end
 115     else
 116       HzFlag=1;
 117     end
 118
 119     if isempty(nfft), nfft = 512; end
 120     if ~strcmp(whole,'whole'), nfft = 2*nfft; end
 121     w = Fs*[0:nfft-1]/nfft;
 122   end
 123
 124   if ~HzFlag, w = w * 2*pi; end
 125
 126   oa = length(a)-1;             % order of a(z)
 127   if oa<0, a=1; oa=0; end       % a can be []
 128   ob = length(b)-1;             % order of b(z)
 129   if ob<0, b=1; ob=0; end       % b can be [] as well
 130   oc = oa + ob;                 % order of c(z)
 131
 132   c = conv(b,fliplr(conj(a)));  % c(z) = b(z)*conj(a)(1/z)*z^(-oa)
 133   cr = c.*[0:oc];               % cr(z) = derivative of c wrt 1/z
 134   num = fft(cr,nfft);
 135   den = fft(c,nfft);
 136   minmag = 10*eps;
 137   polebins = find(abs(den)<minmag);
 138   for b=polebins
 139     warning('grpdelay: setting group delay to 0 at singularity');
 140     num(b) = 0;
 141     den(b) = 1;
 142     % try to preserve angle:
 143     % db = den(b);
 144     % den(b) = minmag*abs(num(b))*exp(j*atan2(imag(db),real(db)));
 145     % warning(sprintf('grpdelay: den(b) changed from %f to %f',db,den(b)));
 146   end
 147   gd = real(num ./ den) - oa;
 148
 149   if strcmp(whole,'whole')==0
 150     ns = nfft/2; % Matlab convention ... should be nfft/2 + 1
 151     gd = gd(1:ns);
 152     w = w(1:ns);
 153   else
 154     ns = nfft; % used in plot below
 155   end
 156
 157   % compatibility
 158   gd = gd(:); w = w(:);
 159
 160   if nargout==0
 161     unwind_protect
 162       grid('on'); % grid() should return its previous state
 163       if HzFlag
 164         funits = 'Hz';
 165       else
 166         funits = 'radian/sample';
 167       end
 168       xlabel(['Frequency (', funits, ')']);
 169       ylabel('Group delay (samples)');
 170       plot(w(1:ns), gd(1:ns), ';;');
 171     unwind_protect_cleanup
 172       grid('on');
 173     end_unwind_protect
 174   end
 175 end
 176
 177 % ------------------------ DEMOS -----------------------
 178
 179 %!demo % 1
 180 %! %--------------------------------------------------------------
 181 %! % From Oppenheim and Schafer, a single zero of radius r=0.9 at
 182 %! % angle pi should have a group delay of about -9 at 1 and 1/2
 183 %! % at zero and 2*pi.
 184 %! %--------------------------------------------------------------
 185 %! title ('Zero at z = -0.9');
 186 %! grpdelay([1 0.9],[],512,'whole',1);
 187 %! hold on;
 188 %! xlabel('Normalized Frequency (cycles/sample)');
 189 %! stem([0, 0.5, 1],[0.5, -9, 0.5],'*b;target;');
 190 %! hold off;
 191 %!
 192 %!demo % 2
 193 %! %--------------------------------------------------------------
 194 %! % confirm the group delays approximately meet the targets
 195 %! % don't worry that it is not exact, as I have not entered
 196 %! % the exact targets.
 197 %! %--------------------------------------------------------------
 198 %! b = poly([1/0.9*exp(1i*pi*0.2), 0.9*exp(1i*pi*0.6)]);
 199 %! a = poly([0.9*exp(-1i*pi*0.6), 1/0.9*exp(-1i*pi*0.2)]);
 200 %! title ('Two Zeros and Two Poles');
 201 %! grpdelay(b,a,512,'whole',1);
 202 %! hold on;
 203 %! xlabel('Normalized Frequency (cycles/sample)');
 204 %! stem([0.1, 0.3, 0.7, 0.9], [9, -9, 9, -9],'*b;target;');
 205 %! hold off;
 206
 207 %!demo % 3
 208 %! %--------------------------------------------------------------
 209 %! % fir lowpass order 40 with cutoff at w=0.3 and details of
 210 %! % the transition band [.3, .5]
 211 %! %--------------------------------------------------------------
 212 %! subplot(211);
 213 %! Fs = 8000;     % sampling rate
 214 %! Fc = 0.3*Fs/2; % lowpass cut-off frequency
 215 %! nb = 40;
 216 %! b = fir1(nb,2*Fc/Fs); % matlab freq normalization: 1=Fs/2
 217 %! [H,f] = freqz(b,1,[],1);
 218 %! [gd,f] = grpdelay(b,1,[],1);
 219 %! title(sprintf('b = fir1(%d,2*%d/%d);',nb,Fc,Fs));
 220 %! xlabel('Normalized Frequency (cycles/sample)');
 221 %! ylabel('Amplitude Response (dB)');
 222 %! grid('on');
 223 %! plot(f,20*log10(abs(H)));
 224 %! subplot(212);
 225 %! del = nb/2; % should equal this
 226 %! title(sprintf('Group Delay in Pass-Band (Expect %d samples)',del));
 227 %! ylabel('Group Delay (samples)');
 228 %! axis([0, 0.2, del-1, del+1]);
 229 %! plot(f,gd);
 230 %! axis();
 231
 232 %!demo % 4
 233 %! %--------------------------------------------------------------
 234 %! % IIR bandstop filter has delays at [1000, 3000]
 235 %! %--------------------------------------------------------------
 236 %! Fs = 8000;
 237 %! [b, a] = cheby1(3, 3, 2*[1000, 3000]/Fs, 'stop');
 238 %! [H,f] = freqz(b,a,[],Fs);
 239 %! [gd,f] = grpdelay(b,a,[],Fs);
 240 %! subplot(211);
 241 %! title('[b,a] = cheby1(3, 3, 2*[1000, 3000]/Fs, \'stop\');');
 242 %! xlabel('Frequency (Hz)');
 243 %! ylabel('Amplitude Response');
 244 %! grid('on');
 245 %! plot(f,abs(H));
 246 %! subplot(212);
 247 %! title('[gd,f] = grpdelay(b,a,[],Fs);');
 248 %! ylabel('Group Delay (samples)');
 249 %! plot(f,gd);
 250
 251
 252 % ------------------------ TESTS -----------------------
 253
 254 %!test % 00
 255 %! [gd1,w] = grpdelay([0,1]);
 256 %! [gd2,w] = grpdelay([0,1],1);
 257 %! assert(gd1,gd2,10*eps);
 258
 259 %!test % 0A
 260 %! [gd,w] = grpdelay([0,1],1,4);
 261 %! assert(gd,[1;1;1;1]);
 262 %! assert(w,pi/4*[0:3]',10*eps);
 263
 264 %!test % 0B
 265 %! [gd,w] = grpdelay([0,1],1,4,'whole');
 266 %! assert(gd,[1;1;1;1]);
 267 %! assert(w,pi/2*[0:3]',10*eps);
 268
 269 %!test % 0C
 270 %! [gd,f] = grpdelay([0,1],1,4,0.5);
 271 %! assert(gd,[1;1;1;1]);
 272 %! assert(f,1/16*[0:3]',10*eps);
 273
 274 %!test % 0D
 275 %! [gd,w] = grpdelay([0,1],1,4,'whole',1);
 276 %! assert(gd,[1;1;1;1]);
 277 %! assert(w,1/4*[0:3]',10*eps);
 278
 279 %!test % 0E
 280 %! [gd,f] = grpdelay([1 -0.9j],[],4,'whole',1);
 281 %! gd0 = 0.447513812154696; gdm1 =0.473684210526316;
 282 %! assert(gd,[gd0;-9;gd0;gdm1],20*eps);
 283 %! assert(f,1/4*[0:3]',10*eps);
 284
 285 %!test % 1A:
 286 %! gd= grpdelay(1,[1,.9],2*pi*[0,0.125,0.25,0.375]);
 287 %! assert(gd, [-0.47368;-0.46918;-0.44751;-0.32316],1e-5);
 288
 289 %!test % 1B:
 290 %! gd= grpdelay(1,[1,.9],[0,0.125,0.25,0.375],1);
 291 %! assert(gd, [-0.47368;-0.46918;-0.44751;-0.32316],1e-5);
 292
 293 %!test % 2:
 294 %! gd = grpdelay([1,2],[1,0.5,.9],4);
 295 %! assert(gd,[-0.29167;-0.24218;0.53077;0.40658],1e-5);
 296
 297 %!test % 3
 298 %! b1=[1,2];a1f=[0.25,0.5,1];a1=fliplr(a1f);
 299 %! % gd1=grpdelay(b1,a1,4);
 300 %! gd=grpdelay(conv(b1,a1f),1,4)-2;
 301 %! assert(gd, [0.095238;0.239175;0.953846;1.759360],1e-5);
 302
 303 %!test % 4
 304 %! Fs = 8000;
 305 %! [b, a] = cheby1(3, 3, 2*[1000, 3000]/Fs, 'stop');
 306 %! [h, w] = grpdelay(b, a, 256, 'half', Fs);
 307 %! [h2, w2] = grpdelay(b, a, 512, 'whole', Fs);
 308 %! assert (size(h), size(w));
 309 %! assert (length(h), 256);
 310 %! assert (size(h2), size(w2));
 311 %! assert (length(h2), 512);
 312 %! assert (h, h2(1:256));
 313 %! assert (w, w2(1:256));
 314
 315 %!test % 5
 316 %! a = [1 0 0.9];
 317 %! b = [0.9 0 1];
 318 %! [dh, wf] = grpdelay(b, a, 512, 'whole');
 319 %! [da, wa] = grpdelay(1, a, 512, 'whole');
 320 %! [db, wb] = grpdelay(b, 1, 512, 'whole');
 321 %! assert(dh,db+da,1e-5);