Update of the license file + update of headers
[scilab.git] / scilab / modules / signal_processing / macros / group.sci
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
9
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 funtion 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 //!
33
34
35 //get frequency grid values in [0,.5)
36  
37    fr=(0:.5/npts:.5*(npts-1)/npts);
38  
39 //get the of arguments used to called group
40  
41    [ns,ne]=argn(0);
42  
43 //if the number of arguments is 2 then
44 //the case may be non-cascade
45  
46    hcs=1;
47    if ne==2 then
48  
49 //get the type of h and the variable name
50  
51       h=a1i;
52       ht=type(h);
53  
54 //if ht==1 then h is a vector containing filter coefficients
55  
56       if ht==1 then
57  
58 //make h a rational polynomial
59  
60          hs=maxi(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,
66  
67 //if ht==16 then h is a rational polynomial
68 //(perhaps in cascade form)
69  
70       //-compat ht==15 retained for list/tlist compatibility
71       if ht==15|ht==16 then
72          z=varn(h(3));
73          hcs=maxi(size(h(2)));
74       end,
75  
76 //if the rational polynomial is not in cascade form then
77  
78       if hcs==1 then
79  
80 //if ht==2 then h is a regular polynomial
81  
82          if ht==2 then
83             z=varn(h);
84          end,
85  
86 //get the derivative of h(z)
87  
88          hzd=derivat(h);
89  
90 //get the group delay of h(z)
91  
92          z=poly(0,z);
93          tgz=-z*hzd/h;
94  
95 //evaluate tg
96  
97          rfr=exp(2*%pi*%i*fr);
98          tg=real(freq(tgz(2),tgz(3),rfr));
99  
100 //done with non-cascade calculation of group delay
101  
102       end,
103    end,
104  
105 //re-organize if h is in cascade form
106  
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,
115  
116 //if implementation is in cascade form
117  
118    if ne==5 then
119  
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)
123  
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));
128  
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;
135  
136       aterm=(a1+2*ejw)./(a1+ejw+a2.*emjw);
137       bterm=(b1+2*ejw)./(b1+ejw+b2.*emjw);
138  
139       tgi=real(bterm-aterm);
140       tg=ones(a1i)*tgi;
141  
142 //done with cascade calculation of group delay
143 end
144 endfunction