Reorganize and improve the reverb effect

Code supplied by Christopher Fitzgerald. This update also implements the echo
and modulation parameters.

1 files changed, 895 insertions, 354 deletions

diff --git a/Alc/alcReverb.c b/Alc/alcReverb.c
index 8a5a1149..631d1b94 100644
--- a/Alc/alcReverb.c
+++ b/Alc/alcReverb.c
@@ -20,8 +20,9 @@
 
 #include "config.h"
 
-#include <math.h>
+#include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 
 #include "AL/al.h"
 #include "AL/alc.h"
@@ -33,8 +34,8 @@
 
 typedef struct DelayLine
 {
-    // The delay lines use sample 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 the
+    // use of bit-masking instead of a modulus for wrapping.
     ALuint   Mask;
     ALfloat *Line;
 } DelayLine;
@@ -46,29 +47,48 @@ typedef struct ALverbState {
     // All delay lines are allocated as a single buffer to reduce memory
     // fragmentation and management code.
     ALfloat  *SampleBuffer;
-    ALuint    TotalLength;
+    ALuint    TotalSamples;
     // Master effect low-pass filter (2 chained 1-pole filters).
     FILTER    LpFilter;
     ALfloat   LpHistory[2];
-    // Initial effect delay and decorrelation.
+    struct {
+        // Modulator delay line.
+        DelayLine Delay;
+        // The vibrato time is tracked with an index over a modulus-wrapped
+        // range (in samples).
+        ALuint    Index;
+        ALuint    Range;
+        // The depth of frequency change (also in samples) and its filter.
+        ALfloat   Depth;
+        ALfloat   Coeff;
+        ALfloat   Filter;
+    } Mod;
+    // Initial effect delay.
     DelayLine Delay;
     // The tap points for the initial delay.  First tap goes to early
-    // reflections, the last four decorrelate to late reverb.
-    ALuint    Tap[5];
+    // reflections, the last to late reverb.
+    ALuint    DelayTap[2];
     struct {
-        // Total gain for early reflections.
+        // Output gain for early reflections.
         ALfloat   Gain;
         // Early reflections are done with 4 delay lines.
         ALfloat   Coeff[4];
         DelayLine Delay[4];
         ALuint    Offset[4];
-        // The gain for each output channel based on 3D panning.
+        // The gain for each output channel based on 3D panning (only for the
+        // EAX path).
         ALfloat   PanGain[OUTPUTCHANNELS];
     } Early;
+    // Decorrelator delay line.
+    DelayLine Decorrelator;
+    // There are actually 4 decorrelator taps, but the first occurs at the
+    // initial sample.
+    ALuint    DecoTap[3];
     struct {
-        // Total gain for late reverb.
+        // Output gain for late reverb.
         ALfloat   Gain;
-        // Attenuation to compensate for modal density and decay rate.
+        // Attenuation to compensate for the modal density and decay rate of
+        // the late lines.
         ALfloat   DensityGain;
         // The feed-back and feed-forward all-pass coefficient.
         ALfloat   ApFeedCoeff;
@@ -85,13 +105,58 @@ typedef struct ALverbState {
         // The cyclical delay lines are 1-pole low-pass filtered.
         ALfloat   LpCoeff[4];
         ALfloat   LpSample[4];
-        // The gain for each output channel based on 3D panning.
+        // The gain for each output channel based on 3D panning (only for the
+        // EAX path).
         ALfloat   PanGain[OUTPUTCHANNELS];
     } Late;
+    struct {
+        // Attenuation to compensate for the modal density and decay rate of
+        // the echo line.
+        ALfloat   DensityGain;
+        // Echo delay and all-pass lines.
+        DelayLine Delay;
+        DelayLine ApDelay;
+        ALfloat   Coeff;
+        ALfloat   ApFeedCoeff;
+        ALfloat   ApCoeff;
+        ALuint    Offset;
+        ALuint    ApOffset;
+        // The echo line is 1-pole low-pass filtered.
+        ALfloat   LpCoeff;
+        ALfloat   LpSample;
+        // Echo mixing coefficients.
+        ALfloat   MixCoeff[2];
+    } Echo;
     // The current read offset for all delay lines.
     ALuint Offset;
 } ALverbState;
 
+/* This coefficient is used to define the maximum frequency range controlled
+ * by the modulation depth.  The current value of 0.1 will allow it to swing
+ * from 0.9x to 1.1x.  This value must be below 1.  At 1 it will cause the
+ * sampler to stall on the downswing, and above 1 it will cause it to sample
+ * backwards.
+ */
+static const ALfloat MODULATION_DEPTH_COEFF = 0.1f;
+
+/* A filter is used to avoid the terrible distortion caused by changing
+ * modulation time and/or depth.  To be consistent across different sample
+ * rates, the coefficient must be raised to a constant divided by the sample
+ * rate:  coeff^(constant / rate).
+ */
+static const ALfloat MODULATION_FILTER_COEFF = 0.048f;
+static const ALfloat MODULATION_FILTER_CONST = 100000.0f;
+
+// When diffusion is above 0, an all-pass filter is used to take the edge off
+// the echo effect.  It uses the following line length (in seconds).
+static const ALfloat ECHO_ALLPASS_LENGTH = 0.0133f;
+
+// Input into the late reverb is decorrelated between four channels.  Their
+// timings are dependent on a fraction and multiplier.  See the
+// UpdateDecorrelator() routine for the calculations involved.
+static const ALfloat DECO_FRACTION = 0.15f;
+static const ALfloat DECO_MULTIPLIER = 2.0f;
+
 // All delay line lengths are specified in seconds.
 
 // The lengths of the early delay lines.
@@ -116,54 +181,479 @@ static const ALfloat LATE_LINE_LENGTH[4] =
 // 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;
-
-static ALuint CalcLengths(ALuint length[13], ALuint frequency)
+// Calculate the length of a delay line and store its mask and offset.
+static ALuint CalcLineLength(ALfloat length, ALuint offset, ALuint frequency, DelayLine *Delay)
 {
-    ALuint samples, totalLength, index;
+    ALuint samples;
 
     // All line lengths are powers of 2, calculated from their lengths, with
     // an additional sample in case of rounding errors.
+    samples = NextPowerOf2((ALuint)(length * frequency) + 1);
+    // All lines share a single sample buffer.
+    Delay->Mask = samples - 1;
+    Delay->Line = (ALfloat*)offset;
+    // Return the sample count for accumulation.
+    return samples;
+}
+
+// Given the allocated sample buffer, this function updates each delay line
+// offset.
+static __inline ALvoid RealizeLineOffset(ALfloat * sampleBuffer, DelayLine *Delay)
+{
+    Delay->Line = &sampleBuffer[(ALuint)Delay->Line];
+}
+
+/* Calculates the delay line metrics and allocates the shared sample buffer
+ * for all lines given a flag indicating whether or not to allocate the EAX-
+ * related delays (eaxFlag) and the sample rate (frequency).  If an
+ * allocation failure occurs, it returns AL_FALSE.
+ */
+static ALboolean AllocLines(ALboolean eaxFlag, ALuint frequency, ALverbState *State)
+{
+    ALuint totalSamples, index;
+    ALfloat length;
+    ALfloat *newBuffer = NULL;
+
+    // All delay line lengths are calculated to accomodate the full range of
+    // lengths given their respective paramters.
+    totalSamples = 0;
+    if(eaxFlag)
+    {
+        /* The modulator's line length is calculated from the maximum
+         * modulation time and depth coefficient, and halfed for the low-to-
+         * high frequency swing.  An additional sample is added to keep it
+         * stable when there is no modulation.
+         */
+        length = (AL_EAXREVERB_MAX_MODULATION_TIME * MODULATION_DEPTH_COEFF /
+                  2.0f) + (1.0f / frequency);
+        totalSamples += CalcLineLength(length, totalSamples, frequency,
+                                       &State->Mod.Delay);
+    }
+
+    // The initial delay is the sum of the reflections and late reverb
+    // delays.
+    if(eaxFlag)
+        length = AL_EAXREVERB_MAX_REFLECTIONS_DELAY +
+                 AL_EAXREVERB_MAX_LATE_REVERB_DELAY;
+    else
+        length = AL_REVERB_MAX_REFLECTIONS_DELAY +
+                 AL_REVERB_MAX_LATE_REVERB_DELAY;
+    totalSamples += CalcLineLength(length, totalSamples, frequency,
+                                   &State->Delay);
+
+    // The early reflection lines.
+    for(index = 0;index < 4;index++)
+        totalSamples += CalcLineLength(EARLY_LINE_LENGTH[index], totalSamples,
+                                       frequency, &State->Early.Delay[index]);
+
+    // The decorrelator line is calculated from the lowest reverb density (a
+    // parameter value of 1).
+    length = (DECO_FRACTION * DECO_MULTIPLIER * DECO_MULTIPLIER) *
+             LATE_LINE_LENGTH[0] * (1.0f + LATE_LINE_MULTIPLIER);
+    totalSamples += CalcLineLength(length, totalSamples, frequency,
+                                   &State->Decorrelator);
+
+    // The late all-pass lines.
+    for(index = 0;index < 4;index++)
+        totalSamples += CalcLineLength(ALLPASS_LINE_LENGTH[index], totalSamples,
+                                       frequency, &State->Late.ApDelay[index]);
+
+    // 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,
+                                       &State->Late.Delay[index]);
+    }
 
