AudioResamplerFirGen.h revision 86eae0e5931103e040ac2cdd023ef5db252e09f6
186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Copyright (C) 2013 The Android Open Source Project 386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Licensed under the Apache License, Version 2.0 (the "License"); 586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * you may not use this file except in compliance with the License. 686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * You may obtain a copy of the License at 786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * http://www.apache.org/licenses/LICENSE-2.0 986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 1086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Unless required by applicable law or agreed to in writing, software 1186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * distributed under the License is distributed on an "AS IS" BASIS, 1286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 1386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * See the License for the specific language governing permissions and 1486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * limitations under the License. 1586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 1686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 1786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung#ifndef ANDROID_AUDIO_RESAMPLER_FIR_GEN_H 1886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung#define ANDROID_AUDIO_RESAMPLER_FIR_GEN_H 1986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 2086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungnamespace android { 2186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 2286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 2386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Sinc function is the traditional variant. 2486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 2586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * TODO: Investigate optimizations (regular sampling grid, NEON vector accelerations) 2686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * TODO: Remove comparison at 0 and trap at a higher level. 2786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 2886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 2986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 3086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double sinc(double x) { 3186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (fabs(x) < FLT_MIN) { 3286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return 1.; 3386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 3486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return sin(x) / x; 3586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 3686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 3786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double sqr(double x) { 3886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return x * x; 3986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 4086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 4186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 4286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * rounds a double to the nearest integer for FIR coefficients. 4386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 4486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * One variant uses noise shaping, which must keep error history 4586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * to work (the err parameter, initialized to 0). 4686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * The other variant is a non-noise shaped version for 4786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * S32 coefficients (noise shaping doesn't gain much). 4886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 4986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Caution: No bounds saturation is applied, but isn't needed in 5086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * this case. 5186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 5286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param x is the value to round. 5386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 5486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param maxval is the maximum integer scale factor expressed as an int64 (for headroom). 5586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Typically this may be the maximum positive integer+1 (using the fact that double precision 5686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * FIR coefficients generated here are never that close to 1.0 to pose an overflow condition). 5786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 5886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param err is the previous error (actual - rounded) for the previous rounding op. 5986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 6086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 6186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 6286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline int64_t toint(double x, int64_t maxval, double& err) { 6386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double val = x * maxval; 6486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double ival = floor(val + 0.5 + err*0.17); 6586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung err = val - ival; 6686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return static_cast<int64_t>(ival); 6786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 6886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 6986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline int64_t toint(double x, int64_t maxval) { 7086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return static_cast<int64_t>(floor(x * maxval + 0.5)); 7186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 7286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 7386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 7486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Modified Bessel function of the first kind 7586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * http://en.wikipedia.org/wiki/Bessel_function 7686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 7786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * The formulas are taken from Abramowitz and Stegun: 7886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 7986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * http://people.math.sfu.ca/~cbm/aands/page_375.htm 8086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * http://people.math.sfu.ca/~cbm/aands/page_378.htm 8186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 8286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * http://dlmf.nist.gov/10.25 8386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * http://dlmf.nist.gov/10.40 8486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 8586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Note we assume x is nonnegative (the function is symmetric, 8686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * pass in the absolute value as needed). 8786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 8886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Constants are compile time derived with templates I0Term<> and 8986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * I0ATerm<> to the precision of the compiler. The series can be expanded 9086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * to any precision needed, but currently set around 24b precision. 9186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 9286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * We use a bit of template math here, constexpr would probably be 9386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * more appropriate for a C++11 compiler. 9486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 9586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 9686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 9786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <int N> 9886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstruct I0Term { 9986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung static const double value = I0Term<N-1>::value/ (4. * N * N); 10086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung}; 10186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 10286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <> 10386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstruct I0Term<0> { 10486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung static const double value = 1.; 10586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung}; 10686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 10786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <int N> 10886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstruct I0ATerm { 10986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung static const double value = I0ATerm<N-1>::value * (2.*N-1.) * (2.*N-1.) / (8. * N); 11086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung}; 11186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 11286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <> 11386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstruct I0ATerm<0> { // 1/sqrt(2*PI); 11486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung static const double value = 0.398942280401432677939946059934381868475858631164934657665925; 11586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung}; 11686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 11786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double I0(double x) { 11886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (x < 3.75) { // TODO: Estrin's method instead of Horner's method? 11986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung x *= x; 12086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return I0Term<0>::value + x*( 12186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0Term<1>::value + x*( 12286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0Term<2>::value + x*( 12386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0Term<3>::value + x*( 12486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0Term<4>::value + x*( 12586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0Term<5>::value + x*( 12686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0Term<6>::value)))))); // e < 1.6e-7 12786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 12886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // a bit ugly here - perhaps we expand the top series 12986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // to permit computation to x < 20 (a reasonable range) 13086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double y = 1./x; 13186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return exp(x) * sqrt(y) * ( 13286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // note: reciprocal squareroot may be easier! 13386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // http://en.wikipedia.org/wiki/Fast_inverse_square_root 13486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<0>::value + y*( 13586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<1>::value + y*( 13686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<2>::value + y*( 13786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<3>::value + y*( 13886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<4>::value + y*( 13986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<5>::value + y*( 14086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<6>::value + y*( 14186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<7>::value + y*( 14286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung I0ATerm<8>::value))))))))); // (... e) < 1.9e-7 14386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 14486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 14586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 14686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * calculates the transition bandwidth for a Kaiser filter 14786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 14886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Formula 3.2.8, Multirate Systems and Filter Banks, PP Vaidyanathan, pg. 48 14986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 15086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param halfNumCoef is half the number of coefficients per filter phase. 15186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param stopBandAtten is the stop band attenuation desired. 15286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @return the transition bandwidth in normalized frequency (0 <= f <= 0.5) 15386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 15486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double firKaiserTbw(int halfNumCoef, double stopBandAtten) { 15586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return (stopBandAtten - 7.95)/(2.*14.36*halfNumCoef); 15686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 15786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 15886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 15986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * calculates the fir transfer response. 16086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 16186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * calculates the transfer coefficient H(w) for 0 <= w <= PI. 16286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Be careful be careful to consider the fact that this is an interpolated filter 16386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * of length L, so normalizing H(w)/L is probably what you expect. 16486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 16586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <typename T> 16686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline double firTransfer(const T* coef, int L, int halfNumCoef, double w) { 16786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double accum = static_cast<double>(coef[0])*0.5; 16886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung coef += halfNumCoef; // skip first row. 16986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung for (int i=1 ; i<=L ; ++i) { 17086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung for (int j=0, ix=i ; j<halfNumCoef ; ++j, ix+=L) { 17186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung accum += cos(ix*w)*static_cast<double>(*coef++); 17286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 17386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 17486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung return accum*2.; 17586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 17686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 17786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 17886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * returns the minimum and maximum |H(f)| bounds 17986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 18086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param coef is the designed polyphase filter banks 18186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 18286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param L is the number of phases (for interpolation) 18386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 18486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param halfNumCoef should be half the number of coefficients for a single 18586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * polyphase. 18686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 18786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param fstart is the normalized frequency start. 18886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 18986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param fend is the normalized frequency end. 19086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 19186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param steps is the number of steps to take (sampling) between frequency start and end 19286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 19386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param firMin returns the minimum transfer |H(f)| found 19486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 19586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param firMax returns the maximum transfer |H(f)| found 19686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 19786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 0 <= f <= 0.5. 19886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * This is used to test passband and stopband performance. 19986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 20086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <typename T> 20186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic void testFir(const T* coef, int L, int halfNumCoef, 20286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double fstart, double fend, int steps, double &firMin, double &firMax) { 20386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double wstart = fstart*(2.*M_PI); 20486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double wend = fend*(2.*M_PI); 20586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double wstep = (wend - wstart)/steps; 20686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double fmax, fmin; 20786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double trf = firTransfer(coef, L, halfNumCoef, wstart); 20886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (trf<0) { 20986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung trf = -trf; 21086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 21186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung fmin = fmax = trf; 21286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung wstart += wstep; 21386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung for (int i=1; i<steps; ++i) { 21486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung trf = firTransfer(coef, L, halfNumCoef, wstart); 21586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (trf<0) { 21686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung trf = -trf; 21786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 21886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (trf>fmax) { 21986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung fmax = trf; 22086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 22186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung else if (trf<fmin) { 22286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung fmin = trf; 22386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 22486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung wstart += wstep; 22586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 22686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // renormalize - this is only needed for integer filter types 22786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double norm = 1./((1ULL<<(sizeof(T)*8-1))*L); 22886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 22986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung firMin = fmin * norm; 23086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung firMax = fmax * norm; 23186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 23286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 23386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung/* 23486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * Calculates the polyphase filter banks based on a windowed sinc function. 23586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 23686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * The windowed sinc is an odd length symmetric filter of exactly L*halfNumCoef*2+1 23786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * taps for the entire kernel. This is then decomposed into L+1 polyphase filterbanks. 23886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * The last filterbank is used for interpolation purposes (and is mostly composed 23986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * of the first bank shifted by one sample), and is unnecessary if one does 24086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * not do interpolation. 24186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 24286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param coef is the caller allocated space for coefficients. This should be 24386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * exactly (L+1)*halfNumCoef in size. 24486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 24586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param L is the number of phases (for interpolation) 24686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 24786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param halfNumCoef should be half the number of coefficients for a single 24886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * polyphase. 24986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 25086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param stopBandAtten is the stopband value, should be >50dB. 25186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 25286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param fcr is cutoff frequency/sampling rate (<0.5). At this point, the energy 25386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * should be 6dB less. (fcr is where the amplitude drops by half). Use the 25486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * firKaiserTbw() to calculate the transition bandwidth. fcr is the midpoint 25586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * between the stop band and the pass band (fstop+fpass)/2. 25686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * 25786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung * @param atten is the attenuation (generally slightly less than 1). 25886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung */ 25986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 26086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungtemplate <typename T> 26186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hungstatic inline void firKaiserGen(T* coef, int L, int halfNumCoef, 26286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double stopBandAtten, double fcr, double atten) { 26386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 26486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // Formula 3.2.5, 3.2.7, Multirate Systems and Filter Banks, PP Vaidyanathan, pg. 48 26586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 26686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // See also: http://melodi.ee.washington.edu/courses/ee518/notes/lec17.pdf 26786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 26886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // Kaiser window and beta parameter 26986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 27086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // | 0.1102*(A - 8.7) A > 50 27186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // beta = | 0.5842*(A - 21)^0.4 + 0.07886*(A - 21) 21 <= A <= 50 27286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // | 0. A < 21 27386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 27486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // with A is the desired stop-band attenuation in dBFS 27586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 27686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 30 dB 2.210 27786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 40 dB 3.384 27886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 50 dB 4.538 27986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 60 dB 5.658 28086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 70 dB 6.764 28186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 80 dB 7.865 28286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 90 dB 8.960 28386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // 100 dB 10.056 28486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 28586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung const int N = L * halfNumCoef; // non-negative half 28686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung const double beta = 0.1102 * (stopBandAtten - 8.7); // >= 50dB always 28786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung const double yscale = 2. * atten * fcr / I0(beta); 28886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung const double xstep = 2. * M_PI * fcr / L; 28986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung const double xfrac = 1. / N; 29086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double err = 0; // for noise shaping on int16_t coefficients 29186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung for (int i=0 ; i<=L ; ++i) { // generate an extra set of coefs for interpolation 29286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung for (int j=0, ix=i ; j<halfNumCoef ; ++j, ix+=L) { 29386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung double y = I0(beta * sqrt(1.0 - sqr(ix * xfrac))) * sinc(ix * xstep) * yscale; 29486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 29586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung // (caution!) float version does not need rounding 29686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung if (is_same<T, int16_t>::value) { // int16_t needs noise shaping 29786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung *coef++ = static_cast<T>(toint(y, 1ULL<<(sizeof(T)*8-1), err)); 29886eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } else { 29986eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung *coef++ = static_cast<T>(toint(y, 1ULL<<(sizeof(T)*8-1))); 30086eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 30186eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 30286eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung } 30386eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung} 30486eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 30586eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung}; // namespace android 30686eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung 30786eae0e5931103e040ac2cdd023ef5db252e09f6Andy Hung#endif /*ANDROID_AUDIO_RESAMPLER_FIR_GEN_H*/ 308