sys/tools/sound/feeder_rate_mkfilter.awk

   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 }