-    // 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)))) *
-               frequency) + 1;
-    length[0] = NextPowerOf2(samples);
-    totalLength = length[0];
+    if(eaxFlag)
+    {
+        // The echo all-pass and delay lines.
+        totalSamples += CalcLineLength(ECHO_ALLPASS_LENGTH, totalSamples,
+                                       frequency, &State->Echo.ApDelay);
+        totalSamples += CalcLineLength(AL_EAXREVERB_MAX_ECHO_TIME, totalSamples,
+                                       frequency, &State->Echo.Delay);
+    }
+
+    if(totalSamples != State->TotalSamples)
+    {
+        newBuffer = realloc(State->SampleBuffer, sizeof(ALfloat) * totalSamples);
+        if(newBuffer == NULL)
+            return AL_FALSE;
+        State->SampleBuffer = newBuffer;
+        State->TotalSamples = totalSamples;
+    }
+
+    // Update all delays to reflect the new sample buffer.
+    RealizeLineOffset(State->SampleBuffer, &State->Delay);
+    RealizeLineOffset(State->SampleBuffer, &State->Decorrelator);
     for(index = 0;index < 4;index++)
     {
-        samples = (ALuint)(EARLY_LINE_LENGTH[index] * frequency) + 1;
-        length[1 + index] = NextPowerOf2(samples);
-        totalLength += length[1 + index];
+        RealizeLineOffset(State->SampleBuffer, &State->Early.Delay[index]);
+        RealizeLineOffset(State->SampleBuffer, &State->Late.ApDelay[index]);
+        RealizeLineOffset(State->SampleBuffer, &State->Late.Delay[index]);
+    }
+    if(eaxFlag)
+    {
+        RealizeLineOffset(State->SampleBuffer, &State->Mod.Delay);
+        RealizeLineOffset(State->SampleBuffer, &State->Echo.ApDelay);
+        RealizeLineOffset(State->SampleBuffer, &State->Echo.Delay);
+    }
+
+    // Clear the sample buffer.
+    for(index = 0;index < State->TotalSamples;index++)
+        State->SampleBuffer[index] = 0.0f;
+
+    return AL_TRUE;
+}
+
+// Calculate a decay coefficient given the length of each cycle and the time
+// until the decay reaches -60 dB.
+static __inline ALfloat CalcDecayCoeff(ALfloat length, ALfloat decayTime)
+{
+    return pow(10.0f, length / decayTime * -60.0f / 20.0f);
+}
+
+// Calculate a decay length from a coefficient and the time until the decay
+// reaches -60 dB.
+static __inline ALfloat CalcDecayLength(ALfloat coeff, ALfloat decayTime)
+{
+    return log10(coeff) / -60.0 * 20.0f * decayTime;
+}
+
+// Calculate the high frequency parameter for the I3DL2 coefficient
+// calculation.
+static __inline ALfloat CalcI3DL2HFreq(ALfloat hfRef, ALuint frequency)
+{
+    return cos(2.0f * M_PI * hfRef / frequency);
+}
+
+/* Calculate the I3DL2 coefficient given the gain and frequency parameters.
+ * To allow for optimization when using multiple chained filters, the gain
+ * is not squared in this function.  Callers using a single filter should
+ * square it to produce the correct coefficient.  Those using multiple
+ * filters should find its N-1 root (where N is the number of chained
+ * filters).
+ */
+static __inline ALfloat CalcI3DL2Coeff(ALfloat g, ALfloat cw)
+{
+    ALfloat coeff;
+
+    coeff = 0.0f;
+    if(g < 0.9999f) // 1-epsilon
+        coeff = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
+
+    return coeff;
+}
+
+// Calculate an attenuation to be applied to the input of any echo models to
+// compensate for modal density and decay time.
+static __inline ALfloat CalcDensityGain(ALfloat a)
+{
+    /* The energy of a signal can be obtained by finding the area under the
+     * squared signal.  This takes the form of Sum(x_n^2), where x is the
+     * amplitude for the sample n.
+     *
+     * Decaying feedback matches exponential decay of the form Sum(a^n),
+     * where a is the attenuation coefficient, and n is the sample.  The area
+     * under this decay curve can be calculated as:  1 / (1 - a).
+     *
+     * Modifying the above equation to find the squared area under the curve
+     * (for energy) yields:  1 / (1 - a^2).  Input attenuation can then be
+     * calculated by inverting the square root of this approximation,
+     * yielding:  1 / sqrt(1 / (1 - a^2)), simplified to: sqrt(1 - a^2).
+     */
+    return aluSqrt(1.0f - (a * a));
+}
+
+// Calculate the mixing matrix coefficients given a diffusion factor.
+static __inline ALvoid CalcMatrixCoeffs(ALfloat diffusion, ALfloat *x, ALfloat *y)
+{
+    ALfloat n, t;
+
+    // The matrix is of order 4, so n is sqrt (4 - 1).
+    n = aluSqrt(3.0f);
+    t = diffusion * atan(n);
+
+    // Calculate the first mixing matrix coefficient.
+    *x = cos(t);
+    // Calculate the second mixing matrix coefficient.
+    *y = sin(t) / n;
+}
+
+// Calculate the limited HF ratio for use with the late reverb low-pass
+// filters.
+static __inline ALfloat CalcLimitedHfRatio(ALfloat hfRatio, ALfloat airAbsorptionGainHF, ALfloat decayTime)
+{
+    ALfloat 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);
+    // Need to limit the result to a minimum of 0.1, just like the HF ratio
+    // parameter.
+    limitRatio = __max(limitRatio, 0.1f);
+
+    // Using the limit calculated above, apply the upper bound to the HF
+    // ratio.
+    return __min(hfRatio, limitRatio);
+}
+
+// Calculate the coefficient for a HF (and eventually LF) decay damping
+// filter.
+static __inline ALfloat CalcDampingCoeff(ALfloat hfRatio, ALfloat length, ALfloat decayTime, ALfloat decayCoeff, ALfloat cw)
+{
+    ALfloat coeff, g;
+
+    // Eventually this should boost the high frequencies when the ratio
+    // exceeds 1.
+    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;
+        g  = __max(g, 0.1f);
+
+        // Damping is done with a 1-pole filter, so g needs to be squared.
+        g *= g;
+        coeff = CalcI3DL2Coeff(g, cw);
+
+        // Very low decay times will produce minimal output, so apply an
+        // upper bound to the coefficient.
+        coeff = __min(coeff, 0.98f);
+    }
+    return coeff;
+}
+
+// Update the EAX modulation index, range, and depth.  Keep in mind that this
+// kind of vibrato is additive and not multiplicative as one may expect.  The
+// downswing will sound stronger than the upswing.
+static ALvoid UpdateModulator(ALfloat modTime, ALfloat modDepth, ALuint frequency, ALverbState *State)
+{
+    ALfloat length;
+
+    /* Modulation is calculated in two parts.
+     *
+     * The modulation time effects the sinus applied to the change in
+     * frequency.  An index out of the current time range (both in samples)
+     * is incremented each sample.  The range is bound to a reasonable
+     * minimum (1 sample) and when the timing changes, the index is rescaled
+     * to the new range (to keep the sinus consistent).
+     */
+    length = modTime * frequency;
+    if (length >= 1.0f) {
+       State->Mod.Index = (ALuint)(State->Mod.Index * length /
+                                   State->Mod.Range);
+       State->Mod.Range = (ALuint)length;
+    } else {
+       State->Mod.Index = 0;
+       State->Mod.Range = 1;
     }
+
+    /* The modulation depth effects the amount of frequency change over the
+     * range of the sinus.  It needs to be scaled by the modulation time so
+     * that a given depth produces a consistent change in frequency over all
+     * ranges of time.  Since the depth is applied to a sinus value, it needs
+     * to be halfed once for the sinus range and again for the sinus swing
+     * in time (half of it is spent decreasing the frequency, half is spent
+     * increasing it).
+     */
+    State->Mod.Depth = modDepth * MODULATION_DEPTH_COEFF * modTime / 2.0f /
+                       2.0f * frequency;
+}
+
+// Update the offsets for the initial effect delay line.
+static ALvoid UpdateDelayLine(ALfloat earlyDelay, ALfloat lateDelay, ALuint frequency, ALverbState *State)
+{
+    // Calculate the initial delay taps.
+    State->DelayTap[0] = (ALuint)(earlyDelay * frequency);
+    State->DelayTap[1] = (ALuint)((earlyDelay + lateDelay) * frequency);
+}
+
+// Update the early reflections gain and line coefficients.
+static ALvoid UpdateEarlyLines(ALfloat reverbGain, ALfloat earlyGain, ALfloat lateDelay, ALverbState *State)
+{
+    ALuint index;
+
+    // Calculate the early reflections gain (from the master effect gain, and
+    // reflections gain parameters) with a constant attenuation of 0.5.
+    State->Early.Gain = 0.5f * reverbGain * earlyGain;
+
+    // Calculate the gain (coefficient) for each early delay line using the
+    // late delay time.  This expands the early reflections to the start of
+    // the late reverb.
     for(index = 0;index < 4;index++)
+        State->Early.Coeff[index] = CalcDecayCoeff(EARLY_LINE_LENGTH[index],
+                                                   lateDelay);
+}
+
+// Update the offsets for the decorrelator line.
+static ALvoid UpdateDecorrelator(ALfloat density, ALuint frequency, ALverbState *State)
+{
+    ALuint index;
+    ALfloat length;
+
+    /* The late reverb inputs are decorrelated to smooth the reverb tail and
+     * reduce harsh echos.  The first tap occurs immediately, while the
+     * remaining taps are delayed by multiples of a fraction of the smallest
+     * cyclical delay time.
+     *
+     * offset[index] = (FRACTION (MULTIPLIER^index)) smallest_delay
+     */
+    for(index = 0;index < 3;index++)
     {
-        samples = (ALuint)(ALLPASS_LINE_LENGTH[index] * frequency) + 1;
-        length[5 + index] = NextPowerOf2(samples);
-        totalLength += length[5 + index];
+        length = (DECO_FRACTION * pow(DECO_MULTIPLIER, (ALfloat)index)) *
+                 LATE_LINE_LENGTH[0] * (1.0f + (density * LATE_LINE_MULTIPLIER));
+        State->DecoTap[index] = (ALuint)(length * frequency);
     }
+}
+
+// Update the late reverb gains, line lengths, and line coefficients.
+static ALvoid UpdateLateLines(ALfloat reverbGain, ALfloat lateGain, ALfloat xMix, ALfloat density, ALfloat decayTime, ALfloat diffusion, ALfloat hfRatio, ALfloat cw, ALuint frequency, ALverbState *State)
+{
+    ALfloat length;
+    ALuint 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));
+
+    // Calculate the all-pass feed-back and feed-forward coefficient.
+    State->Late.ApFeedCoeff = 0.5f * pow(diffusion, 2.0f);
+
     for(index = 0;index < 4;index++)
     {
-        samples = (ALuint)(LATE_LINE_LENGTH[index] *
-                           (1.0f + LATE_LINE_MULTIPLIER) * frequency) + 1;
-        length[9 + index] = NextPowerOf2(samples);
-        totalLength += length[9 + index];
+        // Calculate the gain (coefficient) for each all-pass line.
+        State->Late.ApCoeff[index] = CalcDecayCoeff(ALLPASS_LINE_LENGTH[index],
+                                                    decayTime);
+
+        // 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] = (ALuint)(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;
+    }
+}
+
+// Update the echo gain, line offset, line coefficients, and mixing
+// coefficients.
+static ALvoid UpdateEchoLine(ALfloat reverbGain, ALfloat lateGain, ALfloat echoTime, ALfloat decayTime, ALfloat diffusion, ALfloat echoDepth, ALfloat hfRatio, ALfloat cw, ALuint frequency, ALverbState *State)
+{
+    // Calculate the energy-based attenuation coefficient for the echo delay
+    // line.
+    State->Echo.DensityGain = CalcDensityGain(CalcDecayCoeff(echoTime,
+                                                             decayTime));
+
+    // Update the offset and coefficient for the echo delay line.
+    State->Echo.Offset = (ALuint)(echoTime * frequency);
+
+    // Calculate the decay coefficient for the echo line.
+    State->Echo.Coeff = CalcDecayCoeff(echoTime, decayTime);
+
+    // Calculate the echo all-pass feed coefficient.
+    State->Echo.ApFeedCoeff = 0.5f * pow(diffusion, 2.0f);
+
+    // Calculate the echo all-pass attenuation coefficient.
+    State->Echo.ApCoeff = CalcDecayCoeff(ECHO_ALLPASS_LENGTH, decayTime);
+
+    // Calculate the damping coefficient for each low-pass filter.
+    State->Echo.LpCoeff = CalcDampingCoeff(hfRatio, echoTime, decayTime,
+                                           State->Echo.Coeff, cw);
+
+    /* Calculate the echo mixing coefficients.  The first is applied to the
+     * echo itself.  The second is used to attenuate the late reverb when
+     * echo depth is high and diffusion is low, so the echo is slightly
+     * stronger than the decorrelated echos in the reverb tail.
+     */
+    State->Echo.MixCoeff[0] = reverbGain * lateGain * echoDepth;
+    State->Echo.MixCoeff[1] = 1.0f - (echoDepth * 0.5f * (1.0f - diffusion));
+}
+
+// Update the early and late 3D panning gains.
+static ALvoid Update3DPanning(const ALfloat *ReflectionsPan, const ALfloat *LateReverbPan, ALfloat *PanningLUT, ALverbState *State)
+{
+    ALfloat length;
+    ALfloat earlyPan[3] = { ReflectionsPan[0], ReflectionsPan[1],
+                            ReflectionsPan[2] };
+    ALfloat latePan[3] = { LateReverbPan[0], LateReverbPan[1],
+                           LateReverbPan[2] };
+    ALint pos;
+    ALfloat *speakerGain, dirGain, ambientGain;
+    ALuint index;
+
+    // Calculate the 3D-panning gains for the early reflections and late
+    // reverb.
+    length = earlyPan[0]*earlyPan[0] + earlyPan[1]*earlyPan[1] + earlyPan[2]*earlyPan[2];
+    if(length > 1.0f)
+    {
+        length = 1.0f / aluSqrt(length);
+        earlyPan[0] *= length;
+        earlyPan[1] *= length;
+        earlyPan[2] *= length;
+    }
+    length = latePan[0]*latePan[0] + latePan[1]*latePan[1] + latePan[2]*latePan[2];
+    if(length > 1.0f)
+    {
+        length = 1.0f / aluSqrt(length);
+        latePan[0] *= length;
+        latePan[1] *= length;
+        latePan[2] *= length;
     }
 
-    return totalLength;
+    /* This code applies directional reverb just like the mixer applies
+     * directional sources.  It diffuses the sound toward all speakers as the
+     * magnitude of the panning vector drops, which is only a rough
+     * approximation of the expansion of sound across the speakers from the
+     * panning direction.
+     */
+    pos = aluCart2LUTpos(earlyPan[2], earlyPan[0]);
+    speakerGain = &PanningLUT[OUTPUTCHANNELS * pos];
+    dirGain = aluSqrt((earlyPan[0] * earlyPan[0]) + (earlyPan[2] * earlyPan[2]));
+    ambientGain = (1.0 - dirGain);
+    for(index = 0;index < OUTPUTCHANNELS;index++)
+         State->Early.PanGain[index] = dirGain * speakerGain[index] + ambientGain;
+
+    pos = aluCart2LUTpos(latePan[2], latePan[0]);
+    speakerGain = &PanningLUT[OUTPUTCHANNELS * pos];
+    dirGain = aluSqrt((latePan[0] * latePan[0]) + (latePan[2] * latePan[2]));
+    ambientGain = (1.0 - dirGain);
+    for(index = 0;index < OUTPUTCHANNELS;index++)
+         State->Late.PanGain[index] = dirGain * speakerGain[index] + ambientGain;
 }
 
 // Basic delay line input/output routines.
