update openal-soft

sync point: master-ac5d40e40a0155351fe1be4aab30017b6a13a859
This commit is contained in:
AzaezelX 2021-01-26 13:01:35 -06:00
parent 762a84550f
commit 3603188b7f
365 changed files with 76053 additions and 53126 deletions

View file

@ -0,0 +1,375 @@
#include "config.h"
#include "ambdec.h"
#include <algorithm>
#include <cctype>
#include <cstddef>
#include <iterator>
#include <sstream>
#include <string>
#include "alfstream.h"
#include "core/logging.h"
namespace {
template<typename T, std::size_t N>
constexpr inline std::size_t size(const T(&)[N]) noexcept
{ return N; }
int readline(std::istream &f, std::string &output)
{
while(f.good() && f.peek() == '\n')
f.ignore();
return std::getline(f, output) && !output.empty();
}
bool read_clipped_line(std::istream &f, std::string &buffer)
{
while(readline(f, buffer))
{
std::size_t pos{0};
while(pos < buffer.length() && std::isspace(buffer[pos]))
pos++;
buffer.erase(0, pos);
std::size_t cmtpos{buffer.find_first_of('#')};
if(cmtpos < buffer.length())
buffer.resize(cmtpos);
while(!buffer.empty() && std::isspace(buffer.back()))
buffer.pop_back();
if(!buffer.empty())
return true;
}
return false;
}
std::string read_word(std::istream &f)
{
std::string ret;
f >> ret;
return ret;
}
bool is_at_end(const std::string &buffer, std::size_t endpos)
{
while(endpos < buffer.length() && std::isspace(buffer[endpos]))
++endpos;
return !(endpos < buffer.length());
}
al::optional<std::string> load_ambdec_speakers(AmbDecConf::SpeakerConf *spkrs,
const std::size_t num_speakers, std::istream &f, std::string &buffer)
{
size_t cur_speaker{0};
while(cur_speaker < num_speakers)
{
std::istringstream istr{buffer};
std::string cmd{read_word(istr)};
if(cmd.empty())
{
if(!read_clipped_line(f, buffer))
return al::make_optional<std::string>("Unexpected end of file");
continue;
}
if(cmd == "add_spkr")
{
AmbDecConf::SpeakerConf &spkr = spkrs[cur_speaker++];
const size_t spkr_num{cur_speaker};
istr >> spkr.Name;
if(istr.fail()) WARN("Name not specified for speaker %zu\n", spkr_num);
istr >> spkr.Distance;
if(istr.fail()) WARN("Distance not specified for speaker %zu\n", spkr_num);
istr >> spkr.Azimuth;
if(istr.fail()) WARN("Azimuth not specified for speaker %zu\n", spkr_num);
istr >> spkr.Elevation;
if(istr.fail()) WARN("Elevation not specified for speaker %zu\n", spkr_num);
istr >> spkr.Connection;
if(istr.fail()) TRACE("Connection not specified for speaker %zu\n", spkr_num);
}
else
return al::make_optional("Unexpected speakers command: "+cmd);
istr.clear();
const auto endpos = static_cast<std::size_t>(istr.tellg());
if(!is_at_end(buffer, endpos))
return al::make_optional("Extra junk on line: " + buffer.substr(endpos));
buffer.clear();
}
return al::nullopt;
}
al::optional<std::string> load_ambdec_matrix(float (&gains)[MaxAmbiOrder+1],
AmbDecConf::CoeffArray *matrix, const std::size_t maxrow, std::istream &f, std::string &buffer)
{
bool gotgains{false};
std::size_t cur{0u};
while(cur < maxrow)
{
std::istringstream istr{buffer};
std::string cmd{read_word(istr)};
if(cmd.empty())
{
if(!read_clipped_line(f, buffer))
return al::make_optional<std::string>("Unexpected end of file");
continue;
}
if(cmd == "order_gain")
{
std::size_t curgain{0u};
float value;
while(istr.good())
{
istr >> value;
if(istr.fail()) break;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk on gain "+std::to_string(curgain+1)+": "+
buffer.substr(static_cast<std::size_t>(istr.tellg())));
if(curgain < size(gains))
gains[curgain++] = value;
}
std::fill(std::begin(gains)+curgain, std::end(gains), 0.0f);
gotgains = true;
}
else if(cmd == "add_row")
{
AmbDecConf::CoeffArray &mtxrow = matrix[cur++];
std::size_t curidx{0u};
float value{};
while(istr.good())
{
istr >> value;
if(istr.fail()) break;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk on matrix element "+
std::to_string(curidx)+"x"+std::to_string(cur-1)+": "+
buffer.substr(static_cast<std::size_t>(istr.tellg())));
if(curidx < mtxrow.size())
mtxrow[curidx++] = value;
}
std::fill(mtxrow.begin()+curidx, mtxrow.end(), 0.0f);
}
else
return al::make_optional("Unexpected matrix command: "+cmd);
istr.clear();
const auto endpos = static_cast<std::size_t>(istr.tellg());
if(!is_at_end(buffer, endpos))
return al::make_optional("Extra junk on line: " + buffer.substr(endpos));
buffer.clear();
}
if(!gotgains)
return al::make_optional<std::string>("Matrix order_gain not specified");
return al::nullopt;
}
} // namespace
al::optional<std::string> AmbDecConf::load(const char *fname) noexcept
{
al::ifstream f{fname};
if(!f.is_open())
return al::make_optional<std::string>("Failed to open file");
bool speakers_loaded{false};
bool matrix_loaded{false};
bool lfmatrix_loaded{false};
std::string buffer;
while(read_clipped_line(f, buffer))
{
std::istringstream istr{buffer};
std::string command{read_word(istr)};
if(command.empty())
return al::make_optional("Malformed line: "+buffer);
if(command == "/description")
istr >> Description;
else if(command == "/version")
{
istr >> Version;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk after version: " +
buffer.substr(static_cast<std::size_t>(istr.tellg())));
if(Version != 3)
return al::make_optional("Unsupported version: "+std::to_string(Version));
}
else if(command == "/dec/chan_mask")
{
if(ChanMask)
return al::make_optional<std::string>("Duplicate chan_mask definition");
istr >> std::hex >> ChanMask >> std::dec;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk after mask: " +
buffer.substr(static_cast<std::size_t>(istr.tellg())));
if(!ChanMask)
return al::make_optional("Invalid chan_mask: "+std::to_string(ChanMask));
}
else if(command == "/dec/freq_bands")
{
if(FreqBands)
return al::make_optional<std::string>("Duplicate freq_bands");
istr >> FreqBands;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk after freq_bands: " +
buffer.substr(static_cast<std::size_t>(istr.tellg())));
if(FreqBands != 1 && FreqBands != 2)
return al::make_optional("Invalid freq_bands: "+std::to_string(FreqBands));
}
else if(command == "/dec/speakers")
{
if(NumSpeakers)
return al::make_optional<std::string>("Duplicate speakers");
istr >> NumSpeakers;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk after speakers: " +
buffer.substr(static_cast<std::size_t>(istr.tellg())));
if(!NumSpeakers)
return al::make_optional("Invalid speakers: "+std::to_string(NumSpeakers));
Speakers = std::make_unique<SpeakerConf[]>(NumSpeakers);
}
else if(command == "/dec/coeff_scale")
{
std::string scale = read_word(istr);
if(scale == "n3d") CoeffScale = AmbDecScale::N3D;
else if(scale == "sn3d") CoeffScale = AmbDecScale::SN3D;
else if(scale == "fuma") CoeffScale = AmbDecScale::FuMa;
else
return al::make_optional("Unexpected coeff_scale: "+scale);
}
else if(command == "/opt/xover_freq")
{
istr >> XOverFreq;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk after xover_freq: " +
buffer.substr(static_cast<std::size_t>(istr.tellg())));
}
else if(command == "/opt/xover_ratio")
{
istr >> XOverRatio;
if(!istr.eof() && !std::isspace(istr.peek()))
return al::make_optional("Extra junk after xover_ratio: " +
buffer.substr(static_cast<std::size_t>(istr.tellg())));
}
else if(command == "/opt/input_scale" || command == "/opt/nfeff_comp" ||
command == "/opt/delay_comp" || command == "/opt/level_comp")
{
/* Unused */
read_word(istr);
}
else if(command == "/speakers/{")
{
if(!NumSpeakers)
return al::make_optional<std::string>("Speakers defined without a count");
const auto endpos = static_cast<std::size_t>(istr.tellg());
if(!is_at_end(buffer, endpos))
return al::make_optional("Extra junk on line: " + buffer.substr(endpos));
buffer.clear();
if(auto err = load_ambdec_speakers(Speakers.get(), NumSpeakers, f, buffer))
return err;
speakers_loaded = true;
if(!read_clipped_line(f, buffer))
return al::make_optional<std::string>("Unexpected end of file");
std::istringstream istr2{buffer};
std::string endmark{read_word(istr2)};
if(endmark != "/}")
return al::make_optional("Expected /} after speaker definitions, got "+endmark);
istr.swap(istr2);
}
else if(command == "/lfmatrix/{" || command == "/hfmatrix/{" || command == "/matrix/{")
{
if(!NumSpeakers)
return al::make_optional<std::string>("Matrix defined without a count");
const auto endpos = static_cast<std::size_t>(istr.tellg());
if(!is_at_end(buffer, endpos))
return al::make_optional("Extra junk on line: " + buffer.substr(endpos));
buffer.clear();
if(!Matrix)
{
Matrix = std::make_unique<CoeffArray[]>(NumSpeakers * FreqBands);
LFMatrix = Matrix.get();
HFMatrix = LFMatrix + NumSpeakers*(FreqBands-1);
}
if(FreqBands == 1)
{
if(command != "/matrix/{")
return al::make_optional(
"Unexpected \""+command+"\" type for a single-band decoder");
if(auto err = load_ambdec_matrix(HFOrderGain, HFMatrix, NumSpeakers, f, buffer))
return err;
matrix_loaded = true;
}
else
{
if(command == "/lfmatrix/{")
{
if(auto err=load_ambdec_matrix(LFOrderGain, LFMatrix, NumSpeakers, f, buffer))
return err;
lfmatrix_loaded = true;
}
else if(command == "/hfmatrix/{")
{
if(auto err=load_ambdec_matrix(HFOrderGain, HFMatrix, NumSpeakers, f, buffer))
return err;
matrix_loaded = true;
}
else
return al::make_optional(
"Unexpected \""+command+"\" type for a dual-band decoder");
}
if(!read_clipped_line(f, buffer))
return al::make_optional<std::string>("Unexpected end of file");
std::istringstream istr2{buffer};
std::string endmark{read_word(istr2)};
if(endmark != "/}")
return al::make_optional("Expected /} after matrix definitions, got "+endmark);
istr.swap(istr2);
}
else if(command == "/end")
{
const auto endpos = static_cast<std::size_t>(istr.tellg());
if(!is_at_end(buffer, endpos))
return al::make_optional("Extra junk on end: " + buffer.substr(endpos));
if(!speakers_loaded || !matrix_loaded || (FreqBands == 2 && !lfmatrix_loaded))
return al::make_optional<std::string>("No decoder defined");
return al::nullopt;
}
else
return al::make_optional("Unexpected command: " + command);
istr.clear();
const auto endpos = static_cast<std::size_t>(istr.tellg());
if(!is_at_end(buffer, endpos))
return al::make_optional("Extra junk on line: " + buffer.substr(endpos));
buffer.clear();
}
return al::make_optional<std::string>("Unexpected end of file");
}

View file

@ -0,0 +1,52 @@
#ifndef CORE_AMBDEC_H
#define CORE_AMBDEC_H
#include <array>
#include <memory>
#include <string>
#include "aloptional.h"
#include "core/ambidefs.h"
/* Helpers to read .ambdec configuration files. */
enum class AmbDecScale {
N3D,
SN3D,
FuMa,
};
struct AmbDecConf {
std::string Description;
int Version{0}; /* Must be 3 */
unsigned int ChanMask{0u};
unsigned int FreqBands{0u}; /* Must be 1 or 2 */
AmbDecScale CoeffScale{};
float XOverFreq{0.0f};
float XOverRatio{0.0f};
struct SpeakerConf {
std::string Name;
float Distance{0.0f};
float Azimuth{0.0f};
float Elevation{0.0f};
std::string Connection;
};
size_t NumSpeakers{0};
std::unique_ptr<SpeakerConf[]> Speakers;
using CoeffArray = std::array<float,MaxAmbiChannels>;
std::unique_ptr<CoeffArray[]> Matrix;
/* Unused when FreqBands == 1 */
float LFOrderGain[MaxAmbiOrder+1]{};
CoeffArray *LFMatrix;
float HFOrderGain[MaxAmbiOrder+1]{};
CoeffArray *HFMatrix;
al::optional<std::string> load(const char *fname) noexcept;
};
#endif /* CORE_AMBDEC_H */

View file

@ -0,0 +1,171 @@
#ifndef CORE_AMBIDEFS_H
#define CORE_AMBIDEFS_H
#include <array>
#include <stddef.h>
#include <stdint.h>
using uint = unsigned int;
/* The maximum number of Ambisonics channels. For a given order (o), the size
* needed will be (o+1)**2, thus zero-order has 1, first-order has 4, second-
* order has 9, third-order has 16, and fourth-order has 25.
*/
constexpr uint8_t MaxAmbiOrder{3};
constexpr inline size_t AmbiChannelsFromOrder(size_t order) noexcept
{ return (order+1) * (order+1); }
constexpr size_t MaxAmbiChannels{AmbiChannelsFromOrder(MaxAmbiOrder)};
/* A bitmask of ambisonic channels for 0 to 4th order. This only specifies up
* to 4th order, which is the highest order a 32-bit mask value can specify (a
* 64-bit mask could handle up to 7th order).
*/
constexpr uint Ambi0OrderMask{0x00000001};
constexpr uint Ambi1OrderMask{0x0000000f};
constexpr uint Ambi2OrderMask{0x000001ff};
constexpr uint Ambi3OrderMask{0x0000ffff};
constexpr uint Ambi4OrderMask{0x01ffffff};
/* A bitmask of ambisonic channels with height information. If none of these
* channels are used/needed, there's no height (e.g. with most surround sound
* speaker setups). This is ACN ordering, with bit 0 being ACN 0, etc.
*/
constexpr uint AmbiPeriphonicMask{0xfe7ce4};
/* The maximum number of ambisonic channels for 2D (non-periphonic)
* representation. This is 2 per each order above zero-order, plus 1 for zero-
* order. Or simply, o*2 + 1.
*/
constexpr inline size_t Ambi2DChannelsFromOrder(size_t order) noexcept
{ return order*2 + 1; }
constexpr size_t MaxAmbi2DChannels{Ambi2DChannelsFromOrder(MaxAmbiOrder)};
/* NOTE: These are scale factors as applied to Ambisonics content. Decoder
* coefficients should be divided by these values to get proper scalings.
*/
struct AmbiScale {
static auto& FromN3D() noexcept
{
static constexpr const std::array<float,MaxAmbiChannels> ret{{
1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f,
1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f
}};
return ret;
}
static auto& FromSN3D() noexcept
{
static constexpr const std::array<float,MaxAmbiChannels> ret{{
1.000000000f, /* ACN 0, sqrt(1) */
1.732050808f, /* ACN 1, sqrt(3) */
1.732050808f, /* ACN 2, sqrt(3) */
1.732050808f, /* ACN 3, sqrt(3) */
2.236067978f, /* ACN 4, sqrt(5) */
2.236067978f, /* ACN 5, sqrt(5) */
2.236067978f, /* ACN 6, sqrt(5) */
2.236067978f, /* ACN 7, sqrt(5) */
2.236067978f, /* ACN 8, sqrt(5) */
2.645751311f, /* ACN 9, sqrt(7) */
2.645751311f, /* ACN 10, sqrt(7) */
2.645751311f, /* ACN 11, sqrt(7) */
2.645751311f, /* ACN 12, sqrt(7) */
2.645751311f, /* ACN 13, sqrt(7) */
2.645751311f, /* ACN 14, sqrt(7) */
2.645751311f, /* ACN 15, sqrt(7) */
}};
return ret;
}
static auto& FromFuMa() noexcept
{
static constexpr const std::array<float,MaxAmbiChannels> ret{{
1.414213562f, /* ACN 0 (W), sqrt(2) */
1.732050808f, /* ACN 1 (Y), sqrt(3) */
1.732050808f, /* ACN 2 (Z), sqrt(3) */
1.732050808f, /* ACN 3 (X), sqrt(3) */
1.936491673f, /* ACN 4 (V), sqrt(15)/2 */
1.936491673f, /* ACN 5 (T), sqrt(15)/2 */
2.236067978f, /* ACN 6 (R), sqrt(5) */
1.936491673f, /* ACN 7 (S), sqrt(15)/2 */
1.936491673f, /* ACN 8 (U), sqrt(15)/2 */
2.091650066f, /* ACN 9 (Q), sqrt(35/8) */
1.972026594f, /* ACN 10 (O), sqrt(35)/3 */
2.231093404f, /* ACN 11 (M), sqrt(224/45) */
2.645751311f, /* ACN 12 (K), sqrt(7) */
2.231093404f, /* ACN 13 (L), sqrt(224/45) */
1.972026594f, /* ACN 14 (N), sqrt(35)/3 */
2.091650066f, /* ACN 15 (P), sqrt(35/8) */
}};
return ret;
}
};
struct AmbiIndex {
static auto& FromFuMa() noexcept
{
static constexpr const std::array<uint8_t,MaxAmbiChannels> ret{{
0, /* W */
3, /* X */
1, /* Y */
2, /* Z */
6, /* R */
7, /* S */
5, /* T */
8, /* U */
4, /* V */
12, /* K */
13, /* L */
11, /* M */
14, /* N */
10, /* O */
15, /* P */
9, /* Q */
}};
return ret;
}
static auto& FromFuMa2D() noexcept
{
static constexpr const std::array<uint8_t,MaxAmbi2DChannels> ret{{
0, /* W */
3, /* X */
1, /* Y */
8, /* U */
4, /* V */
15, /* P */
9, /* Q */
}};
return ret;
}
static auto& FromACN() noexcept
{
static constexpr const std::array<uint8_t,MaxAmbiChannels> ret{{
0, 1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15
}};
return ret;
}
static auto& FromACN2D() noexcept
{
static constexpr const std::array<uint8_t,MaxAmbi2DChannels> ret{{
0, 1,3, 4,8, 9,15
}};
return ret;
}
static auto& OrderFromChannel() noexcept
{
static constexpr const std::array<uint8_t,MaxAmbiChannels> ret{{
0, 1,1,1, 2,2,2,2,2, 3,3,3,3,3,3,3,
}};
return ret;
}
static auto& OrderFrom2DChannel() noexcept
{
static constexpr const std::array<uint8_t,MaxAmbi2DChannels> ret{{
0, 1,1, 2,2, 3,3,
}};
return ret;
}
};
#endif /* CORE_AMBIDEFS_H */

View file

