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 #if defined(__ANDROID__) && !(defined(_ISOC99_SOURCE) || (defined(_POSIX_C_SOURCE) && _POSIX_C_SOURCE >= 200112L))
  45 #define log2(x)  (log(x) / log(2.0))
  46 #endif
  47
  48 // The number of distinct scale and phase intervals within the filter table.
  49 #define BSINC_SCALE_COUNT (16)
  50 #define BSINC_PHASE_COUNT (16)
  51
  52 #define BSINC_REJECTION (60.0)
  53 #define BSINC_POINTS_MIN (12)
  54 #define BSINC_ORDER (BSINC_POINTS_MIN - 1)
  55
  56 static double MinDouble(double a, double b)
  57 { return (a <= b) ? a : b; }
  58
  59 static double MaxDouble(double a, double b)
  60 { return (a >= b) ? a : b; }
  61
  62 /* NOTE: This is the normalized (instead of just sin(x)/x) cardinal sine
  63  *       function.
  64  *       2 f_t sinc(2 f_t x)
  65  *       f_t -- normalized transition frequency (0.5 is nyquist)
  66  *       x   -- sample index (-N to N)
  67  */
  68 static double Sinc(const double x)
  69 {
  70     if(fabs(x) < 1e-15)
  71         return 1.0;
  72     return sin(M_PI * x) / (M_PI * x);
  73 }
  74
  75 static double BesselI_0(const double x)
  76 {
  77     double term, sum, last_sum, x2, y;
  78     int i;
  79
  80     term = 1.0;
  81     sum = 1.0;
  82     x2 = x / 2.0;
  83     i = 1;
  84
  85     do {
  86         y = x2 / i;
  87         i++;
  88         last_sum = sum;
  89         term *= y * y;
  90         sum += term;
  91     } while(sum != last_sum);
  92
  93     return sum;
  94 }
  95
  96 /* NOTE: k is assumed normalized (-1 to 1)
  97  *       beta is equivalent to 2 alpha
  98  */
  99 static double Kaiser(const double b, const double k)
 100 {
 101     if(!(k >= -1.0 && k <= 1.0))
 102         return 0.0;
 103     return BesselI_0(b * sqrt(1.0 - k*k)) / 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][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][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_ORDER);
 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(floor(BSINC_POINTS_MIN / (2.0 * scale)),
 160                                    BSINC_POINTS_MIN);
 161         int m = 2 * (int)a;
 162
 163         mt[si] = m;
 164         at[si] = a;
 165     }
 166
 167     /* Calculate the Kaiser-windowed Sinc filter coefficients for each scale
 168        and phase.
 169     */
 170     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 171     {
 172         const int m = mt[si];
 173         const int o = BSINC_POINTS_MIN - (m / 2);
 174         const int l = (m / 2) - 1;
 175         const double a = at[si];
 176         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 177         const double cutoff = (0.5 * scale) - (scaleBase * MaxDouble(0.5, scale));
 178
 179         for(pi = 0; pi <= BSINC_PHASE_COUNT; pi++)
 180         {
 181             const double phase = l + ((double)pi / BSINC_PHASE_COUNT);
 182
 183             for(i = 0; i < m; i++)
 184             {
 185                 const double x = i - phase;
 186                 filter[si][pi][o + i] = Kaiser(beta, x / a) * 2.0 * cutoff * Sinc(2.0 * cutoff * x);
 187             }
 188         }
 189     }
 190
 191     /* Linear interpolation between scales is simplified by pre-calculating
 192        the delta (b - a) in: x = a + f (b - a)
 193
 194        Given a difference in points between scales, the destination points
 195        will be 0, thus: x = a + f (-a)
 196     */
 197     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 198     {
 199         const int m = mt[si];
 200         const int o = BSINC_POINTS_MIN - (m / 2);
 201
 202         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 203         {
 204             for(i = 0; i < m; i++)
 205                 scDeltas[si][pi][o + i] = filter[si + 1][pi][o + i] - filter[si][pi][o + i];
 206         }
 207     }
 208
 209     // Linear interpolation between phases is also simplified.
 210     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 211     {
 212         const int m = mt[si];
 213         const int o = BSINC_POINTS_MIN - (m / 2);
 214
 215         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 216         {
 217             for(i = 0; i < m; i++)
 218                 phDeltas[si][pi][o + i] = filter[si][pi + 1][o + i] - filter[si][pi][o + i];
 219         }
 220     }
 221
 222     /* This last simplification is done to complete the bilinear equation for
 223        the combination of scale and phase.
 224     */
 225     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 226     {
 227         const int m = mt[si];
 228         const int o = BSINC_POINTS_MIN - (m / 2);
 229
 230         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 231         {
 232             for(i = 0; i < m; i++)
 233                 spDeltas[si][pi][o + i] = phDeltas[si + 1][pi][o + i] - phDeltas[si][pi][o + i];
 234         }
 235     }
 236
 237     // Make sure the number of points is a multiple of 4 (for SIMD).
 238     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 239         mt[si] = (mt[si]+3) & ~3;
 240
 241     // Calculate the table size.
 242     i = 0;
 243     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 244         i += 4 * BSINC_PHASE_COUNT * mt[si];
 245
 246     fprintf(stdout, "/* Generated by bsincgen, do not edit! */\n\n"
 247 "/* Table of windowed sinc coefficients and deltas.  This %dth order filter\n"
 248 " * has a rejection of -%gdB, yielding a transition width of ~%.3f\n"
 249 " * (normalized frequency).  Order increases when downsampling to a limit of\n"
 250 " * one octave, after which the quality of the filter (transition width)\n"
 251 " * suffers to reduce the CPU cost.  The bandlimiting will cut all sound after\n"
 252 " * downsampling by ~%.2f octaves.\n"
 253 " */\n"
 254 "static const struct {\n"
 255 "    alignas(16) const float Tab[%d];\n"
 256 "    const float scaleBase, scaleRange;\n"
 257 "    const int m[BSINC_SCALE_COUNT];\n"
 258 "    const int filterOffset[BSINC_SCALE_COUNT];\n"
 259 "} bsinc = {\n", BSINC_ORDER, BSINC_REJECTION, width, log2(1.0/scaleBase), i);
 260
 261     fprintf(stdout, "    /* Tab */ {\n");
 262     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 263     {
 264         const int m = mt[si];
 265         const int o = BSINC_POINTS_MIN - (m / 2);
 266
 267         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 268         {
 269             fprintf(stdout, "        /* %2d,%2d (%d) */", si, pi, m);
 270             fprintf(stdout, "\n       ");
 271             for(i = 0; i < m; i++)
 272                 fprintf(stdout, " %+14.9ef,", filter[si][pi][o + i]);
 273             fprintf(stdout, "\n       ");
 274             for(i = 0; i < m; i++)
 275                 fprintf(stdout, " %+14.9ef,", scDeltas[si][pi][o + i]);
 276             fprintf(stdout, "\n       ");
 277             for(i = 0; i < m; i++)
 278                 fprintf(stdout, " %+14.9ef,", phDeltas[si][pi][o + i]);
 279             fprintf(stdout, "\n       ");
 280             for(i = 0; i < m; i++)
 281                 fprintf(stdout, " %+14.9ef,", spDeltas[si][pi][o + i]);
 282             fprintf(stdout, "\n");
 283         }
 284     }
 285     fprintf(stdout, "    },\n\n");
 286
 287     /* The scaleBase is calculated from the Kaiser window transition width.
 288        It represents the absolute limit to the filter before it fully cuts
 289        the signal.  The limit in octaves can be calculated by taking the
 290        base-2 logarithm of its inverse: log_2(1 / scaleBase)
 291     */
 292     fprintf(stdout, "    /* scaleBase */ %.9ef, /* scaleRange */ %.9ef,\n", scaleBase, 1.0 / scaleRange);
 293
 294     fprintf(stdout, "    /* m */ {");
 295     fprintf(stdout, " %d", mt[0]);
 296     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 297         fprintf(stdout, ", %d", mt[si]);
 298     fprintf(stdout, " },\n");
 299
 300     fprintf(stdout, "    /* filterOffset */ {");
 301     fprintf(stdout, " %d", 0);
 302     i = mt[0]*4*BSINC_PHASE_COUNT;
 303     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 304     {
 305         fprintf(stdout, ", %d", i);
 306         i += mt[si]*4*BSINC_PHASE_COUNT;
 307     }
 308
 309     fprintf(stdout, " }\n};\n\n");
 310 }
 311
 312
 313 /* These methods generate a much simplified 4-point sinc interpolator using a
 314  * Kaiser window. This is much simpler to process at run-time, but has notably
 315  * more aliasing noise.
 316  */
 317
 318 /* Same as in alu.h! */
 319 #define FRACTIONBITS (12)
 320 #define FRACTIONONE  (1<<FRACTIONBITS)
 321
 322 static void Sinc4GenerateTables(void)
 323 {
 324     static double filter[FRACTIONONE][4];
 325
 326     const double width = CalcKaiserWidth(BSINC_REJECTION, 3);
 327     const double beta = CalcKaiserBeta(BSINC_REJECTION);
 328     const double scaleBase = width / 2.0;
 329     const double scaleRange = 1.0 - scaleBase;
 330     const double scale = scaleBase + scaleRange;
 331     const double a = MinDouble(4.0, floor(4.0 / (2.0*scale)));
 332     const int m = 2 * (int)a;
 333     const int l = (m/2) - 1;
 334     int pi;
 335     for(pi = 0;pi < FRACTIONONE;pi++)
 336     {
 337         const double phase = l + ((double)pi / FRACTIONONE);
 338         int i;
 339
 340         for(i = 0;i < m;i++)
 341         {
 342             double x = i - phase;
 343             filter[pi][i] = Kaiser(beta, x / a) * Sinc(x);
 344         }
 345     }
 346
 347     fprintf(stdout, "alignas(16) static const float sinc4Tab[FRACTIONONE][4] = {\n");
 348     for(pi = 0;pi < FRACTIONONE;pi++)
 349         fprintf(stdout, "    { %+14.9ef, %+14.9ef, %+14.9ef, %+14.9ef },\n",
 350                 filter[pi][0], filter[pi][1], filter[pi][2], filter[pi][3]);
 351     fprintf(stdout, "};\n\n");
 352 }
 353
 354 int main(void)
 355 {
 356     BsiGenerateTables();
 357     Sinc4GenerateTables();
 358     return 0;
 359 }