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 // 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 32
  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 phDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][BSINC_POINTS_MAX];
 141     static double scDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT    ][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(phDeltas, 0, sizeof(phDeltas));
 151     memset(scDeltas, 0, sizeof(scDeltas));
 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 phases is simplified by pre-calculating the
 199      * delta (b - a) in: x = a + f (b - a)
 200      */
 201     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 202     {
 203         const int m = mt[si];
 204         const int o = num_points_min - (m / 2);
 205
 206         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 207         {
 208             for(i = 0; i < m; i++)
 209                 phDeltas[si][pi][o + i] = filter[si][pi + 1][o + i] - filter[si][pi][o + i];
 210         }
 211     }
 212
 213     /* Linear interpolation between scales is also simplified.
 214      *
 215      * Given a difference in points between scales, the destination points will
 216      * be 0, thus: x = a + f (-a)
 217      */
 218     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 219     {
 220         const int m = mt[si];
 221         const int o = num_points_min - (m / 2);
 222
 223         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 224         {
 225             for(i = 0; i < m; i++)
 226                 scDeltas[si][pi][o + i] = filter[si + 1][pi][o + i] - filter[si][pi][o + i];
 227         }
 228     }
 229
 230     /* This last simplification is done to complete the bilinear equation for
 231      * the combination of phase and scale.
 232      */
 233     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 234     {
 235         const int m = mt[si];
 236         const int o = num_points_min - (m / 2);
 237
 238         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 239         {
 240             for(i = 0; i < m; i++)
 241                 spDeltas[si][pi][o + i] = phDeltas[si + 1][pi][o + i] - phDeltas[si][pi][o + i];
 242         }
 243     }
 244
 245     /* Make sure the number of points is a multiple of 4 (for SIMD). */
 246     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 247         mt[si] = (mt[si]+3) & ~3;
 248
 249     /* Calculate the table size. */
 250     i = 0;
 251     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 252         i += 4 * BSINC_PHASE_COUNT * mt[si];
 253
 254     fprintf(output,
 255 "/* This %d%s order filter has a rejection of -%.0fdB, yielding a transition width\n"
 256 " * of ~%.3f (normalized frequency). Order increases when downsampling to a\n"
 257 " * limit of one octave, after which the quality of the filter (transition\n"
 258 " * width) suffers to reduce the CPU cost. The bandlimiting will cut all sound\n"
 259 " * after downsampling by ~%.2f octaves.\n"
 260 " */\n"
 261 "alignas(16) constexpr float %s_tab[%d] = {\n",
 262             order, (((order%100)/10) == 1) ? "th" :
 263                    ((order%10) == 1) ? "st" :
 264                    ((order%10) == 2) ? "nd" :
 265                    ((order%10) == 3) ? "rd" : "th",
 266             rejection, width, log2(1.0/scaleBase), tabname, i);
 267     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 268     {
 269         const int m = mt[si];
 270         const int o = num_points_min - (m / 2);
 271
 272         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 273         {
 274             fprintf(output, "    /* %2d,%2d (%d) */", si, pi, m);
 275             fprintf(output, "\n   ");
 276             for(i = 0; i < m; i++)
 277                 fprintf(output, " %+14.9ef,", filter[si][pi][o + i]);
 278             fprintf(output, "\n   ");
 279             for(i = 0; i < m; i++)
 280                 fprintf(output, " %+14.9ef,", phDeltas[si][pi][o + i]);
 281             fprintf(output, "\n   ");
 282             for(i = 0; i < m; i++)
 283                 fprintf(output, " %+14.9ef,", scDeltas[si][pi][o + i]);
 284             fprintf(output, "\n   ");
 285             for(i = 0; i < m; i++)
 286                 fprintf(output, " %+14.9ef,", spDeltas[si][pi][o + i]);
 287             fprintf(output, "\n");
 288         }
 289     }
 290     fprintf(output, "};\nconstexpr BSincTable %s = {\n", tabname);
 291
 292     /* The scaleBase is calculated from the Kaiser window transition width.
 293        It represents the absolute limit to the filter before it fully cuts
 294        the signal.  The limit in octaves can be calculated by taking the
 295        base-2 logarithm of its inverse: log_2(1 / scaleBase)
 296     */
 297     fprintf(output, "    /* scaleBase */ %.9ef, /* scaleRange */ %.9ef,\n", scaleBase, 1.0 / scaleRange);
 298
 299     fprintf(output, "    /* m */ {");
 300     fprintf(output, " %du", mt[0]);
 301     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 302         fprintf(output, ", %du", mt[si]);
 303     fprintf(output, " },\n");
 304
 305     fprintf(output, "    /* filterOffset */ {");
 306     fprintf(output, " %du", 0);
 307     i = mt[0]*4*BSINC_PHASE_COUNT;
 308     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 309     {
 310         fprintf(output, ", %du", i);
 311         i += mt[si]*4*BSINC_PHASE_COUNT;
 312     }
 313
 314     fprintf(output, " },\n");
 315     fprintf(output, "    %s_tab\n", tabname);
 316     fprintf(output, "};\n\n");
 317 }
 318
 319
 320 int main(int argc, char *argv[])
 321 {
 322     FILE *output;
 323
 324     GET_UNICODE_ARGS(&argc, &argv);
 325
 326     if(argc > 2)
 327     {
 328         fprintf(stderr, "Usage: %s [output file]\n", argv[0]);
 329         return 1;
 330     }
 331
 332     if(argc == 2)
 333     {
 334         output = fopen(argv[1], "wb");
 335         if(!output)
 336         {
 337             fprintf(stderr, "Failed to open %s for writing\n", argv[1]);
 338             return 1;
 339         }
 340     }
 341     else
 342         output = stdout;
 343
 344     fprintf(output, "/* Generated by bsincgen, do not edit! */\n"
 345 "#pragma once\n\n"
 346 "static_assert(BSINC_SCALE_COUNT == %d, \"Unexpected BSINC_SCALE_COUNT value!\");\n"
 347 "static_assert(BSINC_PHASE_COUNT == %d, \"Unexpected BSINC_PHASE_COUNT value!\");\n\n"
 348 "namespace {\n\n"
 349 "struct BSincTable {\n"
 350 "    const float scaleBase, scaleRange;\n"
 351 "    const unsigned int m[BSINC_SCALE_COUNT];\n"
 352 "    const unsigned int filterOffset[BSINC_SCALE_COUNT];\n"
 353 "    const float *Tab;\n"
 354 "};\n\n", BSINC_SCALE_COUNT, BSINC_PHASE_COUNT);
 355     /* A 23rd order filter with a -60dB drop at nyquist. */
 356     BsiGenerateTables(output, "bsinc24", 60.0, 23);
 357     /* An 11th order filter with a -60dB drop at nyquist. */
 358     BsiGenerateTables(output, "bsinc12", 60.0, 11);
 359
 360     fprintf(output, "} // namespace\n");
 361     if(output != stdout)
 362         fclose(output);
 363     output = NULL;
 364
 365     return 0;
 366 }