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
|
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 funtion is
//formulated by the following expression:
//
// h(z)=prod(a1i+a2i*z+z**2)/prod(b1i+b2i*z+z^2)
//
//!
//author: C. Bunks date: 2 March 1988
// Copyright INRIA
//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 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 ht==1 then h is a vector containing filter coefficients
if ht==1 then
//make h a rational polynomial
hs=maxi(size(h));
z=poly(0,'z');
h=poly(h,'z','c');
h=gtild(h,'d')*(1/z^(hs-1));
ht=16;
end,
//if ht==16 then h is a rational polynomial
//(perhaps in cascade form)
//-compat ht==15 retained for list/tlist compatibility
if ht==15|ht==16 then
z=varn(h(3));
hcs=maxi(size(h(2)));
end,
//if the rational polynomial is not in cascade form then
if hcs==1 then
//if ht==2 then h is a regular polynomial
if ht==2 then
z=varn(h);
end,
//get the derivative of h(z)
hzd=derivat(h);
//get the group delay of h(z)
z=poly(0,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,
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,
//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
|
