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 (11)
  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     double k2;
 102
 103     if((k < -1.0) || (k > 1.0))
 104         return 0.0;
 105
 106     k2 = MaxDouble(1.0 - (k * k), 0.0);
 107
 108     return BesselI_0(b * sqrt(k2)) / BesselI_0(b);
 109 }
 110
 111 /* NOTE: Calculates the transition width of the Kaiser window.  Rejection is
 112  *       in dB.
 113  */
 114 static double CalcKaiserWidth(const double rejection, const int order)
 115 {
 116     double w_t = 2.0 * M_PI;
 117
 118     if(rejection > 21.0)
 119        return (rejection - 7.95) / (order * 2.285 * w_t);
 120
 121     return 5.79 / (order * w_t);
 122 }
 123
 124 static double CalcKaiserBeta(const double rejection)
 125 {
 126     if(rejection > 50.0)
 127        return 0.1102 * (rejection - 8.7);
 128     else if(rejection >= 21.0)
 129        return (0.5842 * pow(rejection - 21.0, 0.4)) +
 130               (0.07886 * (rejection - 21.0));
 131     return 0.0;
 132 }
 133
 134 /* Generates the coefficient, delta, and index tables required by the bsinc resampler */
 135 static void BsiGenerateTables()
 136 {
 137     static double   filter[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][2 * BSINC_POINTS_MIN];
 138     static double scDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT    ][2 * BSINC_POINTS_MIN];
 139     static double phDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT + 1][2 * BSINC_POINTS_MIN];
 140     static double spDeltas[BSINC_SCALE_COUNT][BSINC_PHASE_COUNT    ][2 * BSINC_POINTS_MIN];
 141     static int mt[BSINC_SCALE_COUNT];
 142     static double at[BSINC_SCALE_COUNT];
 143     double width, beta, scaleBase, scaleRange;
 144     int si, pi, i;
 145
 146     memset(filter, 0, sizeof(filter));
 147     memset(scDeltas, 0, sizeof(scDeltas));
 148     memset(phDeltas, 0, sizeof(phDeltas));
 149     memset(spDeltas, 0, sizeof(spDeltas));
 150
 151     /* Calculate windowing parameters.  The width describes the transition
 152        band, but it may vary due to the linear interpolation between scales
 153        of the filter.
 154     */
 155     width = CalcKaiserWidth(BSINC_REJECTION, BSINC_ORDER);
 156     beta = CalcKaiserBeta(BSINC_REJECTION);
 157     scaleBase = width / 2.0;
 158     scaleRange = 1.0 - scaleBase;
 159
 160     // Determine filter scaling.
 161     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 162     {
 163         const double scale = scaleBase + (scaleRange * si / (BSINC_SCALE_COUNT - 1));
 164         const double a = MinDouble(BSINC_POINTS_MIN, BSINC_POINTS_MIN / (2.0 * scale));
 165         int m = 2 * (int)floor(a);
 166
 167         // Make sure the number of points is a multiple of 4 (for SSE).
 168         m += ~(m - 1) & 3;
 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 = BSINC_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 = BSINC_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 = BSINC_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 = BSINC_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     // Calculate the table size.
 245     i = 0;
 246     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 247         i += 4 * BSINC_PHASE_COUNT * mt[si];
 248
 249     fprintf(stdout, "/* Generated by bsincgen, do not edit! */\n\n"
 250 "/* Table of windowed sinc coefficients and deltas.  This %dth order filter\n"
 251 " * has a rejection of -%gdB, yielding a transition width of ~%.3f\n"
 252 " * (normalized frequency).  Order increases when downsampling to a limit of\n"
 253 " * one octave, after which the quality of the filter (transition width)\n"
 254 " * suffers to reduce the CPU cost.  The bandlimiting will cut all sound after\n"
 255 " * downsampling by ~%.2f octaves.\n"
 256 " */\n"
 257 "static const struct {\n"
 258 "    alignas(16) const float Tab[%d];\n"
 259 "    const float scaleBase, scaleRange;\n"
 260 "    const int m[BSINC_SCALE_COUNT];\n"
 261 "    const int filterOffset[BSINC_SCALE_COUNT];\n"
 262 "} bsinc = {\n", BSINC_ORDER, BSINC_REJECTION, width, log2(1.0/scaleBase), i);
 263
 264     fprintf(stdout, "    /* Tab */ {\n");
 265     for(si = 0; si < BSINC_SCALE_COUNT; si++)
 266     {
 267         const int m = mt[si];
 268         const int o = BSINC_POINTS_MIN - (m / 2);
 269
 270         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
 271         {
 272             fprintf(stdout, "        /* %2d,%2d (%d) */", si, pi, m);
 273             fprintf(stdout, "\n       ");
 274             for(i = 0; i < m; i++)
 275                 fprintf(stdout, " %+14.9ef,", filter[si][pi][o + i]);
 276             fprintf(stdout, "\n       ");
 277             for(i = 0; i < m; i++)
 278                 fprintf(stdout, " %+14.9ef,", scDeltas[si][pi][o + i]);
 279             fprintf(stdout, "\n       ");
 280             for(i = 0; i < m; i++)
 281                 fprintf(stdout, " %+14.9ef,", phDeltas[si][pi][o + i]);
 282             fprintf(stdout, "\n       ");
 283             for(i = 0; i < m; i++)
 284                 fprintf(stdout, " %+14.9ef,", spDeltas[si][pi][o + i]);
 285             fprintf(stdout, "\n");
 286         }
 287     }
 288     fprintf(stdout, "    },\n\n");
 289
 290     /* The scaleBase is calculated from the Kaiser window transition width.
 291        It represents the absolute limit to the filter before it fully cuts
 292        the signal.  The limit in octaves can be calculated by taking the
 293        base-2 logarithm of its inverse: log_2(1 / scaleBase)
 294     */
 295     fprintf(stdout, "    /* scaleBase */ %.9ef, /* scaleRange */ %.9ef,\n", scaleBase, 1.0 / scaleRange);
 296
 297     fprintf(stdout, "    /* m */ {");
 298     fprintf(stdout, " %d", mt[0]);
 299     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 300         fprintf(stdout, ", %d", mt[si]);
 301     fprintf(stdout, " },\n");
 302
 303     fprintf(stdout, "    /* filterOffset */ {");
 304     fprintf(stdout, " %d", 0);
 305     i = mt[0]*4*BSINC_PHASE_COUNT;
 306     for(si = 1; si < BSINC_SCALE_COUNT; si++)
 307     {
 308         fprintf(stdout, ", %d", i);
 309         i += mt[si]*4*BSINC_PHASE_COUNT;
 310     }
 311
 312     fprintf(stdout, " }\n};\n\n");
 313 }
 314
 315
 316 /* These methods generate a much simplified 4-point sinc interpolator using a
 317  * Kaiser window. This is much simpler to process at run-time, but has notably
 318  * more aliasing noise.
 319  */
 320
 321 /* Same as in alu.h! */
 322 #define FRACTIONBITS (12)
 323 #define FRACTIONONE  (1<<FRACTIONBITS)
 324
 325 static void Sinc4GenerateTables(void)
 326 {
 327     static double filter[FRACTIONONE][4];
 328
 329     const double width = CalcKaiserWidth(BSINC_REJECTION, 3);
 330     const double beta = CalcKaiserBeta(BSINC_REJECTION);
 331     const double scaleBase = width / 2.0;
 332     const double scaleRange = 1.0 - scaleBase;
 333     const double scale = scaleBase + scaleRange;
 334     const double a = MinDouble(4.0, 4.0 / (2.0*scale));
 335     const int m = 2 * (int)floor(a);
 336     const int l = (m/2) - 1;
 337     int pi;
 338     for(pi = 0;pi < FRACTIONONE;pi++)
 339     {
 340         const double phase = l + ((double)pi / FRACTIONONE);
 341         int i;
 342
 343         for(i = 0;i < m;i++)
 344         {
 345             double x = i - phase;
 346             filter[pi][i] = Kaiser(beta, x / a) * Sinc(x);
 347         }
 348     }
 349
 350     fprintf(stdout, "alignas(16) static const float sinc4Tab[FRACTIONONE][4] = {\n");
 351     for(pi = 0;pi < FRACTIONONE;pi++)
 352         fprintf(stdout, "    { %+14.9ef, %+14.9ef, %+14.9ef, %+14.9ef },\n",
 353                 filter[pi][0], filter[pi][1], filter[pi][2], filter[pi][3]);
 354     fprintf(stdout, "};\n\n");
 355 }
 356
 357 int main(void)
 358 {
 359     BsiGenerateTables();
 360     Sinc4GenerateTables();
 361     return 0;
 362 }