14b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian/* 24cb4f7ce9184de9a221239c28afcf912e7e1ed43Dan Bornstein * Copyright (C) 2007 The Android Open Source Project 34b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * 44b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * Licensed under the Apache License, Version 2.0 (the "License"); 54b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * you may not use this file except in compliance with the License. 64b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * You may obtain a copy of the License at 74b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * 84b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * http://www.apache.org/licenses/LICENSE-2.0 94b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * 104b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * Unless required by applicable law or agreed to in writing, software 114b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * distributed under the License is distributed on an "AS IS" BASIS, 124b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 134b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * See the License for the specific language governing permissions and 144b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian * limitations under the License. 154b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian */ 164b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 174b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian#include <math.h> 184b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian#include <stdio.h> 1973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger#include <unistd.h> 2073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger#include <stdlib.h> 2173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger#include <string.h> 224b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 2386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double sinc(double x) { 244b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian if (fabs(x) == 0.0f) return 1.0f; 254b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian return sin(x) / x; 264b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian} 274b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 2886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double sqr(double x) { 294b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian return x*x; 304b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian} 314b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 3286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline int64_t toint(double x, int64_t maxval) { 3386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung int64_t v; 3486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 35523476a40b40857ca44cb6fab933469cd3598425Andy Hung v = static_cast<int64_t>(floor(x * maxval + 0.5)); 3686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (v >= maxval) { 3786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return maxval - 1; // error! 3886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 3986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return v; 4086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 4186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 424b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopianstatic double I0(double x) { 434b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian // from the Numerical Recipes in C p. 237 444b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian double ax,ans,y; 454b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian ax=fabs(x); 464b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian if (ax < 3.75) { 474b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian y=x/3.75; 484b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian y*=y; 494b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492 5073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2))))); 514b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian } else { 524b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian y=3.75/ax; 534b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1 5473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2 5573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1 5673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger +y*0.392377e-2)))))))); 574b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian } 584b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian return ans; 594b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian} 604b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 6173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflingerstatic double kaiser(int k, int N, double beta) { 624b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian if (k < 0 || k > N) 634b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian return 0; 6473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger return I0(beta * sqrt(1.0 - sqr((2.0*k)/N - 1.0))) / I0(beta); 6573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger} 6673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 6773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflingerstatic void usage(char* name) { 6873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger fprintf(stderr, 696506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten "usage: %s [-h] [-d] [-D] [-s sample_rate] [-c cut-off_frequency] [-n half_zero_crossings]" 7086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung " [-f {float|fixed|fixed16}] [-b beta] [-v dBFS] [-l lerp]\n" 716506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten " %s [-h] [-d] [-D] [-s sample_rate] [-c cut-off_frequency] [-n half_zero_crossings]" 7286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung " [-f {float|fixed|fixed16}] [-b beta] [-v dBFS] -p M/N\n" 7373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -h this help message\n" 7473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -d debug, print comma-separated coefficient table\n" 756506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten " -D generate extra declarations\n" 7673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -p generate poly-phase filter coefficients, with sample increment M/N\n" 7773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -s sample rate (48000)\n" 7873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -c cut-off frequency (20478)\n" 7973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -n number of zero-crossings on one side (8)\n" 8073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -l number of lerping bits (4)\n" 8186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung " -m number of polyphases (related to -l, default 16)\n" 82904e982b8ec64df6646a61823199ce22c0674826Glenn Kasten " -f output format, can be fixed, fixed16, or float (fixed)\n" 8373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -b kaiser window parameter beta (7.865 [-80dB])\n" 8473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger " -v attenuation in dBFS (0)\n", 8573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger name, name 8673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger ); 8773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger exit(0); 884b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian} 894b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 904b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopianint main(int argc, char** argv) 914b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian{ 924b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian // nc is the number of bits to store the coefficients 9386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung int nc = 32; 9473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger bool polyphase = false; 9573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger unsigned int polyM = 160; 9673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger unsigned int polyN = 147; 9773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger bool debug = false; 9873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger double Fs = 48000; 99b4b75b47c2a4248e60bbc3229d6acc4d5f872431Mathias Agopian double Fc = 20478; 10073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger double atten = 1; 1016506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten int format = 0; // 0=fixed, 1=float 102bf5d0666b803921541e6a5d6abae05316e533855Glenn Kasten bool declarations = false; 1034b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 10473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // in order to keep the errors associated with the linear 10573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // interpolation of the coefficients below the quantization error 10673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // we must satisfy: 10773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 2^nz >= 2^(nc/2) 10873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 10973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // for 16 bit coefficients that would be 256 11073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 11173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // note that increasing nz only increases memory requirements, 11273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // but doesn't increase the amount of computation to do. 11373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 11473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 11573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // see: 11673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // Smith, J.O. Digital Audio Resampling Home Page 11773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // https://ccrma.stanford.edu/~jos/resample/, 2011-03-29 11873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 11973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 12073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // | 0.1102*(A - 8.7) A > 50 12173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // beta = | 0.5842*(A - 21)^0.4 + 0.07886*(A - 21) 21 <= A <= 50 12273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // | 0 A < 21 12373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // with A is the desired stop-band attenuation in dBFS 12473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 12573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // for eg: 12673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 1272a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 30 dB 2.210 1282a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 40 dB 3.384 1292a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 50 dB 4.538 1302a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 60 dB 5.658 1312a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 70 dB 6.764 1322a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 80 dB 7.865 1332a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 90 dB 8.960 1342a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian // 100 dB 10.056 13573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger double beta = 7.865; 13673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 13773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 2*nzc = (A - 8) / (2.285 * dw) 13873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // with dw the transition width = 2*pi*dF/Fs 13973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // 14073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger int nzc = 8; 14173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 142a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten /* 143a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * Example: 144a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * 44.1 KHz to 48 KHz resampling 145a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * 100 dB rejection above 28 KHz 146a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * (the spectrum will fold around 24 KHz and we want 100 dB rejection 147a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * at the point where the folding reaches 20 KHz) 148a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * ...___|_____ 149a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * | \| 150a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * | ____/|\____ 151a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * |/alias| \ 152a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * ------/------+------\---------> KHz 153a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * 20 24 28 154a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * 155a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * Transition band 8 KHz, or dw = 1.0472 156a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * 157a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * beta = 10.056 158a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten * nzc = 20 159a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten */ 16073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 16186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung int M = 1 << 4; // number of phases for interpolation 16273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger int ch; 1636506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten while ((ch = getopt(argc, argv, ":hds:c:n:f:l:m:b:p:v:z:D")) != -1) { 16473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger switch (ch) { 16573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'd': 16673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger debug = true; 16773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 1686506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten case 'D': 1696506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten declarations = true; 1706506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten break; 17173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'p': 17273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger if (sscanf(optarg, "%u/%u", &polyM, &polyN) != 2) { 17373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger usage(argv[0]); 17473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 17573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger polyphase = true; 17673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 17773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 's': 17873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger Fs = atof(optarg); 17973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 18073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'c': 18173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger Fc = atof(optarg); 18273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 18373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'n': 18473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger nzc = atoi(optarg); 18573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 18686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung case 'm': 18786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung M = atoi(optarg); 18886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung break; 18973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'l': 19086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung M = 1 << atoi(optarg); 19173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 19273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'f': 19386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (!strcmp(optarg, "fixed")) { 19486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung format = 0; 19586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 19686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung else if (!strcmp(optarg, "fixed16")) { 19786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung format = 0; 19886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung nc = 16; 19986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 20086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung else if (!strcmp(optarg, "float")) { 20186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung format = 1; 20286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 20386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung else { 20486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung usage(argv[0]); 20586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 20673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 20773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'b': 20873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger beta = atof(optarg); 20973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 21073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'v': 21173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger atten = pow(10, -fabs(atof(optarg))*0.05 ); 21273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 21373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger case 'h': 21473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger default: 21573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger usage(argv[0]); 21673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger break; 21773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 21873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 21973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 22073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // cut off frequency ratio Fc/Fs 22173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger double Fcr = Fc / Fs; 22273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 22373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // total number of coefficients (one side) 22486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 22546afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian const int N = M * nzc; 2264b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian 22786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // lerp (which is most useful if M is a power of 2) 22886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 22986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung int nz = 0; // recalculate nz as the bits needed to represent M 23086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung for (int i = M-1 ; i; i>>=1, nz++); 2314b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian // generate the right half of the filter 23273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger if (!debug) { 233bf5d0666b803921541e6a5d6abae05316e533855Glenn Kasten printf("// cmd-line:"); 234393d4d2082351b0f63ff19fbe82308481e508b3eGlenn Kasten for (int i=0 ; i<argc ; i++) { 235bf5d0666b803921541e6a5d6abae05316e533855Glenn Kasten printf(" %s", argv[i]); 23673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 23773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger printf("\n"); 2386506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten if (declarations) { 2396506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten if (!polyphase) { 2406506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("const int32_t RESAMPLE_FIR_SIZE = %d;\n", N); 2416506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("const int32_t RESAMPLE_FIR_INT_PHASES = %d;\n", M); 2426506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("const int32_t RESAMPLE_FIR_NUM_COEF = %d;\n", nzc); 2436506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten } else { 2446506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("const int32_t RESAMPLE_FIR_SIZE = %d;\n", 2*nzc*polyN); 2456506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("const int32_t RESAMPLE_FIR_NUM_COEF = %d;\n", 2*nzc); 2466506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten } 2476506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten if (!format) { 2486506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("const int32_t RESAMPLE_FIR_COEF_BITS = %d;\n", nc); 2496506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten } 2506506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("\n"); 2516506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("static %s resampleFIR[] = {", !format ? "int32_t" : "float"); 25273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 25373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 25473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 25573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger if (!polyphase) { 25646afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian for (int i=0 ; i<=M ; i++) { // an extra set of coefs for interpolation 25746afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian for (int j=0 ; j<nzc ; j++) { 25846afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian int ix = j*M + i; 25986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double x = (2.0 * M_PI * ix * Fcr) / M; 26046afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian double y = kaiser(ix+N, 2*N, beta) * sinc(x) * 2.0 * Fcr; 26146afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian y *= atten; 26273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 26346afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian if (!debug) { 26446afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian if (j == 0) 265e1984fc72d9f814461992b9f5b1e5f9d45ce7afaGlenn Kasten printf("\n "); 26646afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian } 26746afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian if (!format) { 26886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung int64_t yi = toint(y, 1ULL<<(nc-1)); 26986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (nc > 16) { 2709c380ac229c5815437819b346fb468655a423747Glenn Kasten printf("0x%08x,", int32_t(yi)); 27186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } else { 2729c380ac229c5815437819b346fb468655a423747Glenn Kasten printf("0x%04x,", int32_t(yi)&0xffff); 27386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 27446afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian } else { 2759c380ac229c5815437819b346fb468655a423747Glenn Kasten printf("%.9g%s", y, debug ? "," : "f,"); 2769c380ac229c5815437819b346fb468655a423747Glenn Kasten } 2779c380ac229c5815437819b346fb468655a423747Glenn Kasten if (j != nzc-1) { 2789c380ac229c5815437819b346fb468655a423747Glenn Kasten printf(" "); 27946afbec3743f1d799f185273ff897d1f8e0175ddMathias Agopian } 28073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 28173e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 28273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } else { 283a4ebf1324ac7d8fa1a5ec77bd7f27d124ede67d6Glenn Kasten for (unsigned int j=0 ; j<polyN ; j++) { 28473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // calculate the phase 28573e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger double p = ((polyM*j) % polyN) / double(polyN); 28673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger if (!debug) printf("\n "); 28773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger else printf("\n"); 28873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger // generate a FIR per phase 28973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger for (int i=-nzc ; i<nzc ; i++) { 29073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger double x = 2.0 * M_PI * Fcr * (i + p); 291d88a051aff15fdf5c57e1e5a4083bbd9635af3adMathias Agopian double y = kaiser(i+N, 2*N, beta) * sinc(x) * 2.0 * Fcr;; 29273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger y *= atten; 29373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger if (!format) { 29486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung int64_t yi = toint(y, 1ULL<<(nc-1)); 29586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (nc > 16) { 2966506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("0x%08x,", int32_t(yi)); 29786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } else { 2986506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("0x%04x,", int32_t(yi)&0xffff); 29986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 30073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } else { 3016506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf("%.9g%s", y, debug ? "," : "f,"); 30273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 30373e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 3046506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten if (i != nzc-1) { 3056506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten printf(" "); 30673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 30773e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 30873e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 30973e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger } 31073e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger 3116506d19004ce9d6943b80721c37e26532c6e6bdfGlenn Kasten if (!debug && declarations) { 31273e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger printf("\n};"); 3134b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian } 31473e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger printf("\n"); 3154b8a3d8a89814dc3fb365f18d01733e26eb495a1Mathias Agopian return 0; 31673e90268adf4c9638b8d820a802e5e9a8ebe6597Pixelflinger} 3172a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian 3182a967b3fff07b8711aef41f839ad7521323bb64dMathias Agopian// http://www.csee.umbc.edu/help/sound/AFsp-V2R1/html/audio/ResampAudio.html 319