From 86eae0e5931103e040ac2cdd023ef5db252e09f6 Mon Sep 17 00:00:00 2001 From: Andy Hung Date: Mon, 9 Dec 2013 12:12:46 -0800 Subject: Audio resampler update to add S16 filters This does not affect the existing resamplers. New resampler accessed through additional quality settings: DYN_LOW_QUALITY = 5 DYN_MED_QUALITY = 6 DYN_HIGH_QUALITY = 7 Change-Id: Iebbd31871e808a4a6dee3f3abfd7e9dcf77c48e1 Signed-off-by: Andy Hung --- services/audioflinger/AudioResamplerFirGen.h | 307 +++++++++++++++++++++++++++ 1 file changed, 307 insertions(+) create mode 100644 services/audioflinger/AudioResamplerFirGen.h (limited to 'services/audioflinger/AudioResamplerFirGen.h') diff --git a/services/audioflinger/AudioResamplerFirGen.h b/services/audioflinger/AudioResamplerFirGen.h new file mode 100644 index 0000000..fa6ee3e --- /dev/null +++ b/services/audioflinger/AudioResamplerFirGen.h @@ -0,0 +1,307 @@ +/* + * Copyright (C) 2013 The Android Open Source Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef ANDROID_AUDIO_RESAMPLER_FIR_GEN_H +#define ANDROID_AUDIO_RESAMPLER_FIR_GEN_H + +namespace android { + +/* + * Sinc function is the traditional variant. + * + * TODO: Investigate optimizations (regular sampling grid, NEON vector accelerations) + * TODO: Remove comparison at 0 and trap at a higher level. + * + */ + +static inline double sinc(double x) { + if (fabs(x) < FLT_MIN) { + return 1.; + } + return sin(x) / x; +} + +static inline double sqr(double x) { + return x * x; +} + +/* + * rounds a double to the nearest integer for FIR coefficients. + * + * One variant uses noise shaping, which must keep error history + * to work (the err parameter, initialized to 0). + * The other variant is a non-noise shaped version for + * S32 coefficients (noise shaping doesn't gain much). + * + * Caution: No bounds saturation is applied, but isn't needed in + * this case. + * + * @param x is the value to round. + * + * @param maxval is the maximum integer scale factor expressed as an int64 (for headroom). + * Typically this may be the maximum positive integer+1 (using the fact that double precision + * FIR coefficients generated here are never that close to 1.0 to pose an overflow condition). + * + * @param err is the previous error (actual - rounded) for the previous rounding op. + * + */ + +static inline int64_t toint(double x, int64_t maxval, double& err) { + double val = x * maxval; + double ival = floor(val + 0.5 + err*0.17); + err = val - ival; + return static_cast(ival); +} + +static inline int64_t toint(double x, int64_t maxval) { + return static_cast(floor(x * maxval + 0.5)); +} + +/* + * Modified Bessel function of the first kind + * http://en.wikipedia.org/wiki/Bessel_function + * + * The formulas are taken from Abramowitz and Stegun: + * + * http://people.math.sfu.ca/~cbm/aands/page_375.htm + * http://people.math.sfu.ca/~cbm/aands/page_378.htm + * + * http://dlmf.nist.gov/10.25 + * http://dlmf.nist.gov/10.40 + * + * Note we assume x is nonnegative (the function is symmetric, + * pass in the absolute value as needed). + * + * Constants are compile time derived with templates I0Term<> and + * I0ATerm<> to the precision of the compiler. The series can be expanded + * to any precision needed, but currently set around 24b precision. + * + * We use a bit of template math here, constexpr would probably be + * more appropriate for a C++11 compiler. + * + */ + +template +struct I0Term { + static const double value = I0Term::value/ (4. * N * N); +}; + +template <> +struct I0Term<0> { + static const double value = 1.; +}; + +template +struct I0ATerm { + static const double value = I0ATerm::value * (2.*N-1.) * (2.*N-1.) / (8. * N); +}; + +template <> +struct I0ATerm<0> { // 1/sqrt(2*PI); + static const double value = 0.398942280401432677939946059934381868475858631164934657665925; +}; + +static inline double I0(double x) { + if (x < 3.75) { // TODO: Estrin's method instead of Horner's method? + x *= x; + return I0Term<0>::value + x*( + I0Term<1>::value + x*( + I0Term<2>::value + x*( + I0Term<3>::value + x*( + I0Term<4>::value + x*( + I0Term<5>::value + x*( + I0Term<6>::value)))))); // e < 1.6e-7 + } + // a bit ugly here - perhaps we expand the top series + // to permit computation to x < 20 (a reasonable range) + double y = 1./x; + return exp(x) * sqrt(y) * ( + // note: reciprocal squareroot may be easier! + // http://en.wikipedia.org/wiki/Fast_inverse_square_root + I0ATerm<0>::value + y*( + I0ATerm<1>::value + y*( + I0ATerm<2>::value + y*( + I0ATerm<3>::value + y*( + I0ATerm<4>::value + y*( + I0ATerm<5>::value + y*( + I0ATerm<6>::value + y*( + I0ATerm<7>::value + y*( + I0ATerm<8>::value))))))))); // (... e) < 1.9e-7 +} + +/* + * calculates the transition bandwidth for a Kaiser filter + * + * Formula 3.2.8, Multirate Systems and Filter Banks, PP Vaidyanathan, pg. 48 + * + * @param halfNumCoef is half the number of coefficients per filter phase. + * @param stopBandAtten is the stop band attenuation desired. + * @return the transition bandwidth in normalized frequency (0 <= f <= 0.5) + */ +static inline double firKaiserTbw(int halfNumCoef, double stopBandAtten) { + return (stopBandAtten - 7.95)/(2.*14.36*halfNumCoef); +} + +/* + * calculates the fir transfer response. + * + * calculates the transfer coefficient H(w) for 0 <= w <= PI. + * Be careful be careful to consider the fact that this is an interpolated filter + * of length L, so normalizing H(w)/L is probably what you expect. + */ +template +static inline double firTransfer(const T* coef, int L, int halfNumCoef, double w) { + double accum = static_cast(coef[0])*0.5; + coef += halfNumCoef; // skip first row. + for (int i=1 ; i<=L ; ++i) { + for (int j=0, ix=i ; j(*coef++); + } + } + return accum*2.; +} + +/* + * returns the minimum and maximum |H(f)| bounds + * + * @param coef is the designed polyphase filter banks + * + * @param L is the number of phases (for interpolation) + * + * @param halfNumCoef should be half the number of coefficients for a single + * polyphase. + * + * @param fstart is the normalized frequency start. + * + * @param fend is the normalized frequency end. + * + * @param steps is the number of steps to take (sampling) between frequency start and end + * + * @param firMin returns the minimum transfer |H(f)| found + * + * @param firMax returns the maximum transfer |H(f)| found + * + * 0 <= f <= 0.5. + * This is used to test passband and stopband performance. + */ +template +static void testFir(const T* coef, int L, int halfNumCoef, + double fstart, double fend, int steps, double &firMin, double &firMax) { + double wstart = fstart*(2.*M_PI); + double wend = fend*(2.*M_PI); + double wstep = (wend - wstart)/steps; + double fmax, fmin; + double trf = firTransfer(coef, L, halfNumCoef, wstart); + if (trf<0) { + trf = -trf; + } + fmin = fmax = trf; + wstart += wstep; + for (int i=1; ifmax) { + fmax = trf; + } + else if (trf50dB. + * + * @param fcr is cutoff frequency/sampling rate (<0.5). At this point, the energy + * should be 6dB less. (fcr is where the amplitude drops by half). Use the + * firKaiserTbw() to calculate the transition bandwidth. fcr is the midpoint + * between the stop band and the pass band (fstop+fpass)/2. + * + * @param atten is the attenuation (generally slightly less than 1). + */ + +template +static inline void firKaiserGen(T* coef, int L, int halfNumCoef, + double stopBandAtten, double fcr, double atten) { + // + // Formula 3.2.5, 3.2.7, Multirate Systems and Filter Banks, PP Vaidyanathan, pg. 48 + // + // See also: http://melodi.ee.washington.edu/courses/ee518/notes/lec17.pdf + // + // Kaiser window and beta parameter + // + // | 0.1102*(A - 8.7) A > 50 + // beta = | 0.5842*(A - 21)^0.4 + 0.07886*(A - 21) 21 <= A <= 50 + // | 0. A < 21 + // + // with A is the desired stop-band attenuation in dBFS + // + // 30 dB 2.210 + // 40 dB 3.384 + // 50 dB 4.538 + // 60 dB 5.658 + // 70 dB 6.764 + // 80 dB 7.865 + // 90 dB 8.960 + // 100 dB 10.056 + + const int N = L * halfNumCoef; // non-negative half + const double beta = 0.1102 * (stopBandAtten - 8.7); // >= 50dB always + const double yscale = 2. * atten * fcr / I0(beta); + const double xstep = 2. * M_PI * fcr / L; + const double xfrac = 1. / N; + double err = 0; // for noise shaping on int16_t coefficients + for (int i=0 ; i<=L ; ++i) { // generate an extra set of coefs for interpolation + for (int j=0, ix=i ; j::value) { // int16_t needs noise shaping + *coef++ = static_cast(toint(y, 1ULL<<(sizeof(T)*8-1), err)); + } else { + *coef++ = static_cast(toint(y, 1ULL<<(sizeof(T)*8-1))); + } + } + } +} + +}; // namespace android + +#endif /*ANDROID_AUDIO_RESAMPLER_FIR_GEN_H*/ -- cgit v1.1