sys/tools/sound/feeder_rate_mkfilter.awk

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