← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to macros/signal/group.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/signal/group.sci
 
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
 
17
//
 
18
//In the coefficient polynomial form the tranfer funtion is
 
19
//formulated by the following expression:
 
20
//
 
21
//       h(z)=prod(a1i+a2i*z+z**2)/prod(b1i+b2i*z+z^2)
 
22
//
 
23
//!
 
24
//author: C. Bunks  date: 2 March 1988
 
25
// Copyright INRIA
 
26
 
 
27
//get frequency grid values in [0,.5)
 
28
 
 
29
   fr=(0:.5/npts:.5*(npts-1)/npts);
 
30
 
 
31
//get the of arguments used to called group
 
32
 
 
33
   [ns,ne]=argn(0);
 
34
 
 
35
//if the number of arguments is 2 then
 
36
//the case may be non-cascade
 
37
 
 
38
   hcs=1;
 
39
   if ne==2 then
 
40
 
 
41
//get the type of h and the variable name
 
42
 
 
43
      h=a1i;
 
44
      ht=type(h);
 
45
 
 
46
//if ht==1 then h is a vector containing filter coefficients
 
47
 
 
48
      if ht==1 then
 
49
 
 
50
//make h a rational polynomial
 
51
 
 
52
         hs=maxi(size(h));
 
53
         z=poly(0,'z');
 
54
         h=poly(h,'z','c');
 
55
         h=gtild(h,'d')*(1/z^(hs-1));
 
56
         ht=16;
 
57
      end,
 
58
 
 
59
//if ht==16 then h is a rational polynomial
 
60
//(perhaps in cascade form)
 
61
 
 
62
      //-compat ht==15 retained for list/tlist compatibility
 
63
      if ht==15|ht==16 then
 
64
         z=varn(h(3));
 
65
         hcs=maxi(size(h(2)));
 
66
      end,
 
67
 
 
68
//if the rational polynomial is not in cascade form then
 
69
 
 
70
      if hcs==1 then
 
71
 
 
72
//if ht==2 then h is a regular polynomial
 
73
 
 
74
         if ht==2 then
 
75
            z=varn(h);
 
76
         end,
 
77
 
 
78
//get the derivative of h(z)
 
79
 
 
80
         hzd=derivat(h);
 
81
 
 
82
//get the group delay of h(z)
 
83
 
 
84
         z=poly(0,z);
 
85
         tgz=-z*hzd/h;
 
86
 
 
87
//evaluate tg
 
88
 
 
89
         rfr=exp(2*%pi*%i*fr);
 
90
         tg=real(freq(tgz(2),tgz(3),rfr));
 
91
 
 
92
//done with non-cascade calculation of group delay
 
93
 
 
94
      end,
 
95
   end,
 
96
 
 
97
//re-organize if h is in cascade form
 
98
 
 
99
   if hcs>1 then
 
100
      xc=[coeff(h(2)),coeff(h(3))];
 
101
      a2i=xc(1:hcs);
 
102
      a1i=xc(hcs+1:2*hcs);
 
103
      b2i=xc(3*hcs+1:4*hcs);
 
104
      b1i=xc(4*hcs+1:5*hcs);
 
105
      ne=5;
 
106
   end,
 
107
 
 
108
//if implementation is in cascade form
 
109
 
 
110
   if ne==5 then
 
111
 
 
112
//a1i,a2i,b1i,b2i are the coefficients of a
 
113
//second order decomposition of the filter
 
114
//(i.e. in cascade form, see Deczky)
 
115
 
 
116
      phi=2*%pi*fr;
 
117
      z=poly(0,'z');
 
118
      ejw=freq(z,1,exp(%i*phi));
 
119
      emjw=freq(z,1,exp(-%i*phi));
 
120
 
 
121
      a1=a1i'*ones(phi);
 
122
      b1=b1i'*ones(phi);
 
123
      a2=a2i'*ones(phi);
 
124
      b2=b2i'*ones(phi);
 
125
      ejw=ones(a1i)'*ejw;
 
126
      emjw=ones(a1i)'*emjw;
 
127
 
 
128
      aterm=(a1+2*ejw)./(a1+ejw+a2.*emjw);
 
129
      bterm=(b1+2*ejw)./(b1+ejw+b2.*emjw);
 
130
 
 
131
      tgi=real(bterm-aterm);
 
132
      tg=ones(a1i)*tgi;
 
133
 
 
134
//done with cascade calculation of group delay
 
135
end