@@ -177,16 +667,75 @@ static __inline ALvoid DelayLineIn(DelayLine *Delay, ALuint offset, ALfloat in)
     Delay->Line[offset&Delay->Mask] = in;
 }
 
+// Attenuated delay line output routine.
+static __inline ALfloat AttenuatedDelayLineOut(DelayLine *Delay, ALuint offset, ALfloat coeff)
+{
+    return coeff * Delay->Line[offset&Delay->Mask];
+}
+
+// Basic attenuated all-pass input/output routine.
+static __inline ALfloat AllpassInOut(DelayLine *Delay, ALuint outOffset, ALuint inOffset, ALfloat in, ALfloat feedCoeff, ALfloat coeff)
+{
+    ALfloat out, feed;
+
+    out = DelayLineOut(Delay, outOffset);
+    feed = feedCoeff * in;
+    DelayLineIn(Delay, inOffset, (feedCoeff * (out - feed)) + in);
+
+    // The time-based attenuation is only applied to the delay output to
+    // keep it from affecting the feed-back path (which is already controlled
+    // by the all-pass feed coefficient).
+    return (coeff * out) - feed;
+}
+
+// Given an input sample, this function produces modulation for the late
+// reverb.
+static __inline ALfloat EAXModulation(ALverbState *State, ALfloat in)
+{
+    ALfloat sinus, frac;
+    ALuint offset;
+    ALfloat out0, out1;
+
+    // Calculate the sinus rythm (dependent on modulation time and the
+    // sampling rate).  The center of the sinus is moved to reduce the delay
+    // of the effect when the time or depth are low.
+    sinus = 1.0f + sin(2.0f * M_PI * State->Mod.Index / State->Mod.Range);
+
+    // The depth determines the range over which to read the input samples
+    // from, so it must be filtered to reduce the distortion caused by even
+    // small parameter changes.
+    State->Mod.Filter += (State->Mod.Depth - State->Mod.Filter) *
+                         State->Mod.Coeff;
+
+    // Calculate the read offset and fraction between it and the next sample.
+    frac   = (1.0f + (State->Mod.Filter * sinus));
+    offset = (ALuint)frac;
+    frac  -= offset;
+
+    // Get the two samples crossed by the offset, and feed the delay line
+    // with the next input sample.
+    out0 = DelayLineOut(&State->Mod.Delay, State->Offset - offset);
+    out1 = DelayLineOut(&State->Mod.Delay, State->Offset - offset - 1);
+    DelayLineIn(&State->Mod.Delay, State->Offset, in);
+
+    // Step the modulation index forward, keeping it bound to its range.
+    State->Mod.Index = (State->Mod.Index + 1) % State->Mod.Range;
+
+    // The output is obtained by linearly interpolating the two samples that
+    // were acquired above.
+    return out0 + ((out1 - out0) * frac);
+}
+
 // Delay line output routine for early reflections.
 static __inline ALfloat EarlyDelayLineOut(ALverbState *State, ALuint index)
 {
-    return State->Early.Coeff[index] *
-           DelayLineOut(&State->Early.Delay[index],
-                        State->Offset - State->Early.Offset[index]);
+    return AttenuatedDelayLineOut(&State->Early.Delay[index],
+                                  State->Offset - State->Early.Offset[index],
+                                  State->Early.Coeff[index]);
 }
 
