Merge pull request #12479 from kencooke/audio-nearfield-hrtf

Near-field HRTF
This commit is contained in:
Ken Cooke 2018-02-28 07:44:24 -08:00 committed by GitHub
commit 328b5a1e8b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 260 additions and 62 deletions

View file

@ -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;
}

View file

@ -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;
}

View file

@ -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) {

View file

@ -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: