Apply the ambisonic HF scaling in real-time with HRTF

Rather than applying the HF scale to the IRs necessitating them to be truncated
along with increasing the IR size, it can be applied to the input signal for
the same results. Consequently, the IR size can be notably shortened while
avoiding the extra truncation. In its place, the delayed reversed all-pass
technique can still be used on the input for maintaining phase when applying
the bandsplit/hfscalar filter to the input signal.

4 files changed, 102 insertions, 23 deletions

diff --git a/alc/alc.cpp b/alc/alc.cpp
index 707fc34d..74fed67a 100644
--- a/alc/alc.cpp
+++ b/alc/alc.cpp
@@ -2098,11 +2098,16 @@ static ALCenum UpdateDeviceParams(ALCdevice *device, const int *attrList)
         device->SourcesMax, device->NumMonoSources, device->NumStereoSources,
         device->AuxiliaryEffectSlotMax, device->NumAuxSends);
 
-    if(Uhj2Encoder *uhj{device->Uhj_Encoder.get()})
+    if(device->Uhj_Encoder)
     {
+        /* NOTE: Don't know why this has to be "copied" into a local constexpr
+         * variable to avoid a reference on Uhj2Encoder::sFilterSize...
+         */
         constexpr size_t filter_len{Uhj2Encoder::sFilterSize};
         device->FixedLatency += nanoseconds{seconds{filter_len}} / device->Frequency;
     }
+    if(device->mHrtfState)
+        device->FixedLatency += nanoseconds{seconds{HRTF_DIRECT_DELAY}} / device->Frequency;
 
     /* Enable the stablizer only for formats that have front-left, front-right,
      * and front-center outputs.
@@ -2125,10 +2130,6 @@ static ALCenum UpdateDeviceParams(ALCdevice *device, const int *attrList)
             stablizer->MidFilter.init(5000.0f / static_cast<float>(device->Frequency));
 
             device->Stablizer = std::move(stablizer);
-            /* NOTE: Don't know why this has to be "copied" into a local static
-             * constexpr variable to avoid a reference on
-             * FrontStablizer::DelayLength...
-             */
             constexpr size_t StablizerDelay{FrontStablizer::DelayLength};
             device->FixedLatency += nanoseconds{seconds{StablizerDelay}} / device->Frequency;
         }
