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