-// Given an input sample, this function produces stereo output for early
-// reflections.
+// Given an input sample, this function produces four-channel output for the
+// early reflections.
 static __inline ALvoid EarlyReflection(ALverbState *State, ALfloat in, ALfloat *out)
 {
     ALfloat d[4], v, f[4];
@@ -200,7 +749,7 @@ 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 feedback delay network (FDN).
+     * be considered a simple feed-back delay network (FDN).
      *          N
      *         ---
      *         \
@@ -218,13 +767,13 @@ static __inline ALvoid EarlyReflection(ALverbState *State, ALfloat in, ALfloat *
     f[2] = v - d[2];
     f[3] = v - d[3];
 
-    // Refeed the delay lines.
+    // Re-feed the delay lines.
     DelayLineIn(&State->Early.Delay[0], State->Offset, f[0]);
     DelayLineIn(&State->Early.Delay[1], State->Offset, f[1]);
     DelayLineIn(&State->Early.Delay[2], State->Offset, f[2]);
     DelayLineIn(&State->Early.Delay[3], State->Offset, f[3]);
 
-    // Output the results of the junction for all four lines.
+    // Output the results of the junction for all four channels.
     out[0] = State->Early.Gain * f[0];
     out[1] = State->Early.Gain * f[1];
     out[2] = State->Early.Gain * f[2];
@@ -234,23 +783,18 @@ static __inline ALvoid EarlyReflection(ALverbState *State, ALfloat in, ALfloat *
 // 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;
+    return AllpassInOut(&State->Late.ApDelay[index],
+                        State->Offset - State->Late.ApOffset[index],
+                        State->Offset, in, State->Late.ApFeedCoeff,
+                        State->Late.ApCoeff[index]);
 }
 
 // Delay line output routine for late reverb.
 static __inline ALfloat LateDelayLineOut(ALverbState *State, ALuint index)
 {
-    return State->Late.Coeff[index] *
-           DelayLineOut(&State->Late.Delay[index],
-                        State->Offset - State->Late.Offset[index]);
+    return AttenuatedDelayLineOut(&State->Late.Delay[index],
+                                  State->Offset - State->Late.Offset[index],
+                                  State->Late.Coeff[index]);
 }
 
 // Low-pass filter input/output routine for late reverb.
@@ -261,33 +805,32 @@ static __inline ALfloat LateLowPassInOut(ALverbState *State, ALuint index, ALflo
     return State->Late.LpSample[index];
 }
 
-// Given four decorrelated input samples, this function produces stereo
-// output for late reverb.
+// Given four decorrelated input samples, this function produces four-channel
+// output for the late reverb.
 static __inline ALvoid LateReverb(ALverbState *State, ALfloat *in, ALfloat *out)
 {
     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));
+    // 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(State, 2, in[2] + LateDelayLineOut(State, 2));
+    d[1] = LateLowPassInOut(State, 0, in[0] + LateDelayLineOut(State, 0));
+    d[2] = LateLowPassInOut(State, 3, in[3] + LateDelayLineOut(State, 3));
+    d[3] = LateLowPassInOut(State, 1, in[1] + LateDelayLineOut(State, 1));
 
     // 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)
+    // When there is no diffusion, the shortest all-pass filter will feed the
+    // shortest delay line.
+    d[0] = LateAllPassInOut(State, 0, d[0]);
+    d[1] = LateAllPassInOut(State, 1, d[1]);
+    d[2] = LateAllPassInOut(State, 2, d[2]);
+    d[3] = LateAllPassInOut(State, 3, d[3]);
+
+    /* 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
@@ -305,10 +848,10 @@ static __inline ALvoid LateReverb(ALverbState *State, ALfloat *in, ALfloat *out)
      * 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 ]
+     *  [  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
@@ -319,52 +862,129 @@ static __inline ALvoid LateReverb(ALverbState *State, ALfloat *in, ALfloat *out)
     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 the results of the matrix for all four cyclical delay lines,
-    // attenuated by the late reverb gain (which is attenuated by the 'x'
-    // mix coefficient).
+    // 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] = State->Late.Gain * f[0];
     out[1] = State->Late.Gain * f[1];
     out[2] = State->Late.Gain * f[2];
     out[3] = State->Late.Gain * f[3];
 
-    // 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]);
+    // Re-feed the cyclical delay lines.
+    DelayLineIn(&State->Late.Delay[0], State->Offset, f[0]);
+    DelayLineIn(&State->Late.Delay[1], State->Offset, f[1]);
+    DelayLineIn(&State->Late.Delay[2], State->Offset, f[2]);
+    DelayLineIn(&State->Late.Delay[3], State->Offset, f[3]);
+}
+
+// Given an input sample, this function mixes echo into the four-channel late
+// reverb.
+static __inline ALvoid EAXEcho(ALverbState *State, ALfloat in, ALfloat *late)
+{
+    ALfloat out, feed;
+
+    // Get the latest attenuated echo sample for output.
+    feed = AttenuatedDelayLineOut(&State->Echo.Delay,
+                                  State->Offset - State->Echo.Offset,
+                                  State->Echo.Coeff);
+
+    // Mix the output into the late reverb channels.
+    out = State->Echo.MixCoeff[0] * feed;
+    late[0] = (State->Echo.MixCoeff[1] * late[0]) + out;
+    late[1] = (State->Echo.MixCoeff[1] * late[1]) + out;
+    late[2] = (State->Echo.MixCoeff[1] * late[2]) + out;
+    late[3] = (State->Echo.MixCoeff[1] * late[3]) + out;
+
+    // Mix the energy-attenuated input with the output and pass it through
+    // the echo low-pass filter.
+    feed += State->Echo.DensityGain * in;
+    feed += ((State->Echo.LpSample - feed) * State->Echo.LpCoeff);
+    State->Echo.LpSample = feed;
+
+    // Then the echo all-pass filter.
+    feed = AllpassInOut(&State->Echo.ApDelay,
+                       State->Offset - State->Echo.ApOffset,
+                       State->Offset, feed, State->Echo.ApFeedCoeff,
+                       State->Echo.ApCoeff);
+
+    // Feed the delay with the mixed and filtered sample.
+    DelayLineIn(&State->Echo.Delay, State->Offset, feed);
+}
+
+// Perform the non-EAX reverb pass on a given input sample, resulting in
+// four-channel output.
+static __inline ALvoid VerbPass(ALverbState *State, ALfloat in, ALfloat *early, ALfloat *late)
+{
+    ALfloat taps[4];
+
+    // Low-pass filter the incoming sample.
+    in = lpFilter2P(&State->LpFilter, 0, in);
+
+    // Feed the initial delay line.
+    DelayLineIn(&State->Delay, State->Offset, in);
+
+    // Calculate the early reflection from the first delay tap.
+    in = DelayLineOut(&State->Delay, State->Offset - State->DelayTap[0]);
+    EarlyReflection(State, in, early);
+
+    // Feed the decorrelator from the energy-attenuated output of the second
+    // delay tap.
+    in = DelayLineOut(&State->Delay, State->Offset - State->DelayTap[1]);
+    in *= State->Late.DensityGain;
+    DelayLineIn(&State->Decorrelator, State->Offset, in);
+
+    // Calculate the late reverb from the decorrelator taps.
+    taps[0] = in;
+    taps[1] = DelayLineOut(&State->Decorrelator, State->Offset - State->DecoTap[0]);
+    taps[2] = DelayLineOut(&State->Decorrelator, State->Offset - State->DecoTap[1]);
+    taps[3] = DelayLineOut(&State->Decorrelator, State->Offset - State->DecoTap[2]);
+    LateReverb(State, taps, late);
+
+    // Step all delays forward one sample.
+    State->Offset++;
 }
 
-// Process the reverb for a given input sample, resulting in separate four-
-// channel output for both early reflections and late reverb.
-static __inline ALvoid ReverbInOut(ALverbState *State, ALfloat in, ALfloat *early, ALfloat *late)
+// Perform the EAX reverb pass on a given input sample, resulting in four-
+// channel output.
+static __inline ALvoid EAXVerbPass(ALverbState *State, ALfloat in, ALfloat *early, ALfloat *late)
 {
     ALfloat taps[4];
 
     // Low-pass filter the incoming sample.
     in = lpFilter2P(&State->LpFilter, 0, in);
 
+    // Perform any modulation on the input.
+    in = EAXModulation(State, in);
+
     // Feed the initial delay line.
     DelayLineIn(&State->Delay, State->Offset, in);
 
     // Calculate the early reflection from the first delay tap.
-    in = DelayLineOut(&State->Delay, State->Offset - State->Tap[0]);
+    in = DelayLineOut(&State->Delay, State->Offset - State->DelayTap[0]);
     EarlyReflection(State, in, early);
 
-    // Calculate the late reverb from the last four delay taps.
-    taps[0] = DelayLineOut(&State->Delay, State->Offset - State->Tap[1]);
-    taps[1] = DelayLineOut(&State->Delay, State->Offset - State->Tap[2]);
-    taps[2] = DelayLineOut(&State->Delay, State->Offset - State->Tap[3]);
-    taps[3] = DelayLineOut(&State->Delay, State->Offset - State->Tap[4]);
+    // Feed the decorrelator from the energy-attenuated output of the second
+    // delay tap.
+    in = DelayLineOut(&State->Delay, State->Offset - State->DelayTap[1]);
+    in *= State->Late.DensityGain;
+    DelayLineIn(&State->Decorrelator, State->Offset, in);
+
+    // Calculate the late reverb from the decorrelator taps.
+    taps[0] = in;
+    taps[1] = DelayLineOut(&State->Decorrelator, State->Offset - State->DecoTap[0]);
+    taps[2] = DelayLineOut(&State->Decorrelator, State->Offset - State->DecoTap[1]);
+    taps[3] = DelayLineOut(&State->Decorrelator, State->Offset - State->DecoTap[2]);
     LateReverb(State, taps, late);
 
+    // Calculate and mix in any echo.
+    EAXEcho(State, in, late);
+
     // Step all delays forward one sample.
     State->Offset++;
 }
 
 // This destroys the reverb state.  It should be called only when the effect
 // slot has a different (or no) effect loaded over the reverb effect.
-ALvoid VerbDestroy(ALeffectState *effect)
+static ALvoid VerbDestroy(ALeffectState *effect)
 {
     ALverbState *State = (ALverbState*)effect;
     if(State)
@@ -377,283 +997,174 @@ ALvoid VerbDestroy(ALeffectState *effect)
 
 // This updates the device-dependant reverb state.  This is called on
 // initialization and any time the device parameters (eg. playback frequency,
-// format) have been changed.
-ALboolean VerbDeviceUpdate(ALeffectState *effect, ALCdevice *Device)
+// or format) have been changed.
+static ALboolean VerbDeviceUpdate(ALeffectState *effect, ALCdevice *Device)
 {
     ALverbState *State = (ALverbState*)effect;
-    ALuint length[13], totalLength;
-    ALuint index;
+    ALuint frequency = Device->Frequency, index;
 
-    totalLength = CalcLengths(length, Device->Frequency);
-    if(totalLength != State->TotalLength)
+    // Allocate the delay lines.
+    if(!AllocLines(AL_FALSE, frequency, State))
     {
-        void *temp;
-
-        temp = realloc(State->SampleBuffer, totalLength * sizeof(ALfloat));
-        if(!temp)
-        {
-            alSetError(AL_OUT_OF_MEMORY);
-            return AL_FALSE;
-        }
-        State->TotalLength = totalLength;
-        State->SampleBuffer = temp;
-
-        // All lines share a single sample buffer
-        State->Delay.Mask = length[0] - 1;
-        State->Delay.Line = &State->SampleBuffer[0];
-        totalLength = length[0];
-        for(index = 0;index < 4;index++)
-        {
-            State->Early.Delay[index].Mask = length[1 + index] - 1;
-            State->Early.Delay[index].Line = &State->SampleBuffer[totalLength];
-            totalLength += length[1 + index];
-        }
-        for(index = 0;index < 4;index++)
-        {
-            State->Late.ApDelay[index].Mask = length[5 + index] - 1;
-            State->Late.ApDelay[index].Line = &State->SampleBuffer[totalLength];
-            totalLength += length[5 + index];
-        }
-        for(index = 0;index < 4;index++)
-        {
-            State->Late.Delay[index].Mask = length[9 + index] - 1;
-            State->Late.Delay[index].Line = &State->SampleBuffer[totalLength];
-            totalLength += length[9 + index];
-        }
+        alSetError(AL_OUT_OF_MEMORY);
+        return AL_FALSE;
     }
 
+    // The early reflection and late all-pass filter line lengths are static,
+    // so their offsets only need to be calculated once.
     for(index = 0;index < 4;index++)
     {
         State->Early.Offset[index] = (ALuint)(EARLY_LINE_LENGTH[index] *
-                                              Device->Frequency);
+                                              frequency);
         State->Late.ApOffset[index] = (ALuint)(ALLPASS_LINE_LENGTH[index] *
-                                               Device->Frequency);
+                                               frequency);
     }
 
-    for(index = 0;index < State->TotalLength;index++)
-        State->SampleBuffer[index] = 0.0f;
-
     return AL_TRUE;
 }
 
-// This updates the reverb state.  This is called any time the reverb effect
-// is loaded into a slot.
-ALvoid VerbUpdate(ALeffectState *effect, ALCcontext *Context, const ALeffect *Effect)
+// This updates the device-dependant EAX reverb state.  This is called on
+// initialization and any time the device parameters (eg. playback frequency,
+// format) have been changed.
+static ALboolean EAXVerbDeviceUpdate(ALeffectState *effect, ALCdevice *Device)
 {
     ALverbState *State = (ALverbState*)effect;
-    ALuint frequency = Context->Device->Frequency;
-    ALuint index;
-    ALfloat length, mixCoeff, cw, g, coeff;
-    ALfloat hfRatio = Effect->Reverb.DecayHFRatio;
-
-    // Calculate the master low-pass filter (from the master effect HF gain).
-    cw = cos(2.0*M_PI * Effect->Reverb.HFReference / frequency);
-    g = __max(Effect->Reverb.GainHF, 0.0001f);
-    State->LpFilter.coeff = 0.0f;
-    if(g < 0.9999f) // 1-epsilon
-        State->LpFilter.coeff = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
+    ALuint frequency = Device->Frequency, index;
 
-    // Calculate the initial delay taps.
-    length = Effect->Reverb.ReflectionsDelay;
-    State->Tap[0] = (ALuint)(length * frequency);
+    // Allocate the delay lines.
+    if(!AllocLines(AL_TRUE, frequency, State))
+    {
+        alSetError(AL_OUT_OF_MEMORY);
+        return AL_FALSE;
+    }
 
-    length += Effect->Reverb.LateReverbDelay;
+    // Calculate the modulation filter coefficient.  Notice that the exponent
+    // is calculated given the current sample rate.  This ensures that the
+    // resulting filter response over time is consistent across all sample
+    // rates.
+    State->Mod.Coeff = pow(MODULATION_FILTER_COEFF, MODULATION_FILTER_CONST /
+                                                    frequency);
 
-    /* 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
-     */
+    // The early reflection and late all-pass filter line lengths are static,
+    // so their offsets only need to be calculated once.
     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 * frequency);
