alc/uhjfilter.cpp

   1
   2 #include "config.h"
   3
   4 #include "uhjfilter.h"
   5
   6 #include <algorithm>
   7 #include <iterator>
   8
   9 #include "AL/al.h"
  10
  11 #include "alnumeric.h"
  12 #include "opthelpers.h"
  13
  14
  15 namespace {
  16
  17 /* This is the maximum number of samples processed for each inner loop
  18  * iteration. */
  19 #define MAX_UPDATE_SAMPLES  128
  20
  21
  22 constexpr ALfloat Filter1CoeffSqr[4] = {
  23     0.479400865589f, 0.876218493539f, 0.976597589508f, 0.997499255936f
  24 };
  25 constexpr ALfloat Filter2CoeffSqr[4] = {
  26     0.161758498368f, 0.733028932341f, 0.945349700329f, 0.990599156685f
  27 };
  28
  29 void allpass_process(AllPassState *state, ALfloat *dst, const ALfloat *src, const ALfloat aa,
  30     const size_t todo)
  31 {
  32     ALfloat z1{state->z[0]};
  33     ALfloat z2{state->z[1]};
  34     auto proc_sample = [aa,&z1,&z2](const ALfloat input) noexcept -> ALfloat
  35     {
  36         const ALfloat output{input*aa + z1};
  37         z1 = z2; z2 = output*aa - input;
  38         return output;
  39     };
  40     std::transform(src, src+todo, dst, proc_sample);
  41     state->z[0] = z1;
  42     state->z[1] = z2;
  43 }
  44
  45 } // namespace
  46
  47
  48 /* NOTE: There seems to be a bit of an inconsistency in how this encoding is
  49  * supposed to work. Some references, such as
  50  *
  51  * http://members.tripod.com/martin_leese/Ambisonic/UHJ_file_format.html
  52  *
  53  * specify a pre-scaling of sqrt(2) on the W channel input, while other
  54  * references, such as
  55  *
  56  * https://en.wikipedia.org/wiki/Ambisonic_UHJ_format#Encoding.5B1.5D
  57  * and
  58  * https://wiki.xiph.org/Ambisonics#UHJ_format
  59  *
  60  * do not. The sqrt(2) scaling is in line with B-Format decoder coefficients
  61  * which include such a scaling for the W channel input, however the original
  62  * source for this equation is a 1985 paper by Michael Gerzon, which does not
  63  * apparently include the scaling. Applying the extra scaling creates a louder
  64  * result with a narrower stereo image compared to not scaling, and I don't
  65  * know which is the intended result.
  66  */
  67
  68 void Uhj2Encoder::encode(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
  69     FloatBufferLine *InSamples, const size_t SamplesToDo)
  70 {
  71     alignas(16) ALfloat D[MAX_UPDATE_SAMPLES], S[MAX_UPDATE_SAMPLES];
  72     alignas(16) ALfloat temp[MAX_UPDATE_SAMPLES];
  73
  74     ASSUME(SamplesToDo > 0);
  75
  76     auto winput = InSamples[0].cbegin();
  77     auto xinput = InSamples[1].cbegin();
  78     auto yinput = InSamples[2].cbegin();
  79     for(size_t base{0};base < SamplesToDo;)
  80     {
  81         const size_t todo{minz(SamplesToDo - base, MAX_UPDATE_SAMPLES)};
  82         ASSUME(todo > 0);
  83
  84         /* D = 0.6554516*Y */
  85         std::transform(yinput, yinput+todo, std::begin(temp),
  86             [](const float y) noexcept -> float { return 0.6554516f*y; });
  87         allpass_process(&mFilter1_Y[0], temp, temp, Filter1CoeffSqr[0], todo);
  88         allpass_process(&mFilter1_Y[1], temp, temp, Filter1CoeffSqr[1], todo);
  89         allpass_process(&mFilter1_Y[2], temp, temp, Filter1CoeffSqr[2], todo);
  90         allpass_process(&mFilter1_Y[3], temp, temp, Filter1CoeffSqr[3], todo);
  91         /* NOTE: Filter1 requires a 1 sample delay for the final output, so
  92          * take the last processed sample from the previous run as the first
  93          * output sample.
  94          */
  95         D[0] = mLastY;
  96         for(size_t i{1};i < todo;i++)
  97             D[i] = temp[i-1];
  98         mLastY = temp[todo-1];
  99
 100         /* D += j(-0.3420201*W + 0.5098604*X) */
 101         std::transform(winput, winput+todo, xinput, std::begin(temp),
 102             [](const float w, const float x) noexcept -> float
 103             { return -0.3420201f*w + 0.5098604f*x; });
 104         allpass_process(&mFilter2_WX[0], temp, temp, Filter2CoeffSqr[0], todo);
 105         allpass_process(&mFilter2_WX[1], temp, temp, Filter2CoeffSqr[1], todo);
 106         allpass_process(&mFilter2_WX[2], temp, temp, Filter2CoeffSqr[2], todo);
 107         allpass_process(&mFilter2_WX[3], temp, temp, Filter2CoeffSqr[3], todo);
 108         for(size_t i{0};i < todo;i++)
 109             D[i] += temp[i];
 110
 111         /* S = 0.9396926*W + 0.1855740*X */
 112         std::transform(winput, winput+todo, xinput, std::begin(temp),
 113             [](const float w, const float x) noexcept -> float
 114             { return 0.9396926f*w + 0.1855740f*x; });
 115         allpass_process(&mFilter1_WX[0], temp, temp, Filter1CoeffSqr[0], todo);
 116         allpass_process(&mFilter1_WX[1], temp, temp, Filter1CoeffSqr[1], todo);
 117         allpass_process(&mFilter1_WX[2], temp, temp, Filter1CoeffSqr[2], todo);
 118         allpass_process(&mFilter1_WX[3], temp, temp, Filter1CoeffSqr[3], todo);
 119         S[0] = mLastWX;
 120         for(size_t i{1};i < todo;i++)
 121             S[i] = temp[i-1];
 122         mLastWX = temp[todo-1];
 123
 124         /* Left = (S + D)/2.0 */
 125         ALfloat *RESTRICT left = al::assume_aligned<16>(LeftOut.data()+base);
 126         for(size_t i{0};i < todo;i++)
 127             left[i] += (S[i] + D[i]) * 0.5f;
 128         /* Right = (S - D)/2.0 */
 129         ALfloat *RESTRICT right = al::assume_aligned<16>(RightOut.data()+base);
 130         for(size_t i{0};i < todo;i++)
 131             right[i] += (S[i] - D[i]) * 0.5f;
 132
 133         winput += todo;
 134         xinput += todo;
 135         yinput += todo;
 136         base += todo;
 137     }
 138 }