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
18 //In the coefficient polynomial form the tranfer funtion is
19 //formulated by the following expression:
21 // h(z)=prod(a1i+a2i*z+z**2)/prod(b1i+b2i*z+z^2)
24 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
25 // Copyright (C) INRIA - 1988 - C. Bunks
27 // This file must be used under the terms of the CeCILL.
28 // This source file is licensed as described in the file COPYING, which
29 // you should have received as part of this distribution. The terms
30 // are also available at
31 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
33 //get frequency grid values in [0,.5)
35 fr=(0:.5/npts:.5*(npts-1)/npts);
37 //get the of arguments used to called group
41 //if the number of arguments is 2 then
42 //the case may be non-cascade
47 //get the type of h and the variable name
52 //if ht==1 then h is a vector containing filter coefficients
56 //make h a rational polynomial
61 h=gtild(h,'d')*(1/z^(hs-1));
65 //if ht==16 then h is a rational polynomial
66 //(perhaps in cascade form)
68 //-compat ht==15 retained for list/tlist compatibility
74 //if the rational polynomial is not in cascade form then
78 //if ht==2 then h is a regular polynomial
84 //get the derivative of h(z)
88 //get the group delay of h(z)
96 tg=real(freq(tgz(2),tgz(3),rfr));
98 //done with non-cascade calculation of group delay
103 //re-organize if h is in cascade form
106 xc=[coeff(h(2)),coeff(h(3))];
109 b2i=xc(3*hcs+1:4*hcs);
110 b1i=xc(4*hcs+1:5*hcs);
114 //if implementation is in cascade form
118 //a1i,a2i,b1i,b2i are the coefficients of a
119 //second order decomposition of the filter
120 //(i.e. in cascade form, see Deczky)
124 ejw=freq(z,1,exp(%i*phi));
125 emjw=freq(z,1,exp(-%i*phi));
132 emjw=ones(a1i)'*emjw;
134 aterm=(a1+2*ejw)./(a1+ejw+a2.*emjw);
135 bterm=(b1+2*ejw)./(b1+ejw+b2.*emjw);
137 tgi=real(bterm-aterm);
140 //done with cascade calculation of group delay