Alc/uhjfilter.c

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