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 // The number of distinct scale and phase intervals within the filter table.
  45 #define BSINC_SCALE_COUNT (16)
  46 #define BSINC_PHASE_COUNT (16)
  47
  48 #define BSINC_REJECTION (60.0)
  49 #define BSINC_POINTS_MIN (12)
  50
  51 static double MinDouble(double a, double b)
  52 { return (a <= b) ? a : b; }
  53
  54 static double MaxDouble(double a, double b)
  55 { return (a >= b) ? a : b; }
  56
  57 /* NOTE: This is the normalized (instead of just sin(x)/x) cardinal sine
  58  *       function.
  59  *       2 f_t sinc(2 f_t x)
  60  *       f_t -- normalized transition frequency (0.5 is nyquist)
  61  *       x   -- sample index (-N to N)
  62  */
  63 static double Sinc(const double x)
  64 {
  65     if(fabs(x) < 1e-15)
  66         return 1.0;
  67     return sin(M_PI * x) / (M_PI * x);
  68 }
  69
  70 static double BesselI_0(const double x)
  71 {
  72     double term, sum, last_sum, x2, y;
  73     int i;
  74
  75     term = 1.0;
  76     sum = 1.0;
  77     x2 = x / 2.0;
  78     i = 1;
  79
  80     do {
  81         y = x2 / i;
  82         i++;
  83         last_sum = sum;
  84         term *= y * y;
  85         sum += term;
  86     } while(sum != last_sum);
  87
  88     return sum;
  89 }
  90
  91 /* NOTE: k is assumed normalized (-1 to 1)
  92  *       beta is equivalent to 2 alpha
  93  */
  94 static double Kaiser(const double b, const double k)
  95 {
  96     double k2;
  97
  98     if((k < -1.0) || (k > 1.0))
  99         return 0.0;
 100
 101     k2 = MaxDouble(1.0 - (k * k), 0.0);
 102
 103     return BesselI_0(b * sqrt(k2)) / 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 - 1][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 - 1][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_POINTS_MIN);
 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(BSINC_POINTS_MIN, BSINC_POINTS_MIN / (2.0 * scale));
 160         int m = 2 * (int)floor(a);
 161
 162         // Make sure the number of points is a multiple of 4 (for SSE).
 163         m += ~(m - 1) & 3;
 164
 165         mt[si] = m;
 166         at[si] = a;
 167     }
 168
 169     /* Calculate the Kaiser-windowed Sinc filter coefficients for each scale
 170        and phase.
 171     */
 172     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 173     {
 174         const int m = mt[si];
 175         const int o = BSINC_POINTS_MIN - (m / 2);
 176         const int l = (m / 2) - 1;
 177         const double a = at[si];
 178         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 179         const double cutoff = (0.5 * scale) - (scaleBase * MaxDouble(0.5, scale));
 180
 181         for(pi = 0; pi <= BSINC_PHASE_COUNT; pi++)
 182         {
 183             const double phase = l + ((double)pi / BSINC_PHASE_COUNT);
 184
 185             for(i = 0; i < m; i++)
 186             {
 187                 const double x = i - phase;
 188                 filter[si][pi][o + i] = Kaiser(beta, x / a) * 2.0 * cutoff * Sinc(2.0 * cutoff * x);
 189             }
 190         }
 191     }
 192
 193     /* Linear interpolation between scales is simplified by pre-calculating
 194        the delta (b - a) in: x = a + f (b - a)
 195
 196        Given a difference in points between scales, the destination points
 197        will be 0, thus: x = a + f (-a)
 198     */
 199     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 200     {
 201         const int m = mt[si];
 202         const int o = BSINC_POINTS_MIN - (m / 2);
 203
 204         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 205         {
 206             for(i = 0; i < m; i++)
 207                 scDeltas[si][pi][o + i] = filter[si + 1][pi][o + i] - filter[si][pi][o + i];
 208         }
 209     }
 210
 211     // Linear interpolation between phases is also simplified.
 212     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 213     {
 214         const int m = mt[si];
 215         const int o = BSINC_POINTS_MIN - (m / 2);
 216
 217         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 218         {
 219             for(i = 0; i < m; i++)
 220                 phDeltas[si][pi][o + i] = filter[si][pi + 1][o + i] - filter[si][pi][o + i];
 221         }
 222     }
 223
 224     /* This last simplification is done to complete the bilinear equation for
 225        the combination of scale and phase.
 226     */
 227     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 228     {
 229         const int m = mt[si];
 230         const int o = BSINC_POINTS_MIN - (m / 2);
 231
 232         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 233         {
 234             for(i = 0; i < m; i++)
 235                 spDeltas[si][pi][o + i] = phDeltas[si + 1][pi][o + i] - phDeltas[si][pi][o + i];
 236         }
 237     }
 238
 239     // Calculate the table size.
 240     i = mt[0];
 241     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 242         i += BSINC_PHASE_COUNT * mt[si];
 243     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 244         i += 2 * BSINC_PHASE_COUNT * mt[si];
 245     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 246         i += BSINC_PHASE_COUNT * mt[si];
 247
 248     fprintf(stdout, "static const float bsincTab[%d] =\n{\n", i);
 249
 250     /* Only output enough coefficients for the first (cut) scale as needed to
 251        perform interpolation without extra branching.
 252     */
 253     fprintf(stdout, "    /* %2d,%2d */", mt[0], 0);
 254     for(i = 0; i < mt[0]; i++)
 255         fprintf(stdout, " %+14.9ef,", filter[0][0][i]);
 256     fprintf(stdout, "\n\n");
 257
 258     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 259     {
 260         const int m = mt[si];
 261         const int o = BSINC_POINTS_MIN - (m / 2);
 262
 263         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 264         {
 265             fprintf(stdout, "    /* %2d,%2d */", m, pi);
 266             for(i = 0; i < m; i++)
 267                 fprintf(stdout, " %+14.9ef,", filter[si][pi][o + i]);
 268             fprintf(stdout, "\n");
 269         }
 270     }
 271     fprintf(stdout, "\n");
 272
 273     // There are N-1 scale deltas for N scales.
 274     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 275     {
 276         const int m = mt[si];
 277         const int o = BSINC_POINTS_MIN - (m / 2);
 278
 279         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 280         {
 281             fprintf(stdout, "    /* %2d,%2d */", m, pi);
 282             for(i = 0; i < m; i++)
 283                 fprintf(stdout, " %+14.9ef,", scDeltas[si][pi][o + i]);
 284             fprintf(stdout, "\n");
 285         }
 286     }
 287     fprintf(stdout, "\n");
 288
 289     // Exclude phases for the first (cut) scale.
 290     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 291     {
 292         const int m = mt[si];
 293         const int o = BSINC_POINTS_MIN - (m / 2);
 294
 295         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 296         {
 297             fprintf(stdout, "    /* %2d,%2d */", m, pi);
 298             for(i = 0; i < m; i++)
 299                 fprintf(stdout, " %+14.9ef,", phDeltas[si][pi][o + i]);
 300             fprintf(stdout, "\n");
 301         }
 302     }
 303     fprintf(stdout, "\n");
 304
 305     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 306     {
 307         const int m = mt[si];
 308         const int o = BSINC_POINTS_MIN - (m / 2);
 309
 310         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 311         {
 312             fprintf(stdout, "    /* %2d,%2d */", m, pi);
 313             for(i = 0; i < m; i++)
 314                 fprintf(stdout, " %+14.9ef,", spDeltas[si][pi][o + i]);
 315             fprintf(stdout, "\n");
 316         }
 317     }
 318     fprintf(stdout, "};\n\n");
 319
 320     /* The scaleBase is calculated from the Kaiser window transition width.
 321        It represents the absolute limit to the filter before it fully cuts
 322        the signal.  The limit in octaves can be calculated by taking the
 323        base-2 logarithm of its inverse: log_2(1 / scaleBase)
 324     */
 325     fprintf(stdout, "    static const ALfloat scaleBase = %.9ef, scaleRange = %.9ef;\n", scaleBase, 1.0 / scaleRange);
 326     fprintf(stdout, "    static const ALuint m[BSINC_SCALE_COUNT] = {");
 327
 328     fprintf(stdout, " %d", mt[0]);
 329     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 330         fprintf(stdout, ", %d", mt[si]);
 331
 332     fprintf(stdout, " };\n");
 333     fprintf(stdout, "    static const ALuint to[4][BSINC_SCALE_COUNT] =\n    {\n        { 0");
 334
 335     i = mt[0];
 336     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 337     {
 338         fprintf(stdout, ", %d", i);
 339         i += BSINC_PHASE_COUNT * mt[si];
 340     }
 341     fprintf(stdout, " },\n        {");
 342     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 343     {
 344         fprintf(stdout, " %d,", i);
 345         i += BSINC_PHASE_COUNT * mt[si];
 346     }
 347     fprintf(stdout, " 0 },\n        { 0");
 348     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 349     {
 350         fprintf(stdout, ", %d", i);
 351         i += BSINC_PHASE_COUNT * mt[si];
 352     }
 353     fprintf (stdout, " },\n        {");
 354     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 355     {
 356         fprintf(stdout, " %d,", i);
 357         i += BSINC_PHASE_COUNT * mt[si];
 358     }
 359     fprintf(stdout, " 0 }\n    };\n");
 360
 361     fprintf(stdout, "    static const ALuint tm[2][BSINC_SCALE_COUNT] = \n    {\n        { 0");
 362     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 363         fprintf(stdout, ", %d", mt[si]);
 364     fprintf(stdout, " },\n        {");
 365     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
 366         fprintf(stdout, " %d,", mt[si]);
 367     fprintf(stdout, " 0 }\n    };\n\n");
 368 }
 369
 370
 371 /* These methods generate a much simplified 4-point sinc interpolator using a
 372  * Kaiser windows. This is much simpler to process at run-time, but has notably
 373  * more aliasing noise.
 374  */
 375
 376 /* Same as in alu.h! */
 377 #define FRACTIONBITS (12)
 378 #define FRACTIONONE  (1<<FRACTIONBITS)
 379
 380 static double SincKaiser(double r, double x)
 381 {
 382     /* Limit rippling to -60dB. */
 383     return Kaiser(CalcKaiserBeta(60.0), x / r) * Sinc(x);
 384 }
 385
 386 static void Sinc4GenerateTables(void)
 387 {
 388     static double filter[FRACTIONONE][4];
 389
 390     int i;
 391     for(i = 0;i < FRACTIONONE;i++)
 392     {
 393         double mu = (double)i / FRACTIONONE;
 394         filter[i][0] = SincKaiser(2.0, mu - -1.0);
 395         filter[i][1] = SincKaiser(2.0, mu -  0.0);
 396         filter[i][2] = SincKaiser(2.0, mu -  1.0);
 397         filter[i][3] = SincKaiser(2.0, mu -  2.0);
 398     }
 399
 400     fprintf(stdout, "static const float sinc4Tab[%d][4] =\n{\n", FRACTIONONE);
 401     for(i = 0;i < FRACTIONONE;i++)
 402         fprintf(stdout, "    { %+14.9ef, %+14.9ef, %+14.9ef, %+14.9ef },\n",
 403                 filter[i][0], filter[i][1], filter[i][2], filter[i][3]);
 404     fprintf(stdout, "};\n\n");
 405 }
 406
 407 int main(void)
 408 {
 409     BsiGenerateTables();
 410     Sinc4GenerateTables();
 411     return 0;
 412 }