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