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