Apply biquad and T60 filters using transposed direct form II

5 files changed, 62 insertions, 74 deletions

diff --git a/Alc/effects/echo.c b/Alc/effects/echo.c
index a3a4e120..282aa1d7 100644
--- a/Alc/effects/echo.c
+++ b/Alc/effects/echo.c
@@ -155,14 +155,12 @@ static ALvoid ALechoState_process(ALechoState *state, ALsizei SamplesToDo, const
     const ALsizei tap2 = state->Tap[1].delay;
     ALfloat *restrict delaybuf = state->SampleBuffer;
     ALsizei offset = state->Offset;
-    ALfloat x[2], y[2], in, out;
+    ALfloat z1, z2, in, out;
     ALsizei base;
     ALsizei c, i;
 
-    x[0] = state->Filter.x[0];
-    x[1] = state->Filter.x[1];
-    y[0] = state->Filter.y[0];
-    y[1] = state->Filter.y[1];
+    z1 = state->Filter.z1;
+    z2 = state->Filter.z2;
     for(base = 0;base < SamplesToDo;)
     {
         alignas(16) ALfloat temps[2][128];
@@ -182,11 +180,9 @@ static ALvoid ALechoState_process(ALechoState *state, ALsizei SamplesToDo, const
              * feedback attenuation.
              */
             in = temps[1][i];
-            out = in*state->Filter.b0 +
-                  x[0]*state->Filter.b1 + x[1]*state->Filter.b2 -
-                  y[0]*state->Filter.a1 - y[1]*state->Filter.a2;
-            x[1] = x[0]; x[0] = in;
-            y[1] = y[0]; y[0] = out;
+            out = in*state->Filter.b0 + z1;
+            z1 = in*state->Filter.b1 - out*state->Filter.a1 + z2;
+            z2 = in*state->Filter.b2 - out*state->Filter.a2;
 
             delaybuf[offset&mask] += out * state->FeedGain;
             offset++;
@@ -198,10 +194,8 @@ static ALvoid ALechoState_process(ALechoState *state, ALsizei SamplesToDo, const
 
         base += td;
     }
-    state->Filter.x[0] = x[0];
-    state->Filter.x[1] = x[1];
-    state->Filter.y[0] = y[0];
-    state->Filter.y[1] = y[1];
+    state->Filter.z1 = z1;
+    state->Filter.z2 = z2;
 
     state->Offset = offset;
 }
diff --git a/Alc/effects/reverb.c b/Alc/effects/reverb.c
index 0d486e63..831e4a3a 100644
--- a/Alc/effects/reverb.c
+++ b/Alc/effects/reverb.c
@@ -234,11 +234,9 @@ typedef struct T60Filter {
     ALfloat HFCoeffs[3];
     ALfloat LFCoeffs[3];
 
-    /* The HF and LF filters each keep a state of the last input and last
-     * output sample.
-     */
-    ALfloat HFState[2];
-    ALfloat LFState[2];
+    /* The HF and LF filters each keep a delay component. */
+    ALfloat HFState;
+    ALfloat LFState;
 } T60Filter;
 
 typedef struct EarlyReflections {
@@ -407,10 +405,8 @@ static void ALreverbState_Construct(ALreverbState *state)
             state->Late.T60[i].HFCoeffs[j] = 0.0f;
             state->Late.T60[i].LFCoeffs[j] = 0.0f;
         }
-        state->Late.T60[i].HFState[0] = 0.0f;
-        state->Late.T60[i].HFState[1] = 0.0f;
-        state->Late.T60[i].LFState[0] = 0.0f;
-        state->Late.T60[i].LFState[1] = 0.0f;
+        state->Late.T60[i].HFState = 0.0f;
+        state->Late.T60[i].LFState = 0.0f;
     }
 
     for(i = 0;i < NUM_LINES;i++)
@@ -1443,9 +1439,8 @@ DECL_TEMPLATE(Faded)
 static inline ALfloat FirstOrderFilter(const ALfloat in, const ALfloat *restrict coeffs,
                                        ALfloat *restrict state)
 {
-    ALfloat out = coeffs[0]*in + coeffs[1]*state[0] + coeffs[2]*state[1];
-    state[0] = in;
-    state[1] = out;
+    ALfloat out = coeffs[0]*in + *state;
+    *state = coeffs[1]*in + coeffs[2]*out;
     return out;
 }
 
@@ -1456,8 +1451,8 @@ static inline void LateT60Filter(ALfloat *restrict out, const ALfloat *restrict
     ALsizei i;
     for(i = 0;i < NUM_LINES;i++)
         out[i] = FirstOrderFilter(
-            FirstOrderFilter(in[i], filter[i].HFCoeffs, filter[i].HFState),
-            filter[i].LFCoeffs, filter[i].LFState
+            FirstOrderFilter(in[i], filter[i].HFCoeffs, &filter[i].HFState),
+            filter[i].LFCoeffs, &filter[i].LFState
         );
 }
 
diff --git a/Alc/filters/defs.h b/Alc/filters/defs.h
index c26427fb..8257cec7 100644
--- a/Alc/filters/defs.h
+++ b/Alc/filters/defs.h
@@ -31,10 +31,10 @@ typedef enum BiquadType {
 } BiquadType;
 
 typedef struct BiquadState {
-    ALfloat x[2]; /* History of two last input samples  */
-    ALfloat y[2]; /* History of two last output samples */
-    ALfloat b0, b1, b2; /* Transfer function coefficients "b" */
-    ALfloat a1, a2; /* Transfer function coefficients "a" (a0 is pre-applied) */
+    ALfloat z1, z2; /* Last two delayed components for direct form II. */
+    ALfloat b0, b1, b2; /* Transfer function coefficients "b" (numerator) */
+    ALfloat a1, a2; /* Transfer function coefficients "a" (denominator; a0 is
+                     * pre-applied). */
 } BiquadState;
 /* Currently only a C-based filter process method is implemented. */
 #define BiquadState_process BiquadState_processC
@@ -63,10 +63,8 @@ inline ALfloat calc_rcpQ_from_bandwidth(ALfloat f0norm, ALfloat bandwidth)
 
 inline void BiquadState_clear(BiquadState *filter)
 {
-    filter->x[0] = 0.0f;
-    filter->x[1] = 0.0f;
-    filter->y[0] = 0.0f;
-    filter->y[1] = 0.0f;
+    filter->z1 = 0.0f;
+    filter->z2 = 0.0f;
 }
 
 /**
@@ -97,21 +95,17 @@ inline void BiquadState_copyParams(BiquadState *restrict dst, const BiquadState
 
 void BiquadState_processC(BiquadState *filter, ALfloat *restrict dst, const ALfloat *restrict src, ALsizei numsamples);
 
-inline void BiquadState_processPassthru(BiquadState *filter, const ALfloat *restrict src, ALsizei numsamples)
+inline void BiquadState_processPassthru(BiquadState *filter, ALsizei numsamples)
 {
-    if(numsamples >= 2)
+    if(LIKELY(numsamples >= 2))
     {
-        filter->x[1] = src[numsamples-2];
-        filter->x[0] = src[numsamples-1];
-        filter->y[1] = src[numsamples-2];
-        filter->y[0] = src[numsamples-1];
+        filter->z1 = 0.0f;
+        filter->z2 = 0.0f;
     }
     else if(numsamples == 1)
     {
-        filter->x[1] = filter->x[0];
-        filter->x[0] = src[0];
-        filter->y[1] = filter->y[0];
-        filter->y[0] = src[0];
+        filter->z1 = filter->z2;
+        filter->z2 = 0.0f;
     }
 }
 
diff --git a/Alc/filters/filter.c b/Alc/filters/filter.c
index 4d757be7..fe59977b 100644
--- a/Alc/filters/filter.c
+++ b/Alc/filters/filter.c
@@ -9,7 +9,7 @@
 
 extern inline void BiquadState_clear(BiquadState *filter);
 extern inline void BiquadState_copyParams(BiquadState *restrict dst, const BiquadState *restrict src);
-extern inline void BiquadState_processPassthru(BiquadState *filter, const ALfloat *restrict src, ALsizei numsamples);
+extern inline void BiquadState_processPassthru(BiquadState *filter, ALsizei numsamples);
 extern inline ALfloat calc_rcpQ_from_slope(ALfloat gain, ALfloat slope);
 extern inline ALfloat calc_rcpQ_from_bandwidth(ALfloat f0norm, ALfloat bandwidth);
 
@@ -96,38 +96,43 @@ void BiquadState_setParams(BiquadState *filter, BiquadType type, ALfloat gain, A
 
 void BiquadState_processC(BiquadState *filter, ALfloat *restrict dst, const ALfloat *restrict src, ALsizei numsamples)
 {
-    ALsizei i;
     if(LIKELY(numsamples > 1))
     {
-        ALfloat x0 = filter->x[0];
-        ALfloat x1 = filter->x[1];
-        ALfloat y0 = filter->y[0];
-        ALfloat y1 = filter->y[1];
+        const ALfloat a1 = filter->a1;
+        const ALfloat a2 = filter->a2;
+        const ALfloat b0 = filter->b0;
+        const ALfloat b1 = filter->b1;
+        const ALfloat b2 = filter->b2;
+        ALfloat z1 = filter->z1;
+        ALfloat z2 = filter->z2;
+        ALsizei i;
 
+        /* Processing loop is transposed direct form II. This requires less
+         * storage versus direct form I (only two delay components, instead of
+         * a four-sample history; the last two inputs and outputs), and works
+         * better for floating-point which favors summing similarly-sized
+         * values while being less bothered by overflow.
+         *
+         * See: http://www.earlevel.com/main/2003/02/28/biquads/
+         */
         for(i = 0;i < numsamples;i++)
         {
-            dst[i] = filter->b0* src[i] +
-                     filter->b1*x0 + filter->b2*x1 -
-                     filter->a1*y0 - filter->a2*y1;
-            y1 = y0; y0 = dst[i];
-            x1 = x0; x0 = src[i];
+            ALfloat input = src[i];
+            ALfloat output = input*b0 + z1;
+            z1 = input*b1 - output*a1 + z2;
+            z2 = input*b2 - output*a2;
+            dst[i] = output;
         }
 
-        filter->x[0] = x0;
-        filter->x[1] = x1;
-        filter->y[0] = y0;
-        filter->y[1] = y1;
+        filter->z1 = z1;
+        filter->z2 = z2;
     }
     else if(numsamples == 1)
     {
-        dst[0] = filter->b0 * src[0] +
-                 filter->b1 * filter->x[0] +
-                 filter->b2 * filter->x[1] -
-                 filter->a1 * filter->y[0] -
-                 filter->a2 * filter->y[1];
-        filter->x[1] = filter->x[0];
-        filter->x[0] = src[0];
-        filter->y[1] = filter->y[0];
-        filter->y[0] = dst[0];
+        ALfloat input = *src;
+        ALfloat output = input*filter->b0 + filter->z1;
+        filter->z1 = input*filter->b1 - output*filter->a1 + filter->z2;
+        filter->z2 = input*filter->b2 - output*filter->a2;
+        *dst = output;
     }
 }
diff --git a/Alc/mixvoice.c b/Alc/mixvoice.c
index d01ca781..82d1ce72 100644
--- a/Alc/mixvoice.c
+++ b/Alc/mixvoice.c
@@ -271,16 +271,16 @@ static const ALfloat *DoFilters(BiquadState *lpfilter, BiquadState *hpfilter,
     switch(type)
     {
         case AF_None:
-            BiquadState_processPassthru(lpfilter, src, numsamples);
-            BiquadState_processPassthru(hpfilter, src, numsamples);
+            BiquadState_processPassthru(lpfilter, numsamples);
+            BiquadState_processPassthru(hpfilter, numsamples);
             break;
 
         case AF_LowPass:
             BiquadState_process(lpfilter, dst, src, numsamples);
-            BiquadState_processPassthru(hpfilter, dst, numsamples);
+            BiquadState_processPassthru(hpfilter, numsamples);
             return dst;
         case AF_HighPass:
-            BiquadState_processPassthru(lpfilter, src, numsamples);
+            BiquadState_processPassthru(lpfilter, numsamples);
             BiquadState_process(hpfilter, dst, src, numsamples);
             return dst;