@ -0,0 +1,183 @@
/*-
* Copyright (c) 2005 Boris Mikhaylov
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "config.h"
#include <algorithm>
#include <cmath>
#include <iterator>
#include "bs2b.h"
#include "math_defs.h"
/* Set up all data. */
static void init(struct bs2b *bs2b)
{
float Fc_lo, Fc_hi;
float G_lo, G_hi;
float x, g;
switch(bs2b->level)
{
case BS2B_LOW_CLEVEL: /* Low crossfeed level */
Fc_lo = 360.0f;
Fc_hi = 501.0f;
G_lo = 0.398107170553497f;
G_hi = 0.205671765275719f;
break;
case BS2B_MIDDLE_CLEVEL: /* Middle crossfeed level */
Fc_lo = 500.0f;
Fc_hi = 711.0f;
G_lo = 0.459726988530872f;
G_hi = 0.228208484414988f;
break;
case BS2B_HIGH_CLEVEL: /* High crossfeed level (virtual speakers are closer to itself) */
Fc_lo = 700.0f;
Fc_hi = 1021.0f;
G_lo = 0.530884444230988f;
G_hi = 0.250105790667544f;
break;
case BS2B_LOW_ECLEVEL: /* Low easy crossfeed level */
Fc_lo = 360.0f;
Fc_hi = 494.0f;
G_lo = 0.316227766016838f;
G_hi = 0.168236228897329f;
break;
case BS2B_MIDDLE_ECLEVEL: /* Middle easy crossfeed level */
Fc_lo = 500.0f;
Fc_hi = 689.0f;
G_lo = 0.354813389233575f;
G_hi = 0.187169483835901f;
break;
default: /* High easy crossfeed level */
bs2b->level = BS2B_HIGH_ECLEVEL;
Fc_lo = 700.0f;
Fc_hi = 975.0f;
G_lo = 0.398107170553497f;
G_hi = 0.205671765275719f;
break;
} /* switch */
g = 1.0f / (1.0f - G_hi + G_lo);
/* $fc = $Fc / $s;
* $d = 1 / 2 / pi / $fc;
* $x = exp(-1 / $d);
*/
x = std::exp(-al::MathDefs<float>::Tau() * Fc_lo / static_cast<float>(bs2b->srate));
bs2b->b1_lo = x;
bs2b->a0_lo = G_lo * (1.0f - x) * g;
x = std::exp(-al::MathDefs<float>::Tau() * Fc_hi / static_cast<float>(bs2b->srate));
bs2b->b1_hi = x;
bs2b->a0_hi = (1.0f - G_hi * (1.0f - x)) * g;
bs2b->a1_hi = -x * g;
} /* init */
/* Exported functions.
* See descriptions in "bs2b.h"
*/
void bs2b_set_params(struct bs2b *bs2b, int level, int srate)
{
if(srate <= 0) srate = 1;
bs2b->level = level;
bs2b->srate = srate;
init(bs2b);
} /* bs2b_set_params */
int bs2b_get_level(struct bs2b *bs2b)
{
return bs2b->level;
} /* bs2b_get_level */
int bs2b_get_srate(struct bs2b *bs2b)
{
return bs2b->srate;
} /* bs2b_get_srate */
void bs2b_clear(struct bs2b *bs2b)
{
std::fill(std::begin(bs2b->history), std::end(bs2b->history), bs2b::t_last_sample{});
} /* bs2b_clear */
void bs2b_cross_feed(struct bs2b *bs2b, float *Left, float *Right, size_t SamplesToDo)
{
const float a0_lo{bs2b->a0_lo};
const float b1_lo{bs2b->b1_lo};
const float a0_hi{bs2b->a0_hi};
const float a1_hi{bs2b->a1_hi};
const float b1_hi{bs2b->b1_hi};
float lsamples[128][2];
float rsamples[128][2];
for(size_t base{0};base < SamplesToDo;)
{
const size_t todo{std::min<size_t>(128, SamplesToDo-base)};
/* Process left input */
float z_lo{bs2b->history[0].lo};
float z_hi{bs2b->history[0].hi};
for(size_t i{0};i < todo;i++)
{
lsamples[i][0] = a0_lo*Left[i] + z_lo;
z_lo = b1_lo*lsamples[i][0];
lsamples[i][1] = a0_hi*Left[i] + z_hi;
z_hi = a1_hi*Left[i] + b1_hi*lsamples[i][1];
}
bs2b->history[0].lo = z_lo;
bs2b->history[0].hi = z_hi;
/* Process right input */
z_lo = bs2b->history[1].lo;
z_hi = bs2b->history[1].hi;
for(size_t i{0};i < todo;i++)
{
rsamples[i][0] = a0_lo*Right[i] + z_lo;
z_lo = b1_lo*rsamples[i][0];
rsamples[i][1] = a0_hi*Right[i] + z_hi;
z_hi = a1_hi*Right[i] + b1_hi*rsamples[i][1];
}
bs2b->history[1].lo = z_lo;
bs2b->history[1].hi = z_hi;
/* Crossfeed */
for(size_t i{0};i < todo;i++)
*(Left++) = lsamples[i][1] + rsamples[i][0];
for(size_t i{0};i < todo;i++)
*(Right++) = rsamples[i][1] + lsamples[i][0];
base += todo;
}
} /* bs2b_cross_feed */

View file

@ -0,0 +1,89 @@
/*-
* Copyright (c) 2005 Boris Mikhaylov
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef CORE_BS2B_H
#define CORE_BS2B_H
#include "almalloc.h"
/* Number of crossfeed levels */
#define BS2B_CLEVELS 3
/* Normal crossfeed levels */
#define BS2B_HIGH_CLEVEL 3
#define BS2B_MIDDLE_CLEVEL 2
#define BS2B_LOW_CLEVEL 1
/* Easy crossfeed levels */
#define BS2B_HIGH_ECLEVEL BS2B_HIGH_CLEVEL + BS2B_CLEVELS
#define BS2B_MIDDLE_ECLEVEL BS2B_MIDDLE_CLEVEL + BS2B_CLEVELS
#define BS2B_LOW_ECLEVEL BS2B_LOW_CLEVEL + BS2B_CLEVELS
/* Default crossfeed levels */
#define BS2B_DEFAULT_CLEVEL BS2B_HIGH_ECLEVEL
/* Default sample rate (Hz) */
#define BS2B_DEFAULT_SRATE 44100
struct bs2b {
int level; /* Crossfeed level */
int srate; /* Sample rate (Hz) */
/* Lowpass IIR filter coefficients */
float a0_lo;
float b1_lo;
/* Highboost IIR filter coefficients */
float a0_hi;
float a1_hi;
float b1_hi;
/* Buffer of filter history
* [0] - first channel, [1] - second channel
*/
struct t_last_sample {
float lo;
float hi;
} history[2];
DEF_NEWDEL(bs2b)
};
/* Clear buffers and set new coefficients with new crossfeed level and sample
* rate values.
* level - crossfeed level of *LEVEL values.
* srate - sample rate by Hz.
*/
void bs2b_set_params(bs2b *bs2b, int level, int srate);
/* Return current crossfeed level value */
int bs2b_get_level(bs2b *bs2b);
/* Return current sample rate value */
int bs2b_get_srate(bs2b *bs2b);
/* Clear buffer */
void bs2b_clear(bs2b *bs2b);
void bs2b_cross_feed(bs2b *bs2b, float *Left, float *Right, size_t SamplesToDo);
#endif /* CORE_BS2B_H */

View file

@ -0,0 +1,16 @@
#ifndef CORE_BSINC_DEFS_H
#define CORE_BSINC_DEFS_H
/* The number of distinct scale and phase intervals within the filter table. */
constexpr unsigned int BSincScaleBits{4};
constexpr unsigned int BSincScaleCount{1 << BSincScaleBits};
constexpr unsigned int BSincPhaseBits{5};
constexpr unsigned int BSincPhaseCount{1 << BSincPhaseBits};
/* The maximum number of sample points for the bsinc filters. The max points
* includes the doubling for downsampling, so the maximum number of base sample
* points is 24, which is 23rd order.
*/
constexpr unsigned int BSincPointsMax{48};
#endif /* CORE_BSINC_DEFS_H */

View file

