RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
bilinear.cpp
Go to the documentation of this file.
1 //
2 // File = bilinear.cpp
3 //
4 
5 #include "bilinear.h"
6 #include "complex.h"
7 #include "misdefs.h"
8 #include "typedefs.h"
9 #include <fstream>
10 #include <iostream>
11 #include <stdlib.h>
12 #ifdef _DEBUG
13 extern std::ofstream DebugFile;
14 #endif
15 
16 //====================================================
17 //
18 //----------------------------------------------------
20 BilinearTransf(FilterTransFunc* analog_filter, double sampling_interval)
21 
22 {
23  int max_poles;
24  int num_poles, num_zeros;
25  int j, m, n;
26  int max_coeff, num_numer_coeff;
27  double h_const, h_sub_zero, denom_mu_zero;
28  complex *pole, *zero;
29  complex* mu;
30  complex alpha;
31  complex beta;
32  complex gamma;
33  complex delta;
34  complex eta;
35  complex work;
36  complex c_two;
37  double *a, *b;
38  IirFilterDesign* iir_filter;
39 
40  pole = analog_filter->GetPoles(&num_poles);
41  zero = analog_filter->GetZeros(&num_zeros);
42  h_sub_zero = analog_filter->GetHSubZero();
43 #ifdef _DEBUG
44  DebugFile << "num analog poles = " << num_poles << std::endl;
45  DebugFile << "num analog zeros = " << num_zeros << std::endl;
46 #endif
47 
48  if (num_poles > num_zeros) {
49  max_poles = num_poles;
50  } else {
51  max_poles = num_zeros;
52  }
53 
54  //--------------------------------------------
55  // allocate and initialize working storage
56 
57  mu = new complex[max_poles + 1];
58  a = new double[max_poles + 1];
59  b = new double[max_poles + 1];
60 
61  for (j = 0; j <= max_poles; j++) {
62  mu[j] = complex(0.0, 0.0);
63  a[j] = 0.0;
64  b[j] = 0.0;
65  }
66 
67  //-------------------------------------------
68  // compute constant gain factor
69 
70  h_const = 1.0;
71  work = complex(1.0, 0.0);
72  c_two = complex(2.0, 0.0);
73 
74  for (n = 1; n <= num_poles; n++) {
75  work = work * (c_two - (sampling_interval * pole[n]));
76 #ifdef _DEBUG
77  DebugFile << "work = " << work << std::endl;
78 #endif
79  h_const = h_const * sampling_interval;
80  }
81 #ifdef _DEBUG
82  DebugFile << "T**2 = " << h_const << std::endl;
83 #endif
84  h_const = h_sub_zero * h_const / real(work);
85 
86 #ifdef _DEBUG
87  DebugFile << "in BilinearTransf, h_const = " << h_const << std::endl;
88  DebugFile << "work = " << work << std::endl;
89 #endif
90 
91  //--------------------------------------------------
92  // compute denominator coefficients
93 
94  mu[0] = complex(1.0, 0.0);
95 
96  for (n = 1; n <= num_poles; n++) {
97 #ifdef _DEBUG
98  DebugFile << "in BilinearTransf, pole [" << n << "] = " << pole[n]
99  << std::endl;
100 #endif
101 
102  gamma = complex((2.0 / sampling_interval), 0.0) - pole[n];
103  delta = complex((-2.0 / sampling_interval), 0.0) - pole[n];
104 
105 #ifdef _DEBUG
106  DebugFile << "gamma = " << gamma << std::endl;
107  DebugFile << "delta = " << delta << std::endl;
108 #endif
109 
110  for (j = n; j >= 1; j--) {
111  mu[j] = gamma * mu[j] + (delta * mu[j - 1]);
112  }
113  mu[0] = gamma * mu[0];
114  }
115 #ifdef _DEBUG
116  DebugFile << "for denom, mu[0] = " << mu[0] << std::endl;
117 #endif
118  denom_mu_zero = real(mu[0]);
119  for (j = 1; j <= num_poles; j++) {
120  a[j] = -1.0 * real(mu[j]) / denom_mu_zero;
121 #ifdef _DEBUG
122  DebugFile << "a[" << j << "] = " << a[j] << std::endl;
123  DebugFile << "imag(mu[" << j << "]) = " << imag(mu[j]) << std::endl;
124 #endif
125  }
126  //-----------------------------------------------------
127  // compute numerator coeffcients
128 
129  mu[0] = complex(1.0, 0.0);
130  for (n = 1; n <= max_poles; n++) {
131  mu[n] = complex(0.0, 0.0);
132  }
133 
134  max_coeff = 0;
135 
136  //- - - - - - - - - - - - - - - - - - - - -
137  // compute (1+z**(-1)) ** (N-M)
138 
139  for (m = 1; m <= (num_poles - num_zeros); m++) {
140  max_coeff++;
141  for (j = max_coeff; j >= 1; j--) {
142  mu[j] = mu[j] + mu[j - 1];
143  }
144  }
145  for (m = 1; m <= num_zeros; m++) {
146  max_coeff++;
147 #ifdef _DEBUG
148  DebugFile << "zero[" << m << "] = " << zero[m] << std::endl;
149 #endif
150  alpha = complex((2.0 / sampling_interval), 0.0) - zero[m];
151  beta = complex((-2.0 / sampling_interval), 0.0) - zero[m];
152 
153  for (j = max_coeff; j >= 1; j--) {
154  mu[j] = alpha * mu[j] + (beta * mu[j - 1]);
155  }
156  mu[0] = alpha * mu[0];
157  }
158  num_numer_coeff = max_coeff + 1;
159  for (j = 0; j < num_numer_coeff; j++) {
160  b[j] = h_sub_zero * real(mu[j]) / denom_mu_zero;
161 #ifdef _DEBUG
162  DebugFile << "b[" << j << "] = " << b[j] << std::endl;
163  DebugFile << "imag(mu[" << j << "]) = " << imag(mu[j]) << std::endl;
164 #endif
165  }
166 
167  delete[] mu;
168  iir_filter = new IirFilterDesign(num_numer_coeff, num_poles, b, a);
169  iir_filter->SetSamplingInterval(sampling_interval);
170  return (iir_filter);
171 }
complex * GetPoles(int *num_poles)
Definition: filtfunc.cpp:58
complex * GetZeros(int *num_zeros)
Definition: filtfunc.cpp:238
void SetSamplingInterval(double sampling_interval)
Definition: iir_dsgn.cpp:241
std::ofstream DebugFile
IirFilterDesign * BilinearTransf(FilterTransFunc *analog_filter, double sampling_interval)
Definition: bilinear.cpp:20
float GetHSubZero(void)
Definition: filtfunc.cpp:253
double imag(const complex &_z)
Definition: complex.h:103
double real(const complex &_z)
Definition: complex.h:98