From e367167709faaf00324924a145dab6c4c34604ad Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Thu, 22 Feb 2018 15:53:28 -0800 Subject: [PATCH 1/6] New HRTF with near-field audio spatialization --- libraries/audio/src/AudioHRTF.cpp | 221 ++++++++++++++++++++++++++---- libraries/audio/src/AudioHRTF.h | 4 + 2 files changed, 195 insertions(+), 30 deletions(-) diff --git a/libraries/audio/src/AudioHRTF.cpp b/libraries/audio/src/AudioHRTF.cpp index f2f0235ccb..8ccb79b2a5 100644 --- a/libraries/audio/src/AudioHRTF.cpp +++ b/libraries/audio/src/AudioHRTF.cpp @@ -71,6 +71,28 @@ ALIGN32 static const float crossfadeTable[HRTF_BLOCK] = { 0.0029059408f, 0.0022253666f, 0.0016352853f, 0.0011358041f, 0.0007270137f, 0.0004089886f, 0.0001817865f, 0.0000454487f, }; +// +// 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 +190,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 +740,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 +771,151 @@ 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 nearFieldAzimuth(float azimuth, float distance, float& azimuthL, float& azimuthR) { + + // FIXME: optimize + float dx = distance * cosf(azimuth); + float dy = distance * sinf(azimuth); + + azimuthL = atan2f(dy + HRTF_HEAD_RADIUS, dx); + azimuthR = atan2f(dy - HRTF_HEAD_RADIUS, dx); +} + +// +// Approximate the near-field DC gain correction at each ear. +// +static void nearFieldGain(float azimuth, float distance, float& gainL, float& gainR) { + + // normalized distance factor = [0,1] as distance = [1,HRTF_HEAD_RADIUS] + assert(distance < 1.0f); + assert(distance > HRTF_HEAD_RADIUS); + float d = (1.0f - distance) * ( 1.0f / (1.0f - 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: this must converge to 1.0 when distance = 1.0m at all azimuth + gainL = 1.0f - d * cL; + gainR = 1.0f - 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_DISTANCE_MIN); + + // compute the azimuth correction at each ear + float azimuthL, azimuthR; + nearFieldAzimuth(azimuth, distance, azimuthL, azimuthR); + + // compute the DC gain correction at each ear + float gainL = 1.0f; + float gainR = 1.0f; + if (distance < 1.0f) { + nearFieldGain(azimuth, distance, gainL, gainR); } - // table parameters - int az0 = (int)azimuth; - int az1 = (az0 + 1) % HRTF_AZIMUTHS; - float frac = azimuth - (float)az0; + // parameters for table interpolation + int azL0, azR0, az0; + int azL1, azR1, az1; + float fracL, fracR, frac; - assert((az0 >= 0) && (az0 < HRTF_AZIMUTHS)); - assert((az1 >= 0) && (az1 < HRTF_AZIMUTHS)); - assert((frac >= 0.0f) && (frac < 1.0f)); + 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 +974,40 @@ 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 < 1.0f) { + nearFieldFilter(gainL, b0, b1, a1); - 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] = b0; + bqCoef[1][channel+4] = b1; + bqCoef[2][channel+4] = 0.0f; + bqCoef[3][channel+4] = a1; + bqCoef[4][channel+4] = 0.0f; - 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; + 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..b8422b7c73 100644 --- a/libraries/audio/src/AudioHRTF.h +++ b/libraries/audio/src/AudioHRTF.h @@ -23,6 +23,10 @@ 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_DISTANCE_MIN = 0.125f; +static const float HRTF_HEAD_RADIUS = 0.0875f; // average human head + class AudioHRTF { public: From f25b8e8df01c6256e2618dcc073eed89639b4650 Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Thu, 22 Feb 2018 16:05:01 -0800 Subject: [PATCH 2/6] Refactor distance attenuation for near-field HRTF. Adds near-field attenuation from 0.125m to 1m, calibrated to 0dB at 1m. Overload will trigger dynamic range compression at the peak-limiter. --- .../src/audio/AudioMixerSlave.cpp | 23 +++++-------------- libraries/audio-client/src/AudioClient.cpp | 13 +++-------- 2 files changed, 9 insertions(+), 27 deletions(-) diff --git a/assignment-client/src/audio/AudioMixerSlave.cpp b/assignment-client/src/audio/AudioMixerSlave.cpp index abc2f7220d..e5ad9a681a 100644 --- a/assignment-client/src/audio/AudioMixerSlave.cpp +++ b/assignment-client/src/audio/AudioMixerSlave.cpp @@ -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; @@ -536,23 +533,15 @@ 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 *= exp2f(log2f(g) * log2f(std::max(distance, HRTF_DISTANCE_MIN))); + gain = std::min(gain, 1.0f / HRTF_DISTANCE_MIN); return gain; } diff --git a/libraries/audio-client/src/AudioClient.cpp b/libraries/audio-client/src/AudioClient.cpp index 579910d9f7..40a55c272c 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_DISTANCE_MIN); return gain; } From bb2f3cac2cc7140d4e971e414a930a5e1a56fabe Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Tue, 27 Feb 2018 08:51:55 -0800 Subject: [PATCH 3/6] Use named constant for the onset of near-field corrections (1 meter) --- assignment-client/src/audio/AudioMixerSlave.cpp | 4 ++-- libraries/audio-client/src/AudioClient.cpp | 2 +- libraries/audio/src/AudioHRTF.cpp | 12 ++++++------ libraries/audio/src/AudioHRTF.h | 5 +++-- 4 files changed, 12 insertions(+), 11 deletions(-) diff --git a/assignment-client/src/audio/AudioMixerSlave.cpp b/assignment-client/src/audio/AudioMixerSlave.cpp index e5ad9a681a..04b818adee 100644 --- a/assignment-client/src/audio/AudioMixerSlave.cpp +++ b/assignment-client/src/audio/AudioMixerSlave.cpp @@ -540,8 +540,8 @@ float computeGain(const AudioMixerClientData& listenerNodeData, const AvatarAudi // calculate the attenuation using the distance to this node // reference attenuation of 0dB at distance = 1.0m - gain *= exp2f(log2f(g) * log2f(std::max(distance, HRTF_DISTANCE_MIN))); - gain = std::min(gain, 1.0f / HRTF_DISTANCE_MIN); + gain *= exp2f(log2f(g) * log2f(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 40a55c272c..8379151d1c 100644 --- a/libraries/audio-client/src/AudioClient.cpp +++ b/libraries/audio-client/src/AudioClient.cpp @@ -1826,7 +1826,7 @@ float AudioClient::gainForSource(float distance, float volume) { // attenuation = -6dB * log2(distance) // reference attenuation of 0dB at distance = 1.0m - float gain = volume / std::max(distance, HRTF_DISTANCE_MIN); + 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 8ccb79b2a5..a6451b1771 100644 --- a/libraries/audio/src/AudioHRTF.cpp +++ b/libraries/audio/src/AudioHRTF.cpp @@ -790,10 +790,10 @@ static void nearFieldAzimuth(float azimuth, float distance, float& azimuthL, flo // static void nearFieldGain(float azimuth, float distance, float& gainL, float& gainR) { - // normalized distance factor = [0,1] as distance = [1,HRTF_HEAD_RADIUS] - assert(distance < 1.0f); + // 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 = (1.0f - distance) * ( 1.0f / (1.0f - 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; @@ -816,7 +816,7 @@ static void nearFieldGain(float azimuth, float distance, float& gainL, float& ga float cR = ((-0.000452339132f * angleR - 0.00173192444f) * angleR + 0.162476536f) * angleR; // approximate the gain correction - // NOTE: this must converge to 1.0 when distance = 1.0m at all azimuth + // NOTE: this must converge to 1.0 when distance = HRTF_NEARFIELD_MAX at all azimuth gainL = 1.0f - d * cL; gainR = 1.0f - d * cR; } @@ -891,7 +891,7 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela assert(azimuth >= -PI); assert(azimuth <= +PI); - distance = MAX(distance, HRTF_DISTANCE_MIN); + distance = MAX(distance, HRTF_NEARFIELD_MIN); // compute the azimuth correction at each ear float azimuthL, azimuthR; @@ -900,7 +900,7 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela // compute the DC gain correction at each ear float gainL = 1.0f; float gainR = 1.0f; - if (distance < 1.0f) { + if (distance < HRTF_NEARFIELD_MAX) { nearFieldGain(azimuth, distance, gainL, gainR); } diff --git a/libraries/audio/src/AudioHRTF.h b/libraries/audio/src/AudioHRTF.h index b8422b7c73..f0f4dd849b 100644 --- a/libraries/audio/src/AudioHRTF.h +++ b/libraries/audio/src/AudioHRTF.h @@ -24,8 +24,9 @@ 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_DISTANCE_MIN = 0.125f; -static const float HRTF_HEAD_RADIUS = 0.0875f; // average human head +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 { From cc7b9a254b63cc50579ab8541cfbeb48a65c4c63 Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Tue, 27 Feb 2018 09:02:03 -0800 Subject: [PATCH 4/6] Optimized compute of distance attenuation --- assignment-client/src/audio/AudioMixerSlave.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/assignment-client/src/audio/AudioMixerSlave.cpp b/assignment-client/src/audio/AudioMixerSlave.cpp index 04b818adee..4ec4a0f67f 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; @@ -496,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,11 +536,9 @@ float computeGain(const AudioMixerClientData& listenerNodeData, const AvatarAudi // translate the zone setting to gain per log2(distance) float g = glm::clamp(1.0f - attenuationPerDoublingInDistance, EPSILON, 1.0f); - float distance = glm::length(relativePosition); - // calculate the attenuation using the distance to this node // reference attenuation of 0dB at distance = 1.0m - gain *= exp2f(log2f(g) * log2f(std::max(distance, HRTF_NEARFIELD_MIN))); + gain *= fastExp2f(fastLog2f(g) * fastLog2f(std::max(distance, HRTF_NEARFIELD_MIN))); gain = std::min(gain, 1.0f / HRTF_NEARFIELD_MIN); return gain; From 41ccf6afd8640d6f4547d4f72ce46804caf55267 Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Tue, 27 Feb 2018 12:49:32 -0800 Subject: [PATCH 5/6] Compensate the near-field azimuth correction for the parallax baked into the HRTF dataset --- libraries/audio/src/AudioHRTF.cpp | 37 +++++++++++++++++++++---------- libraries/audio/src/AudioHRTF.h | 1 + 2 files changed, 26 insertions(+), 12 deletions(-) diff --git a/libraries/audio/src/AudioHRTF.cpp b/libraries/audio/src/AudioHRTF.cpp index a6451b1771..f0679a3684 100644 --- a/libraries/audio/src/AudioHRTF.cpp +++ b/libraries/audio/src/AudioHRTF.cpp @@ -775,20 +775,28 @@ static void distanceBiquad(float distance, float& b0, float& b1, float& b2, floa // 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 nearFieldAzimuth(float azimuth, float distance, float& azimuthL, float& azimuthR) { +static void nearFieldAzimuthCorrection(float azimuth, float distance, float& azimuthL, float& azimuthR) { - // FIXME: optimize float dx = distance * cosf(azimuth); float dy = distance * sinf(azimuth); - azimuthL = atan2f(dy + HRTF_HEAD_RADIUS, dx); - azimuthR = atan2f(dy - HRTF_HEAD_RADIUS, dx); + // at reference distance, the azimuth parallax is correct + float dx0 = HRTF_AZIMUTH_REF * cosf(azimuth); + float dy0 = HRTF_AZIMUTH_REF * sinf(azimuth); + + float deltaL = atan2f(dy + HRTF_HEAD_RADIUS, dx) - atan2f(dy0 + HRTF_HEAD_RADIUS, dx0); + float deltaR = atan2f(dy - HRTF_HEAD_RADIUS, dx) - atan2f(dy0 - HRTF_HEAD_RADIUS, dx0); + + // add the azimuth correction + // NOTE: must converge to 0 when distance = HRTF_AZIMUTH_REF + azimuthL += deltaL; + azimuthR += deltaR; } // // Approximate the near-field DC gain correction at each ear. // -static void nearFieldGain(float azimuth, float distance, float& gainL, float& gainR) { +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); @@ -816,9 +824,9 @@ static void nearFieldGain(float azimuth, float distance, float& gainL, float& ga float cR = ((-0.000452339132f * angleR - 0.00173192444f) * angleR + 0.162476536f) * angleR; // approximate the gain correction - // NOTE: this must converge to 1.0 when distance = HRTF_NEARFIELD_MAX at all azimuth - gainL = 1.0f - d * cL; - gainR = 1.0f - d * cR; + // NOTE: must converge to 0 when distance = HRTF_NEARFIELD_MAX + gainL -= d * cL; + gainR -= d * cR; } // @@ -894,14 +902,17 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela distance = MAX(distance, HRTF_NEARFIELD_MIN); // compute the azimuth correction at each ear - float azimuthL, azimuthR; - nearFieldAzimuth(azimuth, distance, azimuthL, azimuthR); + float azimuthL = azimuth; + float azimuthR = azimuth; + if (distance < HRTF_AZIMUTH_REF) { + nearFieldAzimuthCorrection(azimuth, distance, azimuthL, azimuthR); + } // compute the DC gain correction at each ear float gainL = 1.0f; float gainR = 1.0f; if (distance < HRTF_NEARFIELD_MAX) { - nearFieldGain(azimuth, distance, gainL, gainR); + nearFieldGainCorrection(azimuth, distance, gainL, gainR); } // parameters for table interpolation @@ -976,7 +987,8 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela // // Second biquad implements the near-field or distance filter. // - if (distance < 1.0f) { + if (distance < HRTF_NEARFIELD_MAX) { + nearFieldFilter(gainL, b0, b1, a1); bqCoef[0][channel+4] = b0; @@ -994,6 +1006,7 @@ static void setFilters(float firCoef[4][HRTF_TAPS], float bqCoef[5][8], int dela bqCoef[4][channel+5] = 0.0f; } else { + distanceBiquad(distance, b0, b1, b2, a1, a2); bqCoef[0][channel+4] = b0; diff --git a/libraries/audio/src/AudioHRTF.h b/libraries/audio/src/AudioHRTF.h index f0f4dd849b..8993842d6e 100644 --- a/libraries/audio/src/AudioHRTF.h +++ b/libraries/audio/src/AudioHRTF.h @@ -24,6 +24,7 @@ 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 From bc7b983f4871c5b4bf2e7e28a57525762eed2c20 Mon Sep 17 00:00:00 2001 From: Ken Cooke Date: Wed, 28 Feb 2018 04:44:20 -0800 Subject: [PATCH 6/6] Fast approximation of the azimuth parallax correction --- libraries/audio/src/AudioHRTF.cpp | 50 +++++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/libraries/audio/src/AudioHRTF.cpp b/libraries/audio/src/AudioHRTF.cpp index f0679a3684..aa951210a6 100644 --- a/libraries/audio/src/AudioHRTF.cpp +++ b/libraries/audio/src/AudioHRTF.cpp @@ -71,6 +71,25 @@ 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. // @@ -777,20 +796,39 @@ static void distanceBiquad(float distance, float& b0, float& b1, float& b2, floa // 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); - // at reference distance, the azimuth parallax is correct float dx0 = HRTF_AZIMUTH_REF * cosf(azimuth); float dy0 = HRTF_AZIMUTH_REF * sinf(azimuth); - float deltaL = atan2f(dy + HRTF_HEAD_RADIUS, dx) - atan2f(dy0 + HRTF_HEAD_RADIUS, dx0); - float deltaR = atan2f(dy - HRTF_HEAD_RADIUS, dx) - atan2f(dy0 - HRTF_HEAD_RADIUS, dx0); + 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; - // add the azimuth correction + 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 += deltaL; - azimuthR += deltaR; + azimuthL += fx0 * fy; + azimuthR -= fx1 * fy; +#endif } //