@ -0,0 +1,288 @@
#include "bsinc_tables.h"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <limits>
#include <memory>
#include <stdexcept>
#include "math_defs.h"
namespace {
using uint = unsigned int;
/* This is the normalized cardinal sine (sinc) function.
*
* sinc(x) = { 1, x = 0
* { sin(pi x) / (pi x), otherwise.
*/
constexpr double Sinc(const double x)
{
if(!(x > 1e-15 || x < -1e-15))
return 1.0;
return std::sin(al::MathDefs<double>::Pi()*x) / (al::MathDefs<double>::Pi()*x);
}
/* The zero-order modified Bessel function of the first kind, used for the
* Kaiser window.
*
* I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
* = sum_{k=0}^inf ((x / 2)^k / k!)^2
*/
constexpr double BesselI_0(const double x)
{
/* Start at k=1 since k=0 is trivial. */
const double x2{x / 2.0};
double term{1.0};
double sum{1.0};
double last_sum{};
int k{1};
/* Let the integration converge until the term of the sum is no longer
* significant.
*/
do {
const double y{x2 / k};
++k;
last_sum = sum;
term *= y * y;
sum += term;
} while(sum != last_sum);
return sum;
}
/* Calculate a Kaiser window from the given beta value and a normalized k
* [-1, 1].
*
* w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B), -1 <= k <= 1
* { 0, elsewhere.
*
* Where k can be calculated as:
*
* k = i / l, where -l <= i <= l.
*
* or:
*
* k = 2 i / M - 1, where 0 <= i <= M.
*/
constexpr double Kaiser(const double beta, const double k, const double besseli_0_beta)
{
if(!(k >= -1.0 && k <= 1.0))
return 0.0;
return BesselI_0(beta * std::sqrt(1.0 - k*k)) / besseli_0_beta;
}
/* Calculates the (normalized frequency) transition width of the Kaiser window.
* Rejection is in dB.
*/
constexpr double CalcKaiserWidth(const double rejection, const uint order)
{
if(rejection > 21.19)
return (rejection - 7.95) / (order * 2.285 * al::MathDefs<double>::Tau());
/* This enforces a minimum rejection of just above 21.18dB */
return 5.79 / (order * al::MathDefs<double>::Tau());
}
/* Calculates the beta value of the Kaiser window. Rejection is in dB. */
constexpr double CalcKaiserBeta(const double rejection)
{
if(rejection > 50.0)
return 0.1102 * (rejection-8.7);
else if(rejection >= 21.0)
return (0.5842 * std::pow(rejection-21.0, 0.4)) + (0.07886 * (rejection-21.0));
return 0.0;
}
struct BSincHeader {
double width{};
double beta{};
double scaleBase{};
double scaleRange{};
double besseli_0_beta{};
uint a[BSincScaleCount]{};
uint total_size{};
constexpr BSincHeader(uint Rejection, uint Order) noexcept
{
width = CalcKaiserWidth(Rejection, Order);
beta = CalcKaiserBeta(Rejection);
scaleBase = width / 2.0;
scaleRange = 1.0 - scaleBase;
besseli_0_beta = BesselI_0(beta);
uint num_points{Order+1};
for(uint si{0};si < BSincScaleCount;++si)
{
const double scale{scaleBase + (scaleRange * si / (BSincScaleCount-1))};
const uint a_{std::min(static_cast<uint>(num_points / 2.0 / scale), num_points)};
const uint m{2 * a_};
a[si] = a_;
total_size += 4 * BSincPhaseCount * ((m+3) & ~3u);
}
}
};
/* 11th and 23rd order filters (12 and 24-point respectively) with a 60dB drop
* at nyquist. Each filter will scale up the order when downsampling, to 23rd
* and 47th order respectively.
*/
constexpr BSincHeader bsinc12_hdr{60, 11};
constexpr BSincHeader bsinc24_hdr{60, 23};
/* NOTE: GCC 5 has an issue with BSincHeader objects being in an anonymous
* namespace while also being used as non-type template parameters.
*/
#if !defined(__clang__) && defined(__GNUC__) && __GNUC__ < 6
template<size_t total_size>
struct BSincFilterArray {
alignas(16) std::array<float, total_size> mTable;
BSincFilterArray(const BSincHeader &hdr)
#else
template<const BSincHeader &hdr>
struct BSincFilterArray {
alignas(16) std::array<float, hdr.total_size> mTable;
BSincFilterArray()
#endif
{
using filter_type = double[][BSincPhaseCount+1][BSincPointsMax];
auto filter = std::make_unique<filter_type>(BSincScaleCount);
/* Calculate the Kaiser-windowed Sinc filter coefficients for each
* scale and phase index.
*/
for(uint si{0};si < BSincScaleCount;++si)
{
const uint m{hdr.a[si] * 2};
const size_t o{(BSincPointsMax-m) / 2};
const double scale{hdr.scaleBase + (hdr.scaleRange * si / (BSincScaleCount-1))};
const double cutoff{scale - (hdr.scaleBase * std::max(0.5, scale) * 2.0)};
const auto a = static_cast<double>(hdr.a[si]);
const double l{a - 1.0};
/* Do one extra phase index so that the phase delta has a proper
* target for its last index.
*/
for(uint pi{0};pi <= BSincPhaseCount;++pi)
{
const double phase{l + (pi/double{BSincPhaseCount})};
for(uint i{0};i < m;++i)
{
const double x{i - phase};
filter[si][pi][o+i] = Kaiser(hdr.beta, x/a, hdr.besseli_0_beta) * cutoff *
Sinc(cutoff*x);
}
}
}
size_t idx{0};
for(size_t si{0};si < BSincScaleCount-1;++si)
{
const size_t m{((hdr.a[si]*2) + 3) & ~3u};
const size_t o{(BSincPointsMax-m) / 2};
for(size_t pi{0};pi < BSincPhaseCount;++pi)
{
/* Write out the filter. Also calculate and write out the phase
* and scale deltas.
*/
for(size_t i{0};i < m;++i)
mTable[idx++] = static_cast<float>(filter[si][pi][o+i]);
/* Linear interpolation between phases is simplified by pre-
* calculating the delta (b - a) in: x = a + f (b - a)
*/
for(size_t i{0};i < m;++i)
{
const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
mTable[idx++] = static_cast<float>(phDelta);
}
/* Linear interpolation between scales is also simplified.
*
* Given a difference in points between scales, the destination
* points will be 0, thus: x = a + f (-a)
*/
for(size_t i{0};i < m;++i)
{
const double scDelta{filter[si+1][pi][o+i] - filter[si][pi][o+i]};
mTable[idx++] = static_cast<float>(scDelta);
}
/* This last simplification is done to complete the bilinear
* equation for the combination of phase and scale.
*/
for(size_t i{0};i < m;++i)
{
const double spDelta{(filter[si+1][pi+1][o+i] - filter[si+1][pi][o+i]) -
(filter[si][pi+1][o+i] - filter[si][pi][o+i])};
mTable[idx++] = static_cast<float>(spDelta);
}
}
}
{
/* The last scale index doesn't have any scale or scale-phase
* deltas.
*/
constexpr size_t si{BSincScaleCount-1};
const size_t m{((hdr.a[si]*2) + 3) & ~3u};
const size_t o{(BSincPointsMax-m) / 2};
for(size_t pi{0};pi < BSincPhaseCount;++pi)
{
for(size_t i{0};i < m;++i)
mTable[idx++] = static_cast<float>(filter[si][pi][o+i]);
for(size_t i{0};i < m;++i)
{
const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
mTable[idx++] = static_cast<float>(phDelta);
}
for(size_t i{0};i < m;++i)
mTable[idx++] = 0.0f;
for(size_t i{0};i < m;++i)
mTable[idx++] = 0.0f;
}
}
assert(idx == hdr.total_size);
}
};
#if !defined(__clang__) && defined(__GNUC__) && __GNUC__ < 6
const BSincFilterArray<bsinc12_hdr.total_size> bsinc12_filter{bsinc12_hdr};
const BSincFilterArray<bsinc24_hdr.total_size> bsinc24_filter{bsinc24_hdr};
#else
const BSincFilterArray<bsinc12_hdr> bsinc12_filter{};
const BSincFilterArray<bsinc24_hdr> bsinc24_filter{};
#endif
constexpr BSincTable GenerateBSincTable(const BSincHeader &hdr, const float *tab)
{
BSincTable ret{};
ret.scaleBase = static_cast<float>(hdr.scaleBase);
ret.scaleRange = static_cast<float>(1.0 / hdr.scaleRange);
for(size_t i{0};i < BSincScaleCount;++i)
ret.m[i] = ((hdr.a[i]*2) + 3) & ~3u;
ret.filterOffset[0] = 0;
for(size_t i{1};i < BSincScaleCount;++i)
ret.filterOffset[i] = ret.filterOffset[i-1] + ret.m[i-1]*4*BSincPhaseCount;
ret.Tab = tab;
return ret;
}
} // namespace
const BSincTable bsinc12{GenerateBSincTable(bsinc12_hdr, &bsinc12_filter.mTable.front())};
const BSincTable bsinc24{GenerateBSincTable(bsinc24_hdr, &bsinc24_filter.mTable.front())};

View file

@ -0,0 +1,17 @@
#ifndef CORE_BSINC_TABLES_H
#define CORE_BSINC_TABLES_H
#include "bsinc_defs.h"
struct BSincTable {
float scaleBase, scaleRange;
unsigned int m[BSincScaleCount];
unsigned int filterOffset[BSincScaleCount];
const float *Tab;
};
extern const BSincTable bsinc12;
extern const BSincTable bsinc24;
#endif /* CORE_BSINC_TABLES_H */

View file

@ -0,0 +1,14 @@
#ifndef CORE_BUFFERLINE_H
#define CORE_BUFFERLINE_H
#include <array>
/* Size for temporary storage of buffer data, in floats. Larger values need
* more memory and are harder on cache, while smaller values may need more
* iterations for mixing.
*/
constexpr int BufferLineSize{1024};
using FloatBufferLine = std::array<float,BufferLineSize>;
#endif /* CORE_BUFFERLINE_H */

View file

@ -0,0 +1,142 @@
#include "config.h"
#include "cpu_caps.h"
#if defined(_WIN32) && (defined(_M_ARM) || defined(_M_ARM64))
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
#ifndef PF_ARM_NEON_INSTRUCTIONS_AVAILABLE
#define PF_ARM_NEON_INSTRUCTIONS_AVAILABLE 19
#endif
#endif
#ifdef HAVE_INTRIN_H
#include <intrin.h>
#endif
#ifdef HAVE_CPUID_H
#include <cpuid.h>
#endif
#include <array>
#include <cctype>
#include <string>
int CPUCapFlags{0};
namespace {
#if defined(HAVE_GCC_GET_CPUID) \
&& (defined(__i386__) || defined(__x86_64__) || defined(_M_IX86) || defined(_M_X64))
using reg_type = unsigned int;
inline std::array<reg_type,4> get_cpuid(unsigned int f)
{
std::array<reg_type,4> ret{};
__get_cpuid(f, &ret[0], &ret[1], &ret[2], &ret[3]);
return ret;
}
#define CAN_GET_CPUID
#elif defined(HAVE_CPUID_INTRINSIC) \
&& (defined(__i386__) || defined(__x86_64__) || defined(_M_IX86) || defined(_M_X64))
using reg_type = int;
inline std::array<reg_type,4> get_cpuid(unsigned int f)
{
std::array<reg_type,4> ret{};
(__cpuid)(ret.data(), f);
return ret;
}
#define CAN_GET_CPUID
#endif
} // namespace
al::optional<CPUInfo> GetCPUInfo()
{
CPUInfo ret;
#ifdef CAN_GET_CPUID
auto cpuregs = get_cpuid(0);
if(cpuregs[0] == 0)
return al::nullopt;
const reg_type maxfunc{cpuregs[0]};
cpuregs = get_cpuid(0x80000000);
const reg_type maxextfunc{cpuregs[0]};
ret.mVendor.append(reinterpret_cast<char*>(&cpuregs[1]), 4);
ret.mVendor.append(reinterpret_cast<char*>(&cpuregs[3]), 4);
ret.mVendor.append(reinterpret_cast<char*>(&cpuregs[2]), 4);
auto iter_end = std::remove(ret.mVendor.begin(), ret.mVendor.end(), '\0');
iter_end = std::unique(ret.mVendor.begin(), iter_end,
[](auto&& c0, auto&& c1) { return std::isspace(c0) && std::isspace(c1); });
ret.mVendor.erase(iter_end, ret.mVendor.end());
if(!ret.mVendor.empty() && std::isspace(ret.mVendor.back()))
ret.mVendor.pop_back();
if(!ret.mVendor.empty() && std::isspace(ret.mVendor.front()))
ret.mVendor.erase(ret.mVendor.begin());
if(maxextfunc >= 0x80000004)
{
cpuregs = get_cpuid(0x80000002);
ret.mName.append(reinterpret_cast<char*>(cpuregs.data()), 16);
cpuregs = get_cpuid(0x80000003);
ret.mName.append(reinterpret_cast<char*>(cpuregs.data()), 16);
cpuregs = get_cpuid(0x80000004);
ret.mName.append(reinterpret_cast<char*>(cpuregs.data()), 16);
iter_end = std::remove(ret.mName.begin(), ret.mName.end(), '\0');
iter_end = std::unique(ret.mName.begin(), iter_end,
[](auto&& c0, auto&& c1) { return std::isspace(c0) && std::isspace(c1); });
ret.mName.erase(iter_end, ret.mName.end());
if(!ret.mName.empty() && std::isspace(ret.mName.back()))
ret.mName.pop_back();
if(!ret.mName.empty() && std::isspace(ret.mName.front()))
ret.mName.erase(ret.mName.begin());
}
if(maxfunc >= 1)
{
cpuregs = get_cpuid(1);
if((cpuregs[3]&(1<<25)))
ret.mCaps |= CPU_CAP_SSE;
if((ret.mCaps&CPU_CAP_SSE) && (cpuregs[3]&(1<<26)))
ret.mCaps |= CPU_CAP_SSE2;
if((ret.mCaps&CPU_CAP_SSE2) && (cpuregs[2]&(1<<0)))
ret.mCaps |= CPU_CAP_SSE3;
if((ret.mCaps&CPU_CAP_SSE3) && (cpuregs[2]&(1<<19)))
ret.mCaps |= CPU_CAP_SSE4_1;
}
#else
/* Assume support for whatever's supported if we can't check for it */
#if defined(HAVE_SSE4_1)
#warning "Assuming SSE 4.1 run-time support!"
ret.mCaps |= CPU_CAP_SSE | CPU_CAP_SSE2 | CPU_CAP_SSE3 | CPU_CAP_SSE4_1;
#elif defined(HAVE_SSE3)
#warning "Assuming SSE 3 run-time support!"
ret.mCaps |= CPU_CAP_SSE | CPU_CAP_SSE2 | CPU_CAP_SSE3;
#elif defined(HAVE_SSE2)
#warning "Assuming SSE 2 run-time support!"
ret.mCaps |= CPU_CAP_SSE | CPU_CAP_SSE2;
#elif defined(HAVE_SSE)
#warning "Assuming SSE run-time support!"
ret.mCaps |= CPU_CAP_SSE;
#endif
#endif /* CAN_GET_CPUID */
#ifdef HAVE_NEON
#ifdef __ARM_NEON
ret.mCaps |= CPU_CAP_NEON;
#elif defined(_WIN32) && (defined(_M_ARM) || defined(_M_ARM64))
if(IsProcessorFeaturePresent(PF_ARM_NEON_INSTRUCTIONS_AVAILABLE))
ret.mCaps |= CPU_CAP_NEON;
#else
#warning "Assuming NEON run-time support!"
ret.mCaps |= CPU_CAP_NEON;
#endif
#endif
return al::make_optional(ret);
}

View file

@ -0,0 +1,26 @@
#ifndef CORE_CPU_CAPS_H
#define CORE_CPU_CAPS_H
#include <string>
#include "aloptional.h"
extern int CPUCapFlags;
enum {
CPU_CAP_SSE = 1<<0,
CPU_CAP_SSE2 = 1<<1,
CPU_CAP_SSE3 = 1<<2,
CPU_CAP_SSE4_1 = 1<<3,
CPU_CAP_NEON = 1<<4,
};
struct CPUInfo {
std::string mVendor;
std::string mName;
int mCaps{0};
};
al::optional<CPUInfo> GetCPUInfo();
#endif /* CORE_CPU_CAPS_H */

View file

@ -0,0 +1,65 @@
#include "config.h"
#include "devformat.h"
uint BytesFromDevFmt(DevFmtType type) noexcept
{
switch(type)
{
case DevFmtByte: return sizeof(int8_t);
case DevFmtUByte: return sizeof(uint8_t);
case DevFmtShort: return sizeof(int16_t);
case DevFmtUShort: return sizeof(uint16_t);
case DevFmtInt: return sizeof(int32_t);
case DevFmtUInt: return sizeof(uint32_t);
case DevFmtFloat: return sizeof(float);
}
return 0;
}
uint ChannelsFromDevFmt(DevFmtChannels chans, uint ambiorder) noexcept
{
switch(chans)
{
case DevFmtMono: return 1;
case DevFmtStereo: return 2;
case DevFmtQuad: return 4;
case DevFmtX51: return 6;
case DevFmtX51Rear: return 6;
case DevFmtX61: return 7;
case DevFmtX71: return 8;
case DevFmtAmbi3D: return (ambiorder+1) * (ambiorder+1);
}
return 0;
}
const char *DevFmtTypeString(DevFmtType type) noexcept
{
switch(type)
{
case DevFmtByte: return "Int8";
case DevFmtUByte: return "UInt8";
case DevFmtShort: return "Int16";
case DevFmtUShort: return "UInt16";
case DevFmtInt: return "Int32";
case DevFmtUInt: return "UInt32";
case DevFmtFloat: return "Float32";
}
return "(unknown type)";
}
const char *DevFmtChannelsString(DevFmtChannels chans) noexcept
{
switch(chans)
{
case DevFmtMono: return "Mono";
case DevFmtStereo: return "Stereo";
case DevFmtQuad: return "Quadraphonic";
case DevFmtX51: return "5.1 Surround";
case DevFmtX51Rear: return "5.1 Surround (Rear)";
case DevFmtX61: return "6.1 Surround";
case DevFmtX71: return "7.1 Surround";
case DevFmtAmbi3D: return "Ambisonic 3D";
}
return "(unknown channels)";
}

View file

@ -0,0 +1,106 @@
#ifndef CORE_DEVFORMAT_H
#define CORE_DEVFORMAT_H
#include <cstdint>
using uint = unsigned int;
enum Channel : unsigned char {
FrontLeft = 0,
FrontRight,
FrontCenter,
LFE,
BackLeft,
BackRight,
BackCenter,
SideLeft,
SideRight,
TopFrontLeft,
TopFrontCenter,
TopFrontRight,
TopCenter,
TopBackLeft,
TopBackCenter,
TopBackRight,
MaxChannels
};
/* Device formats */
enum DevFmtType : unsigned char {
DevFmtByte,
DevFmtUByte,
DevFmtShort,
DevFmtUShort,
DevFmtInt,
DevFmtUInt,
DevFmtFloat,
DevFmtTypeDefault = DevFmtFloat
};
enum DevFmtChannels : unsigned char {
DevFmtMono,
DevFmtStereo,
DevFmtQuad,
DevFmtX51,
DevFmtX61,
DevFmtX71,
DevFmtAmbi3D,
/* Similar to 5.1, except using rear channels instead of sides */
DevFmtX51Rear,
DevFmtChannelsDefault = DevFmtStereo
};
#define MAX_OUTPUT_CHANNELS 16
/* DevFmtType traits, providing the type, etc given a DevFmtType. */
template<DevFmtType T>
struct DevFmtTypeTraits { };
template<>
struct DevFmtTypeTraits<DevFmtByte> { using Type = int8_t; };
template<>
struct DevFmtTypeTraits<DevFmtUByte> { using Type = uint8_t; };
template<>
struct DevFmtTypeTraits<DevFmtShort> { using Type = int16_t; };
template<>
struct DevFmtTypeTraits<DevFmtUShort> { using Type = uint16_t; };
template<>
struct DevFmtTypeTraits<DevFmtInt> { using Type = int32_t; };
template<>
struct DevFmtTypeTraits<DevFmtUInt> { using Type = uint32_t; };
template<>
struct DevFmtTypeTraits<DevFmtFloat> { using Type = float; };
template<DevFmtType T>
using DevFmtType_t = typename DevFmtTypeTraits<T>::Type;
uint BytesFromDevFmt(DevFmtType type) noexcept;
uint ChannelsFromDevFmt(DevFmtChannels chans, uint ambiorder) noexcept;
inline uint FrameSizeFromDevFmt(DevFmtChannels chans, DevFmtType type, uint ambiorder) noexcept
{ return ChannelsFromDevFmt(chans, ambiorder) * BytesFromDevFmt(type); }
const char *DevFmtTypeString(DevFmtType type) noexcept;
const char *DevFmtChannelsString(DevFmtChannels chans) noexcept;
enum class DevAmbiLayout : bool {
FuMa,
ACN,
Default = ACN
};
enum class DevAmbiScaling : unsigned char {
FuMa,
SN3D,
N3D,
Default = SN3D
};
#endif /* CORE_DEVFORMAT_H */

View file

@ -0,0 +1,31 @@
#include "config.h"
#include "except.h"
#include <cstdio>
#include <cstdarg>
#include "opthelpers.h"
namespace al {
/* Defined here to avoid inlining it. */
base_exception::~base_exception() { }
void base_exception::setMessage(const char* msg, std::va_list args)
{
std::va_list args2;
va_copy(args2, args);
int msglen{std::vsnprintf(nullptr, 0, msg, args)};
if LIKELY(msglen > 0)
{
mMessage.resize(static_cast<size_t>(msglen)+1);
std::vsnprintf(&mMessage[0], mMessage.length(), msg, args2);
mMessage.pop_back();
}
va_end(args2);
}
} // namespace al

View file

@ -0,0 +1,31 @@
#ifndef CORE_EXCEPT_H
#define CORE_EXCEPT_H
#include <cstdarg>
#include <exception>
#include <string>
#include <utility>
namespace al {
class base_exception : public std::exception {
std::string mMessage;
protected:
base_exception() = default;
virtual ~base_exception();
void setMessage(const char *msg, std::va_list args);
public:
const char *what() const noexcept override { return mMessage.c_str(); }
};
} // namespace al
#define START_API_FUNC try
#define END_API_FUNC catch(...) { std::terminate(); }
#endif /* CORE_EXCEPT_H */

View file

@ -0,0 +1,163 @@
#include "config.h"
#include "biquad.h"
#include <algorithm>
#include <cassert>
#include <cmath>
#include "opthelpers.h"
template<typename Real>
void BiquadFilterR<Real>::setParams(BiquadType type, Real f0norm, Real gain, Real rcpQ)
{
// Limit gain to -100dB
assert(gain > 0.00001f);
const Real w0{al::MathDefs<Real>::Tau() * f0norm};
const Real sin_w0{std::sin(w0)};
const Real cos_w0{std::cos(w0)};
const Real alpha{sin_w0/2.0f * rcpQ};
Real sqrtgain_alpha_2;
Real a[3]{ 1.0f, 0.0f, 0.0f };
Real b[3]{ 1.0f, 0.0f, 0.0f };
/* Calculate filter coefficients depending on filter type */
switch(type)
{
case BiquadType::HighShelf:
sqrtgain_alpha_2 = 2.0f * std::sqrt(gain) * alpha;
b[0] = gain*((gain+1.0f) + (gain-1.0f)*cos_w0 + sqrtgain_alpha_2);
b[1] = -2.0f*gain*((gain-1.0f) + (gain+1.0f)*cos_w0 );
b[2] = gain*((gain+1.0f) + (gain-1.0f)*cos_w0 - sqrtgain_alpha_2);
a[0] = (gain+1.0f) - (gain-1.0f)*cos_w0 + sqrtgain_alpha_2;
a[1] = 2.0f* ((gain-1.0f) - (gain+1.0f)*cos_w0 );
a[2] = (gain+1.0f) - (gain-1.0f)*cos_w0 - sqrtgain_alpha_2;
break;
case BiquadType::LowShelf:
sqrtgain_alpha_2 = 2.0f * std::sqrt(gain) * alpha;
b[0] = gain*((gain+1.0f) - (gain-1.0f)*cos_w0 + sqrtgain_alpha_2);
b[1] = 2.0f*gain*((gain-1.0f) - (gain+1.0f)*cos_w0 );
b[2] = gain*((gain+1.0f) - (gain-1.0f)*cos_w0 - sqrtgain_alpha_2);
a[0] = (gain+1.0f) + (gain-1.0f)*cos_w0 + sqrtgain_alpha_2;
a[1] = -2.0f* ((gain-1.0f) + (gain+1.0f)*cos_w0 );
a[2] = (gain+1.0f) + (gain-1.0f)*cos_w0 - sqrtgain_alpha_2;
break;
case BiquadType::Peaking:
b[0] = 1.0f + alpha * gain;
b[1] = -2.0f * cos_w0;
b[2] = 1.0f - alpha * gain;
a[0] = 1.0f + alpha / gain;
a[1] = -2.0f * cos_w0;
a[2] = 1.0f - alpha / gain;
break;
case BiquadType::LowPass:
b[0] = (1.0f - cos_w0) / 2.0f;
b[1] = 1.0f - cos_w0;
b[2] = (1.0f - cos_w0) / 2.0f;
a[0] = 1.0f + alpha;
a[1] = -2.0f * cos_w0;
a[2] = 1.0f - alpha;
break;
case BiquadType::HighPass:
b[0] = (1.0f + cos_w0) / 2.0f;
b[1] = -(1.0f + cos_w0);
b[2] = (1.0f + cos_w0) / 2.0f;
a[0] = 1.0f + alpha;
a[1] = -2.0f * cos_w0;
a[2] = 1.0f - alpha;
break;
case BiquadType::BandPass:
b[0] = alpha;
b[1] = 0.0f;
b[2] = -alpha;
a[0] = 1.0f + alpha;
a[1] = -2.0f * cos_w0;
a[2] = 1.0f - alpha;
break;
}
mA1 = a[1] / a[0];
mA2 = a[2] / a[0];
mB0 = b[0] / a[0];
mB1 = b[1] / a[0];
mB2 = b[2] / a[0];
}
template<typename Real>
void BiquadFilterR<Real>::process(const al::span<const Real> src, Real *dst)
{
const Real b0{mB0};
const Real b1{mB1};
const Real b2{mB2};
const Real a1{mA1};
const Real a2{mA2};
Real z1{mZ1};
Real z2{mZ2};
/* Processing loop is Transposed Direct Form II. This requires less storage
* compared to Direct Form I (only two delay components, instead of a four-
* sample history; the last two inputs and outputs), and works better for
* floating-point which favors summing similarly-sized values while being
* less bothered by overflow.
*
* See: http://www.earlevel.com/main/2003/02/28/biquads/
*/
auto proc_sample = [b0,b1,b2,a1,a2,&z1,&z2](Real input) noexcept -> Real
{
const Real output{input*b0 + z1};
z1 = input*b1 - output*a1 + z2;
z2 = input*b2 - output*a2;
return output;
};
std::transform(src.cbegin(), src.cend(), dst, proc_sample);
mZ1 = z1;
mZ2 = z2;
}
template<typename Real>
void BiquadFilterR<Real>::dualProcess(BiquadFilterR &other, const al::span<const Real> src,
Real *dst)
{
const Real b00{mB0};
const Real b01{mB1};
const Real b02{mB2};
const Real a01{mA1};
const Real a02{mA2};
const Real b10{other.mB0};
const Real b11{other.mB1};
const Real b12{other.mB2};
const Real a11{other.mA1};
const Real a12{other.mA2};
Real z01{mZ1};
Real z02{mZ2};
Real z11{other.mZ1};
Real z12{other.mZ2};
auto proc_sample = [b00,b01,b02,a01,a02,b10,b11,b12,a11,a12,&z01,&z02,&z11,&z12](Real input) noexcept -> Real
{
const Real tmpout{input*b00 + z01};
z01 = input*b01 - tmpout*a01 + z02;
z02 = input*b02 - tmpout*a02;
input = tmpout;
const Real output{input*b10 + z11};
z11 = input*b11 - output*a11 + z12;
z12 = input*b12 - output*a12;
return output;
};
std::transform(src.cbegin(), src.cend(), dst, proc_sample);
mZ1 = z01;
mZ2 = z02;
other.mZ1 = z11;
other.mZ2 = z12;
}
template class BiquadFilterR<float>;
template class BiquadFilterR<double>;

View file

@ -0,0 +1,144 @@
#ifndef CORE_FILTERS_BIQUAD_H
#define CORE_FILTERS_BIQUAD_H
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <utility>
#include "alspan.h"
#include "math_defs.h"
/* Filters implementation is based on the "Cookbook formulae for audio
* EQ biquad filter coefficients" by Robert Bristow-Johnson
* http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
*/
/* Implementation note: For the shelf and peaking filters, the specified gain
* is for the centerpoint of the transition band. This better fits EFX filter
* behavior, which expects the shelf's reference frequency to reach the given
* gain. To set the gain for the shelf or peak itself, use the square root of
* the desired linear gain (or halve the dB gain).
*/
enum class BiquadType {
/** EFX-style low-pass filter, specifying a gain and reference frequency. */
HighShelf,
/** EFX-style high-pass filter, specifying a gain and reference frequency. */
LowShelf,
/** Peaking filter, specifying a gain and reference frequency. */
Peaking,
/** Low-pass cut-off filter, specifying a cut-off frequency. */
LowPass,
/** High-pass cut-off filter, specifying a cut-off frequency. */
HighPass,
/** Band-pass filter, specifying a center frequency. */
BandPass,
};
template<typename Real>
class BiquadFilterR {
/* Last two delayed components for direct form II. */
Real mZ1{0.0f}, mZ2{0.0f};
/* Transfer function coefficients "b" (numerator) */
Real mB0{1.0f}, mB1{0.0f}, mB2{0.0f};
/* Transfer function coefficients "a" (denominator; a0 is pre-applied). */
Real mA1{0.0f}, mA2{0.0f};
void setParams(BiquadType type, Real f0norm, Real gain, Real rcpQ);
/**
* Calculates the rcpQ (i.e. 1/Q) coefficient for shelving filters, using
* the reference gain and shelf slope parameter.
* \param gain 0 < gain
* \param slope 0 < slope <= 1
*/
static Real rcpQFromSlope(Real gain, Real slope)
{ return std::sqrt((gain + 1.0f/gain)*(1.0f/slope - 1.0f) + 2.0f); }
/**
* Calculates the rcpQ (i.e. 1/Q) coefficient for filters, using the
* normalized reference frequency and bandwidth.
* \param f0norm 0 < f0norm < 0.5.
* \param bandwidth 0 < bandwidth
*/
static Real rcpQFromBandwidth(Real f0norm, Real bandwidth)
{
const Real w0{al::MathDefs<Real>::Tau() * f0norm};
return 2.0f*std::sinh(std::log(Real{2.0f})/2.0f*bandwidth*w0/std::sin(w0));
}
public:
void clear() noexcept { mZ1 = mZ2 = 0.0f; }
/**
* Sets the filter state for the specified filter type and its parameters.
*
* \param type The type of filter to apply.
* \param f0norm The normalized reference frequency (ref / sample_rate).
* This is the center point for the Shelf, Peaking, and BandPass filter
* types, or the cutoff frequency for the LowPass and HighPass filter
* types.
* \param gain The gain for the reference frequency response. Only used by
* the Shelf and Peaking filter types.
* \param slope Slope steepness of the transition band.
*/
void setParamsFromSlope(BiquadType type, Real f0norm, Real gain, Real slope)
{
gain = std::max<Real>(gain, 0.001f); /* Limit -60dB */
setParams(type, f0norm, gain, rcpQFromSlope(gain, slope));
}
/**
* Sets the filter state for the specified filter type and its parameters.
*
* \param type The type of filter to apply.
* \param f0norm The normalized reference frequency (ref / sample_rate).
* This is the center point for the Shelf, Peaking, and BandPass filter
* types, or the cutoff frequency for the LowPass and HighPass filter
* types.
* \param gain The gain for the reference frequency response. Only used by
* the Shelf and Peaking filter types.
* \param bandwidth Normalized bandwidth of the transition band.
*/
void setParamsFromBandwidth(BiquadType type, Real f0norm, Real gain, Real bandwidth)
{ setParams(type, f0norm, gain, rcpQFromBandwidth(f0norm, bandwidth)); }
void copyParamsFrom(const BiquadFilterR &other)
{
mB0 = other.mB0;
mB1 = other.mB1;
mB2 = other.mB2;
mA1 = other.mA1;
mA2 = other.mA2;
}
void process(const al::span<const Real> src, Real *dst);
/** Processes this filter and the other at the same time. */
void dualProcess(BiquadFilterR &other, const al::span<const Real> src, Real *dst);
/* Rather hacky. It's just here to support "manual" processing. */
std::pair<Real,Real> getComponents() const noexcept { return {mZ1, mZ2}; }
void setComponents(Real z1, Real z2) noexcept { mZ1 = z1; mZ2 = z2; }
Real processOne(const Real in, Real &z1, Real &z2) const noexcept
{
const Real out{in*mB0 + z1};
z1 = in*mB1 - out*mA1 + z2;
z2 = in*mB2 - out*mA2;
return out;
}
};
template<typename Real>
struct DualBiquadR {
BiquadFilterR<Real> &f0, &f1;
void process(const al::span<const Real> src, Real *dst)
{ f0.dualProcess(f1, src, dst); }
};
using BiquadFilter = BiquadFilterR<float>;
using DualBiquad = DualBiquadR<float>;
#endif /* CORE_FILTERS_BIQUAD_H */

View file

@ -0,0 +1,383 @@
#include "config.h"
#include "nfc.h"
#include <algorithm>
#include "opthelpers.h"
/* Near-field control filters are the basis for handling the near-field effect.
* The near-field effect is a bass-boost present in the directional components
* of a recorded signal, created as a result of the wavefront curvature (itself
* a function of sound distance). Proper reproduction dictates this be
* compensated for using a bass-cut given the playback speaker distance, to
* avoid excessive bass in the playback.
*
* For real-time rendered audio, emulating the near-field effect based on the
* sound source's distance, and subsequently compensating for it at output
* based on the speaker distances, can create a more realistic perception of
* sound distance beyond a simple 1/r attenuation.
*
* These filters do just that. Each one applies a low-shelf filter, created as
* the combination of a bass-boost for a given sound source distance (near-
* field emulation) along with a bass-cut for a given control/speaker distance
* (near-field compensation).
*
* Note that it is necessary to apply a cut along with the boost, since the
* boost alone is unstable in higher-order ambisonics as it causes an infinite
* DC gain (even first-order ambisonics requires there to be no DC offset for
* the boost to work). Consequently, ambisonics requires a control parameter to
* be used to avoid an unstable boost-only filter. NFC-HOA defines this control
* as a reference delay, calculated with:
*
* reference_delay = control_distance / speed_of_sound
*
* This means w0 (for input) or w1 (for output) should be set to:
*
* wN = 1 / (reference_delay * sample_rate)
*
* when dealing with NFC-HOA content. For FOA input content, which does not
* specify a reference_delay variable, w0 should be set to 0 to apply only
* near-field compensation for output. It's important that w1 be a finite,
* positive, non-0 value or else the bass-boost will become unstable again.
* Also, w0 should not be too large compared to w1, to avoid excessively loud
* low frequencies.
*/
namespace {
constexpr float B[5][4] = {
{ 0.0f },
{ 1.0f },
{ 3.0f, 3.0f },
{ 3.6778f, 6.4595f, 2.3222f },
{ 4.2076f, 11.4877f, 5.7924f, 9.1401f }
};
NfcFilter1 NfcFilterCreate1(const float w0, const float w1) noexcept
{
NfcFilter1 nfc{};
float b_00, g_0;
float r;
nfc.base_gain = 1.0f;
nfc.gain = 1.0f;
/* Calculate bass-boost coefficients. */
r = 0.5f * w0;
b_00 = B[1][0] * r;
g_0 = 1.0f + b_00;
nfc.gain *= g_0;
nfc.b1 = 2.0f * b_00 / g_0;
/* Calculate bass-cut coefficients. */
r = 0.5f * w1;
b_00 = B[1][0] * r;
g_0 = 1.0f + b_00;
nfc.base_gain /= g_0;
nfc.gain /= g_0;
nfc.a1 = 2.0f * b_00 / g_0;
return nfc;
}
void NfcFilterAdjust1(NfcFilter1 *nfc, const float w0) noexcept
{
const float r{0.5f * w0};
const float b_00{B[1][0] * r};
const float g_0{1.0f + b_00};
nfc->gain = nfc->base_gain * g_0;
nfc->b1 = 2.0f * b_00 / g_0;
}
NfcFilter2 NfcFilterCreate2(const float w0, const float w1) noexcept
{
NfcFilter2 nfc{};
float b_10, b_11, g_1;
float r;
nfc.base_gain = 1.0f;
nfc.gain = 1.0f;
/* Calculate bass-boost coefficients. */
r = 0.5f * w0;
b_10 = B[2][0] * r;
b_11 = B[2][1] * r * r;
g_1 = 1.0f + b_10 + b_11;
nfc.gain *= g_1;
nfc.b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc.b2 = 4.0f * b_11 / g_1;
/* Calculate bass-cut coefficients. */
r = 0.5f * w1;
b_10 = B[2][0] * r;
b_11 = B[2][1] * r * r;
g_1 = 1.0f + b_10 + b_11;
nfc.base_gain /= g_1;
nfc.gain /= g_1;
nfc.a1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc.a2 = 4.0f * b_11 / g_1;
return nfc;
}
void NfcFilterAdjust2(NfcFilter2 *nfc, const float w0) noexcept
{
const float r{0.5f * w0};
const float b_10{B[2][0] * r};
const float b_11{B[2][1] * r * r};
const float g_1{1.0f + b_10 + b_11};
nfc->gain = nfc->base_gain * g_1;
nfc->b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc->b2 = 4.0f * b_11 / g_1;
}
NfcFilter3 NfcFilterCreate3(const float w0, const float w1) noexcept
{
NfcFilter3 nfc{};
float b_10, b_11, g_1;
float b_00, g_0;
float r;
nfc.base_gain = 1.0f;
nfc.gain = 1.0f;
/* Calculate bass-boost coefficients. */
r = 0.5f * w0;
b_10 = B[3][0] * r;
b_11 = B[3][1] * r * r;
b_00 = B[3][2] * r;
g_1 = 1.0f + b_10 + b_11;
g_0 = 1.0f + b_00;
nfc.gain *= g_1 * g_0;
nfc.b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc.b2 = 4.0f * b_11 / g_1;
nfc.b3 = 2.0f * b_00 / g_0;
/* Calculate bass-cut coefficients. */
r = 0.5f * w1;
b_10 = B[3][0] * r;
b_11 = B[3][1] * r * r;
b_00 = B[3][2] * r;
g_1 = 1.0f + b_10 + b_11;
g_0 = 1.0f + b_00;
nfc.base_gain /= g_1 * g_0;
nfc.gain /= g_1 * g_0;
nfc.a1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc.a2 = 4.0f * b_11 / g_1;
nfc.a3 = 2.0f * b_00 / g_0;
return nfc;
}
void NfcFilterAdjust3(NfcFilter3 *nfc, const float w0) noexcept
{
const float r{0.5f * w0};
const float b_10{B[3][0] * r};
const float b_11{B[3][1] * r * r};
const float b_00{B[3][2] * r};
const float g_1{1.0f + b_10 + b_11};
const float g_0{1.0f + b_00};
nfc->gain = nfc->base_gain * g_1 * g_0;
nfc->b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc->b2 = 4.0f * b_11 / g_1;
nfc->b3 = 2.0f * b_00 / g_0;
}
NfcFilter4 NfcFilterCreate4(const float w0, const float w1) noexcept
{
NfcFilter4 nfc{};
float b_10, b_11, g_1;
float b_00, b_01, g_0;
float r;
nfc.base_gain = 1.0f;
nfc.gain = 1.0f;
/* Calculate bass-boost coefficients. */
r = 0.5f * w0;
b_10 = B[4][0] * r;
b_11 = B[4][1] * r * r;
b_00 = B[4][2] * r;
b_01 = B[4][3] * r * r;
g_1 = 1.0f + b_10 + b_11;
g_0 = 1.0f + b_00 + b_01;
nfc.gain *= g_1 * g_0;
nfc.b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc.b2 = 4.0f * b_11 / g_1;
nfc.b3 = (2.0f*b_00 + 4.0f*b_01) / g_0;
nfc.b4 = 4.0f * b_01 / g_0;
/* Calculate bass-cut coefficients. */
r = 0.5f * w1;
b_10 = B[4][0] * r;
b_11 = B[4][1] * r * r;
b_00 = B[4][2] * r;
b_01 = B[4][3] * r * r;
g_1 = 1.0f + b_10 + b_11;
g_0 = 1.0f + b_00 + b_01;
nfc.base_gain /= g_1 * g_0;
nfc.gain /= g_1 * g_0;
nfc.a1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc.a2 = 4.0f * b_11 / g_1;
nfc.a3 = (2.0f*b_00 + 4.0f*b_01) / g_0;
nfc.a4 = 4.0f * b_01 / g_0;
return nfc;
}
void NfcFilterAdjust4(NfcFilter4 *nfc, const float w0) noexcept
{
const float r{0.5f * w0};
const float b_10{B[4][0] * r};
const float b_11{B[4][1] * r * r};
const float b_00{B[4][2] * r};
const float b_01{B[4][3] * r * r};
const float g_1{1.0f + b_10 + b_11};
const float g_0{1.0f + b_00 + b_01};
nfc->gain = nfc->base_gain * g_1 * g_0;
nfc->b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
nfc->b2 = 4.0f * b_11 / g_1;
nfc->b3 = (2.0f*b_00 + 4.0f*b_01) / g_0;
nfc->b4 = 4.0f * b_01 / g_0;
}
} // namespace
void NfcFilter::init(const float w1) noexcept
{
first = NfcFilterCreate1(0.0f, w1);
second = NfcFilterCreate2(0.0f, w1);
third = NfcFilterCreate3(0.0f, w1);
fourth = NfcFilterCreate4(0.0f, w1);
}
void NfcFilter::adjust(const float w0) noexcept
{
NfcFilterAdjust1(&first, w0);
NfcFilterAdjust2(&second, w0);
NfcFilterAdjust3(&third, w0);
NfcFilterAdjust4(&fourth, w0);
}
void NfcFilter::process1(const al::span<const float> src, float *RESTRICT dst)
{
const float gain{first.gain};
const float b1{first.b1};
const float a1{first.a1};
float z1{first.z[0]};
auto proc_sample = [gain,b1,a1,&z1](const float in) noexcept -> float
{
const float y{in*gain - a1*z1};
const float out{y + b1*z1};
z1 += y;
return out;
};
std::transform(src.cbegin(), src.cend(), dst, proc_sample);
first.z[0] = z1;
}
void NfcFilter::process2(const al::span<const float> src, float *RESTRICT dst)
{
const float gain{second.gain};
const float b1{second.b1};
const float b2{second.b2};
const float a1{second.a1};
const float a2{second.a2};
float z1{second.z[0]};
float z2{second.z[1]};
auto proc_sample = [gain,b1,b2,a1,a2,&z1,&z2](const float in) noexcept -> float
{
const float y{in*gain - a1*z1 - a2*z2};
const float out{y + b1*z1 + b2*z2};
z2 += z1;
z1 += y;
return out;
};
std::transform(src.cbegin(), src.cend(), dst, proc_sample);
second.z[0] = z1;
second.z[1] = z2;
}
void NfcFilter::process3(const al::span<const float> src, float *RESTRICT dst)
{
const float gain{third.gain};
const float b1{third.b1};
const float b2{third.b2};
const float b3{third.b3};
const float a1{third.a1};
const float a2{third.a2};
const float a3{third.a3};
float z1{third.z[0]};
float z2{third.z[1]};
float z3{third.z[2]};
auto proc_sample = [gain,b1,b2,b3,a1,a2,a3,&z1,&z2,&z3](const float in) noexcept -> float
{
float y{in*gain - a1*z1 - a2*z2};
float out{y + b1*z1 + b2*z2};
z2 += z1;
z1 += y;
y = out - a3*z3;
out = y + b3*z3;
z3 += y;
return out;
};
std::transform(src.cbegin(), src.cend(), dst, proc_sample);
third.z[0] = z1;
third.z[1] = z2;
third.z[2] = z3;
}
void NfcFilter::process4(const al::span<const float> src, float *RESTRICT dst)
{
const float gain{fourth.gain};
const float b1{fourth.b1};
const float b2{fourth.b2};
const float b3{fourth.b3};
const float b4{fourth.b4};
const float a1{fourth.a1};
const float a2{fourth.a2};
const float a3{fourth.a3};
const float a4{fourth.a4};
float z1{fourth.z[0]};
float z2{fourth.z[1]};
float z3{fourth.z[2]};
float z4{fourth.z[3]};
auto proc_sample = [gain,b1,b2,b3,b4,a1,a2,a3,a4,&z1,&z2,&z3,&z4](const float in) noexcept -> float
{
float y{in*gain - a1*z1 - a2*z2};
float out{y + b1*z1 + b2*z2};
z2 += z1;
z1 += y;
y = out - a3*z3 - a4*z4;
out = y + b3*z3 + b4*z4;
z4 += z3;
z3 += y;
return out;
};
std::transform(src.cbegin(), src.cend(), dst, proc_sample);
fourth.z[0] = z1;
fourth.z[1] = z2;
fourth.z[2] = z3;
fourth.z[3] = z4;
}

View file

@ -0,0 +1,63 @@
#ifndef CORE_FILTERS_NFC_H
#define CORE_FILTERS_NFC_H
#include <cstddef>
#include "alspan.h"
struct NfcFilter1 {
float base_gain, gain;
float b1, a1;
float z[1];
};
struct NfcFilter2 {
float base_gain, gain;
float b1, b2, a1, a2;
float z[2];
};
struct NfcFilter3 {
float base_gain, gain;
float b1, b2, b3, a1, a2, a3;
float z[3];
};
struct NfcFilter4 {
float base_gain, gain;
float b1, b2, b3, b4, a1, a2, a3, a4;
float z[4];
};
class NfcFilter {
NfcFilter1 first;
NfcFilter2 second;
NfcFilter3 third;
NfcFilter4 fourth;
public:
/* NOTE:
* w0 = speed_of_sound / (source_distance * sample_rate);
* w1 = speed_of_sound / (control_distance * sample_rate);
*
* Generally speaking, the control distance should be approximately the
* average speaker distance, or based on the reference delay if outputing
* NFC-HOA. It must not be negative, 0, or infinite. The source distance
* should not be too small relative to the control distance.
*/
void init(const float w1) noexcept;
void adjust(const float w0) noexcept;
/* Near-field control filter for first-order ambisonic channels (1-3). */
void process1(const al::span<const float> src, float *RESTRICT dst);
/* Near-field control filter for second-order ambisonic channels (4-8). */
void process2(const al::span<const float> src, float *RESTRICT dst);
/* Near-field control filter for third-order ambisonic channels (9-15). */
void process3(const al::span<const float> src, float *RESTRICT dst);
/* Near-field control filter for fourth-order ambisonic channels (16-24). */
void process4(const al::span<const float> src, float *RESTRICT dst);
};
#endif /* CORE_FILTERS_NFC_H */

View file

@ -0,0 +1,113 @@
#include "config.h"
#include "splitter.h"
#include <algorithm>
#include <cmath>
#include <limits>
#include "math_defs.h"
#include "opthelpers.h"
template<typename Real>
void BandSplitterR<Real>::init(Real f0norm)
{
const Real w{f0norm * al::MathDefs<Real>::Tau()};
const Real cw{std::cos(w)};
if(cw > std::numeric_limits<float>::epsilon())
mCoeff = (std::sin(w) - 1.0f) / cw;
else
mCoeff = cw * -0.5f;
mLpZ1 = 0.0f;
mLpZ2 = 0.0f;
mApZ1 = 0.0f;
}
template<typename Real>
void BandSplitterR<Real>::process(const al::span<const Real> input, Real *hpout, Real *lpout)
{
const Real ap_coeff{mCoeff};
const Real lp_coeff{mCoeff*0.5f + 0.5f};
Real lp_z1{mLpZ1};
Real lp_z2{mLpZ2};
Real ap_z1{mApZ1};
auto proc_sample = [ap_coeff,lp_coeff,&lp_z1,&lp_z2,&ap_z1,&lpout](const Real in) noexcept -> Real
{
/* Low-pass sample processing. */
Real d{(in - lp_z1) * lp_coeff};
Real lp_y{lp_z1 + d};
lp_z1 = lp_y + d;
d = (lp_y - lp_z2) * lp_coeff;
lp_y = lp_z2 + d;
lp_z2 = lp_y + d;
*(lpout++) = lp_y;
/* All-pass sample processing. */
Real ap_y{in*ap_coeff + ap_z1};
ap_z1 = in - ap_y*ap_coeff;
/* High-pass generated from removing low-passed output. */
return ap_y - lp_y;
};
std::transform(input.cbegin(), input.cend(), hpout, proc_sample);
mLpZ1 = lp_z1;
mLpZ2 = lp_z2;
mApZ1 = ap_z1;
}
template<typename Real>
void BandSplitterR<Real>::processHfScale(const al::span<Real> samples, const Real hfscale)
{
const Real ap_coeff{mCoeff};
const Real lp_coeff{mCoeff*0.5f + 0.5f};
Real lp_z1{mLpZ1};
Real lp_z2{mLpZ2};
Real ap_z1{mApZ1};
auto proc_sample = [hfscale,ap_coeff,lp_coeff,&lp_z1,&lp_z2,&ap_z1](const Real in) noexcept -> Real
{
/* Low-pass sample processing. */
Real d{(in - lp_z1) * lp_coeff};
Real lp_y{lp_z1 + d};
lp_z1 = lp_y + d;
d = (lp_y - lp_z2) * lp_coeff;
lp_y = lp_z2 + d;
lp_z2 = lp_y + d;
/* All-pass sample processing. */
Real ap_y{in*ap_coeff + ap_z1};
ap_z1 = in - ap_y*ap_coeff;
/* High-pass generated by removing the low-passed signal, which is then
* scaled and added back to the low-passed signal.
*/
return (ap_y-lp_y)*hfscale + lp_y;
};
std::transform(samples.begin(), samples.end(), samples.begin(), proc_sample);
mLpZ1 = lp_z1;
mLpZ2 = lp_z2;
mApZ1 = ap_z1;
}
template<typename Real>
void BandSplitterR<Real>::applyAllpass(const al::span<Real> samples) const
{
const Real coeff{mCoeff};
Real z1{0.0f};
auto proc_sample = [coeff,&z1](const Real in) noexcept -> Real
{
const Real out{in*coeff + z1};
z1 = in - out*coeff;
return out;
};
std::transform(samples.begin(), samples.end(), samples.begin(), proc_sample);
}
template class BandSplitterR<float>;
template class BandSplitterR<double>;

View file

@ -0,0 +1,36 @@
#ifndef CORE_FILTERS_SPLITTER_H
#define CORE_FILTERS_SPLITTER_H
#include <cstddef>
#include "alspan.h"
/* Band splitter. Splits a signal into two phase-matching frequency bands. */
template<typename Real>
class BandSplitterR {
Real mCoeff{0.0f};
Real mLpZ1{0.0f};
Real mLpZ2{0.0f};
Real mApZ1{0.0f};
public:
BandSplitterR() = default;
BandSplitterR(const BandSplitterR&) = default;
BandSplitterR(Real f0norm) { init(f0norm); }
void init(Real f0norm);
void clear() noexcept { mLpZ1 = mLpZ2 = mApZ1 = 0.0f; }
void process(const al::span<const Real> input, Real *hpout, Real *lpout);
void processHfScale(const al::span<Real> samples, const Real hfscale);
/* The all-pass portion of the band splitter. Applies the same phase shift
* without splitting the signal. Note that each use of this method is
* indepedent, it does not track history between calls.
*/
void applyAllpass(const al::span<Real> samples) const;
};
using BandSplitter = BandSplitterR<float>;
#endif /* CORE_FILTERS_SPLITTER_H */

View file

@ -0,0 +1,79 @@
#include "config.h"
#include "fmt_traits.h"
namespace al {
const int16_t muLawDecompressionTable[256] = {
-32124,-31100,-30076,-29052,-28028,-27004,-25980,-24956,
-23932,-22908,-21884,-20860,-19836,-18812,-17788,-16764,
-15996,-15484,-14972,-14460,-13948,-13436,-12924,-12412,
-11900,-11388,-10876,-10364, -9852, -9340, -8828, -8316,
-7932, -7676, -7420, -7164, -6908, -6652, -6396, -6140,
-5884, -5628, -5372, -5116, -4860, -4604, -4348, -4092,
-3900, -3772, -3644, -3516, -3388, -3260, -3132, -3004,
-2876, -2748, -2620, -2492, -2364, -2236, -2108, -1980,
-1884, -1820, -1756, -1692, -1628, -1564, -1500, -1436,
-1372, -1308, -1244, -1180, -1116, -1052, -988, -924,
-876, -844, -812, -780, -748, -716, -684, -652,
-620, -588, -556, -524, -492, -460, -428, -396,
-372, -356, -340, -324, -308, -292, -276, -260,
-244, -228, -212, -196, -180, -164, -148, -132,
-120, -112, -104, -96, -88, -80, -72, -64,
-56, -48, -40, -32, -24, -16, -8, 0,
32124, 31100, 30076, 29052, 28028, 27004, 25980, 24956,
23932, 22908, 21884, 20860, 19836, 18812, 17788, 16764,
15996, 15484, 14972, 14460, 13948, 13436, 12924, 12412,
11900, 11388, 10876, 10364, 9852, 9340, 8828, 8316,
7932, 7676, 7420, 7164, 6908, 6652, 6396, 6140,
5884, 5628, 5372, 5116, 4860, 4604, 4348, 4092,
3900, 3772, 3644, 3516, 3388, 3260, 3132, 3004,
2876, 2748, 2620, 2492, 2364, 2236, 2108, 1980,
1884, 1820, 1756, 1692, 1628, 1564, 1500, 1436,
1372, 1308, 1244, 1180, 1116, 1052, 988, 924,
876, 844, 812, 780, 748, 716, 684, 652,
620, 588, 556, 524, 492, 460, 428, 396,
372, 356, 340, 324, 308, 292, 276, 260,
244, 228, 212, 196, 180, 164, 148, 132,
120, 112, 104, 96, 88, 80, 72, 64,
56, 48, 40, 32, 24, 16, 8, 0
};
const int16_t aLawDecompressionTable[256] = {
-5504, -5248, -6016, -5760, -4480, -4224, -4992, -4736,
-7552, -7296, -8064, -7808, -6528, -6272, -7040, -6784,
-2752, -2624, -3008, -2880, -2240, -2112, -2496, -2368,
-3776, -3648, -4032, -3904, -3264, -3136, -3520, -3392,
-22016,-20992,-24064,-23040,-17920,-16896,-19968,-18944,
-30208,-29184,-32256,-31232,-26112,-25088,-28160,-27136,
-11008,-10496,-12032,-11520, -8960, -8448, -9984, -9472,
-15104,-14592,-16128,-15616,-13056,-12544,-14080,-13568,
-344, -328, -376, -360, -280, -264, -312, -296,
-472, -456, -504, -488, -408, -392, -440, -424,
-88, -72, -120, -104, -24, -8, -56, -40,
-216, -200, -248, -232, -152, -136, -184, -168,
-1376, -1312, -1504, -1440, -1120, -1056, -1248, -1184,
-1888, -1824, -2016, -1952, -1632, -1568, -1760, -1696,
-688, -656, -752, -720, -560, -528, -624, -592,
-944, -912, -1008, -976, -816, -784, -880, -848,
5504, 5248, 6016, 5760, 4480, 4224, 4992, 4736,
7552, 7296, 8064, 7808, 6528, 6272, 7040, 6784,
2752, 2624, 3008, 2880, 2240, 2112, 2496, 2368,
3776, 3648, 4032, 3904, 3264, 3136, 3520, 3392,
22016, 20992, 24064, 23040, 17920, 16896, 19968, 18944,
30208, 29184, 32256, 31232, 26112, 25088, 28160, 27136,
11008, 10496, 12032, 11520, 8960, 8448, 9984, 9472,
15104, 14592, 16128, 15616, 13056, 12544, 14080, 13568,
344, 328, 376, 360, 280, 264, 312, 296,
472, 456, 504, 488, 408, 392, 440, 424,
88, 72, 120, 104, 24, 8, 56, 40,
216, 200, 248, 232, 152, 136, 184, 168,
1376, 1312, 1504, 1440, 1120, 1056, 1248, 1184,
1888, 1824, 2016, 1952, 1632, 1568, 1760, 1696,
688, 656, 752, 720, 560, 528, 624, 592,
944, 912, 1008, 976, 816, 784, 880, 848
};
} // namespace al

View file

@ -0,0 +1,81 @@
#ifndef CORE_FMT_TRAITS_H
#define CORE_FMT_TRAITS_H
#include <stddef.h>
#include <stdint.h>
#include "albyte.h"
#include "buffer_storage.h"
namespace al {
extern const int16_t muLawDecompressionTable[256];
extern const int16_t aLawDecompressionTable[256];
template<FmtType T>
struct FmtTypeTraits { };
template<>
struct FmtTypeTraits<FmtUByte> {
using Type = uint8_t;
template<typename OutT>
static constexpr inline OutT to(const Type val) noexcept
{ return val*OutT{1.0/128.0} - OutT{1.0}; }
};
template<>
struct FmtTypeTraits<FmtShort> {
using Type = int16_t;
template<typename OutT>
static constexpr inline OutT to(const Type val) noexcept { return val*OutT{1.0/32768.0}; }
};
template<>
struct FmtTypeTraits<FmtFloat> {
using Type = float;
template<typename OutT>
static constexpr inline OutT to(const Type val) noexcept { return val; }
};
template<>
struct FmtTypeTraits<FmtDouble> {
using Type = double;
template<typename OutT>
static constexpr inline OutT to(const Type val) noexcept { return static_cast<OutT>(val); }
};
template<>
struct FmtTypeTraits<FmtMulaw> {
using Type = uint8_t;
template<typename OutT>
static constexpr inline OutT to(const Type val) noexcept
{ return muLawDecompressionTable[val] * OutT{1.0/32768.0}; }
};
template<>
struct FmtTypeTraits<FmtAlaw> {
using Type = uint8_t;
template<typename OutT>
static constexpr inline OutT to(const Type val) noexcept
{ return aLawDecompressionTable[val] * OutT{1.0/32768.0}; }
};
template<FmtType SrcType, typename DstT>
inline void LoadSampleArray(DstT *RESTRICT dst, const al::byte *src, const size_t srcstep,
const size_t samples) noexcept
{
using TypeTraits = FmtTypeTraits<SrcType>;
using SampleType = typename TypeTraits::Type;
const SampleType *RESTRICT ssrc{reinterpret_cast<const SampleType*>(src)};
for(size_t i{0u};i < samples;i++)
dst[i] = TypeTraits::template to<DstT>(ssrc[i*srcstep]);
}
} // namespace al
#endif /* CORE_FMT_TRAITS_H */

View file

@ -0,0 +1,56 @@
#include "config.h"
#include "fpu_ctrl.h"
#ifdef HAVE_INTRIN_H
#include <intrin.h>
#endif
#ifdef HAVE_SSE_INTRINSICS
#include <xmmintrin.h>
#endif
#include "cpu_caps.h"
void FPUCtl::enter() noexcept
{
if(this->in_mode) return;
#if defined(HAVE_SSE_INTRINSICS)
this->sse_state = _mm_getcsr();
unsigned int sseState{this->sse_state};
sseState |= 0x8000; /* set flush-to-zero */
sseState |= 0x0040; /* set denormals-are-zero */
_mm_setcsr(sseState);
#elif defined(__GNUC__) && defined(HAVE_SSE)
if((CPUCapFlags&CPU_CAP_SSE))
{
__asm__ __volatile__("stmxcsr %0" : "=m" (*&this->sse_state));
unsigned int sseState{this->sse_state};
sseState |= 0x8000; /* set flush-to-zero */
if((CPUCapFlags&CPU_CAP_SSE2))
sseState |= 0x0040; /* set denormals-are-zero */
__asm__ __volatile__("ldmxcsr %0" : : "m" (*&sseState));
}
#endif
this->in_mode = true;
}
void FPUCtl::leave() noexcept
{
if(!this->in_mode) return;
#if defined(HAVE_SSE_INTRINSICS)
_mm_setcsr(this->sse_state);
#elif defined(__GNUC__) && defined(HAVE_SSE)
if((CPUCapFlags&CPU_CAP_SSE))
__asm__ __volatile__("ldmxcsr %0" : : "m" (*&this->sse_state));
#endif
this->in_mode = false;
}

View file

@ -0,0 +1,21 @@
#ifndef CORE_FPU_CTRL_H
#define CORE_FPU_CTRL_H
class FPUCtl {
#if defined(HAVE_SSE_INTRINSICS) || (defined(__GNUC__) && defined(HAVE_SSE))
unsigned int sse_state{};
#endif
bool in_mode{};
public:
FPUCtl() noexcept { enter(); in_mode = true; };
~FPUCtl() { if(in_mode) leave(); }
FPUCtl(const FPUCtl&) = delete;
FPUCtl& operator=(const FPUCtl&) = delete;
void enter() noexcept;
void leave() noexcept;
};
#endif /* CORE_FPU_CTRL_H */

View file

@ -0,0 +1,95 @@
#include "config.h"
#include "logging.h"
#include <cstdarg>
#include <cstdio>
#include <string>
#include "strutils.h"
#include "vector.h"
#ifdef _WIN32
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
void al_print(LogLevel level, FILE *logfile, const char *fmt, ...)
{
al::vector<char> dynmsg;
char stcmsg[256];
char *str{stcmsg};
std::va_list args, args2;
va_start(args, fmt);
va_copy(args2, args);
int msglen{std::vsnprintf(str, sizeof(stcmsg), fmt, args)};
if UNLIKELY(msglen >= 0 && static_cast<size_t>(msglen) >= sizeof(stcmsg))
{
dynmsg.resize(static_cast<size_t>(msglen) + 1u);
str = dynmsg.data();
msglen = std::vsnprintf(str, dynmsg.size(), fmt, args2);
}
va_end(args2);
va_end(args);
std::wstring wstr{utf8_to_wstr(str)};
if(gLogLevel >= level)
{
fputws(wstr.c_str(), logfile);
fflush(logfile);
}
OutputDebugStringW(wstr.c_str());
}
#else
#ifdef __ANDROID__
#include <android/log.h>
#endif
void al_print(LogLevel level, FILE *logfile, const char *fmt, ...)
{
al::vector<char> dynmsg;
char stcmsg[256];
char *str{stcmsg};
std::va_list args, args2;
va_start(args, fmt);
va_copy(args2, args);
int msglen{std::vsnprintf(str, sizeof(stcmsg), fmt, args)};
if UNLIKELY(msglen >= 0 && static_cast<size_t>(msglen) >= sizeof(stcmsg))
{
dynmsg.resize(static_cast<size_t>(msglen) + 1u);
str = dynmsg.data();
msglen = std::vsnprintf(str, dynmsg.size(), fmt, args2);
}
va_end(args2);
va_end(args);
if(gLogLevel >= level)
{
std::fputs(str, logfile);
std::fflush(logfile);
}
#ifdef __ANDROID__
auto android_severity = [](LogLevel l) noexcept
{
switch(l)
{
case LogLevel::Trace: return ANDROID_LOG_DEBUG;
case LogLevel::Warning: return ANDROID_LOG_WARN;
case LogLevel::Error: return ANDROID_LOG_ERROR;
/* Should not happen. */
case LogLevel::Disable:
break;
}
return ANDROID_LOG_ERROR;
};
__android_log_print(android_severity(level), "openal", "%s", str);
#endif
}
#endif

View file

@ -0,0 +1,47 @@
#ifndef CORE_LOGGING_H
#define CORE_LOGGING_H
#include <stdio.h>
#include "opthelpers.h"
enum class LogLevel {
Disable,
Error,
Warning,
Trace
};
extern LogLevel gLogLevel;
extern FILE *gLogFile;
#if !defined(_WIN32) && !defined(__ANDROID__)
#define TRACE(...) do { \
if UNLIKELY(gLogLevel >= LogLevel::Trace) \
fprintf(gLogFile, "[ALSOFT] (II) " __VA_ARGS__); \
} while(0)
#define WARN(...) do { \
if UNLIKELY(gLogLevel >= LogLevel::Warning) \
fprintf(gLogFile, "[ALSOFT] (WW) " __VA_ARGS__); \
} while(0)
#define ERR(...) do { \
if UNLIKELY(gLogLevel >= LogLevel::Error) \
fprintf(gLogFile, "[ALSOFT] (EE) " __VA_ARGS__); \
} while(0)
#else
[[gnu::format(printf,3,4)]] void al_print(LogLevel level, FILE *logfile, const char *fmt, ...);
#define TRACE(...) al_print(LogLevel::Trace, gLogFile, "[ALSOFT] (II) " __VA_ARGS__)
#define WARN(...) al_print(LogLevel::Warning, gLogFile, "[ALSOFT] (WW) " __VA_ARGS__)
#define ERR(...) al_print(LogLevel::Error, gLogFile, "[ALSOFT] (EE) " __VA_ARGS__)
#endif
#endif /* CORE_LOGGING_H */

