From fa55fc84f5942c4e6b4dcac523c7f28d816ad22e Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Thu, 14 Jul 2016 12:04:05 -0700 Subject: [PATCH] Optimized compute of distance filters using log-quantized lookup tables. Magnitude error < 0.25dB for entire parameter space. --- libraries/audio/src/AudioHRTF.cpp | 253 +++++++++++++++++++----------- 1 file changed, 159 insertions(+), 94 deletions(-) diff --git a/libraries/audio/src/AudioHRTF.cpp b/libraries/audio/src/AudioHRTF.cpp index 658af01408..378e2154c1 100644 --- a/libraries/audio/src/AudioHRTF.cpp +++ b/libraries/audio/src/AudioHRTF.cpp @@ -16,6 +16,13 @@ #include "AudioHRTF.h" #include "AudioHRTFData.h" +#ifndef MAX +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#endif +#ifndef MIN +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#endif + // // Equal-gain crossfade // @@ -58,6 +65,103 @@ static const float crossfadeTable[HRTF_BLOCK] = { 0.0024846123f, 0.0019026510f, 0.0013981014f, 0.0009710421f, 0.0006215394f, 0.0003496476f, 0.0001554090f, 0.0000388538f, }; +// +// Model the frequency-dependent attenuation of sound propogation in air. +// +// Fit using linear regression to a log-log model of lowpass cutoff frequency vs distance, +// loosely based on data from Handbook of Acoustics. Only the onset of significant +// attenuation is modelled, not the filter slope. +// +// 1m -> -3dB @ 55kHz +// 10m -> -3dB @ 12kHz +// 100m -> -3dB @ 2.5kHz +// 1km -> -3dB @ 0.6kHz +// 10km -> -3dB @ 0.1kHz +// +static const int NLOWPASS = 64; +static const float lowpassTable[NLOWPASS][5] = { // { b0, b1, b2, a1, a2 } + // distance = 1 + { 0.999772371f, 1.399489756f, 0.454495527f, 1.399458985f, 0.454298669f }, + { 0.999631480f, 1.357609808f, 0.425210203f, 1.357549905f, 0.424901586f }, + { 0.999405154f, 1.311503050f, 0.394349994f, 1.311386830f, 0.393871368f }, + { 0.999042876f, 1.260674595f, 0.361869089f, 1.260450057f, 0.361136504f }, + // distance = 2 + { 0.998465222f, 1.204646525f, 0.327757118f, 1.204214978f, 0.326653886f }, + { 0.997548106f, 1.143019308f, 0.292064663f, 1.142195387f, 0.290436690f }, + { 0.996099269f, 1.075569152f, 0.254941286f, 1.074009405f, 0.252600301f }, + { 0.993824292f, 1.002389610f, 0.216688640f, 0.999469185f, 0.213433357f }, + // distance = 4 + { 0.990280170f, 0.924075266f, 0.177827150f, 0.918684864f, 0.173497723f }, + { 0.984818279f, 0.841917936f, 0.139164195f, 0.832151968f, 0.133748443f }, + { 0.976528670f, 0.758036513f, 0.101832398f, 0.740761682f, 0.095635899f }, + { 0.964216485f, 0.675305244f, 0.067243474f, 0.645654855f, 0.061110348f }, + // distance = 8 + { 0.946463038f, 0.596943020f, 0.036899688f, 0.547879974f, 0.032425772f }, + { 0.921823868f, 0.525770189f, 0.012060451f, 0.447952111f, 0.011702396f }, + { 0.890470015f, 0.463334299f, -0.001227816f, 0.347276405f, 0.005300092f }, + { 0.851335343f, 0.407521164f, -0.009353968f, 0.241900234f, 0.007602305f }, + // distance = 16 + { 0.804237360f, 0.358139558f, -0.014293332f, 0.130934213f, 0.017149373f }, + { 0.750073259f, 0.314581568f, -0.016625381f, 0.014505388f, 0.033524057f }, + { 0.690412072f, 0.275936128f, -0.017054561f, -0.106682490f, 0.055976129f }, + { 0.627245545f, 0.241342015f, -0.016246850f, -0.231302564f, 0.083643275f }, + // distance = 32 + { 0.562700627f, 0.210158533f, -0.014740899f, -0.357562697f, 0.115680957f }, + { 0.498787849f, 0.181982455f, -0.012925406f, -0.483461730f, 0.151306628f }, + { 0.437224055f, 0.156585449f, -0.011055180f, -0.607042210f, 0.189796534f }, + { 0.379336998f, 0.133834032f, -0.009281617f, -0.726580065f, 0.230469477f }, + // distance = 64 + { 0.326040627f, 0.113624970f, -0.007683443f, -0.840693542f, 0.272675696f }, + { 0.277861727f, 0.095845793f, -0.006291936f, -0.948380091f, 0.315795676f }, + { 0.234997480f, 0.080357656f, -0.005109519f, -1.049001190f, 0.359246807f }, + { 0.197386484f, 0.066993521f, -0.004122547f, -1.142236313f, 0.402493771f }, + // distance = 128 + { 0.164780457f, 0.055564709f, -0.003309645f, -1.228023442f, 0.445058962f }, + { 0.136808677f, 0.045870650f, -0.002646850f, -1.306498037f, 0.486530514f }, + { 0.113031290f, 0.037708627f, -0.002110591f, -1.377937457f, 0.526566783f }, + { 0.092980475f, 0.030881892f, -0.001679255f, -1.442713983f, 0.564897095f }, + // distance = 256 + { 0.076190239f, 0.025205585f, -0.001333863f, -1.501257246f, 0.601319206f }, + { 0.062216509f, 0.020510496f, -0.001058229f, -1.554025452f, 0.635694228f }, + { 0.050649464f, 0.016644994f, -0.000838826f, -1.601484205f, 0.667939837f }, + { 0.041120009f, 0.013475547f, -0.000664513f, -1.644091518f, 0.698022561f }, + // distance = 512 + { 0.033302044f, 0.010886252f, -0.000526217f, -1.682287704f, 0.725949783f }, + { 0.026911868f, 0.008777712f, -0.000416605f, -1.716488979f, 0.751761953f }, + { 0.021705773f, 0.007065551f, -0.000329788f, -1.747083800f, 0.775525335f }, + { 0.017476603f, 0.005678758f, -0.000261057f, -1.774431204f, 0.797325509f }, + // distance = 1024 + { 0.014049828f, 0.004558012f, -0.000206658f, -1.798860530f, 0.817261711f }, + { 0.011279504f, 0.003654067f, -0.000163610f, -1.820672082f, 0.835442043f }, + { 0.009044384f, 0.002926264f, -0.000129544f, -1.840138412f, 0.851979516f }, + { 0.007244289f, 0.002341194f, -0.000102586f, -1.857505967f, 0.866988864f }, + // distance = 2048 + { 0.005796846f, 0.001871515f, -0.000081250f, -1.872996926f, 0.880584038f }, + { 0.004634607f, 0.001494933f, -0.000064362f, -1.886811124f, 0.892876302f }, + { 0.003702543f, 0.001193324f, -0.000050993f, -1.899127955f, 0.903972829f }, + { 0.002955900f, 0.000951996f, -0.000040407f, -1.910108223f, 0.913975712f }, + // distance = 4096 + { 0.002358382f, 0.000759068f, -0.000032024f, -1.919895894f, 0.922981321f }, + { 0.001880626f, 0.000604950f, -0.000025383f, -1.928619738f, 0.931079931f }, + { 0.001498926f, 0.000481920f, -0.000020123f, -1.936394836f, 0.938355560f }, + { 0.001194182f, 0.000383767f, -0.000015954f, -1.943323983f, 0.944885977f }, + // distance = 8192 + { 0.000951028f, 0.000305502f, -0.000012651f, -1.949498943f, 0.950742822f }, + { 0.000757125f, 0.000243126f, -0.000010033f, -1.955001608f, 0.955991826f }, + { 0.000602572f, 0.000193434f, -0.000007957f, -1.959905036f, 0.960693085f }, + { 0.000479438f, 0.000153861f, -0.000006312f, -1.964274383f, 0.964901371f }, + // distance = 16384 + { 0.000381374f, 0.000122359f, -0.000005007f, -1.968167752f, 0.968666478f }, + { 0.000303302f, 0.000097288f, -0.000003972f, -1.971636944f, 0.972033562f }, + { 0.000241166f, 0.000077342f, -0.000003151f, -1.974728138f, 0.975043493f }, + { 0.000191726f, 0.000061475f, -0.000002500f, -1.977482493f, 0.977733194f }, + // distance = 32768 + { 0.000152399f, 0.000048857f, -0.000001984f, -1.979936697f, 0.980135969f }, + { 0.000121122f, 0.000038825f, -0.000001574f, -1.982123446f, 0.982281818f }, + { 0.000096252f, 0.000030849f, -0.000001249f, -1.984071877f, 0.984197728f }, + { 0.000076480f, 0.000024509f, -0.000000991f, -1.985807957f, 0.985907955f }, +}; + static const float TWOPI = 6.283185307f; // @@ -578,80 +682,58 @@ static void ThiranBiquad(float f, float& b0, float& b1, float& b2, float& a1, fl b2 = 1.0f; } -// returns the gain of analog (s-plane) lowpass evaluated at w -static double analogFilter(double w0, double w) { - double w0sq, wsq; - double num, den; +// split x into exponent and fraction (0.0f to 1.0f) +static void splitf(float x, int& expn, float& frac) { - w0sq = w0 * w0; - wsq = w * w; + union { float f; int i; } mant, bits = { x }; + const int IEEE754_MANT_BITS = 23; + const int IEEE754_EXPN_BIAS = 127; - num = w0sq * w0sq; - den = wsq * wsq + w0sq * w0sq; - - return sqrt(num / den); + mant.i = bits.i & ((1 << IEEE754_MANT_BITS) - 1); + mant.i |= (IEEE754_EXPN_BIAS << IEEE754_MANT_BITS); + + frac = mant.f - 1.0f; + expn = (bits.i >> IEEE754_MANT_BITS) - IEEE754_EXPN_BIAS; } -// design a lowpass biquad using analog matching -static void LowpassBiquad(double coef[5], double w0) { - double G1; - double wpi, wn, wd; - double wna, wda; - double gn, gd, gnsq, gdsq; - double num, den; - double Wnsq, Wdsq, B, A; - double b0, b1, b2, a0, a1, a2; - double temp, scale; - const double PI = 3.14159265358979323846; +static void distanceBiquad(float distance, float& b0, float& b1, float& b2, float& a1, float& a2) { - // compute the Nyquist gain - wpi = w0 + 2.8 * (1.0 - w0/PI); // minimax-like error - wpi = (wpi > PI) ? PI : wpi; - G1 = analogFilter(w0, wpi); + // + // Computed from a lookup table quantized to distance = 2^(N/4) + // and reconstructed by piecewise linear interpolation. + // Approximation error < 0.25dB + // - // approximate wn and wd - wd = 0.5 * w0; - wn = wd * sqrt(1.0/G1); // down G1 at pi, instead of zeros + float x = distance; + x = MIN(MAX(x, 1.0f), 1<<30); + x *= x; + x *= x; // x = distance^4 - Wnsq = wn * wn; - Wdsq = wd * wd; + // split x into e and frac, such that x = 2^(e+0) + frac * (2^(e+1) - 2^(e+0)) + int e; + float frac; + splitf(x, e, frac); - // analog freqs of wn and wd - wna = 2.0 * atan(wn); - wda = 2.0 * atan(wd); + // clamp to table limits + if (e < 0) { + e = 0; + frac = 0.0f; + } + if (e > NLOWPASS-2) { + e = NLOWPASS-2; + frac = 1.0f; + } + assert(frac >= 0.0f); + assert(frac <= 1.0f); + assert(e+0 >= 0); + assert(e+1 < NLOWPASS); - // normalized analog gains at wna and wda - temp = 1.0 / G1; - gn = temp * analogFilter(w0, wna); - gd = temp * analogFilter(w0, wda); - gnsq = gn * gn; - gdsq = gd * gd; - - // compute B, matching gains at wn and wd - temp = 1.0 / (wn * wd); - den = fabs(gnsq - gdsq); - num = gnsq * (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gdsq * Wdsq); - B = temp * sqrt(num / den); - - // compute A, matching gains at wn and wd - num = (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gnsq * Wdsq); - A = temp * sqrt(num / den); - - // design digital filter via bilinear transform - b0 = G1 * (1.0 + B + Wnsq); - b1 = G1 * 2.0 * (Wnsq - 1.0); - b2 = G1 * (1.0 - B + Wnsq); - a0 = 1.0 + A + Wdsq; - a1 = 2.0 * (Wdsq - 1.0); - a2 = 1.0 - A + Wdsq; - - // normalize - scale = 1.0 / a0; - coef[0] = b0 * scale; - coef[1] = b1 * scale; - coef[2] = b2 * scale; - coef[3] = a1 * scale; - coef[4] = a2 * scale; + // piecewise linear interpolation + b0 = lowpassTable[e+0][0] + frac * (lowpassTable[e+1][0] - lowpassTable[e+0][0]); + b1 = lowpassTable[e+0][1] + frac * (lowpassTable[e+1][1] - lowpassTable[e+0][1]); + b2 = lowpassTable[e+0][2] + frac * (lowpassTable[e+1][2] - lowpassTable[e+0][2]); + a1 = lowpassTable[e+0][3] + frac * (lowpassTable[e+1][3] - lowpassTable[e+0][3]); + a2 = lowpassTable[e+0][4] + frac * (lowpassTable[e+1][4] - lowpassTable[e+0][4]); } // compute new filters for a given azimuth, distance and gain @@ -739,38 +821,21 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela } // - // Model the frequency-dependent attenuation of sound propogation in air. - // Fit using linear regression to a log-log model of lowpass cutoff frequency vs distance, - // loosely based on data from Handbook of Acoustics. Only the onset of significant - // attenuation is modelled, not the filter slope. + // Second biquad implements the distance filter. // - // 1m -> -3dB @ 55kHz - // 10m -> -3dB @ 12kHz - // 100m -> -3dB @ 2.5kHz - // 1km -> -3dB @ 0.6kHz - // 10km -> -3dB @ 0.1kHz - // - distance = (distance < 1.0f) ? 1.0f : distance; - double freq = exp2(-0.666 * log2(distance) + 15.75); - double coef[5]; - LowpassBiquad(coef, (double)TWOPI * freq / 24000); + distanceBiquad(distance, b0, b1, b2, a1, a2); - // TESTING: compute attn at w=pi - //double num = coef[0] - coef[1] + coef[2]; - //double den = 1.0 - coef[3] + coef[4]; - //double mag = 10 * log10((num * num) / (den * den)); + bqCoef[0][channel+4] = b0; + bqCoef[1][channel+4] = b1; + bqCoef[2][channel+4] = b2; + bqCoef[3][channel+4] = a1; + bqCoef[4][channel+4] = a2; - bqCoef[0][channel+4] = (float)coef[0]; - bqCoef[1][channel+4] = (float)coef[1]; - bqCoef[2][channel+4] = (float)coef[2]; - bqCoef[3][channel+4] = (float)coef[3]; - bqCoef[4][channel+4] = (float)coef[4]; - - bqCoef[0][channel+5] = (float)coef[0]; - bqCoef[1][channel+5] = (float)coef[1]; - bqCoef[2][channel+5] = (float)coef[2]; - bqCoef[3][channel+5] = (float)coef[3]; - bqCoef[4][channel+5] = (float)coef[4]; + bqCoef[0][channel+5] = b0; + bqCoef[1][channel+5] = b1; + bqCoef[2][channel+5] = b2; + bqCoef[3][channel+5] = a1; + bqCoef[4][channel+5] = a2; } void AudioHRTF::render(int16_t* input, float* output, int index, float azimuth, float distance, float gain, int numFrames) {