scilab/modules/signal_processing/macros/group.sci

   1 function [tg,fr]=group(npts,a1i,a2i,b1i,b2i)
   2 //Calculate the group delay of a digital filter
   3 //with transfer function h(z).
   4 //The filter specification can be in coefficient form,
   5 //polynomial form, rational polynomial form, cascade
   6 //polynomial form, or in coefficient polynomial form.
   7 //  npts :Number of points desired in calculation of group delay
   8 //  a1i  :In coefficient, polynomial, rational polynomial, or
   9 //       :cascade polynomial form this variable is the transfer
  10 //       :function of the filter.  In coefficient polynomial
  11 //       :form this is a vector of coefficients (see below).
  12 //  a2i  :In coeff poly form this is a vector of coeffs
  13 //  b1i  :In coeff poly form this is a vector of coeffs
  14 //  b2i  :In coeff poly form this is a vector of coeffs
  15 //  tg   :Values of group delay evaluated on the grid fr
  16 //  fr   :Grid of frequency values where group delay is evaluated
  17 //
  18 //In the coefficient polynomial form the tranfer funtion is
  19 //formulated by the following expression:
  20 //
  21 //       h(z)=prod(a1i+a2i*z+z**2)/prod(b1i+b2i*z+z^2)
  22 //
  23 //!
  24 //author: C. Bunks  date: 2 March 1988
  25 // Copyright INRIA
  26
  27 //get frequency grid values in [0,.5)
  28
  29    fr=(0:.5/npts:.5*(npts-1)/npts);
  30
  31 //get the of arguments used to called group
  32
  33    [ns,ne]=argn(0);
  34
  35 //if the number of arguments is 2 then
  36 //the case may be non-cascade
  37
  38    hcs=1;
  39    if ne==2 then
  40
  41 //get the type of h and the variable name
  42
  43       h=a1i;
  44       ht=type(h);
  45
  46 //if ht==1 then h is a vector containing filter coefficients
  47
  48       if ht==1 then
  49
  50 //make h a rational polynomial
  51
  52          hs=maxi(size(h));
  53          z=poly(0,'z');
  54          h=poly(h,'z','c');
  55          h=gtild(h,'d')*(1/z^(hs-1));
  56          ht=16;
  57       end,
  58
  59 //if ht==16 then h is a rational polynomial
  60 //(perhaps in cascade form)
  61
  62       //-compat ht==15 retained for list/tlist compatibility
  63       if ht==15|ht==16 then
  64          z=varn(h(3));
  65          hcs=maxi(size(h(2)));
  66       end,
  67
  68 //if the rational polynomial is not in cascade form then
  69
  70       if hcs==1 then
  71
  72 //if ht==2 then h is a regular polynomial
  73
  74          if ht==2 then
  75             z=varn(h);
  76          end,
  77
  78 //get the derivative of h(z)
  79
  80          hzd=derivat(h);
  81
  82 //get the group delay of h(z)
  83
  84          z=poly(0,z);
  85          tgz=-z*hzd/h;
  86
  87 //evaluate tg
  88
  89          rfr=exp(2*%pi*%i*fr);
  90          tg=real(freq(tgz(2),tgz(3),rfr));
  91
  92 //done with non-cascade calculation of group delay
  93
  94       end,
  95    end,
  96
  97 //re-organize if h is in cascade form
  98
  99    if hcs>1 then
 100       xc=[coeff(h(2)),coeff(h(3))];
 101       a2i=xc(1:hcs);
 102       a1i=xc(hcs+1:2*hcs);
 103       b2i=xc(3*hcs+1:4*hcs);
 104       b1i=xc(4*hcs+1:5*hcs);
 105       ne=5;
 106    end,
 107
 108 //if implementation is in cascade form
 109
 110    if ne==5 then
 111
 112 //a1i,a2i,b1i,b2i are the coefficients of a
 113 //second order decomposition of the filter
 114 //(i.e. in cascade form, see Deczky)
 115
 116       phi=2*%pi*fr;
 117       z=poly(0,'z');
 118       ejw=freq(z,1,exp(%i*phi));
 119       emjw=freq(z,1,exp(-%i*phi));
 120
 121       a1=a1i'*ones(phi);
 122       b1=b1i'*ones(phi);
 123       a2=a2i'*ones(phi);
 124       b2=b2i'*ones(phi);
 125       ejw=ones(a1i)'*ejw;
 126       emjw=ones(a1i)'*emjw;
 127
 128       aterm=(a1+2*ejw)./(a1+ejw+a2.*emjw);
 129       bterm=(b1+2*ejw)./(b1+ejw+b2.*emjw);
 130
 131       tgi=real(bterm-aterm);
 132       tg=ones(a1i)*tgi;
 133
 134 //done with cascade calculation of group delay
 135 end
 136 endfunction