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