Fade between HRTF coefficients, to reduce noise from sudden changes

3 files changed, 254 insertions, 24 deletions

diff --git a/Alc/ALu.c b/Alc/ALu.c
index 3dc1a1ff..8a81aad1 100644
--- a/Alc/ALu.c
+++ b/Alc/ALu.c
@@ -292,13 +292,17 @@ ALvoid CalcNonAttnSourceParams(ALsource *ALSource, const ALCcontext *ALContext)
                     ALSource->Params.HrtfCoeffs[c][i][0] = 0.0f;
                     ALSource->Params.HrtfCoeffs[c][i][1] = 0.0f;
                 }
-                continue;
             }
-
-            GetLerpedHrtfCoeffs(0.0, angles[c] * (M_PI/180.0),
-                                DryGain*ListenerGain,
-                                ALSource->Params.HrtfCoeffs[c],
-                                ALSource->Params.HrtfDelay[c]);
+            else
+            {
+                /* Get the static HRIR coefficients and delays for this
+                 * channel. */
+                GetLerpedHrtfCoeffs(0.0, angles[c] * (M_PI/180.0),
+                                    DryGain*ListenerGain,
+                                    ALSource->Params.HrtfCoeffs[c],
+                                    ALSource->Params.HrtfDelay[c]);
+            }
+            ALSource->HrtfCounter = 0;
         }
     }
     else
