Group Delay Computation: grpdelay.m

Figure J.11 gives a listing of a matlab program for computing the group delay of an IIR digital filter $ H(z)=B(z)/A(z)$ using the method described in §7.6.6.

In Matlab with the Signal Processing Toolbox installed, (or Octave with the Octave Forge package installed), say 'help grpdelay' for usage documentation, and say 'type grpdelay' to additionally see test, demo, and plotting code. Here, we include only the code relevant to computation of the group delay itself.

Figure J.11: Matlab/Octave function for computing the group delay of a digital filter.

function [gd,w] = grpdelay(b,a,nfft,whole,Fs)

  if (nargin<1 || nargin>5)
    usage("[g,w]=grpdelay(b [, a [, n [,'whole'[,Fs]]]])");
  if nargin<5
    Fs=0; % return w in radians per sample
    if nargin<4, whole='';
    elseif ~isstr(whole)
      Fs = whole;
      whole = '';
    if nargin<3, nfft=512; end
    if nargin<2, a=1; end

  if strcmp(whole,'whole')==0, nfft = 2*nfft; end

  w = 2*pi*[0:nfft-1]/nfft;
  if Fs>0, w = Fs*w/(2*pi); end

  oa = length(a)-1;             % order of a(z)
  oc = oa + length(b)-1;        % order of c(z)
  c = conv(b,fliplr(a));	% c(z) = b(z)*a(1/z)*z^(-oa)
  cr = c.*[0:oc];               % derivative of c wrt 1/z
  num = fft(cr,nfft);
  den = fft(c,nfft);
  minmag = 10*eps;
  polebins = find(abs(den)<minmag);
  for b=polebins
    disp('*** grpdelay: group delay singular! setting to 0')
    num(b) = 0;
    den(b) = 1;
  gd = real(num ./ den) - oa;

  if strcmp(whole,'whole')==0
    ns = nfft/2; % Matlab convention - should be nfft/2 + 1
    gd = gd(1:ns);
    w = w(1:ns);

  w = w'; % Matlab returns column vectors
  gd = gd';

Next Section:
Matlab listing: fold.m
Previous Section:
Parallel SOS to Transfer Function: psos2tf.m