From 4dd53ab9429f0ca561dfd45f8f6022ae9871455a Mon Sep 17 00:00:00 2001
From: Chris Robinson <chris.kcat@gmail.com>
Date: Tue, 15 Aug 2017 04:15:50 -0700
Subject: Keep bsinc info together in a struct

---
 utils/bsincgen.c | 98 +++++++++++++++++++++++++++++++++++---------------------
 1 file changed, 61 insertions(+), 37 deletions(-)

(limited to 'utils')

diff --git a/utils/bsincgen.c b/utils/bsincgen.c
index 33783a27..fcc5ec85 100644
--- a/utils/bsincgen.c
+++ b/utils/bsincgen.c
@@ -245,16 +245,32 @@ static void BsiGenerateTables()
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
         i += BSINC_PHASE_COUNT * mt[si];
 
-    fprintf(stdout, "static const float bsincTab[%d] =\n{\n", i);
-
+    fprintf(stdout, "/* Generated by bsincgen, do not edit! */\n\n"
+"/* Table of windowed sinc coefficients and deltas.  This 11th order filter\n"
+" * has a rejection of -60 dB, yielding a transition width of ~0.302\n"
+" * (normalized frequency).  Order increases when downsampling to a limit of\n"
+" * one octave, after which the quality of the filter (transition width)\n"
+" * suffers to reduce the CPU cost.  The bandlimiting will cut all sound after\n"
+" * downsampling by ~2.73 octaves.\n"
+" */\n"
+"static const struct {\n"
+"    alignas(16) const float Tab[%d];\n"
+"    const float scaleBase, scaleRange;\n"
+"    const int m[BSINC_SCALE_COUNT];\n"
+"    const int to[4][BSINC_SCALE_COUNT];\n"
+"    const int tm[2][BSINC_SCALE_COUNT];\n"
+"} bsinc = {\n", i);
+
+    fprintf(stdout, "    /* Tab */ {\n");
     /* Only output enough coefficients for the first (cut) scale as needed to
        perform interpolation without extra branching.
     */
-    fprintf(stdout, "    /* %2d,%2d */", mt[0], 0);
+    fprintf(stdout, "        /* %2d,%2d */", mt[0], 0);
     for(i = 0; i < mt[0]; i++)
         fprintf(stdout, " %+14.9ef,", filter[0][0][i]);
     fprintf(stdout, "\n\n");
 
