From 73e90268adf4c9638b8d820a802e5e9a8ebe6597 Mon Sep 17 00:00:00 2001 From: Pixelflinger Date: Thu, 25 Oct 2012 19:43:49 -0700 Subject: 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. --- tools/resampler_tools/fir.cpp | 263 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 212 insertions(+), 51 deletions(-) (limited to 'tools') 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 #include +#include +#include +#include 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= (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= (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= (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 - + -- cgit v1.1