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 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 double sqr(double x) {
29 return x*x;
30 }
31
I0(double x)32 static 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
kaiser(int k,int N,double beta)51 static 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
usage(char * name)58 static 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
main(int argc,char ** argv)77 int 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