1// Copyright (C) 2022 The Qt Company Ltd.
2// SPDX-License-Identifier: LicenseRef-Qt-Commercial OR LGPL-3.0-only OR GPL-3.0-only
3
4#include "qambisonicdecoder_p.h"
5
6#include <QtCore/qdebug.h>
7#include <QtSpatialAudio/private/qambisonicdecoderdata_p.h>
8
9#include <cmath>
10
11QT_BEGIN_NAMESPACE
12
13// Ambisonic decoding is described in detail in https://ambisonics.dreamhosters.com/BLaH3.pdf.
14// We're using a phase matched band splitting filter to split the ambisonic signal into a low
15// and high frequency component and apply matrix conversions to those components individually
16// as described in the document.
17//
18// We are currently not using a near field compensation filter, something that could potentially
19// improve sound quality further.
20//
21// For mono and stereo decoding, we use a simpler algorithm to avoid artificially dampening signals
22// coming from the back, as we do not have any speakers in that direction and the calculations
23// through matlab would give us audible 'holes'.
24
25struct QAmbisonicDecoderData
26{
27 QAudioFormat::ChannelConfig config;
28 const float *lf[3];
29 const float *hf[3];
30 const float *reverb;
31};
32
33constexpr float reverb_x_0[] = {
34 1.f, 0.f, // L
35 0.f, 1.f, // R
36 .7f, .7f, // C
37 1.f, 0.f, // Ls
38 0.f, 1.f, // Rs
39 1.f, 0.f, // Lb
40 0.f, 1.f, // Rb
41};
42
43constexpr float reverb_x_1[] = {
44 1.f, 0.f, // L
45 0.f, 1.f, // R
46 .7f, .7f, // C
47 .0f, .0f, // LFE
48 1.f, 0.f, // Ls
49 0.f, 1.f, // Rs
50 1.f, 0.f, // Lb
51 0.f, 1.f, // Rb
52};
53
54static constexpr QAmbisonicDecoderData decoderMap[] = {
55 { .config: QAudioFormat::ChannelConfigSurround5Dot0,
56 .lf: { decoderMatrix_5dot0_1_lf, decoderMatrix_5dot0_2_lf, decoderMatrix_5dot0_3_lf },
57 .hf: { decoderMatrix_5dot0_1_hf, decoderMatrix_5dot0_2_hf, decoderMatrix_5dot0_3_hf },
58 .reverb: reverb_x_0 },
59 { .config: QAudioFormat::ChannelConfigSurround5Dot1,
60 .lf: { decoderMatrix_5dot1_1_lf, decoderMatrix_5dot1_2_lf, decoderMatrix_5dot1_3_lf },
61 .hf: { decoderMatrix_5dot1_1_hf, decoderMatrix_5dot1_2_hf, decoderMatrix_5dot1_3_hf },
62 .reverb: reverb_x_1 },
63 { .config: QAudioFormat::ChannelConfigSurround7Dot0,
64 .lf: { decoderMatrix_7dot0_1_lf, decoderMatrix_7dot0_2_lf, decoderMatrix_7dot0_3_lf },
65 .hf: { decoderMatrix_7dot0_1_hf, decoderMatrix_7dot0_2_hf, decoderMatrix_7dot0_3_hf },
66 .reverb: reverb_x_0 },
67 { .config: QAudioFormat::ChannelConfigSurround7Dot1,
68 .lf: { decoderMatrix_7dot1_1_lf, decoderMatrix_7dot1_2_lf, decoderMatrix_7dot1_3_lf },
69 .hf: { decoderMatrix_7dot1_1_hf, decoderMatrix_7dot1_2_hf, decoderMatrix_7dot1_3_hf },
70 .reverb: reverb_x_1 }
71};
72
73// Implements a split second order IIR filter
74// The audio data is split into a phase synced low and high frequency part
75// This allows us to apply different factors to both parts for better sound
76// localization when converting from ambisonic formats
77//
78// Details are described in https://ambisonics.dreamhosters.com/BLaH3.pdf, Appendix A.2.
79class QAmbisonicDecoderFilter
80{
81public:
82 QAmbisonicDecoderFilter() = default;
83 void configure(float sampleRate, float cutoffFrequency = 380)
84 {
85 double k = tan(M_PI*cutoffFrequency/sampleRate);
86 a1 = float(2.*(k*k - 1.)/(k*k + 2*k + 1.));
87 a2 = float((k*k - 2*k + 1.)/(k*k + 2*k + 1.));
88
89 b0_lf = float(k*k/(k*k + 2*k + 1));
90 b1_lf = 2.f*b0_lf;
91
92 b0_hf = float(1./(k*k + 2*k + 1));
93 b1_hf = -2.f*b0_hf;
94 }
95
96 struct Output
97 {
98 float lf;
99 float hf;
100 };
101
102 Output next(float x)
103 {
104 float r_lf = x*b0_lf +
105 prevX[0]*b1_lf +
106 prevX[1]*b0_lf -
107 prevR_lf[0]*a1 -
108 prevR_lf[1]*a2;
109 float r_hf = x*b0_hf +
110 prevX[0]*b1_hf +
111 prevX[1]*b0_hf -
112 prevR_hf[0]*a1 -
113 prevR_hf[1]*a2;
114 prevX[1] = prevX[0];
115 prevX[0] = x;
116 prevR_lf[1] = prevR_lf[0];
117 prevR_lf[0] = r_lf;
118 prevR_hf[1] = prevR_hf[0];
119 prevR_hf[0] = r_hf;
120 return { .lf: r_lf, .hf: r_hf };
121 }
122
123private:
124 float a1 = 0.;
125 float a2 = 0.;
126
127 float b0_hf = 0.;
128 float b1_hf = 0.;
129
130 float b0_lf = 0.;
131 float b1_lf = 0.;
132
133 float prevX[2] = {};
134 float prevR_lf[2] = {};
135 float prevR_hf[2] = {};
136};
137
138QAmbisonicDecoder::QAmbisonicDecoder(AmbisonicOrder ambisonicOrder, const QAudioFormat &format)
139 : order(ambisonicOrder)
140{
141 Q_ASSERT(order > 0 && order <= 3);
142 inputChannels = (order + 1) * (order + 1);
143 outputChannels = format.channelCount();
144
145 channelConfig = format.channelConfig();
146 if (channelConfig == QAudioFormat::ChannelConfigUnknown)
147 channelConfig = QAudioFormat::defaultChannelConfigForChannelCount(channelCount: format.channelCount());
148
149 if (channelConfig == QAudioFormat::ChannelConfigMono ||
150 channelConfig == QAudioFormat::ChannelConfigStereo ||
151 channelConfig == QAudioFormat::ChannelConfig2Dot1 ||
152 channelConfig == QAudioFormat::ChannelConfig3Dot0 ||
153 channelConfig == QAudioFormat::ChannelConfig3Dot1) {
154 // these are non surround configs and handled manually to avoid
155 // audible holes for sounds coming from behing
156 //
157 // We use a simpler decoding process here, only taking first order
158 // ambisonics into account
159 //
160 // Left and right channels get 50% W and 50% X
161 // Center gets 50% W and 50% Y
162 // LFE gets 50% W
163 simpleDecoderFactors = new float[4*outputChannels];
164 float *r = new float[2*outputChannels]; // reverb output is in stereo
165 float *f = simpleDecoderFactors;
166 reverbFactors = r;
167 if (channelConfig & QAudioFormat::channelConfig(channels: QAudioFormat::FrontLeft)) {
168 f[0] = 0.5f; f[1] = 0.5f; f[2] = 0.; f[3] = 0.f;
169 f += 4;
170 r[0] = 1.; r[1] = 0.;
171 r += 2;
172 }
173 if (channelConfig & QAudioFormat::channelConfig(channels: QAudioFormat::FrontRight)) {
174 f[0] = 0.5f; f[1] = -0.5f; f[2] = 0.; f[3] = 0.f;
175 f += 4;
176 r[0] = 0.; r[1] = 1.;
177 r += 2;
178 }
179 if (channelConfig & QAudioFormat::channelConfig(channels: QAudioFormat::FrontCenter)) {
180 f[0] = 0.5f; f[1] = -0.f; f[2] = 0.; f[3] = 0.5f;
181 f += 4;
182 r[0] = .5; r[1] = .5;
183 r += 2;
184 }
185 if (channelConfig & QAudioFormat::channelConfig(channels: QAudioFormat::LFE)) {
186 f[0] = 0.5f; f[1] = -0.f; f[2] = 0.; f[3] = 0.0f;
187 f += 4;
188 r[0] = 0.; r[1] = 0.;
189 r += 2;
190 }
191 Q_UNUSED(f);
192 Q_UNUSED(r);
193 Q_ASSERT((f - simpleDecoderFactors) == 4*outputChannels);
194 Q_ASSERT((r - reverbFactors) == 2*outputChannels);
195
196 return;
197 }
198
199 for (const auto &d : decoderMap) {
200 if (d.config == channelConfig) {
201 decoderData = &d;
202 reverbFactors = decoderData->reverb;
203 break;
204 }
205 }
206 if (!decoderData) {
207 // can't handle this,
208 outputChannels = 0;
209 return;
210 }
211
212 filters = new QAmbisonicDecoderFilter[inputChannels];
213 for (int i = 0; i < inputChannels; ++i)
214 filters[i].configure(sampleRate: format.sampleRate());
215}
216
217QAmbisonicDecoder::~QAmbisonicDecoder()
218{
219 if (simpleDecoderFactors) {
220 delete[] simpleDecoderFactors;
221 delete[] reverbFactors;
222 }
223}
224
225void QAmbisonicDecoder::processBuffer(const float *input[], float *output, int nSamples)
226{
227 float *o = output;
228 memset(s: o, c: 0, n: nSamples*outputChannels*sizeof(float));
229
230 if (simpleDecoderFactors) {
231 for (int i = 0; i < nSamples; ++i) {
232 for (int j = 0; j < 4; ++j) {
233 for (int k = 0; k < outputChannels; ++k)
234 o[k] += simpleDecoderFactors[k*4 + j]*input[j][i];
235 }
236 o += outputChannels;
237 }
238 return;
239 }
240
241 const float *matrix_hi = decoderData->hf[order - 1];
242 const float *matrix_lo = decoderData->lf[order - 1];
243 for (int i = 0; i < nSamples; ++i) {
244 QAmbisonicDecoderFilter::Output buf[maxAmbisonicChannels];
245 for (int j = 0; j < inputChannels; ++j)
246 buf[j] = filters[j].next(x: input[j][i]);
247 for (int j = 0; j < inputChannels; ++j) {
248 for (int k = 0; k < outputChannels; ++k)
249 o[k] += matrix_lo[k*inputChannels + j]*buf[j].lf + matrix_hi[k*inputChannels + j]*buf[j].hf;
250 }
251 o += outputChannels;
252 }
253}
254
255void QAmbisonicDecoder::processBuffer(const float *input[], short *output, int nSamples)
256{
257 const float *reverb[] = { nullptr, nullptr };
258 return processBufferWithReverb(input, reverb, output, nSamples);
259}
260
261void QAmbisonicDecoder::processBufferWithReverb(const float *input[], const float *reverb[2], short *output, int nSamples)
262{
263 if (simpleDecoderFactors) {
264 for (int i = 0; i < nSamples; ++i) {
265 float o[4] = {};
266 for (int k = 0; k < outputChannels; ++k) {
267 for (int j = 0; j < 4; ++j)
268 o[k] += simpleDecoderFactors[k*4 + j]*input[j][i];
269 }
270 if (reverb[0]) {
271 for (int k = 0; k < outputChannels; ++k) {
272 o[k] += reverb[0][i]*reverbFactors[2*k] + reverb[1][i]*reverbFactors[2*k+1];
273 }
274 }
275
276 for (int k = 0; k < outputChannels; ++k)
277 output[k] = static_cast<short>(o[k]*32768.);
278 output += outputChannels;
279 }
280 return;
281 }
282
283 // qDebug() << "XXX" << inputChannels << outputChannels;
284 const float *matrix_hi = decoderData->hf[order - 1];
285 const float *matrix_lo = decoderData->lf[order - 1];
286 for (int i = 0; i < nSamples; ++i) {
287 QAmbisonicDecoderFilter::Output buf[maxAmbisonicChannels];
288 for (int j = 0; j < inputChannels; ++j)
289 buf[j] = filters[j].next(x: input[j][i]);
290 float o[32] = {}; // we can't support more than 32 channels from our API
291 for (int j = 0; j < inputChannels; ++j) {
292 for (int k = 0; k < outputChannels; ++k)
293 o[k] += matrix_lo[k*inputChannels + j]*buf[j].lf + matrix_hi[k*inputChannels + j]*buf[j].hf;
294 }
295 if (reverb[0]) {
296 for (int k = 0; k < outputChannels; ++k) {
297 o[k] += reverb[0][i]*reverbFactors[2*k] + reverb[1][i]*reverbFactors[2*k+1];
298 }
299 }
300 for (int k = 0; k < outputChannels; ++k)
301 output[k] = static_cast<short>(o[k]*32768.);
302 output += outputChannels;
303 }
304
305}
306
307QT_END_NAMESPACE
308
309

source code of qtmultimedia/src/spatialaudio/qambisonicdecoder.cpp