Apply biquad and T60 filters using transposed direct form II

This commit is contained in:
Chris Robinson 2018-04-03 10:15:35 -07:00
parent 414b56edec
commit 869637af2e
5 changed files with 62 additions and 74 deletions

View File

@ -155,14 +155,12 @@ static ALvoid ALechoState_process(ALechoState *state, ALsizei SamplesToDo, const
const ALsizei tap2 = state->Tap[1].delay; const ALsizei tap2 = state->Tap[1].delay;
ALfloat *restrict delaybuf = state->SampleBuffer; ALfloat *restrict delaybuf = state->SampleBuffer;
ALsizei offset = state->Offset; ALsizei offset = state->Offset;
ALfloat x[2], y[2], in, out; ALfloat z1, z2, in, out;
ALsizei base; ALsizei base;
ALsizei c, i; ALsizei c, i;
x[0] = state->Filter.x[0]; z1 = state->Filter.z1;
x[1] = state->Filter.x[1]; z2 = state->Filter.z2;
y[0] = state->Filter.y[0];
y[1] = state->Filter.y[1];
for(base = 0;base < SamplesToDo;) for(base = 0;base < SamplesToDo;)
{ {
alignas(16) ALfloat temps[2][128]; alignas(16) ALfloat temps[2][128];
@ -182,11 +180,9 @@ static ALvoid ALechoState_process(ALechoState *state, ALsizei SamplesToDo, const
* feedback attenuation. * feedback attenuation.
*/ */
in = temps[1][i]; in = temps[1][i];
out = in*state->Filter.b0 + out = in*state->Filter.b0 + z1;
x[0]*state->Filter.b1 + x[1]*state->Filter.b2 - z1 = in*state->Filter.b1 - out*state->Filter.a1 + z2;
y[0]*state->Filter.a1 - y[1]*state->Filter.a2; z2 = in*state->Filter.b2 - out*state->Filter.a2;
x[1] = x[0]; x[0] = in;
y[1] = y[0]; y[0] = out;
delaybuf[offset&mask] += out * state->FeedGain; delaybuf[offset&mask] += out * state->FeedGain;
offset++; offset++;
@ -198,10 +194,8 @@ static ALvoid ALechoState_process(ALechoState *state, ALsizei SamplesToDo, const
base += td; base += td;
} }
state->Filter.x[0] = x[0]; state->Filter.z1 = z1;
state->Filter.x[1] = x[1]; state->Filter.z2 = z2;
state->Filter.y[0] = y[0];
state->Filter.y[1] = y[1];
state->Offset = offset; state->Offset = offset;
} }

View File

@ -234,11 +234,9 @@ typedef struct T60Filter {
ALfloat HFCoeffs[3]; ALfloat HFCoeffs[3];
ALfloat LFCoeffs[3]; ALfloat LFCoeffs[3];
/* The HF and LF filters each keep a state of the last input and last /* The HF and LF filters each keep a delay component. */
* output sample. ALfloat HFState;
*/ ALfloat LFState;
ALfloat HFState[2];
ALfloat LFState[2];
} T60Filter; } T60Filter;
typedef struct EarlyReflections { 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].HFCoeffs[j] = 0.0f;
state->Late.T60[i].LFCoeffs[j] = 0.0f; state->Late.T60[i].LFCoeffs[j] = 0.0f;
} }
state->Late.T60[i].HFState[0] = 0.0f; state->Late.T60[i].HFState = 0.0f;
state->Late.T60[i].HFState[1] = 0.0f; state->Late.T60[i].LFState = 0.0f;
state->Late.T60[i].LFState[0] = 0.0f;
state->Late.T60[i].LFState[1] = 0.0f;
} }
for(i = 0;i < NUM_LINES;i++) 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, static inline ALfloat FirstOrderFilter(const ALfloat in, const ALfloat *restrict coeffs,
ALfloat *restrict state) ALfloat *restrict state)
{ {
ALfloat out = coeffs[0]*in + coeffs[1]*state[0] + coeffs[2]*state[1]; ALfloat out = coeffs[0]*in + *state;
state[0] = in; *state = coeffs[1]*in + coeffs[2]*out;
state[1] = out;
return out; return out;
} }
@ -1456,8 +1451,8 @@ static inline void LateT60Filter(ALfloat *restrict out, const ALfloat *restrict
ALsizei i; ALsizei i;
for(i = 0;i < NUM_LINES;i++) for(i = 0;i < NUM_LINES;i++)
out[i] = FirstOrderFilter( out[i] = FirstOrderFilter(
FirstOrderFilter(in[i], filter[i].HFCoeffs, filter[i].HFState), FirstOrderFilter(in[i], filter[i].HFCoeffs, &filter[i].HFState),
filter[i].LFCoeffs, filter[i].LFState filter[i].LFCoeffs, &filter[i].LFState
); );
} }

View File