+        State->Early.Offset[index] = (ALuint)(EARLY_LINE_LENGTH[index] *
+                                              frequency);
+        State->Late.ApOffset[index] = (ALuint)(ALLPASS_LINE_LENGTH[index] *
+                                               frequency);
     }
 
-    // Calculate the early reflections gain (from the master effect gain, and
-    // reflections gain parameters).
-    State->Early.Gain = Effect->Reverb.Gain * Effect->Reverb.ReflectionsGain;
+    // The echo all-pass filter line length is static, so its offset only
+    // needs to be calculated once.
+    State->Echo.ApOffset = (ALuint)(ECHO_ALLPASS_LENGTH * frequency);
 
-    // Calculate the gain (coefficient) for each early delay line.
-    for(index = 0;index < 4;index++)
-        State->Early.Coeff[index] = pow(10.0f, EARLY_LINE_LENGTH[index] /
-                                               Effect->Reverb.LateReverbDelay *
-                                               -60.0f / 20.0f);
+    return AL_TRUE;
+}
 
-    // Calculate the first mixing matrix coefficient (x).
-    mixCoeff = 1.0f - (0.5f * pow(Effect->Reverb.Diffusion, 3.0f));
+// This updates the reverb state.  This is called any time the reverb effect
+// is loaded into a slot.
+static ALvoid VerbUpdate(ALeffectState *effect, ALCcontext *Context, const ALeffect *Effect)
+{
+    ALverbState *State = (ALverbState*)effect;
+    ALuint frequency = Context->Device->Frequency;
+    ALfloat cw, g, x, y, hfRatio;
 
-    // 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 delay line coefficients, this gain needs to be
-    // attenuated by the 'x' mix coefficient from above.
-    State->Late.Gain = Effect->Reverb.Gain * Effect->Reverb.LateReverbGain * mixCoeff;
+    // Calculate the master low-pass filter (from the master effect HF gain).
+    cw = CalcI3DL2HFreq(Effect->Reverb.HFReference, frequency);
+    g = __max(Effect->Reverb.GainHF, 0.0001f);
+    // This is done with 2 chained 1-pole filters, so no need to square g.
+    State->LpFilter.coeff = CalcI3DL2Coeff(g, cw);
 
-    /* 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).
-     */
-    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 * 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 * length));
-    State->Late.DensityGain = __min(aluSqrt(g / length), 1.0f);
+    // Update the initial effect delay.
+    UpdateDelayLine(Effect->Reverb.ReflectionsDelay,
+                    Effect->Reverb.LateReverbDelay, frequency, State);
 