diff --git a/alc/hrtf.cpp b/alc/hrtf.cpp
index 26e08cba..89e07efc 100644
--- a/alc/hrtf.cpp
+++ b/alc/hrtf.cpp
@@ -295,7 +295,15 @@ void DirectHrtfState::build(const HrtfStore *Hrtf, const al::span<const AngularP
      * calculation of the new IR length to deal with the head and tail
      * generated by the HF scaling.
      */
-    static constexpr bool DualBand{true};
+    constexpr bool DualBand{false};
+    const double xover_norm{400.0 / Hrtf->sampleRate};
+
+    for(size_t i{0};i < mChannels.size();++i)
+    {
+        const size_t order{AmbiIndex::OrderFromChannel[i]};
+        mChannels[i].mSplitter.init(static_cast<float>(xover_norm));
+        mChannels[i].mHfScale = AmbiOrderHFGain[order];
+    }
 
     ALuint min_delay{HRTF_HISTORY_LENGTH*HRIR_DELAY_FRACONE};
     ALuint max_delay{0};
@@ -342,11 +350,10 @@ void DirectHrtfState::build(const HrtfStore *Hrtf, const al::span<const AngularP
     /* For dual-band processing, add a 16-sample delay to compensate for the HF
      * scale on the minimum-phase response.
      */
-    static constexpr ALuint base_delay{DualBand ? 16 : 0};
-    const double xover_norm{400.0 / Hrtf->sampleRate};
+    constexpr ALuint base_delay{DualBand ? 16 : 0};
     BandSplitterR<double> splitter{xover_norm};
 
-    auto tmpres = al::vector<std::array<double2,HRIR_LENGTH>>(mCoeffs.size());
+    auto tmpres = al::vector<std::array<double2,HRIR_LENGTH>>(mChannels.size());
     auto tmpflt = al::vector<std::array<double,HRIR_LENGTH*4>>(3);
     const al::span<double,HRIR_LENGTH*4> tempir{tmpflt[2]};
     for(size_t c{0u};c < AmbiPoints.size();++c)
@@ -357,11 +364,9 @@ void DirectHrtfState::build(const HrtfStore *Hrtf, const al::span<const AngularP
 
         if /*constexpr*/(!DualBand)
         {
-            /* For single-band decoding, apply the HF scale to the response. */
-            for(size_t i{0u};i < mCoeffs.size();++i)
+            for(size_t i{0u};i < mChannels.size();++i)
             {
-                const size_t order{AmbiIndex::OrderFromChannel[i]};
-                const double mult{double{AmbiOrderHFGain[order]} * AmbiMatrix[c][i]};
+                const double mult{AmbiMatrix[c][i]};
                 const ALuint numirs{HRIR_LENGTH - maxu(ldelay, rdelay)};
                 ALuint lidx{ldelay}, ridx{rdelay};
                 for(ALuint j{0};j < numirs;++j)
@@ -399,7 +404,7 @@ void DirectHrtfState::build(const HrtfStore *Hrtf, const al::span<const AngularP
         splitter.process(tempir, tmpflt[0].data(), tmpflt[1].data());
 
         /* Apply left ear response with delay and HF scale. */
-        for(size_t i{0u};i < mCoeffs.size();++i)
+        for(size_t i{0u};i < mChannels.size();++i)
         {
             const double mult{AmbiMatrix[c][i]};
             const double hfgain{AmbiOrderHFGain[AmbiIndex::OrderFromChannel[i]]};
@@ -419,7 +424,7 @@ void DirectHrtfState::build(const HrtfStore *Hrtf, const al::span<const AngularP
         splitter.clear();
         splitter.process(tempir, tmpflt[0].data(), tmpflt[1].data());
 
-        for(size_t i{0u};i < mCoeffs.size();++i)
+        for(size_t i{0u};i < mChannels.size();++i)
         {
             const double mult{AmbiMatrix[c][i]};
             const double hfgain{AmbiOrderHFGain[AmbiIndex::OrderFromChannel[i]]};
@@ -431,11 +436,12 @@ void DirectHrtfState::build(const HrtfStore *Hrtf, const al::span<const AngularP
     tmpflt.clear();
     impres.clear();
 
-    for(size_t i{0u};i < mCoeffs.size();++i)
+    for(size_t i{0u};i < mChannels.size();++i)
     {
         auto copy_arr = [](const double2 &in) noexcept -> float2
         { return float2{{static_cast<float>(in[0]), static_cast<float>(in[1])}}; };
-        std::transform(tmpres[i].cbegin(), tmpres[i].cend(), mCoeffs[i].begin(), copy_arr);
+        std::transform(tmpres[i].cbegin(), tmpres[i].cend(), mChannels[i].mCoeffs.begin(),
+            copy_arr);
     }
     tmpres.clear();
 
diff --git a/alc/hrtf.h b/alc/hrtf.h
index 5b709e30..7876265b 100644
--- a/alc/hrtf.h
+++ b/alc/hrtf.h
@@ -12,6 +12,8 @@
 #include "alspan.h"
 #include "ambidefs.h"
 #include "atomic.h"
+#include "bufferline.h"
+#include "filters/splitter.h"
 #include "intrusive_ptr.h"
 #include "vector.h"
 
@@ -77,12 +79,24 @@ struct AngularPoint {
     AzRadians Azim;
 };
 
+#define HRTF_DIRECT_DELAY 128
 struct DirectHrtfState {
+    struct ChannelData {
+        std::array<float,HRTF_DIRECT_DELAY> mDelay{};
+        BandSplitter mSplitter;
+        float mHfScale{};
+        alignas(16) HrirArray mCoeffs{};
+    };
+
+    std::array<float,HRTF_DIRECT_DELAY> mLeftDelay{};
+    std::array<float,HRTF_DIRECT_DELAY> mRightDelay{};
+    std::array<float,HRTF_DIRECT_DELAY+BUFFERSIZE> mTemp;
+
     /* HRTF filter state for dry buffer content */
     ALuint mIrSize{0};
-    al::FlexArray<HrirArray,16> mCoeffs;
+    al::FlexArray<ChannelData> mChannels;
 
-    DirectHrtfState(size_t numchans) : mCoeffs{numchans} { }
+    DirectHrtfState(size_t numchans) : mChannels{numchans} { }
     /**
      * Produces HRTF filter coefficients for decoding B-Format, given a set of
      * virtual speaker positions, a matching decoding matrix, and per-order
@@ -95,7 +109,7 @@ struct DirectHrtfState {
 
     static std::unique_ptr<DirectHrtfState> Create(size_t num_chans);
 
-    DEF_FAM_NEWDEL(DirectHrtfState, mCoeffs)
+    DEF_FAM_NEWDEL(DirectHrtfState, mChannels)
 };
 
 
diff --git a/alc/mixer/hrtfbase.h b/alc/mixer/hrtfbase.h
index 5bef1ddd..eea63efa 100644
--- a/alc/mixer/hrtfbase.h
+++ b/alc/mixer/hrtfbase.h
@@ -87,16 +87,74 @@ inline void MixDirectHrtfBase(FloatBufferLine &LeftOut, FloatBufferLine &RightOu
 
     const uint_fast32_t IrSize{State->mIrSize};
 
-    auto coeff_iter = State->mCoeffs.begin();
+    auto chan_iter = State->mChannels.begin();
     for(const FloatBufferLine &input : InSamples)
     {
-        const auto &Coeffs = *(coeff_iter++);
+        /* For dual-band processing, the signal needs extra scaling applied to
+         * the high frequency response. The band-splitter alone creates a
+         * frequency-dependent phase shift, which is not ideal. To counteract
+         * it, combine it with a backwards phase-shift.
+         */
+
+        /* Load the input signal backwards, into a temp buffer with delay
+         * padding. The delay serves to reduce the error caused by IIR filter's
+         * phase shift on a partial input.
+         */
+        al::span<float> tempbuf{State->mTemp.data(), HRTF_DIRECT_DELAY+BufferSize};
+        auto tmpiter = std::reverse_copy(input.begin(), input.begin()+BufferSize, tempbuf.begin());
+        std::copy(chan_iter->mDelay.cbegin(), chan_iter->mDelay.cend(), tmpiter);
+
+        /* Save the unfiltered newest input samples for next time. */
+        std::copy_n(tempbuf.begin(), HRTF_DIRECT_DELAY, chan_iter->mDelay.begin());
+
+        /* 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).
+         */
+        chan_iter->mSplitter.applyAllpass(tempbuf);
+        tempbuf = tempbuf.subspan<HRTF_DIRECT_DELAY>();
+        std::reverse(tempbuf.begin(), tempbuf.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.
+         */
+        chan_iter->mSplitter.applyHfScale(tempbuf, chan_iter->mHfScale);
+
+        /* Now apply the HRIR coefficients to this channel. */
+        const auto &Coeffs = chan_iter->mCoeffs;
+        ++chan_iter;
+
         for(size_t i{0u};i < BufferSize;++i)
         {
-            const float insample{input[i]};
+            const float insample{tempbuf[i]};
             ApplyCoeffs(AccumSamples+i, IrSize, Coeffs, insample, insample);
         }
     }
+
+    /* Apply a delay to the existing signal to align with the input delay. */
+    auto &ldelay = State->mLeftDelay;
+    auto &rdelay = State->mRightDelay;
+    if LIKELY(BufferSize >= HRTF_DIRECT_DELAY)
+    {
+        auto buffer_end = LeftOut.begin() + BufferSize;
+        auto delay_end = std::rotate(LeftOut.begin(), buffer_end - HRTF_DIRECT_DELAY, buffer_end);
+        std::swap_ranges(LeftOut.begin(), delay_end, ldelay.begin());
+
+        buffer_end = RightOut.begin() + BufferSize;
+        delay_end = std::rotate(RightOut.begin(), buffer_end - HRTF_DIRECT_DELAY, buffer_end);
+        std::swap_ranges(RightOut.begin(), delay_end, rdelay.begin());
+    }
+    else
+    {
+        auto buffer_end = LeftOut.begin() + BufferSize;
+        auto delay_start = std::swap_ranges(LeftOut.begin(), buffer_end, ldelay.begin());
+        std::rotate(ldelay.begin(), delay_start, ldelay.end());
+
+        buffer_end = RightOut.begin() + BufferSize;
+        delay_start = std::swap_ranges(RightOut.begin(), buffer_end, rdelay.begin());
+        std::rotate(rdelay.begin(), delay_start, rdelay.end());
+    }
     for(size_t i{0u};i < BufferSize;++i)
         LeftOut[i]  += AccumSamples[i][0];
     for(size_t i{0u};i < BufferSize;++i)