Fix UHJ encoding/decoding factors

Classic B-Format uses scaling factors W=1, X=sqrt(2), Y=sqrt(2), and Z=sqrt(2),
which is +3dB louder than FuMa. The base factors are designed assuming classic
scaling, so encoding a 0dBFS FuMa signal without accounting for this would
result in the UHJ signal peaking at about -3dBFS. Similarly, decoding UHJ to
FuMa B-Format would be +3dB louder than intended.

So encoding needs to implicitly boost the signal by +3dB, and decoding needs to
attenuate by -3dB.

3 files changed, 37 insertions, 17 deletions

diff --git a/core/uhjfilter.cpp b/core/uhjfilter.cpp
index fad04d20..5a236d38 100644
--- a/core/uhjfilter.cpp
+++ b/core/uhjfilter.cpp
@@ -43,6 +43,9 @@ const PhaseShifterT<UhjEncoder::sFilterDelay*2> PShift{};
 void UhjEncoder::encode(const FloatBufferSpan LeftOut, const FloatBufferSpan RightOut,
     const FloatBufferLine *InSamples, const size_t SamplesToDo)
 {
+    /* Given FuMa input, a +3dB boost is needed for the expected levels. */
+    constexpr float sqrt2{1.41421356237f};
+
     ASSUME(SamplesToDo > 0);
 
     float *RESTRICT left{al::assume_aligned<16>(LeftOut.data())};
@@ -60,14 +63,14 @@ void UhjEncoder::encode(const FloatBufferSpan LeftOut, const FloatBufferSpan Rig
     auto miditer = mS.begin() + sFilterDelay;
     std::transform(winput, winput+SamplesToDo, xinput, miditer,
         [](const float w, const float x) noexcept -> float
-        { return 0.9396926f*w + 0.1855740f*x; });
+        { return 0.9396926f*sqrt2*w + 0.1855740f*sqrt2*x; });
     for(size_t i{0};i < SamplesToDo;++i,++miditer)
         *miditer += left[i] + right[i];
 
     /* D = 0.6554516*Y */
     auto sideiter = mD.begin() + sFilterDelay;
     std::transform(yinput, yinput+SamplesToDo, sideiter,
-        [](const float y) noexcept -> float { return 0.6554516f*y; });
+        [](const float y) noexcept -> float { return 0.6554516f*sqrt2*y; });
     for(size_t i{0};i < SamplesToDo;++i,++sideiter)
         *sideiter += left[i] - right[i];
 
@@ -75,7 +78,7 @@ void UhjEncoder::encode(const FloatBufferSpan LeftOut, const FloatBufferSpan Rig
     auto tmpiter = std::copy(mWXHistory.cbegin(), mWXHistory.cend(), mTemp.begin());
     std::transform(winput, winput+SamplesToDo, xinput, tmpiter,
         [](const float w, const float x) noexcept -> float
-        { return -0.3420201f*w + 0.5098604f*x; });
+        { return -0.3420201f*sqrt2*w + 0.5098604f*sqrt2*x; });
     std::copy_n(mTemp.cbegin()+SamplesToDo, mWXHistory.size(), mWXHistory.begin());
     PShift.processAccum({mD.data(), SamplesToDo}, mTemp.data());
 
