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