-    // Calculate the all-pass feed-back and feed-forward coefficient.
-    State->Late.ApFeedCoeff = 0.6f * pow(Effect->Reverb.Diffusion, 3.0f);
+    // Update the early lines.
+    UpdateEarlyLines(Effect->Reverb.Gain, Effect->Reverb.ReflectionsGain,
+                     Effect->Reverb.LateReverbDelay, State);
 
-    // Calculate the mixing matrix coefficient (y / x).
-    g = aluSqrt((1.0f - (mixCoeff * mixCoeff)) / 3.0f);
-    State->Late.MixCoeff = g / mixCoeff;
+    // Update the decorrelator.
+    UpdateDecorrelator(Effect->Reverb.Density, frequency, State);
 
-    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);
-    }
+    // Get the mixing matrix coefficients (x and y).
+    CalcMatrixCoeffs(Effect->Reverb.Diffusion, &x, &y);
+    // Then divide x into y to simplify the matrix calculation.
+    State->Late.MixCoeff = y / x;
 
     // If the HF limit parameter is flagged, calculate an appropriate limit
     // based on the air absorption parameter.
+    hfRatio = Effect->Reverb.DecayHFRatio;
     if(Effect->Reverb.DecayHFLimit && Effect->Reverb.AirAbsorptionGainHF < 1.0f)
