diff --git a/assignment-client/src/audio/AudioMixerSlave.cpp b/assignment-client/src/audio/AudioMixerSlave.cpp index b975ba8bfe..ccdd2fca4f 100644 --- a/assignment-client/src/audio/AudioMixerSlave.cpp +++ b/assignment-client/src/audio/AudioMixerSlave.cpp @@ -49,7 +49,7 @@ void sendEnvironmentPacket(const SharedNodePointer& node, AudioMixerClientData& inline float approximateGain(const AvatarAudioStream& listeningNodeStream, const PositionalAudioStream& streamToAdd, const glm::vec3& relativePosition); inline float computeGain(const AudioMixerClientData& listenerNodeData, const AvatarAudioStream& listeningNodeStream, - const PositionalAudioStream& streamToAdd, const glm::vec3& relativePosition, bool isEcho); + const PositionalAudioStream& streamToAdd, const glm::vec3& relativePosition, float distance, bool isEcho); inline float computeAzimuth(const AvatarAudioStream& listeningNodeStream, const PositionalAudioStream& streamToAdd, const glm::vec3& relativePosition); @@ -276,7 +276,7 @@ void AudioMixerSlave::addStream(AudioMixerClientData& listenerNodeData, const QU glm::vec3 relativePosition = streamToAdd.getPosition() - listeningNodeStream.getPosition(); float distance = glm::max(glm::length(relativePosition), EPSILON); - float gain = computeGain(listenerNodeData, listeningNodeStream, streamToAdd, relativePosition, isEcho); + float gain = computeGain(listenerNodeData, listeningNodeStream, streamToAdd, relativePosition, distance, isEcho); float azimuth = isEcho ? 0.0f : computeAzimuth(listeningNodeStream, listeningNodeStream, relativePosition); const int HRTF_DATASET_INDEX = 1; @@ -489,9 +489,6 @@ float approximateGain(const AvatarAudioStream& listeningNodeStream, const Positi // avatar: skip attenuation - it is too costly to approximate // distance attenuation: approximate, ignore zone-specific attenuations - // this is a good approximation for streams further than ATTENUATION_START_DISTANCE - // those streams closer will be amplified; amplifying close streams is acceptable - // when throttling, as close streams are expected to be heard by a user float distance = glm::length(relativePosition); return gain / distance; @@ -499,7 +496,7 @@ float approximateGain(const AvatarAudioStream& listeningNodeStream, const Positi } float computeGain(const AudioMixerClientData& listenerNodeData, const AvatarAudioStream& listeningNodeStream, - const PositionalAudioStream& streamToAdd, const glm::vec3& relativePosition, bool isEcho) { + const PositionalAudioStream& streamToAdd, const glm::vec3& relativePosition, float distance, bool isEcho) { float gain = 1.0f; // injector: apply attenuation @@ -536,23 +533,13 @@ float computeGain(const AudioMixerClientData& listenerNodeData, const AvatarAudi break; } } + // translate the zone setting to gain per log2(distance) + float g = glm::clamp(1.0f - attenuationPerDoublingInDistance, EPSILON, 1.0f); - // distance attenuation - const float ATTENUATION_START_DISTANCE = 1.0f; - float distance = glm::length(relativePosition); - assert(ATTENUATION_START_DISTANCE > EPSILON); - if (distance >= ATTENUATION_START_DISTANCE) { - - // translate the zone setting to gain per log2(distance) - float g = 1.0f - attenuationPerDoublingInDistance; - g = glm::clamp(g, EPSILON, 1.0f); - - // calculate the distance coefficient using the distance to this node - float distanceCoefficient = fastExp2f(fastLog2f(g) * fastLog2f(distance/ATTENUATION_START_DISTANCE)); - - // multiply the current attenuation coefficient by the distance coefficient - gain *= distanceCoefficient; - } + // calculate the attenuation using the distance to this node + // reference attenuation of 0dB at distance = 1.0m + gain *= fastExp2f(fastLog2f(g) * fastLog2f(std::max(distance, HRTF_NEARFIELD_MIN))); + gain = std::min(gain, 1.0f / HRTF_NEARFIELD_MIN); return gain; } diff --git a/libraries/audio-client/src/AudioClient.cpp b/libraries/audio-client/src/AudioClient.cpp index 579910d9f7..8379151d1c 100644 --- a/libraries/audio-client/src/AudioClient.cpp +++ b/libraries/audio-client/src/AudioClient.cpp @@ -1824,16 +1824,9 @@ float AudioClient::azimuthForSource(const glm::vec3& relativePosition) { float AudioClient::gainForSource(float distance, float volume) { - const float ATTENUATION_BEGINS_AT_DISTANCE = 1.0f; - - // I'm assuming that the AudioMixer's getting of the stream's attenuation - // factor is basically same as getting volume - float gain = volume; - - // attenuate based on distance - if (distance >= ATTENUATION_BEGINS_AT_DISTANCE) { - gain /= (distance/ATTENUATION_BEGINS_AT_DISTANCE); // attenuation = -6dB * log2(distance) - } + // attenuation = -6dB * log2(distance) + // reference attenuation of 0dB at distance = 1.0m + float gain = volume / std::max(distance, HRTF_NEARFIELD_MIN); return gain; } diff --git a/libraries/audio/src/AudioHRTF.cpp b/libraries/audio/src/AudioHRTF.cpp index f2f0235ccb..aa951210a6 100644 --- a/libraries/audio/src/AudioHRTF.cpp +++ b/libraries/audio/src/AudioHRTF.cpp @@ -71,6 +71,47 @@ ALIGN32 static const float crossfadeTable[HRTF_BLOCK] = { 0.0029059408f, 0.0022253666f, 0.0016352853f, 0.0011358041f, 0.0007270137f, 0.0004089886f, 0.0001817865f, 0.0000454487f, }; +// +// Fast approximation of the azimuth parallax correction, +// for azimuth = [-pi,pi] and distance = [0.125,2]. +// +// Correction becomes 0 at distance = 2. +// Correction is continuous and smooth. +// +static const int NAZIMUTH = 8; +static const float azimuthTable[NAZIMUTH][3] = { + { 0.018719007f, 0.097263971f, 0.080748954f }, // [-4pi/4,-3pi/4] + { 0.066995833f, 0.319754290f, 0.336963269f }, // [-3pi/4,-2pi/4] + { -0.066989851f, -0.101178847f, 0.006359474f }, // [-2pi/4,-1pi/4] + { -0.018727343f, -0.020357568f, 0.040065626f }, // [-1pi/4,-0pi/4] + { -0.005519051f, -0.018744412f, 0.040065629f }, // [ 0pi/4, 1pi/4] + { -0.001201296f, -0.025103593f, 0.042396711f }, // [ 1pi/4, 2pi/4] + { 0.001198959f, -0.032642381f, 0.048316220f }, // [ 2pi/4, 3pi/4] + { 0.005519640f, -0.053424870f, 0.073296888f }, // [ 3pi/4, 4pi/4] +}; + +// +// Model the normalized near-field Distance Variation Filter. +// +// This version is parameterized by the DC gain correction, instead of directly by azimuth and distance. +// A first-order shelving filter is used to minimize the disturbance in ITD. +// +// Loosely based on data from S. Spagnol, "Distance rendering and perception of nearby virtual sound sources +// with a near-field filter model,” Applied Acoustics (2017) +// +static const int NNEARFIELD = 9; +static const float nearFieldTable[NNEARFIELD][3] = { // { b0, b1, a1 } + { 0.008410604f, -0.000262748f, -0.991852144f }, // gain = 1/256 + { 0.016758914f, -0.000529590f, -0.983770676f }, // gain = 1/128 + { 0.033270607f, -0.001075350f, -0.967804743f }, // gain = 1/64 + { 0.065567740f, -0.002213762f, -0.936646021f }, // gain = 1/32 + { 0.127361554f, -0.004667324f, -0.877305769f }, // gain = 1/16 + { 0.240536414f, -0.010201827f, -0.769665412f }, // gain = 1/8 + { 0.430858205f, -0.023243052f, -0.592384847f }, // gain = 1/4 + { 0.703238106f, -0.054157913f, -0.350919807f }, // gain = 1/2 + { 1.000000000f, -0.123144711f, -0.123144711f }, // gain = 1/1 +}; + // // Model the frequency-dependent attenuation of sound propogation in air. // @@ -168,6 +209,8 @@ static const float lowpassTable[NLOWPASS][5] = { // { b0, b1, b2, a1, a2 } { 0.000076480f, 0.000024509f, -0.000000991f, -1.985807957f, 0.985907955f }, }; +static const float HALFPI = 1.570796327f; +static const float PI = 3.141592654f; static const float TWOPI = 6.283185307f; // @@ -716,7 +759,7 @@ static void distanceBiquad(float distance, float& b0, float& b1, float& b2, floa // float x = distance; - x = MIN(MAX(x, 1.0f), 1<<30); + x = MIN(x, 1<<30); x *= x; x *= x; // x = distance^4 @@ -747,33 +790,181 @@ static void distanceBiquad(float distance, float& b0, float& b1, float& b2, floa a2 = lowpassTable[e+0][4] + frac * (lowpassTable[e+1][4] - lowpassTable[e+0][4]); } +// +// Geometric correction of the azimuth, to the angle at each ear. +// D. Brungart, "Auditory parallax effects in the HRTF for nearby sources," IEEE WASPAA (1999). +// +static void nearFieldAzimuthCorrection(float azimuth, float distance, float& azimuthL, float& azimuthR) { + +#ifdef HRTF_AZIMUTH_EXACT + float dx = distance * cosf(azimuth); + float dy = distance * sinf(azimuth); + + float dx0 = HRTF_AZIMUTH_REF * cosf(azimuth); + float dy0 = HRTF_AZIMUTH_REF * sinf(azimuth); + + azimuthL += atan2f(dy + HRTF_HEAD_RADIUS, dx) - atan2f(dy0 + HRTF_HEAD_RADIUS, dx0); + azimuthR += atan2f(dy - HRTF_HEAD_RADIUS, dx) - atan2f(dy0 - HRTF_HEAD_RADIUS, dx0); +#else + // at reference distance, the azimuth parallax is correct + float fy = (HRTF_AZIMUTH_REF - distance) / distance; + + float x0 = +azimuth; + float x1 = -azimuth; // compute using symmetry + + const float RAD_TO_INDEX = 1.2732394f; // 8/(2*pi), rounded down + int k0 = (int)(RAD_TO_INDEX * x0 + 4.0f); + int k1 = (NAZIMUTH-1) - k0; + assert(k0 >= 0); + assert(k1 >= 0); + assert(k0 < NAZIMUTH); + assert(k1 < NAZIMUTH); + + // piecewise polynomial approximation over azimuth=[-pi,pi] + float fx0 = (azimuthTable[k0][0] * x0 + azimuthTable[k0][1]) * x0 + azimuthTable[k0][2]; + float fx1 = (azimuthTable[k1][0] * x1 + azimuthTable[k1][1]) * x1 + azimuthTable[k1][2]; + + // approximate the azimuth correction + // NOTE: must converge to 0 when distance = HRTF_AZIMUTH_REF + azimuthL += fx0 * fy; + azimuthR -= fx1 * fy; +#endif +} + +// +// Approximate the near-field DC gain correction at each ear. +// +static void nearFieldGainCorrection(float azimuth, float distance, float& gainL, float& gainR) { + + // normalized distance factor = [0,1] as distance = [HRTF_NEARFIELD_MAX,HRTF_HEAD_RADIUS] + assert(distance < HRTF_NEARFIELD_MAX); + assert(distance > HRTF_HEAD_RADIUS); + float d = (HRTF_NEARFIELD_MAX - distance) * ( 1.0f / (HRTF_NEARFIELD_MAX - HRTF_HEAD_RADIUS)); + + // angle of incidence at each ear + float angleL = azimuth + HALFPI; + float angleR = azimuth - HALFPI; + if (angleL > +PI) { + angleL -= TWOPI; + } + if (angleR < -PI) { + angleR += TWOPI; + } + assert(angleL >= -PI); + assert(angleL <= +PI); + assert(angleR >= -PI); + assert(angleR <= +PI); + + // normalized occlusion factor = [0,1] as angle of incidence = [0,pi] + angleL *= angleL; + angleR *= angleR; + float cL = ((-0.000452339132f * angleL - 0.00173192444f) * angleL + 0.162476536f) * angleL; + float cR = ((-0.000452339132f * angleR - 0.00173192444f) * angleR + 0.162476536f) * angleR; + + // approximate the gain correction + // NOTE: must converge to 0 when distance = HRTF_NEARFIELD_MAX + gainL -= d * cL; + gainR -= d * cR; +} + +// +// Approximate the normalized near-field Distance Variation Function at each ear. +// A. Kan, "Distance Variation Function for simulation of near-field virtual auditory space," IEEE ICASSP (2006) +// +static void nearFieldFilter(float gain, float& b0, float& b1, float& a1) { + + // + // Computed from a lookup table quantized to gain = 2^(-N) + // and reconstructed by piecewise linear interpolation. + // + + // split gain into e and frac, such that gain = 2^(e+0) + frac * (2^(e+1) - 2^(e+0)) + int e; + float frac; + splitf(gain, e, frac); + + // clamp to table limits + e += NNEARFIELD-1; + if (e < 0) { + e = 0; + frac = 0.0f; + } + if (e > NNEARFIELD-2) { + e = NNEARFIELD-2; + frac = 1.0f; + } + assert(frac >= 0.0f); + assert(frac <= 1.0f); + assert(e+0 >= 0); + assert(e+1 < NNEARFIELD); + + // piecewise linear interpolation + b0 = nearFieldTable[e+0][0] + frac * (nearFieldTable[e+1][0] - nearFieldTable[e+0][0]); + b1 = nearFieldTable[e+0][1] + frac * (nearFieldTable[e+1][1] - nearFieldTable[e+0][1]); + a1 = nearFieldTable[e+0][2] + frac * (nearFieldTable[e+1][2] - nearFieldTable[e+0][2]); +} + +static void azimuthToIndex(float azimuth, int& index0, int& index1, float& frac) { + + // convert from radians to table units + azimuth *= (HRTF_AZIMUTHS / TWOPI); + + if (azimuth < 0.0f) { + azimuth += HRTF_AZIMUTHS; + } + + // table parameters + index0 = (int)azimuth; + index1 = index0 + 1; + frac = azimuth - (float)index0; + + if (index1 >= HRTF_AZIMUTHS) { + index1 -= HRTF_AZIMUTHS; + } + + assert((index0 >= 0) && (index0 < HRTF_AZIMUTHS)); + assert((index1 >= 0) && (index1 < HRTF_AZIMUTHS)); + assert((frac >= 0.0f) && (frac < 1.0f)); +} + // compute new filters for a given azimuth, distance and gain static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int delay[4], int index, float azimuth, float distance, float gain, int channel) { - // convert from radians to table units - azimuth *= HRTF_AZIMUTHS / TWOPI; - - // wrap to principle value - while (azimuth < 0.0f) { - azimuth += HRTF_AZIMUTHS; + if (azimuth > PI) { + azimuth -= TWOPI; } - while (azimuth >= HRTF_AZIMUTHS) { - azimuth -= HRTF_AZIMUTHS; + assert(azimuth >= -PI); + assert(azimuth <= +PI); + + distance = MAX(distance, HRTF_NEARFIELD_MIN); + + // compute the azimuth correction at each ear + float azimuthL = azimuth; + float azimuthR = azimuth; + if (distance < HRTF_AZIMUTH_REF) { + nearFieldAzimuthCorrection(azimuth, distance, azimuthL, azimuthR); } - // table parameters - int az0 = (int)azimuth; - int az1 = (az0 + 1) % HRTF_AZIMUTHS; - float frac = azimuth - (float)az0; + // compute the DC gain correction at each ear + float gainL = 1.0f; + float gainR = 1.0f; + if (distance < HRTF_NEARFIELD_MAX) { + nearFieldGainCorrection(azimuth, distance, gainL, gainR); + } - assert((az0 >= 0) && (az0 < HRTF_AZIMUTHS)); - assert((az1 >= 0) && (az1 < HRTF_AZIMUTHS)); - assert((frac >= 0.0f) && (frac < 1.0f)); + // parameters for table interpolation + int azL0, azR0, az0; + int azL1, azR1, az1; + float fracL, fracR, frac; + + azimuthToIndex(azimuthL, azL0, azL1, fracL); + azimuthToIndex(azimuthR, azR0, azR1, fracR); + azimuthToIndex(azimuth, az0, az1, frac); // interpolate FIR - interpolate(firCoef[channel+0], ir_table_table[index][az0][0], ir_table_table[index][az1][0], frac, gain); - interpolate(firCoef[channel+1], ir_table_table[index][az0][1], ir_table_table[index][az1][1], frac, gain); + interpolate(firCoef[channel+0], ir_table_table[index][azL0][0], ir_table_table[index][azL1][0], fracL, gain * gainL); + interpolate(firCoef[channel+1], ir_table_table[index][azR0][1], ir_table_table[index][azR1][1], fracR, gain * gainR); // interpolate ITD float itd = (1.0f - frac) * itd_table_table[index][az0] + frac * itd_table_table[index][az1]; @@ -832,21 +1023,42 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela } // - // Second biquad implements the distance filter. + // Second biquad implements the near-field or distance filter. // - distanceBiquad(distance, b0, b1, b2, a1, a2); + if (distance < HRTF_NEARFIELD_MAX) { - 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; + nearFieldFilter(gainL, b0, b1, a1); - 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; + bqCoef[0][channel+4] = b0; + bqCoef[1][channel+4] = b1; + bqCoef[2][channel+4] = 0.0f; + bqCoef[3][channel+4] = a1; + bqCoef[4][channel+4] = 0.0f; + + nearFieldFilter(gainR, b0, b1, a1); + + bqCoef[0][channel+5] = b0; + bqCoef[1][channel+5] = b1; + bqCoef[2][channel+5] = 0.0f; + bqCoef[3][channel+5] = a1; + bqCoef[4][channel+5] = 0.0f; + + } else { + + distanceBiquad(distance, b0, b1, b2, a1, a2); + + 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+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) { diff --git a/libraries/audio/src/AudioHRTF.h b/libraries/audio/src/AudioHRTF.h index a197264994..8993842d6e 100644 --- a/libraries/audio/src/AudioHRTF.h +++ b/libraries/audio/src/AudioHRTF.h @@ -23,6 +23,12 @@ static const int HRTF_BLOCK = 240; // block processing size static const float HRTF_GAIN = 1.0f; // HRTF global gain adjustment +// Near-field HRTF +static const float HRTF_AZIMUTH_REF = 2.0f; // IRCAM Listen HRTF was recorded at 2 meters +static const float HRTF_NEARFIELD_MAX = 1.0f; // distance in meters +static const float HRTF_NEARFIELD_MIN = 0.125f; // distance in meters +static const float HRTF_HEAD_RADIUS = 0.0875f; // average human head in meters + class AudioHRTF { public: