« back to all changes in this revision
Viewing changes to mixxx/lib/fidlib-0.9.10/fidmkf.h
- browse files at revision 2841
- compare with another revision
- download diff
- download tarball
- view history from revision 2841
- Committer: Alex Barker
- Date: 2011-07-24 21:06:29 UTC
- Revision ID: alex@1stleg.com-20110724210629-zyr2khq2o0nv8iru
Updated libsoundtouch and fidlib to latest versions
- files added:
- mixxx/lib/fidlib-0.9.10
- mixxx/lib/soundtouch-1.6.0
- files removed:
- mixxx/lib/fidlib-0.9.9
- mixxx/lib/soundtouch-1.5.0
- files modified:
- mixxx/build/depends.py
added
removed
mixxx/lib/fidlib-0.9.10/fidmkf.h
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
146
my_sqrt(double aa) {
147
return aa <= 0.0 ? 0.0 : sqrt(aa);
148
}
149
150
// 'csqrt' clashes with builtin in GCC 4, so call it 'c_sqrt'
151
STATIC_INLINE void
152
c_sqrt(double *aa) {
153
double mag= hypot(aa[0], aa[1]);
154
double rr= my_sqrt((mag + aa[0]) * 0.5);
155
double ii= my_sqrt((mag - aa[0]) * 0.5);
156
if (aa[1] < 0.0) ii= -ii;
157
aa[0]= rr;
158
aa[1]= ii;
159
}
160
161
//
162
// Complex imaginary exponent: aa= e^i.theta
163
//
164
165
STATIC_INLINE void
166
cexpj(double *aa, double theta) {
167
aa[0]= cos(theta);
168
aa[1]= sin(theta);
169
}
170
171
//
172
// Complex exponent: aa= e^aa
173
//
174
175
// 'cexp' clashes with builtin in GCC 4, so call it 'c_exp'
176
STATIC_INLINE void
177
c_exp(double *aa) {
178
double mag= exp(aa[0]);
179
aa[0]= mag * cos(aa[1]);
180
aa[1]= mag * sin(aa[1]);
181
}
182
183
//
184
// Global temp buffer for generating filters. *NOT THREAD SAFE*
185
//
186
// Note that the poles and zeros are stored in a strange way.
187
// Rather than storing both a pole (or zero) and its complex
188
// conjugate, I'm storing just one of the pair. Also, for real
189
// poles, I'm not storing the imaginary part (which is zero).
190
// This results in a list of numbers exactly half the length you
191
// might otherwise expect. However, since some of these numbers
192
// are in pairs, and some are single, we need a separate flag
193
// array to indicate which is which. poltyp[] serves this
194
// purpose. An entry is 1 if the corresponding offset is a real
195
// pole, or 2 if it is the first of a pair of values making up a
196
// complex pole. The second value of the pair has an entry of 0
197
// attached. (Similarly for zeros in zertyp[])
198
//
199
200
#define MAXPZ 64
201
202
static int n_pol; // Number of poles
203
static double pol[MAXPZ]; // Pole values (see above)
204
static char poltyp[MAXPZ]; // Pole value types: 1 real, 2 first of complex pair, 0 second
205
static int n_zer; // Same for zeros ...
206
static double zer[MAXPZ];
207
static char zertyp[MAXPZ];
208
209
210
//
211
// Pre-warp a frequency
212
//
213
214
STATIC_INLINE double
215
prewarp(double val) {
216
return tan(val * M_PI) / M_PI;
217
}
218
219
220
//
221
// Bessel poles; final one is a real value for odd numbers of
222
// poles
223
//
224
225
static double bessel_1[]= {
226
-1.00000000000e+00
227
};
228
229
static double bessel_2[]= {
230
-1.10160133059e+00, 6.36009824757e-01,
231
};
232
233
static double bessel_3[]= {
234
-1.04740916101e+00, 9.99264436281e-01,
235
-1.32267579991e+00,
236
};
237
238
static double bessel_4[]= {
239
-9.95208764350e-01, 1.25710573945e+00,
240
-1.37006783055e+00, 4.10249717494e-01,
241
};
242
243
static double bessel_5[]= {
244
-9.57676548563e-01, 1.47112432073e+00,
245
-1.38087732586e+00, 7.17909587627e-01,
246
-1.50231627145e+00,
247
};
248
249
static double bessel_6[]= {
250
-9.30656522947e-01, 1.66186326894e+00,
251
-1.38185809760e+00, 9.71471890712e-01,
252
-1.57149040362e+00, 3.20896374221e-01,
253
};
254
255
static double bessel_7[]= {
256
-9.09867780623e-01, 1.83645135304e+00,
257
-1.37890321680e+00, 1.19156677780e+00,
258
-1.61203876622e+00, 5.89244506931e-01,
259
-1.68436817927e+00,
260
};
261
262
static double bessel_8[]= {
263
-8.92869718847e-01, 1.99832584364e+00,
264
-1.37384121764e+00, 1.38835657588e+00,
265
-1.63693941813e+00, 8.22795625139e-01,
266
-1.75740840040e+00, 2.72867575103e-01,
267
};
268
269
static double bessel_9[]= {
270
-8.78399276161e-01, 2.14980052431e+00,
271
-1.36758830979e+00, 1.56773371224e+00,
272
-1.65239648458e+00, 1.03138956698e+00,
273
-1.80717053496e+00, 5.12383730575e-01,
274
-1.85660050123e+00,
275
};
276
277
static double bessel_10[]= {
278
-8.65756901707e-01, 2.29260483098e+00,
279
-1.36069227838e+00, 1.73350574267e+00,
280
-1.66181024140e+00, 1.22110021857e+00,
281
-1.84219624443e+00, 7.27257597722e-01,
282
-1.92761969145e+00, 2.41623471082e-01,
283
};
284
285
static double *bessel_poles[10]= {
286
bessel_1, bessel_2, bessel_3, bessel_4, bessel_5,
287
bessel_6, bessel_7, bessel_8, bessel_9, bessel_10
288
};
289
290
//
291
// Generate Bessel poles for the given order.
292
//
293
294
static void
295
bessel(int order) {
296
int a;
297
298
if (order > 10) error("Maximum Bessel order is 10");
299
n_pol= order;
300
memcpy(pol, bessel_poles[order-1], n_pol * sizeof(double));
301
302
for (a= 0; a<order-1; ) {
303
poltyp[a++]= 2;
304
poltyp[a++]= 0;
305
}
306
if (a < order)
307
poltyp[a++]= 1;
308
}
309
310
//
311
// Generate Butterworth poles for the given order. These are
312
// regularly-spaced points on the unit circle to the left of the
313
// real==0 line.
314
//
315
316
static void
317
butterworth(int order) {
318
int a;
319
if (order > MAXPZ)
320
error("Maximum butterworth/chebyshev order is %d", MAXPZ);
321
n_pol= order;
322
for (a= 0; a<order-1; a += 2) {
323
poltyp[a]= 2;
324
poltyp[a+1]= 0;
325
cexpj(pol+a, M_PI - (order-a-1) * 0.5 * M_PI / order);
326
}
327
if (a < order) {
328
poltyp[a]= 1;
329
pol[a]= -1.0;
330
}
331
}
332
333
//
334
// Generate Chebyshev poles for the given order and ripple.
335
//
336
337
static void
338
chebyshev(int order, double ripple) {
339
double eps, y;
340
double sh, ch;
341
int a;
342
343
butterworth(order);
344
if (ripple >= 0.0) error("Chebyshev ripple in dB should be -ve");
345
346
eps= sqrt(-1.0 + pow(10.0, -0.1 * ripple));
347
y= asinh(1.0 / eps) / order;
348
if (y <= 0.0) error("Internal error; chebyshev y-value <= 0.0: %g", y);
349
sh= sinh(y);
350
ch= cosh(y);
351
352
for (a= 0; a<n_pol; ) {
353
if (poltyp[a] == 1)
354
pol[a++] *= sh;
355
else {
356
pol[a++] *= sh;
357
pol[a++] *= ch;
358
}
359
}
360
}
361
362
363
//
364
// Adjust raw poles to LP filter
365
//
366
367
static void
368
lowpass(double freq) {
369
int a;
370
371
// Adjust poles
372
freq *= TWOPI;
373
for (a= 0; a<n_pol; a++)
374
pol[a] *= freq;
375
376
// Add zeros
377
n_zer= n_pol;
378
for (a= 0; a<n_zer; a++) {
379
zer[a]= -INF;
380
zertyp[a]= 1;
381
}
382
}
383
384
//
385
// Adjust raw poles to HP filter
386
//
387
388
static void
389
highpass(double freq) {
390
int a;
391
392
// Adjust poles
393
freq *= TWOPI;
394
for (a= 0; a<n_pol; ) {
395
if (poltyp[a] == 1) {
396
pol[a]= freq / pol[a];
397
a++;
398
} else {
399
crecip(pol + a);
400
pol[a++] *= freq;
401
pol[a++] *= freq;
402
}
403
}
404
405
// Add zeros
406
n_zer= n_pol;
407
for (a= 0; a<n_zer; a++) {
408
zer[a]= 0.0;
409
zertyp[a]= 1;
410
}
411
}
412
413
//
414
// Adjust raw poles to BP filter. The number of poles is
415
// doubled.
416
//
417
418
static void
419
bandpass(double freq1, double freq2) {
420
double w0= TWOPI * sqrt(freq1*freq2);
421
double bw= 0.5 * TWOPI * (freq2-freq1);
422
int a, b;
423
424
if (n_pol * 2 > MAXPZ)
425
error("Maximum order for bandpass filters is %d", MAXPZ/2);
426
427
// Run through the list backwards, expanding as we go
428
for (a= n_pol, b= n_pol*2; a>0; ) {
429
// hba= pole * bw;
430
// temp= c_sqrt(1.0 - square(w0 / hba));
431
// pole1= hba * (1.0 + temp);
432
// pole2= hba * (1.0 - temp);
433
434
if (poltyp[a-1] == 1) {
435
double hba;
436
a--; b -= 2;
437
poltyp[b]= 2; poltyp[b+1]= 0;
438
hba= pol[a] * bw;
439
cassz(pol+b, 1.0 - (w0 / hba) * (w0 / hba), 0.0);
440
c_sqrt(pol+b);
441
caddz(pol+b, 1.0, 0.0);
442
cmulr(pol+b, hba);
443
} else { // Assume poltyp[] data is valid
444
double hba[2];
445
a -= 2; b -= 4;
446
poltyp[b]= 2; poltyp[b+1]= 0;
447
poltyp[b+2]= 2; poltyp[b+3]= 0;
448
cass(hba, pol+a);
449
cmulr(hba, bw);
450
cass(pol+b, hba);
451
crecip(pol+b);
452
cmulr(pol+b, w0);
453
csqu(pol+b);
454
cneg(pol+b);
455
caddz(pol+b, 1.0, 0.0);
456
c_sqrt(pol+b);
457
cmul(pol+b, hba);
458
cass(pol+b+2, pol+b);
459
cneg(pol+b+2);
460
cadd(pol+b, hba);
461
cadd(pol+b+2, hba);
462
}
463
}
464
n_pol *= 2;
465
466
// Add zeros
467
n_zer= n_pol;
468
for (a= 0; a<n_zer; a++) {
469
zertyp[a]= 1;
470
zer[a]= (a<n_zer/2) ? 0.0 : -INF;
471
}
472
}
473
474
//
475
// Adjust raw poles to BS filter. The number of poles is
476
// doubled.
477
//
478
479
static void
480
bandstop(double freq1, double freq2) {
481
double w0= TWOPI * sqrt(freq1*freq2);
482
double bw= 0.5 * TWOPI * (freq2-freq1);
483
int a, b;
484
485
if (n_pol * 2 > MAXPZ)
486
error("Maximum order for bandstop filters is %d", MAXPZ/2);
487
488
// Run through the list backwards, expanding as we go
489
for (a= n_pol, b= n_pol*2; a>0; ) {
490
// hba= bw / pole;
491
// temp= c_sqrt(1.0 - square(w0 / hba));
492
// pole1= hba * (1.0 + temp);
493
// pole2= hba * (1.0 - temp);
494
495
if (poltyp[a-1] == 1) {
496
double hba;
497
a--; b -= 2;
498
poltyp[b]= 2; poltyp[b+1]= 0;
499
hba= bw / pol[a];
500
cassz(pol+b, 1.0 - (w0 / hba) * (w0 / hba), 0.0);
501
c_sqrt(pol+b);
502
caddz(pol+b, 1.0, 0.0);
503
cmulr(pol+b, hba);
504
} else { // Assume poltyp[] data is valid
505
double hba[2];
506
a -= 2; b -= 4;
507
poltyp[b]= 2; poltyp[b+1]= 0;
508
poltyp[b+2]= 2; poltyp[b+3]= 0;
509
cass(hba, pol+a);
510
crecip(hba);
511
cmulr(hba, bw);
512
cass(pol+b, hba);
513
crecip(pol+b);
514
cmulr(pol+b, w0);
515
csqu(pol+b);
516
cneg(pol+b);
517
caddz(pol+b, 1.0, 0.0);
518
c_sqrt(pol+b);
519
cmul(pol+b, hba);
520
cass(pol+b+2, pol+b);
521
cneg(pol+b+2);
522
cadd(pol+b, hba);
523
cadd(pol+b+2, hba);
524
}
525
}
526
n_pol *= 2;
527
528
// Add zeros
529
n_zer= n_pol;
530
for (a= 0; a<n_zer; a+=2) {
531
zertyp[a]= 2; zertyp[a+1]= 0;
532
zer[a]= 0.0; zer[a+1]= w0;
533
}
534
}
535
536
//
537
// Convert list of poles+zeros from S to Z using bilinear
538
// transform
539
//
540
541
static void
542
s2z_bilinear() {
543
int a;
544
for (a= 0; a<n_pol; ) {
545
// Calculate (2 + val) / (2 - val)
546
if (poltyp[a] == 1) {
547
if (pol[a] == -INF)
548
pol[a]= -1.0;
549
else
550
pol[a]= (2 + pol[a]) / (2 - pol[a]);
551
a++;
552
} else {
553
double val[2];
554
cass(val, pol+a);
555
cneg(val);
556
caddz(val, 2, 0);
557
caddz(pol+a, 2, 0);
558
cdiv(pol+a, val);
559
a += 2;
560
}
561
}
562
for (a= 0; a<n_zer; ) {
563
// Calculate (2 + val) / (2 - val)
564
if (zertyp[a] == 1) {
565
if (zer[a] == -INF)
566
zer[a]= -1.0;
567
else
568
zer[a]= (2 + zer[a]) / (2 - zer[a]);
569
a++;
570
} else {
571
double val[2];
572
cass(val, zer+a);
573
cneg(val);
574
caddz(val, 2, 0);
575
caddz(zer+a, 2, 0);
576
cdiv(zer+a, val);
577
a += 2;
578
}
579
}
580
}
581
582
//
583
// Convert S to Z using matched-Z transform
584
//
585
586
static void
587
s2z_matchedZ() {
588
int a;
589
590
for (a= 0; a<n_pol; ) {
591
// Calculate cexp(val)
592
if (poltyp[a] == 1) {
593
if (pol[a] == -INF)
594
pol[a]= 0.0;
595
else
596
pol[a]= exp(pol[a]);
597
a++;
598
} else {
599
c_exp(pol+a);
600
a += 2;
601
}
602
}
603
604
for (a= 0; a<n_zer; ) {
605
// Calculate cexp(val)
606
if (zertyp[a] == 1) {
607
if (zer[a] == -INF)
608
zer[a]= 0.0;
609
else
610
zer[a]= exp(zer[a]);
611
a++;
612
} else {
613
c_exp(zer+a);
614
a += 2;
615
}
616
}
617
}
618
619
//
620
// Generate a FidFilter for the current set of poles and zeros.
621
// The given gain is inserted at the start of the FidFilter as a
622
// one-coefficient FIR filter. This is positioned to be easily
623
// adjusted later to correct the filter gain.
624
//
625
// 'cbm' should be a bitmap indicating which FIR coefficients are
626
// constants for this filter type. Normal values are ~0 for all
627
// constant, or 0 for none constant, or some other bitmask for a
628
// mixture. Filters generated with lowpass(), highpass() and
629
// bandpass() above should pass ~0, but bandstop() requires 0x5.
630
//
631
// This routine requires that any lone real poles/zeros are at
632
// the end of the list. All other poles/zeros are handled in
633
// pairs (whether pairs of real poles/zeros, or conjugate pairs).
634
//
635
636
static FidFilter*
637
z2fidfilter(double gain, int cbm) {
638
int n_head, n_val;
639
int a;
640
FidFilter *rv;
641
FidFilter *ff;
642
643
n_head= 1 + n_pol + n_zer; // Worst case: gain + 2-element IIR/FIR
644
n_val= 1 + 2 * (n_pol+n_zer); // for each pole/zero
645
646
rv= ff= FFALLOC(n_head, n_val);
647
648
ff->typ= 'F';
649
ff->len= 1;
650
ff->val[0]= gain;
651
ff= FFNEXT(ff);
652
653
// Output as much as possible as 2x2 IIR/FIR filters
654
for (a= 0; a <= n_pol-2 && a <= n_zer-2; a += 2) {
655
// Look for a pair of values for an IIR
656
if (poltyp[a] == 1 && poltyp[a+1] == 1) {
657
// Two real values
658
ff->typ= 'I';
659
ff->len= 3;
660
ff->val[0]= 1;
661
ff->val[1]= -(pol[a] + pol[a+1]);
662
ff->val[2]= pol[a] * pol[a+1];
663
ff= FFNEXT(ff);
664
} else if (poltyp[a] == 2) {
665
// A complex value and its conjugate pair
666
ff->typ= 'I';
667
ff->len= 3;
668
ff->val[0]= 1;
669
ff->val[1]= -2 * pol[a];
670
ff->val[2]= pol[a] * pol[a] + pol[a+1] * pol[a+1];
671
ff= FFNEXT(ff);
672
} else error("Internal error -- bad poltyp[] values for z2fidfilter()");
673
674
// Look for a pair of values for an FIR
675
if (zertyp[a] == 1 && zertyp[a+1] == 1) {
676
// Two real values
677
// Skip if constant and 0/0
678
if (!cbm || zer[a] != 0.0 || zer[a+1] != 0.0) {
679
ff->typ= 'F';
680
ff->cbm= cbm;
681
ff->len= 3;
682
ff->val[0]= 1;
683
ff->val[1]= -(zer[a] + zer[a+1]);
684
ff->val[2]= zer[a] * zer[a+1];
685
ff= FFNEXT(ff);
686
}
687
} else if (zertyp[a] == 2) {
688
// A complex value and its conjugate pair
689
// Skip if constant and 0/0
690
if (!cbm || zer[a] != 0.0 || zer[a+1] != 0.0) {
691
ff->typ= 'F';
692
ff->cbm= cbm;
693
ff->len= 3;
694
ff->val[0]= 1;
695
ff->val[1]= -2 * zer[a];
696
ff->val[2]= zer[a] * zer[a] + zer[a+1] * zer[a+1];
697
ff= FFNEXT(ff);
698
}
699
} else error("Internal error -- bad zertyp[] values");
700
}
701
702
// Clear up any remaining bits and pieces. Should only be a 1x1
703
// IIR/FIR.
704
if (n_pol-a == 0 && n_zer-a == 0)
705
;
706
else if (n_pol-a == 1 && n_zer-a == 1) {
707
if (poltyp[a] != 1 || zertyp[a] != 1)
708
error("Internal error; bad poltyp or zertyp for final pole/zero");
709
ff->typ= 'I';
710
ff->len= 2;
711
ff->val[0]= 1;
712
ff->val[1]= -pol[a];
713
ff= FFNEXT(ff);
714
715
// Skip FIR if it is constant and zero
716
if (!cbm || zer[a] != 0.0) {
717
ff->typ= 'F';
718
ff->cbm= cbm;
719
ff->len= 2;
720
ff->val[0]= 1;
721
ff->val[1]= -zer[a];
722
ff= FFNEXT(ff);
723
}
724
} else
725
error("Internal error: unexpected poles/zeros at end of list");
726
727
// End of list
728
ff->typ= 0;
729
ff->len= 0;
730
ff= FFNEXT(ff);
731
732
rv= realloc(rv, ((char*)ff)-((char*)rv));
733
if (!rv) error("Out of memory");
734
return rv;
735
}
736
737
//
738
// Setup poles/zeros for a band-pass resonator. 'qfact' gives
739
// the Q-factor; 0 is a special value indicating +infinity,
740
// giving an oscillator.
741
//
742
743
static void
744
bandpass_res(double freq, double qfact) {
745
double mag;
746
double th0, th1, th2;
747
double theta= freq * TWOPI;
748
double val[2];
749
double tmp1[2], tmp2[2], tmp3[2], tmp4[2];
750
int cnt;
751
752
n_pol= 2;
753
poltyp[0]= 2; poltyp[1]= 0;
754
n_zer= 2;
755
zertyp[0]= 1; zertyp[1]= 1;
756
zer[0]= 1; zer[1]= -1;
757
758
if (qfact == 0.0) {
759
cexpj(pol, theta);
760
return;
761
}
762
763
// Do a full binary search, rather than seeding it as Tony Fisher does
764
cexpj(val, theta);
765
mag= exp(-theta / (2.0 * qfact));
766
th0= 0; th2= M_PI;
767
for (cnt= 60; cnt > 0; cnt--) {
768
th1= 0.5 * (th0 + th2);
769
cexpj(pol, th1);
770
cmulr(pol, mag);
771
772
// Evaluate response of filter for Z= val
773
memcpy(tmp1, val, 2*sizeof(double));
774
memcpy(tmp2, val, 2*sizeof(double));
775
memcpy(tmp3, val, 2*sizeof(double));
776
memcpy(tmp4, val, 2*sizeof(double));
777
csubz(tmp1, 1, 0);
778
csubz(tmp2, -1, 0);
779
cmul(tmp1, tmp2);
780
csub(tmp3, pol); cconj(pol);
781
csub(tmp4, pol); cconj(pol);
782
cmul(tmp3, tmp4);
783
cdiv(tmp1, tmp3);
784
785
if (fabs(tmp1[1] / tmp1[0]) < 1e-10) break;
786
787
//printf("%-24.16g%-24.16g -> %-24.16g%-24.16g\n", th0, th2, tmp1[0], tmp1[1]);
788
789
if (tmp1[1] > 0.0) th2= th1;
790
else th0= th1;
791
}
792
793
if (cnt <= 0) fprintf(stderr, "Resonator binary search failed to converge");
794
}
795
796
//
797
// Setup poles/zeros for a bandstop resonator
798
//
799
800
static void
801
bandstop_res(double freq, double qfact) {
802
bandpass_res(freq, qfact);
803
zertyp[0]= 2; zertyp[1]= 0;
804
cexpj(zer, TWOPI * freq);
805
}
806
807
//
808
// Setup poles/zeros for an allpass resonator
809
//
810
811
static void
812
allpass_res(double freq, double qfact) {
813
bandpass_res(freq, qfact);
814
zertyp[0]= 2; zertyp[1]= 0;
815
memcpy(zer, pol, 2*sizeof(double));
816
cmulr(zer, 1.0 / (zer[0]*zer[0] + zer[1]*zer[1]));
817
}
818
819
//
820
// Setup poles/zeros for a proportional-integral filter
821
//
822
823
static void
824
prop_integral(double freq) {
825
n_pol= 1;
826
poltyp[0]= 1;
827
pol[0]= 0.0;
828
n_zer= 1;
829
zertyp[0]= 1;
830
zer[0]= -TWOPI * freq;
831
}
832
833
// END //
Loggerhead is a web-based interface for Breezy
Version: 2.0.1