1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - 1988 - C. Bunks
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 function [tg,fr]=group(npts,a1i,a2i,b1i,b2i)
11     //Calculate the group delay of a digital filter
12     //with transfer function h(z).
13     //The filter specification can be in coefficient form,
14     //polynomial form, rational polynomial form, cascade
15     //polynomial form, or in coefficient polynomial form.
16     //  npts :Number of points desired in calculation of group delay
17     //  a1i  :In coefficient, polynomial, rational polynomial, or
18     //       :cascade polynomial form this variable is the transfer
19     //       :function of the filter.  In coefficient polynomial
20     //       :form this is a vector of coefficients (see below).
21     //  a2i  :In coeff poly form this is a vector of coeffs
22     //  b1i  :In coeff poly form this is a vector of coeffs
23     //  b2i  :In coeff poly form this is a vector of coeffs
24     //  tg   :Values of group delay evaluated on the grid fr
25     //  fr   :Grid of frequency values where group delay is evaluated
26     //
27     //In the coefficient polynomial form the tranfer function is
28     //formulated by the following expression:
29     //
30     //       h(z)=prod(a1i+a2i*z+z**2)/prod(b1i+b2i*z+z^2)
31     //
32     //!
35     //get frequency grid values in [0,.5)
37     fr=(0:.5/npts:.5*(npts-1)/npts);
39     //get the of arguments used to called group
41     [ns,ne]=argn(0);
43     //if the number of arguments is 2 then
44     //the case may be non-cascade
46     hcs=1;
47     if ne==2 then
49         //get the type of h and the variable name
51         h=a1i;
52         ht=type(h);
54         //if ht==1 then h is a vector containing filter coefficients
56         if ht==1 then
58             //make h a rational polynomial
60             hs=max(size(h));
61             z=poly(0,"z");
62             h=poly(h,"z","c");
63             h=gtild(h,"d")*(1/z^(hs-1));
64             ht=16;
65         end,
67         //if ht==16 then h is a rational polynomial
68         //(perhaps in cascade form)
70         //-compat ht==15 retained for list/tlist compatibility
71         if ht==15|ht==16 then
72             z=varn(h(3));
73             hcs=max(size(h(2)));
74         end,
76         //if the rational polynomial is not in cascade form then
78         if hcs==1 then
80             //if ht==2 then h is a regular polynomial
82             if ht==2 then
83                 z=varn(h);
84             end,
86             //get the derivative of h(z)
88             hzd=derivat(h);
90             //get the group delay of h(z)
92             z=poly(0,z);
93             tgz=-z*hzd/h;
95             //evaluate tg
97             rfr=exp(2*%pi*%i*fr);
98             tg=real(freq(tgz(2),tgz(3),rfr));
100             //done with non-cascade calculation of group delay
102         end,
103     end,
105     //re-organize if h is in cascade form
107     if hcs>1 then
108         xc=[coeff(h(2)),coeff(h(3))];
109         a2i=xc(1:hcs);
110         a1i=xc(hcs+1:2*hcs);
111         b2i=xc(3*hcs+1:4*hcs);
112         b1i=xc(4*hcs+1:5*hcs);
113         ne=5;
114     end,
116     //if implementation is in cascade form
118     if ne==5 then
120         //a1i,a2i,b1i,b2i are the coefficients of a
121         //second order decomposition of the filter
122         //(i.e. in cascade form, see Deczky)
124         phi=2*%pi*fr;
125         z=poly(0,"z");
126         ejw=freq(z,1,exp(%i*phi));
127         emjw=freq(z,1,exp(-%i*phi));
129         a1=a1i'*ones(phi);
130         b1=b1i'*ones(phi);
131         a2=a2i'*ones(phi);
132         b2=b2i'*ones(phi);
133         ejw=ones(a1i)'*ejw;
134         emjw=ones(a1i)'*emjw;
136         aterm=(a1+2*ejw)./(a1+ejw+a2.*emjw);
137         bterm=(b1+2*ejw)./(b1+ejw+b2.*emjw);
139         tgi=real(bterm-aterm);
140         tg=ones(a1i)*tgi;
142         //done with cascade calculation of group delay
143     end
144 endfunction