View file

@ -0,0 +1,441 @@
#include "config.h"
#include "mastering.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <limits>
#include <new>
#include "almalloc.h"
#include "alnumeric.h"
#include "alspan.h"
#include "opthelpers.h"
/* These structures assume BufferLineSize is a power of 2. */
static_assert((BufferLineSize & (BufferLineSize-1)) == 0, "BufferLineSize is not a power of 2");
struct SlidingHold {
alignas(16) float mValues[BufferLineSize];
uint mExpiries[BufferLineSize];
uint mLowerIndex;
uint mUpperIndex;
uint mLength;
};
namespace {
using namespace std::placeholders;
/* This sliding hold follows the input level with an instant attack and a
* fixed duration hold before an instant release to the next highest level.
* It is a sliding window maximum (descending maxima) implementation based on
* Richard Harter's ascending minima algorithm available at:
*
* http://www.richardhartersworld.com/cri/2001/slidingmin.html
*/
float UpdateSlidingHold(SlidingHold *Hold, const uint i, const float in)
{
static constexpr uint mask{BufferLineSize - 1};
const uint length{Hold->mLength};
float (&values)[BufferLineSize] = Hold->mValues;
uint (&expiries)[BufferLineSize] = Hold->mExpiries;
uint lowerIndex{Hold->mLowerIndex};
uint upperIndex{Hold->mUpperIndex};
if(i >= expiries[upperIndex])
upperIndex = (upperIndex + 1) & mask;
if(in >= values[upperIndex])
{
values[upperIndex] = in;
expiries[upperIndex] = i + length;
lowerIndex = upperIndex;
}
else
{
do {
do {
if(!(in >= values[lowerIndex]))
goto found_place;
} while(lowerIndex--);
lowerIndex = mask;
} while(1);
found_place:
lowerIndex = (lowerIndex + 1) & mask;
values[lowerIndex] = in;
expiries[lowerIndex] = i + length;
}
Hold->mLowerIndex = lowerIndex;
Hold->mUpperIndex = upperIndex;
return values[upperIndex];
}
void ShiftSlidingHold(SlidingHold *Hold, const uint n)
{
auto exp_begin = std::begin(Hold->mExpiries) + Hold->mUpperIndex;
auto exp_last = std::begin(Hold->mExpiries) + Hold->mLowerIndex;
if(exp_last-exp_begin < 0)
{
std::transform(exp_begin, std::end(Hold->mExpiries), exp_begin,
std::bind(std::minus<>{}, _1, n));
exp_begin = std::begin(Hold->mExpiries);
}
std::transform(exp_begin, exp_last+1, exp_begin, std::bind(std::minus<>{}, _1, n));
}
/* Multichannel compression is linked via the absolute maximum of all
* channels.
*/
void LinkChannels(Compressor *Comp, const uint SamplesToDo, const FloatBufferLine *OutBuffer)
{
const size_t numChans{Comp->mNumChans};
ASSUME(SamplesToDo > 0);
ASSUME(numChans > 0);
auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
std::fill(side_begin, side_begin+SamplesToDo, 0.0f);
auto fill_max = [SamplesToDo,side_begin](const FloatBufferLine &input) -> void
{
const float *RESTRICT buffer{al::assume_aligned<16>(input.data())};
auto max_abs = std::bind(maxf, _1, std::bind(static_cast<float(&)(float)>(std::fabs), _2));
std::transform(side_begin, side_begin+SamplesToDo, buffer, side_begin, max_abs);
};
std::for_each(OutBuffer, OutBuffer+numChans, fill_max);
}
/* This calculates the squared crest factor of the control signal for the
* basic automation of the attack/release times. As suggested by the paper,
* it uses an instantaneous squared peak detector and a squared RMS detector
* both with 200ms release times.
*/
static void CrestDetector(Compressor *Comp, const uint SamplesToDo)
{
const float a_crest{Comp->mCrestCoeff};
float y2_peak{Comp->mLastPeakSq};
float y2_rms{Comp->mLastRmsSq};
ASSUME(SamplesToDo > 0);
auto calc_crest = [&y2_rms,&y2_peak,a_crest](const float x_abs) noexcept -> float
{
const float x2{clampf(x_abs * x_abs, 0.000001f, 1000000.0f)};
y2_peak = maxf(x2, lerp(x2, y2_peak, a_crest));
y2_rms = lerp(x2, y2_rms, a_crest);
return y2_peak / y2_rms;
};
auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
std::transform(side_begin, side_begin+SamplesToDo, std::begin(Comp->mCrestFactor), calc_crest);
Comp->mLastPeakSq = y2_peak;
Comp->mLastRmsSq = y2_rms;
}
/* The side-chain starts with a simple peak detector (based on the absolute
* value of the incoming signal) and performs most of its operations in the
* log domain.
*/
void PeakDetector(Compressor *Comp, const uint SamplesToDo)
{
ASSUME(SamplesToDo > 0);
/* Clamp the minimum amplitude to near-zero and convert to logarithm. */
auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
std::transform(side_begin, side_begin+SamplesToDo, side_begin,
[](const float s) -> float { return std::log(maxf(0.000001f, s)); });
}
/* An optional hold can be used to extend the peak detector so it can more
* solidly detect fast transients. This is best used when operating as a
* limiter.
*/
void PeakHoldDetector(Compressor *Comp, const uint SamplesToDo)
{
ASSUME(SamplesToDo > 0);
SlidingHold *hold{Comp->mHold};
uint i{0};
auto detect_peak = [&i,hold](const float x_abs) -> float
{
const float x_G{std::log(maxf(0.000001f, x_abs))};
return UpdateSlidingHold(hold, i++, x_G);
};
auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
std::transform(side_begin, side_begin+SamplesToDo, side_begin, detect_peak);
ShiftSlidingHold(hold, SamplesToDo);
}
/* This is the heart of the feed-forward compressor. It operates in the log
* domain (to better match human hearing) and can apply some basic automation
* to knee width, attack/release times, make-up/post gain, and clipping
* reduction.
*/
void GainCompressor(Compressor *Comp, const uint SamplesToDo)
{
const bool autoKnee{Comp->mAuto.Knee};
const bool autoAttack{Comp->mAuto.Attack};
const bool autoRelease{Comp->mAuto.Release};
const bool autoPostGain{Comp->mAuto.PostGain};
const bool autoDeclip{Comp->mAuto.Declip};
const uint lookAhead{Comp->mLookAhead};
const float threshold{Comp->mThreshold};
const float slope{Comp->mSlope};
const float attack{Comp->mAttack};
const float release{Comp->mRelease};
const float c_est{Comp->mGainEstimate};
const float a_adp{Comp->mAdaptCoeff};
const float *crestFactor{Comp->mCrestFactor};
float postGain{Comp->mPostGain};
float knee{Comp->mKnee};
float t_att{attack};
float t_rel{release - attack};
float a_att{std::exp(-1.0f / t_att)};
float a_rel{std::exp(-1.0f / t_rel)};
float y_1{Comp->mLastRelease};
float y_L{Comp->mLastAttack};
float c_dev{Comp->mLastGainDev};
ASSUME(SamplesToDo > 0);
for(float &sideChain : al::span<float>{Comp->mSideChain, SamplesToDo})
{
if(autoKnee)
knee = maxf(0.0f, 2.5f * (c_dev + c_est));
const float knee_h{0.5f * knee};
/* This is the gain computer. It applies a static compression curve
* to the control signal.
*/
const float x_over{std::addressof(sideChain)[lookAhead] - threshold};
const float y_G{
(x_over <= -knee_h) ? 0.0f :
(std::fabs(x_over) < knee_h) ? (x_over + knee_h) * (x_over + knee_h) / (2.0f * knee) :
x_over};
const float y2_crest{*(crestFactor++)};
if(autoAttack)
{
t_att = 2.0f*attack/y2_crest;
a_att = std::exp(-1.0f / t_att);
}
if(autoRelease)
{
t_rel = 2.0f*release/y2_crest - t_att;
a_rel = std::exp(-1.0f / t_rel);
}
/* Gain smoothing (ballistics) is done via a smooth decoupled peak
* detector. The attack time is subtracted from the release time
* above to compensate for the chained operating mode.
*/
const float x_L{-slope * y_G};
y_1 = maxf(x_L, lerp(x_L, y_1, a_rel));
y_L = lerp(y_1, y_L, a_att);
/* Knee width and make-up gain automation make use of a smoothed
* measurement of deviation between the control signal and estimate.
* The estimate is also used to bias the measurement to hot-start its
* average.
*/
c_dev = lerp(-(y_L+c_est), c_dev, a_adp);
if(autoPostGain)
{
/* Clipping reduction is only viable when make-up gain is being
* automated. It modifies the deviation to further attenuate the
* control signal when clipping is detected. The adaptation time
* is sufficiently long enough to suppress further clipping at the
* same output level.
*/
if(autoDeclip)
c_dev = maxf(c_dev, sideChain - y_L - threshold - c_est);
postGain = -(c_dev + c_est);
}
sideChain = std::exp(postGain - y_L);
}
Comp->mLastRelease = y_1;
Comp->mLastAttack = y_L;
Comp->mLastGainDev = c_dev;
}
/* Combined with the hold time, a look-ahead delay can improve handling of
* fast transients by allowing the envelope time to converge prior to
* reaching the offending impulse. This is best used when operating as a
* limiter.
*/
void SignalDelay(Compressor *Comp, const uint SamplesToDo, FloatBufferLine *OutBuffer)
{
const size_t numChans{Comp->mNumChans};
const uint lookAhead{Comp->mLookAhead};
ASSUME(SamplesToDo > 0);
ASSUME(numChans > 0);
ASSUME(lookAhead > 0);
for(size_t c{0};c < numChans;c++)
{
float *inout{al::assume_aligned<16>(OutBuffer[c].data())};
float *delaybuf{al::assume_aligned<16>(Comp->mDelay[c].data())};
auto inout_end = inout + SamplesToDo;
if LIKELY(SamplesToDo >= lookAhead)
{
auto delay_end = std::rotate(inout, inout_end - lookAhead, inout_end);
std::swap_ranges(inout, delay_end, delaybuf);
}
else
{
auto delay_start = std::swap_ranges(inout, inout_end, delaybuf);
std::rotate(delaybuf, delay_start, delaybuf + lookAhead);
}
}
}
} // namespace
std::unique_ptr<Compressor> Compressor::Create(const size_t NumChans, const float SampleRate,
const bool AutoKnee, const bool AutoAttack, const bool AutoRelease, const bool AutoPostGain,
const bool AutoDeclip, const float LookAheadTime, const float HoldTime, const float PreGainDb,
const float PostGainDb, const float ThresholdDb, const float Ratio, const float KneeDb,
const float AttackTime, const float ReleaseTime)
{
const auto lookAhead = static_cast<uint>(
clampf(std::round(LookAheadTime*SampleRate), 0.0f, BufferLineSize-1));
const auto hold = static_cast<uint>(
clampf(std::round(HoldTime*SampleRate), 0.0f, BufferLineSize-1));
size_t size{sizeof(Compressor)};
if(lookAhead > 0)
{
size += sizeof(*Compressor::mDelay) * NumChans;
/* The sliding hold implementation doesn't handle a length of 1. A 1-
* sample hold is useless anyway, it would only ever give back what was
* just given to it.
*/
if(hold > 1)
size += sizeof(*Compressor::mHold);
}
auto Comp = std::unique_ptr<Compressor>{new (al_calloc(16, size)) Compressor{}};
Comp->mNumChans = NumChans;
Comp->mAuto.Knee = AutoKnee;
Comp->mAuto.Attack = AutoAttack;
Comp->mAuto.Release = AutoRelease;
Comp->mAuto.PostGain = AutoPostGain;
Comp->mAuto.Declip = AutoPostGain && AutoDeclip;
Comp->mLookAhead = lookAhead;
Comp->mPreGain = std::pow(10.0f, PreGainDb / 20.0f);
Comp->mPostGain = PostGainDb * std::log(10.0f) / 20.0f;
Comp->mThreshold = ThresholdDb * std::log(10.0f) / 20.0f;
Comp->mSlope = 1.0f / maxf(1.0f, Ratio) - 1.0f;
Comp->mKnee = maxf(0.0f, KneeDb * std::log(10.0f) / 20.0f);
Comp->mAttack = maxf(1.0f, AttackTime * SampleRate);
Comp->mRelease = maxf(1.0f, ReleaseTime * SampleRate);
/* Knee width automation actually treats the compressor as a limiter. By
* varying the knee width, it can effectively be seen as applying
* compression over a wide range of ratios.
*/
if(AutoKnee)
Comp->mSlope = -1.0f;
if(lookAhead > 0)
{
if(hold > 1)
{
Comp->mHold = ::new (static_cast<void*>(Comp.get() + 1)) SlidingHold{};
Comp->mHold->mValues[0] = -std::numeric_limits<float>::infinity();
Comp->mHold->mExpiries[0] = hold;
Comp->mHold->mLength = hold;
Comp->mDelay = ::new(static_cast<void*>(Comp->mHold + 1)) FloatBufferLine[NumChans];
}
else
{
Comp->mDelay = ::new(static_cast<void*>(Comp.get() + 1)) FloatBufferLine[NumChans];
}
std::fill_n(Comp->mDelay, NumChans, FloatBufferLine{});
}
Comp->mCrestCoeff = std::exp(-1.0f / (0.200f * SampleRate)); // 200ms
Comp->mGainEstimate = Comp->mThreshold * -0.5f * Comp->mSlope;
Comp->mAdaptCoeff = std::exp(-1.0f / (2.0f * SampleRate)); // 2s
return Comp;
}
Compressor::~Compressor()
{
if(mHold)
al::destroy_at(mHold);
mHold = nullptr;
if(mDelay)
al::destroy_n(mDelay, mNumChans);
mDelay = nullptr;
}
void Compressor::process(const uint SamplesToDo, FloatBufferLine *OutBuffer)
{
const size_t numChans{mNumChans};
ASSUME(SamplesToDo > 0);
ASSUME(numChans > 0);
const float preGain{mPreGain};
if(preGain != 1.0f)
{
auto apply_gain = [SamplesToDo,preGain](FloatBufferLine &input) noexcept -> void
{
float *buffer{al::assume_aligned<16>(input.data())};
std::transform(buffer, buffer+SamplesToDo, buffer,
std::bind(std::multiplies<float>{}, _1, preGain));
};
std::for_each(OutBuffer, OutBuffer+numChans, apply_gain);
}
LinkChannels(this, SamplesToDo, OutBuffer);
if(mAuto.Attack || mAuto.Release)
CrestDetector(this, SamplesToDo);
if(mHold)
PeakHoldDetector(this, SamplesToDo);
else
PeakDetector(this, SamplesToDo);
GainCompressor(this, SamplesToDo);
if(mDelay)
SignalDelay(this, SamplesToDo, OutBuffer);
const float (&sideChain)[BufferLineSize*2] = mSideChain;
auto apply_comp = [SamplesToDo,&sideChain](FloatBufferLine &input) noexcept -> void
{
float *buffer{al::assume_aligned<16>(input.data())};
const float *gains{al::assume_aligned<16>(&sideChain[0])};
std::transform(gains, gains+SamplesToDo, buffer, buffer,
std::bind(std::multiplies<float>{}, _1, _2));
};
std::for_each(OutBuffer, OutBuffer+numChans, apply_comp);
auto side_begin = std::begin(mSideChain) + SamplesToDo;
std::copy(side_begin, side_begin+mLookAhead, std::begin(mSideChain));
}

