1/*
2 * Copyright (C) 2007 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *      http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17#include <math.h>
18#include <stdio.h>
19#include <unistd.h>
20#include <stdlib.h>
21#include <string.h>
22
23static double sinc(double x) {
24    if (fabs(x) == 0.0f) return 1.0f;
25    return sin(x) / x;
26}
27
28static double sqr(double x) {
29    return x*x;
30}
31
32static double I0(double x) {
33    // from the Numerical Recipes in C p. 237
34    double ax,ans,y;
35    ax=fabs(x);
36    if (ax < 3.75) {
37        y=x/3.75;
38        y*=y;
39        ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
40                +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
41    } else {
42        y=3.75/ax;
43        ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
44                +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
45                        +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
46                                +y*0.392377e-2))))))));
47    }
48    return ans;
49}
50
51static double kaiser(int k, int N, double beta) {
52    if (k < 0 || k > N)
53        return 0;
54    return I0(beta * sqrt(1.0 - sqr((2.0*k)/N - 1.0))) / I0(beta);
55}
56
57
58static void usage(char* name) {
59    fprintf(stderr,
60            "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"
61            "       %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"
62            "    -h    this help message\n"
63            "    -d    debug, print comma-separated coefficient table\n"
64            "    -p    generate poly-phase filter coefficients, with sample increment M/N\n"
65            "    -s    sample rate (48000)\n"
66            "    -c    cut-off frequency (20478)\n"
67            "    -n    number of zero-crossings on one side (8)\n"
68            "    -l    number of lerping bits (4)\n"
69            "    -f    output format, can be fixed-point or floating-point (fixed)\n"
70            "    -b    kaiser window parameter beta (7.865 [-80dB])\n"
71            "    -v    attenuation in dBFS (0)\n",
72            name, name
73    );
74    exit(0);
75}
76
77int main(int argc, char** argv)
78{
79    // nc is the number of bits to store the coefficients
80    const int nc = 32;
81
82    bool polyphase = false;
83    unsigned int polyM = 160;
84    unsigned int polyN = 147;
85    bool debug = false;
86    double Fs = 48000;
87    double Fc = 20478;
88    double atten = 1;
89    int format = 0;
90
91
92    // in order to keep the errors associated with the linear
93    // interpolation of the coefficients below the quantization error
94    // we must satisfy:
95    //   2^nz >= 2^(nc/2)
96    //
97    // for 16 bit coefficients that would be 256
98    //
99    // note that increasing nz only increases memory requirements,
100    // but doesn't increase the amount of computation to do.
101    //
102    //
103    // see:
104    // Smith, J.O. Digital Audio Resampling Home Page
105    // https://ccrma.stanford.edu/~jos/resample/, 2011-03-29
106    //
107    int nz = 4;
108
109    //         | 0.1102*(A - 8.7)                         A > 50
110    //  beta = | 0.5842*(A - 21)^0.4 + 0.07886*(A - 21)   21 <= A <= 50
111    //         | 0                                        A < 21
112    //   with A is the desired stop-band attenuation in dBFS
113    //
114    // for eg:
115    //
116    //    30 dB    2.210
117    //    40 dB    3.384
118    //    50 dB    4.538
119    //    60 dB    5.658
120    //    70 dB    6.764
121    //    80 dB    7.865
122    //    90 dB    8.960
123    //   100 dB   10.056
124    double beta = 7.865;
125
126
127    // 2*nzc = (A - 8) / (2.285 * dw)
128    //      with dw the transition width = 2*pi*dF/Fs
129    //
130    int nzc = 8;
131
132    //
133    // Example:
134    // 44.1 KHz to 48 KHz resampling
135    // 100 dB rejection above 28 KHz
136    //   (the spectrum will fold around 24 KHz and we want 100 dB rejection
137    //    at the point where the folding reaches 20 KHz)
138    //  ...___|_____
139    //        |     \|
140    //        | ____/|\____
141    //        |/alias|     \
142    //  ------/------+------\---------> KHz
143    //       20     24     28
144
145    // Transition band 8 KHz, or dw = 1.0472
146    //
147    // beta = 10.056
148    // nzc  = 20
149    //
150
151    int ch;
152    while ((ch = getopt(argc, argv, ":hds:c:n:f:l:b:p:v:")) != -1) {
153        switch (ch) {
154            case 'd':
155                debug = true;
156                break;
157            case 'p':
158                if (sscanf(optarg, "%u/%u", &polyM, &polyN) != 2) {
159                    usage(argv[0]);
160                }
161                polyphase = true;
162                break;
163            case 's':
164                Fs = atof(optarg);
165                break;
166            case 'c':
167                Fc = atof(optarg);
168                break;
169            case 'n':
170                nzc = atoi(optarg);
171                break;
172            case 'l':
173                nz = atoi(optarg);
174                break;
175            case 'f':
176                if (!strcmp(optarg,"fixed")) format = 0;
177                else if (!strcmp(optarg,"float")) format = 1;
178                else usage(argv[0]);
179                break;
180            case 'b':
181                beta = atof(optarg);
182                break;
183            case 'v':
184                atten = pow(10, -fabs(atof(optarg))*0.05 );
185                break;
186            case 'h':
187            default:
188                usage(argv[0]);
189                break;
190        }
191    }
192
193    // cut off frequency ratio Fc/Fs
194    double Fcr = Fc / Fs;
195
196
197    // total number of coefficients (one side)
198    const int M = (1 << nz);
199    const int N = M * nzc;
200
201    // generate the right half of the filter
202    if (!debug) {
203        printf("// cmd-line: ");
204        for (int i=1 ; i<argc ; i++) {
205            printf("%s ", argv[i]);
206        }
207        printf("\n");
208        if (!polyphase) {
209            printf("const int32_t RESAMPLE_FIR_SIZE           = %d;\n", N);
210            printf("const int32_t RESAMPLE_FIR_LERP_INT_BITS  = %d;\n", nz);
211            printf("const int32_t RESAMPLE_FIR_NUM_COEF       = %d;\n", nzc);
212        } else {
213            printf("const int32_t RESAMPLE_FIR_SIZE           = %d;\n", 2*nzc*polyN);
214            printf("const int32_t RESAMPLE_FIR_NUM_COEF       = %d;\n", 2*nzc);
215        }
216        if (!format) {
217            printf("const int32_t RESAMPLE_FIR_COEF_BITS      = %d;\n", nc);
218        }
219        printf("\n");
220        printf("static %s resampleFIR[] = {", !format ? "int32_t" : "float");
221    }
222
223    if (!polyphase) {
224        for (int i=0 ; i<=M ; i++) { // an extra set of coefs for interpolation
225            for (int j=0 ; j<nzc ; j++) {
226                int ix = j*M + i;
227                double x = (2.0 * M_PI * ix * Fcr) / (1 << nz);
228                double y = kaiser(ix+N, 2*N, beta) * sinc(x) * 2.0 * Fcr;
229                y *= atten;
230
231                if (!debug) {
232                    if (j == 0)
233                        printf("\n    ");
234                }
235
236                if (!format) {
237                    int64_t yi = floor(y * ((1ULL<<(nc-1))) + 0.5);
238                    if (yi >= (1LL<<(nc-1))) yi = (1LL<<(nc-1))-1;
239                    printf("0x%08x, ", int32_t(yi));
240                } else {
241                    printf("%.9g%s ", y, debug ? "," : "f,");
242                }
243            }
244        }
245    } else {
246        for (int j=0 ; j<polyN ; j++) {
247            // calculate the phase
248            double p = ((polyM*j) % polyN) / double(polyN);
249            if (!debug) printf("\n    ");
250            else        printf("\n");
251            // generate a FIR per phase
252            for (int i=-nzc ; i<nzc ; i++) {
253                double x = 2.0 * M_PI * Fcr * (i + p);
254                double y = kaiser(i+N, 2*N, beta) * sinc(x) * 2.0 * Fcr;;
255                y *= atten;
256                if (!format) {
257                    int64_t yi = floor(y * ((1ULL<<(nc-1))) + 0.5);
258                    if (yi >= (1LL<<(nc-1))) yi = (1LL<<(nc-1))-1;
259                    printf("0x%08x", int32_t(yi));
260                } else {
261                    printf("%.9g%s", y, debug ? "" : "f");
262                }
263
264                if (debug && (i==nzc-1)) {
265                } else {
266                    printf(", ");
267                }
268            }
269        }
270    }
271
272    if (!debug) {
273        printf("\n};");
274    }
275    printf("\n");
276    return 0;
277}
278
279// http://www.csee.umbc.edu/help/sound/AFsp-V2R1/html/audio/ResampAudio.html
280
281
282