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 23 static inline double sinc(double x) { 24 if (fabs(x) == 0.0f) return 1.0f; 25 return sin(x) / x; 26 } 27 28 static inline double sqr(double x) { 29 return x*x; 30 } 31 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 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 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 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 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