View file

@ -0,0 +1,104 @@
#ifndef CORE_MASTERING_H
#define CORE_MASTERING_H
#include <memory>
#include "almalloc.h"
#include "bufferline.h"
struct SlidingHold;
using uint = unsigned int;
/* General topology and basic automation was based on the following paper:
*
* D. Giannoulis, M. Massberg and J. D. Reiss,
* "Parameter Automation in a Dynamic Range Compressor,"
* Journal of the Audio Engineering Society, v61 (10), Oct. 2013
*
* Available (along with supplemental reading) at:
*
* http://c4dm.eecs.qmul.ac.uk/audioengineering/compressors/
*/
struct Compressor {
size_t mNumChans{0u};
struct {
bool Knee : 1;
bool Attack : 1;
bool Release : 1;
bool PostGain : 1;
bool Declip : 1;
} mAuto{};
uint mLookAhead{0};
float mPreGain{0.0f};
float mPostGain{0.0f};
float mThreshold{0.0f};
float mSlope{0.0f};
float mKnee{0.0f};
float mAttack{0.0f};
float mRelease{0.0f};
alignas(16) float mSideChain[2*BufferLineSize]{};
alignas(16) float mCrestFactor[BufferLineSize]{};
SlidingHold *mHold{nullptr};
FloatBufferLine *mDelay{nullptr};
float mCrestCoeff{0.0f};
float mGainEstimate{0.0f};
float mAdaptCoeff{0.0f};
float mLastPeakSq{0.0f};
float mLastRmsSq{0.0f};
float mLastRelease{0.0f};
float mLastAttack{0.0f};
float mLastGainDev{0.0f};
~Compressor();
void process(const uint SamplesToDo, FloatBufferLine *OutBuffer);
int getLookAhead() const noexcept { return static_cast<int>(mLookAhead); }
DEF_PLACE_NEWDEL()
/**
* The compressor is initialized with the following settings:
*
* \param NumChans Number of channels to process.
* \param SampleRate Sample rate to process.
* \param AutoKnee Whether to automate the knee width parameter.
* \param AutoAttack Whether to automate the attack time parameter.
* \param AutoRelease Whether to automate the release time parameter.
* \param AutoPostGain Whether to automate the make-up (post) gain
* parameter.
* \param AutoDeclip Whether to automate clipping reduction. Ignored
* when not automating make-up gain.
* \param LookAheadTime Look-ahead time (in seconds).
* \param HoldTime Peak hold-time (in seconds).
* \param PreGainDb Gain applied before detection (in dB).
* \param PostGainDb Make-up gain applied after compression (in dB).
* \param ThresholdDb Triggering threshold (in dB).
* \param Ratio Compression ratio (x:1). Set to INFINIFTY for true
* limiting. Ignored when automating knee width.
* \param KneeDb Knee width (in dB). Ignored when automating knee
* width.
* \param AttackTime Attack time (in seconds). Acts as a maximum when
* automating attack time.
* \param ReleaseTime Release time (in seconds). Acts as a maximum when
* automating release time.
*/
static std::unique_ptr<Compressor> Create(const size_t NumChans, const float SampleRate,
const bool AutoKnee, const bool AutoAttack, const bool AutoRelease,
const bool AutoPostGain, const bool AutoDeclip, const float LookAheadTime,
const float HoldTime, const float PreGainDb, const float PostGainDb,
const float ThresholdDb, const float Ratio, const float KneeDb, const float AttackTime,
const float ReleaseTime);
};
#endif /* CORE_MASTERING_H */

