Add a function to encode 2-channel UHJ from B-Format

2 files changed, 140 insertions, 0 deletions

diff --git a/Alc/uhjfilter.c b/Alc/uhjfilter.c
new file mode 100644
index 00000000..63a88aa7
--- /dev/null
+++ b/Alc/uhjfilter.c
@@ -0,0 +1,91 @@
+
+#include "config.h"
+
+#include "alu.h"
+#include "uhjfilter.h"
+
+/* This is the maximum number of samples processed for each inner loop
+ * iteration. */
+#define MAX_UPDATE_SAMPLES  256
+
+
+static const ALfloat Filter1Coeff[4] = {
+    0.6923878f, 0.9360654322959f, 0.9882295226860f, 0.9987488452737f
+};
+static const ALfloat Filter2Coeff[4] = {
+    0.4021921162426f, 0.8561710882420f, 0.9722909545651f, 0.9952884791278f
+};
+
+void EncodeUhj2(Uhj2Encoder *enc, ALfloat (*restrict OutBuffer)[BUFFERSIZE], ALfloat (*restrict InSamples)[BUFFERSIZE], ALuint SamplesToDo)
+{
+    ALuint base, i, c;
+
+    for(base = 0;base < SamplesToDo;)
+    {
+        ALfloat D[MAX_UPDATE_SAMPLES], S;
+        ALuint todo = minu(SamplesToDo - base, MAX_UPDATE_SAMPLES);
+
+        /* D = 0.6554516*Y */
+        for(i = 0;i < todo;i++)
+        {
+            ALfloat in = 0.6554516f*InSamples[2][base+i];
+            for(c = 0;c < 4;c++)
+            {
+                ALfloat aa = Filter1Coeff[c]*Filter1Coeff[c];
+                ALfloat out = aa*(in + enc->Filter1_Y[c].y[1]) - enc->Filter1_Y[c].x[1];
+                enc->Filter1_Y[c].x[1] = enc->Filter1_Y[c].x[0];
+                enc->Filter1_Y[c].x[0] = in;
+                enc->Filter1_Y[c].y[1] = enc->Filter1_Y[c].y[0];
+                enc->Filter1_Y[c].y[0] = out;
+                in = out;
+            }
+            /* NOTE: Filter1 requires a 1 sample delay for the base output, so
+             * take the sample before the last for output.
+             */
+            D[i] = enc->Filter1_Y[3].y[1];
+        }
+
+        /* D += j(-0.3420201*W' + 0.5098604*X) */
+        for(i = 0;i < todo;i++)
+        {
+            ALfloat in = -0.3420201f*1.414213562f*InSamples[0][base+i] +
+                          0.5098604f*InSamples[1][base+i];
+            for(c = 0;c < 4;c++)
+            {
+                ALfloat aa = Filter2Coeff[c]*Filter2Coeff[c];
+                ALfloat out = aa*(in + enc->Filter2_WX[c].y[1]) - enc->Filter2_WX[c].x[1];
+                enc->Filter2_WX[c].x[1] = enc->Filter2_WX[c].x[0];
+                enc->Filter2_WX[c].x[0] = in;
+                enc->Filter2_WX[c].y[1] = enc->Filter2_WX[c].y[0];
+                enc->Filter2_WX[c].y[0] = out;
+                in = out;
+            }
+            D[i] += enc->Filter2_WX[3].y[0];
+        }
+
+        /* S = 0.9396926*W' + 0.1855740*X
+         * Left = (S + D)/2.0
+         * Right = (S - D)/2.0
+         */
+        for(i = 0;i < todo;i++)
+        {
+            ALfloat in = 0.9396926f*1.414213562f*InSamples[0][base+i] +
+                         0.1855740f*InSamples[1][base+i];
+            for(c = 0;c < 4;c++)
+            {
+                ALfloat aa = Filter1Coeff[c]*Filter1Coeff[c];
+                ALfloat out = aa*(in + enc->Filter1_WX[c].y[1]) - enc->Filter1_WX[c].x[1];
+                enc->Filter1_WX[c].x[1] = enc->Filter1_WX[c].x[0];
+                enc->Filter1_WX[c].x[0] = in;
+                enc->Filter1_WX[c].y[1] = enc->Filter1_WX[c].y[0];
+                enc->Filter1_WX[c].y[0] = out;
+                in = out;
+            }
+            S = enc->Filter1_WX[3].y[1];
+            OutBuffer[0][base + i] += (S + D[i]) * 0.5f;
+            OutBuffer[1][base + i] += (S - D[i]) * 0.5f;
+        }
+
+        base += todo;
+    }
+}
diff --git a/Alc/uhjfilter.h b/Alc/uhjfilter.h
new file mode 100644
index 00000000..5238e202
--- /dev/null
+++ b/Alc/uhjfilter.h
@@ -0,0 +1,49 @@
+#ifndef UHJFILTER_H
+#define UHJFILTER_H
+
+#include "AL/al.h"
+
+#include "alMain.h"
+
+typedef struct AllPassState {
+    ALfloat x[2]; /* Last two input samples */
+    ALfloat y[2]; /* Last two output samples */
+} AllPassState;
+
+/* Encoding 2-channel UHJ from B-Format is done as:
+ *
+ * W' = W * sqrt(2)
+ * S = 0.9396926*W' + 0.1855740*X
+ * D = j(-0.3420201*W' + 0.5098604*X) + 0.6554516*Y
+ *
+ * Left = (S + D)/2.0
+ * Right = (S - D)/2.0
+ *
+ * where j is a wide-band +90 degree phase shift.
+ *
+ * The phase shift is done using a Hilbert transform, described here:
+ * https://web.archive.org/web/20060708031958/http://www.biochem.oulu.fi/~oniemita/dsp/hilbert/
+ * It works using 2 sets of 4 chained filters. The first filter chain produces
+ * a phase shift of varying magnitude over a wide range of frequencies, while
+ * the second filter chain produces a phase shift 90 degrees ahead of the
+ * first over the same range.
+ *
+ * Combining these two stages requires the use of three filter chains. S-
+ * channel output uses a Filter1 chain on the W and X channel mix, while the D-
+ * channel output uses a Filter1 chain on the Y channel plus a Filter2 chain on
+ * the W and X channel mix. This results in the W and X input mix on the D-
+ * channel output having the required +90 degree phase shift relative to the
+ * other inputs.
+ */
+
+typedef struct Uhj2Encoder {
+    AllPassState Filter1_WX[4];
+    AllPassState Filter1_Y[4];
+    AllPassState Filter2_WX[4];
+} Uhj2Encoder;
+
+/* Encodes a 2-channel UHJ (stereo-compatible) signal from a B-Format input
+ * signal. */
+void EncodeUhj2(Uhj2Encoder *enc, ALfloat (*restrict OutBuffer)[BUFFERSIZE], ALfloat (*restrict InSamples)[BUFFERSIZE], ALuint SamplesToDo);
+
+#endif /* UHJFILTER_H */