RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
elipfunc.cpp
Go to the documentation of this file.
1 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 //
3 // File = elipfunc.cpp
4 //
5 // Elliptical Filter Function
6 //
7 
8 #include "elipfunc.h"
9 #include "complex.h"
10 #include "misdefs.h"
11 #include <math.h>
12 #include <stdlib.h>
13 #ifdef _DEBUG
14 extern std::ofstream DebugFile;
15 #endif
16 
17 //======================================================
18 // constructor
19 
20 EllipticalTransFunc::EllipticalTransFunc(int order, double passband_ripple,
21  double stopband_ripple,
22  double passband_edge,
23  double stopband_edge,
24  int upper_summation_limit)
25  : FilterTransFunc(order)
26 {
27  int m;
28  int min_order;
29  double u, u4, x;
30  double modular_const;
31  double selec_factor;
32  double double_temp;
33  double discrim_factor;
34  double term, sum;
35  double numer, denom;
36  double vv, ww, xx, yy;
37  double p_sub_zero;
38  double mu;
39  int i, i_mirror, r;
40  double aa, bb, cc;
41  complex cmplx_work, work2;
42 
43  //------------------------------------
44  // Check viability of parameter set
45 
46  selec_factor = passband_edge / stopband_edge;
47 
48  x = pow((1. - selec_factor * selec_factor), 0.25);
49  u = (1.0 - x) / (2.0 * (1 + x));
50  u4 = u * u * u * u;
51 
52  // compute:
53  // modular_const = u + 2u**5 + 15u**9 + 150u**13
54 
55  modular_const =
56  u * (1. + (2 * u * u4 * (1. + (7.5 * u4 * (1 + (10. * u4))))));
57 
58  discrim_factor = (pow(10.0, stopband_ripple / 10.0) - 1.0) /
59  (pow(10.0, passband_ripple / 10.0) - 1.0);
60 
61  min_order =
62  (int)ceil(log10(16.0 * discrim_factor) / log10(1.0 / modular_const));
63 
64  if (order < min_order) {
65  std::cout << "Fatal error -- minimum order of " << min_order << " required"
66  << std::endl;
67  exit(1);
68  }
69 
70  //------------------------------------------------------
71  // compute transfer function
72 
73  Num_Prototype_Poles = order;
74  Prototype_Pole_Locs = new complex[order + 1];
75 
76  if (order % 2) // order is odd
77  {
78  Num_Prototype_Zeros = order - 1;
79  } else // order is even
80  {
81  Num_Prototype_Zeros = order;
82  }
84 
85  //--------------------------------------------------
86  // step 7 Algorith 5.4
87 
88  r = (order - (order % 2)) / 2;
89  Num_Biquad_Sects = r;
90 
91  A_Biquad_Coef = new double[r + 1];
92  B_Biquad_Coef = new double[r + 1];
93  C_Biquad_Coef = new double[r + 1];
94 
95  //-------------------------------------------------------
96  // Eq. (5.28)
97 
98  numer = pow(10.0, passband_ripple / 20.0) + 1.0;
99 
100  vv = log(numer / (pow(10.0, passband_ripple / 20.0) - 1.)) / (2. * order);
101 
102  //-------------------------------------------------------
103  // Eq. (5.29)
104 
105  sum = 0.0;
106  for (m = 0; m < upper_summation_limit; m++) {
107  term = ipow(-1.0, m);
108  term *= ipow(modular_const, m * (m + 1));
109  term *= sinh((2. * m + 1) * vv);
110  sum = sum + term;
111  }
112  // numer = 2.0 * sum * sqrt(sqrt(modular_const));
113  numer = sum * sqrt(sqrt(modular_const));
114 
115  sum = 0.0;
116  for (m = 1; m < upper_summation_limit; m++) {
117  term = ipow(-1.0, m);
118  term *= ipow(modular_const, m * m);
119  term *= cosh(2.0 * m * vv);
120  sum = sum + term;
121  }
122  // p_sub_zero = fabs(numer/(1.+2.*sum));
123  p_sub_zero = fabs(numer / (0.5 + sum));
124 
125  //------------------------------------------
126  // Eq. (5.30)
127 
128  ww = 1.0 + selec_factor * p_sub_zero * p_sub_zero;
129  ww = sqrt(ww * (1.0 + p_sub_zero * p_sub_zero / selec_factor));
130 
131  //---------------------------------------
132  // loop for steps 8, 9, 10, of Alg 5.4
133 
134  H_Sub_Zero = 1.0;
135 
136  for (i = 1; i <= r; i++) {
137  if (order % 2) // if order is odd
138  {
139  mu = i;
140  } else // order is even
141  {
142  mu = i - 0.5;
143  }
144  //------------------------------
145  // Eq. (5.31) numerator
146 
147  sum = 0.0;
148  for (m = 0; m < upper_summation_limit; m++) {
149  term = ipow(-1.0, m);
150  term *= ipow(modular_const, m * (m + 1));
151  term *= sin((2 * m + 1) * PI * mu / order);
152  sum += term;
153  }
154  numer = 2.0 * sum * sqrt(sqrt(modular_const));
155 
156  //---------------------------------------
157  // Eq. (5.31) denominator and finish
158 
159  sum = 0.0;
160  for (m = 1; m < upper_summation_limit; m++) {
161  term = ipow(-1.0, m);
162  term *= ipow(modular_const, m * m);
163  term *= cos(2.0 * PI * m * mu / order);
164  sum += term;
165  }
166  xx = numer / (1. + 2. * sum);
167 
168  //----------------------------------------
169  // Eq. (5.32)
170 
171  yy = 1.0 - selec_factor * xx * xx;
172  yy = sqrt(yy * (1.0 - (xx * xx / selec_factor)));
173 
174  //-----------------------------------------
175  // Eq. (5.33)
176 
177  aa = 1.0 / (xx * xx);
178  aa *= (passband_edge * passband_edge);
179  A_Biquad_Coef[i] = aa;
180 
181  //-----------------------------------------
182  // Eq. (5.34)
183 
184  denom = 1.0 + ipow(p_sub_zero * xx, 2);
185  bb = 2.0 * p_sub_zero * yy / denom;
186  bb *= passband_edge;
187  B_Biquad_Coef[i] = bb;
188 
189  //-----------------------------------------
190  // Eq. (5.35)
191 
192  denom = ipow(denom, 2);
193  numer = ipow(p_sub_zero * yy, 2) + ipow(xx * ww, 2);
194  cc = numer / denom;
195  cc *= (passband_edge * passband_edge);
196  C_Biquad_Coef[i] = cc;
197 
198  //--------------------------------------------
199  //
200  H_Sub_Zero *= (cc / aa);
201 
202  //-------------------------------------------
203  // compute pair of pole locations
204  // by finding roots of s**2 + bb*s + cc = 0
205 
206  //------------------------------------------------------
207  // we need to compute:
208  // cmplx_work = sqrt((complex)(bb*bb - 4.*cc));
209  //
210  // work around for missing sqrt(complex) function
211  //
212  // (bb*bb - 4.0*cc) will always be real and negative
213  // so sqrt(bb*bb -4.0*cc) will always be pure imaginary
214  // equal to sqrt(-1)*sqrt(4.0*cc - bb*bb)
215  // therefore:
216 
217  double_temp = sqrt(4.0 * cc - bb * bb);
218  cmplx_work = complex(0.0, double_temp);
219 
220  Prototype_Pole_Locs[i] = (complex(-bb, 0.0) - cmplx_work) / 2.0;
221 
222 #ifdef _DEBUG
223  DebugFile << "in ellip response, pole[" << i
224  << "] = " << Prototype_Pole_Locs[i] << std::endl;
225 #endif
226 
227  Prototype_Pole_Locs[order + 1 - i] = (complex(-bb, 0.0) + cmplx_work) / 2.0;
228  //-----------------------------------------------------------
229  // compute pair of zero locations
230  // by finding roots of s**2 + a = 0
231  //
232  // roots = +/- sqrt(-a)
233  //
234 
235  if (order % 2) {
236  i_mirror = order - i;
237  } else {
238  i_mirror = order + 1 - i;
239  }
240  if (aa < 0.0) {
241  double_temp = sqrt(-aa);
242  Prototype_Zero_Locs[i] = complex(double_temp, 0.0);
243  Prototype_Zero_Locs[i_mirror] = complex((-double_temp), 0.0);
244  } else {
245  double_temp = sqrt(aa);
246  Prototype_Zero_Locs[i] = complex(0.0, double_temp);
247  Prototype_Zero_Locs[i_mirror] = complex(0.0, (-double_temp));
248  }
249 #ifdef _DEBUG
250  DebugFile << "in ellip response, zero[" << i
251  << "] = " << Prototype_Zero_Locs[i] << std::endl;
252 #endif
253  }
254  //---------------------------
255  // Finish up Ho
256 
257  if (order % 2) {
258  // p_sub_zero *= stopband_edge;
259  p_sub_zero *= passband_edge;
260  H_Sub_Zero *= p_sub_zero;
261  Prototype_Pole_Locs[(order + 1) / 2] = complex(-p_sub_zero, 0.0);
262 #ifdef _DEBUG
263  DebugFile << "p_sub_zero = " << p_sub_zero << std::endl;
264  DebugFile << "in ellip, H_Sub_Zero = " << H_Sub_Zero << std::endl;
265 #endif
266  } else {
267  H_Sub_Zero *= pow(10.0, passband_ripple / (-20.0));
268  }
269  return;
270 };
271 
272 //-----------------------------------------------------
273 // raise double to an integer power
274 
275 double
276 ipow(double x, int m)
277 {
278  int i;
279  double result;
280  result = 1.0;
281  for (i = 1; i <= m; i++) {
282  result *= x;
283  }
284  return (result);
285 }
int Num_Biquad_Sects
Definition: filtfunc.h:59
double * C_Biquad_Coef
Definition: filtfunc.h:65
double * B_Biquad_Coef
Definition: filtfunc.h:64
EllipticalTransFunc(int order, double passband_ripple, double stopband_ripple, double passband_edge, double stopband_edge, int upper_summation_limit)
Definition: elipfunc.cpp:20
double * A_Biquad_Coef
Definition: filtfunc.h:63
#define PI
Definition: misdefs.h:9
complex sqrt(const complex _z)
Definition: complex.h:148
double ipow(double x, int m)
Definition: elipfunc.cpp:276
double H_Sub_Zero
Definition: filtfunc.h:66
int Num_Prototype_Zeros
Definition: filtfunc.h:58
std::ofstream DebugFile
complex * Prototype_Pole_Locs
Definition: filtfunc.h:67
int Num_Prototype_Poles
Definition: filtfunc.h:57
complex * Prototype_Zero_Locs
Definition: filtfunc.h:68