dsound: Implement EAX late reverb.

dlls/dsound/dsound_eax.h

@@ -154,6 +154,19 @@ typedef struct {
     DelayLine Decorrelator;
     unsigned int DecoTap[3];
 
+    struct {
+        float Gain;
+        float DensityGain;
+        float MixCoeff;
+
+        float Coeff[4];
+        DelayLine Delay[4];
+        unsigned int Offset[4];
+
+        float LpCoeff[4];
+        float LpSample[4];
+    } Late;
+
     unsigned int Offset;
 } eax_buffer_info;
 

dlls/dsound/eax.c

@@ -107,6 +107,8 @@ static const float LATE_LINE_LENGTH[4] =
 
 static const float LATE_LINE_MULTIPLIER = 4.0f;
 
+#define SPEEDOFSOUNDMETRESPERSEC 343.3f
+
 static float lpFilter2P(FILTER *iir, unsigned int offset, float input)
 {
     float *history = &iir->history[offset*2];
@@ -121,6 +123,11 @@ static float lpFilter2P(FILTER *iir, unsigned int offset, float input)
     return output;
 }
 
+static float lerp(float val1, float val2, float mu)
+{
+    return val1 + (val2-val1)*mu;
+}
+
 static void DelayLineIn(DelayLine *Delay, unsigned int offset, float in)
 {
     Delay->Line[offset&Delay->Mask] = in;
@@ -187,8 +194,85 @@ static void EarlyReflection(IDirectSoundBufferImpl* dsb, float in, float *out)
     out[3] = dsb->eax.Early.Gain * f[3];
 }
 