-    {
-        ALfloat limitRatio;
-
-        // 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
-        // calculated once for all lines.
-        limitRatio = 1.0f / (log10(Effect->Reverb.AirAbsorptionGainHF) *
-                             SPEEDOFSOUNDMETRESPERSEC *
-                             Effect->Reverb.DecayTime / -60.0f * 20.0f);
-        // Need to limit the result to a minimum of 0.1, just like the HF
-        // ratio parameter.
-        limitRatio = __max(limitRatio, 0.1f);
-
-        // Using the limit calculated above, apply the upper bound to the
-        // HF ratio.
-        hfRatio = __min(hfRatio, limitRatio);
-    }
+        hfRatio = CalcLimitedHfRatio(hfRatio, Effect->Reverb.AirAbsorptionGainHF,
+                                     Effect->Reverb.DecayTime);
 
-    // Calculate the low-pass filter frequency.
-    cw = cos(2.0*M_PI * Effect->Reverb.HFReference / frequency);
+    // Update the late lines.
+    UpdateLateLines(Effect->Reverb.Gain, Effect->Reverb.LateReverbGain,
+                    x, Effect->Reverb.Density, Effect->Reverb.DecayTime,
+                    Effect->Reverb.Diffusion, hfRatio, cw, frequency, State);
+}
 
