1 #!/usr/bin/awk -f
2 #
3 # SPDX-License-Identifier: BSD-2-Clause-FreeBSD
4 #
5 # Copyright (c) 2007-2009 Ariff Abdullah <ariff@FreeBSD.org>
6 # All rights reserved.
7 #
8 # Redistribution and use in source and binary forms, with or without
9 # modification, are permitted provided that the following conditions
10 # are met:
11 # 1. Redistributions of source code must retain the above copyright
12 # notice, this list of conditions and the following disclaimer.
13 # 2. Redistributions in binary form must reproduce the above copyright
14 # notice, this list of conditions and the following disclaimer in the
15 # documentation and/or other materials provided with the distribution.
16 #
17 # THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
18 # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 # ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
21 # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22 # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23 # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26 # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27 # SUCH DAMAGE.
28 #
29 # $FreeBSD$
30 #
31
32 #
33 # FIR filter design by windowing method. This might become one of the
34 # funniest joke I've ever written due to too many tricks being applied to
35 # ensure maximum precision (well, in fact this is already have the same
36 # precision granularity compared to its C counterpart). Nevertheless, it's
37 # working, precise, dynamically tunable based on "presets".
38 #
39 # XXX EXPECT TOTAL REWRITE! DON'T ARGUE!
40 #
41 # TODO: Using ultraspherical window might be a good idea.
42 #
43 # Based on:
44 #
45 # "Digital Audio Resampling" by Julius O. Smith III
46 #
47 # - http://ccrma.stanford.edu/~jos/resample/
48 #
49
50 #
51 # Some basic Math functions.
52 #
53 function abs(x)
54 {
55 return (((x < 0) ? -x : x) + 0);
56 }
57
58 function fabs(x)
59 {
60 return (((x < 0.0) ? -x : x) + 0.0);
61 }
62
63 function ceil(x, r)
64 {
65 r = int(x);
66 if (r < x)
67 r++;
68 return (r + 0);
69 }
70
71 function floor(x, r)
72 {
73 r = int(x);
74 if (r > x)
75 r--;
76 return (r + 0);
77 }
78
79 function pow(x, y)
80 {
81 return (exp(1.0 * y * log(1.0 * x)));
82 }
83
84 #
85 # What the hell...
86 #
87 function shl(x, y)
88 {
89 while (y > 0) {
90 x *= 2;
91 y--;
92 }
93 return (x);
94 }
95
96 function shr(x, y)
97 {
98 while (y > 0 && x != 0) {
99 x = floor(x / 2);
100 y--;
101 }
102 return (x);
103 }
104
105 function fx_floor(v, o, r)
106 {
107 if (fabs(v) < fabs(smallest))
108 smallest = v;
109 if (fabs(v) > fabs(largest))
110 largest = v;
111
112 r = floor((v * o) + 0.5);
113 if (r < INT32_MIN || r > INT32_MAX)
114 printf("\n#error overflow v=%f, please reduce %d\n", v, o);
115
116 return (r);
117 }
118
119 #
120 # Kaiser linear piecewise functions.
121 #
122 function kaiserAttn2Beta(attn, beta)
123 {
124 if (attn < 0.0)
125 return (Z_KAISER_BETA_DEFAULT);
126
127 if (attn > 50.0)
128 beta = 0.1102 * ((1.0 * attn) - 8.7);
129 else if (attn > 21.0)
130 beta = (0.5842 * pow((1.0 * attn) - 21.0, 0.4)) + \
131 (0.07886 * ((1.0 * attn) - 21.0));
132 else
133 beta = 0.0;
134
135 return (beta);
136 }
137
138 function kaiserBeta2Attn(beta, x, y, i, attn, xbeta)
139 {
140 if (beta < Z_WINDOW_KAISER)
141 return (Z_KAISER_ATTN_DEFAULT);
142
143 if (beta > kaiserAttn2Beta(50.0))
144 attn = ((1.0 * beta) / 0.1102) + 8.7;
145 else {
146 x = 21.0;
147 y = 50.0;
148 attn = 0.5 * (x + y);
149 for (i = 0; i < 128; i++) {
150 xbeta = kaiserAttn2Beta(attn)
151 if (beta == xbeta || \
152 (i > 63 && \
153 fabs(beta - xbeta) < Z_KAISER_EPSILON))
154 break;
155 if (beta > xbeta)
156 x = attn;
157 else
158 y = attn;
159 attn = 0.5 * (x + y);
160 }
161 }
162
163 return (attn);
164 }
165
166 function kaiserRolloff(len, attn)
167 {
168 return (1.0 - (((1.0 * attn) - 7.95) / (((1.0 * len) - 1.0) * 14.36)));
169 }
170
171 #
172 # 0th order modified Bessel function of the first kind.
173 #
174 function I0(x, s, u, n, h, t)
175 {
176 s = n = u = 1.0;
177 h = x * 0.5;
178
179 do {
180 t = h / n;
181 n += 1.0;
182 t *= t;
183 u *= t;
184 s += u;
185 } while (u >= (I0_EPSILON * s));
186
187 return (s);
188 }
189
190 function wname(beta)
191 {
192 if (beta >= Z_WINDOW_KAISER)
193 return ("Kaiser");
194 else if (beta == Z_WINDOW_BLACKMAN_NUTTALL)
195 return ("Blackman - Nuttall");
196 else if (beta == Z_WINDOW_NUTTALL)
197 return ("Nuttall");
198 else if (beta == Z_WINDOW_BLACKMAN_HARRIS)
199 return ("Blackman - Harris");
200 else if (beta == Z_WINDOW_BLACKMAN)
201 return ("Blackman");
202 else if (beta == Z_WINDOW_HAMMING)
203 return ("Hamming");
204 else if (beta == Z_WINDOW_HANN)
205 return ("Hann");
206 else
207 return ("What The Hell !?!?");
208 }
209
210 function rolloff_round(x)
211 {
212 if (x < 0.67)
213 x = 0.67;
214 else if (x > 1.0)
215 x = 1.0;
216
217 return (x);
218 }
219
220 function tap_round(x, y)
221 {
222 y = floor(x + 3);
223 y -= y % 4;
224 return (y);
225 }
226
227 function lpf(imp, n, rolloff, beta, num, i, j, x, nm, ibeta, w)
228 {
229 rolloff = rolloff_round(rolloff + (Z_NYQUIST_HOVER * (1.0 - rolloff)));
230 imp[0] = rolloff;
231
232 #
233 # Generate ideal sinc impulses, locate the last zero-crossing and pad
234 # the remaining with 0.
235 #
236 # Note that there are other (faster) ways of calculating this without
237 # the misery of traversing the entire sinc given the fact that the
238 # distance between each zero crossings is actually the bandwidth of
239 # the impulses, but it seems having 0.0001% chances of failure due to
240 # limited precision.
241 #
242 j = n;
243 for (i = 1; i < n; i++) {
244 x = (M_PI * i) / (1.0 * num);
245 imp[i] = sin(x * rolloff) / x;
246 if (i != 1 && (imp[i] * imp[i - 1]) <= 0.0)
247 j = i;
248 }
249
250 for (i = j; i < n; i++)
251 imp[i] = 0.0;
252
253 nm = 1.0 * (j - 1);
254
255 if (beta >= Z_WINDOW_KAISER)
256 ibeta = I0(beta);
257
258 for (i = 1; i < j; i++) {
259 if (beta >= Z_WINDOW_KAISER) {
260 # Kaiser window...
261 x = (1.0 * i) / nm;
262 w = I0(beta * sqrt(1.0 - (x * x))) / ibeta;
263 } else {
264 # Cosined windows...
265 x = (M_PI * i) / nm;
266 if (beta == Z_WINDOW_BLACKMAN_NUTTALL) {
267 # Blackman - Nuttall
268 w = 0.36335819 + (0.4891775 * cos(x)) + \
269 (0.1365995 * cos(2 * x)) + \
270 (0.0106411 * cos(3 * x));
271 } else if (beta == Z_WINDOW_NUTTALL) {
272 # Nuttall
273 w = 0.355768 + (0.487396 * cos(x)) + \
274 (0.144232 * cos(2 * x)) + \
275 (0.012604 * cos(3 * x));
276 } else if (beta == Z_WINDOW_BLACKMAN_HARRIS) {
277 # Blackman - Harris
278 w = 0.422323 + (0.49755 * cos(x)) + \
279 (0.07922 * cos(2 * x));
280 } else if (beta == Z_WINDOW_BLACKMAN) {
281 # Blackman
282 w = 0.42 + (0.50 * cos(x)) + \
283 (0.08 * cos(2 * x));
284 } else if (beta == Z_WINDOW_HAMMING) {
285 # Hamming
286 w = 0.54 + (0.46 * cos(x));
287 } else if (beta == Z_WINDOW_HANN) {
288 # Hann
289 w = 0.50 + (0.50 * cos(x));
290 } else {
291 # What The Hell !?!?
292 w = 0.0;
293 }
294 }
295 imp[i] *= w;
296 }
297
298 imp["impulse_length"] = j;
299 imp["rolloff"] = rolloff;
300 }
301
302 function mkfilter(imp, nmult, rolloff, beta, num, \
303 nwing, mwing, nrolloff, i, dcgain, v, quality)
304 {
305 nwing = floor((nmult * num) / 2) + 1;
306
307 lpf(imp, nwing, rolloff, beta, num);
308
309 mwing = imp["impulse_length"];
310 nrolloff = imp["rolloff"];
311 quality = imp["quality"];
312
313 dcgain = 0.0;
314 for (i = num; i < mwing; i += num)
315 dcgain += imp[i];
316 dcgain *= 2.0;
317 dcgain += imp[0];
318
319 for (i = 0; i < nwing; i++)
320 imp[i] /= dcgain;
321
322 if (quality > 2)
323 printf("\n");
324 printf("/*\n");
325 printf(" * quality = %d\n", quality);
326 printf(" * window = %s\n", wname(beta));
327 if (beta >= Z_WINDOW_KAISER) {
328 printf(" * beta: %.2f\n", beta);
329 printf(" * stop: -%.2f dB\n", \
330 kaiserBeta2Attn(beta));
331 }
332 printf(" * length = %d\n", nmult);
333 printf(" * bandwidth = %.2f%%", rolloff * 100.0);
334 if (rolloff != nrolloff) {
335 printf(" + %.2f%% = %.2f%% (nyquist hover: %.2f%%)", \
336 (nrolloff - rolloff) * 100.0, nrolloff * 100.0, \
337 Z_NYQUIST_HOVER * 100.0);
338 }
339 printf("\n");
340 printf(" * drift = %d\n", num);
341 printf(" * width = %d\n", mwing);
342 printf(" */\n");
343 printf("static int32_t z_coeff_q%d[%d] = {", \
344 quality, nwing + (Z_COEFF_OFFSET * 2));
345 for (i = 0; i < (nwing + (Z_COEFF_OFFSET * 2)); i++) {
346 if ((i % 5) == 0)
347 printf("\n ");
348 if (i < Z_COEFF_OFFSET)
349 v = fx_floor(imp[Z_COEFF_OFFSET - i], Z_COEFF_ONE);
350 else if ((i - Z_COEFF_OFFSET) >= nwing)
351 v = fx_floor( \
352 imp[nwing + nwing - i + Z_COEFF_OFFSET - 1],\
353 Z_COEFF_ONE);
354 else
355 v = fx_floor(imp[i - Z_COEFF_OFFSET], Z_COEFF_ONE);
356 printf(" %s0x%08x,", (v < 0) ? "-" : " ", abs(v));
357 }
358 printf("\n};\n\n");
359 printf("/*\n");
360 printf(" * interpolated q%d differences.\n", quality);
361 printf(" */\n");
362 printf("static int32_t z_dcoeff_q%d[%d] = {", quality, nwing);
363 for (i = 1; i <= nwing; i++) {
364 if ((i % 5) == 1)
365 printf("\n ");
366 v = -imp[i - 1];
367 if (i != nwing)
368 v += imp[i];
369 v = fx_floor(v, Z_INTERP_COEFF_ONE);
370 if (abs(v) > abs(largest_interp))
371 largest_interp = v;
372 printf(" %s0x%08x,", (v < 0) ? "-" : " ", abs(v));
373 }
374 printf("\n};\n");
375
376 return (nwing);
377 }
378
379 function filter_parse(s, a, i, attn, alen)
380 {
381 split(s, a, ":");
382 alen = length(a);
383
384 if (alen > 0 && a[1] == "OVERSAMPLING_FACTOR") {
385 if (alen != 2)
386 return (-1);
387 init_drift(floor(a[2]));
388 return (-1);
389 }
390
391 if (alen > 0 && a[1] == "COEFFICIENT_BIT") {
392 if (alen != 2)
393 return (-1);
394 init_coeff_bit(floor(a[2]));
395 return (-1);
396 }
397
398 if (alen > 0 && a[1] == "ACCUMULATOR_BIT") {
399 if (alen != 2)
400 return (-1);
401 init_accum_bit(floor(a[2]));
402 return (-1);
403 }
404
405 if (alen > 0 && a[1] == "INTERPOLATOR") {
406 if (alen != 2)
407 return (-1);
408 init_coeff_interpolator(toupper(a[2]));
409 return (-1);
410 }
411
412 if (alen == 1 || alen == 2) {
413 if (a[1] == "NYQUIST_HOVER") {
414 i = 1.0 * a[2];
415 Z_NYQUIST_HOVER = (i > 0.0 && i < 1.0) ? i : 0.0;
416 return (-1);
417 }
418 i = 1;
419 if (alen == 1) {
420 attn = Z_KAISER_ATTN_DEFAULT;
421 Popts["beta"] = Z_KAISER_BETA_DEFAULT;
422 } else if (Z_WINDOWS[a[1]] < Z_WINDOW_KAISER) {
423 Popts["beta"] = Z_WINDOWS[a[1]];
424 i = tap_round(a[2]);
425 Popts["nmult"] = i;
426 if (i < 28)
427 i = 28;
428 i = 1.0 - (6.44 / i);
429 Popts["rolloff"] = rolloff_round(i);
430 return (0);
431 } else {
432 attn = 1.0 * a[i++];
433 Popts["beta"] = kaiserAttn2Beta(attn);
434 }
435 i = tap_round(a[i]);
436 Popts["nmult"] = i;
437 if (i > 7 && i < 28)
438 i = 27;
439 i = kaiserRolloff(i, attn);
440 Popts["rolloff"] = rolloff_round(i);
441
442 return (0);
443 }
444
445 if (!(alen == 3 || alen == 4))
446 return (-1);
447
448 i = 2;
449
450 if (a[1] == "kaiser") {
451 if (alen > 2)
452 Popts["beta"] = 1.0 * a[i++];
453 else
454 Popts["beta"] = Z_KAISER_BETA_DEFAULT;
455 } else if (Z_WINDOWS[a[1]] < Z_WINDOW_KAISER)
456 Popts["beta"] = Z_WINDOWS[a[1]];
457 else if (1.0 * a[1] < Z_WINDOW_KAISER)
458 return (-1);
459 else
460 Popts["beta"] = kaiserAttn2Beta(1.0 * a[1]);
461 Popts["nmult"] = tap_round(a[i++]);
462 if (a[1] == "kaiser" && alen == 3)
463 i = kaiserRolloff(Popts["nmult"], \
464 kaiserBeta2Attn(Popts["beta"]));
465 else
466 i = 1.0 * a[i];
467 Popts["rolloff"] = rolloff_round(i);
468
469 return (0);
470 }
471
472 function genscale(bit, s1, s2, scale)
473 {
474 if ((bit + Z_COEFF_SHIFT) > Z_ACCUMULATOR_BIT)
475 s1 = Z_COEFF_SHIFT - \
476 (32 - (Z_ACCUMULATOR_BIT - Z_COEFF_SHIFT));
477 else
478 s1 = Z_COEFF_SHIFT - (32 - bit);
479
480 s2 = Z_SHIFT + (32 - bit);
481
482 if (s1 == 0)
483 scale = "v";
484 else if (s1 < 0)
485 scale = sprintf("(v) << %d", abs(s1));
486 else
487 scale = sprintf("(v) >> %d", s1);
488
489 scale = sprintf("(%s) * Z_SCALE_CAST(s)", scale);
490
491 if (s2 != 0)
492 scale = sprintf("(%s) >> %d", scale, s2);
493
494 printf("#define Z_SCALE_%d(v, s)\t%s(%s)\n", \
495 bit, (bit < 10) ? "\t" : "", scale);
496 }
497
498 function genlerp(bit, use64, lerp)
499 {
500 if ((bit + Z_LINEAR_SHIFT) <= 32) {
501 lerp = sprintf("(((y) - (x)) * (z)) >> %d", Z_LINEAR_SHIFT);
502 } else if (use64 != 0) {
503 if ((bit + Z_LINEAR_SHIFT) <= 64) {
504 lerp = sprintf( \
505 "(((int64_t)(y) - (x)) * (z)) " \
506 ">> %d", \
507 Z_LINEAR_SHIFT);
508 } else {
509 lerp = sprintf( \
510 "((int64_t)((y) >> %d) - ((x) >> %d)) * ", \
511 "(z)" \
512 bit + Z_LINEAR_SHIFT - 64, \
513 bit + Z_LINEAR_SHIFT - 64);
514 if ((64 - bit) != 0)
515 lerp = sprintf("(%s) >> %d", lerp, 64 - bit);
516 }
517 } else {
518 lerp = sprintf( \
519 "(((y) >> %d) - ((x) >> %d)) * (z)", \
520 bit + Z_LINEAR_SHIFT - 32, \
521 bit + Z_LINEAR_SHIFT - 32);
522 if ((32 - bit) != 0)
523 lerp = sprintf("(%s) >> %d", lerp, 32 - bit);
524 }
525
526 printf("#define Z_LINEAR_INTERPOLATE_%d(z, x, y)" \
527 "\t\t\t\t%s\\\n\t((x) + (%s))\n", \
528 bit, (bit < 10) ? "\t" : "", lerp);
529 }
530
531 function init_drift(drift, xdrift)
532 {
533 xdrift = floor(drift);
534
535 if (Z_DRIFT_SHIFT != -1) {
536 if (xdrift != Z_DRIFT_SHIFT)
537 printf("#error Z_DRIFT_SHIFT reinitialize!\n");
538 return;
539 }
540
541 #
542 # Initialize filter oversampling factor, or in other word
543 # Z_DRIFT_SHIFT.
544 #
545 if (xdrift < 0)
546 xdrift = 1;
547 else if (xdrift > 31)
548 xdrift = 31;
549
550 Z_DRIFT_SHIFT = xdrift;
551 Z_DRIFT_ONE = shl(1, Z_DRIFT_SHIFT);
552
553 Z_SHIFT = Z_FULL_SHIFT - Z_DRIFT_SHIFT;
554 Z_ONE = shl(1, Z_SHIFT);
555 Z_MASK = Z_ONE - 1;
556 }
557
558 function init_coeff_bit(cbit, xcbit)
559 {
560 xcbit = floor(cbit);
561
562 if (Z_COEFF_SHIFT != 0) {
563 if (xcbit != Z_COEFF_SHIFT)
564 printf("#error Z_COEFF_SHIFT reinitialize!\n");
565 return;
566 }
567
568 #
569 # Initialize dynamic range of coefficients.
570 #
571 if (xcbit < 1)
572 xcbit = 1;
573 else if (xcbit > 30)
574 xcbit = 30;
575
576 Z_COEFF_SHIFT = xcbit;
577 Z_COEFF_ONE = shl(1, Z_COEFF_SHIFT);
578 }
579
580 function init_accum_bit(accbit, xaccbit)
581 {
582 xaccbit = floor(accbit);
583
584 if (Z_ACCUMULATOR_BIT != 0) {
585 if (xaccbit != Z_ACCUMULATOR_BIT)
586 printf("#error Z_ACCUMULATOR_BIT reinitialize!\n");
587 return;
588 }
589
590 #
591 # Initialize dynamic range of accumulator.
592 #
593 if (xaccbit > 64)
594 xaccbit = 64;
595 else if (xaccbit < 32)
596 xaccbit = 32;
597
598 Z_ACCUMULATOR_BIT = xaccbit;
599 }
600
601 function init_coeff_interpolator(interp)
602 {
603 #
604 # Validate interpolator type.
605 #
606 if (interp == "ZOH" || interp == "LINEAR" || \
607 interp == "QUADRATIC" || interp == "HERMITE" || \
608 interp == "BSPLINE" || interp == "OPT32X" || \
609 interp == "OPT16X" || interp == "OPT8X" || \
610 interp == "OPT4X" || interp == "OPT2X")
611 Z_COEFF_INTERPOLATOR = interp;
612 }
613
614 BEGIN {
615 I0_EPSILON = 1e-21;
616 M_PI = atan2(0.0, -1.0);
617
618 INT32_MAX = 1 + ((shl(1, 30) - 1) * 2);
619 INT32_MIN = -1 - INT32_MAX;
620
621 Z_COEFF_OFFSET = 5;
622
623 Z_ACCUMULATOR_BIT_DEFAULT = 58;
624 Z_ACCUMULATOR_BIT = 0;
625
626 Z_FULL_SHIFT = 30;
627 Z_FULL_ONE = shl(1, Z_FULL_SHIFT);
628
629 Z_COEFF_SHIFT_DEFAULT = 30;
630 Z_COEFF_SHIFT = 0;
631 Z_COEFF_ONE = 0;
632
633 Z_COEFF_INTERPOLATOR = 0;
634
635 Z_INTERP_COEFF_SHIFT = 24;
636 Z_INTERP_COEFF_ONE = shl(1, Z_INTERP_COEFF_SHIFT);
637
638 Z_LINEAR_FULL_SHIFT = Z_FULL_SHIFT;
639 Z_LINEAR_FULL_ONE = shl(1, Z_LINEAR_FULL_SHIFT);
640 Z_LINEAR_SHIFT = 8;
641 Z_LINEAR_UNSHIFT = Z_LINEAR_FULL_SHIFT - Z_LINEAR_SHIFT;
642 Z_LINEAR_ONE = shl(1, Z_LINEAR_SHIFT)
643
644 Z_DRIFT_SHIFT_DEFAULT = 5;
645 Z_DRIFT_SHIFT = -1;
646 # meehhhh... let it overflow...
647 #Z_SCALE_SHIFT = 31;
648 #Z_SCALE_ONE = shl(1, Z_SCALE_SHIFT);
649
650 Z_WINDOW_KAISER = 0.0;
651 Z_WINDOW_BLACKMAN_NUTTALL = -1.0;
652 Z_WINDOW_NUTTALL = -2.0;
653 Z_WINDOW_BLACKMAN_HARRIS = -3.0;
654 Z_WINDOW_BLACKMAN = -4.0;
655 Z_WINDOW_HAMMING = -5.0;
656 Z_WINDOW_HANN = -6.0;
657
658 Z_WINDOWS["blackman_nuttall"] = Z_WINDOW_BLACKMAN_NUTTALL;
659 Z_WINDOWS["nuttall"] = Z_WINDOW_NUTTALL;
660 Z_WINDOWS["blackman_harris"] = Z_WINDOW_BLACKMAN_HARRIS;
661 Z_WINDOWS["blackman"] = Z_WINDOW_BLACKMAN;
662 Z_WINDOWS["hamming"] = Z_WINDOW_HAMMING;
663 Z_WINDOWS["hann"] = Z_WINDOW_HANN;
664
665 Z_KAISER_2_BLACKMAN_BETA = 8.568611;
666 Z_KAISER_2_BLACKMAN_NUTTALL_BETA = 11.98;
667
668 Z_KAISER_ATTN_DEFAULT = 100;
669 Z_KAISER_BETA_DEFAULT = kaiserAttn2Beta(Z_KAISER_ATTN_DEFAULT);
670
671 Z_KAISER_EPSILON = 1e-21;
672
673 #
674 # This is practically a joke.
675 #
676 Z_NYQUIST_HOVER = 0.0;
677
678 smallest = 10.000000;
679 largest = 0.000010;
680 largest_interp = 0;
681
682 if (ARGC < 2) {
683 ARGC = 1;
684 ARGV[ARGC++] = "100:8:0.85";
685 ARGV[ARGC++] = "100:36:0.92";
686 ARGV[ARGC++] = "100:164:0.97";
687 #ARGV[ARGC++] = "100:8";
688 #ARGV[ARGC++] = "100:16";
689 #ARGV[ARGC++] = "100:32:0.7929";
690 #ARGV[ARGC++] = "100:64:0.8990";
691 #ARGV[ARGC++] = "100:128:0.9499";
692 }
693
694 printf("#ifndef _FEEDER_RATE_GEN_H_\n");
695 printf("#define _FEEDER_RATE_GEN_H_\n\n");
696 printf("/*\n");
697 printf(" * Generated using feeder_rate_mkfilter.awk, heaven, wind and awesome.\n");
698 printf(" *\n");
699 printf(" * DO NOT EDIT!\n");
700 printf(" */\n\n");
701 printf("#define FEEDER_RATE_PRESETS\t\"");
702 for (i = 1; i < ARGC; i++)
703 printf("%s%s", (i == 1) ? "" : " ", ARGV[i]);
704 printf("\"\n\n");
705 imp["quality"] = 2;
706 for (i = 1; i < ARGC; i++) {
707 if (filter_parse(ARGV[i]) == 0) {
708 beta = Popts["beta"];
709 nmult = Popts["nmult"];
710 rolloff = Popts["rolloff"];
711 if (Z_DRIFT_SHIFT == -1)
712 init_drift(Z_DRIFT_SHIFT_DEFAULT);
713 if (Z_COEFF_SHIFT == 0)
714 init_coeff_bit(Z_COEFF_SHIFT_DEFAULT);
715 if (Z_ACCUMULATOR_BIT == 0)
716 init_accum_bit(Z_ACCUMULATOR_BIT_DEFAULT);
717 ztab[imp["quality"] - 2] = \
718 mkfilter(imp, nmult, rolloff, beta, Z_DRIFT_ONE);
719 imp["quality"]++;
720 }
721 }
722
723 printf("\n");
724 #
725 # XXX
726 #
727 #if (length(ztab) > 0) {
728 # j = 0;
729 # for (i = 0; i < length(ztab); i++) {
730 # if (ztab[i] > j)
731 # j = ztab[i];
732 # }
733 # printf("static int32_t z_coeff_zero[%d] = {", j);
734 # for (i = 0; i < j; i++) {
735 # if ((i % 19) == 0)
736 # printf("\n");
737 # printf(" 0,");
738 # }
739 # printf("\n};\n\n");
740 #}
741 #
742 # XXX
743 #
744 printf("static const struct {\n");
745 printf("\tint32_t len;\n");
746 printf("\tint32_t *coeff;\n");
747 printf("\tint32_t *dcoeff;\n");
748 printf("} z_coeff_tab[] = {\n");
749 if (length(ztab) > 0) {
750 j = 0;
751 for (i = 0; i < length(ztab); i++) {
752 if (ztab[i] > j)
753 j = ztab[i];
754 }
755 j = length(sprintf("%d", j));
756 lfmt = sprintf("%%%dd", j);
757 j = length(sprintf("z_coeff_q%d", length(ztab) + 1));
758 zcfmt = sprintf("%%-%ds", j);
759 zdcfmt = sprintf("%%-%ds", j + 1);
760
761 for (i = 0; i < length(ztab); i++) {
762 l = sprintf(lfmt, ztab[i]);
763 zc = sprintf("z_coeff_q%d", i + 2);
764 zc = sprintf(zcfmt, zc);
765 zdc = sprintf("z_dcoeff_q%d", i + 2);
766 zdc = sprintf(zdcfmt, zdc);
767 printf("\t{ %s, %s, %s },\n", l, zc, zdc);
768 }
769 } else
770 printf("\t{ 0, NULL, NULL }\n");
771 printf("};\n\n");
772
773 #Z_UNSHIFT = 0;
774 #v = shr(Z_ONE - 1, Z_UNSHIFT) * abs(largest_interp);
775 #while (v < 0 || v > INT32_MAX) {
776 # Z_UNSHIFT += 1;
777 # v = shr(Z_ONE - 1, Z_UNSHIFT) * abs(largest_interp);
778 #}
779 v = ((Z_ONE - 1) * abs(largest_interp)) / INT32_MAX;
780 Z_UNSHIFT = ceil(log(v) / log(2.0));
781 Z_INTERP_SHIFT = Z_SHIFT - Z_UNSHIFT + Z_INTERP_COEFF_SHIFT;
782
783 Z_INTERP_UNSHIFT = (Z_SHIFT - Z_UNSHIFT) + Z_INTERP_COEFF_SHIFT \
784 - Z_COEFF_SHIFT;
785
786 printf("#define Z_COEFF_TAB_SIZE\t\t\t\t\t\t\\\n");
787 printf("\t((int32_t)(sizeof(z_coeff_tab) /");
788 printf(" sizeof(z_coeff_tab[0])))\n\n");
789 printf("#define Z_COEFF_OFFSET\t\t%d\n\n", Z_COEFF_OFFSET);
790 printf("#define Z_RSHIFT(x, y)\t\t(((x) + " \
791 "(1 << ((y) - 1))) >> (y))\n");
792 printf("#define Z_RSHIFT_L(x, y)\t(((x) + " \
793 "(1LL << ((y) - 1))) >> (y))\n\n");
794 printf("#define Z_FULL_SHIFT\t\t%d\n", Z_FULL_SHIFT);
795 printf("#define Z_FULL_ONE\t\t0x%08x%s\n", Z_FULL_ONE, \
796 (Z_FULL_ONE > INT32_MAX) ? "U" : "");
797 printf("\n");
798 printf("#define Z_DRIFT_SHIFT\t\t%d\n", Z_DRIFT_SHIFT);
799 #printf("#define Z_DRIFT_ONE\t\t0x%08x\n", Z_DRIFT_ONE);
800 printf("\n");
801 printf("#define Z_SHIFT\t\t\t%d\n", Z_SHIFT);
802 printf("#define Z_ONE\t\t\t0x%08x\n", Z_ONE);
803 printf("#define Z_MASK\t\t\t0x%08x\n", Z_MASK);
804 printf("\n");
805 printf("#define Z_COEFF_SHIFT\t\t%d\n", Z_COEFF_SHIFT);
806 zinterphp = "(z) * (d)";
807 zinterpunshift = Z_SHIFT + Z_INTERP_COEFF_SHIFT - Z_COEFF_SHIFT;
808 if (zinterpunshift > 0) {
809 v = (Z_ONE - 1) * abs(largest_interp);
810 if (v < INT32_MIN || v > INT32_MAX)
811 zinterphp = sprintf("(int64_t)%s", zinterphp);
812 zinterphp = sprintf("(%s) >> %d", zinterphp, zinterpunshift);
813 } else if (zinterpunshift < 0)
814 zinterphp = sprintf("(%s) << %d", zinterphp, \
815 abs(zinterpunshift));
816 if (Z_UNSHIFT == 0)
817 zinterp = "z";
818 else
819 zinterp = sprintf("(z) >> %d", Z_UNSHIFT);
820 zinterp = sprintf("(%s) * (d)", zinterp);
821 if (Z_INTERP_UNSHIFT < 0)
822 zinterp = sprintf("(%s) << %d", zinterp, \
823 abs(Z_INTERP_UNSHIFT));
824 else if (Z_INTERP_UNSHIFT > 0)
825 zinterp = sprintf("(%s) >> %d", zinterp, Z_INTERP_UNSHIFT);
826 if (zinterphp != zinterp) {
827 printf("\n#ifdef SND_FEEDER_RATE_HP\n");
828 printf("#define Z_COEFF_INTERPOLATE(z, c, d)" \
829 "\t\t\t\t\t\\\n\t((c) + (%s))\n", zinterphp);
830 printf("#else\n");
831 printf("#define Z_COEFF_INTERPOLATE(z, c, d)" \
832 "\t\t\t\t\t\\\n\t((c) + (%s))\n", zinterp);
833 printf("#endif\n");
834 } else
835 printf("#define Z_COEFF_INTERPOLATE(z, c, d)" \
836 "\t\t\t\t\t\\\n\t((c) + (%s))\n", zinterp);
837 #printf("\n");
838 #printf("#define Z_SCALE_SHIFT\t\t%d\n", Z_SCALE_SHIFT);
839 #printf("#define Z_SCALE_ONE\t\t0x%08x%s\n", Z_SCALE_ONE, \
840 # (Z_SCALE_ONE > INT32_MAX) ? "U" : "");
841 printf("\n");
842 printf("#define Z_SCALE_CAST(s)\t\t((uint32_t)(s))\n");
843 genscale(8);
844 genscale(16);
845 genscale(24);
846 genscale(32);
847 printf("\n");
848 printf("#define Z_ACCUMULATOR_BIT\t%d\n\n", Z_ACCUMULATOR_BIT)
849 for (i = 8; i <= 32; i += 8) {
850 gbit = ((i + Z_COEFF_SHIFT) > Z_ACCUMULATOR_BIT) ? \
851 (i - (Z_ACCUMULATOR_BIT - Z_COEFF_SHIFT)) : 0;
852 printf("#define Z_GUARD_BIT_%d\t\t%d\n", i, gbit);
853 if (gbit == 0)
854 printf("#define Z_NORM_%d(v)\t\t(v)\n\n", i);
855 else
856 printf("#define Z_NORM_%d(v)\t\t" \
857 "((v) >> Z_GUARD_BIT_%d)\n\n", i, i);
858 }
859 printf("\n");
860 printf("#define Z_LINEAR_FULL_ONE\t0x%08xU\n", Z_LINEAR_FULL_ONE);
861 printf("#define Z_LINEAR_SHIFT\t\t%d\n", Z_LINEAR_SHIFT);
862 printf("#define Z_LINEAR_UNSHIFT\t%d\n", Z_LINEAR_UNSHIFT);
863 printf("#define Z_LINEAR_ONE\t\t0x%08x\n", Z_LINEAR_ONE);
864 printf("\n");
865 printf("#ifdef SND_PCM_64\n");
866 genlerp(8, 1);
867 genlerp(16, 1);
868 genlerp(24, 1);
869 genlerp(32, 1);
870 printf("#else\t/* !SND_PCM_64 */\n");
871 genlerp(8, 0);
872 genlerp(16, 0);
873 genlerp(24, 0);
874 genlerp(32, 0);
875 printf("#endif\t/* SND_PCM_64 */\n");
876 printf("\n");
877 printf("#define Z_QUALITY_ZOH\t\t0\n");
878 printf("#define Z_QUALITY_LINEAR\t1\n");
879 printf("#define Z_QUALITY_SINC\t\t%d\n", \
880 floor((length(ztab) - 1) / 2) + 2);
881 printf("\n");
882 printf("#define Z_QUALITY_MIN\t\t0\n");
883 printf("#define Z_QUALITY_MAX\t\t%d\n", length(ztab) + 1);
884 if (Z_COEFF_INTERPOLATOR != 0)
885 printf("\n#define Z_COEFF_INTERP_%s\t1\n", \
886 Z_COEFF_INTERPOLATOR);
887 printf("\n/*\n * smallest: %.32f\n * largest: %.32f\n *\n", \
888 smallest, largest);
889 printf(" * z_unshift=%d, z_interp_shift=%d\n *\n", \
890 Z_UNSHIFT, Z_INTERP_SHIFT);
891 v = shr(Z_ONE - 1, Z_UNSHIFT) * abs(largest_interp);
892 printf(" * largest interpolation multiplication: %d\n */\n", v);
893 if (v < INT32_MIN || v > INT32_MAX) {
894 printf("\n#ifndef SND_FEEDER_RATE_HP\n");
895 printf("#error interpolation overflow, please reduce" \
896 " Z_INTERP_SHIFT\n");
897 printf("#endif\n");
898 }
899
900 printf("\n#endif /* !_FEEDER_RATE_GEN_H_ */\n");
901 }
Cache object: e925a0d842442fa2d5f31d7f9fb2bb3b
|