+static float LateDelayLineOut(IDirectSoundBufferImpl* dsb, unsigned int index)
+{
+    return AttenuatedDelayLineOut(&dsb->eax.Late.Delay[index],
+                                  dsb->eax.Offset - dsb->eax.Late.Offset[index],
+                                  dsb->eax.Late.Coeff[index]);
+}
+
+static float LateLowPassInOut(IDirectSoundBufferImpl* dsb, unsigned int index, float in)
+{
+    in = lerp(in, dsb->eax.Late.LpSample[index], dsb->eax.Late.LpCoeff[index]);
+    dsb->eax.Late.LpSample[index] = in;
+    return in;
+}
+
+static void LateReverb(IDirectSoundBufferImpl* dsb, const float *in, float *out)
+{
+    float d[4], f[4];
+
+    /* Obtain the decayed results of the cyclical delay lines, and add the
+     * corresponding input channels.  Then pass the results through the
+     * low-pass filters. */
+
+    /* This is where the feed-back cycles from line 0 to 1 to 3 to 2 and back
+     * to 0. */
+    d[0] = LateLowPassInOut(dsb, 2, in[2] + LateDelayLineOut(dsb, 2));
+    d[1] = LateLowPassInOut(dsb, 0, in[0] + LateDelayLineOut(dsb, 0));
+    d[2] = LateLowPassInOut(dsb, 3, in[3] + LateDelayLineOut(dsb, 3));
+    d[3] = LateLowPassInOut(dsb, 1, in[1] + LateDelayLineOut(dsb, 1));
+
+    /* Late reverb is done with a modified feed-back delay network (FDN)
+     * topology.  Four input lines are each fed through their own all-pass
+     * filter and then into the mixing matrix.  The four outputs of the
+     * mixing matrix are then cycled back to the inputs.  Each output feeds
+     * a different input to form a circlular feed cycle.
+     *
+     * The mixing matrix used is a 4D skew-symmetric rotation matrix derived
+     * using a single unitary rotational parameter:
+     *
+     *  [  d,  a,  b,  c ]          1 = a^2 + b^2 + c^2 + d^2
+     *  [ -a,  d,  c, -b ]
+     *  [ -b, -c,  d,  a ]
+     *  [ -c,  b, -a,  d ]
+     *
+     * The rotation is constructed from the effect's diffusion parameter,
+     * yielding:  1 = x^2 + 3 y^2; where a, b, and c are the coefficient y
+     * with differing signs, and d is the coefficient x.  The matrix is thus:
+     *
+     *  [  x,  y, -y,  y ]          n = sqrt(matrix_order - 1)
+     *  [ -y,  x,  y,  y ]          t = diffusion_parameter * atan(n)
+     *  [  y, -y,  x,  y ]          x = cos(t)
+     *  [ -y, -y, -y,  x ]          y = sin(t) / n
+     *
+     * To reduce the number of multiplies, the x coefficient is applied with
+     * the cyclical delay line coefficients.  Thus only the y coefficient is
+     * applied when mixing, and is modified to be:  y / x.
+     */
+    f[0] = d[0] + (dsb->eax.Late.MixCoeff * (         d[1] + -d[2] + d[3]));
+    f[1] = d[1] + (dsb->eax.Late.MixCoeff * (-d[0]         +  d[2] + d[3]));
+    f[2] = d[2] + (dsb->eax.Late.MixCoeff * ( d[0] + -d[1]         + d[3]));
+    f[3] = d[3] + (dsb->eax.Late.MixCoeff * (-d[0] + -d[1] + -d[2]       ));
+
+    /* Output the results of the matrix for all four channels, attenuated by
+     * the late reverb gain (which is attenuated by the 'x' mix coefficient). */
+    out[0] = dsb->eax.Late.Gain * f[0];
+    out[1] = dsb->eax.Late.Gain * f[1];
+    out[2] = dsb->eax.Late.Gain * f[2];
+    out[3] = dsb->eax.Late.Gain * f[3];
+
+    /* Re-feed the cyclical delay lines. */
+    DelayLineIn(&dsb->eax.Late.Delay[0], dsb->eax.Offset, f[0]);
+    DelayLineIn(&dsb->eax.Late.Delay[1], dsb->eax.Offset, f[1]);
+    DelayLineIn(&dsb->eax.Late.Delay[2], dsb->eax.Offset, f[2]);
+    DelayLineIn(&dsb->eax.Late.Delay[3], dsb->eax.Offset, f[3]);
+}
+
 static void VerbPass(IDirectSoundBufferImpl* dsb, float in, float* out)
 {
+    float feed, late[4], taps[4];
+
     /* Low-pass filter the incoming sample. */
     in = lpFilter2P(&dsb->eax.LpFilter, 0, in);
 
@@ -199,6 +283,25 @@ static void VerbPass(IDirectSoundBufferImpl* dsb, float in, float* out)
     in = DelayLineOut(&dsb->eax.Delay, dsb->eax.Offset - dsb->eax.DelayTap[0]);
     EarlyReflection(dsb, in, out);
 
+    /* Feed the decorrelator from the energy-attenuated output of the second
+     * delay tap. */
+    in = DelayLineOut(&dsb->eax.Delay, dsb->eax.Offset - dsb->eax.DelayTap[1]);
+    feed = in * dsb->eax.Late.DensityGain;
+    DelayLineIn(&dsb->eax.Decorrelator, dsb->eax.Offset, feed);
+
+    /* Calculate the late reverb from the decorrelator taps. */
+    taps[0] = feed;
+    taps[1] = DelayLineOut(&dsb->eax.Decorrelator, dsb->eax.Offset - dsb->eax.DecoTap[0]);
+    taps[2] = DelayLineOut(&dsb->eax.Decorrelator, dsb->eax.Offset - dsb->eax.DecoTap[1]);
+    taps[3] = DelayLineOut(&dsb->eax.Decorrelator, dsb->eax.Offset - dsb->eax.DecoTap[2]);
+    LateReverb(dsb, taps, late);
+
+    /* Mix early reflections and late reverb. */
+    out[0] += late[0];
+    out[1] += late[1];
+    out[2] += late[2];
+    out[3] += late[3];
+
     /* Step all delays forward one sample. */
     dsb->eax.Offset++;
 }
@@ -308,6 +411,95 @@ static float CalcI3DL2HFreq(float hfRef, unsigned int frequency)
     return cosf(M_PI*2.0f * hfRef / frequency);
 }
 
