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.
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 @@
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 @@
// 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;
- // total number of coefficients
+
+ // 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 (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");
}
- printf("\n};\n");
- return 0;
- }
-// http://www.dsptutor.freeuk.com/KaiserFilterDesign/KaiserFilterDesign.html
+ 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");
+ return 0;
+}
+
// http://www.csee.umbc.edu/help/sound/AFsp-V2R1/html/audio/ResampAudio.html
-
+