Band-split the HRIRs when building the ambisonic decoder filters

This allows each HRIR to contribute a frequency-dependent response, essentially
acting like a dual-band decoder playing over the cube speaker array.

3 files changed, 68 insertions, 31 deletions

diff --git a/Alc/bformatdec.c b/Alc/bformatdec.c
index 3c175b33..356423ce 100644
--- a/Alc/bformatdec.c
+++ b/Alc/bformatdec.c
@@ -10,14 +10,7 @@
 #include "almalloc.h"
 
 
-typedef struct BandSplitter {
-    ALfloat coeff;
-    ALfloat lp_z1;
-    ALfloat lp_z2;
-    ALfloat hp_z1;
-} BandSplitter;
-
-static void bandsplit_init(BandSplitter *splitter, ALfloat freq_mult)
+void bandsplit_init(BandSplitter *splitter, ALfloat freq_mult)
 {
     ALfloat w = freq_mult * F_TAU;
     ALfloat cw = cosf(w);
@@ -31,8 +24,15 @@ static void bandsplit_init(BandSplitter *splitter, ALfloat freq_mult)
     splitter->hp_z1 = 0.0f;
 }
 
-static void bandsplit_process(BandSplitter *splitter, ALfloat *restrict hpout, ALfloat *restrict lpout,
-                              const ALfloat *input, ALuint count)
+void bandsplit_clear(BandSplitter *splitter)
+{
+    splitter->lp_z1 = 0.0f;
+    splitter->lp_z2 = 0.0f;
+    splitter->hp_z1 = 0.0f;
+}
+
+void bandsplit_process(BandSplitter *splitter, ALfloat *restrict hpout, ALfloat *restrict lpout,
+                       const ALfloat *input, ALuint count)
 {
     ALfloat coeff, d, x;
     ALfloat z1, z2;
diff --git a/Alc/bformatdec.h b/Alc/bformatdec.h
index be5a69f6..433603dd 100644
--- a/Alc/bformatdec.h
+++ b/Alc/bformatdec.h
@@ -22,6 +22,7 @@ void bformatdec_process(struct BFormatDec *dec, ALfloat (*restrict OutBuffer)[BU
 /* Up-samples a first-order input to the decoder's configuration. */
 void bformatdec_upSample(struct BFormatDec *dec, ALfloat (*restrict OutBuffer)[BUFFERSIZE], ALfloat (*restrict InSamples)[BUFFERSIZE], ALuint InChannels, ALuint SamplesToDo);
 
+
 /* Stand-alone first-order upsampler. Kept here because it shares some stuff
  * with bformatdec.
  */
@@ -31,4 +32,18 @@ void ambiup_reset(struct AmbiUpsampler *ambiup, const ALCdevice *device);
 
 void ambiup_process(struct AmbiUpsampler *ambiup, ALfloat (*restrict OutBuffer)[BUFFERSIZE], ALuint OutChannels, ALfloat (*restrict InSamples)[BUFFERSIZE], ALuint SamplesToDo);
 
+
+/* Band splitter. Splits a signal into two phase-matching frequency bands. */
+typedef struct BandSplitter {
+    ALfloat coeff;
+    ALfloat lp_z1;
+    ALfloat lp_z2;
+    ALfloat hp_z1;
+} BandSplitter;
+
+void bandsplit_init(BandSplitter *splitter, ALfloat freq_mult);
+void bandsplit_clear(BandSplitter *splitter);
+void bandsplit_process(BandSplitter *splitter, ALfloat *restrict hpout, ALfloat *restrict lpout,
+                       const ALfloat *input, ALuint count);
+
 #endif /* BFORMATDEC_H */
diff --git a/Alc/hrtf.c b/Alc/hrtf.c
index 3be30d95..b11fb9a8 100644
--- a/Alc/hrtf.c
+++ b/Alc/hrtf.c
@@ -28,6 +28,7 @@
 #include "alMain.h"
 #include "alSource.h"
 #include "alu.h"
+#include "bformatdec.h"
 #include "hrtf.h"
 
 #include "compat.h"
@@ -187,20 +188,22 @@ ALuint BuildBFormatHrtf(const struct Hrtf *Hrtf, ALfloat (*coeffs)[HRIR_LENGTH][
         { DEG2RAD(-35.0f), DEG2RAD(-135.0f) },
         { DEG2RAD(-35.0f), DEG2RAD( 135.0f) },
     };
-    static const ALfloat CubeMatrix[8][MAX_AMBI_COEFFS] = {
-        { 0.25f,  0.14425f,  0.14425f,  0.14425f },
-        { 0.25f, -0.14425f,  0.14425f,  0.14425f },
-        { 0.25f,  0.14425f,  0.14425f, -0.14425f },
-        { 0.25f, -0.14425f,  0.14425f, -0.14425f },
-        { 0.25f,  0.14425f, -0.14425f,  0.14425f },
-        { 0.25f, -0.14425f, -0.14425f,  0.14425f },
-        { 0.25f,  0.14425f, -0.14425f, -0.14425f },
-        { 0.25f, -0.14425f, -0.14425f, -0.14425f },
+    static const ALfloat CubeMatrix[8][2][MAX_AMBI_COEFFS] = {
+        { { 0.25f,  0.14425f,  0.14425f,  0.14425f }, { 0.125f,  0.125f,  0.125f,  0.125f } },
+        { { 0.25f, -0.14425f,  0.14425f,  0.14425f }, { 0.125f, -0.125f,  0.125f,  0.125f } },
+        { { 0.25f,  0.14425f,  0.14425f, -0.14425f }, { 0.125f,  0.125f,  0.125f, -0.125f } },
+        { { 0.25f, -0.14425f,  0.14425f, -0.14425f }, { 0.125f, -0.125f,  0.125f, -0.125f } },
+        { { 0.25f,  0.14425f, -0.14425f,  0.14425f }, { 0.125f,  0.125f, -0.125f,  0.125f } },
+        { { 0.25f, -0.14425f, -0.14425f,  0.14425f }, { 0.125f, -0.125f, -0.125f,  0.125f } },
+        { { 0.25f,  0.14425f, -0.14425f, -0.14425f }, { 0.125f,  0.125f, -0.125f, -0.125f } },
+        { { 0.25f, -0.14425f, -0.14425f, -0.14425f }, { 0.125f, -0.125f, -0.125f, -0.125f } },
     };
+    BandSplitter splitter;
+    ALfloat temps[3][HRIR_LENGTH];
     ALuint lidx[8], ridx[8];
     ALuint min_delay = HRTF_HISTORY_LENGTH;
     ALuint max_length = 0;
-    ALuint i, j, c;
+    ALuint i, j, c, b;
 
     assert(NumChannels == 4);
 
@@ -229,33 +232,52 @@ ALuint BuildBFormatHrtf(const struct Hrtf *Hrtf, ALfloat (*coeffs)[HRIR_LENGTH][
         min_delay = minu(min_delay, minu(Hrtf->delays[lidx[c]], Hrtf->delays[ridx[c]]));
     }
 
+    memset(temps, 0, sizeof(temps));
+    bandsplit_init(&splitter, 400.0f / (ALfloat)Hrtf->sampleRate);
     for(c = 0;c < 8;c++)
     {
         const ALshort *fir;
-        ALuint length;
         ALuint delay;
 
+        /* Band-split left HRIR into low and high frequency responses. */
+        bandsplit_clear(&splitter);
         fir = &Hrtf->coeffs[lidx[c] * Hrtf->irSize];
+        for(i = 0;i < Hrtf->irSize;i++)
+            temps[2][i] = fir[i] / 32767.0f;
+        bandsplit_process(&splitter, temps[0], temps[1], temps[2], HRIR_LENGTH);
+
+        /* Add to the left output coefficients with the specified delay. */
         delay = Hrtf->delays[lidx[c]] - min_delay;
-        length = minu(delay + Hrtf->irSize, HRIR_LENGTH);
         for(i = 0;i < NumChannels;++i)
         {
-            ALuint k = 0;
-            for(j = delay;j < length;++j)
-                coeffs[i][j][0] += fir[k++]/32767.0f * CubeMatrix[c][i];
+            for(b = 0;b < 2;b++)
+            {
+                ALuint k = 0;
+                for(j = delay;j < HRIR_LENGTH;++j)
+                    coeffs[i][j][0] += temps[b][k++] * CubeMatrix[c][b][i];
+            }
         }
-        max_length = maxu(max_length, length);
+        max_length = maxu(max_length, minu(delay + Hrtf->irSize, HRIR_LENGTH));
 
+        /* Band-split right HRIR into low and high frequency responses. */
+        bandsplit_clear(&splitter);
         fir = &Hrtf->coeffs[ridx[c] * Hrtf->irSize];
+        for(i = 0;i < Hrtf->irSize;i++)
+            temps[2][i] = fir[i] / 32767.0f;
+        bandsplit_process(&splitter, temps[0], temps[1], temps[2], HRIR_LENGTH);
+
+        /* Add to the right output coefficients with the specified delay. */
         delay = Hrtf->delays[ridx[c]] - min_delay;
-        length = minu(delay + Hrtf->irSize, HRIR_LENGTH);
         for(i = 0;i < NumChannels;++i)
         {
-            ALuint k = 0;
-            for(j = delay;j < length;++j)
-                coeffs[i][j][1] += fir[k++]/32767.0f * CubeMatrix[c][i];
+            for(b = 0;b < 2;b++)
+            {
+                ALuint k = 0;
+                for(j = delay;j < HRIR_LENGTH;++j)
+                    coeffs[i][j][1] += temps[b][k++] * CubeMatrix[c][b][i];
+            }
         }
-        max_length = maxu(max_length, length);
+        max_length = maxu(max_length, minu(delay + Hrtf->irSize, HRIR_LENGTH));
     }
     TRACE("Skipped min delay: %u, new combined length: %u\n", min_delay, max_length);