+    fprintf(stdout, "        /* Filters */\n");
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
     {
         const int m = mt[si];
@@ -262,7 +278,7 @@ static void BsiGenerateTables()
 
         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
         {
-            fprintf(stdout, "    /* %2d,%2d */", m, pi);
+            fprintf(stdout, "        /* %2d,%2d */", m, pi);
             for(i = 0; i < m; i++)
                 fprintf(stdout, " %+14.9ef,", filter[si][pi][o + i]);
             fprintf(stdout, "\n");
@@ -271,6 +287,7 @@ static void BsiGenerateTables()
     fprintf(stdout, "\n");
 
     // There are N-1 scale deltas for N scales.
+    fprintf(stdout, "        /* Scale deltas */\n");
     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
     {
         const int m = mt[si];
@@ -278,7 +295,7 @@ static void BsiGenerateTables()
 
         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
         {
-            fprintf(stdout, "    /* %2d,%2d */", m, pi);
+            fprintf(stdout, "        /* %2d,%2d */", m, pi);
             for(i = 0; i < m; i++)
                 fprintf(stdout, " %+14.9ef,", scDeltas[si][pi][o + i]);
             fprintf(stdout, "\n");
@@ -287,6 +304,7 @@ static void BsiGenerateTables()
     fprintf(stdout, "\n");
 
     // Exclude phases for the first (cut) scale.
+    fprintf(stdout, "        /* Phase deltas */\n");
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
     {
         const int m = mt[si];
@@ -294,7 +312,7 @@ static void BsiGenerateTables()
 
         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
         {
-            fprintf(stdout, "    /* %2d,%2d */", m, pi);
+            fprintf(stdout, "        /* %2d,%2d */", m, pi);
             for(i = 0; i < m; i++)
                 fprintf(stdout, " %+14.9ef,", phDeltas[si][pi][o + i]);
             fprintf(stdout, "\n");
@@ -302,6 +320,7 @@ static void BsiGenerateTables()
     }
     fprintf(stdout, "\n");
 
+    fprintf(stdout, "        /* Scale phase deltas */\n");
     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
     {
         const int m = mt[si];
@@ -309,67 +328,67 @@ static void BsiGenerateTables()
 
         for(pi = 0; pi < BSINC_PHASE_COUNT; pi++)
         {
-            fprintf(stdout, "    /* %2d,%2d */", m, pi);
+            fprintf(stdout, "        /* %2d,%2d */", m, pi);
             for(i = 0; i < m; i++)
                 fprintf(stdout, " %+14.9ef,", spDeltas[si][pi][o + i]);
             fprintf(stdout, "\n");
         }
     }
-    fprintf(stdout, "};\n\n");
+    fprintf(stdout, "    },\n\n");
 
     /* The scaleBase is calculated from the Kaiser window transition width.
        It represents the absolute limit to the filter before it fully cuts
        the signal.  The limit in octaves can be calculated by taking the
        base-2 logarithm of its inverse: log_2(1 / scaleBase)
     */
-    fprintf(stdout, "    static const ALfloat scaleBase = %.9ef, scaleRange = %.9ef;\n", scaleBase, 1.0 / scaleRange);
-    fprintf(stdout, "    static const ALuint m[BSINC_SCALE_COUNT] = {");
+    fprintf(stdout, "    /* scaleBase */ %.9ef, /* scaleRange */ %.9ef,\n", scaleBase, 1.0 / scaleRange);
+    fprintf(stdout, "    /* m */ {");
 
     fprintf(stdout, " %d", mt[0]);
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
         fprintf(stdout, ", %d", mt[si]);
 
-    fprintf(stdout, " };\n");
-    fprintf(stdout, "    static const ALuint to[4][BSINC_SCALE_COUNT] =\n    {\n        { 0");
+    fprintf(stdout, " },\n");
+    fprintf(stdout, "    /* to */ {\n        { %5d", 0);
 
     i = mt[0];
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
     {
-        fprintf(stdout, ", %d", i);
+        fprintf(stdout, ", %5d", i);
         i += BSINC_PHASE_COUNT * mt[si];
     }
     fprintf(stdout, " },\n        {");
     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
     {
-        fprintf(stdout, " %d,", i);
+        fprintf(stdout, " %5d,", i);
         i += BSINC_PHASE_COUNT * mt[si];
     }
-    fprintf(stdout, " 0 },\n        { 0");
+    fprintf(stdout, " %5d },\n        { %5d", 0, 0);
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
     {
-        fprintf(stdout, ", %d", i);
+        fprintf(stdout, ", %5d", i);
         i += BSINC_PHASE_COUNT * mt[si];
     }
     fprintf (stdout, " },\n        {");
     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
     {
-        fprintf(stdout, " %d,", i);
+        fprintf(stdout, " %5d,", i);
         i += BSINC_PHASE_COUNT * mt[si];
     }
-    fprintf(stdout, " 0 }\n    };\n");
+    fprintf(stdout, " %5d }\n    },\n", 0);
 
-    fprintf(stdout, "    static const ALuint tm[2][BSINC_SCALE_COUNT] = \n    {\n        { 0");
+    fprintf(stdout, "    /* tm */ {\n        {  0");
     for(si = 1; si < BSINC_SCALE_COUNT; si++)
         fprintf(stdout, ", %d", mt[si]);
     fprintf(stdout, " },\n        {");
     for(si = 0; si < (BSINC_SCALE_COUNT - 1); si++)
         fprintf(stdout, " %d,", mt[si]);
-    fprintf(stdout, " 0 }\n    };\n\n");
+    fprintf(stdout, "  0 }\n    }\n};\n\n");
 }
 
 
 /* These methods generate a much simplified 4-point sinc interpolator using a
- * Kaiser windows. This is much simpler to process at run-time, but has notably
+ * Kaiser window. This is much simpler to process at run-time, but has notably
  * more aliasing noise.
  */
 
@@ -377,30 +396,35 @@ static void BsiGenerateTables()
 #define FRACTIONBITS (12)
 #define FRACTIONONE  (1<<FRACTIONBITS)
 
-static double SincKaiser(double r, double x)
-{
-    /* Limit rippling to -60dB. */
-    return Kaiser(CalcKaiserBeta(60.0), x / r) * Sinc(x);
-}
-
 static void Sinc4GenerateTables(void)
 {
     static double filter[FRACTIONONE][4];
 
-    int i;
-    for(i = 0;i < FRACTIONONE;i++)
+    const double width = CalcKaiserWidth(BSINC_REJECTION, 4);
+    const double beta = CalcKaiserBeta(BSINC_REJECTION);
+    const double scaleBase = width / 2.0;
+    const double scaleRange = 1.0 - scaleBase;
+    const double scale = scaleBase + scaleRange;
+    const double a = MinDouble(4.0, 4.0 / (2.0*scale));
+    const int m = 2 * (int)floor(a);
+    const int l = (m/2) - 1;
+    int pi;
+    for(pi = 0;pi < FRACTIONONE;pi++)
     {
-        double mu = (double)i / FRACTIONONE;
-        filter[i][0] = SincKaiser(2.0, mu - -1.0);
-        filter[i][1] = SincKaiser(2.0, mu -  0.0);
-        filter[i][2] = SincKaiser(2.0, mu -  1.0);
-        filter[i][3] = SincKaiser(2.0, mu -  2.0);
+        const double phase = l + ((double)pi / FRACTIONONE);
+        int i;
+
+        for(i = 0;i < m;i++)
+        {
+            double x = i - phase;
+            filter[pi][i] = Kaiser(beta, x / a) * Sinc(x);
+        }
     }
 
-    fprintf(stdout, "static const float sinc4Tab[%d][4] =\n{\n", FRACTIONONE);
-    for(i = 0;i < FRACTIONONE;i++)
+    fprintf(stdout, "alignas(16) const float sinc4Tab[FRACTIONONE][4] = {\n");
+    for(pi = 0;pi < FRACTIONONE;pi++)
         fprintf(stdout, "    { %+14.9ef, %+14.9ef, %+14.9ef, %+14.9ef },\n",
-                filter[i][0], filter[i][1], filter[i][2], filter[i][3]);
+                filter[pi][0], filter[pi][1], filter[pi][2], filter[pi][3]);
     fprintf(stdout, "};\n\n");
 }
 
-- 
cgit v1.2.3