Alc/uhjfilter.cpp

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