utils/bsincgen.c

   1 /*
   2  * Sinc interpolator coefficient and delta generator for the OpenAL Soft
   3  * cross platform audio library.
   4  *
   5  * Copyright (C) 2015 by Christopher Fitzgerald.
   6  *
   7  * This library is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Library General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2 of the License, or (at your option) any later version.
  11  *
  12  * This library is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Library General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Library General Public
  18  * License along with this library; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20  * MA 02110-1301  USA
  21  *
  22  * Or visit:  http://www.gnu.org/licenses/old-licenses/lgpl-2.0.html
  23  *
  24  * --------------------------------------------------------------------------
  25  *
  26  * This is a modified version of the bandlimited windowed sinc interpolator
  27  * algorithm presented here:
  28  *
  29  *   Smith, J.O. "Windowed Sinc Interpolation", in
  30  *   Physical Audio Signal Processing,
  31  *   https://ccrma.stanford.edu/~jos/pasp/Windowed_Sinc_Interpolation.html,
  32  *   online book,
  33  *   accessed October 2012.
  34  */
  35
  36 #include <stdio.h>
  37 #include <math.h>
  38 #include <string.h>
  39
  40 #ifndef M_PI
  41 #define M_PI                         (3.14159265358979323846)
  42 #endif
  43
  44 // The number of distinct scale and phase intervals within the filter table.
  45 #define BSINC_SCALE_COUNT (16)
  46 #define BSINC_PHASE_COUNT (16)
  47
  48 #define BSINC_REJECTION (60.0)
  49 #define BSINC_POINTS_MIN (12)
  50
  51 static double MinDouble(double a, double b)
  52 { return (a <= b) ? a : b; }
  53
  54 static double MaxDouble(double a, double b)
  55 { return (a >= b) ? a : b; }
  56
  57 /* NOTE: This is the normalized (instead of just sin(x)/x) cardinal sine
  58  *       function.
  59  *       2 f_t sinc(2 f_t x)
  60  *       f_t -- normalized transition frequency (0.5 is nyquist)
  61  *       x   -- sample index (-N to N)
  62  */
  63 static double Sinc(const double x)
  64 {
  65     if(fabs(x) < 1e-15)
  66         return 1.0;
  67     return sin(M_PI * x) / (M_PI * x);
  68 }
  69
  70 static double BesselI_0(const double x)
  71 {
  72     double term, sum, last_sum, x2, y;
  73     int i;
  74
  75     term = 1.0;
  76     sum = 1.0;
  77     x2 = x / 2.0;
  78     i = 1;
  79
  80     do {
  81         y = x2 / i;
  82         i++;
  83         last_sum = sum;
  84         term *= y * y;
  85         sum += term;
  86     } while(sum != last_sum);
  87
  88     return sum;
  89 }
  90
  91 /* NOTE: k is assumed normalized (-1 to 1)
  92  *       beta is equivalent to 2 alpha
  93  */
  94 static double Kaiser(const double b, const double k)
  95 {
  96     double k2;
  97
  98     if((k < -1.0) || (k > 1.0))
  99         return 0.0;
 100
 101     k2 = MaxDouble(1.0 - (k * k), 0.0);
 102
 103     return BesselI_0(b * sqrt(k2)) / BesselI_0(b);
 104 }
 105
 106 /* NOTE: Calculates the transition width of the Kaiser window.  Rejection is
 107  *       in dB.
 108  */
 109 static double CalcKaiserWidth(const double rejection, const int order)
 110 {
 111     double w_t = 2.0 * M_PI;
 112
 113     if(rejection > 21.0)
 114        return (rejection - 7.95) / (order * 2.285 * w_t);
 115
 116     return 5.79 / (order * w_t);
 117 }
 118
 119 static double CalcKaiserBeta(const double rejection)
 120 {
 121     if(rejection > 50.0)
 122        return 0.1102 * (rejection - 8.7);
 123     else if(rejection >= 21.0)
 124        return (0.5842 * pow(rejection - 21.0, 0.4)) +
 125               (0.07886 * (rejection - 21.0));
 126     return 0.0;
 127 }
 128
 129 /* Generates the coefficient, delta, and index tables required by the bsinc resampler */
 130 static void BsiGenerateTables()
 131 {
 132     static double filter[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][2 * BSINC_POINTS_MIN];
 133     static double scDeltas[BSINC_SCALE_COUNT - 1][BSINC_PHASE_COUNT][2 * BSINC_POINTS_MIN];
 134     static double phDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][2 * BSINC_POINTS_MIN];
 135     static double spDeltas[BSINC_SCALE_COUNT - 1][BSINC_PHASE_COUNT][2 * BSINC_POINTS_MIN];
 136     static int mt[BSINC_SCALE_COUNT];
 137     static double at[BSINC_SCALE_COUNT];
 138     double width, beta, scaleBase, scaleRange;
 139     int si, pi, i;
 140
 141     memset(filter, 0, sizeof(filter));
 142     memset(scDeltas, 0, sizeof(scDeltas));
 143     memset(phDeltas, 0, sizeof(phDeltas));
 144     memset(spDeltas, 0, sizeof(spDeltas));
 145
 146     /* Calculate windowing parameters.  The width describes the transition
 147        band, but it may vary due to the linear interpolation between scales
 148        of the filter.
 149     */
 150     width = CalcKaiserWidth(BSINC_REJECTION, BSINC_POINTS_MIN);
 151     beta = CalcKaiserBeta(BSINC_REJECTION);
 152     scaleBase = width / 2.0;
 153     scaleRange = 1.0 - scaleBase;
 154
 155     // Determine filter scaling.
 156     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 157     {
 158         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 159         const double a = MinDouble(BSINC_POINTS_MIN, BSINC_POINTS_MIN / (2.0 * scale));
 160         int m = 2 * (int)floor(a);
 161
 162         // Make sure the number of points is a multiple of 4 (for SSE).
 163         m += ~(m - 1) & 3;
 164
 165         mt[si] = m;
 166         at[si] = a;
 167     }
 168
 169     /* Calculate the Kaiser-windowed Sinc filter coefficients for each scale
 170        and phase.
 171     */
 172     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 173     {
 174         const int m = mt[si];
 175         const int o = BSINC_POINTS_MIN - (m / 2);
 176         const int l = (m / 2) - 1;
 177         const double a = at[si];
 178         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 179         const double cutoff = (0.5 * scale) - (scaleBase * MaxDouble(0.5, scale));
 180
 181         for(pi = 0; pi <= BSINC_PHASE_COUNT; pi++)
 182         {
 183             const double phase = l + ((double)pi / BSINC_PHASE_COUNT);
 184
 185             for(i = 0; i < m; i++)
 186             {
 187                 const double x = i - phase;
 188                 filter[si][pi][o + i] = Kaiser(beta, x / a) * 2.0 * cutoff * Sinc(2.0 * cutoff * x);
 189             }
 190         }
 191     }
 192
 193     /* Linear interpolation between scales is simplified by pre-calculating
 194        the delta (b - a) in: x = a + f (b - a)
 195
 196        Given a difference in points between scales, the destination points
 197        will be 0, thus: x = a + f (-a)
 198     */
 199     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 200     {
 201         const int m = mt[si];
 202         const int o = BSINC_POINTS_MIN - (m / 2);
 203
 204         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 205         {
 206             for(i = 0; i < m; i++)
 207                 scDeltas[si][pi][o + i] = filter[si + 1][pi][o + i] - filter[si][pi][o + i];
 208         }
 209     }
 210
 211     // Linear interpolation between phases is also simplified.
 212     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 213     {
 214         const int m = mt[si];
 215         const int o = BSINC_POINTS_MIN - (m / 2);
 216
 217         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 218         {
 219             for(i = 0; i < m; i++)
 220                 phDeltas[si][pi][o + i] = filter[si][pi + 1][o + i] - filter[si][pi][o + i];
 221         }
 222     }
 223
 224     /* This last simplification is done to complete the bilinear equation for
 225        the combination of scale and phase.
 226     */
 227     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 228     {
 229         const int m = mt[si];
 230         const int o = BSINC_POINTS_MIN - (m / 2);
 231
 232         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 233         {
 234             for(i = 0; i < m; i++)
 235                 spDeltas[si][pi][o + i] = phDeltas[si + 1][pi][o + i] - phDeltas[si][pi][o + i];
 236         }
 237     }
 238
 239     // Calculate the table size.
 240     i = mt[0];
 241     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 242         i += BSINC_PHASE_COUNT * mt[si];
 243     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 244         i += 2 * BSINC_PHASE_COUNT * mt[si];
 245     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 246         i += BSINC_PHASE_COUNT * mt[si];
 247
 248     fprintf(stdout, "/* Generated by bsincgen, do not edit! */\n\n"
 249 "/* Table of windowed sinc coefficients and deltas.  This 11th order filter\n"
 250 " * has a rejection of -60 dB, yielding a transition width of ~0.302\n"
 251 " * (normalized frequency).  Order increases when downsampling to a limit of\n"
 252 " * one octave, after which the quality of the filter (transition width)\n"
 253 " * suffers to reduce the CPU cost.  The bandlimiting will cut all sound after\n"
 254 " * downsampling by ~2.73 octaves.\n"
 255 " */\n"
 256 "static const struct {\n"
 257 "    alignas(16) const float Tab[%d];\n"
 258 "    const float scaleBase, scaleRange;\n"
 259 "    const int m[BSINC_SCALE_COUNT];\n"
 260 "    const int to[4][BSINC_SCALE_COUNT];\n"
 261 "    const int tm[2][BSINC_SCALE_COUNT];\n"
 262 "} bsinc = {\n", i);
 263
 264     fprintf(stdout, "    /* Tab */ {\n");
 265     /* Only output enough coefficients for the first (cut) scale as needed to
 266        perform interpolation without extra branching.
 267     */
 268     fprintf(stdout, "        /* %2d,%2d */", mt[0], 0);
 269     for(i = 0; i < mt[0]; i++)
 270         fprintf(stdout, " %+14.9ef,", filter[0][0][i]);
 271     fprintf(stdout, "\n\n");
 272
 273     fprintf(stdout, "        /* Filters */\n");
 274     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 275     {
 276         const int m = mt[si];
 277         const int o = BSINC_POINTS_MIN - (m / 2);
 278
 279         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 280         {
 281             fprintf(stdout, "        /* %2d,%2d */", m, pi);
 282             for(i = 0; i < m; i++)
 283                 fprintf(stdout, " %+14.9ef,", filter[si][pi][o + i]);
 284             fprintf(stdout, "\n");
 285         }
 286     }
 287     fprintf(stdout, "\n");
 288
 289     // There are N-1 scale deltas for N scales.
 290     fprintf(stdout, "        /* Scale deltas */\n");
 291     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 292     {
 293         const int m = mt[si];
 294         const int o = BSINC_POINTS_MIN - (m / 2);
 295
 296         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 297         {
 298             fprintf(stdout, "        /* %2d,%2d */", m, pi);
 299             for(i = 0; i < m; i++)
 300                 fprintf(stdout, " %+14.9ef,", scDeltas[si][pi][o + i]);
 301             fprintf(stdout, "\n");
 302         }
 303     }
 304     fprintf(stdout, "\n");
 305
 306     // Exclude phases for the first (cut) scale.
 307     fprintf(stdout, "        /* Phase deltas */\n");
 308     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 309     {
 310         const int m = mt[si];
 311         const int o = BSINC_POINTS_MIN - (m / 2);
 312
 313         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 314         {
 315             fprintf(stdout, "        /* %2d,%2d */", m, pi);
 316             for(i = 0; i < m; i++)
 317                 fprintf(stdout, " %+14.9ef,", phDeltas[si][pi][o + i]);
 318             fprintf(stdout, "\n");
 319         }
 320     }
 321     fprintf(stdout, "\n");
 322
 323     fprintf(stdout, "        /* Scale phase deltas */\n");
 324     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 325     {
 326         const int m = mt[si];
 327         const int o = BSINC_POINTS_MIN - (m / 2);
 328
 329         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 330         {
 331             fprintf(stdout, "        /* %2d,%2d */", m, pi);
 332             for(i = 0; i < m; i++)
 333                 fprintf(stdout, " %+14.9ef,", spDeltas[si][pi][o + i]);
 334             fprintf(stdout, "\n");
 335         }
 336     }
 337     fprintf(stdout, "    },\n\n");
 338
 339     /* The scaleBase is calculated from the Kaiser window transition width.
 340        It represents the absolute limit to the filter before it fully cuts
 341        the signal.  The limit in octaves can be calculated by taking the
 342        base-2 logarithm of its inverse: log_2(1 / scaleBase)
 343     */
 344     fprintf(stdout, "    /* scaleBase */ %.9ef, /* scaleRange */ %.9ef,\n", scaleBase, 1.0 / scaleRange);
 345     fprintf(stdout, "    /* m */ {");
 346
 347     fprintf(stdout, " %d", mt[0]);
 348     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 349         fprintf(stdout, ", %d", mt[si]);
 350
 351     fprintf(stdout, " },\n");
 352     fprintf(stdout, "    /* to */ {\n        { %5d", 0);
 353
 354     i = mt[0];
 355     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 356     {
 357         fprintf(stdout, ", %5d", i);
 358         i += BSINC_PHASE_COUNT * mt[si];
 359     }
 360     fprintf(stdout, " },\n        {");
 361     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 362     {
 363         fprintf(stdout, " %5d,", i);
 364         i += BSINC_PHASE_COUNT * mt[si];
 365     }
 366     fprintf(stdout, " %5d },\n        { %5d", 0, 0);
 367     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 368     {
 369         fprintf(stdout, ", %5d", i);
 370         i += BSINC_PHASE_COUNT * mt[si];
 371     }
 372     fprintf (stdout, " },\n        {");
 373     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 374     {
 375         fprintf(stdout, " %5d,", i);
 376         i += BSINC_PHASE_COUNT * mt[si];
 377     }
 378     fprintf(stdout, " %5d }\n    },\n", 0);
 379
 380     fprintf(stdout, "    /* tm */ {\n        {  0");
 381     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 382         fprintf(stdout, ", %d", mt[si]);
 383     fprintf(stdout, " },\n        {");
 384     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 385         fprintf(stdout, " %d,", mt[si]);
 386     fprintf(stdout, "  0 }\n    }\n};\n\n");
 387 }
 388
 389
 390 /* These methods generate a much simplified 4-point sinc interpolator using a
 391  * Kaiser window. This is much simpler to process at run-time, but has notably
 392  * more aliasing noise.
 393  */
 394
 395 /* Same as in alu.h! */
 396 #define FRACTIONBITS (12)
 397 #define FRACTIONONE  (1<<FRACTIONBITS)
 398
 399 static void Sinc4GenerateTables(void)
 400 {
 401     static double filter[FRACTIONONE][4];
 402
 403     const double width = CalcKaiserWidth(BSINC_REJECTION, 4);
 404     const double beta = CalcKaiserBeta(BSINC_REJECTION);
 405     const double scaleBase = width / 2.0;
 406     const double scaleRange = 1.0 - scaleBase;
 407     const double scale = scaleBase + scaleRange;
 408     const double a = MinDouble(4.0, 4.0 / (2.0*scale));
 409     const int m = 2 * (int)floor(a);
 410     const int l = (m/2) - 1;
 411     int pi;
 412     for(pi = 0;pi < FRACTIONONE;pi++)
 413     {
 414         const double phase = l + ((double)pi / FRACTIONONE);
 415         int i;
 416
 417         for(i = 0;i < m;i++)
 418         {
 419             double x = i - phase;
 420             filter[pi][i] = Kaiser(beta, x / a) * Sinc(x);
 421         }
 422     }
 423
 424     fprintf(stdout, "alignas(16) const float sinc4Tab[FRACTIONONE][4] = {\n");
 425     for(pi = 0;pi < FRACTIONONE;pi++)
 426         fprintf(stdout, "    { %+14.9ef, %+14.9ef, %+14.9ef, %+14.9ef },\n",
 427                 filter[pi][0], filter[pi][1], filter[pi][2], filter[pi][3]);
 428     fprintf(stdout, "};\n\n");
 429 }
 430
 431 int main(void)
 432 {
 433     BsiGenerateTables();
 434     Sinc4GenerateTables();
 435     return 0;
 436 }