Decode directly from B-Format to HRTF instead of a cube

Last time this attempted to average the HRIRs according to their contribution
to a given B-Format channel as if they were loudspeakers, as well as averaging
the HRIR delays. The latter part resulted in the loss of the ITD (inter-aural
time delay), a key component of HRTF.

This time, the HRIRs are averaged similar to above, except instead of averaging
the delays, they're applied to the resulting coefficients (for example, a delay
of 8 would apply the HRIR starting at the 8th sample of the target HRIR). This
does roughly double the IR length, as the largest delay is about 35 samples
while the filter is normally 32 samples. However, this is still smaller the
original data set IR (which was 256 samples), it also only needs to be applied
to 4 channels for first-order ambisonics, rather than the 8-channel cube. So
it's doing twice as much work per sample, but only working on half the number
of samples.

Additionally, since the resulting HRIRs no longer rely on an extra delay line,
a more efficient HRTF mixing function can be made that doesn't use one. Such a
function can also avoid the per-sample stepping parameters the original uses.

4 files changed, 106 insertions, 40 deletions

diff --git a/Alc/ALu.c b/Alc/ALu.c
index 8b885121..41377eec 100644
--- a/Alc/ALu.c
+++ b/Alc/ALu.c
@@ -1505,7 +1505,7 @@ ALvoid aluMixData(ALCdevice *device, ALvoid *buffer, ALsizei size)
             if(lidx != -1 && ridx != -1)
             {
                 HrtfMixerFunc HrtfMix = SelectHrtfMixer();
-                ALuint irsize = device->Hrtf->irSize;
+                ALuint irsize = device->Hrtf_IrSize;
                 MixHrtfParams hrtfparams;
                 memset(&hrtfparams, 0, sizeof(hrtfparams));
                 for(c = 0;c < device->Dry.NumChannels;c++)
diff --git a/Alc/hrtf.c b/Alc/hrtf.c
index 9e927e4f..ed99554c 100644
--- a/Alc/hrtf.c
+++ b/Alc/hrtf.c
@@ -172,6 +172,82 @@ void GetLerpedHrtfCoeffs(const struct Hrtf *Hrtf, ALfloat elevation, ALfloat azi
 }
 
 
+ALuint BuildBFormatHrtf(const struct Hrtf *Hrtf, ALfloat (*coeffs)[HRIR_LENGTH][2], ALuint NumChannels)
+{
+    ALuint total_hrirs = 0;
+    ALuint max_length = 0;
+    ALuint eidx, aidx, i, j;
+    ALfloat scale;
+
+    assert(NumChannels == 4);
+
+    for(eidx = 0;eidx < Hrtf->evCount;++eidx)
+    {
+        const ALfloat elev = (ALfloat)eidx/(ALfloat)(Hrtf->evCount-1)*F_PI - F_PI_2;
+        const ALuint evoffset = Hrtf->evOffset[eidx];
+        const ALuint azcount = Hrtf->azCount[eidx];
+
+        for(aidx = 0;aidx < azcount;++aidx)
+        {
+            ALfloat ambcoeffs[4];
+            const ALshort *fir;
+            ALuint length, delay;
+            ALuint lidx, ridx;
+            ALfloat x, y, z;
+            ALfloat azi;
+
+            lidx = evoffset + aidx;
+            ridx = evoffset + ((azcount - aidx) % azcount);
+
+            azi = (ALfloat)aidx/(ALfloat)azcount * -F_TAU;
+            x = cosf(azi) * cosf(elev);
+            y = sinf(azi) * cosf(elev);
+            z = sinf(elev);
+
+            ambcoeffs[0] = 1.0f;
+            ambcoeffs[1] = y / 1.732050808f;
+            ambcoeffs[2] = z / 1.732050808f;
+            ambcoeffs[3] = x / 1.732050808f;
+
+            /* Apply left ear response */
+            delay = Hrtf->delays[lidx];
+            fir = &Hrtf->coeffs[lidx * Hrtf->irSize];
+            length = minu(delay + Hrtf->irSize, HRIR_LENGTH);
+            for(i = 0;i < NumChannels;++i)
+            {
+                for(j = delay;j < length;++j)
+                    coeffs[i][j][0] += fir[j-delay]/32767.0f * ambcoeffs[i];
+            }
+            max_length = maxu(max_length, length);
+
+            /* Apply right ear response */
+            delay = Hrtf->delays[ridx];
+            fir = &Hrtf->coeffs[ridx * Hrtf->irSize];
+            length = minu(delay + Hrtf->irSize, HRIR_LENGTH);
+            for(i = 0;i < NumChannels;++i)
+            {
+                for(j = delay;j < length;++j)
+                    coeffs[i][j][1] += fir[j-delay]/32767.0f * ambcoeffs[i];
+            }
+            max_length = maxu(max_length, length);
+
+            total_hrirs++;
+        }
+    }
+
+    scale = (ALfloat)total_hrirs;
+    for(i = 0;i < NumChannels;++i)
+    {
+        for(j = 0;j < max_length;++j)
+        {
+            coeffs[i][j][0] /= scale;
+            coeffs[i][j][1] /= scale;
+        }
+    }
+    return max_length;
+}
+
+
 static struct Hrtf *LoadHrtf00(FILE *f, const_al_string filename)
 {
     const ALubyte maxDelay = HRTF_HISTORY_LENGTH-1;
diff --git a/Alc/hrtf.h b/Alc/hrtf.h
index e9f04d49..d3ecfc3c 100644
--- a/Alc/hrtf.h
+++ b/Alc/hrtf.h
@@ -42,4 +42,11 @@ void FreeHrtfList(vector_HrtfEntry *list);
 
 void GetLerpedHrtfCoeffs(const struct Hrtf *Hrtf, ALfloat elevation, ALfloat azimuth, ALfloat spread, ALfloat gain, ALfloat (*coeffs)[2], ALuint *delays);
 
+/* Produces HRTF filter coefficients for decoding B-Format. The result will
+ * have ACN ordering with N3D normalization. NumChannels must currently be 4,
+ * for first-order. Returns the maximum impulse-response length of the
+ * generated coefficients.
+ */
+ALuint BuildBFormatHrtf(const struct Hrtf *Hrtf, ALfloat (*coeffs)[HRIR_LENGTH][2], ALuint NumChannels);
+
 #endif /* ALC_HRTF_H */
diff --git a/Alc/panning.c b/Alc/panning.c
index 2d7502ec..65506303 100644
--- a/Alc/panning.c
+++ b/Alc/panning.c
@@ -758,52 +758,35 @@ static void InitHQPanning(ALCdevice *device, const AmbDecConf *conf, const ALuin
 
 static void InitHrtfPanning(ALCdevice *device)
 {
-    static const enum Channel CubeChannels[MAX_OUTPUT_CHANNELS] = {
-        UpperFrontLeft, UpperFrontRight, UpperBackLeft, UpperBackRight,
-        LowerFrontLeft, LowerFrontRight, LowerBackLeft, LowerBackRight,
-        InvalidChannel, InvalidChannel, InvalidChannel, InvalidChannel,
-        InvalidChannel, InvalidChannel, InvalidChannel, InvalidChannel
-    };
-    static const ChannelMap Cube8Cfg[8] = {
-        { UpperFrontLeft,  { 0.176776695f,  0.072168784f,  0.072168784f,  0.072168784f } },
-        { UpperFrontRight, { 0.176776695f,  0.072168784f, -0.072168784f,  0.072168784f } },
-        { UpperBackLeft,   { 0.176776695f, -0.072168784f,  0.072168784f,  0.072168784f } },
-        { UpperBackRight,  { 0.176776695f, -0.072168784f, -0.072168784f,  0.072168784f } },
-        { LowerFrontLeft,  { 0.176776695f,  0.072168784f,  0.072168784f, -0.072168784f } },
-        { LowerFrontRight, { 0.176776695f,  0.072168784f, -0.072168784f, -0.072168784f } },
-        { LowerBackLeft,   { 0.176776695f, -0.072168784f,  0.072168784f, -0.072168784f } },
-        { LowerBackRight,  { 0.176776695f, -0.072168784f, -0.072168784f, -0.072168784f } },
-    };
-    static const struct {
-        enum Channel Channel;
-        ALfloat Angle;
-        ALfloat Elevation;
-    } CubeInfo[8] = {
-        { UpperFrontLeft,  DEG2RAD( -45.0f), DEG2RAD( 45.0f) },
-        { UpperFrontRight, DEG2RAD(  45.0f), DEG2RAD( 45.0f) },
-        { UpperBackLeft,   DEG2RAD(-135.0f), DEG2RAD( 45.0f) },
-        { UpperBackRight,  DEG2RAD( 135.0f), DEG2RAD( 45.0f) },
-        { LowerFrontLeft,  DEG2RAD( -45.0f), DEG2RAD(-45.0f) },
-        { LowerFrontRight, DEG2RAD(  45.0f), DEG2RAD(-45.0f) },
-        { LowerBackLeft,   DEG2RAD(-135.0f), DEG2RAD(-45.0f) },
-        { LowerBackRight,  DEG2RAD( 135.0f), DEG2RAD(-45.0f) },
-    };
-    const ChannelMap *chanmap = Cube8Cfg;
-    size_t count = COUNTOF(Cube8Cfg);
-    ALuint i;
+    ALfloat hrtf_coeffs[4][HRIR_LENGTH][2];
+    size_t count = 4;
+    ALuint i, j;
 
-    SetChannelMap(CubeChannels, device->Dry.Ambi.Coeffs, chanmap, count,
-                  &device->Dry.NumChannels, AL_TRUE);
-    device->Dry.CoeffCount = 4;
+    for(i = 0;i < count;i++)
+    {
+        device->Dry.Ambi.Map[i].Scale = 1.0f;
+        device->Dry.Ambi.Map[i].Index = i;
+    }
+    device->Dry.CoeffCount = 0;
+    device->Dry.NumChannels = count;
 
     device->FOAOut.Ambi = device->Dry.Ambi;
     device->FOAOut.CoeffCount = device->Dry.CoeffCount;
 
+    memset(hrtf_coeffs, 0, sizeof(hrtf_coeffs));
+    device->Hrtf_IrSize = BuildBFormatHrtf(device->Hrtf, hrtf_coeffs, device->Dry.NumChannels);
+
+    /* Round up to the nearest multiple of 8 */
+    device->Hrtf_IrSize = (device->Hrtf_IrSize+7)&~7;
     for(i = 0;i < device->Dry.NumChannels;i++)
     {
-        int chan = GetChannelIndex(CubeChannels, CubeInfo[i].Channel);
-        GetLerpedHrtfCoeffs(device->Hrtf, CubeInfo[i].Elevation, CubeInfo[i].Angle, 1.0f, 0.0f,
-                            device->Hrtf_Params[chan].Coeffs, device->Hrtf_Params[chan].Delay);
+        for(j = 0;j < HRIR_LENGTH;j++)
+        {
+            device->Hrtf_Params[i].Coeffs[j][0] = hrtf_coeffs[i][j][0];
+            device->Hrtf_Params[i].Coeffs[j][1] = hrtf_coeffs[i][j][1];
+        }
+        device->Hrtf_Params[i].Delay[0] = 0;
+        device->Hrtf_Params[i].Delay[1] = 0;
     }
 }