@@ -709,24 +713,55 @@ ALvoid CalcSourceParams(ALsource *ALSource, const ALCcontext *ALContext)
     if((Device->Flags&DEVICE_USE_HRTF))
     {
         // Use a binaural HRTF algorithm for stereo headphone playback
+        ALfloat delta, ev = 0.0f, az = 0.0f;
+
         if(Distance > 0.0f)
         {
             ALfloat invlen = 1.0f/Distance;
             Position[0] *= invlen;
             Position[1] *= invlen;
             Position[2] *= invlen;
-            GetLerpedHrtfCoeffs(asin(Position[1]),
-                                atan2(Position[0], -Position[2]*ZScale),
-                                DryGain, ALSource->Params.HrtfCoeffs[0],
-                                ALSource->Params.HrtfDelay[0]);
+
+            // Calculate elevation and azimuth only when the source is not at
+            // the listener.  This prevents +0 and -0 Z from producing
+            // inconsistent panning.
+            ev = asin(Position[1]);
+            az = atan2(Position[0], -Position[2]*ZScale);
+        }
+
+        // Check to see if the HRIR is already moving.
+        if(ALSource->HrtfMoving)
+        {
+            // Calculate the normalized HRTF transition factor (delta).
+            delta = CalcHrtfDelta(ALSource->Params.HrtfGain, DryGain,
+                                  ALSource->Params.HrtfDir, Position);
+            // If the delta is large enough, get the moving HRIR target
+            // coefficients, target delays, steppping values, and counter.
+            if(delta > 0.001f)
+            {
+                ALSource->HrtfCounter = GetMovingHrtfCoeffs(ev, az, DryGain,
+                                          delta, ALSource->HrtfCounter,
+                                          ALSource->Params.HrtfCoeffs[0],
+                                          ALSource->Params.HrtfDelay[0],
+                                          ALSource->Params.HrtfCoeffStep,
+                                          ALSource->Params.HrtfDelayStep);
+                ALSource->Params.HrtfGain = DryGain;
+                ALSource->Params.HrtfDir[0] = Position[0];
+                ALSource->Params.HrtfDir[1] = Position[1];
+                ALSource->Params.HrtfDir[2] = Position[2];
+            }
         }
         else
         {
-            /* Force front-centered for sounds that comes from the listener,
-             * to prevent +0 and -0 Z from producing inconsistent panning */
-            GetLerpedHrtfCoeffs(0.0f, 0.0f, DryGain,
+            // Get the initial (static) HRIR coefficients and delays.
+            GetLerpedHrtfCoeffs(ev, az, DryGain,
                                 ALSource->Params.HrtfCoeffs[0],
                                 ALSource->Params.HrtfDelay[0]);
+            ALSource->HrtfCounter = 0;
+            ALSource->Params.HrtfGain = DryGain;
+            ALSource->Params.HrtfDir[0] = Position[0];
+            ALSource->Params.HrtfDir[1] = Position[1];
+            ALSource->Params.HrtfDir[2] = Position[2];
         }
     }
     else
diff --git a/Alc/hrtf.c b/Alc/hrtf.c
index 4e03e198..6edc37e7 100644
--- a/Alc/hrtf.c
+++ b/Alc/hrtf.c
@@ -80,6 +80,37 @@ static void CalcAzIndices(ALuint evidx, ALfloat az, ALuint *azidx, ALfloat *azmu
     *azmu = az - (ALuint)az;
 }
 
+// Calculates the normalized HRTF transition factor (delta) from the changes
+// in gain and listener to source angle between updates.  The result is a
+// normalized delta factor than can be used to calculate moving HRIR stepping
+// values.
+ALfloat CalcHrtfDelta(ALfloat oldGain, ALfloat newGain, const ALfloat olddir[3], const ALfloat newdir[3])
+{
+    ALfloat gainChange, angleChange, delta;
+
+    // Calculate the normalized dB gain change.
+    gainChange = aluFabs((log10(__max(newGain, 0.0001f)) -
+                          log10(__max(oldGain, 0.0001f))) / log10(0.0001f));
+
+    // Calculate the normalized listener to source angle change when there is
+    // enough gain to notice it.
+    angleChange = 0.0f;
+    if(gainChange > 0.0001f || newGain > 0.0001f)
+    {
+        // No angle change when the directions are equal or degenerate (when
+        // both have zero length).
+        if(newdir[0]-olddir[0] || newdir[1]-olddir[1] || newdir[2]-olddir[2])
+            angleChange = aluAcos(olddir[0]*newdir[0] +
+                                  olddir[1]*newdir[1] +
+                                  olddir[2]*newdir[2]) / M_PI;
+
+    }
+
+    // Use the largest of the two changes for the delta factor.
+    delta = __max(gainChange, angleChange) * 2.0f;
+    return __min(delta, 1.0f);
+}
+
 // Calculates static HRIR coefficients and delays for the given polar
 // elevation and azimuth in radians.  Linear interpolation is used to
 // increase the apparent resolution of the HRIR dataset.  The coefficients
@@ -144,10 +175,115 @@ void GetLerpedHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALflo
     // Calculate the HRIR delays using linear interpolation.
     delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
                               lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
-                              mu[2]) + 0.5f);
+                              mu[2]) * 65536.0f);
+    delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
+                              lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
+                              mu[2]) * 65536.0f);
+}
+
+// Calculates the moving HRIR target coefficients, target delays, and
+// stepping values for the given polar elevation and azimuth in radians.
+// Linear interpolation is used to increase the apparent resolution of the
+// HRIR dataset.  The coefficients are also normalized and attenuated by the
+// specified gain.  Stepping resolution and count is determined using the
+// given delta factor between 0.0 and 1.0.
+ALint GetMovingHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat delta, ALint counter, ALfloat (*coeffs)[2], ALuint *delays, ALfloat (*coeffStep)[2], ALint *delayStep)
+{
+    ALfloat step;
+    ALuint evidx[2], azidx[2];
+    ALfloat mu[3];
+    ALuint lidx[4], ridx[4];
+    ALuint i;
+    ALfloat left, right;
+
+    // Claculate elevation indices and interpolation factor.
+    CalcEvIndices(elevation, evidx, &mu[2]);
+
+    // Calculate azimuth indices and interpolation factor for the first
+    // elevation.
+    CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
+
+    // Calculate the first set of linear HRIR indices for left and right
+    // channels.
+    lidx[0] = evOffset[evidx[0]] + azidx[0];
+    lidx[1] = evOffset[evidx[0]] + azidx[1];
+    ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
+    ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
+
+    // Calculate azimuth indices and interpolation factor for the second
+    // elevation.
+    CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
+
+    // Calculate the second set of linear HRIR indices for left and right
+    // channels.
+    lidx[2] = evOffset[evidx[1]] + azidx[0];
+    lidx[3] = evOffset[evidx[1]] + azidx[1];
+    ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
+    ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
+
+    // Calculate the stepping parameters.
+    delta = __max(0.015f * delta * Hrtf.sampleRate, 1.0f);
+    step = 1.0f / delta;
+
+    // Calculate the normalized and attenuated target HRIR coefficients using
+    // linear interpolation when there is enough gain to warrant it.  Zero
+    // the target coefficients if gain is too low.  Then calculate the
+    // coefficient stepping values using the target and previous running
+    // coefficients.
+    if(gain > 0.0001f)
+    {
+        ALdouble scale = gain * (1.0/32767.0);
+        for(i = 0;i < HRIR_LENGTH;i++)
+        {
+            left = coeffs[i][0] - (coeffStep[i][0] * counter);
+            right = coeffs[i][1] - (coeffStep[i][1] * counter);
+
+            coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
+                                lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
+                                mu[2]) * scale;
+            coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
+                                lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
+                                mu[2]) * scale;
+
+            coeffStep[i][0] = step * (coeffs[i][0] - left);
+            coeffStep[i][1] = step * (coeffs[i][1] - right);
+        }
+    }
+    else
+    {
+        for(i = 0;i < HRIR_LENGTH;i++)
+        {
+            left = coeffs[i][0] - (coeffStep[i][0] * counter);
+            right = coeffs[i][1] - (coeffStep[i][1] * counter);
+
+            coeffs[i][0] = 0.0f;
+            coeffs[i][1] = 0.0f;
+
+            coeffStep[i][0] = step * -left;
+            coeffStep[i][1] = step * -right;
+        }
+    }
+
+    // Calculate the HRIR delays using linear interpolation.  Then calculate
+    // the delay stepping values using the target and previous running
+    // delays.
+    left = delays[0] - (delayStep[0] * counter);
+    right = delays[1] - (delayStep[1] * counter);
+
+    delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
+                              lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
+                              mu[2]) * 65536.0f);
     delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
                               lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