View file

@ -0,0 +1,101 @@
#ifndef CORE_MIXER_DEFS_H
#define CORE_MIXER_DEFS_H
#include <array>
#include <stdlib.h>
#include "alspan.h"
#include "core/bufferline.h"
struct HrtfChannelState;
struct HrtfFilter;
struct MixHrtfFilter;
using uint = unsigned int;
using float2 = std::array<float,2>;
constexpr int MixerFracBits{12};
constexpr int MixerFracOne{1 << MixerFracBits};
constexpr int MixerFracMask{MixerFracOne - 1};
/* Maximum number of samples to pad on the ends of a buffer for resampling.
* Note that the padding is symmetric (half at the beginning and half at the
* end)!
*/
constexpr int MaxResamplerPadding{48};
constexpr float GainSilenceThreshold{0.00001f}; /* -100dB */
enum class Resampler {
Point,
Linear,
Cubic,
FastBSinc12,
BSinc12,
FastBSinc24,
BSinc24,
Max = BSinc24
};
/* Interpolator state. Kind of a misnomer since the interpolator itself is
* stateless. This just keeps it from having to recompute scale-related
* mappings for every sample.
*/
struct BsincState {
float sf; /* Scale interpolation factor. */
uint m; /* Coefficient count. */
uint l; /* Left coefficient offset. */
/* Filter coefficients, followed by the phase, scale, and scale-phase
* delta coefficients. Starting at phase index 0, each subsequent phase
* index follows contiguously.
*/
const float *filter;
};
union InterpState {
BsincState bsinc;
};
using ResamplerFunc = float*(*)(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst);
ResamplerFunc PrepareResampler(Resampler resampler, uint increment, InterpState *state);
template<typename TypeTag, typename InstTag>
float *Resample_(const InterpState *state, float *RESTRICT src, uint frac, uint increment,
const al::span<float> dst);
template<typename InstTag>
void Mix_(const al::span<const float> InSamples, const al::span<FloatBufferLine> OutBuffer,
float *CurrentGains, const float *TargetGains, const size_t Counter, const size_t OutPos);
template<typename InstTag>
void MixHrtf_(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const MixHrtfFilter *hrtfparams, const size_t BufferSize);
template<typename InstTag>
void MixHrtfBlend_(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const HrtfFilter *oldparams, const MixHrtfFilter *newparams, const size_t BufferSize);
template<typename InstTag>
void MixDirectHrtf_(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const al::span<const FloatBufferLine> InSamples, float2 *AccumSamples,
float *TempBuf, HrtfChannelState *ChanState, const size_t IrSize, const size_t BufferSize);
/* Vectorized resampler helpers */
template<size_t N>
inline void InitPosArrays(uint frac, uint increment, uint (&frac_arr)[N], uint (&pos_arr)[N])
{
pos_arr[0] = 0;
frac_arr[0] = frac;
for(size_t i{1};i < N;i++)
{
const uint frac_tmp{frac_arr[i-1] + increment};
pos_arr[i] = pos_arr[i-1] + (frac_tmp>>MixerFracBits);
frac_arr[i] = frac_tmp&MixerFracMask;
}
}
#endif /* CORE_MIXER_DEFS_H */

View file

@ -0,0 +1,159 @@
#ifndef CORE_MIXER_HRTFBASE_H
#define CORE_MIXER_HRTFBASE_H
#include <algorithm>
#include <cmath>
#include "almalloc.h"
#include "hrtfdefs.h"
#include "opthelpers.h"
using uint = unsigned int;
using ApplyCoeffsT = void(&)(float2 *RESTRICT Values, const size_t irSize,
const HrirArray &Coeffs, const float left, const float right);
template<ApplyCoeffsT ApplyCoeffs>
inline void MixHrtfBase(const float *InSamples, float2 *RESTRICT AccumSamples, const size_t IrSize,
const MixHrtfFilter *hrtfparams, const size_t BufferSize)
{
ASSUME(BufferSize > 0);
const HrirArray &Coeffs = *hrtfparams->Coeffs;
const float gainstep{hrtfparams->GainStep};
const float gain{hrtfparams->Gain};
size_t ldelay{HrtfHistoryLength - hrtfparams->Delay[0]};
size_t rdelay{HrtfHistoryLength - hrtfparams->Delay[1]};
float stepcount{0.0f};
for(size_t i{0u};i < BufferSize;++i)
{
const float g{gain + gainstep*stepcount};
const float left{InSamples[ldelay++] * g};
const float right{InSamples[rdelay++] * g};
ApplyCoeffs(AccumSamples+i, IrSize, Coeffs, left, right);
stepcount += 1.0f;
}
}
template<ApplyCoeffsT ApplyCoeffs>
inline void MixHrtfBlendBase(const float *InSamples, float2 *RESTRICT AccumSamples,
const size_t IrSize, const HrtfFilter *oldparams, const MixHrtfFilter *newparams,
const size_t BufferSize)
{
ASSUME(BufferSize > 0);
const auto &OldCoeffs = oldparams->Coeffs;
const float oldGainStep{oldparams->Gain / static_cast<float>(BufferSize)};
const auto &NewCoeffs = *newparams->Coeffs;
const float newGainStep{newparams->GainStep};
if LIKELY(oldparams->Gain > GainSilenceThreshold)
{
size_t ldelay{HrtfHistoryLength - oldparams->Delay[0]};
size_t rdelay{HrtfHistoryLength - oldparams->Delay[1]};
auto stepcount = static_cast<float>(BufferSize);
for(size_t i{0u};i < BufferSize;++i)
{
const float g{oldGainStep*stepcount};
const float left{InSamples[ldelay++] * g};
const float right{InSamples[rdelay++] * g};
ApplyCoeffs(AccumSamples+i, IrSize, OldCoeffs, left, right);
stepcount -= 1.0f;
}
}
if LIKELY(newGainStep*static_cast<float>(BufferSize) > GainSilenceThreshold)
{
size_t ldelay{HrtfHistoryLength+1 - newparams->Delay[0]};
size_t rdelay{HrtfHistoryLength+1 - newparams->Delay[1]};
float stepcount{1.0f};
for(size_t i{1u};i < BufferSize;++i)
{
const float g{newGainStep*stepcount};
const float left{InSamples[ldelay++] * g};
const float right{InSamples[rdelay++] * g};
ApplyCoeffs(AccumSamples+i, IrSize, NewCoeffs, left, right);
stepcount += 1.0f;
}
}
}
template<ApplyCoeffsT ApplyCoeffs>
inline void MixDirectHrtfBase(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const al::span<const FloatBufferLine> InSamples, float2 *RESTRICT AccumSamples,
float *TempBuf, HrtfChannelState *ChanState, const size_t IrSize, const size_t BufferSize)
{
ASSUME(BufferSize > 0);
/* Add the existing signal directly to the accumulation buffer, unfiltered,
* and with a delay to align with the input delay.
*/
for(size_t i{0};i < BufferSize;++i)
{
AccumSamples[HrtfDirectDelay+i][0] += LeftOut[i];
AccumSamples[HrtfDirectDelay+i][1] += RightOut[i];
}
for(const FloatBufferLine &input : InSamples)
{
/* For dual-band processing, the signal needs extra scaling applied to
* the high frequency response. The band-splitter alone creates a
* frequency-dependent phase shift, which is not ideal. To counteract
* it, combine it with a backwards phase shift.
*/
/* Load the input signal backwards, into a temp buffer with delay
* padding. The delay serves to reduce the error caused by the IIR
* filter's phase shift on a partial input.
*/
al::span<float> tempbuf{al::assume_aligned<16>(TempBuf), HrtfDirectDelay+BufferSize};
auto tmpiter = std::reverse_copy(input.begin(), input.begin()+BufferSize, tempbuf.begin());
std::copy(ChanState->mDelay.cbegin(), ChanState->mDelay.cend(), tmpiter);
/* Save the unfiltered newest input samples for next time. */
std::copy_n(tempbuf.begin(), ChanState->mDelay.size(), ChanState->mDelay.begin());
/* Apply the all-pass on the reversed signal and reverse the resulting
* sample array. This produces the forward response with a backwards
* phase shift (+n degrees becomes -n degrees).
*/
ChanState->mSplitter.applyAllpass(tempbuf);
tempbuf = tempbuf.subspan<HrtfDirectDelay>();
std::reverse(tempbuf.begin(), tempbuf.end());
/* Now apply the HF scale with the band-splitter. This applies the
* forward phase shift, which cancels out with the backwards phase
* shift to get the original phase on the scaled signal.
*/
ChanState->mSplitter.processHfScale(tempbuf, ChanState->mHfScale);
/* Now apply the HRIR coefficients to this channel. */
const auto &Coeffs = ChanState->mCoeffs;
for(size_t i{0u};i < BufferSize;++i)
{
const float insample{tempbuf[i]};
ApplyCoeffs(AccumSamples+i, IrSize, Coeffs, insample, insample);
}
++ChanState;
}
for(size_t i{0u};i < BufferSize;++i)
LeftOut[i] = AccumSamples[i][0];
for(size_t i{0u};i < BufferSize;++i)
RightOut[i] = AccumSamples[i][1];
/* Copy the new in-progress accumulation values to the front and clear the
* following samples for the next mix.
*/
auto accum_iter = std::copy_n(AccumSamples+BufferSize, HrirLength+HrtfDirectDelay,
AccumSamples);
std::fill_n(accum_iter, BufferSize, float2{});
}
#endif /* CORE_MIXER_HRTFBASE_H */

View file

@ -0,0 +1,53 @@
#ifndef CORE_MIXER_HRTFDEFS_H
#define CORE_MIXER_HRTFDEFS_H
#include <array>
#include "core/ambidefs.h"
#include "core/bufferline.h"
#include "core/filters/splitter.h"
using float2 = std::array<float,2>;
using ubyte = unsigned char;
using ubyte2 = std::array<ubyte,2>;
using ushort = unsigned short;
using uint = unsigned int;
using uint2 = std::array<uint,2>;
constexpr uint HrtfHistoryBits{6};
constexpr uint HrtfHistoryLength{1 << HrtfHistoryBits};
constexpr uint HrtfHistoryMask{HrtfHistoryLength - 1};
constexpr uint HrirBits{7};
constexpr uint HrirLength{1 << HrirBits};
constexpr uint HrirMask{HrirLength - 1};
constexpr uint MinIrLength{8};
constexpr uint HrtfDirectDelay{256};
using HrirArray = std::array<float2,HrirLength>;
struct MixHrtfFilter {
const HrirArray *Coeffs;
uint2 Delay;
float Gain;
float GainStep;
};
struct HrtfFilter {
alignas(16) HrirArray Coeffs;
uint2 Delay;
float Gain;
};
struct HrtfChannelState {
std::array<float,HrtfDirectDelay> mDelay{};
BandSplitter mSplitter;
float mHfScale{};
alignas(16) HrirArray mCoeffs{};
};
#endif /* CORE_MIXER_HRTFDEFS_H */

View file

@ -0,0 +1,198 @@
#include "config.h"
#include <cassert>
#include <cmath>
#include <limits>
#include "alnumeric.h"
#include "core/bsinc_tables.h"
#include "defs.h"
#include "hrtfbase.h"
struct CTag;
struct CopyTag;
struct PointTag;
struct LerpTag;
struct CubicTag;
struct BSincTag;
struct FastBSincTag;
namespace {
constexpr uint FracPhaseBitDiff{MixerFracBits - BSincPhaseBits};
constexpr uint FracPhaseDiffOne{1 << FracPhaseBitDiff};
inline float do_point(const InterpState&, const float *RESTRICT vals, const uint)
{ return vals[0]; }
inline float do_lerp(const InterpState&, const float *RESTRICT vals, const uint frac)
{ return lerp(vals[0], vals[1], static_cast<float>(frac)*(1.0f/MixerFracOne)); }
inline float do_cubic(const InterpState&, const float *RESTRICT vals, const uint frac)
{ return cubic(vals[0], vals[1], vals[2], vals[3], static_cast<float>(frac)*(1.0f/MixerFracOne)); }
inline float do_bsinc(const InterpState &istate, const float *RESTRICT vals, const uint frac)
{
const size_t m{istate.bsinc.m};
// Calculate the phase index and factor.
const uint pi{frac >> FracPhaseBitDiff};
const float pf{static_cast<float>(frac & (FracPhaseDiffOne-1)) * (1.0f/FracPhaseDiffOne)};
const float *fil{istate.bsinc.filter + m*pi*4};
const float *phd{fil + m};
const float *scd{phd + m};
const float *spd{scd + m};
// Apply the scale and phase interpolated filter.
float r{0.0f};
for(size_t j_f{0};j_f < m;j_f++)
r += (fil[j_f] + istate.bsinc.sf*scd[j_f] + pf*(phd[j_f] + istate.bsinc.sf*spd[j_f])) * vals[j_f];
return r;
}
inline float do_fastbsinc(const InterpState &istate, const float *RESTRICT vals, const uint frac)
{
const size_t m{istate.bsinc.m};
// Calculate the phase index and factor.
const uint pi{frac >> FracPhaseBitDiff};
const float pf{static_cast<float>(frac & (FracPhaseDiffOne-1)) * (1.0f/FracPhaseDiffOne)};
const float *fil{istate.bsinc.filter + m*pi*4};
const float *phd{fil + m};
// Apply the phase interpolated filter.
float r{0.0f};
for(size_t j_f{0};j_f < m;j_f++)
r += (fil[j_f] + pf*phd[j_f]) * vals[j_f];
return r;
}
using SamplerT = float(&)(const InterpState&, const float*RESTRICT, const uint);
template<SamplerT Sampler>
float *DoResample(const InterpState *state, float *RESTRICT src, uint frac, uint increment,
const al::span<float> dst)
{
const InterpState istate{*state};
for(float &out : dst)
{
out = Sampler(istate, src, frac);
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
}
return dst.data();
}
inline void ApplyCoeffs(float2 *RESTRICT Values, const size_t IrSize, const HrirArray &Coeffs,
const float left, const float right)
{
ASSUME(IrSize >= MinIrLength);
for(size_t c{0};c < IrSize;++c)
{
Values[c][0] += Coeffs[c][0] * left;
Values[c][1] += Coeffs[c][1] * right;
}
}
} // namespace
template<>
float *Resample_<CopyTag,CTag>(const InterpState*, float *RESTRICT src, uint, uint,
const al::span<float> dst)
{
#if defined(HAVE_SSE) || defined(HAVE_NEON)
/* Avoid copying the source data if it's aligned like the destination. */
if((reinterpret_cast<intptr_t>(src)&15) == (reinterpret_cast<intptr_t>(dst.data())&15))
return src;
#endif
std::copy_n(src, dst.size(), dst.begin());
return dst.data();
}
template<>
float *Resample_<PointTag,CTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{ return DoResample<do_point>(state, src, frac, increment, dst); }
template<>
float *Resample_<LerpTag,CTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{ return DoResample<do_lerp>(state, src, frac, increment, dst); }
template<>
float *Resample_<CubicTag,CTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{ return DoResample<do_cubic>(state, src-1, frac, increment, dst); }
template<>
float *Resample_<BSincTag,CTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{ return DoResample<do_bsinc>(state, src-state->bsinc.l, frac, increment, dst); }
template<>
float *Resample_<FastBSincTag,CTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{ return DoResample<do_fastbsinc>(state, src-state->bsinc.l, frac, increment, dst); }
template<>
void MixHrtf_<CTag>(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const MixHrtfFilter *hrtfparams, const size_t BufferSize)
{ MixHrtfBase<ApplyCoeffs>(InSamples, AccumSamples, IrSize, hrtfparams, BufferSize); }
template<>
void MixHrtfBlend_<CTag>(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const HrtfFilter *oldparams, const MixHrtfFilter *newparams, const size_t BufferSize)
{
MixHrtfBlendBase<ApplyCoeffs>(InSamples, AccumSamples, IrSize, oldparams, newparams,
BufferSize);
}
template<>
void MixDirectHrtf_<CTag>(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const al::span<const FloatBufferLine> InSamples, float2 *AccumSamples,
float *TempBuf, HrtfChannelState *ChanState, const size_t IrSize, const size_t BufferSize)
{
MixDirectHrtfBase<ApplyCoeffs>(LeftOut, RightOut, InSamples, AccumSamples, TempBuf, ChanState,
IrSize, BufferSize);
}
template<>
void Mix_<CTag>(const al::span<const float> InSamples, const al::span<FloatBufferLine> OutBuffer,
float *CurrentGains, const float *TargetGains, const size_t Counter, const size_t OutPos)
{
const float delta{(Counter > 0) ? 1.0f / static_cast<float>(Counter) : 0.0f};
const auto min_len = minz(Counter, InSamples.size());
for(FloatBufferLine &output : OutBuffer)
{
float *RESTRICT dst{al::assume_aligned<16>(output.data()+OutPos)};
float gain{*CurrentGains};
const float step{(*TargetGains-gain) * delta};
size_t pos{0};
if(!(std::abs(step) > std::numeric_limits<float>::epsilon()))
gain = *TargetGains;
else
{
float step_count{0.0f};
for(;pos != min_len;++pos)
{
dst[pos] += InSamples[pos] * (gain + step*step_count);
step_count += 1.0f;
}
if(pos == Counter)
gain = *TargetGains;
else
gain += step*step_count;
}
*CurrentGains = gain;
++CurrentGains;
++TargetGains;
if(!(std::abs(gain) > GainSilenceThreshold))
continue;
for(;pos != InSamples.size();++pos)
dst[pos] += InSamples[pos] * gain;
}
}

