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 #define _UNICODE
  37 #include <stdio.h>
  38 #include <math.h>
  39 #include <string.h>
  40 #include <stdlib.h>
  41
  42 #ifdef _WIN32
  43 #define WIN32_LEAN_AND_MEAN
  44 #include <windows.h>
  45
  46 static char *ToUTF8(const wchar_t *from)
  47 {
  48     char *out = NULL;
  49     int len;
  50     if((len=WideCharToMultiByte(CP_UTF8, 0, from, -1, NULL, 0, NULL, NULL)) > 0)
  51     {
  52         out = calloc(sizeof(*out), len);
  53         WideCharToMultiByte(CP_UTF8, 0, from, -1, out, len, NULL, NULL);
  54         out[len-1] = 0;
  55     }
  56     return out;
  57 }
  58
  59 static WCHAR *FromUTF8(const char *str)
  60 {
  61     WCHAR *out = NULL;
  62     int len;
  63
  64     if((len=MultiByteToWideChar(CP_UTF8, 0, str, -1, NULL, 0)) > 0)
  65     {
  66         out = calloc(sizeof(WCHAR), len);
  67         MultiByteToWideChar(CP_UTF8, 0, str, -1, out, len);
  68         out[len-1] = 0;
  69     }
  70     return out;
  71 }
  72
  73
  74 static FILE *my_fopen(const char *fname, const char *mode)
  75 {
  76     WCHAR *wname=NULL, *wmode=NULL;
  77     FILE *file = NULL;
  78
  79     wname = FromUTF8(fname);
  80     wmode = FromUTF8(mode);
  81     if(!wname)
  82         fprintf(stderr, "Failed to convert UTF-8 filename: \"%s\"\n", fname);
  83     else if(!wmode)
  84         fprintf(stderr, "Failed to convert UTF-8 mode: \"%s\"\n", mode);
  85     else
  86         file = _wfopen(wname, wmode);
  87
  88     free(wname);
  89     free(wmode);
  90
  91     return file;
  92 }
  93 #define fopen my_fopen
  94 #endif
  95
  96
  97 #ifndef M_PI
  98 #define M_PI                         (3.14159265358979323846)
  99 #endif
 100
 101 #if defined(__ANDROID__) && !(defined(_ISOC99_SOURCE) || (defined(_POSIX_C_SOURCE) && _POSIX_C_SOURCE >= 200112L))
 102 #define log2(x)  (log(x) / log(2.0))
 103 #endif
 104
 105 // The number of distinct scale and phase intervals within the filter table.
 106 #define BSINC_SCALE_COUNT (16)
 107 #define BSINC_PHASE_COUNT (16)
 108
 109 /* 24 points includes the doubling for downsampling. So the maximum allowed
 110  * order is 11, which is 12 sample points that multiplied by 2 is 24.
 111  */
 112 #define BSINC_POINTS_MAX (24)
 113
 114 static double MinDouble(double a, double b)
 115 { return (a <= b) ? a : b; }
 116
 117 static double MaxDouble(double a, double b)
 118 { return (a >= b) ? a : b; }
 119
 120 /* NOTE: This is the normalized (instead of just sin(x)/x) cardinal sine
 121  *       function.
 122  *       2 f_t sinc(2 f_t x)
 123  *       f_t -- normalized transition frequency (0.5 is nyquist)
 124  *       x   -- sample index (-N to N)
 125  */
 126 static double Sinc(const double x)
 127 {
 128     if(fabs(x) < 1e-15)
 129         return 1.0;
 130     return sin(M_PI * x) / (M_PI * x);
 131 }
 132
 133 static double BesselI_0(const double x)
 134 {
 135     double term, sum, last_sum, x2, y;
 136     int i;
 137
 138     term = 1.0;
 139     sum = 1.0;
 140     x2 = x / 2.0;
 141     i = 1;
 142
 143     do {
 144         y = x2 / i;
 145         i++;
 146         last_sum = sum;
 147         term *= y * y;
 148         sum += term;
 149     } while(sum != last_sum);
 150
 151     return sum;
 152 }
 153
 154 /* NOTE: k is assumed normalized (-1 to 1)
 155  *       beta is equivalent to 2 alpha
 156  */
 157 static double Kaiser(const double b, const double k)
 158 {
 159     if(!(k >= -1.0 && k <= 1.0))
 160         return 0.0;
 161     return BesselI_0(b * sqrt(1.0 - k*k)) / BesselI_0(b);
 162 }
 163
 164 /* Calculates the (normalized frequency) transition width of the Kaiser window.
 165  * Rejection is in dB.
 166  */
 167 static double CalcKaiserWidth(const double rejection, const int order)
 168 {
 169     double w_t = 2.0 * M_PI;
 170
 171     if(rejection > 21.0)
 172        return (rejection - 7.95) / (order * 2.285 * w_t);
 173     /* This enforces a minimum rejection of just above 21.18dB */
 174     return 5.79 / (order * w_t);
 175 }
 176
 177 static double CalcKaiserBeta(const double rejection)
 178 {
 179     if(rejection > 50.0)
 180        return 0.1102 * (rejection - 8.7);
 181     else if(rejection >= 21.0)
 182        return (0.5842 * pow(rejection - 21.0, 0.4)) +
 183               (0.07886 * (rejection - 21.0));
 184     return 0.0;
 185 }
 186
 187 /* Generates the coefficient, delta, and index tables required by the bsinc resampler */
 188 static void BsiGenerateTables(FILE *output, const char *tabname, const double rejection, const int order)
 189 {
 190     static double   filter[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][BSINC_POINTS_MAX];
 191     static double scDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT    ][BSINC_POINTS_MAX];
 192     static double phDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][BSINC_POINTS_MAX];
 193     static double spDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT    ][BSINC_POINTS_MAX];
 194     static int mt[BSINC_SCALE_COUNT];
 195     static double at[BSINC_SCALE_COUNT];
 196     const int num_points_min = order + 1;
 197     double width, beta, scaleBase, scaleRange;
 198     int si, pi, i;
 199
 200     memset(filter, 0, sizeof(filter));
 201     memset(scDeltas, 0, sizeof(scDeltas));
 202     memset(phDeltas, 0, sizeof(phDeltas));
 203     memset(spDeltas, 0, sizeof(spDeltas));
 204
 205     /* Calculate windowing parameters.  The width describes the transition
 206        band, but it may vary due to the linear interpolation between scales
 207        of the filter.
 208     */
 209     width = CalcKaiserWidth(rejection, order);
 210     beta = CalcKaiserBeta(rejection);
 211     scaleBase = width / 2.0;
 212     scaleRange = 1.0 - scaleBase;
 213
 214     // Determine filter scaling.
 215     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 216     {
 217         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 218         const double a = MinDouble(floor(num_points_min / (2.0 * scale)), num_points_min);
 219         const int m = 2 * (int)a;
 220
 221         mt[si] = m;
 222         at[si] = a;
 223     }
 224
 225     /* Calculate the Kaiser-windowed Sinc filter coefficients for each scale
 226        and phase.
 227     */
 228     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 229     {
 230         const int m = mt[si];
 231         const int o = num_points_min - (m / 2);
 232         const int l = (m / 2) - 1;
 233         const double a = at[si];
 234         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 235         const double cutoff = (0.5 * scale) - (scaleBase * MaxDouble(0.5, scale));
 236
 237         for(pi = 0; pi <= BSINC_PHASE_COUNT; pi++)
 238         {
 239             const double phase = l + ((double)pi / BSINC_PHASE_COUNT);
 240
 241             for(i = 0; i < m; i++)
 242             {
 243                 const double x = i - phase;
 244                 filter[si][pi][o + i] = Kaiser(beta, x / a) * 2.0 * cutoff * Sinc(2.0 * cutoff * x);
 245             }
 246         }
 247     }
 248
 249     /* Linear interpolation between scales is simplified by pre-calculating
 250        the delta (b - a) in: x = a + f (b - a)
 251
 252        Given a difference in points between scales, the destination points
 253        will be 0, thus: x = a + f (-a)
 254     */
 255     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 256     {
 257         const int m = mt[si];
 258         const int o = num_points_min - (m / 2);
 259
 260         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 261         {
 262             for(i = 0; i < m; i++)
 263                 scDeltas[si][pi][o + i] = filter[si + 1][pi][o + i] - filter[si][pi][o + i];
 264         }
 265     }
 266
 267     // Linear interpolation between phases is also simplified.
 268     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 269     {
 270         const int m = mt[si];
 271         const int o = num_points_min - (m / 2);
 272
 273         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 274         {
 275             for(i = 0; i < m; i++)
 276                 phDeltas[si][pi][o + i] = filter[si][pi + 1][o + i] - filter[si][pi][o + i];
 277         }
 278     }
 279
 280     /* This last simplification is done to complete the bilinear equation for
 281        the combination of scale and phase.
 282     */
 283     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 284     {
 285         const int m = mt[si];
 286         const int o = num_points_min - (m / 2);
 287
 288         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 289         {
 290             for(i = 0; i < m; i++)
 291                 spDeltas[si][pi][o + i] = phDeltas[si + 1][pi][o + i] - phDeltas[si][pi][o + i];
 292         }
 293     }
 294
 295     // Make sure the number of points is a multiple of 4 (for SIMD).
 296     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 297         mt[si] = (mt[si]+3) & ~3;
 298
 299     fprintf(output,
 300 "/* This %d%s order filter has a rejection of -%.0fdB, yielding a transition width\n"
 301 " * of ~%.3f (normalized frequency). Order increases when downsampling to a\n"
 302 " * limit of one octave, after which the quality of the filter (transition\n"
 303 " * width) suffers to reduce the CPU cost. The bandlimiting will cut all sound\n"
 304 " * after downsampling by ~%.2f octaves.\n"
 305 " */\n"
 306 "static const BSincTable %s = {\n",
 307             order, (((order%100)/10) == 1) ? "th" :
 308                    ((order%10) == 1) ? "st" :
 309                    ((order%10) == 2) ? "nd" :
 310                    ((order%10) == 3) ? "rd" : "th",
 311             rejection, width, log2(1.0/scaleBase), tabname);
 312
 313     /* The scaleBase is calculated from the Kaiser window transition width.
 314        It represents the absolute limit to the filter before it fully cuts
 315        the signal.  The limit in octaves can be calculated by taking the
 316        base-2 logarithm of its inverse: log_2(1 / scaleBase)
 317     */
 318     fprintf(output, "    /* scaleBase */ %.9ef, /* scaleRange */ %.9ef,\n", scaleBase, 1.0 / scaleRange);
 319
 320     fprintf(output, "    /* m */ {");
 321     fprintf(output, " %d", mt[0]);
 322     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 323         fprintf(output, ", %d", mt[si]);
 324     fprintf(output, " },\n");
 325
 326     fprintf(output, "    /* filterOffset */ {");
 327     fprintf(output, " %d", 0);
 328     i = mt[0]*4*BSINC_PHASE_COUNT;
 329     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 330     {
 331         fprintf(output, ", %d", i);
 332         i += mt[si]*4*BSINC_PHASE_COUNT;
 333     }
 334
 335     fprintf(output, " },\n");
 336
 337     // Calculate the table size.
 338     i = 0;
 339     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 340         i += 4 * BSINC_PHASE_COUNT * mt[si];
 341
 342     fprintf(output, "\n    /* Tab (%d entries) */ {\n", i);
 343     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 344     {
 345         const int m = mt[si];
 346         const int o = num_points_min - (m / 2);
 347
 348         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 349         {
 350             fprintf(output, "        /* %2d,%2d (%d) */", si, pi, m);
 351             fprintf(output, "\n       ");
 352             for(i = 0; i < m; i++)
 353                 fprintf(output, " %+14.9ef,", filter[si][pi][o + i]);
 354             fprintf(output, "\n       ");
 355             for(i = 0; i < m; i++)
 356                 fprintf(output, " %+14.9ef,", scDeltas[si][pi][o + i]);
 357             fprintf(output, "\n       ");
 358             for(i = 0; i < m; i++)
 359                 fprintf(output, " %+14.9ef,", phDeltas[si][pi][o + i]);
 360             fprintf(output, "\n       ");
 361             for(i = 0; i < m; i++)
 362                 fprintf(output, " %+14.9ef,", spDeltas[si][pi][o + i]);
 363             fprintf(output, "\n");
 364         }
 365     }
 366     fprintf(output, "    }\n};\n\n");
 367 }
 368
 369
 370 /* These methods generate a much simplified 4-point sinc interpolator using a
 371  * Kaiser window. This is much simpler to process at run-time, but has notably
 372  * more aliasing noise.
 373  */
 374
 375 /* Same as in alu.h! */
 376 #define FRACTIONBITS (12)
 377 #define FRACTIONONE  (1<<FRACTIONBITS)
 378
 379 static void Sinc4GenerateTables(FILE *output, const double rejection)
 380 {
 381     static double filter[FRACTIONONE][4];
 382
 383     const double width = CalcKaiserWidth(rejection, 3);
 384     const double beta = CalcKaiserBeta(rejection);
 385     const double scaleBase = width / 2.0;
 386     const double scaleRange = 1.0 - scaleBase;
 387     const double scale = scaleBase + scaleRange;
 388     const double a = MinDouble(4.0, floor(4.0 / (2.0*scale)));
 389     const int m = 2 * (int)a;
 390     const int l = (m/2) - 1;
 391     int pi;
 392     for(pi = 0;pi < FRACTIONONE;pi++)
 393     {
 394         const double phase = l + ((double)pi / FRACTIONONE);
 395         int i;
 396
 397         for(i = 0;i < m;i++)
 398         {
 399             double x = i - phase;
 400             filter[pi][i] = Kaiser(beta, x / a) * Sinc(x);
 401         }
 402     }
 403
 404     fprintf(output, "alignas(16) static const float sinc4Tab[FRACTIONONE][4] = {\n");
 405     for(pi = 0;pi < FRACTIONONE;pi++)
 406         fprintf(output, "    { %+14.9ef, %+14.9ef, %+14.9ef, %+14.9ef },\n",
 407                 filter[pi][0], filter[pi][1], filter[pi][2], filter[pi][3]);
 408     fprintf(output, "};\n\n");
 409 }
 410
 411
 412 #ifdef _WIN32
 413 #define main my_main
 414 int main(int argc, char *argv[]);
 415
 416 static char **arglist;
 417 static void cleanup_arglist(void)
 418 {
 419     int i;
 420     for(i = 0;arglist[i];i++)
 421         free(arglist[i]);
 422     free(arglist);
 423 }
 424
 425 int wmain(int argc, const wchar_t *wargv[])
 426 {
 427     int i;
 428
 429     atexit(cleanup_arglist);
 430     arglist = calloc(sizeof(*arglist), argc+1);
 431     for(i = 0;i < argc;i++)
 432         arglist[i] = ToUTF8(wargv[i]);
 433
 434     return main(argc, arglist);
 435 }
 436 #endif
 437
 438 int main(int argc, char *argv[])
 439 {
 440     FILE *output;
 441
 442     if(argc > 2)
 443     {
 444         fprintf(stderr, "Usage: %s [output file]\n", argv[0]);
 445         return 1;
 446     }
 447
 448     if(argc == 2)
 449     {
 450         output = fopen(argv[1], "rb");
 451         if(!output)
 452         {
 453             fprintf(stderr, "Failed to open %s for writing\n", argv[1]);
 454             return 1;
 455         }
 456     }
 457     else
 458         output = stdout;
 459
 460     fprintf(output, "/* Generated by bsincgen, do not edit! */\n\n"
 461 "typedef struct BSincTable {\n"
 462 "    const float scaleBase, scaleRange;\n"
 463 "    const int m[BSINC_SCALE_COUNT];\n"
 464 "    const int filterOffset[BSINC_SCALE_COUNT];\n"
 465 "    alignas(16) const float Tab[];\n"
 466 "} BSincTable;\n\n");
 467     /* An 11th order filter with a -60dB drop at nyquist. */
 468     BsiGenerateTables(output, "bsinc12", 60.0, 11);
 469     Sinc4GenerateTables(output, 60.0);
 470
 471     if(output != stdout)
 472         fclose(output);
 473     output = NULL;
 474
 475     return 0;
 476 }