-    for(index = 0;index < 4;index++)
-    {
-        // 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 * frequency);
+// This updates the EAX reverb state.  This is called any time the EAX reverb
+// effect is loaded into a slot.
+static ALvoid EAXVerbUpdate(ALeffectState *effect, ALCcontext *Context, const ALeffect *Effect)
+{
+    ALverbState *State = (ALverbState*)effect;
+    ALuint frequency = Context->Device->Frequency;
+    ALfloat cw, g, x, y, hfRatio;
 
-        // Calculate the gain (coefficient) for each cyclical line.
-        State->Late.Coeff[index] = pow(10.0f, length / Effect->Reverb.DecayTime *
-                                              -60.0f / 20.0f);
-
-        // Eventually this should boost the high frequencies when the ratio
-        // exceeds 1.
-        coeff = 0.0f;
-        if (hfRatio < 1.0f)
-        {
-            // 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.
-            if(g < 0.9999f) // 1-epsilon
-                coeff = (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.
-            coeff = __min(coeff, 0.98f);
-        }
-        State->Late.LpCoeff[index] = coeff;
+    // Calculate the master low-pass filter (from the master effect HF gain).
+    cw = CalcI3DL2HFreq(Effect->Reverb.HFReference, frequency);
+    g = __max(Effect->Reverb.GainHF, 0.0001f);
+    // This is done with 2 chained 1-pole filters, so no need to square g.
+    State->LpFilter.coeff = CalcI3DL2Coeff(g, cw);
 
-        // Attenuate the cyclical line coefficients by the mixing coefficient
-        // (x).
-        State->Late.Coeff[index] *= mixCoeff;
-    }
+    // Update the modulator line.
+    UpdateModulator(Effect->Reverb.ModulationTime,
+                    Effect->Reverb.ModulationDepth, frequency, State);
 
-    // Calculate the 3D-panning gains for the early reflections and late
-    // reverb (for EAX mode).
-    {
-        ALfloat earlyPan[3] = { Effect->Reverb.ReflectionsPan[0], Effect->Reverb.ReflectionsPan[1], Effect->Reverb.ReflectionsPan[2] };
-        ALfloat latePan[3] = { Effect->Reverb.LateReverbPan[0], Effect->Reverb.LateReverbPan[1], Effect->Reverb.LateReverbPan[2] };
-        ALfloat *speakerGain, dirGain, ambientGain;
-        ALfloat length;
-        ALint pos;
-
-        length = earlyPan[0]*earlyPan[0] + earlyPan[1]*earlyPan[1] + earlyPan[2]*earlyPan[2];
-        if(length > 1.0f)
-        {
-            length = 1.0f / aluSqrt(length);
-            earlyPan[0] *= length;
-            earlyPan[1] *= length;
-            earlyPan[2] *= length;
-        }
-        length = latePan[0]*latePan[0] + latePan[1]*latePan[1] + latePan[2]*latePan[2];
-        if(length > 1.0f)
-        {
-            length = 1.0f / aluSqrt(length);
-            latePan[0] *= length;
-            latePan[1] *= length;
-            latePan[2] *= length;
-        }
-
-        // This code applies directional reverb just like the mixer applies
-        // directional sources.  It diffuses the sound toward all speakers
-        // as the magnitude of the panning vector drops, which is only an
-        // approximation of the expansion of sound across the speakers from
-        // the panning direction.
-        pos = aluCart2LUTpos(earlyPan[2], earlyPan[0]);
-        speakerGain = &Context->PanningLUT[OUTPUTCHANNELS * pos];
-        dirGain = aluSqrt((earlyPan[0] * earlyPan[0]) + (earlyPan[2] * earlyPan[2]));
-        ambientGain = (1.0 - dirGain);
-        for(index = 0;index < OUTPUTCHANNELS;index++)
-             State->Early.PanGain[index] = dirGain * speakerGain[index] + ambientGain;
-
-        pos = aluCart2LUTpos(latePan[2], latePan[0]);
-        speakerGain = &Context->PanningLUT[OUTPUTCHANNELS * pos];
-        dirGain = aluSqrt((latePan[0] * latePan[0]) + (latePan[2] * latePan[2]));
-        ambientGain = (1.0 - dirGain);
-        for(index = 0;index < OUTPUTCHANNELS;index++)
-             State->Late.PanGain[index] = dirGain * speakerGain[index] + ambientGain;
-    }
+    // Update the initial effect delay.
+    UpdateDelayLine(Effect->Reverb.ReflectionsDelay,
+                    Effect->Reverb.LateReverbDelay, frequency, State);
+
+    // Update the early lines.
+    UpdateEarlyLines(Effect->Reverb.Gain, Effect->Reverb.ReflectionsGain,
+                     Effect->Reverb.LateReverbDelay, State);
+
+    // Update the decorrelator.
+    UpdateDecorrelator(Effect->Reverb.Density, frequency, State);
+
+    // Get the mixing matrix coefficients (x and y).
+    CalcMatrixCoeffs(Effect->Reverb.Diffusion, &x, &y);
+    // Then divide x into y to simplify the matrix calculation.
+    State->Late.MixCoeff = y / x;
+
+    // If the HF limit parameter is flagged, calculate an appropriate limit
+    // based on the air absorption parameter.
+    hfRatio = Effect->Reverb.DecayHFRatio;
+    if(Effect->Reverb.DecayHFLimit && Effect->Reverb.AirAbsorptionGainHF < 1.0f)
+        hfRatio = CalcLimitedHfRatio(hfRatio, Effect->Reverb.AirAbsorptionGainHF,
+                                     Effect->Reverb.DecayTime);
+
+    // Update the late lines.
+    UpdateLateLines(Effect->Reverb.Gain, Effect->Reverb.LateReverbGain,
+                    x, Effect->Reverb.Density, Effect->Reverb.DecayTime,
+                    Effect->Reverb.Diffusion, hfRatio, cw, frequency, State);
+
+    // Update the echo line.
+    UpdateEchoLine(Effect->Reverb.Gain, Effect->Reverb.LateReverbGain,
+                   Effect->Reverb.EchoTime, Effect->Reverb.DecayTime,
+                   Effect->Reverb.Diffusion, Effect->Reverb.EchoDepth,
+                   hfRatio, cw, frequency, State);
+
+    // Update early and late 3D panning.
+    Update3DPanning(Effect->Reverb.ReflectionsPan, Effect->Reverb.LateReverbPan,
+                    Context->PanningLUT, State);
 }
 
 // This processes the reverb state, given the input samples and an output
 // buffer.
-ALvoid VerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint SamplesToDo, const ALfloat *SamplesIn, ALfloat (*SamplesOut)[OUTPUTCHANNELS])
+static ALvoid VerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint SamplesToDo, const ALfloat *SamplesIn, ALfloat (*SamplesOut)[OUTPUTCHANNELS])
 {
     ALverbState *State = (ALverbState*)effect;
     ALuint index;
@@ -663,7 +1174,7 @@ ALvoid VerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint Sampl
     for(index = 0;index < SamplesToDo;index++)
     {
         // Process reverb for this sample.
-        ReverbInOut(State, SamplesIn[index], early, late);
+        VerbPass(State, SamplesIn[index], early, late);
 
         // Mix early reflections and late reverb.
         out[0] = (early[0] + late[0]) * gain;
@@ -685,7 +1196,7 @@ ALvoid VerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint Sampl
 
 // This processes the EAX reverb state, given the input samples and an output
 // buffer.
-ALvoid EAXVerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint SamplesToDo, const ALfloat *SamplesIn, ALfloat (*SamplesOut)[OUTPUTCHANNELS])
+static ALvoid EAXVerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint SamplesToDo, const ALfloat *SamplesIn, ALfloat (*SamplesOut)[OUTPUTCHANNELS])
 {
     ALverbState *State = (ALverbState*)effect;
     ALuint index;
@@ -695,7 +1206,7 @@ ALvoid EAXVerbProcess(ALeffectState *effect, const ALeffectslot *Slot, ALuint Sa
     for(index = 0;index < SamplesToDo;index++)
     {
         // Process reverb for this sample.
-        ReverbInOut(State, SamplesIn[index], early, late);
+        EAXVerbPass(State, SamplesIn[index], early, late);
 
         // Unfortunately, while the number and configuration of gains for
         // panning adjust according to OUTPUTCHANNELS, the output from the
@@ -746,20 +1257,25 @@ ALeffectState *VerbCreate(void)
     State->state.Update = VerbUpdate;
     State->state.Process = VerbProcess;
 
-    State->TotalLength = 0;
+    State->TotalSamples = 0;
     State->SampleBuffer = NULL;
 
     State->LpFilter.coeff = 0.0f;
     State->LpFilter.history[0] = 0.0f;
     State->LpFilter.history[1] = 0.0f;
+
+    State->Mod.Delay.Mask = 0;
+    State->Mod.Delay.Line = NULL;
+    State->Mod.Index = 0;
+    State->Mod.Range = 1;
+    State->Mod.Depth = 0.0f;
+    State->Mod.Coeff = 0.0f;
+    State->Mod.Filter = 0.0f;
+
     State->Delay.Mask = 0;
     State->Delay.Line = NULL;
-
-    State->Tap[0] = 0;
-    State->Tap[1] = 0;
-    State->Tap[2] = 0;
-    State->Tap[3] = 0;
-    State->Tap[4] = 0;
+    State->DelayTap[0] = 0;
+    State->DelayTap[1] = 0;
 
     State->Early.Gain = 0.0f;
     for(index = 0;index < 4;index++)
@@ -770,11 +1286,16 @@ ALeffectState *VerbCreate(void)
         State->Early.Offset[index] = 0;
     }
 
+    State->Decorrelator.Mask = 0;
+    State->Decorrelator.Line = NULL;
+    State->DecoTap[0] = 0;
+    State->DecoTap[1] = 0;
+    State->DecoTap[2] = 0;
+
     State->Late.Gain = 0.0f;
     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;
@@ -791,20 +1312,40 @@ ALeffectState *VerbCreate(void)
         State->Late.LpSample[index] = 0.0f;
     }
 
-    // Panning is applied as an independent gain for each output channel.
     for(index = 0;index < OUTPUTCHANNELS;index++)
     {
         State->Early.PanGain[index] = 0.0f;
         State->Late.PanGain[index] = 0.0f;
     }
 
+    State->Echo.DensityGain = 0.0f;
+    State->Echo.Delay.Mask = 0;
+    State->Echo.Delay.Line = NULL;
+    State->Echo.ApDelay.Mask = 0;
+    State->Echo.ApDelay.Line = NULL;
+    State->Echo.Coeff = 0.0f;
+    State->Echo.ApFeedCoeff = 0.0f;
+    State->Echo.ApCoeff = 0.0f;
+    State->Echo.Offset = 0;
+    State->Echo.ApOffset = 0;
+    State->Echo.LpCoeff = 0.0f;
+    State->Echo.LpSample = 0.0f;
+    State->Echo.MixCoeff[0] = 0.0f;
+    State->Echo.MixCoeff[1] = 1.0f;
+
     State->Offset = 0;
+
     return &State->state;
 }
 
 ALeffectState *EAXVerbCreate(void)
 {
     ALeffectState *State = VerbCreate();
-    if(State) State->Process = EAXVerbProcess;
+    if(State)
+    {
+        State->DeviceUpdate = EAXVerbDeviceUpdate;
+        State->Update = EAXVerbUpdate;
+        State->Process = EAXVerbProcess;
+    }
     return State;
 }