Apply phase correction to the ambisonic decoder HRIRs

This preserves the original phase of the HRIR frequencies for decoding the
ambisonic signal. This should help prevent extra coloration from the band-
splitter used to scale the HF response.

1 files changed, 60 insertions, 25 deletions

diff --git a/Alc/hrtf.cpp b/Alc/hrtf.cpp
index 96d5ec43..b6b5ce43 100644
--- a/Alc/hrtf.cpp
+++ b/Alc/hrtf.cpp
@@ -75,10 +75,12 @@ using namespace std::placeholders;
 
 using HrtfHandlePtr = std::unique_ptr<HrtfHandle>;
 
-/* Current data set limits defined by the makehrtf utility. */
+/* Data set limits must be the same as or more flexible than those defined in
+ * the makehrtf utility.
+ */
 #define MIN_IR_SIZE                  (8)
 #define MAX_IR_SIZE                  (512)
-#define MOD_IR_SIZE                  (8)
+#define MOD_IR_SIZE                  (2)
 
 #define MIN_FD_COUNT                 (1)
 #define MAX_FD_COUNT                 (16)
@@ -300,9 +302,9 @@ void BuildBFormatHrtf(const HrtfEntry *Hrtf, DirectHrtfState *state, const ALsiz
     static constexpr int OrderFromChan[MAX_AMBI_COEFFS]{
         0, 1,1,1, 2,2,2,2,2, 3,3,3,3,3,3,3,
     };
-    /* Set this to true for dual-band HRTF processing. May require a higher
-     * quality filter, or better calculation of the new IR length to deal with
-     * the tail generated by the filter.
+    /* Set this to true for dual-band HRTF processing. May require better
+     * calculation of the new IR length to deal with the head and tail
+     * generated by the HF scaling.
      */
     static constexpr bool DualBand{true};
 
@@ -315,7 +317,7 @@ void BuildBFormatHrtf(const HrtfEntry *Hrtf, DirectHrtfState *state, const ALsiz
             std::bind(std::mem_fn(&HrtfEntry::Field::evCount), _2)))};
     ALsizei min_delay{HRTF_HISTORY_LENGTH};
     ALsizei max_delay{0};
