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
sinc(double x)23 static inline double sinc(double x) {
24 if (fabs(x) == 0.0f) return 1.0f;
25 return sin(x) / x;
26 }
27
sqr(double x)28 static inline double sqr(double x) {
29 return x*x;
30 }
31
toint(double x,int64_t maxval)32 static 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
I0(double x)42 static 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
kaiser(int k,int N,double beta)61 static 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
usage(char * name)67 static 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
main(int argc,char ** argv)90 int 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