1 //
2 // mkfilter-derived code
3 // ---------------------
4 //
5 // Copyright (c) 2002-2004 Jim Peters <http://uazu.net/>. This
6 // file is released under the GNU Lesser General Public License
7 // (LGPL) version 2.1 as published by the Free Software
8 // Foundation. See the file COPYING_LIB for details, or visit
9 // <http://www.fsf.org/licenses/licenses.html>.
10 //
11 // This is all code derived from 'mkfilter' by Tony Fisher of the
12 // University of York. I've rewritten it all in C, and given it
13 // a thorough overhaul, so there is actually none of his code
14 // here any more, but it is all still very strongly based on the
15 // algorithms and techniques that he used in 'mkfilter'.
16 //
17 // For those who didn't hear, Tony Fisher died in February 2000
18 // at the age of 43. See his web-site for information and a
19 // tribute:
20 //
21 // http://www-users.cs.york.ac.uk/~fisher/
22 // http://www-users.cs.york.ac.uk/~fisher/tribute.html
23 //
24 // The original C++ sources and the rest of the mkfilter tool-set
25 // are still available from his site:
26 //
27 // http://www-users.cs.york.ac.uk/~fisher/mkfilter/
28 //
29 //
30 // I've made a number of improvements and changes whilst
31 // rewriting the code in C. For example, I halved the
32 // calculations required in designing the filters by storing only
33 // one complex pole/zero out of each conjugate pair. This also
34 // made it much easier to output the filter as a list of
35 // sub-filters without lots of searching around to match up
36 // conjugate pairs later on. Outputting as a list of subfilters
37 // permits greater accuracy in calculation of the response, and
38 // also in the execution of the final filter. Also, some FIR
39 // coefficients can be marked as 'constant', allowing optimised
40 // routines to be generated for whole classes of filters, with
41 // just the variable coefficients filled in at run-time.
42 //
43 // On the down-side, complex numbers are not portably available
44 // in C before C99, so complex calculations here are done on
45 // double[] arrays with inline functions, which ends up looking
46 // more like assembly language than C. Never mind.
47 //
48
49 //
50 // LEGAL STUFF
51 // -----------
52 //
53 // Tony Fisher released his software on his University of York
54 // pages for free use and free download. The software itself has
55 // no licence terms attached, nor copyright messages, just the
56 // author's name, E-mail address and date. Nor are there any
57 // licence terms indicated on the website. I understand that
58 // under the Berne convention copyright does not have to be
59 // claimed explicitly, so these are in fact copyright files by
60 // legal default. However, the intention was obviously that
61 // these files should be used by others.
62 //
63 // None of this really helps, though, if we're going to try to be
64 // 100% legally correct, so I wrote to Anthony Moulds who is the
65 // contact name on Tony Fisher's pages now. I explained what I
66 // planned to do with the code, and he answered as follows:
67 //
68 // (Note that I was planning to use it 'as-is' at that time,
69 // rather than rewrite it as I have done now)
70 //
71 // > To: "Jim Peters" <jim@uazu.net>
72 // > From: "Anthony Moulds" <anthony@cs.york.ac.uk>
73 // > Subject: RE: mkfilter source
74 // > Date: Tue, 29 Oct 2002 15:30:19 -0000
75 // >
76 // > Hi Jim,
77 // >
78 // > Thanks for your email.
79 // >
80 // > The University will be happy to let you use Dr Fisher's mkfilter
81 // > code since your intention is not to profit financially from his work.
82 // >
83 // > It would be nice if in some way you could acknowledge his contribution.
84 // >
85 // > Best wishes and good luck with your work,
86 // >
87 // > Anthony Moulds
88 // > Senior Experimental Officer,
89 // > Computer Science Department, University of York,
90 // > York, England, UK. Tel: 44(0)1904 434758 Fax: 44(0)19042767
91 // > ============================================================
92 // >
93 // >
94 // > > -----Original Message-----
95 // > > From: Jim Peters [mailto:jim@uazu.net]
96 // > > Sent: Monday, October 28, 2002 12:36 PM
97 // > > To: anthony@cs.york.ac.uk
98 // > > Subject: mkfilter source
99 // > >
100 // > >
101 // > > I'm very sorry to hear (rather late, I know) that Tony Fisher died --
102 // > > I've always gone straight to the filter page, rather than through his
103 // > > home page. I hope his work remains available for the future.
104 // > >
105 // > > Anyway, the reason I'm writing is to clarify the status of the
106 // > > mkfilter source code. Because copyright is not claimed on the web
107 // > > page nor in the source distribution, I guess that Tony's intention was
108 // > > that this code should be in the public domain. However, I would like
109 // > > to check this now to avoid complications later.
110 // > >
111 // > > I am using his code, modified, to provide a library of filter-design
112 // > > routines for a GPL'd filter design app, which is not yet released.
113 // > > The library could also be used standalone, permitting apps to design
114 // > > filters at run-time rather than using hard-coded compile-time filters.
115 // > > My interest in filters is as a part of my work on the OpenEEG project
116 //
117 // So this looks pretty clear to me. I am not planning to profit
118 // from the work, so everything is fine with the University. I
119 // guess others might profit from the work, indirectly, as with
120 // any free software release, but so long as I don't, we're fine.
121 //
122 // I hope this is watertight enough for Debian/etc. Otherwise
123 // I'll have to go back to Anthony Moulds for clarification.
124 //
125 // Even though there is no code cut-and-pasted from 'mkfilter'
126 // here, it is all very obviously based on that code, so it
127 // probably counts as a derived work -- although as ever "I Am
128 // Not A Lawyer".
129 //
130
131
132 #ifdef HUGE_VAL
133 #define INF HUGE_VAL
134 #else
135 #define INF (1.0/0.0)
136 #endif
137
138 #define TWOPI (2*M_PI)
139
140
141 //
142 // Complex square root: aa= aa^0.5
143 //
144
145 STATIC_INLINE double
my_sqrt(double aa)146 my_sqrt(double aa) {
147 return aa <= 0.0 ? 0.0 : sqrt(aa);
148 }
149
150 STATIC_INLINE void
csqrt(double * aa)151 csqrt(double *aa) {
152 double mag= hypot(aa[0], aa[1]);
153 double rr= my_sqrt((mag + aa[0]) * 0.5);
154 double ii= my_sqrt((mag - aa[0]) * 0.5);
155 if (aa[1] < 0.0) ii= -ii;
156 aa[0]= rr;
157 aa[1]= ii;
158 }
159
160 //
161 // Complex imaginary exponent: aa= e^i.theta
162 //
163
164 STATIC_INLINE void
cexpj(double * aa,double theta)165 cexpj(double *aa, double theta) {
166 aa[0]= cos(theta);
167 aa[1]= sin(theta);
168 }
169
170 //
171 // Complex exponent: aa= e^aa
172 //
173
174 STATIC_INLINE void
cexp(double * aa)175 cexp(double *aa) {
176 double mag= exp(aa[0]);
177 aa[0]= mag * cos(aa[1]);
178 aa[1]= mag * sin(aa[1]);
179 }
180
181 //
182 // Global temp buffer for generating filters. *NOT THREAD SAFE*
183 //
184 // Note that the poles and zeros are stored in a strange way.
185 // Rather than storing both a pole (or zero) and its complex
186 // conjugate, I'm storing just one of the pair. Also, for real
187 // poles, I'm not storing the imaginary part (which is zero).
188 // This results in a list of numbers exactly half the length you
189 // might otherwise expect. However, since some of these numbers
190 // are in pairs, and some are single, we need a separate flag
191 // array to indicate which is which. poltyp[] serves this
192 // purpose. An entry is 1 if the corresponding offset is a real
193 // pole, or 2 if it is the first of a pair of values making up a
194 // complex pole. The second value of the pair has an entry of 0
195 // attached. (Similarly for zeros in zertyp[])
196 //
197
198 #define MAXPZ 64
199
200 static int n_pol; // Number of poles
201 static double pol[MAXPZ]; // Pole values (see above)
202 static char poltyp[MAXPZ]; // Pole value types: 1 real, 2 first of complex pair, 0 second
203 static int n_zer; // Same for zeros ...
204 static double zer[MAXPZ];
205 static char zertyp[MAXPZ];
206
207
208 //
209 // Pre-warp a frequency
210 //
211
212 STATIC_INLINE double
prewarp(double val)213 prewarp(double val) {
214 return tan(val * M_PI) / M_PI;
215 }
216
217
218 //
219 // Bessel poles; final one is a real value for odd numbers of
220 // poles
221 //
222
223 static double bessel_1[]= {
224 -1.00000000000e+00
225 };
226
227 static double bessel_2[]= {
228 -1.10160133059e+00, 6.36009824757e-01,
229 };
230
231 static double bessel_3[]= {
232 -1.04740916101e+00, 9.99264436281e-01,
233 -1.32267579991e+00,
234 };
235
236 static double bessel_4[]= {
237 -9.95208764350e-01, 1.25710573945e+00,
238 -1.37006783055e+00, 4.10249717494e-01,
239 };
240
241 static double bessel_5[]= {
242 -9.57676548563e-01, 1.47112432073e+00,
243 -1.38087732586e+00, 7.17909587627e-01,
244 -1.50231627145e+00,
245 };
246
247 static double bessel_6[]= {
248 -9.30656522947e-01, 1.66186326894e+00,
249 -1.38185809760e+00, 9.71471890712e-01,
250 -1.57149040362e+00, 3.20896374221e-01,
251 };
252
253 static double bessel_7[]= {
254 -9.09867780623e-01, 1.83645135304e+00,
255 -1.37890321680e+00, 1.19156677780e+00,
256 -1.61203876622e+00, 5.89244506931e-01,
257 -1.68436817927e+00,
258 };
259
260 static double bessel_8[]= {
261 -8.92869718847e-01, 1.99832584364e+00,
262 -1.37384121764e+00, 1.38835657588e+00,
263 -1.63693941813e+00, 8.22795625139e-01,
264 -1.75740840040e+00, 2.72867575103e-01,
265 };
266
267 static double bessel_9[]= {
268 -8.78399276161e-01, 2.14980052431e+00,
269 -1.36758830979e+00, 1.56773371224e+00,
270 -1.65239648458e+00, 1.03138956698e+00,
271 -1.80717053496e+00, 5.12383730575e-01,
272 -1.85660050123e+00,
273 };
274
275 static double bessel_10[]= {
276 -8.65756901707e-01, 2.29260483098e+00,
277 -1.36069227838e+00, 1.73350574267e+00,
278 -1.66181024140e+00, 1.22110021857e+00,
279 -1.84219624443e+00, 7.27257597722e-01,
280 -1.92761969145e+00, 2.41623471082e-01,
281 };
282
283 static double *bessel_poles[10]= {
284 bessel_1, bessel_2, bessel_3, bessel_4, bessel_5,
285 bessel_6, bessel_7, bessel_8, bessel_9, bessel_10
286 };
287
288 //
289 // Generate Bessel poles for the given order.
290 //
291
292 static void
bessel(int order)293 bessel(int order) {
294 int a;
295
296 if (order > 10) error("Maximum Bessel order is 10");
297 n_pol= order;
298 memcpy(pol, bessel_poles[order-1], n_pol * sizeof(double));
299
300 for (a= 0; a<order-1; ) {
301 poltyp[a++]= 2;
302 poltyp[a++]= 0;
303 }
304 if (a < order)
305 poltyp[a++]= 1;
306 }
307
308 //
309 // Generate Butterworth poles for the given order. These are
310 // regularly-spaced points on the unit circle to the left of the
311 // real==0 line.
312 //
313
314 static void
butterworth(int order)315 butterworth(int order) {
316 int a;
317 if (order > MAXPZ)
318 error("Maximum butterworth/chebyshev order is %d", MAXPZ);
319 n_pol= order;
320 for (a= 0; a<order-1; a += 2) {
321 poltyp[a]= 2;
322 poltyp[a+1]= 0;
323 cexpj(pol+a, M_PI - (order-a-1) * 0.5 * M_PI / order);
324 }
325 if (a < order) {
326 poltyp[a]= 1;
327 pol[a]= -1.0;
328 }
329 }
330
331 //
332 // Generate Chebyshev poles for the given order and ripple.
333 //
334
335 static void
chebyshev(int order,double ripple)336 chebyshev(int order, double ripple) {
337 double eps, y;
338 double sh, ch;
339 int a;
340
341 butterworth(order);
342 if (ripple >= 0.0) error("Chebyshev ripple in dB should be -ve");
343
344 eps= sqrt(-1.0 + pow(10.0, -0.1 * ripple));
345 y= asinh(1.0 / eps) / order;
346 if (y <= 0.0) error("Internal error; chebyshev y-value <= 0.0: %g", y);
347 sh= sinh(y);
348 ch= cosh(y);
349
350 for (a= 0; a<n_pol; ) {
351 if (poltyp[a] == 1)
352 pol[a++] *= sh;
353 else {
354 pol[a++] *= sh;
355 pol[a++] *= ch;
356 }
357 }
358 }
359
360
361 //
362 // Adjust raw poles to LP filter
363 //
364
365 static void
lowpass(double freq)366 lowpass(double freq) {
367 int a;
368
369 // Adjust poles
370 freq *= TWOPI;
371 for (a= 0; a<n_pol; a++)
372 pol[a] *= freq;
373
374 // Add zeros
375 n_zer= n_pol;
376 for (a= 0; a<n_zer; a++) {
377 zer[a]= -INF;
378 zertyp[a]= 1;
379 }
380 }
381
382 //
383 // Adjust raw poles to HP filter
384 //
385
386 static void
highpass(double freq)387 highpass(double freq) {
388 int a;
389
390 // Adjust poles
391 freq *= TWOPI;
392 for (a= 0; a<n_pol; ) {
393 if (poltyp[a] == 1) {
394 pol[a]= freq / pol[a];
395 a++;
396 } else {
397 crecip(pol + a);
398 pol[a++] *= freq;
399 pol[a++] *= freq;
400 }
401 }
402
403 // Add zeros
404 n_zer= n_pol;
405 for (a= 0; a<n_zer; a++) {
406 zer[a]= 0.0;
407 zertyp[a]= 1;
408 }
409 }
410
411 //
412 // Adjust raw poles to BP filter. The number of poles is
413 // doubled.
414 //
415
416 static void
bandpass(double freq1,double freq2)417 bandpass(double freq1, double freq2) {
418 double w0= TWOPI * sqrt(freq1*freq2);
419 double bw= 0.5 * TWOPI * (freq2-freq1);
420 int a, b;
421
422 if (n_pol * 2 > MAXPZ)
423 error("Maximum order for bandpass filters is %d", MAXPZ/2);
424
425 // Run through the list backwards, expanding as we go
426 for (a= n_pol, b= n_pol*2; a>0; ) {
427 // hba= pole * bw;
428 // temp= csqrt(1.0 - square(w0 / hba));
429 // pole1= hba * (1.0 + temp);
430 // pole2= hba * (1.0 - temp);
431
432 if (poltyp[a-1] == 1) {
433 double hba;
434 a--; b -= 2;
435 poltyp[b]= 2; poltyp[b+1]= 0;
436 hba= pol[a] * bw;
437 cassz(pol+b, 1.0 - (w0 / hba) * (w0 / hba), 0.0);
438 csqrt(pol+b);
439 caddz(pol+b, 1.0, 0.0);
440 cmulr(pol+b, hba);
441 } else { // Assume poltyp[] data is valid
442 double hba[2];
443 a -= 2; b -= 4;
444 poltyp[b]= 2; poltyp[b+1]= 0;
445 poltyp[b+2]= 2; poltyp[b+3]= 0;
446 cass(hba, pol+a);
447 cmulr(hba, bw);
448 cass(pol+b, hba);
449 crecip(pol+b);
450 cmulr(pol+b, w0);
451 csqu(pol+b);
452 cneg(pol+b);
453 caddz(pol+b, 1.0, 0.0);
454 csqrt(pol+b);
455 cmul(pol+b, hba);
456 cass(pol+b+2, pol+b);
457 cneg(pol+b+2);
458 cadd(pol+b, hba);
459 cadd(pol+b+2, hba);
460 }
461 }
462 n_pol *= 2;
463
464 // Add zeros
465 n_zer= n_pol;
466 for (a= 0; a<n_zer; a++) {
467 zertyp[a]= 1;
468 zer[a]= (a<n_zer/2) ? 0.0 : -INF;
469 }
470 }
471
472 //
473 // Adjust raw poles to BS filter. The number of poles is
474 // doubled.
475 //
476
477 static void
bandstop(double freq1,double freq2)478 bandstop(double freq1, double freq2) {
479 double w0= TWOPI * sqrt(freq1*freq2);
480 double bw= 0.5 * TWOPI * (freq2-freq1);
481 int a, b;
482
483 if (n_pol * 2 > MAXPZ)
484 error("Maximum order for bandstop filters is %d", MAXPZ/2);
485
486 // Run through the list backwards, expanding as we go
487 for (a= n_pol, b= n_pol*2; a>0; ) {
488 // hba= bw / pole;
489 // temp= csqrt(1.0 - square(w0 / hba));
490 // pole1= hba * (1.0 + temp);
491 // pole2= hba * (1.0 - temp);
492
493 if (poltyp[a-1] == 1) {
494 double hba;
495 a--; b -= 2;
496 poltyp[b]= 2; poltyp[b+1]= 0;
497 hba= bw / pol[a];
498 cassz(pol+b, 1.0 - (w0 / hba) * (w0 / hba), 0.0);
499 csqrt(pol+b);
500 caddz(pol+b, 1.0, 0.0);
501 cmulr(pol+b, hba);
502 } else { // Assume poltyp[] data is valid
503 double hba[2];
504 a -= 2; b -= 4;
505 poltyp[b]= 2; poltyp[b+1]= 0;
506 poltyp[b+2]= 2; poltyp[b+3]= 0;
507 cass(hba, pol+a);
508 crecip(hba);
509 cmulr(hba, bw);
510 cass(pol+b, hba);
511 crecip(pol+b);
512 cmulr(pol+b, w0);
513 csqu(pol+b);
514 cneg(pol+b);
515 caddz(pol+b, 1.0, 0.0);
516 csqrt(pol+b);
517 cmul(pol+b, hba);
518 cass(pol+b+2, pol+b);
519 cneg(pol+b+2);
520 cadd(pol+b, hba);
521 cadd(pol+b+2, hba);
522 }
523 }
524 n_pol *= 2;
525
526 // Add zeros
527 n_zer= n_pol;
528 for (a= 0; a<n_zer; a+=2) {
529 zertyp[a]= 2; zertyp[a+1]= 0;
530 zer[a]= 0.0; zer[a+1]= w0;
531 }
532 }
533
534 //
535 // Convert list of poles+zeros from S to Z using bilinear
536 // transform
537 //
538
539 static void
s2z_bilinear()540 s2z_bilinear() {
541 int a;
542 for (a= 0; a<n_pol; ) {
543 // Calculate (2 + val) / (2 - val)
544 if (poltyp[a] == 1) {
545 if (pol[a] == -INF)
546 pol[a]= -1.0;
547 else
548 pol[a]= (2 + pol[a]) / (2 - pol[a]);
549 a++;
550 } else {
551 double val[2];
552 cass(val, pol+a);
553 cneg(val);
554 caddz(val, 2, 0);
555 caddz(pol+a, 2, 0);
556 cdiv(pol+a, val);
557 a += 2;
558 }
559 }
560 for (a= 0; a<n_zer; ) {
561 // Calculate (2 + val) / (2 - val)
562 if (zertyp[a] == 1) {
563 if (zer[a] == -INF)
564 zer[a]= -1.0;
565 else
566 zer[a]= (2 + zer[a]) / (2 - zer[a]);
567 a++;
568 } else {
569 double val[2];
570 cass(val, zer+a);
571 cneg(val);
572 caddz(val, 2, 0);
573 caddz(zer+a, 2, 0);
574 cdiv(zer+a, val);
575 a += 2;
576 }
577 }
578 }
579
580 //
581 // Convert S to Z using matched-Z transform
582 //
583
584 static void
s2z_matchedZ()585 s2z_matchedZ() {
586 int a;
587
588 for (a= 0; a<n_pol; ) {
589 // Calculate cexp(val)
590 if (poltyp[a] == 1) {
591 if (pol[a] == -INF)
592 pol[a]= 0.0;
593 else
594 pol[a]= exp(pol[a]);
595 a++;
596 } else {
597 cexp(pol+a);
598 a += 2;
599 }
600 }
601
602 for (a= 0; a<n_zer; ) {
603 // Calculate cexp(val)
604 if (zertyp[a] == 1) {
605 if (zer[a] == -INF)
606 zer[a]= 0.0;
607 else
608 zer[a]= exp(zer[a]);
609 a++;
610 } else {
611 cexp(zer+a);
612 a += 2;
613 }
614 }
615 }
616
617 //
618 // Generate a FidFilter for the current set of poles and zeros.
619 // The given gain is inserted at the start of the FidFilter as a
620 // one-coefficient FIR filter. This is positioned to be easily
621 // adjusted later to correct the filter gain.
622 //
623 // 'cbm' should be a bitmap indicating which FIR coefficients are
624 // constants for this filter type. Normal values are ~0 for all
625 // constant, or 0 for none constant, or some other bitmask for a
626 // mixture. Filters generated with lowpass(), highpass() and
627 // bandpass() above should pass ~0, but bandstop() requires 0x5.
628 //
629 // This routine requires that any lone real poles/zeros are at
630 // the end of the list. All other poles/zeros are handled in
631 // pairs (whether pairs of real poles/zeros, or conjugate pairs).
632 //
633
634 static FidFilter*
z2fidfilter(double gain,int cbm)635 z2fidfilter(double gain, int cbm) {
636 int n_head, n_val;
637 int a;
638 FidFilter *rv;
639 FidFilter *ff;
640
641 n_head= 1 + n_pol + n_zer; // Worst case: gain + 2-element IIR/FIR
642 n_val= 1 + 2 * (n_pol+n_zer); // for each pole/zero
643
644 rv= ff= FFALLOC(n_head, n_val);
645
646 ff->typ= 'F';
647 ff->len= 1;
648 ff->val[0]= gain;
649 ff= FFNEXT(ff);
650
651 // Output as much as possible as 2x2 IIR/FIR filters
652 for (a= 0; a <= n_pol-2 && a <= n_zer-2; a += 2) {
653 // Look for a pair of values for an IIR
654 if (poltyp[a] == 1 && poltyp[a+1] == 1) {
655 // Two real values
656 ff->typ= 'I';
657 ff->len= 3;
658 ff->val[0]= 1;
659 ff->val[1]= -(pol[a] + pol[a+1]);
660 ff->val[2]= pol[a] * pol[a+1];
661 ff= FFNEXT(ff);
662 } else if (poltyp[a] == 2) {
663 // A complex value and its conjugate pair
664 ff->typ= 'I';
665 ff->len= 3;
666 ff->val[0]= 1;
667 ff->val[1]= -2 * pol[a];
668 ff->val[2]= pol[a] * pol[a] + pol[a+1] * pol[a+1];
669 ff= FFNEXT(ff);
670 } else error("Internal error -- bad poltyp[] values for z2fidfilter()");
671
672 // Look for a pair of values for an FIR
673 if (zertyp[a] == 1 && zertyp[a+1] == 1) {
674 // Two real values
675 // Skip if constant and 0/0
676 if (!cbm || zer[a] != 0.0 || zer[a+1] != 0.0) {
677 ff->typ= 'F';
678 ff->cbm= cbm;
679 ff->len= 3;
680 ff->val[0]= 1;
681 ff->val[1]= -(zer[a] + zer[a+1]);
682 ff->val[2]= zer[a] * zer[a+1];
683 ff= FFNEXT(ff);
684 }
685 } else if (zertyp[a] == 2) {
686 // A complex value and its conjugate pair
687 // Skip if constant and 0/0
688 if (!cbm || zer[a] != 0.0 || zer[a+1] != 0.0) {
689 ff->typ= 'F';
690 ff->cbm= cbm;
691 ff->len= 3;
692 ff->val[0]= 1;
693 ff->val[1]= -2 * zer[a];
694 ff->val[2]= zer[a] * zer[a] + zer[a+1] * zer[a+1];
695 ff= FFNEXT(ff);
696 }
697 } else error("Internal error -- bad zertyp[] values");
698 }
699
700 // Clear up any remaining bits and pieces. Should only be a 1x1
701 // IIR/FIR.
702 if (n_pol-a == 0 && n_zer-a == 0)
703 ;
704 else if (n_pol-a == 1 && n_zer-a == 1) {
705 if (poltyp[a] != 1 || zertyp[a] != 1)
706 error("Internal error; bad poltyp or zertyp for final pole/zero");
707 ff->typ= 'I';
708 ff->len= 2;
709 ff->val[0]= 1;
710 ff->val[1]= -pol[a];
711 ff= FFNEXT(ff);
712
713 // Skip FIR if it is constant and zero
714 if (!cbm || zer[a] != 0.0) {
715 ff->typ= 'F';
716 ff->cbm= cbm;
717 ff->len= 2;
718 ff->val[0]= 1;
719 ff->val[1]= -zer[a];
720 ff= FFNEXT(ff);
721 }
722 } else
723 error("Internal error: unexpected poles/zeros at end of list");
724
725 // End of list
726 ff->typ= 0;
727 ff->len= 0;
728 ff= FFNEXT(ff);
729
730 rv= realloc(rv, ((char*)ff)-((char*)rv));
731 if (!rv) error("Out of memory");
732 return rv;
733 }
734
735 //
736 // Setup poles/zeros for a band-pass resonator. 'qfact' gives
737 // the Q-factor; 0 is a special value indicating +infinity,
738 // giving an oscillator.
739 //
740
741 static void
bandpass_res(double freq,double qfact)742 bandpass_res(double freq, double qfact) {
743 double mag;
744 double th0, th1, th2;
745 double theta= freq * TWOPI;
746 double val[2];
747 double tmp1[2], tmp2[2], tmp3[2], tmp4[2];
748 int cnt;
749
750 n_pol= 2;
751 poltyp[0]= 2; poltyp[1]= 0;
752 n_zer= 2;
753 zertyp[0]= 1; zertyp[1]= 1;
754 zer[0]= 1; zer[1]= -1;
755
756 if (qfact == 0.0) {
757 cexpj(pol, theta);
758 return;
759 }
760
761 // Do a full binary search, rather than seeding it as Tony Fisher does
762 cexpj(val, theta);
763 mag= exp(-theta / (2.0 * qfact));
764 th0= 0; th2= M_PI;
765 for (cnt= 60; cnt > 0; cnt--) {
766 th1= 0.5 * (th0 + th2);
767 cexpj(pol, th1);
768 cmulr(pol, mag);
769
770 // Evaluate response of filter for Z= val
771 memcpy(tmp1, val, 2*sizeof(double));
772 memcpy(tmp2, val, 2*sizeof(double));
773 memcpy(tmp3, val, 2*sizeof(double));
774 memcpy(tmp4, val, 2*sizeof(double));
775 csubz(tmp1, 1, 0);
776 csubz(tmp2, -1, 0);
777 cmul(tmp1, tmp2);
778 csub(tmp3, pol); cconj(pol);
779 csub(tmp4, pol); cconj(pol);
780 cmul(tmp3, tmp4);
781 cdiv(tmp1, tmp3);
782
783 if (fabs(tmp1[1] / tmp1[0]) < 1e-10) break;
784
785 //printf("%-24.16g%-24.16g -> %-24.16g%-24.16g\n", th0, th2, tmp1[0], tmp1[1]);
786
787 if (tmp1[1] > 0.0) th2= th1;
788 else th0= th1;
789 }
790
791 if (cnt <= 0) fprintf(stderr, "Resonator binary search failed to converge");
792 }
793
794 //
795 // Setup poles/zeros for a bandstop resonator
796 //
797
798 static void
bandstop_res(double freq,double qfact)799 bandstop_res(double freq, double qfact) {
800 bandpass_res(freq, qfact);
801 zertyp[0]= 2; zertyp[1]= 0;
802 cexpj(zer, TWOPI * freq);
803 }
804
805 //
806 // Setup poles/zeros for an allpass resonator
807 //
808
809 static void
allpass_res(double freq,double qfact)810 allpass_res(double freq, double qfact) {
811 bandpass_res(freq, qfact);
812 zertyp[0]= 2; zertyp[1]= 0;
813 memcpy(zer, pol, 2*sizeof(double));
814 cmulr(zer, 1.0 / (zer[0]*zer[0] + zer[1]*zer[1]));
815 }
816
817 //
818 // Setup poles/zeros for a proportional-integral filter
819 //
820
821 static void
prop_integral(double freq)822 prop_integral(double freq) {
823 n_pol= 1;
824 poltyp[0]= 1;
825 pol[0]= 0.0;
826 n_zer= 1;
827 zertyp[0]= 1;
828 zer[0]= -TWOPI * freq;
829 }
830
831 // END //
832