@ -31,10 +31,10 @@ typedef enum BiquadType {
} BiquadType; } BiquadType;
typedef struct BiquadState { typedef struct BiquadState {
ALfloat x[2]; /* History of two last input samples */ ALfloat z1, z2; /* Last two delayed components for direct form II. */
ALfloat y[2]; /* History of two last output samples */ ALfloat b0, b1, b2; /* Transfer function coefficients "b" (numerator) */
ALfloat b0, b1, b2; /* Transfer function coefficients "b" */ ALfloat a1, a2; /* Transfer function coefficients "a" (denominator; a0 is
ALfloat a1, a2; /* Transfer function coefficients "a" (a0 is pre-applied) */ * pre-applied). */
} BiquadState; } BiquadState;
/* Currently only a C-based filter process method is implemented. */ /* Currently only a C-based filter process method is implemented. */
#define BiquadState_process BiquadState_processC #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) inline void BiquadState_clear(BiquadState *filter)
{ {
filter->x[0] = 0.0f; filter->z1 = 0.0f;
filter->x[1] = 0.0f; filter->z2 = 0.0f;
filter->y[0] = 0.0f;
filter->y[1] = 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); 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->z1 = 0.0f;
filter->x[0] = src[numsamples-1]; filter->z2 = 0.0f;
filter->y[1] = src[numsamples-2];
filter->y[0] = src[numsamples-1];
} }
else if(numsamples == 1) else if(numsamples == 1)
{ {
filter->x[1] = filter->x[0]; filter->z1 = filter->z2;
filter->x[0] = src[0]; filter->z2 = 0.0f;
filter->y[1] = filter->y[0];
filter->y[0] = src[0];
} }
} }

View File

@ -9,7 +9,7 @@
extern inline void BiquadState_clear(BiquadState *filter); extern inline void BiquadState_clear(BiquadState *filter);
extern inline void BiquadState_copyParams(BiquadState *restrict dst, const BiquadState *restrict src); 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_slope(ALfloat gain, ALfloat slope);
extern inline ALfloat calc_rcpQ_from_bandwidth(ALfloat f0norm, ALfloat bandwidth); 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) void BiquadState_processC(BiquadState *filter, ALfloat *restrict dst, const ALfloat *restrict src, ALsizei numsamples)
{ {
ALsizei i;
if(LIKELY(numsamples > 1)) if(LIKELY(numsamples > 1))
{ {
ALfloat x0 = filter->x[0]; const ALfloat a1 = filter->a1;
ALfloat x1 = filter->x[1]; const ALfloat a2 = filter->a2;
ALfloat y0 = filter->y[0]; const ALfloat b0 = filter->b0;
ALfloat y1 = filter->y[1]; 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++) for(i = 0;i < numsamples;i++)
{ {
dst[i] = filter->b0* src[i] + ALfloat input = src[i];
filter->b1*x0 + filter->b2*x1 - ALfloat output = input*b0 + z1;
filter->a1*y0 - filter->a2*y1; z1 = input*b1 - output*a1 + z2;
y1 = y0; y0 = dst[i]; z2 = input*b2 - output*a2;
x1 = x0; x0 = src[i]; dst[i] = output;
} }
filter->x[0] = x0; filter->z1 = z1;
filter->x[1] = x1; filter->z2 = z2;
filter->y[0] = y0;
filter->y[1] = y1;
} }
else if(numsamples == 1) else if(numsamples == 1)
{ {
dst[0] = filter->b0 * src[0] + ALfloat input = *src;
filter->b1 * filter->x[0] + ALfloat output = input*filter->b0 + filter->z1;
filter->b2 * filter->x[1] - filter->z1 = input*filter->b1 - output*filter->a1 + filter->z2;
filter->a1 * filter->y[0] - filter->z2 = input*filter->b2 - output*filter->a2;
filter->a2 * filter->y[1]; *dst = output;
filter->x[1] = filter->x[0];
filter->x[0] = src[0];
filter->y[1] = filter->y[0];
filter->y[0] = dst[0];
} }
} }

View File

@ -271,16 +271,16 @@ static const ALfloat *DoFilters(BiquadState *lpfilter, BiquadState *hpfilter,
switch(type) switch(type)
{ {
case AF_None: case AF_None:
BiquadState_processPassthru(lpfilter, src, numsamples); BiquadState_processPassthru(lpfilter, numsamples);
BiquadState_processPassthru(hpfilter, src, numsamples); BiquadState_processPassthru(hpfilter, numsamples);
break; break;
case AF_LowPass: case AF_LowPass:
BiquadState_process(lpfilter, dst, src, numsamples); BiquadState_process(lpfilter, dst, src, numsamples);
BiquadState_processPassthru(hpfilter, dst, numsamples); BiquadState_processPassthru(hpfilter, numsamples);
return dst; return dst;
case AF_HighPass: case AF_HighPass:
BiquadState_processPassthru(lpfilter, src, numsamples); BiquadState_processPassthru(lpfilter, numsamples);
BiquadState_process(hpfilter, dst, src, numsamples); BiquadState_process(hpfilter, dst, src, numsamples);
return dst; return dst;