View file

@ -0,0 +1,305 @@
#include "config.h"
#include <arm_neon.h>
#include <cmath>
#include <limits>
#include "alnumeric.h"
#include "core/bsinc_defs.h"
#include "defs.h"
#include "hrtfbase.h"
struct NEONTag;
struct LerpTag;
struct BSincTag;
struct FastBSincTag;
#if defined(__GNUC__) && !defined(__clang__) && !defined(__ARM_NEON)
#pragma GCC target("fpu=neon")
#endif
namespace {
inline float32x4_t set_f4(float l0, float l1, float l2, float l3)
{
float32x4_t ret{vmovq_n_f32(l0)};
ret = vsetq_lane_f32(l1, ret, 1);
ret = vsetq_lane_f32(l2, ret, 2);
ret = vsetq_lane_f32(l3, ret, 3);
return ret;
}
constexpr uint FracPhaseBitDiff{MixerFracBits - BSincPhaseBits};
constexpr uint FracPhaseDiffOne{1 << FracPhaseBitDiff};
inline void ApplyCoeffs(float2 *RESTRICT Values, const size_t IrSize, const HrirArray &Coeffs,
const float left, const float right)
{
float32x4_t leftright4;
{
float32x2_t leftright2{vmov_n_f32(left)};
leftright2 = vset_lane_f32(right, leftright2, 1);
leftright4 = vcombine_f32(leftright2, leftright2);
}
ASSUME(IrSize >= MinIrLength);
for(size_t c{0};c < IrSize;c += 2)
{
float32x4_t vals = vld1q_f32(&Values[c][0]);
float32x4_t coefs = vld1q_f32(&Coeffs[c][0]);
vals = vmlaq_f32(vals, coefs, leftright4);
vst1q_f32(&Values[c][0], vals);
}
}
} // namespace
template<>
float *Resample_<LerpTag,NEONTag>(const InterpState*, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const int32x4_t increment4 = vdupq_n_s32(static_cast<int>(increment*4));
const float32x4_t fracOne4 = vdupq_n_f32(1.0f/MixerFracOne);
const int32x4_t fracMask4 = vdupq_n_s32(MixerFracMask);
alignas(16) uint pos_[4], frac_[4];
int32x4_t pos4, frac4;
InitPosArrays(frac, increment, frac_, pos_);
frac4 = vld1q_s32(reinterpret_cast<int*>(frac_));
pos4 = vld1q_s32(reinterpret_cast<int*>(pos_));
auto dst_iter = dst.begin();
for(size_t todo{dst.size()>>2};todo;--todo)
{
const int pos0{vgetq_lane_s32(pos4, 0)};
const int pos1{vgetq_lane_s32(pos4, 1)};
const int pos2{vgetq_lane_s32(pos4, 2)};
const int pos3{vgetq_lane_s32(pos4, 3)};
const float32x4_t val1{set_f4(src[pos0], src[pos1], src[pos2], src[pos3])};
const float32x4_t val2{set_f4(src[pos0+1], src[pos1+1], src[pos2+1], src[pos3+1])};
/* val1 + (val2-val1)*mu */
const float32x4_t r0{vsubq_f32(val2, val1)};
const float32x4_t mu{vmulq_f32(vcvtq_f32_s32(frac4), fracOne4)};
const float32x4_t out{vmlaq_f32(val1, mu, r0)};
vst1q_f32(dst_iter, out);
dst_iter += 4;
frac4 = vaddq_s32(frac4, increment4);
pos4 = vaddq_s32(pos4, vshrq_n_s32(frac4, MixerFracBits));
frac4 = vandq_s32(frac4, fracMask4);
}
if(size_t todo{dst.size()&3})
{
src += static_cast<uint>(vgetq_lane_s32(pos4, 0));
frac = static_cast<uint>(vgetq_lane_s32(frac4, 0));
do {
*(dst_iter++) = lerp(src[0], src[1], static_cast<float>(frac) * (1.0f/MixerFracOne));
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
} while(--todo);
}
return dst.data();
}
template<>
float *Resample_<BSincTag,NEONTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const float *const filter{state->bsinc.filter};
const float32x4_t sf4{vdupq_n_f32(state->bsinc.sf)};
const size_t m{state->bsinc.m};
src -= state->bsinc.l;
for(float &out_sample : dst)
{
// Calculate the phase index and factor.
const uint pi{frac >> FracPhaseBitDiff};
const float pf{static_cast<float>(frac & (FracPhaseDiffOne-1)) * (1.0f/FracPhaseDiffOne)};
// Apply the scale and phase interpolated filter.
float32x4_t r4{vdupq_n_f32(0.0f)};
{
const float32x4_t pf4{vdupq_n_f32(pf)};
const float *fil{filter + m*pi*4};
const float *phd{fil + m};
const float *scd{phd + m};
const float *spd{scd + m};
size_t td{m >> 2};
size_t j{0u};
do {
/* f = ((fil + sf*scd) + pf*(phd + sf*spd)) */
const float32x4_t f4 = vmlaq_f32(
vmlaq_f32(vld1q_f32(&fil[j]), sf4, vld1q_f32(&scd[j])),
pf4, vmlaq_f32(vld1q_f32(&phd[j]), sf4, vld1q_f32(&spd[j])));
/* r += f*src */
r4 = vmlaq_f32(r4, f4, vld1q_f32(&src[j]));
j += 4;
} while(--td);
}
r4 = vaddq_f32(r4, vrev64q_f32(r4));
out_sample = vget_lane_f32(vadd_f32(vget_low_f32(r4), vget_high_f32(r4)), 0);
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
}
return dst.data();
}
template<>
float *Resample_<FastBSincTag,NEONTag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const float *const filter{state->bsinc.filter};
const size_t m{state->bsinc.m};
src -= state->bsinc.l;
for(float &out_sample : dst)
{
// Calculate the phase index and factor.
const uint pi{frac >> FracPhaseBitDiff};
const float pf{static_cast<float>(frac & (FracPhaseDiffOne-1)) * (1.0f/FracPhaseDiffOne)};
// Apply the phase interpolated filter.
float32x4_t r4{vdupq_n_f32(0.0f)};
{
const float32x4_t pf4{vdupq_n_f32(pf)};
const float *fil{filter + m*pi*4};
const float *phd{fil + m};
size_t td{m >> 2};
size_t j{0u};
do {
/* f = fil + pf*phd */
const float32x4_t f4 = vmlaq_f32(vld1q_f32(&fil[j]), pf4, vld1q_f32(&phd[j]));
/* r += f*src */
r4 = vmlaq_f32(r4, f4, vld1q_f32(&src[j]));
j += 4;
} while(--td);
}
r4 = vaddq_f32(r4, vrev64q_f32(r4));
out_sample = vget_lane_f32(vadd_f32(vget_low_f32(r4), vget_high_f32(r4)), 0);
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
}
return dst.data();
}
template<>
void MixHrtf_<NEONTag>(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const MixHrtfFilter *hrtfparams, const size_t BufferSize)
{ MixHrtfBase<ApplyCoeffs>(InSamples, AccumSamples, IrSize, hrtfparams, BufferSize); }
template<>
void MixHrtfBlend_<NEONTag>(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const HrtfFilter *oldparams, const MixHrtfFilter *newparams, const size_t BufferSize)
{
MixHrtfBlendBase<ApplyCoeffs>(InSamples, AccumSamples, IrSize, oldparams, newparams,
BufferSize);
}
template<>
void MixDirectHrtf_<NEONTag>(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const al::span<const FloatBufferLine> InSamples, float2 *AccumSamples,
float *TempBuf, HrtfChannelState *ChanState, const size_t IrSize, const size_t BufferSize)
{
MixDirectHrtfBase<ApplyCoeffs>(LeftOut, RightOut, InSamples, AccumSamples, TempBuf, ChanState,
IrSize, BufferSize);
}
template<>
void Mix_<NEONTag>(const al::span<const float> InSamples, const al::span<FloatBufferLine> OutBuffer,
float *CurrentGains, const float *TargetGains, const size_t Counter, const size_t OutPos)
{
const float delta{(Counter > 0) ? 1.0f / static_cast<float>(Counter) : 0.0f};
const auto min_len = minz(Counter, InSamples.size());
const auto aligned_len = minz((min_len+3) & ~size_t{3}, InSamples.size()) - min_len;
for(FloatBufferLine &output : OutBuffer)
{
float *RESTRICT dst{al::assume_aligned<16>(output.data()+OutPos)};
float gain{*CurrentGains};
const float step{(*TargetGains-gain) * delta};
size_t pos{0};
if(!(std::abs(step) > std::numeric_limits<float>::epsilon()))
gain = *TargetGains;
else
{
float step_count{0.0f};
/* Mix with applying gain steps in aligned multiples of 4. */
if(size_t todo{(min_len-pos) >> 2})
{
const float32x4_t four4{vdupq_n_f32(4.0f)};
const float32x4_t step4{vdupq_n_f32(step)};
const float32x4_t gain4{vdupq_n_f32(gain)};
float32x4_t step_count4{vdupq_n_f32(0.0f)};
step_count4 = vsetq_lane_f32(1.0f, step_count4, 1);
step_count4 = vsetq_lane_f32(2.0f, step_count4, 2);
step_count4 = vsetq_lane_f32(3.0f, step_count4, 3);
do {
const float32x4_t val4 = vld1q_f32(&InSamples[pos]);
float32x4_t dry4 = vld1q_f32(&dst[pos]);
dry4 = vmlaq_f32(dry4, val4, vmlaq_f32(gain4, step4, step_count4));
step_count4 = vaddq_f32(step_count4, four4);
vst1q_f32(&dst[pos], dry4);
pos += 4;
} while(--todo);
/* NOTE: step_count4 now represents the next four counts after
* the last four mixed samples, so the lowest element
* represents the next step count to apply.
*/
step_count = vgetq_lane_f32(step_count4, 0);
}
/* Mix with applying left over gain steps that aren't aligned multiples of 4. */
for(size_t leftover{min_len&3};leftover;++pos,--leftover)
{
dst[pos] += InSamples[pos] * (gain + step*step_count);
step_count += 1.0f;
}
if(pos == Counter)
gain = *TargetGains;
else
gain += step*step_count;
/* Mix until pos is aligned with 4 or the mix is done. */
for(size_t leftover{aligned_len&3};leftover;++pos,--leftover)
dst[pos] += InSamples[pos] * gain;
}
*CurrentGains = gain;
++CurrentGains;
++TargetGains;
if(!(std::abs(gain) > GainSilenceThreshold))
continue;
if(size_t todo{(InSamples.size()-pos) >> 2})
{
const float32x4_t gain4 = vdupq_n_f32(gain);
do {
const float32x4_t val4 = vld1q_f32(&InSamples[pos]);
float32x4_t dry4 = vld1q_f32(&dst[pos]);
dry4 = vmlaq_f32(dry4, val4, gain4);
vst1q_f32(&dst[pos], dry4);
pos += 4;
} while(--todo);
}
for(size_t leftover{(InSamples.size()-pos)&3};leftover;++pos,--leftover)
dst[pos] += InSamples[pos] * gain;
}
}

View file

@ -0,0 +1,271 @@
#include "config.h"
#include <xmmintrin.h>
#include <cmath>
#include <limits>
#include "alnumeric.h"
#include "core/bsinc_defs.h"
#include "defs.h"
#include "hrtfbase.h"
struct SSETag;
struct BSincTag;
struct FastBSincTag;
/* SSE2 is required for any SSE support. */
#if defined(__GNUC__) && !defined(__clang__) && !defined(__SSE2__)
#pragma GCC target("sse2")
#endif
namespace {
constexpr uint FracPhaseBitDiff{MixerFracBits - BSincPhaseBits};
constexpr uint FracPhaseDiffOne{1 << FracPhaseBitDiff};
#define MLA4(x, y, z) _mm_add_ps(x, _mm_mul_ps(y, z))
inline void ApplyCoeffs(float2 *RESTRICT Values, const size_t IrSize, const HrirArray &Coeffs,
const float left, const float right)
{
const __m128 lrlr{_mm_setr_ps(left, right, left, right)};
ASSUME(IrSize >= MinIrLength);
/* This isn't technically correct to test alignment, but it's true for
* systems that support SSE, which is the only one that needs to know the
* alignment of Values (which alternates between 8- and 16-byte aligned).
*/
if(reinterpret_cast<intptr_t>(Values)&0x8)
{
__m128 imp0, imp1;
__m128 coeffs{_mm_load_ps(&Coeffs[0][0])};
__m128 vals{_mm_loadl_pi(_mm_setzero_ps(), reinterpret_cast<__m64*>(&Values[0][0]))};
imp0 = _mm_mul_ps(lrlr, coeffs);
vals = _mm_add_ps(imp0, vals);
_mm_storel_pi(reinterpret_cast<__m64*>(&Values[0][0]), vals);
size_t td{((IrSize+1)>>1) - 1};
size_t i{1};
do {
coeffs = _mm_load_ps(&Coeffs[i+1][0]);
vals = _mm_load_ps(&Values[i][0]);
imp1 = _mm_mul_ps(lrlr, coeffs);
imp0 = _mm_shuffle_ps(imp0, imp1, _MM_SHUFFLE(1, 0, 3, 2));
vals = _mm_add_ps(imp0, vals);
_mm_store_ps(&Values[i][0], vals);
imp0 = imp1;
i += 2;
} while(--td);
vals = _mm_loadl_pi(vals, reinterpret_cast<__m64*>(&Values[i][0]));
imp0 = _mm_movehl_ps(imp0, imp0);
vals = _mm_add_ps(imp0, vals);
_mm_storel_pi(reinterpret_cast<__m64*>(&Values[i][0]), vals);
}
else
{
for(size_t i{0};i < IrSize;i += 2)
{
const __m128 coeffs{_mm_load_ps(&Coeffs[i][0])};
__m128 vals{_mm_load_ps(&Values[i][0])};
vals = MLA4(vals, lrlr, coeffs);
_mm_store_ps(&Values[i][0], vals);
}
}
}
} // namespace
template<>
float *Resample_<BSincTag,SSETag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const float *const filter{state->bsinc.filter};
const __m128 sf4{_mm_set1_ps(state->bsinc.sf)};
const size_t m{state->bsinc.m};
src -= state->bsinc.l;
for(float &out_sample : dst)
{
// Calculate the phase index and factor.
const uint pi{frac >> FracPhaseBitDiff};
const float pf{static_cast<float>(frac & (FracPhaseDiffOne-1)) * (1.0f/FracPhaseDiffOne)};
// Apply the scale and phase interpolated filter.
__m128 r4{_mm_setzero_ps()};
{
const __m128 pf4{_mm_set1_ps(pf)};
const float *fil{filter + m*pi*4};
const float *phd{fil + m};
const float *scd{phd + m};
const float *spd{scd + m};
size_t td{m >> 2};
size_t j{0u};
do {
/* f = ((fil + sf*scd) + pf*(phd + sf*spd)) */
const __m128 f4 = MLA4(
MLA4(_mm_load_ps(&fil[j]), sf4, _mm_load_ps(&scd[j])),
pf4, MLA4(_mm_load_ps(&phd[j]), sf4, _mm_load_ps(&spd[j])));
/* r += f*src */
r4 = MLA4(r4, f4, _mm_loadu_ps(&src[j]));
j += 4;
} while(--td);
}
r4 = _mm_add_ps(r4, _mm_shuffle_ps(r4, r4, _MM_SHUFFLE(0, 1, 2, 3)));
r4 = _mm_add_ps(r4, _mm_movehl_ps(r4, r4));
out_sample = _mm_cvtss_f32(r4);
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
}
return dst.data();
}
template<>
float *Resample_<FastBSincTag,SSETag>(const InterpState *state, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const float *const filter{state->bsinc.filter};
const size_t m{state->bsinc.m};
src -= state->bsinc.l;
for(float &out_sample : dst)
{
// Calculate the phase index and factor.
const uint pi{frac >> FracPhaseBitDiff};
const float pf{static_cast<float>(frac & (FracPhaseDiffOne-1)) * (1.0f/FracPhaseDiffOne)};
// Apply the phase interpolated filter.
__m128 r4{_mm_setzero_ps()};
{
const __m128 pf4{_mm_set1_ps(pf)};
const float *fil{filter + m*pi*4};
const float *phd{fil + m};
size_t td{m >> 2};
size_t j{0u};
do {
/* f = fil + pf*phd */
const __m128 f4 = MLA4(_mm_load_ps(&fil[j]), pf4, _mm_load_ps(&phd[j]));
/* r += f*src */
r4 = MLA4(r4, f4, _mm_loadu_ps(&src[j]));
j += 4;
} while(--td);
}
r4 = _mm_add_ps(r4, _mm_shuffle_ps(r4, r4, _MM_SHUFFLE(0, 1, 2, 3)));
r4 = _mm_add_ps(r4, _mm_movehl_ps(r4, r4));
out_sample = _mm_cvtss_f32(r4);
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
}
return dst.data();
}
template<>
void MixHrtf_<SSETag>(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const MixHrtfFilter *hrtfparams, const size_t BufferSize)
{ MixHrtfBase<ApplyCoeffs>(InSamples, AccumSamples, IrSize, hrtfparams, BufferSize); }
template<>
void MixHrtfBlend_<SSETag>(const float *InSamples, float2 *AccumSamples, const uint IrSize,
const HrtfFilter *oldparams, const MixHrtfFilter *newparams, const size_t BufferSize)
{
MixHrtfBlendBase<ApplyCoeffs>(InSamples, AccumSamples, IrSize, oldparams, newparams,
BufferSize);
}
template<>
void MixDirectHrtf_<SSETag>(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const al::span<const FloatBufferLine> InSamples, float2 *AccumSamples,
float *TempBuf, HrtfChannelState *ChanState, const size_t IrSize, const size_t BufferSize)
{
MixDirectHrtfBase<ApplyCoeffs>(LeftOut, RightOut, InSamples, AccumSamples, TempBuf, ChanState,
IrSize, BufferSize);
}
template<>
void Mix_<SSETag>(const al::span<const float> InSamples, const al::span<FloatBufferLine> OutBuffer,
float *CurrentGains, const float *TargetGains, const size_t Counter, const size_t OutPos)
{
const float delta{(Counter > 0) ? 1.0f / static_cast<float>(Counter) : 0.0f};
const auto min_len = minz(Counter, InSamples.size());
const auto aligned_len = minz((min_len+3) & ~size_t{3}, InSamples.size()) - min_len;
for(FloatBufferLine &output : OutBuffer)
{
float *RESTRICT dst{al::assume_aligned<16>(output.data()+OutPos)};
float gain{*CurrentGains};
const float step{(*TargetGains-gain) * delta};
size_t pos{0};
if(!(std::abs(step) > std::numeric_limits<float>::epsilon()))
gain = *TargetGains;
else
{
float step_count{0.0f};
/* Mix with applying gain steps in aligned multiples of 4. */
if(size_t todo{(min_len-pos) >> 2})
{
const __m128 four4{_mm_set1_ps(4.0f)};
const __m128 step4{_mm_set1_ps(step)};
const __m128 gain4{_mm_set1_ps(gain)};
__m128 step_count4{_mm_setr_ps(0.0f, 1.0f, 2.0f, 3.0f)};
do {
const __m128 val4{_mm_load_ps(&InSamples[pos])};
__m128 dry4{_mm_load_ps(&dst[pos])};
/* dry += val * (gain + step*step_count) */
dry4 = MLA4(dry4, val4, MLA4(gain4, step4, step_count4));
_mm_store_ps(&dst[pos], dry4);
step_count4 = _mm_add_ps(step_count4, four4);
pos += 4;
} while(--todo);
/* NOTE: step_count4 now represents the next four counts after
* the last four mixed samples, so the lowest element
* represents the next step count to apply.
*/
step_count = _mm_cvtss_f32(step_count4);
}
/* Mix with applying left over gain steps that aren't aligned multiples of 4. */
for(size_t leftover{min_len&3};leftover;++pos,--leftover)
{
dst[pos] += InSamples[pos] * (gain + step*step_count);
step_count += 1.0f;
}
if(pos == Counter)
gain = *TargetGains;
else
gain += step*step_count;
/* Mix until pos is aligned with 4 or the mix is done. */
for(size_t leftover{aligned_len&3};leftover;++pos,--leftover)
dst[pos] += InSamples[pos] * gain;
}
*CurrentGains = gain;
++CurrentGains;
++TargetGains;
if(!(std::abs(gain) > GainSilenceThreshold))
continue;
if(size_t todo{(InSamples.size()-pos) >> 2})
{
const __m128 gain4{_mm_set1_ps(gain)};
do {
const __m128 val4{_mm_load_ps(&InSamples[pos])};
__m128 dry4{_mm_load_ps(&dst[pos])};
dry4 = _mm_add_ps(dry4, _mm_mul_ps(val4, gain4));
_mm_store_ps(&dst[pos], dry4);
pos += 4;
} while(--todo);
}
for(size_t leftover{(InSamples.size()-pos)&3};leftover;++pos,--leftover)
dst[pos] += InSamples[pos] * gain;
}
}