+static float CalcDensityGain(float a)
+{
+    return sqrtf(1.0f - (a * a));
+}
+
+static float CalcDampingCoeff(float hfRatio, float length, float decayTime, float decayCoeff, float cw)
+{
+    float coeff, g;
+
+    coeff = 0.0f;
+    if (hfRatio < 1.0f)
+    {
+        /* Calculate the low-pass coefficient by dividing the HF decay
+         * coefficient by the full decay coefficient. */
+        g = CalcDecayCoeff(length, decayTime * hfRatio) / decayCoeff;
+
+        /* Damping is done with a 1-pole filter, so g needs to be squared. */
+        g *= g;
+        coeff = lpCoeffCalc(g, cw);
+
+        /* Very low decay times will produce minimal output, so apply an
+         * upper bound to the coefficient. */
+        if (coeff > 0.98f) coeff = 0.98f;
+    }
+    return coeff;
+}
+
+static void UpdateLateLines(float reverbGain, float lateGain, float xMix, float density, float decayTime,
+                            float diffusion, float hfRatio, float cw, unsigned int frequency, eax_buffer_info *State)
+{
+    float length;
+    unsigned int index;
+
+    /* Calculate the late reverb gain (from the master effect gain, and late
+     * reverb gain parameters).  Since the output is tapped prior to the
+     * application of the next delay line coefficients, this gain needs to be
+     * attenuated by the 'x' mixing matrix coefficient as well.
+     */
+    State->Late.Gain = reverbGain * lateGain * xMix;
+
+    /* To compensate for changes in modal density and decay time of the late
+     * reverb signal, the input is attenuated based on the maximal energy of
+     * the outgoing signal.  This approximation is used to keep the apparent
+     * energy of the signal equal for all ranges of density and decay time.
+     *
+     * The average length of the cyclcical delay lines is used to calculate
+     * the attenuation coefficient.
+     */
+    length = (LATE_LINE_LENGTH[0] + LATE_LINE_LENGTH[1] +
+              LATE_LINE_LENGTH[2] + LATE_LINE_LENGTH[3]) / 4.0f;
+    length *= 1.0f + (density * LATE_LINE_MULTIPLIER);
+    State->Late.DensityGain = CalcDensityGain(CalcDecayCoeff(length,
+                                                             decayTime));
+
+    for(index = 0; index < 4; index++)
+    {
+        /* Calculate the length (in seconds) of each cyclical delay line. */
+        length = LATE_LINE_LENGTH[index] * (1.0f + (density * LATE_LINE_MULTIPLIER));
+
+        /* Calculate the delay offset for each cyclical delay line. */
+        State->Late.Offset[index] = fastf2u(length * frequency);
+
+        /* Calculate the gain (coefficient) for each cyclical line. */
+        State->Late.Coeff[index] = CalcDecayCoeff(length, decayTime);
+
+        /* Calculate the damping coefficient for each low-pass filter. */
+        State->Late.LpCoeff[index] = CalcDampingCoeff(hfRatio, length, decayTime,
+                             State->Late.Coeff[index], cw);
+
+        /* Attenuate the cyclical line coefficients by the mixing coefficient
+         * (x). */
+        State->Late.Coeff[index] *= xMix;
+    }
+}
+
+static void CalcMatrixCoeffs(float diffusion, float *x, float *y)
+{
+    float n, t;
+
+    /* The matrix is of order 4, so n is sqrt (4 - 1). */
+    n = sqrtf(3.0f);
+    t = diffusion * atanf(n);
+
+    /* Calculate the first mixing matrix coefficient. */
+    *x = cosf(t);
+    /* Calculate the second mixing matrix coefficient. */
+    *y = sinf(t) / n;
+}
+
 static unsigned int NextPowerOf2(unsigned int value)
 {
     if (value > 0)
@@ -370,6 +562,14 @@ static BOOL AllocLines(unsigned int frequency, IDirectSoundBufferImpl *dsb)
     totalSamples += CalcLineLength(length, totalSamples, frequency,
                                    &dsb->eax.Decorrelator);
 
+    /* The late delay lines are calculated from the lowest reverb density. */
+    for (index = 0; index < 4; index++)
+    {
+        length = LATE_LINE_LENGTH[index] * (1.0f + LATE_LINE_MULTIPLIER);
+        totalSamples += CalcLineLength(length, totalSamples, frequency,
+                                       &dsb->eax.Late.Delay[index]);
+    }
+
     if (totalSamples != dsb->eax.TotalSamples)
     {
         TRACE("New reverb buffer length: %u samples (%f sec)\n", totalSamples, totalSamples/(float)frequency);
@@ -391,6 +591,7 @@ static BOOL AllocLines(unsigned int frequency, IDirectSoundBufferImpl *dsb)
     for(index = 0; index < 4; index++)
     {
         RealizeLineOffset(dsb->eax.SampleBuffer, &dsb->eax.Early.Delay[index]);
+        RealizeLineOffset(dsb->eax.SampleBuffer, &dsb->eax.Late.Delay[index]);
     }
 
     /* Clear the sample buffer. */
@@ -400,10 +601,35 @@ static BOOL AllocLines(unsigned int frequency, IDirectSoundBufferImpl *dsb)
     return TRUE;
 }
 
+static inline float CalcDecayLength(float coeff, float decayTime)
+{
+    return log10f(coeff) * decayTime / log10f(0.001f)/*-60 dB*/;
+}
+
+static float CalcLimitedHfRatio(float hfRatio, float airAbsorptionGainHF, float decayTime)
+{
+    float limitRatio;
+
+    /* Find the attenuation due to air absorption in dB (converting delay
+     * time to meters using the speed of sound).  Then reversing the decay
+     * equation, solve for HF ratio.  The delay length is cancelled out of
+     * the equation, so it can be calculated once for all lines.
+     */
+    limitRatio = 1.0f / (CalcDecayLength(airAbsorptionGainHF, decayTime) *
+                         SPEEDOFSOUNDMETRESPERSEC);
+    /* Using the limit calculated above, apply the upper bound to the HF
+     * ratio. Also need to limit the result to a minimum of 0.1, just like the
+     * HF ratio parameter. */
+    if (limitRatio < 0.1f) limitRatio = 0.1f;
+    else if (limitRatio > hfRatio) limitRatio = hfRatio;
+
+    return limitRatio;
+}
+
 static void ReverbUpdate(IDirectSoundBufferImpl *dsb)
 {
     unsigned int index;
-    float cw;
+    float cw, hfRatio, x, y;
 
     /* avoid segfaults in mixing thread when we recalculate the line offsets */
     EnterCriticalSection(&dsb->device->mixlock);
@@ -430,6 +656,20 @@ static void ReverbUpdate(IDirectSoundBufferImpl *dsb)
                     dsb->device->eax.eax_props.flLateReverbDelay, &dsb->eax);
 
     UpdateDecorrelator(dsb->device->eax.eax_props.flDensity, dsb->device->pwfx->nSamplesPerSec, &dsb->eax);
+
+    CalcMatrixCoeffs(dsb->device->eax.eax_props.flDiffusion, &x, &y);
+    dsb->eax.Late.MixCoeff = y / x;
+
+    hfRatio = dsb->device->eax.eax_props.flDecayHFRatio;
+
+    if (dsb->device->eax.eax_props.iDecayHFLimit && dsb->device->eax.eax_props.flAirAbsorptionGainHF < 1.0f) {
+        hfRatio = CalcLimitedHfRatio(hfRatio, dsb->device->eax.eax_props.flAirAbsorptionGainHF,
+                                     dsb->device->eax.eax_props.flDecayTime);
+    }
+
+    UpdateLateLines(dsb->device->eax.eax_props.flGain, dsb->device->eax.eax_props.flLateReverbGain,
+                    x, dsb->device->eax.eax_props.flDensity, dsb->device->eax.eax_props.flDecayTime,
+                    dsb->device->eax.eax_props.flDiffusion, hfRatio, cw, dsb->device->pwfx->nSamplesPerSec, &dsb->eax);
 }
 
 static BOOL ReverbDeviceUpdate(DirectSoundDevice *dev)
@@ -474,6 +714,20 @@ void init_eax_buffer(IDirectSoundBufferImpl *dsb)
     dsb->eax.DecoTap[1] = 0;
     dsb->eax.DecoTap[2] = 0;
 
+    dsb->eax.Late.Gain = 0.0f;
+    dsb->eax.Late.DensityGain = 0.0f;
+    dsb->eax.Late.MixCoeff = 0.0f;
+    for(index = 0; index < 4; index++)
+    {
+        dsb->eax.Late.Coeff[index] = 0.0f;
+        dsb->eax.Late.Delay[index].Mask = 0;
+        dsb->eax.Late.Delay[index].Line = NULL;
+        dsb->eax.Late.Offset[index] = 0;
+
+        dsb->eax.Late.LpCoeff[index] = 0.0f;
+        dsb->eax.Late.LpSample[index] = 0.0f;
+    }
+
     dsb->eax.Offset = 0;
 
     ReverbUpdate(dsb);