Use a modified reverb model that obeys the reverb parameters better

1 files changed, 323 insertions, 222 deletions

diff --git a/Alc/alcReverb.c b/Alc/alcReverb.c
index b4eaec80..0a65920f 100644
--- a/Alc/alcReverb.c
+++ b/Alc/alcReverb.c
@@ -1,6 +1,6 @@
 /**
- * OpenAL cross platform audio library
- * Copyright (C) 2008 by Christopher Fitzgerald.
+ * Reverb for the OpenAL cross platform audio library
+ * Copyright (C) 2008-2009 by Christopher Fitzgerald.
  * This library is free software; you can redistribute it and/or
  *  modify it under the terms of the GNU Library General Public
  *  License as published by the Free Software Foundation; either
@@ -46,8 +46,8 @@
 
 typedef struct DelayLine
 {
-    // The delay lines use lengths that are powers of 2 to allow bitmasking
-    // instead of modulus wrapping.
+    // The delay lines use sample lengths that are powers of 2 to allow
+    // bitmasking instead of modulus wrapping.
     ALuint   Mask;
     ALfloat *Line;
 } DelayLine;
@@ -55,15 +55,15 @@ typedef struct DelayLine
 struct ALverbState
 {
     // All delay lines are allocated as a single buffer to reduce memory
-    // fragmentation and teardown code.
+    // fragmentation and management code.
     ALfloat  *SampleBuffer;
-    // Master reverb gain.
+    // Master effect gain.
     ALfloat   Gain;
-    // Initial reverb delay.
+    // Initial effect delay and decorrelation.
     DelayLine Delay;
     // The tap points for the initial delay.  First tap goes to early
-    // reflections, the second to late reverb.
-    ALuint    Tap[2];
+    // reflections, the last four decorrelate to late reverb.
+    ALuint    Tap[5];
     struct {
         // Gain for early reflections.
         ALfloat   Gain;
@@ -75,42 +75,64 @@ struct ALverbState
     struct {
         // Gain for late reverb.
         ALfloat   Gain;
-        // Diffusion of late reverb.
-        ALfloat   Diffusion;
-        // Late reverb is done with 8 delay lines.
-        ALfloat   Coeff[8];
-        DelayLine Delay[8];
-        ALuint    Offset[8];
-        // The input and last 4 delay lines are low-pass filtered.
-        ALfloat   LpCoeff[5];
-        ALfloat   LpSample[5];
+        // Attenuation to compensate for modal density and decay rate.
+        ALfloat   DensityGain;
+        // The feed-back and feed-forward all-pass coefficient.
+        ALfloat   ApFeedCoeff;
+        // Mixing matrix coefficient.
+        ALfloat   MixCoeff;
+        // Late reverb has 4 parallel all-pass filters.
+        ALfloat   ApCoeff[4];
+        DelayLine ApDelay[4];
+        ALuint    ApOffset[4];
+        // In addition to 4 cyclical delay lines.
+        ALfloat   Coeff[4];
+        DelayLine Delay[4];
+        ALuint    Offset[4];
+        // The cyclical delay lines are low-pass filtered.
+        ALfloat   LpCoeff[4][2];
+        ALfloat   LpSample[4];
     } Late;
+    // The current read offset for all delay lines.
     ALuint Offset;
 };
 
 // All delay line lengths are specified in seconds.
 
-// The length of the initial delay line (a sum of the maximum delay before
-// early reflections and late reverb; 0.3 + 0.1).
-static const ALfloat MASTER_LINE_LENGTH = 0.4000f;
-
 // The lengths of the early delay lines.
 static const ALfloat EARLY_LINE_LENGTH[4] =
 {
     0.0015f, 0.0045f, 0.0135f, 0.0405f
 };
 
-// The lengths of the late delay lines.
-static const ALfloat LATE_LINE_LENGTH[8] =
+// The lengths of the late all-pass delay lines.
+static const ALfloat ALLPASS_LINE_LENGTH[4] =
 {
-    0.0015f, 0.0037f, 0.0093f, 0.0234f,
-    0.0100f, 0.0150f, 0.0225f, 0.0337f
+    0.0151f, 0.0167f, 0.0183f, 0.0200f,
 };
 
-// The last 4 late delay lines have a variable length dependent on the effect
-// density parameter and this multiplier.
-static const ALfloat LATE_LINE_MULTIPLIER = 9.0f;
+// The lengths of the late cyclical delay lines.
+static const ALfloat LATE_LINE_LENGTH[4] =
+{
+    0.0211f, 0.0311f, 0.0461f, 0.0680f
+};
+
+// The late cyclical delay lines have a variable length dependent on the
+// effect's density parameter (inverted for some reason) and this multiplier.
+static const ALfloat LATE_LINE_MULTIPLIER = 4.0f;
+
+// Input into the late reverb is decorrelated between four channels.  Their
+// timings are dependent on a fraction and multiplier.  See VerbUpdate() for
+// the calculations involved.
+static const ALfloat DECO_FRACTION = 1.0f / 32.0f;
+static const ALfloat DECO_MULTIPLIER = 2.0f;
 
+// The maximum length of initial delay for the master delay line (a sum of
+// the maximum early reflection and late reverb delays).
+static const ALfloat MASTER_LINE_LENGTH = 0.3f + 0.1f;
+
+// Find the next power of 2.  Actually, this will return the input value if
+// it is already a power of 2.
 static ALuint NextPowerOf2(ALuint value)
 {
     ALuint powerOf2 = 1;
@@ -146,8 +168,8 @@ static __inline ALfloat EarlyDelayLineOut(ALverbState *State, ALuint index)
                         State->Offset - State->Early.Offset[index]);
 }
 
-// Given an input sample, this function produces a decorrelated stereo output
-// for early reflections.
+// Given an input sample, this function produces stereo output for early
+// reflections.
 static __inline ALvoid EarlyReflection(ALverbState *State, ALfloat in, ALfloat *out)
 {
     ALfloat d[4], v, f[4];
@@ -161,11 +183,11 @@ static __inline ALvoid EarlyReflection(ALverbState *State, ALfloat in, ALfloat *
     /* The following uses a lossless scattering junction from waveguide
      * theory.  It actually amounts to a householder mixing matrix, which
      * will produce a maximally diffuse response, and means this can probably
-     * be considered a simple FDN.
+     * be considered a simple feedback delay network (FDN).
      *          N
      *         ---
      *         \
-     * v = 2/N /   di
+     * v = 2/N /   d_i
      *         ---
      *         i=1
      */