-                              mu[2]) + 0.5f);
+                              mu[2]) * 65536.0f);
+
+    delayStep[0] = (ALint)(step * (delays[0] - left));
+    delayStep[1] = (ALint)(step * (delays[1] - right));
+
+    // The stepping count is the number of samples necessary for the HRIR to
+    // complete its transition.  The mixer will only apply stepping for this
+    // many samples.
+    return (ALuint)delta;
 }
 
 ALCboolean IsHrtfCompatible(ALCdevice *device)
diff --git a/Alc/mixer.c b/Alc/mixer.c
index fd351f19..665ef3f8 100644
--- a/Alc/mixer.c
+++ b/Alc/mixer.c
@@ -76,8 +76,10 @@ static void Mix_Hrtf_##T##_##sampler(ALsource *Source, ALCdevice *Device,     \
 {                                                                             \
     const ALuint NumChannels = Source->NumChannels;                           \
     const T *RESTRICT data = srcdata;                                         \
+    const ALint *RESTRICT DelayStep = Source->Params.HrtfDelayStep;           \
     ALfloat (*RESTRICT DryBuffer)[MAXCHANNELS];                               \
     ALfloat *RESTRICT ClickRemoval, *RESTRICT PendingClicks;                  \
+    ALfloat (*RESTRICT CoeffStep)[2] = Source->Params.HrtfCoeffStep;          \
     ALuint pos, frac;                                                         \
     FILTER *DryFilter;                                                        \
     ALuint BufferIdx;                                                         \
@@ -97,38 +99,85 @@ static void Mix_Hrtf_##T##_##sampler(ALsource *Source, ALCdevice *Device,     \
                                                                               \
     for(i = 0;i < NumChannels;i++)                                            \
     {                                                                         \
-        ALfloat (*RESTRICT Coeffs)[2] = Source->Params.HrtfCoeffs[i];         \
-        const ALuint *RESTRICT Delay = Source->Params.HrtfDelay[i];           \
+        ALfloat (*RESTRICT TargetCoeffs)[2] = Source->Params.HrtfCoeffs[i];   \
+        ALuint *RESTRICT TargetDelay = Source->Params.HrtfDelay[i];           \
         ALfloat *RESTRICT History = Source->HrtfHistory[i];                   \
         ALfloat (*RESTRICT Values)[2] = Source->HrtfValues[i];                \
+        ALint Counter = __max(Source->HrtfCounter, OutPos) - OutPos;          \
         ALuint Offset = Source->HrtfOffset + OutPos;                          \
+        ALfloat Coeffs[HRIR_LENGTH][2];                                       \
+        ALuint Delay[2];                                                      \
         ALfloat left, right;                                                  \
                                                                               \
         pos = 0;                                                              \
         frac = *DataPosFrac;                                                  \
                                                                               \
+        for(c = 0;c < HRIR_LENGTH;c++)                                        \
+        {                                                                     \
+            Coeffs[c][0] = TargetCoeffs[c][0] - (CoeffStep[c][0]*Counter);    \
+            Coeffs[c][1] = TargetCoeffs[c][1] - (CoeffStep[c][1]*Counter);    \
+        }                                                                     \
+                                                                              \
+        Delay[0] = TargetDelay[0] - (DelayStep[0]*Counter) + 32768;           \
+        Delay[1] = TargetDelay[1] - (DelayStep[1]*Counter) + 32768;           \
+                                                                              \
         if(LIKELY(OutPos == 0))                                               \
         {                                                                     \
             value = sampler(data + pos*NumChannels + i, NumChannels, frac);   \
             value = lpFilter2PC(DryFilter, i, value);                         \
                                                                               \
             History[Offset&SRC_HISTORY_MASK] = value;                         \
-            left = History[(Offset-Delay[0])&SRC_HISTORY_MASK];               \
-            right = History[(Offset-Delay[1])&SRC_HISTORY_MASK];              \
+            left = History[(Offset-(Delay[0]>>16))&SRC_HISTORY_MASK];         \
+            right = History[(Offset-(Delay[1]>>16))&SRC_HISTORY_MASK];        \
                                                                               \
             ClickRemoval[FRONT_LEFT]  -= Values[(Offset+1)&HRIR_MASK][0] +    \
                                          Coeffs[0][0] * left;                 \
             ClickRemoval[FRONT_RIGHT] -= Values[(Offset+1)&HRIR_MASK][1] +    \
                                          Coeffs[0][1] * right;                \
         }                                                                     \
-        for(BufferIdx = 0;BufferIdx < BufferSize;BufferIdx++)                 \
+        for(BufferIdx = 0;BufferIdx < BufferSize && Counter > 0;BufferIdx++)  \
+        {                                                                     \
+            value = sampler(data + pos*NumChannels + i, NumChannels, frac);   \
+            value = lpFilter2P(DryFilter, i, value);                          \
+                                                                              \
+            History[Offset&SRC_HISTORY_MASK] = value;                         \
+            left = History[(Offset-(Delay[0]>>16))&SRC_HISTORY_MASK];         \
+            right = History[(Offset-(Delay[1]>>16))&SRC_HISTORY_MASK];        \
+                                                                              \
+            Delay[0] += DelayStep[0];                                         \
+            Delay[1] += DelayStep[1];                                         \
+                                                                              \
+            Values[Offset&HRIR_MASK][0] = 0.0f;                               \
+            Values[Offset&HRIR_MASK][1] = 0.0f;                               \
+            Offset++;                                                         \
+                                                                              \
+            for(c = 0;c < HRIR_LENGTH;c++)                                    \
+            {                                                                 \
+                const ALuint off = (Offset+c)&HRIR_MASK;                      \
+                Values[off][0] += Coeffs[c][0] * left;                        \
+                Values[off][1] += Coeffs[c][1] * right;                       \
+                Coeffs[c][0] += CoeffStep[c][0];                              \
+                Coeffs[c][1] += CoeffStep[c][1];                              \
+            }                                                                 \
+                                                                              \
+            DryBuffer[OutPos][FRONT_LEFT]  += Values[Offset&HRIR_MASK][0];    \
+            DryBuffer[OutPos][FRONT_RIGHT] += Values[Offset&HRIR_MASK][1];    \
+                                                                              \
+            frac += increment;                                                \
+            pos  += frac>>FRACTIONBITS;                                       \
+            frac &= FRACTIONMASK;                                             \
+            OutPos++;                                                         \
+            Counter--;                                                        \
+        }                                                                     \
+                                                                              \
+        for(;BufferIdx < BufferSize;BufferIdx++)                              \
         {                                                                     \
             value = sampler(data + pos*NumChannels + i, NumChannels, frac);   \
             value = lpFilter2P(DryFilter, i, value);                          \
                                                                               \
             History[Offset&SRC_HISTORY_MASK] = value;                         \
-            left = History[(Offset-Delay[0])&SRC_HISTORY_MASK];               \
-            right = History[(Offset-Delay[1])&SRC_HISTORY_MASK];              \
+            left = History[(Offset-(Delay[0]>>16))&SRC_HISTORY_MASK];         \
+            right = History[(Offset-(Delay[1]>>16))&SRC_HISTORY_MASK];        \
                                                                               \
             Values[Offset&HRIR_MASK][0] = 0.0f;                               \
             Values[Offset&HRIR_MASK][1] = 0.0f;                               \
@@ -155,8 +204,8 @@ static void Mix_Hrtf_##T##_##sampler(ALsource *Source, ALCdevice *Device,     \
             value = lpFilter2PC(DryFilter, i, value);                         \
                                                                               \
             History[Offset&SRC_HISTORY_MASK] = value;                         \
-            left = History[(Offset-Delay[0])&SRC_HISTORY_MASK];               \
-            right = History[(Offset-Delay[1])&SRC_HISTORY_MASK];              \
+            left = History[(Offset-(Delay[0]>>16))&SRC_HISTORY_MASK];         \
+            right = History[(Offset-(Delay[1]>>16))&SRC_HISTORY_MASK];        \
                                                                               \
             PendingClicks[FRONT_LEFT]  += Values[(Offset+1)&HRIR_MASK][0] +   \
                                           Coeffs[0][0] * left;                \
@@ -743,4 +792,14 @@ ALvoid MixSource(ALsource *Source, ALCdevice *Device, ALuint SamplesToDo)
     Source->position_fraction = DataPosFrac;
     Source->Buffer            = BufferListItem->buffer;
     Source->HrtfOffset       += OutPos;
+    if(State == AL_PLAYING)
+    {
+        Source->HrtfCounter = __max(Source->HrtfCounter, OutPos) - OutPos;
+        Source->HrtfMoving  = AL_TRUE;
+    }
+    else
+    {
+        Source->HrtfCounter = 0;
+        Source->HrtfMoving  = AL_FALSE;
+    }
 }