-    al::vector<ALsizei> idx(AmbiCount);
+    auto idx = al::vector<ALsizei>(AmbiCount);
     auto calc_idxs = [Hrtf,ebase,&field,&max_delay,&min_delay](const AngularPoint &pt) noexcept -> ALsizei
     {
         /* Calculate elevation index. */
@@ -339,17 +341,23 @@ void BuildBFormatHrtf(const HrtfEntry *Hrtf, DirectHrtfState *state, const ALsiz
     };
     std::transform(AmbiPoints, AmbiPoints+AmbiCount, idx.begin(), calc_idxs);
 
+    /* For dual-band processing, add a 12-sample delay to compensate for the HF
+     * scale on the minimum-phase response.
+     */
+    static constexpr ALsizei base_delay{DualBand ? 12 : 0};
     const ALdouble xover_norm{400.0 / Hrtf->sampleRate};
     BandSplitterR<double> splitter;
+    SplitterAllpassR<double> allpass;
     splitter.init(xover_norm);
+    allpass.init(xover_norm);
 
     auto tmpres = al::vector<HrirArray<ALdouble>>(NumChannels);
-    auto tmpfilt = al::vector<std::array<ALdouble,HRIR_LENGTH>>(3);
+    auto tmpfilt = al::vector<std::array<ALdouble,HRIR_LENGTH*4>>(3);
     for(ALsizei c{0};c < AmbiCount;++c)
     {
         const ALfloat (*fir)[2]{&Hrtf->coeffs[idx[c] * Hrtf->irSize]};
-        const ALsizei ldelay{Hrtf->delays[idx[c]][0] - min_delay};
-        const ALsizei rdelay{Hrtf->delays[idx[c]][1] - min_delay};
+        const ALsizei ldelay{Hrtf->delays[idx[c]][0] - min_delay + base_delay};
+        const ALsizei rdelay{Hrtf->delays[idx[c]][1] - min_delay + base_delay};
 
         if(!DualBand)
         {
@@ -369,35 +377,62 @@ void BuildBFormatHrtf(const HrtfEntry *Hrtf, DirectHrtfState *state, const ALsiz
             continue;
         }
 
-        /* Split the left HRIR into low and high frequency bands. */
-        auto tmpfilt_iter = std::transform(fir, fir+Hrtf->irSize, tmpfilt[2].begin(),
-            [](const ALfloat (&ir)[2]) noexcept { return ir[0]; });
-        std::fill(tmpfilt_iter, tmpfilt[2].end(), 0.0);
+        /* For dual-band processing, the HRIR needs to be split into low and
+         * high frequency responses. The band-splitter alone creates frequency-
+         * dependent phase-shifts, which is not ideal. To counteract it,
+         * combine it with a backwards phase-shift.
+         */
+
+        /* Load the (left) HRIR backwards, into a temp buffer with padding. */
+        std::fill(tmpfilt[2].begin(), tmpfilt[2].end(), 0.0);
+        std::transform(fir, fir+Hrtf->irSize, tmpfilt[2].rbegin() + HRIR_LENGTH*3,
+            [](const ALfloat (&ir)[2]) noexcept -> ALdouble { return ir[0]; });
+
+        /* Apply the all-pass on the reversed signal and reverse the resulting
+         * sample array. This produces the forward response with a backwards
+         * phase-shift (+n degrees becomes -n degrees).
+         */
+        allpass.clear();
+        allpass.process(tmpfilt[2].data(), tmpfilt[2].size());
+        std::reverse(tmpfilt[2].begin(), tmpfilt[2].end());
+
+        /* Now apply the band-splitter. This applies the normal phase-shift,
+         * which cancels out with the backwards phase-shift to get the original
+         * phase on the split signal.
+         */
         splitter.clear();
-        splitter.process(tmpfilt[0].data(), tmpfilt[1].data(), tmpfilt[2].data(), HRIR_LENGTH);
+        splitter.process(tmpfilt[0].data(), tmpfilt[1].data(), tmpfilt[2].data(),
+            tmpfilt[2].size());
 
         /* Apply left ear response with delay and HF scale. */
         for(ALsizei i{0};i < NumChannels;++i)
         {
             const ALdouble mult{AmbiMatrix[c][i]};
             const ALdouble hfgain{AmbiOrderHFGain[OrderFromChan[i]]};
-            for(ALsizei lidx{ldelay},j{0};lidx < HRIR_LENGTH;++lidx,++j)
+            ALsizei j{HRIR_LENGTH*3 - ldelay};
+            for(ALsizei lidx{0};lidx < HRIR_LENGTH;++lidx,++j)
                 tmpres[i][lidx][0] += (tmpfilt[0][j]*hfgain + tmpfilt[1][j]) * mult;
         }
 
-        /* Split the right HRIR into low and high frequency bands. */
-        tmpfilt_iter = std::transform(fir, fir+Hrtf->irSize, tmpfilt[2].begin(),
-            [](const ALfloat (&ir)[2]) noexcept { return ir[1]; });
-        std::fill(tmpfilt_iter, tmpfilt[2].end(), 0.0);
+        /* Now run the same process on the right HRIR. */
+        std::fill(tmpfilt[2].begin(), tmpfilt[2].end(), 0.0);
+        std::transform(fir, fir+Hrtf->irSize, tmpfilt[2].rbegin() + HRIR_LENGTH*3,
+            [](const ALfloat (&ir)[2]) noexcept -> ALdouble { return ir[1]; });
+
+        allpass.clear();
+        allpass.process(tmpfilt[2].data(), tmpfilt[2].size());
+        std::reverse(tmpfilt[2].begin(), tmpfilt[2].end());
+
         splitter.clear();
-        splitter.process(tmpfilt[0].data(), tmpfilt[1].data(), tmpfilt[2].data(), HRIR_LENGTH);
+        splitter.process(tmpfilt[0].data(), tmpfilt[1].data(), tmpfilt[2].data(),
+            tmpfilt[2].size());
 
-        /* Apply right ear response with delay and HF scale. */
         for(ALsizei i{0};i < NumChannels;++i)
         {
             const ALdouble mult{AmbiMatrix[c][i]};
             const ALdouble hfgain{AmbiOrderHFGain[OrderFromChan[i]]};
-            for(ALsizei ridx{rdelay},j{0};ridx < HRIR_LENGTH;++ridx,++j)
+            ALsizei j{HRIR_LENGTH*3 - rdelay};
+            for(ALsizei ridx{0};ridx < HRIR_LENGTH;++ridx,++j)
                 tmpres[i][ridx][1] += (tmpfilt[0][j]*hfgain + tmpfilt[1][j]) * mult;
         }
     }
@@ -414,10 +449,10 @@ void BuildBFormatHrtf(const HrtfEntry *Hrtf, DirectHrtfState *state, const ALsiz
     tmpres.clear();
 
     ALsizei max_length{HRIR_LENGTH};
-    /* Increase the IR size by 2/3rds with dual-band processing to account for
-     * the tail generated by the filter.
+    /* Increase the IR size by 24 samples with dual-band processing to account
+     * for the head and tail from the HF response scale.
      */
-    const ALsizei irsize{DualBand ? mini(Hrtf->irSize*5/3, max_length) : Hrtf->irSize};
+    const ALsizei irsize{DualBand ? mini(Hrtf->irSize + base_delay*2, max_length) : Hrtf->irSize};
     max_length = mini(max_delay-min_delay + irsize, max_length);
 
     /* Round up to the next IR size multiple. */