@@ -194,6 +216,20 @@ static __inline ALvoid EarlyReflection(ALverbState *State, ALfloat in, ALfloat *
     out[1] = State->Early.Gain * f[3];
 }
 
+// All-pass input/output routine for late reverb.
+static __inline ALfloat LateAllPassInOut(ALverbState *State, ALuint index, ALfloat in)
+{
+    ALfloat out;
+
+    out = State->Late.ApCoeff[index] *
+          DelayLineOut(&State->Late.ApDelay[index],
+                       State->Offset - State->Late.ApOffset[index]);
+    out -= (State->Late.ApFeedCoeff * in);
+    DelayLineIn(&State->Late.ApDelay[index], State->Offset,
+                (State->Late.ApFeedCoeff * out) + in);
+    return out;
+}
+
 // Delay line output routine for late reverb.
 static __inline ALfloat LateDelayLineOut(ALverbState *State, ALuint index)
 {
@@ -205,94 +241,81 @@ static __inline ALfloat LateDelayLineOut(ALverbState *State, ALuint index)
 // Low-pass filter input/output routine for late reverb.
 static __inline ALfloat LateLowPassInOut(ALverbState *State, ALuint index, ALfloat in)
 {
-    State->Late.LpSample[index] = in + ((State->Late.LpSample[index] - in) *
-                                        State->Late.LpCoeff[index]);
+    State->Late.LpSample[index] = (State->Late.LpCoeff[index][0] * in) +
+        (State->Late.LpCoeff[index][1] * State->Late.LpSample[index]);
     return State->Late.LpSample[index];
 }
 
-// Given an input sample, this function produces a decorrelated stereo output
-// for late reverb.
-static __inline ALvoid LateReverb(ALverbState *State, ALfloat in, ALfloat *out)
+// Given four decorrelated input samples, this function produces stereo
+// output for late reverb.
+static __inline ALvoid LateReverb(ALverbState *State, ALfloat *in, ALfloat *out)
 {
-    ALfloat din, d[8], v, dv, f[8];
-
-    // Since the input will be sent directly to the output as in the early
-    // reflections function, it needs to take into account some immediate
-    // absorption.
-    in = LateLowPassInOut(State, 0, in);
-
-    // When diffusion is full, no input is directly passed to the variable-
-    // length delay lines (the last 4).
-    din = (1.0f - State->Late.Diffusion) * in;
-
-    // Obtain the decayed results of the fixed-length delay lines.
-    d[0] = LateDelayLineOut(State, 0);
-    d[1] = LateDelayLineOut(State, 1);
-    d[2] = LateDelayLineOut(State, 2);
-    d[3] = LateDelayLineOut(State, 3);
-    // Obtain the decayed and low-pass filtered results of the variable-
-    // length delay lines.
-    d[4] = LateLowPassInOut(State, 1, LateDelayLineOut(State, 4));
-    d[5] = LateLowPassInOut(State, 2, LateDelayLineOut(State, 5));
-    d[6] = LateLowPassInOut(State, 3, LateDelayLineOut(State, 6));
-    d[7] = LateLowPassInOut(State, 4, LateDelayLineOut(State, 7));
-
-    // The waveguide formula used in the early reflections function works
-    // great for high diffusion, but it is not obviously paramerized to allow
-    // a variable diffusion.  With only limited time and resources, what
-    // follows is the best variation of that formula I could come up with.
-    // First, there are 8 delay lines used.  The first 4 are fixed-length and
-    // generate the highest density of the diffuse response.  The last 4 are
-    // variable-length, and are used to smooth out the diffuse response.  The
-    // density effect parameter alters their length.  The inner two delay
-    // lines of each group have their signs reversed (more about this later).
-    v = (d[0] - d[1] - d[2] + d[3] +
-         d[4] - d[5] - d[6] + d[7]) * 0.25f;
-    // Diffusion is applied as a reduction of the junction pressure for all
-    // branches.  This presents two problems.  When the diffusion factor (0
-    // to 1) reaches 0.5, the average feed value is reduced (the junction
-    // becomes lossy).  Thus, at 0.5 the signal decays almost twice as fast
-    // as it should.  The second problem is the introduction of some
-    // resonant frequencies (coloration).  The reversed signs above are used
-    // to help combat some of the coloration by adding variations along the
-    // feed cycle.
-    v *= State->Late.Diffusion;
-    // Load the junction with the input.  To reduce the noticeable echo of
-    // the longer delay lines (the variable-length ones) the input is loaded
-    // with the inverse of the effect diffusion.  So at full diffusion, the
-    // input is not applied to the last 4 delay lines.  Input signs reversed
-    // to balance the equation.
-    dv = v + din;
-    v += in;
-
-    // As with the reversed signs above, to balance the equation the signs
-    // need to be reversed here, too.
-    f[0] = d[0] - v;
-    f[1] = d[1] + v;
-    f[2] = d[2] + v;
-    f[3] = d[3] - v;
-    f[4] = d[4] - dv;
-    f[5] = d[5] + dv;
-    f[6] = d[6] + dv;
-    f[7] = d[7] - dv;
-
-    // Feed the fixed-length delay lines with their own cycle (0 -> 1 -> 3 ->
-    // 2 -> 0...).
+    ALfloat d[4], f[4];
+
+    // Obtain the decayed results of the cyclical delay lines, and add the
+    // corresponding input channels attenuated by density.  Then pass the
+    // results through the low-pass filters.
+    d[0] = LateLowPassInOut(State, 0, (State->Late.DensityGain * in[0]) +
+                                      LateDelayLineOut(State, 0));
+    d[1] = LateLowPassInOut(State, 1, (State->Late.DensityGain * in[1]) +
+                                      LateDelayLineOut(State, 1));
+    d[2] = LateLowPassInOut(State, 2, (State->Late.DensityGain * in[2]) +
+                                      LateDelayLineOut(State, 2));
+    d[3] = LateLowPassInOut(State, 3, (State->Late.DensityGain * in[3]) +
+                                      LateDelayLineOut(State, 3));
+
+    // To help increase diffusion, run each line through an all-pass filter.
+    // The order of the all-pass filters is selected so that the shortest
+    // all-pass filter will feed the shortest delay line.
+    d[0] = LateAllPassInOut(State, 1, d[0]);
+    d[1] = LateAllPassInOut(State, 3, d[1]);
+    d[2] = LateAllPassInOut(State, 0, d[2]);
+    d[3] = LateAllPassInOut(State, 2, d[3]);
+
+    /* Late reverb is done with a modified feedback 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 ]          x = 1 - (0.5 diffusion^3)
+     *  [ -y,  x,  y,  y ]          y = sqrt((1 - x^2) / 3)
+     *  [  y, -y,  x,  y ]
+     *  [ -y, -y, -y,  x ]
+     *
+     * 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] + (State->Late.MixCoeff * ( d[1] - d[2] + d[3]));
+    f[1] = d[1] + (State->Late.MixCoeff * (-d[0] + d[2] + d[3]));
+    f[2] = d[2] + (State->Late.MixCoeff * ( d[0] - d[1] + d[3]));
+    f[3] = d[3] + (State->Late.MixCoeff * (-d[0] - d[1] - d[2]));
+
+    // Output is tapped at the input to the shortest two cyclical delay
+    // lines, attenuated by the late reverb gain (which is attenuated by the
+    // mixing coefficient x).
+    out[0] = State->Late.Gain * f[0];
+    out[1] = State->Late.Gain * f[1];
+
+    // The delay lines are fed circularly in the order:
+    // 0 -> 1 -> 3 -> 2 -> 0 ...
     DelayLineIn(&State->Late.Delay[0], State->Offset, f[2]);
     DelayLineIn(&State->Late.Delay[1], State->Offset, f[0]);
     DelayLineIn(&State->Late.Delay[2], State->Offset, f[3]);
     DelayLineIn(&State->Late.Delay[3], State->Offset, f[1]);
-    // Feed the variable-length delay lines with their cycle (4 -> 6 -> 7 ->
-    // 5 -> 4...).
-    DelayLineIn(&State->Late.Delay[4], State->Offset, f[5]);
-    DelayLineIn(&State->Late.Delay[5], State->Offset, f[7]);
-    DelayLineIn(&State->Late.Delay[6], State->Offset, f[4]);
-    DelayLineIn(&State->Late.Delay[7], State->Offset, f[6]);
-
-    // Output is derived from the values fed to the inner two variable-length
-    // delay lines (5 and 6).
-    out[0] = State->Late.Gain * f[7];
-    out[1] = State->Late.Gain * f[4];
 }
 
 // This creates the reverb state.  It should be called only when the reverb
@@ -300,33 +323,46 @@ static __inline ALvoid LateReverb(ALverbState *State, ALfloat in, ALfloat *out)
 ALverbState *VerbCreate(ALCcontext *Context)
 {
     ALverbState *State = NULL;
-    ALuint length[13], totalLength, index;
+    ALuint samples, length[13], totalLength, index;
 
     State = malloc(sizeof(ALverbState));
     if(!State)
         return NULL;
 
-    // All line lengths are powers of 2, calculated from the line timings and
-    // the addition of an extra sample (for safety).
-    length[0] = NextPowerOf2((ALuint)(MASTER_LINE_LENGTH*Context->Frequency) + 1);
+    // All line lengths are powers of 2, calculated from their lengths, with
+    // an additional sample in case of rounding errors.
+
+    // See VerbUpdate() for an explanation of the additional calculation
+    // added to the master line length.
+    samples = (ALuint)
+              ((MASTER_LINE_LENGTH +
+                (LATE_LINE_LENGTH[0] * (1.0f + LATE_LINE_MULTIPLIER) *
+                 (DECO_FRACTION * ((DECO_MULTIPLIER * DECO_MULTIPLIER *
+                                    DECO_MULTIPLIER) - 1.0f)))) *
+               Context->Frequency) + 1;
+    length[0] = NextPowerOf2(samples);
     totalLength = length[0];
     for(index = 0;index < 4;index++)
     {
-        length[1+index] = NextPowerOf2((ALuint)(EARLY_LINE_LENGTH[index]*Context->Frequency) + 1);
-        totalLength += length[1+index];
+        samples = (ALuint)(EARLY_LINE_LENGTH[index] * Context->Frequency) + 1;
+        length[1 + index] = NextPowerOf2(samples);
+        totalLength += length[1 + index];
     }
     for(index = 0;index < 4;index++)
     {
-        length[5+index] = NextPowerOf2((ALuint)(LATE_LINE_LENGTH[index]*Context->Frequency) + 1);
-        totalLength += length[5+index];
+        samples = (ALuint)(ALLPASS_LINE_LENGTH[index] * Context->Frequency) + 1;
+        length[5 + index] = NextPowerOf2(samples);
+        totalLength += length[5 + index];
     }
-    for(index = 4;index < 8;index++)
+    for(index = 0;index < 4;index++)
     {
-        length[5+index] = NextPowerOf2((ALuint)(LATE_LINE_LENGTH[index]*(1.0f + LATE_LINE_MULTIPLIER)*Context->Frequency) + 1);
-        totalLength += length[5+index];
+        samples = (ALuint)(LATE_LINE_LENGTH[index] *
+                           (1.0f + LATE_LINE_MULTIPLIER) * Context->Frequency) + 1;
+        length[9 + index] = NextPowerOf2(samples);
+        totalLength += length[9 + index];
     }
 
-    // They all share a single sample buffer.
+    // All lines share a single sample buffer.
     State->SampleBuffer = malloc(totalLength * sizeof(ALfloat));
     if(!State->SampleBuffer)
     {
@@ -344,10 +380,11 @@ ALverbState *VerbCreate(ALCcontext *Context)
 
     State->Tap[0] = 0;
     State->Tap[1] = 0;
+    State->Tap[2] = 0;
+    State->Tap[3] = 0;
+    State->Tap[4] = 0;
 
     State->Early.Gain = 0.0f;
-    // All fixed-length delay lines have their read-write offsets calculated
-    // one time.
     for(index = 0;index < 4;index++)
     {
         State->Early.Coeff[index] = 0.0f;
@@ -355,30 +392,40 @@ ALverbState *VerbCreate(ALCcontext *Context)
         State->Early.Delay[index].Line = &State->SampleBuffer[totalLength];
         totalLength += length[1 + index];
 
-        State->Early.Offset[index] = (ALuint)(EARLY_LINE_LENGTH[index] * Context->Frequency);
+        // The early delay lines have their read offsets calculated once.
+        State->Early.Offset[index] = (ALuint)(EARLY_LINE_LENGTH[index] *
+                                              Context->Frequency);
     }
 
     State->Late.Gain = 0.0f;
-    State->Late.Diffusion = 0.0f;
-    for(index = 0;index < 8;index++)
+    State->Late.DensityGain = 0.0f;
+    State->Late.ApFeedCoeff = 0.0f;
+    State->Late.MixCoeff = 0.0f;
+
+    for(index = 0;index < 4;index++)
+    {
+        State->Late.ApCoeff[index] = 0.0f;
+        State->Late.ApDelay[index].Mask = length[5 + index] - 1;
+        State->Late.ApDelay[index].Line = &State->SampleBuffer[totalLength];
+        totalLength += length[5 + index];
+
+        // The late all-pass lines have their read offsets calculated once.
+        State->Late.ApOffset[index] = (ALuint)(ALLPASS_LINE_LENGTH[index] *
+                                               Context->Frequency);
+    }
+
+    for(index = 0;index < 4;index++)
     {
         State->Late.Coeff[index] = 0.0f;
-        State->Late.Delay[index].Mask = length[5 + index] - 1;
+        State->Late.Delay[index].Mask = length[9 + index] - 1;
         State->Late.Delay[index].Line = &State->SampleBuffer[totalLength];
-        totalLength += length[5 + index];
+        totalLength += length[9 + index];
 
         State->Late.Offset[index] = 0;
-        if(index < 4)
-        {
-            State->Late.Offset[index] = (ALuint)(LATE_LINE_LENGTH[index] * Context->Frequency);
-            State->Late.LpCoeff[index] = 0.0f;
-            State->Late.LpSample[index] = 0.0f;
-        }
-        else if(index == 4)
-        {
-            State->Late.LpCoeff[index] = 0.0f;
-            State->Late.LpSample[index] = 0.0f;
-        }
+
+        State->Late.LpCoeff[index][0] = 0.0f;
+        State->Late.LpCoeff[index][1] = 0.0f;
+        State->Late.LpSample[index] = 0.0f;
     }
 
     State->Offset = 0;
@@ -402,24 +449,37 @@ ALvoid VerbDestroy(ALverbState *State)
 ALvoid VerbUpdate(ALCcontext *Context, ALeffectslot *Slot, ALeffect *Effect)
 {
     ALverbState *State = Slot->ReverbState;
-    ALuint index, index2;
-    ALfloat length, lpcoeff, cw, g;
+    ALuint index;
+    ALfloat length, mixCoeff, cw, g, lpCoeff;
     ALfloat hfRatio = Effect->Reverb.DecayHFRatio;
 
-    // Calculate the master gain (from the slot and master reverb gain).
+    // Calculate the master gain (from the slot and master effect gain).
     State->Gain = Slot->Gain * Effect->Reverb.Gain;
 
     // Calculate the initial delay taps.
     length = Effect->Reverb.ReflectionsDelay;
     State->Tap[0] = (ALuint)(length * Context->Frequency);
+
     length += Effect->Reverb.LateReverbDelay;
-    State->Tap[1] = (ALuint)(length * Context->Frequency);
 
-    // Calculate the early reflections gain.  Right now this uses a gain of
-    // 0.75 to compensate for the increase in density.  It should probably
-    // use a power (RMS) based measurement from the resulting distribution of
-    // early delay lines.
-    State->Early.Gain = Effect->Reverb.ReflectionsGain * 0.75f;
+    /* The four inputs to the late reverb are decorrelated to smooth the
+     * initial reverb and reduce harsh echos.  The timings are calculated as
+     * multiples of a fraction of the smallest cyclical delay time. This
+     * result is then adjusted so that the first tap occurs immediately (all
+     * taps are reduced by the shortest fraction).
+     *
+     * offset[index] = ((FRACTION MULTIPLIER^index) - 1) delay
+     */
+    for(index = 0;index < 4;index++)
+    {
+        length += LATE_LINE_LENGTH[0] *
+            (1.0f + (Effect->Reverb.Density * LATE_LINE_MULTIPLIER)) *
+            (DECO_FRACTION * (pow(DECO_MULTIPLIER, (ALfloat)index) - 1.0f));
+        State->Tap[1 + index] = (ALuint)(length * Context->Frequency);
+    }
+
+    // Set the early reflections gain.
+    State->Early.Gain = Effect->Reverb.ReflectionsGain;
 
     // Calculate the gain (coefficient) for each early delay line.
     for(index = 0;index < 4;index++)
@@ -427,29 +487,68 @@ ALvoid VerbUpdate(ALCcontext *Context, ALeffectslot *Slot, ALeffect *Effect)
                                                Effect->Reverb.LateReverbDelay *
                                                -60.0f / 20.0f);
 
-    // Calculate the late reverb gain, adjusted by density, diffusion, and
-    // decay time.  To be accurate, the adjustments should probably use power
-    // measurements for each contribution, but they are not too bad as they
-    // are.
-    State->Late.Gain = Effect->Reverb.LateReverbGain *
-                       (0.45f + (0.55f * Effect->Reverb.Density)) *
-                       (1.0f - (0.25f * Effect->Reverb.Diffusion)) *
-                       (1.0f - (0.025f * Effect->Reverb.DecayTime));
-    State->Late.Diffusion = Effect->Reverb.Diffusion;
-
-    // The EFX specification does not make it clear whether the air
-    // absorption parameter should always take effect.  Both Generic Software
-    // and Generic Hardware only apply it when HF limit is flagged, so that's
-    // what is done here.
+    // Calculate the first mixing matrix coefficient (x).
+    mixCoeff = 1.0f - (0.5f * pow(Effect->Reverb.Diffusion, 3.0f));
+
+    // Set the late reverb gain.  Since the output is tapped prior to the
+    // application of the delay line coefficients, this gain needs to be
+    // attenuated by the mix coefficient from above.
+    State->Late.Gain = Effect->Reverb.LateReverbGain * mixCoeff;
+
+    /* 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 is calculated as the ratio between a
+     * reference value and the current approximation of energy for the output
+     * signal.
+     *
+     * Reverb output matches exponential decay of the form Sum(a^n), where a
+     * is the attenuation coefficient, and n is the sample ranging from 0 to
+     * infinity.  The signal energy can thus be approximated using the area
+     * under this curve, calculated as:  1 / (1 - a).
+     *
+     * The reference energy is calculated from a signal at the lowest (effect
+     * at 1.0) density with a decay time of one second.
+     *
+     * The coefficient is calculated as the average length of the cyclical
+     * delay lines.  This produces a better result than calculating the gain
+     * for each line individually (most likely a side effect of diffusion).
+     *
+     * The final result is the square root of the ratio bound to a maximum
+     * value of 1 (no amplification) and attenuated by 1 / sqrt(2) to
+     * compensate for the four decorrelated inputs.
+     */
+    length = (LATE_LINE_LENGTH[0] + LATE_LINE_LENGTH[1] +
+              LATE_LINE_LENGTH[2] + LATE_LINE_LENGTH[3]);
+    g = length * (1.0f + LATE_LINE_MULTIPLIER) * 0.25f;
+    g = pow(10.0f, g * -60.0f / 20.0f);
+    g = 1.0f / (1.0f - g);
+    length *= 1.0f + (Effect->Reverb.Density * LATE_LINE_MULTIPLIER) * 0.25f;
+    length = pow(10.0f, length / Effect->Reverb.DecayTime * -60.0f / 20.0f);
+    length = 1.0f / (1.0f - length);
+    State->Late.DensityGain = 0.707106f * __min(aluSqrt(g / length), 1.0f);
+
+    // Calculate the all-pass feed-back and feed-forward coefficient.
+    State->Late.ApFeedCoeff = 0.6f * pow(Effect->Reverb.Diffusion, 3.0f);
+
+    // Calculate the mixing matrix coefficient (y / x).
+    g = aluSqrt((1.0f - (mixCoeff * mixCoeff)) / 3.0f);
+    State->Late.MixCoeff = g / mixCoeff;
+
+    for(index = 0;index < 4;index++)
+    {
+        // Calculate the gain (coefficient) for each all-pass line.
+        State->Late.ApCoeff[index] = pow(10.0f, ALLPASS_LINE_LENGTH[index] /
+                                                Effect->Reverb.DecayTime *
+                                                -60.0f / 20.0f);
+    }
+
     // If the HF limit parameter is flagged, calculate an appropriate limit
     // based on the air absorption parameter.
     if(Effect->Reverb.DecayHFLimit && Effect->Reverb.AirAbsorptionGainHF < 1.0f)
     {
         ALfloat limitRatio;
 
-        // The following is my best guess at how to limit the HF ratio by the
-        // air absorption parameter.
-        // For each of the last 4 delays, find the attenuation due to air
+        // For each of the cyclical delays, 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
@@ -466,59 +565,58 @@ ALvoid VerbUpdate(ALCcontext *Context, ALeffectslot *Slot, ALeffect *Effect)
         hfRatio = __min(hfRatio, limitRatio);
     }
 
-    cw = cos(2.0f*3.141592654f * LOWPASSFREQCUTOFF / Context->Frequency);
+    // Calculate the filter frequency for low-pass or high-pass depending on
+    // whether the HF ratio is above 1.
+    cw = 2.0f * M_PI * LOWPASSFREQCUTOFF / Context->Frequency;
+    if(hfRatio > 1.0f)
+        cw = M_PI - cw;
+    cw = cos(cw);
 
-    for(index = 0;index < 8;index++)
+    for(index = 0;index < 4;index++)
     {
-        // Calculate the length (in seconds) of each delay line.
-        length = LATE_LINE_LENGTH[index];
-        if(index >= 4)
-        {
-            // Calculate the delay offset for the variable-length delay
-            // lines.
-            length *= 1.0f + (Effect->Reverb.Density * LATE_LINE_MULTIPLIER);
-            State->Late.Offset[index] = (ALuint)(length * Context->Frequency);
-        }
-        // Calculate the gain (coefficient) for each line.
+        // Calculate the length (in seconds) of each cyclical delay line.
+        length = LATE_LINE_LENGTH[index] * (1.0f + (Effect->Reverb.Density *
+                                                    LATE_LINE_MULTIPLIER));
+        // Calculate the delay offset for the cyclical delay lines.
+        State->Late.Offset[index] = (ALuint)(length * Context->Frequency);
+
+        // Calculate the gain (coefficient) for each cyclical line.
         State->Late.Coeff[index] = pow(10.0f, length / Effect->Reverb.DecayTime *
                                               -60.0f / 20.0f);
-        if(index >= 4)
-        {
-            index2 = index - 3;
-
-            // Calculate the decay equation for each low-pass filter.
-            g = pow(10.0f, length / (Effect->Reverb.DecayTime * hfRatio) *
-                           -60.0f / 20.0f) /
-                State->Late.Coeff[index];
-            g  = __max(g, 0.1f);
-            g *= g;
-            // Calculate the gain (coefficient) for each low-pass filter.
-            lpcoeff = 0.0f;
-            if(g < 0.9999f) // 1-epsilon
-                lpcoeff = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
-
-            // Very low decay times will produce minimal output, so apply an
-            // upper bound to the coefficient.
-            State->Late.LpCoeff[index2] = __min(lpcoeff, 0.98f);
-        }
-    }
-
-    // This just calculates the coefficient for the late reverb input low-
-    // pass filter.  It is calculated based the average (hence -30 instead
-    // of -60) length of the inner two variable-length delay lines.
-    length = LATE_LINE_LENGTH[5] * (1.0f + Effect->Reverb.Density * LATE_LINE_MULTIPLIER) +
-             LATE_LINE_LENGTH[6] * (1.0f + Effect->Reverb.Density * LATE_LINE_MULTIPLIER);
-
-    g = pow(10.0f, ((length / (Effect->Reverb.DecayTime * hfRatio))-
-                    (length / Effect->Reverb.DecayTime)) * -30.0f / 20.0f);
-    g  = __max(g, 0.1f);
-    g *= g;
 
-    lpcoeff = 0.0f;
-    if(g < 0.9999f) // 1-epsilon
-        lpcoeff = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
+        // Calculate the decay equation for each low-pass filter.
+        g = pow(10.0f, length / (Effect->Reverb.DecayTime * hfRatio) *
+                       -60.0f / 20.0f);
+        if (hfRatio > 1.0f)
+            g = State->Late.Coeff[index] / g;
+        else
+            g = g / State->Late.Coeff[index];
+        g  = __max(g, 0.1f);
+        g *= g;
+
+        // Calculate the gain (coefficient) for each low-pass filter.
+        lpCoeff = 0.0f;
+        if(g < 0.9999f) // 1-epsilon
+            lpCoeff = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
+
+        // Very low decay times will produce minimal output, so apply an
+        // upper bound to the coefficient.
+        lpCoeff = __min(lpCoeff, 0.98f);
+
+        // Calculate the filter coefficients for high-pass or low-pass
+        // dependent on HF ratio being above 1.
+        if(hfRatio > 1.0f) {
+            State->Late.LpCoeff[index][0] = 1.0f + lpCoeff;
+            State->Late.LpCoeff[index][1] = -lpCoeff;
+        } else {
+            State->Late.LpCoeff[index][0] = 1.0f - lpCoeff;
+            State->Late.LpCoeff[index][1] = lpCoeff;
+        }
 
-    State->Late.LpCoeff[0] = __min(lpcoeff, 0.98f);
+        // Attenuate the cyclical line coefficients by the mixing coefficient
+        // (x).
+        State->Late.Coeff[index] *= mixCoeff;
+    }
 }
 
 // This processes the reverb state, given the input samples and an output
