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