| /* | |
| * Copyright (C) 2007 The Android Open Source Project | |
| * | |
| * Licensed under the Apache License, Version 2.0 (the "License"); | |
| * you may not use this file except in compliance with the License. | |
| * You may obtain a copy of the License at | |
| * | |
| * http://www.apache.org/licenses/LICENSE-2.0 | |
| * | |
| * Unless required by applicable law or agreed to in writing, software | |
| * distributed under the License is distributed on an "AS IS" BASIS, | |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | |
| * See the License for the specific language governing permissions and | |
| * limitations under the License. | |
| */ | |
| #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; | |
| return sin(x) / x; | |
| } | |
| static double sqr(double x) { | |
| return x*x; | |
| } | |
| static double I0(double x) { | |
| // from the Numerical Recipes in C p. 237 | |
| double ax,ans,y; | |
| ax=fabs(x); | |
| if (ax < 3.75) { | |
| 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))))); | |
| } 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)))))))); | |
| } | |
| return ans; | |
| } | |
| static double kaiser(int k, int N, double beta) { | |
| if (k < 0 || k > N) | |
| return 0; | |
| 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 | |
| const int nc = 32; | |
| bool polyphase = false; | |
| unsigned int polyM = 160; | |
| unsigned int polyN = 147; | |
| bool debug = false; | |
| double Fs = 48000; | |
| double Fc = 20478; | |
| double atten = 1; | |
| int format = 0; | |
| // 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 | |
| // 60 dB 5.658 | |
| // 70 dB 6.764 | |
| // 80 dB 7.865 | |
| // 90 dB 8.960 | |
| // 100 dB 10.056 | |
| 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 (one side) | |
| const int M = (1 << nz); | |
| const int N = M * nzc; | |
| // generate the right half of the filter | |
| 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<=M ; i++) { // an extra set of coefs for interpolation | |
| for (int j=0 ; j<nzc ; j++) { | |
| int ix = j*M + i; | |
| double x = (2.0 * M_PI * ix * Fcr) / (1 << nz); | |
| double y = kaiser(ix+N, 2*N, beta) * sinc(x) * 2.0 * Fcr; | |
| y *= atten; | |
| if (!debug) { | |
| if (j == 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) * 2.0 * Fcr;; | |
| 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) { | |
| printf("\n};"); | |
| } | |
| printf("\n"); | |
| return 0; | |
| } | |
| // http://www.csee.umbc.edu/help/sound/AFsp-V2R1/html/audio/ResampAudio.html | |