View file

@ -0,0 +1,89 @@
/**
* OpenAL cross platform audio library
* Copyright (C) 2014 by Timothy Arceri <t_arceri@yahoo.com.au>.
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
* Or go to http://www.gnu.org/copyleft/lgpl.html
*/
#include "config.h"
#include <xmmintrin.h>
#include <emmintrin.h>
#include "alnumeric.h"
#include "defs.h"
struct SSE2Tag;
struct LerpTag;
#if defined(__GNUC__) && !defined(__clang__) && !defined(__SSE2__)
#pragma GCC target("sse2")
#endif
template<>
float *Resample_<LerpTag,SSE2Tag>(const InterpState*, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const __m128i increment4{_mm_set1_epi32(static_cast<int>(increment*4))};
const __m128 fracOne4{_mm_set1_ps(1.0f/MixerFracOne)};
const __m128i fracMask4{_mm_set1_epi32(MixerFracMask)};
alignas(16) uint pos_[4], frac_[4];
InitPosArrays(frac, increment, frac_, pos_);
__m128i frac4{_mm_setr_epi32(static_cast<int>(frac_[0]), static_cast<int>(frac_[1]),
static_cast<int>(frac_[2]), static_cast<int>(frac_[3]))};
__m128i pos4{_mm_setr_epi32(static_cast<int>(pos_[0]), static_cast<int>(pos_[1]),
static_cast<int>(pos_[2]), static_cast<int>(pos_[3]))};
auto dst_iter = dst.begin();
for(size_t todo{dst.size()>>2};todo;--todo)
{
const int pos0{_mm_cvtsi128_si32(_mm_shuffle_epi32(pos4, _MM_SHUFFLE(0, 0, 0, 0)))};
const int pos1{_mm_cvtsi128_si32(_mm_shuffle_epi32(pos4, _MM_SHUFFLE(1, 1, 1, 1)))};
const int pos2{_mm_cvtsi128_si32(_mm_shuffle_epi32(pos4, _MM_SHUFFLE(2, 2, 2, 2)))};
const int pos3{_mm_cvtsi128_si32(_mm_shuffle_epi32(pos4, _MM_SHUFFLE(3, 3, 3, 3)))};
const __m128 val1{_mm_setr_ps(src[pos0 ], src[pos1 ], src[pos2 ], src[pos3 ])};
const __m128 val2{_mm_setr_ps(src[pos0+1], src[pos1+1], src[pos2+1], src[pos3+1])};
/* val1 + (val2-val1)*mu */
const __m128 r0{_mm_sub_ps(val2, val1)};
const __m128 mu{_mm_mul_ps(_mm_cvtepi32_ps(frac4), fracOne4)};
const __m128 out{_mm_add_ps(val1, _mm_mul_ps(mu, r0))};
_mm_store_ps(dst_iter, out);
dst_iter += 4;
frac4 = _mm_add_epi32(frac4, increment4);
pos4 = _mm_add_epi32(pos4, _mm_srli_epi32(frac4, MixerFracBits));
frac4 = _mm_and_si128(frac4, fracMask4);
}
if(size_t todo{dst.size()&3})
{
src += static_cast<uint>(_mm_cvtsi128_si32(pos4));
frac = static_cast<uint>(_mm_cvtsi128_si32(frac4));
do {
*(dst_iter++) = lerp(src[0], src[1], static_cast<float>(frac) * (1.0f/MixerFracOne));
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
} while(--todo);
}
return dst.data();
}

View file

@ -0,0 +1,94 @@
/**
* OpenAL cross platform audio library
* Copyright (C) 2014 by Timothy Arceri <t_arceri@yahoo.com.au>.
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
* Or go to http://www.gnu.org/copyleft/lgpl.html
*/
#include "config.h"
#include <xmmintrin.h>
#include <emmintrin.h>
#include <smmintrin.h>
#include "alnumeric.h"
#include "defs.h"
struct SSE4Tag;
struct LerpTag;
#if defined(__GNUC__) && !defined(__clang__) && !defined(__SSE4_1__)
#pragma GCC target("sse4.1")
#endif
template<>
float *Resample_<LerpTag,SSE4Tag>(const InterpState*, float *RESTRICT src, uint frac,
uint increment, const al::span<float> dst)
{
const __m128i increment4{_mm_set1_epi32(static_cast<int>(increment*4))};
const __m128 fracOne4{_mm_set1_ps(1.0f/MixerFracOne)};
const __m128i fracMask4{_mm_set1_epi32(MixerFracMask)};
alignas(16) uint pos_[4], frac_[4];
InitPosArrays(frac, increment, frac_, pos_);
__m128i frac4{_mm_setr_epi32(static_cast<int>(frac_[0]), static_cast<int>(frac_[1]),
static_cast<int>(frac_[2]), static_cast<int>(frac_[3]))};
__m128i pos4{_mm_setr_epi32(static_cast<int>(pos_[0]), static_cast<int>(pos_[1]),
static_cast<int>(pos_[2]), static_cast<int>(pos_[3]))};
auto dst_iter = dst.begin();
for(size_t todo{dst.size()>>2};todo;--todo)
{
const int pos0{_mm_extract_epi32(pos4, 0)};
const int pos1{_mm_extract_epi32(pos4, 1)};
const int pos2{_mm_extract_epi32(pos4, 2)};
const int pos3{_mm_extract_epi32(pos4, 3)};
const __m128 val1{_mm_setr_ps(src[pos0 ], src[pos1 ], src[pos2 ], src[pos3 ])};
const __m128 val2{_mm_setr_ps(src[pos0+1], src[pos1+1], src[pos2+1], src[pos3+1])};
/* val1 + (val2-val1)*mu */
const __m128 r0{_mm_sub_ps(val2, val1)};
const __m128 mu{_mm_mul_ps(_mm_cvtepi32_ps(frac4), fracOne4)};
const __m128 out{_mm_add_ps(val1, _mm_mul_ps(mu, r0))};
_mm_store_ps(dst_iter, out);
dst_iter += 4;
frac4 = _mm_add_epi32(frac4, increment4);
pos4 = _mm_add_epi32(pos4, _mm_srli_epi32(frac4, MixerFracBits));
frac4 = _mm_and_si128(frac4, fracMask4);
}
if(size_t todo{dst.size()&3})
{
/* NOTE: These four elements represent the position *after* the last
* four samples, so the lowest element is the next position to
* resample.
*/
src += static_cast<uint>(_mm_cvtsi128_si32(pos4));
frac = static_cast<uint>(_mm_cvtsi128_si32(frac4));
do {
*(dst_iter++) = lerp(src[0], src[1], static_cast<float>(frac) * (1.0f/MixerFracOne));
frac += increment;
src += frac>>MixerFracBits;
frac &= MixerFracMask;
} while(--todo);
}
return dst.data();
}

View file

@ -0,0 +1,275 @@
#include "config.h"
#include "uhjfilter.h"
#ifdef HAVE_SSE_INTRINSICS
#include <xmmintrin.h>
#elif defined(HAVE_NEON)
#include <arm_neon.h>
#endif
#include <algorithm>
#include <iterator>
#include "alcomplex.h"
#include "alnumeric.h"
#include "opthelpers.h"
namespace {
using complex_d = std::complex<double>;
struct PhaseShifterT {
alignas(16) std::array<float,Uhj2Encoder::sFilterSize> Coeffs;
/* Some notes on this filter construction.
*
* A wide-band phase-shift filter needs a delay to maintain linearity. A
* dirac impulse in the center of a time-domain buffer represents a filter
* passing all frequencies through as-is with a pure delay. Converting that
* to the frequency domain, adjusting the phase of each frequency bin by
* +90 degrees, then converting back to the time domain, results in a FIR
* filter that applies a +90 degree wide-band phase-shift.
*
* A particularly notable aspect of the time-domain filter response is that
* every other coefficient is 0. This allows doubling the effective size of
* the filter, by storing only the non-0 coefficients and double-stepping
* over the input to apply it.
*
* Additionally, the resulting filter is independent of the sample rate.
* The same filter can be applied regardless of the device's sample rate
* and achieve the same effect.
*/
PhaseShifterT()
{
constexpr size_t fft_size{Uhj2Encoder::sFilterSize * 2};
constexpr size_t half_size{fft_size / 2};
/* Generate a frequency domain impulse with a +90 degree phase offset.
* Reconstruct the mirrored frequencies to convert to the time domain.
*/
auto fftBuffer = std::make_unique<complex_d[]>(fft_size);
std::fill_n(fftBuffer.get(), fft_size, complex_d{});
fftBuffer[half_size] = 1.0;
forward_fft({fftBuffer.get(), fft_size});
for(size_t i{0};i < half_size+1;++i)
fftBuffer[i] = complex_d{-fftBuffer[i].imag(), fftBuffer[i].real()};
for(size_t i{half_size+1};i < fft_size;++i)
fftBuffer[i] = std::conj(fftBuffer[fft_size - i]);
inverse_fft({fftBuffer.get(), fft_size});
/* Reverse the filter for simpler processing, and store only the non-0
* coefficients.
*/
auto fftiter = fftBuffer.get() + half_size + (Uhj2Encoder::sFilterSize-1);
for(float &coeff : Coeffs)
{
coeff = static_cast<float>(fftiter->real() / double{fft_size});
fftiter -= 2;
}
}
};
const PhaseShifterT PShift{};
void allpass_process(al::span<float> dst, const float *RESTRICT src)
{
#ifdef HAVE_SSE_INTRINSICS
size_t pos{0};
if(size_t todo{dst.size()>>1})
{
do {
__m128 r04{_mm_setzero_ps()};
__m128 r14{_mm_setzero_ps()};
for(size_t j{0};j < PShift.Coeffs.size();j+=4)
{
const __m128 coeffs{_mm_load_ps(&PShift.Coeffs[j])};
const __m128 s0{_mm_loadu_ps(&src[j*2])};
const __m128 s1{_mm_loadu_ps(&src[j*2 + 4])};
__m128 s{_mm_shuffle_ps(s0, s1, _MM_SHUFFLE(2, 0, 2, 0))};
r04 = _mm_add_ps(r04, _mm_mul_ps(s, coeffs));
s = _mm_shuffle_ps(s0, s1, _MM_SHUFFLE(3, 1, 3, 1));
r14 = _mm_add_ps(r14, _mm_mul_ps(s, coeffs));
}
r04 = _mm_add_ps(r04, _mm_shuffle_ps(r04, r04, _MM_SHUFFLE(0, 1, 2, 3)));
r04 = _mm_add_ps(r04, _mm_movehl_ps(r04, r04));
dst[pos++] += _mm_cvtss_f32(r04);
r14 = _mm_add_ps(r14, _mm_shuffle_ps(r14, r14, _MM_SHUFFLE(0, 1, 2, 3)));
r14 = _mm_add_ps(r14, _mm_movehl_ps(r14, r14));
dst[pos++] += _mm_cvtss_f32(r14);
src += 2;
} while(--todo);
}
if((dst.size()&1))
{
__m128 r4{_mm_setzero_ps()};
for(size_t j{0};j < PShift.Coeffs.size();j+=4)
{
const __m128 coeffs{_mm_load_ps(&PShift.Coeffs[j])};
/* NOTE: This could alternatively be done with two unaligned loads
* and a shuffle. Which would be better?
*/
const __m128 s{_mm_setr_ps(src[j*2], src[j*2 + 2], src[j*2 + 4], src[j*2 + 6])};
r4 = _mm_add_ps(r4, _mm_mul_ps(s, coeffs));
}
r4 = _mm_add_ps(r4, _mm_shuffle_ps(r4, r4, _MM_SHUFFLE(0, 1, 2, 3)));
r4 = _mm_add_ps(r4, _mm_movehl_ps(r4, r4));
dst[pos] += _mm_cvtss_f32(r4);
}
#elif defined(HAVE_NEON)
size_t pos{0};
if(size_t todo{dst.size()>>1})
{
/* There doesn't seem to be NEON intrinsics to do this kind of stipple
* shuffling, so there's two custom methods for it.
*/
auto shuffle_2020 = [](float32x4_t a, float32x4_t b)
{
float32x4_t ret{vmovq_n_f32(vgetq_lane_f32(a, 0))};
ret = vsetq_lane_f32(vgetq_lane_f32(a, 2), ret, 1);
ret = vsetq_lane_f32(vgetq_lane_f32(b, 0), ret, 2);
ret = vsetq_lane_f32(vgetq_lane_f32(b, 2), ret, 3);
return ret;
};
auto shuffle_3131 = [](float32x4_t a, float32x4_t b)
{
float32x4_t ret{vmovq_n_f32(vgetq_lane_f32(a, 1))};
ret = vsetq_lane_f32(vgetq_lane_f32(a, 3), ret, 1);
ret = vsetq_lane_f32(vgetq_lane_f32(b, 1), ret, 2);
ret = vsetq_lane_f32(vgetq_lane_f32(b, 3), ret, 3);
return ret;
};
do {
float32x4_t r04{vdupq_n_f32(0.0f)};
float32x4_t r14{vdupq_n_f32(0.0f)};
for(size_t j{0};j < PShift.Coeffs.size();j+=4)
{
const float32x4_t coeffs{vld1q_f32(&PShift.Coeffs[j])};
const float32x4_t s0{vld1q_f32(&src[j*2])};
const float32x4_t s1{vld1q_f32(&src[j*2 + 4])};
r04 = vmlaq_f32(r04, shuffle_2020(s0, s1), coeffs);
r14 = vmlaq_f32(r14, shuffle_3131(s0, s1), coeffs);
}
r04 = vaddq_f32(r04, vrev64q_f32(r04));
dst[pos++] = vget_lane_f32(vadd_f32(vget_low_f32(r04), vget_high_f32(r04)), 0);
r14 = vaddq_f32(r14, vrev64q_f32(r14));
dst[pos++] = vget_lane_f32(vadd_f32(vget_low_f32(r14), vget_high_f32(r14)), 0);
src += 2;
} while(--todo);
}
if((dst.size()&1))
{
auto load4 = [](float32_t a, float32_t b, float32_t c, float32_t d)
{
float32x4_t ret{vmovq_n_f32(a)};
ret = vsetq_lane_f32(b, ret, 1);
ret = vsetq_lane_f32(c, ret, 2);
ret = vsetq_lane_f32(d, ret, 3);
return ret;
};
float32x4_t r4{vdupq_n_f32(0.0f)};
for(size_t j{0};j < PShift.Coeffs.size();j+=4)
{
const float32x4_t coeffs{vld1q_f32(&PShift.Coeffs[j])};
const float32x4_t s{load4(src[j*2], src[j*2 + 2], src[j*2 + 4], src[j*2 + 6])};
r4 = vmlaq_f32(r4, s, coeffs);
}
r4 = vaddq_f32(r4, vrev64q_f32(r4));
dst[pos] = vget_lane_f32(vadd_f32(vget_low_f32(r4), vget_high_f32(r4)), 0);
}
#else
for(float &output : dst)
{
float ret{0.0f};
for(size_t j{0};j < PShift.Coeffs.size();++j)
ret += src[j*2] * PShift.Coeffs[j];
output += ret;
++src;
}
#endif
}
} // namespace
/* Encoding 2-channel UHJ from B-Format is done as:
*
* S = 0.9396926*W + 0.1855740*X
* D = j(-0.3420201*W + 0.5098604*X) + 0.6554516*Y
*
* Left = (S + D)/2.0
* Right = (S - D)/2.0
*
* where j is a wide-band +90 degree phase shift.
*
* The phase shift is done using a FIR filter derived from an FFT'd impulse
* with the desired shift.
*/
void Uhj2Encoder::encode(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const FloatBufferLine *InSamples, const size_t SamplesToDo)
{
ASSUME(SamplesToDo > 0);
float *RESTRICT left{al::assume_aligned<16>(LeftOut.data())};
float *RESTRICT right{al::assume_aligned<16>(RightOut.data())};
const float *RESTRICT winput{al::assume_aligned<16>(InSamples[0].data())};
const float *RESTRICT xinput{al::assume_aligned<16>(InSamples[1].data())};
const float *RESTRICT yinput{al::assume_aligned<16>(InSamples[2].data())};
/* Combine the previously delayed mid/side signal with the input. */
/* S = 0.9396926*W + 0.1855740*X */
auto miditer = std::copy(mMidDelay.cbegin(), mMidDelay.cend(), mMid.begin());
std::transform(winput, winput+SamplesToDo, xinput, miditer,
[](const float w, const float x) noexcept -> float
{ return 0.9396926f*w + 0.1855740f*x; });
/* D = 0.6554516*Y */
auto sideiter = std::copy(mSideDelay.cbegin(), mSideDelay.cend(), mSide.begin());
std::transform(yinput, yinput+SamplesToDo, sideiter,
[](const float y) noexcept -> float { return 0.6554516f*y; });
/* Include any existing direct signal in the mid/side buffers. */
for(size_t i{0};i < SamplesToDo;++i,++miditer)
*miditer += left[i] + right[i];
for(size_t i{0};i < SamplesToDo;++i,++sideiter)
*sideiter += left[i] - right[i];
/* Copy the future samples back to the delay buffers for next time. */
std::copy_n(mMid.cbegin()+SamplesToDo, mMidDelay.size(), mMidDelay.begin());
std::copy_n(mSide.cbegin()+SamplesToDo, mSideDelay.size(), mSideDelay.begin());
/* Now add the all-passed signal into the side signal. */
/* D += j(-0.3420201*W + 0.5098604*X) */
auto tmpiter = std::copy(mSideHistory.cbegin(), mSideHistory.cend(), mTemp.begin());
std::transform(winput, winput+SamplesToDo, xinput, tmpiter,
[](const float w, const float x) noexcept -> float
{ return -0.3420201f*w + 0.5098604f*x; });
std::copy_n(mTemp.cbegin()+SamplesToDo, mSideHistory.size(), mSideHistory.begin());
allpass_process({mSide.data(), SamplesToDo}, mTemp.data());
/* Left = (S + D)/2.0 */
for(size_t i{0};i < SamplesToDo;i++)
left[i] = (mMid[i] + mSide[i]) * 0.5f;
/* Right = (S - D)/2.0 */
for(size_t i{0};i < SamplesToDo;i++)
right[i] = (mMid[i] - mSide[i]) * 0.5f;
}

View file

@ -0,0 +1,39 @@
#ifndef CORE_UHJFILTER_H
#define CORE_UHJFILTER_H
#include <array>
#include "almalloc.h"
#include "bufferline.h"
struct Uhj2Encoder {
/* A particular property of the filter allows it to cover nearly twice its
* length, so the filter size is also the effective delay (despite being
* center-aligned).
*/
constexpr static size_t sFilterSize{128};
/* Delays for the unfiltered signal. */
alignas(16) std::array<float,sFilterSize> mMidDelay{};
alignas(16) std::array<float,sFilterSize> mSideDelay{};
alignas(16) std::array<float,BufferLineSize+sFilterSize> mMid{};
alignas(16) std::array<float,BufferLineSize+sFilterSize> mSide{};
/* History for the FIR filter. */
alignas(16) std::array<float,sFilterSize*2 - 1> mSideHistory{};
alignas(16) std::array<float,BufferLineSize + sFilterSize*2> mTemp{};
/**
* Encodes a 2-channel UHJ (stereo-compatible) signal from a B-Format input
* signal. The input must use FuMa channel ordering and scaling.
*/
void encode(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const FloatBufferLine *InSamples, const size_t SamplesToDo);
DEF_NEWDEL(Uhj2Encoder)
};
#endif /* CORE_UHJFILTER_H */