@@ -526,7 +624,7 @@ ALvoid VerbUpdate(ALCcontext *Context, ALeffectslot *Slot, ALeffect *Effect)
 ALvoid VerbProcess(ALverbState *State, ALuint SamplesToDo, const ALfloat *SamplesIn, ALfloat (*SamplesOut)[OUTPUTCHANNELS])
 {
     ALuint index;
-    ALfloat in, early[2], late[2], out[2];
+    ALfloat in[4], early[2], late[2], out[2];
 
     for(index = 0;index < SamplesToDo;index++)
     {
@@ -534,11 +632,14 @@ ALvoid VerbProcess(ALverbState *State, ALuint SamplesToDo, const ALfloat *Sample
         DelayLineIn(&State->Delay, State->Offset, SamplesIn[index]);
 
         // Calculate the early reflection from the first delay tap.
-        in = DelayLineOut(&State->Delay, State->Offset - State->Tap[0]);
-        EarlyReflection(State, in, early);
-
-        // Calculate the late reverb from the second delay tap.
-        in = DelayLineOut(&State->Delay, State->Offset - State->Tap[1]);
+        in[0] = DelayLineOut(&State->Delay, State->Offset - State->Tap[0]);
+        EarlyReflection(State, in[0], early);
+
+        // Calculate the late reverb from the last four delay taps.
+        in[0] = DelayLineOut(&State->Delay, State->Offset - State->Tap[1]);
+        in[1] = DelayLineOut(&State->Delay, State->Offset - State->Tap[2]);
+        in[2] = DelayLineOut(&State->Delay, State->Offset - State->Tap[3]);
+        in[3] = DelayLineOut(&State->Delay, State->Offset - State->Tap[4]);
         LateReverb(State, in, late);
 
         // Mix early reflections and late reverb.