1#!/usr/bin/awk -f 2# 3# SPDX-License-Identifier: BSD-2-Clause 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# 53function abs(x) 54{ 55 return (((x < 0) ? -x : x) + 0); 56} 57 58function fabs(x) 59{ 60 return (((x < 0.0) ? -x : x) + 0.0); 61} 62 63function ceil(x, r) 64{ 65 r = int(x); 66 if (r < x) 67 r++; 68 return (r + 0); 69} 70 71function floor(x, r) 72{ 73 r = int(x); 74 if (r > x) 75 r--; 76 return (r + 0); 77} 78 79function pow(x, y) 80{ 81 return (exp(1.0 * y * log(1.0 * x))); 82} 83 84# 85# What the hell... 86# 87function shl(x, y) 88{ 89 while (y > 0) { 90 x *= 2; 91 y--; 92 } 93 return (x); 94} 95 96function shr(x, y) 97{ 98 while (y > 0 && x != 0) { 99 x = floor(x / 2); 100 y--; 101 } 102 return (x); 103} 104 105function 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# 122function 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 138function 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 166function 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# 174function 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 190function 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 210function 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 220function tap_round(x, y) 221{ 222 y = floor(x + 3); 223 y -= y % 4; 224 return (y); 225} 226 227function 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 302function 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 379function 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 472function 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 498function 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 531function 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 558function 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 580function 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 601function 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 614BEGIN { 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} 902