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  128
  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 static void allpass_process(AllPassState *state, ALfloat *restrict dst, const ALfloat *restrict src, const ALfloat aa, ALsizei todo)
  20 {
  21     ALsizei i;
  22
  23     if(LIKELY(todo > 1))
  24     {
  25         ALfloat x0 = state->x[0];
  26         ALfloat x1 = state->x[1];
  27         ALfloat y0 = state->y[0];
  28         ALfloat y1 = state->y[1];
  29
  30         for(i = 0;i < todo;i++)
  31         {
  32             dst[i] = aa*(src[i] + y1) - x1;
  33             y1 = y0; y0 = dst[i];
  34             x1 = x0; x0 = src[i];
  35         }
  36
  37         state->x[0] = x0;
  38         state->x[1] = x1;
  39         state->y[0] = y0;
  40         state->y[1] = y1;
  41     }
  42     else if(todo == 1)
  43     {
  44         dst[0] = aa*(src[0] + state->y[1]) - state->x[1];
  45         state->x[1] = state->x[0];
  46         state->x[0] = src[0];
  47         state->y[1] = state->y[0];
  48         state->y[0] = dst[0];
  49     }
  50 }
  51
  52
  53 /* NOTE: There seems to be a bit of an inconsistency in how this encoding is
  54  * supposed to work. Some references, such as
  55  *
  56  * http://members.tripod.com/martin_leese/Ambisonic/UHJ_file_format.html
  57  *
  58  * specify a pre-scaling of sqrt(2) on the W channel input, while other
  59  * references, such as
  60  *
  61  * https://en.wikipedia.org/wiki/Ambisonic_UHJ_format#Encoding.5B1.5D
  62  * and
  63  * https://wiki.xiph.org/Ambisonics#UHJ_format
  64  *
  65  * do not. The sqrt(2) scaling is in line with B-Format decoder coefficients
  66  * which include such a scaling for the W channel input, however the original
  67  * source for this equation is a 1985 paper by Michael Gerzon, which does not
  68  * apparently include the scaling. Applying the extra scaling creates a louder
  69  * result with a narrower stereo image compared to not scaling, and I don't
  70  * know which is the intended result.
  71  */
  72
  73 void EncodeUhj2(Uhj2Encoder *enc, ALfloat *restrict LeftOut, ALfloat *restrict RightOut, ALfloat (*restrict InSamples)[BUFFERSIZE], ALsizei SamplesToDo)
  74 {
  75     ALfloat D[MAX_UPDATE_SAMPLES], S[MAX_UPDATE_SAMPLES];
  76     ALfloat temp[2][MAX_UPDATE_SAMPLES];
  77     ALsizei base, i;
  78
  79     for(base = 0;base < SamplesToDo;)
  80     {
  81         ALsizei todo = mini(SamplesToDo - base, MAX_UPDATE_SAMPLES);
  82
  83         /* D = 0.6554516*Y */
  84         for(i = 0;i < todo;i++)
  85             temp[0][i] = 0.6554516f*InSamples[2][base+i];
  86         allpass_process(&enc->Filter1_Y[0], temp[1], temp[0],
  87                         Filter1Coeff[0]*Filter1Coeff[0], todo);
  88         allpass_process(&enc->Filter1_Y[1], temp[0], temp[1],
  89                         Filter1Coeff[1]*Filter1Coeff[1], todo);
  90         allpass_process(&enc->Filter1_Y[2], temp[1], temp[0],
  91                         Filter1Coeff[2]*Filter1Coeff[2], todo);
  92         /* NOTE: Filter1 requires a 1 sample delay for the final output, so
  93          * take the last processed sample from the previous run as the first
  94          * output sample.
  95          */
  96         D[0] = enc->Filter1_Y[3].y[0];
  97         allpass_process(&enc->Filter1_Y[3], temp[0], temp[1],
  98                         Filter1Coeff[3]*Filter1Coeff[3], todo);
  99         for(i = 1;i < todo;i++)
 100             D[i] = temp[0][i-1];
 101
 102         /* D += j(-0.3420201*W + 0.5098604*X) */
 103         for(i = 0;i < todo;i++)
 104             temp[0][i] = -0.3420201f*InSamples[0][base+i] +
 105                           0.5098604f*InSamples[1][base+i];
 106         allpass_process(&enc->Filter2_WX[0], temp[1], temp[0],
 107                         Filter2Coeff[0]*Filter2Coeff[0], todo);
 108         allpass_process(&enc->Filter2_WX[1], temp[0], temp[1],
 109                         Filter2Coeff[1]*Filter2Coeff[1], todo);
 110         allpass_process(&enc->Filter2_WX[2], temp[1], temp[0],
 111                         Filter2Coeff[2]*Filter2Coeff[2], todo);
 112         allpass_process(&enc->Filter2_WX[3], temp[0], temp[1],
 113                         Filter2Coeff[3]*Filter2Coeff[3], todo);
 114         for(i = 0;i < todo;i++)
 115             D[i] += temp[0][i];
 116
 117         /* S = 0.9396926*W + 0.1855740*X */
 118         for(i = 0;i < todo;i++)
 119             temp[0][i] = 0.9396926f*InSamples[0][base+i] +
 120                          0.1855740f*InSamples[1][base+i];
 121         allpass_process(&enc->Filter1_WX[0], temp[1], temp[0],
 122                         Filter1Coeff[0]*Filter1Coeff[0], todo);
 123         allpass_process(&enc->Filter1_WX[1], temp[0], temp[1],
 124                         Filter1Coeff[1]*Filter1Coeff[1], todo);
 125         allpass_process(&enc->Filter1_WX[2], temp[1], temp[0],
 126                         Filter1Coeff[2]*Filter1Coeff[2], todo);
 127         S[0] = enc->Filter1_WX[3].y[0];
 128         allpass_process(&enc->Filter1_WX[3], temp[0], temp[1],
 129                         Filter1Coeff[3]*Filter1Coeff[3], todo);
 130         for(i = 1;i < todo;i++)
 131             S[i] = temp[0][i-1];
 132
 133         /* Left = (S + D)/2.0 */
 134         for(i = 0;i < todo;i++)
 135             *(LeftOut++) += (S[i] + D[i]) * 0.5f;
 136         /* Right = (S - D)/2.0 */
 137         for(i = 0;i < todo;i++)
 138             *(RightOut++) += (S[i] - D[i]) * 0.5f;
 139
 140         base += todo;
 141     }
 142 }