improve fir tool: cleanup, better default, bug fixes

- all parameters can be changed on the command-line
- added float output
- added debug option
- added an option to generate a polyphase filter coefficients
- added an attenuation option in dBFS
- added a lot of comments and references
- fixed kaiser window parameter

also the default should generate a filter with 80 dB rejection
(of the 24 KHz aliasing) above 20 KHz and a 15 KHz transition
band around ~20 KHz (for 48 KHz sampling rate).
It's not very good but corresponds to the current code.

1 files changed, 212 insertions, 51 deletions

diff --git a/tools/resampler_tools/fir.cpp b/tools/resampler_tools/fir.cpp
index 377814f..14707d1 100644
--- a/tools/resampler_tools/fir.cpp
+++ b/tools/resampler_tools/fir.cpp
@@ -16,6 +16,9 @@
 
 #include <math.h>
 #include <stdio.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
 
 static double sinc(double x) {
     if (fabs(x) == 0.0f) return 1.0f;
@@ -34,44 +37,82 @@ static double I0(double x) {
         y=x/3.75;
         y*=y;
         ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
-            +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
+                +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
     } else {
         y=3.75/ax;
         ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
-            +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
-            +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
-            +y*0.392377e-2))))))));
+                +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
+                        +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
+                                +y*0.392377e-2))))))));
     }
     return ans;
 }
 
