Revision control
Copy as Markdown
// Copyright (C) 2022 The Qt Company Ltd.
// SPDX-License-Identifier: LicenseRef-Qt-Commercial OR LGPL-3.0-only OR GPL-3.0-only
#include "qambisonicdecoder_p.h"
#include <QtSpatialAudio/private/qambisonicdecoderdata_p.h>
#include <QtMultimedia/private/qaudio_qspan_support_p.h>
#include <QtCore/qdebug.h>
#include <algorithm>
#include <array>
#include <cmath>
QT_BEGIN_NAMESPACE
// We're using a phase matched band splitting filter to split the ambisonic signal into a low
// and high frequency component and apply matrix conversions to those components individually
// as described in the document.
//
// We are currently not using a near field compensation filter, something that could potentially
// improve sound quality further.
//
// For mono and stereo decoding, we use a simpler algorithm to avoid artificially dampening signals
// coming from the back, as we do not have any speakers in that direction and the calculations
// through matlab would give us audible 'holes'.
struct QAmbisonicDecoderData
{
const QAudioFormat::ChannelConfig config;
const std::array<const float *, 3> lf;
const std::array<const float *, 3> hf;
const float *const reverb;
};
constexpr float reverb_x_0[] = {
1.f, 0.f, // L
0.f, 1.f, // R
.7f, .7f, // C
1.f, 0.f, // Ls
0.f, 1.f, // Rs
1.f, 0.f, // Lb
0.f, 1.f, // Rb
};
constexpr float reverb_x_1[] = {
1.f, 0.f, // L
0.f, 1.f, // R
.7f, .7f, // C
.0f, .0f, // LFE
1.f, 0.f, // Ls
0.f, 1.f, // Rs
1.f, 0.f, // Lb
0.f, 1.f, // Rb
};
static constexpr QAmbisonicDecoderData decoderMap[] = {
{ QAudioFormat::ChannelConfigSurround5Dot0,
{{ decoderMatrix_5dot0_1_lf, decoderMatrix_5dot0_2_lf, decoderMatrix_5dot0_3_lf }},
{{ decoderMatrix_5dot0_1_hf, decoderMatrix_5dot0_2_hf, decoderMatrix_5dot0_3_hf }},
reverb_x_0 },
{ QAudioFormat::ChannelConfigSurround5Dot1,
{{ decoderMatrix_5dot1_1_lf, decoderMatrix_5dot1_2_lf, decoderMatrix_5dot1_3_lf }},
{{ decoderMatrix_5dot1_1_hf, decoderMatrix_5dot1_2_hf, decoderMatrix_5dot1_3_hf }},
reverb_x_1 },
{ QAudioFormat::ChannelConfigSurround7Dot0,
{{ decoderMatrix_7dot0_1_lf, decoderMatrix_7dot0_2_lf, decoderMatrix_7dot0_3_lf }},
{{ decoderMatrix_7dot0_1_hf, decoderMatrix_7dot0_2_hf, decoderMatrix_7dot0_3_hf }},
reverb_x_0 },
{ QAudioFormat::ChannelConfigSurround7Dot1,
{{ decoderMatrix_7dot1_1_lf, decoderMatrix_7dot1_2_lf, decoderMatrix_7dot1_3_lf }},
{{ decoderMatrix_7dot1_1_hf, decoderMatrix_7dot1_2_hf, decoderMatrix_7dot1_3_hf }},
reverb_x_1 }
};
// Implements a split second order IIR filter (TDF-II)
// The audio data is split into a phase synced low and high frequency part
// This allows us to apply different factors to both parts for better sound
// localization when converting from ambisonic formats
//
class QAmbisonicDecoderFilter
{
public:
QAmbisonicDecoderFilter() = default;
void configure(float sampleRate, float cutoffFrequency = 380)
{
double k = std::tan(M_PI * cutoffFrequency / sampleRate);
double denom = k * k + 2.0 * k + 1.0;
a1 = 2.0 * (k * k - 1.0) / denom;
a2 = (k * k - 2.0 * k + 1.0) / denom;
b0_lf = (k * k) / denom;
b1_lf = 2.0 * b0_lf;
b0_hf = 1.0 / denom;
b1_hf = -2.0 * b0_hf;
}
struct Output
{
float lf;
float hf;
};
Output next(float x)
{
#ifdef Q_PROCESSOR_X86
// tiny DC offset to prevent denormals, which can cause severe performance degradation on x86 CPUs
x += 1e-18f;
#endif
// Process LF
double r_lf = x * b0_lf + s1_lf;
s1_lf = x * b1_lf - r_lf * a1 + s2_lf;
s2_lf = x * b0_lf - r_lf * a2;
// Process HF
double r_hf = x * b0_hf + s1_hf;
s1_hf = x * b1_hf - r_hf * a1 + s2_hf;
s2_hf = x * b0_hf - r_hf * a2;
return Output{
float(r_lf),
float(r_hf),
};
}
private:
// coefficients
double a1 = 0.;
double a2 = 0.;
double b0_hf = 0.;
double b1_hf = 0.;
double b0_lf = 0.;
double b1_lf = 0.;
// state
double s1_lf = 0.0, s2_lf = 0.0;
double s1_hf = 0.0, s2_hf = 0.0;
};
namespace {
int inputChannelsForAmbisonicOrder(QAmbisonicDecoder::AmbisonicOrder ambisonicOrder)
{
auto order = qToUnderlying(ambisonicOrder);
Q_ASSERT(order > 0 && order <= 3);
return (order + 1) * (order + 1);
}
QAudioFormat::ChannelConfig ambisonicDecoderChannelConfig(QAudioFormat::ChannelConfig channelConfig,
int numberOfOutputChannels)
{
if (channelConfig == QAudioFormat::ChannelConfigUnknown)
channelConfig = QAudioFormat::defaultChannelConfigForChannelCount(numberOfOutputChannels);
return channelConfig;
}
} // namespace
QAmbisonicDecoder::QAmbisonicDecoder(AmbisonicOrder ambisonicOrder, int sampleRate,
int numberOfOutputChannels,
QAudioFormat::ChannelConfig channelCfg)
: channelConfig(ambisonicDecoderChannelConfig(channelCfg, numberOfOutputChannels)),
order(ambisonicOrder),
inputChannels(inputChannelsForAmbisonicOrder(ambisonicOrder)),
outputChannels(numberOfOutputChannels)
{
if (channelConfig == QAudioFormat::ChannelConfigMono ||
channelConfig == QAudioFormat::ChannelConfigStereo ||
channelConfig == QAudioFormat::ChannelConfig2Dot1 ||
channelConfig == QAudioFormat::ChannelConfig3Dot0 ||
channelConfig == QAudioFormat::ChannelConfig3Dot1) {
// these are non surround configs and handled manually to avoid
// audible holes for sounds coming from behing
//
// We use a simpler decoding process here, only taking first order
// ambisonics into account
//
// Left and right channels get 50% W and 50% X
// Center gets 50% W and 50% Y
// LFE gets 50% W
simpleDecoderFactors = std::make_unique<float[]>(4 * outputChannels);
m_reverbFactorsOwned = std::make_unique<float[]>(2 * outputChannels); // reverb output is in stereo
reverbFactors = m_reverbFactorsOwned.get();
float *f = simpleDecoderFactors.get();
float *r = m_reverbFactorsOwned.get();
if (channelConfig & QAudioFormat::channelConfig(QAudioFormat::FrontLeft)) {
f[0] = 0.5f; f[1] = 0.5f; f[2] = 0.; f[3] = 0.f;
f += 4;
r[0] = 1.; r[1] = 0.;
r += 2;
}
if (channelConfig & QAudioFormat::channelConfig(QAudioFormat::FrontRight)) {
f[0] = 0.5f; f[1] = -0.5f; f[2] = 0.; f[3] = 0.f;
f += 4;
r[0] = 0.; r[1] = 1.;
r += 2;
}
if (channelConfig & QAudioFormat::channelConfig(QAudioFormat::FrontCenter)) {
f[0] = 0.5f; f[1] = -0.f; f[2] = 0.; f[3] = 0.5f;
f += 4;
r[0] = .5; r[1] = .5;
r += 2;
}
if (channelConfig & QAudioFormat::channelConfig(QAudioFormat::LFE)) {
f[0] = 0.5f; f[1] = -0.f; f[2] = 0.; f[3] = 0.0f;
f += 4;
r[0] = 0.; r[1] = 0.;
r += 2;
}
Q_UNUSED(f);
Q_UNUSED(r);
Q_ASSERT((f - simpleDecoderFactors.get()) == 4 * outputChannels);
Q_ASSERT((r - reverbFactors) == 2*outputChannels);
return;
}
for (const auto &d : decoderMap) {
if (d.config == channelConfig) {
decoderData = &d;
reverbFactors = decoderData->reverb;
break;
}
}
if (!decoderData)
return;
filters = std::make_unique<QAmbisonicDecoderFilter[]>(inputChannels);
for (int i = 0; i < inputChannels; ++i)
filters[i].configure(sampleRate);
}
QAmbisonicDecoder::~QAmbisonicDecoder() = default;
void QAmbisonicDecoder::processBuffer(QSpan<const float *> input, QSpan<float> output)
{
const int nSamples = int(output.size()) / outputChannels;
std::fill(output.begin(), output.end(), 0.f);
float *o = output.data();
if (simpleDecoderFactors) {
for (int i = 0; i < nSamples; ++i) {
for (int j = 0; j < 4; ++j) {
for (int k = 0; k < outputChannels; ++k)
o[k] += simpleDecoderFactors[k*4 + j]*input[j][i];
}
o += outputChannels;
}
return;
}
const float *matrix_hi = decoderData->hf[order - 1];
const float *matrix_lo = decoderData->lf[order - 1];
for (int i = 0; i < nSamples; ++i) {
std::array<QAmbisonicDecoderFilter::Output, maxAmbisonicChannels> buf;
for (int j = 0; j < inputChannels; ++j)
buf[j] = filters[j].next(input[j][i]);
for (int j = 0; j < inputChannels; ++j) {
for (int k = 0; k < outputChannels; ++k)
o[k] += matrix_lo[k*inputChannels + j]*buf[j].lf + matrix_hi[k*inputChannels + j]*buf[j].hf;
}
o += outputChannels;
}
}
void QAmbisonicDecoder::processBufferWithReverb(QSpan<const float *> input,
QSpan<const float *, 2> reverb, QSpan<float> output)
{
Q_ASSERT(outputChannels > 0);
Q_ASSERT(int(output.size()) % outputChannels == 0);
const int nSamples = int(output.size()) / outputChannels;
std::fill(output.begin(), output.end(), 0.f);
float *o = output.data();
if (simpleDecoderFactors) {
for (int i = 0; i < nSamples; ++i) {
for (int k = 0; k < outputChannels; ++k) {
for (int j = 0; j < 4; ++j)
o[k] += simpleDecoderFactors[k*4 + j]*input[j][i];
}
if (reverb[0]) {
for (int k = 0; k < outputChannels; ++k)
o[k] += reverb[0][i] * reverbFactors[2 * k]
+ reverb[1][i] * reverbFactors[2 * k + 1];
}
o += outputChannels;
}
return;
}
Q_ASSERT(filters);
const float *matrix_hi = decoderData->hf[order - 1];
const float *matrix_lo = decoderData->lf[order - 1];
for (int i = 0; i < nSamples; ++i) {
std::array<QAmbisonicDecoderFilter::Output, maxAmbisonicChannels> buf;
for (int j = 0; j < inputChannels; ++j)
buf[j] = filters[j].next(input[j][i]);
for (int j = 0; j < inputChannels; ++j) {
for (int k = 0; k < outputChannels; ++k)
o[k] += matrix_lo[k*inputChannels + j]*buf[j].lf + matrix_hi[k*inputChannels + j]*buf[j].hf;
}
if (reverb[0]) {
for (int k = 0; k < outputChannels; ++k)
o[k] += reverb[0][i]*reverbFactors[2*k] + reverb[1][i]*reverbFactors[2*k+1];
}
o += outputChannels;
}
}
QT_END_NAMESPACE