7 /* This is the maximum number of samples processed for each inner loop
9 #define MAX_UPDATE_SAMPLES 128
12 static const ALfloat Filter1CoeffSqr
[4] = {
13 0.479400865589f
, 0.876218493539f
, 0.976597589508f
, 0.997499255936f
15 static const ALfloat Filter2CoeffSqr
[4] = {
16 0.161758498368f
, 0.733028932341f
, 0.945349700329f
, 0.990599156685f
19 static void allpass_process(AllPassState
*state
, ALfloat
*restrict dst
, const ALfloat
*restrict src
, const ALfloat aa
, ALsizei todo
)
21 ALfloat z1
= state
->z
[0];
22 ALfloat z2
= state
->z
[1];
25 for(i
= 0;i
< todo
;i
++)
27 ALfloat input
= src
[i
];
28 ALfloat output
= input
*aa
+ z1
;
29 z1
= z2
; z2
= output
*aa
- input
;
38 /* NOTE: There seems to be a bit of an inconsistency in how this encoding is
39 * supposed to work. Some references, such as
41 * http://members.tripod.com/martin_leese/Ambisonic/UHJ_file_format.html
43 * specify a pre-scaling of sqrt(2) on the W channel input, while other
46 * https://en.wikipedia.org/wiki/Ambisonic_UHJ_format#Encoding.5B1.5D
48 * https://wiki.xiph.org/Ambisonics#UHJ_format
50 * do not. The sqrt(2) scaling is in line with B-Format decoder coefficients
51 * which include such a scaling for the W channel input, however the original
52 * source for this equation is a 1985 paper by Michael Gerzon, which does not
53 * apparently include the scaling. Applying the extra scaling creates a louder
54 * result with a narrower stereo image compared to not scaling, and I don't
55 * know which is the intended result.
58 void EncodeUhj2(Uhj2Encoder
*enc
, ALfloat
*restrict LeftOut
, ALfloat
*restrict RightOut
, ALfloat (*restrict InSamples
)[BUFFERSIZE
], ALsizei SamplesToDo
)
60 ALfloat D
[MAX_UPDATE_SAMPLES
], S
[MAX_UPDATE_SAMPLES
];
61 ALfloat temp
[2][MAX_UPDATE_SAMPLES
];
64 ASSUME(SamplesToDo
> 0);
66 for(base
= 0;base
< SamplesToDo
;)
68 ALsizei todo
= mini(SamplesToDo
- base
, MAX_UPDATE_SAMPLES
);
72 for(i
= 0;i
< todo
;i
++)
73 temp
[0][i
] = 0.6554516f
*InSamples
[2][base
+i
];
74 allpass_process(&enc
->Filter1_Y
[0], temp
[1], temp
[0], Filter1CoeffSqr
[0], todo
);
75 allpass_process(&enc
->Filter1_Y
[1], temp
[0], temp
[1], Filter1CoeffSqr
[1], todo
);
76 allpass_process(&enc
->Filter1_Y
[2], temp
[1], temp
[0], Filter1CoeffSqr
[2], todo
);
77 allpass_process(&enc
->Filter1_Y
[3], temp
[0], temp
[1], Filter1CoeffSqr
[3], todo
);
78 /* NOTE: Filter1 requires a 1 sample delay for the final output, so
79 * take the last processed sample from the previous run as the first
83 for(i
= 1;i
< todo
;i
++)
85 enc
->LastY
= temp
[0][i
-1];
87 /* D += j(-0.3420201*W + 0.5098604*X) */
88 for(i
= 0;i
< todo
;i
++)
89 temp
[0][i
] = -0.3420201f
*InSamples
[0][base
+i
] +
90 0.5098604f
*InSamples
[1][base
+i
];
91 allpass_process(&enc
->Filter2_WX
[0], temp
[1], temp
[0], Filter2CoeffSqr
[0], todo
);
92 allpass_process(&enc
->Filter2_WX
[1], temp
[0], temp
[1], Filter2CoeffSqr
[1], todo
);
93 allpass_process(&enc
->Filter2_WX
[2], temp
[1], temp
[0], Filter2CoeffSqr
[2], todo
);
94 allpass_process(&enc
->Filter2_WX
[3], temp
[0], temp
[1], Filter2CoeffSqr
[3], todo
);
95 for(i
= 0;i
< todo
;i
++)
98 /* S = 0.9396926*W + 0.1855740*X */
99 for(i
= 0;i
< todo
;i
++)
100 temp
[0][i
] = 0.9396926f
*InSamples
[0][base
+i
] +
101 0.1855740f
*InSamples
[1][base
+i
];
102 allpass_process(&enc
->Filter1_WX
[0], temp
[1], temp
[0], Filter1CoeffSqr
[0], todo
);
103 allpass_process(&enc
->Filter1_WX
[1], temp
[0], temp
[1], Filter1CoeffSqr
[1], todo
);
104 allpass_process(&enc
->Filter1_WX
[2], temp
[1], temp
[0], Filter1CoeffSqr
[2], todo
);
105 allpass_process(&enc
->Filter1_WX
[3], temp
[0], temp
[1], Filter1CoeffSqr
[3], todo
);
107 for(i
= 1;i
< todo
;i
++)
109 enc
->LastWX
= temp
[0][i
-1];
111 /* Left = (S + D)/2.0 */
112 for(i
= 0;i
< todo
;i
++)
113 *(LeftOut
++) += (S
[i
] + D
[i
]) * 0.5f
;
114 /* Right = (S - D)/2.0 */
115 for(i
= 0;i
< todo
;i
++)
116 *(RightOut
++) += (S
[i
] - D
[i
]) * 0.5f
;