-static double kaiser(int k, int N, double alpha) {
+static double kaiser(int k, int N, double beta) {
     if (k < 0 || k > N)
         return 0;
-    return I0(M_PI*alpha * sqrt(1.0 - sqr((2.0*k)/N - 1.0))) / I0(M_PI*alpha);
+    return I0(beta * sqrt(1.0 - sqr((2.0*k)/N - 1.0))) / I0(beta);
+}
+
+
+static void usage(char* name) {
+    fprintf(stderr,
+            "usage: %s [-h] [-d] [-s sample_rate] [-c cut-off_frequency] [-n half_zero_crossings] [-f {float|fixed}] [-b beta] [-v dBFS] [-l lerp]\n"
+            "       %s [-h] [-d] [-s sample_rate] [-c cut-off_frequency] [-n half_zero_crossings] [-f {float|fixed}] [-b beta] [-v dBFS] -p M/N\n"
+            "    -h    this help message\n"
+            "    -d    debug, print comma-separated coefficient table\n"
+            "    -p    generate poly-phase filter coefficients, with sample increment M/N\n"
+            "    -s    sample rate (48000)\n"
+            "    -c    cut-off frequency (20478)\n"
+            "    -n    number of zero-crossings on one side (8)\n"
+            "    -l    number of lerping bits (4)\n"
+            "    -f    output format, can be fixed-point or floating-point (fixed)\n"
+            "    -b    kaiser window parameter beta (7.865 [-80dB])\n"
+            "    -v    attenuation in dBFS (0)\n",
+            name, name
+    );
+    exit(0);
 }
 
 int main(int argc, char** argv)
 {
     // nc is the number of bits to store the coefficients
-    int nc = 32;
+    const int nc = 32;
 
-    // ni is the minimum number of bits needed for interpolation
-    // (not used for generating the coefficients)
-    const int ni = nc / 2;
+    bool polyphase = false;
+    unsigned int polyM = 160;
+    unsigned int polyN = 147;
+    bool debug = false;
+    double Fs = 48000;
+    double Fc = 24000;
+    double atten = 1;
+    int format = 0;
 
-    // cut off frequency ratio Fc/Fs
-    // The bigger the stop-band, the less coefficients we'll need.
-    double Fcr = 20000.0 / 48000.0;
 
-    // nzc is the number of zero-crossing on one half of the filter
-    int nzc = 8;
-    
-    // alpha parameter of the kaiser window
-    // Larger numbers reduce ripples in the rejection band but increase
-    // the width of the transition band. 
-    // the table below gives some value of alpha for a given
-    // stop-band attenuation.
+    // in order to keep the errors associated with the linear
+    // interpolation of the coefficients below the quantization error
+    // we must satisfy:
+    //   2^nz >= 2^(nc/2)
+    //
+    // for 16 bit coefficients that would be 256
+    //
+    // note that increasing nz only increases memory requirements,
+    // but doesn't increase the amount of computation to do.
+    //
+    //
+    // see:
+    // Smith, J.O. Digital Audio Resampling Home Page
+    // https://ccrma.stanford.edu/~jos/resample/, 2011-03-29
+    //
+    int nz = 4;
+
+    //         | 0.1102*(A - 8.7)                         A > 50
+    //  beta = | 0.5842*(A - 21)^0.4 + 0.07886*(A - 21)   21 <= A <= 50
+    //         | 0                                        A < 21
+    //   with A is the desired stop-band attenuation in dBFS
+    //
+    // for eg:
+    //
     //    30 dB    2.210
     //    40 dB    3.384
     //    50 dB    4.538
@@ -80,42 +121,162 @@ int main(int argc, char** argv)
     //    80 dB    7.865
     //    90 dB    8.960
     //   100 dB   10.056
-    double alpha = 7.865;	// -80dB stop-band attenuation
-    
-    // 2^nz is the number coefficients per zero-crossing
-    // (int theory this should be 1<<(nc/2))
-    const int nz = 4;
+    double beta = 7.865;
+
+
+    // 2*nzc = (A - 8) / (2.285 * dw)
+    //      with dw the transition width = 2*pi*dF/Fs
+    //
+    int nzc = 8;
+
+    //
+    // Example:
+    // 44.1 KHz to 48 KHz resampling
+    // 100 dB rejection above 28 KHz
+    //   (the spectrum will fold around 24 KHz and we want 100 dB rejection
+    //    at the point where the folding reaches 20 KHz)
+    //  ...___|_____
+    //        |     \|
+    //        | ____/|\____
+    //        |/alias|     \
+    //  ------/------+------\---------> KHz
+    //       20     24     28
+
+    // Transition band 8 KHz, or dw = 1.0472
+    //
+    // beta = 10.056
+    // nzc  = 20
+    //
+
+    int ch;
+    while ((ch = getopt(argc, argv, ":hds:c:n:f:l:b:p:v:")) != -1) {
+        switch (ch) {
+            case 'd':
+                debug = true;
+                break;
+            case 'p':
+                if (sscanf(optarg, "%u/%u", &polyM, &polyN) != 2) {
+                    usage(argv[0]);
+                }
+                polyphase = true;
+                break;
+            case 's':
+                Fs = atof(optarg);
+                break;
+            case 'c':
+                Fc = atof(optarg);
+                break;
+            case 'n':
+                nzc = atoi(optarg);
+                break;
+            case 'l':
+                nz = atoi(optarg);
+                break;
+            case 'f':
+                if (!strcmp(optarg,"fixed")) format = 0;
+                else if (!strcmp(optarg,"float")) format = 1;
+                else usage(argv[0]);
+                break;
+            case 'b':
+                beta = atof(optarg);
+                break;
+            case 'v':
+                atten = pow(10, -fabs(atof(optarg))*0.05 );
+                break;
+            case 'h':
+            default:
+                usage(argv[0]);
+                break;
+        }
+    }
+
+    // cut off frequency ratio Fc/Fs
+    double Fcr = Fc / Fs;
+
 
-    // total number of coefficients
+    // total number of coefficients (one side)
     const int N = (1 << nz) * nzc;
 
     // generate the right half of the filter
-    printf("const int32_t RESAMPLE_FIR_SIZE           = %d;\n", N);
-    printf("const int32_t RESAMPLE_FIR_NUM_COEF       = %d;\n", nzc);
-    printf("const int32_t RESAMPLE_FIR_COEF_BITS      = %d;\n", nc);
-    printf("const int32_t RESAMPLE_FIR_LERP_FRAC_BITS = %d;\n", ni);
-    printf("const int32_t RESAMPLE_FIR_LERP_INT_BITS  = %d;\n", nz);
-    printf("\n");
-    printf("static int16_t resampleFIR[%d] = {", N);
-    for (int i=0 ; i<N ; i++)
-    {
-        double x = (2.0 * M_PI * i * Fcr) / (1 << nz);
-        double y = kaiser(i+N, 2*N, alpha) * sinc(x);
-
-        long yi = floor(y * ((1ULL<<(nc-1))) + 0.5);
-        if (yi >= (1LL<<(nc-1))) yi = (1LL<<(nc-1))-1;        
-
-        if ((i % (1 << 4)) == 0) printf("\n    ");
-        if (nc > 16)
-        	printf("0x%08x, ", int(yi));
-        else 
-        	printf("0x%04x, ", int(yi)&0xFFFF);
+    if (!debug) {
+        printf("// cmd-line: ");
+        for (int i=1 ; i<argc ; i++) {
+            printf("%s ", argv[i]);
+        }
+        printf("\n");
+        if (!polyphase) {
+            printf("const int32_t RESAMPLE_FIR_SIZE           = %d;\n", N);
+            printf("const int32_t RESAMPLE_FIR_LERP_INT_BITS  = %d;\n", nz);
+            printf("const int32_t RESAMPLE_FIR_NUM_COEF       = %d;\n", nzc);
+        } else {
+            printf("const int32_t RESAMPLE_FIR_SIZE           = %d;\n", 2*nzc*polyN);
+            printf("const int32_t RESAMPLE_FIR_NUM_COEF       = %d;\n", 2*nzc);
+        }
+        if (!format) {
+            printf("const int32_t RESAMPLE_FIR_COEF_BITS      = %d;\n", nc);
+        }
+        printf("\n");
+        printf("static %s resampleFIR[] = {", !format ? "int32_t" : "float");
+    }
+
+    if (!polyphase) {
+        for (int i=0 ; i<N ; i++) {
+            double x = (2.0 * M_PI * i * Fcr) / (1 << nz);
+            double y = kaiser(i+N, 2*N, beta) * sinc(x);
+            y *= atten;
+
+            if (!debug) {
+                if ((i % (1<<nz)) == 0)
+                    printf("\n    ");
+            }
+
+            if (!format) {
+                int64_t yi = floor(y * ((1ULL<<(nc-1))) + 0.5);
+                if (yi >= (1LL<<(nc-1))) yi = (1LL<<(nc-1))-1;
+                printf("0x%08x, ", int32_t(yi));
+            } else {
+                printf("%.9g%s ", y, debug ? "," : "f,");
+            }
+        }
+    } else {
+        for (int j=0 ; j<polyN ; j++) {
+            // calculate the phase
+            double p = ((polyM*j) % polyN) / double(polyN);
+            if (!debug) printf("\n    ");
+            else        printf("\n");
+            // generate a FIR per phase
+            for (int i=-nzc ; i<nzc ; i++) {
+                double x = 2.0 * M_PI * Fcr * (i + p);
+                double y = kaiser(i+N, 2*N, beta) * sinc(x);
+                y *= atten;
+                if (!format) {
+                    int64_t yi = floor(y * ((1ULL<<(nc-1))) + 0.5);
+                    if (yi >= (1LL<<(nc-1))) yi = (1LL<<(nc-1))-1;
+                    printf("0x%08x", int32_t(yi));
+                } else {
+                    printf("%.9g%s", y, debug ? "" : "f");
+                }
+
+                if (debug && (i==nzc-1)) {
+                } else {
+                    printf(", ");
+                }
+            }
+        }
+    }
+
+    if (!debug) {
+        if (!format) {
+            printf("\n    0x%08x ", 0);
+        } else {
+            printf("\n    %.9g ", 0.0f);
+        }
+        printf("\n};");
     }
-    printf("\n};\n");
+    printf("\n");
     return 0;
- }
+}
 
-// http://www.dsptutor.freeuk.com/KaiserFilterDesign/KaiserFilterDesign.html
 // http://www.csee.umbc.edu/help/sound/AFsp-V2R1/html/audio/ResampAudio.html
 
- 
+