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