@@ -110,23 +113,32 @@ void UhjEncoder::encode(const FloatBufferSpan LeftOut, const FloatBufferSpan Rig
 void UhjDecoder::decode(const al::span<BufferLine> samples, const size_t offset,
     const size_t samplesToDo, const size_t forwardSamples)
 {
+    /* A -3dB attenuation is needed for FuMa output. */
+    constexpr float sqrt1_2{0.707106781187f};
+
     ASSUME(samplesToDo > 0);
 
-    /* S = Left + Right */
-    for(size_t i{0};i < samplesToDo+sFilterDelay;++i)
-        mS[i] = samples[0][offset+i] + samples[1][offset+i];
+    {
+        const float *RESTRICT left{al::assume_aligned<16>(samples[0].data() + offset)};
+        const float *RESTRICT right{al::assume_aligned<16>(samples[1].data() + offset)};
+        const float *RESTRICT t{al::assume_aligned<16>(samples[2].data() + offset)};
 
-    /* D = Left - Right */
-    for(size_t i{0};i < samplesToDo+sFilterDelay;++i)
-        mD[i] = samples[0][offset+i] - samples[1][offset+i];
+        /* S = Left + Right */
+        for(size_t i{0};i < samplesToDo+sFilterDelay;++i)
+            mS[i] = (left[i] + right[i]) * sqrt1_2;
 
-    /* T */
-    for(size_t i{0};i < samplesToDo+sFilterDelay;++i)
-        mT[i] = samples[2][offset+i];
+        /* D = Left - Right */
+        for(size_t i{0};i < samplesToDo+sFilterDelay;++i)
+            mD[i] = (left[i] - right[i]) * sqrt1_2;
+
+        /* T */
+        for(size_t i{0};i < samplesToDo+sFilterDelay;++i)
+            mT[i] = t[i] * sqrt1_2;
+    }
 
-    float *RESTRICT woutput{samples[0].data() + offset};
-    float *RESTRICT xoutput{samples[1].data() + offset};
-    float *RESTRICT youtput{samples[2].data() + offset};
+    float *RESTRICT woutput{al::assume_aligned<16>(samples[0].data() + offset)};
+    float *RESTRICT xoutput{al::assume_aligned<16>(samples[1].data() + offset)};
+    float *RESTRICT youtput{al::assume_aligned<16>(samples[2].data() + offset)};
 
     /* Precompute j(0.828347*D + 0.767835*T) and store in xoutput. */
     auto tmpiter = std::copy(mDTHistory.cbegin(), mDTHistory.cend(), mTemp.begin());
@@ -157,6 +169,6 @@ void UhjDecoder::decode(const al::span<BufferLine> samples, const size_t offset,
         float *RESTRICT zoutput{samples[3].data() + offset};
         /* Z = 1.023332*Q */
         for(size_t i{0};i < samplesToDo;++i)
-            zoutput[i] = 1.023332f*zoutput[i];
+            zoutput[i] = 1.023332f*sqrt1_2*zoutput[i];
     }
 }
diff --git a/core/uhjfilter.h b/core/uhjfilter.h
index c28b9510..897de2f8 100644
--- a/core/uhjfilter.h
+++ b/core/uhjfilter.h
@@ -49,6 +49,12 @@ struct UhjDecoder {
 
     alignas(16) std::array<float,BufferLineSize+MaxResamplerEdge + sFilterDelay*2> mTemp{};
 
+    /**
+     * Decodes a 3- or 4-channel UHJ signal into a B-Format signal with FuMa
+     * channel ordering and scaling. For 3-channel, the 3rd channel may be
+     * attenuated by 'n', where 0 <= n <= 1. So 2-channel UHJ can be decoded by
+     * leaving the 3rd channel silent (n=0).
+     */
     void decode(const al::span<BufferLine> samples, const size_t offset, const size_t samplesToDo,
         const size_t forwardSamples);
 
diff --git a/utils/uhjdecoder.cpp b/utils/uhjdecoder.cpp
index f9db4895..cda16a40 100644
--- a/utils/uhjdecoder.cpp
+++ b/utils/uhjdecoder.cpp
@@ -492,8 +492,10 @@ int main(int argc, char **argv)
                 decoder->decode2(inmem.get(), decmem, got);
             for(size_t i{0};i < got;++i)
             {
+                /* Attenuate by -3dB for FuMa output levels. */
+                constexpr float sqrt1_2{0.707106781187f};
                 for(size_t j{0};j < outchans;++j)
-                    outmem[i*outchans + j] = f32AsLEBytes(decmem[j][i]);
+                    outmem[i*outchans + j] = f32AsLEBytes(decmem[j][i] * sqrt1_2);
             }
 
             size_t wrote{fwrite(outmem.get(), sizeof(byte4)*outchans, got, outfile.get())};