Init foo_input_joc
build / windows (push) Has been cancelled

This commit is contained in:
2026-09-25 17:07:04 +08:00
commit ad62a914e6
127 changed files with 38545 additions and 0 deletions
+401
View File
@@ -0,0 +1,401 @@
#include "hrtf/jochrtf.h"
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include "foundation/mini_json.h"
#include "foundation/sha256.h"
#include "io/npy.h"
#include "io/zip_reader.h"
namespace joc::hrtf {
namespace {
std::string to_upper(std::string text) {
for (char& c : text) {
if (c >= 'a' && c <= 'z') {
c = static_cast<char>(c - 'a' + 'A');
}
}
return text;
}
bool is_sha256_hex(const std::string& text) {
if (text.size() != 64) {
return false;
}
for (const char c : text) {
const bool digit = c >= '0' && c <= '9';
const bool upper = c >= 'A' && c <= 'F';
if (!digit && !upper) {
return false;
}
}
return true;
}
// json.dumps(list(shape)) as the reference writes it, e.g. "[36, 2, 77]".
std::string shape_json(const std::vector<std::int64_t>& shape) {
std::string text = "[";
for (std::size_t i = 0; i < shape.size(); ++i) {
text += (i == 0 ? "" : ", ");
text += std::to_string(shape[i]);
}
text += "]";
return text;
}
Status hrtf_fail(const std::string& message) {
return Status::fail(JOC_ERR_HRTF_FORMAT, stage::kRender, message);
}
std::string payload_sha256(const std::vector<double>& centers,
const std::vector<double>& coefficients,
const std::vector<double>& delay_coefficients,
const std::vector<double>& delay_bounds) {
crypto::Sha256 hash;
const char prefix[] = "JOC-HRTF-CACHE-PAYLOAD-V1";
hash.update(prefix, sizeof(prefix) - 1);
const std::uint8_t zero = 0;
hash.update(&zero, 1);
struct Entry {
const char* name;
const char* dtype;
const std::vector<double>* values;
std::vector<std::int64_t> shape;
};
const Entry entries[4] = {
{"band_center_frequencies_hz", "<f8", &centers, {kHybridBands}},
{"coefficients", "<c16", &coefficients, {kShTerms, kEars, kHybridBands}},
{"delay_coefficients", "<f8", &delay_coefficients, {kShTerms, kEars}},
{"delay_bounds", "<f8", &delay_bounds, {2, 2}},
};
for (const Entry& entry : entries) {
const std::string name(entry.name);
const std::string dtype(entry.dtype);
const std::string shape = shape_json(entry.shape);
hash.update(name.data(), name.size());
hash.update(&zero, 1);
hash.update(dtype.data(), dtype.size());
hash.update(&zero, 1);
hash.update(shape.data(), shape.size());
hash.update(&zero, 1);
hash.update(entry.values->data(), entry.values->size() * sizeof(double));
}
return hash.finish_hex();
}
} // namespace
Status load_jochrtf(const std::string& path, Field* out) {
if (out == nullptr) {
return Status::fail(JOC_ERR_INVALID_ARGUMENT, stage::kRender, "null field");
}
io::ZipArchive archive;
std::string error;
if (!archive.open(path, &error)) {
return Status::fail(JOC_ERR_HRTF_NOT_FOUND, stage::kRender,
"cannot read compiled HRTF " + path + ": " + error);
}
// Member set must be exactly the five expected names.
static const char* kMembers[5] = {"metadata_json.npy", "band_center_frequencies_hz.npy",
"coefficients.npy", "delay_coefficients.npy",
"delay_bounds.npy"};
if (archive.entries().size() != 5u) {
return hrtf_fail("compiled HRTF cache has an invalid member set (" +
std::to_string(archive.entries().size()) + " members)");
}
for (const char* name : kMembers) {
if (archive.find(name) == nullptr) {
return hrtf_fail(std::string("compiled HRTF cache is missing ") + name);
}
}
auto read_member = [&](const char* name, std::vector<std::uint8_t>* raw,
io::NpyArray* array) -> Status {
if (!archive.read_member(name, raw, &error)) {
return hrtf_fail(std::string("compiled HRTF member ") + name + ": " + error);
}
if (!io::parse_npy(raw->data(), raw->size(), array, &error)) {
return hrtf_fail(std::string("compiled HRTF member ") + name + ": " + error);
}
if (array->fortran_order) {
return hrtf_fail(std::string("compiled HRTF member must be C-contiguous: ") + name);
}
return Status::success();
};
std::vector<std::uint8_t> raw;
io::NpyArray array;
Status status = read_member("metadata_json.npy", &raw, &array);
if (!status.ok()) {
return status;
}
std::string metadata_text;
if (!io::npy_unicode_to_utf8(array, &metadata_text, &error)) {
return hrtf_fail("compiled HRTF metadata: " + error);
}
if (metadata_text.size() > 64u * 1024u) {
return hrtf_fail("compiled HRTF metadata is too large");
}
status = read_member("band_center_frequencies_hz.npy", &raw, &array);
if (!status.ok()) {
return status;
}
if (array.descr != "<f8" || !io::npy_shape_is(array, {kHybridBands})) {
return hrtf_fail("band_center_frequencies_hz must be <f8(77,)");
}
std::vector<double> centers;
io::npy_to_double(array, &centers, &error);
status = read_member("coefficients.npy", &raw, &array);
if (!status.ok()) {
return status;
}
if (array.descr != "<c16" || !io::npy_shape_is(array, {kShTerms, kEars, kHybridBands})) {
return hrtf_fail("coefficients must be <c16(36, 2, 77)");
}
std::vector<double> coefficients;
if (!io::npy_to_double(array, &coefficients, &error)) {
return hrtf_fail("coefficients: " + error);
}
status = read_member("delay_coefficients.npy", &raw, &array);
if (!status.ok()) {
return status;
}
if (array.descr != "<f8" || !io::npy_shape_is(array, {kShTerms, kEars})) {
return hrtf_fail("delay_coefficients must be <f8(36, 2)");
}
std::vector<double> delay_coefficients;
io::npy_to_double(array, &delay_coefficients, &error);
status = read_member("delay_bounds.npy", &raw, &array);
if (!status.ok()) {
return status;
}
if (array.descr != "<f8" || !io::npy_shape_is(array, {2, 2})) {
return hrtf_fail("delay_bounds must be <f8(2, 2)");
}
std::vector<double> delay_bounds;
io::npy_to_double(array, &delay_bounds, &error);
std::vector<json::Member> members;
if (!json::parse_object(metadata_text, &members, &error)) {
return hrtf_fail("compiled HRTF metadata: " + error);
}
auto require_string = [&](const char* key, std::string* value) -> Status {
const json::Member* member = json::find(members, key);
if (member == nullptr || !json::as_string(*member, value)) {
return hrtf_fail(std::string("compiled HRTF metadata is missing ") + key);
}
return Status::success();
};
std::string magic;
std::string schema;
std::string source_sha256;
std::string cache_key;
std::string payload_hash;
std::string delay_source;
status = require_string("magic", &magic);
if (!status.ok()) { return status; }
status = require_string("cache_schema", &schema);
if (!status.ok()) { return status; }
status = require_string("source_sha256", &source_sha256);
if (!status.ok()) { return status; }
status = require_string("cache_key", &cache_key);
if (!status.ok()) { return status; }
status = require_string("payload_sha256", &payload_hash);
if (!status.ok()) { return status; }
status = require_string("delay_source", &delay_source);
if (!status.ok()) { return status; }
if (magic != kMagic) {
return hrtf_fail("compiled HRTF magic mismatch: " + magic);
}
if (schema != kCacheSchema) {
return hrtf_fail("compiled HRTF cache schema mismatch: " + schema);
}
const json::Member* version_member = json::find(members, "format_version");
long long version = -1;
if (version_member == nullptr || !json::as_integer(*version_member, &version)) {
return hrtf_fail("compiled HRTF metadata is missing format_version");
}
if (version != kFormatVersion) {
return Status::fail(JOC_ERR_HRTF_VERSION, stage::kRender,
"unsupported .jochrtf version " + std::to_string(version) +
"; rebuild it from the source SOFA");
}
out->source_sha256 = to_upper(source_sha256);
out->cache_key = to_upper(cache_key);
if (!is_sha256_hex(out->source_sha256)) {
return hrtf_fail("compiled HRTF source_sha256 is not a 64-digit digest");
}
if (!is_sha256_hex(out->cache_key)) {
return hrtf_fail("compiled HRTF cache_key is not a 64-digit digest");
}
const std::string expected = payload_sha256(centers, coefficients, delay_coefficients,
delay_bounds);
if (to_upper(payload_hash) != to_upper(expected)) {
return Status::fail(JOC_ERR_HRTF_HASH, stage::kRender,
"compiled HRTF payload hash mismatch");
}
out->payload_sha256 = to_upper(payload_hash);
const json::Member* radius_member = json::find(members, "measurement_radius_m");
double radius = 0.0;
if (radius_member == nullptr || !json::as_number(*radius_member, &radius) || radius <= 0.0) {
return hrtf_fail("compiled HRTF measurement_radius_m must be a positive number");
}
out->measurement_radius_m = radius;
const json::Member* order_member = json::find(members, "order");
long long order = 0;
if (order_member == nullptr || !json::as_integer(*order_member, &order) || order <= 0 ||
order * order > kShTerms) {
return hrtf_fail("compiled HRTF order is out of range");
}
out->order = order;
for (const double value : coefficients) {
if (!std::isfinite(value)) {
return hrtf_fail("compiled HRTF coefficients contain non-finite values");
}
}
for (const double value : delay_coefficients) {
if (!std::isfinite(value) || std::abs(value) > 48000.0 * 64.0) {
return hrtf_fail("compiled HRTF delay coefficients are out of range");
}
}
for (const double value : delay_bounds) {
if (!std::isfinite(value)) {
return hrtf_fail("compiled HRTF delay bounds contain non-finite values");
}
}
if (delay_bounds.size() == 4u && delay_bounds[0] > delay_bounds[1]) {
return hrtf_fail("compiled HRTF delay bounds are inverted");
}
if (const json::Member* member = json::find(members, "compiler_version")) {
json::as_string(*member, &out->compiler_version);
}
if (const json::Member* member = json::find(members, "phase_policy_version")) {
json::as_string(*member, &out->phase_policy_version);
}
if (const json::Member* member = json::find(members, "sh_convention")) {
json::as_string(*member, &out->sh_convention);
}
if (const json::Member* member = json::find(members, "source_display_name")) {
json::as_string(*member, &out->source_display_name);
}
if (const json::Member* member = json::find(members, "projection_ridge")) {
json::as_number(*member, &out->projection_ridge);
}
if (const json::Member* member = json::find(members, "spherical_harmonic_ridge")) {
json::as_number(*member, &out->spherical_harmonic_ridge);
}
if (const json::Member* member = json::find(members, "fit_report")) {
out->fit_report_json = member->raw;
}
if (const json::Member* member = json::find(members, "filterbank")) {
out->filterbank_json = member->raw;
}
out->delay_source = delay_source;
out->metadata_json = metadata_text;
out->coefficients = std::move(coefficients);
out->delay_coefficients = std::move(delay_coefficients);
out->delay_bounds = std::move(delay_bounds);
out->band_centers_hz = std::move(centers);
return Status::success();
}
Status load_kernels(const std::string& npz_path, Kernels* out) {
if (out == nullptr) {
return Status::fail(JOC_ERR_INVALID_ARGUMENT, stage::kRender, "null kernels");
}
io::ZipArchive archive;
std::string error;
if (!archive.open(npz_path, &error)) {
return Status::fail(JOC_ERR_HRTF_NOT_FOUND, stage::kRender,
"cannot read kernel tables " + npz_path + ": " + error);
}
struct Request {
const char* member;
const char* shape_text;
std::vector<std::int64_t> shape;
};
const Request requests[6] = {
{"qmf_analysis_coefficients.npy", "<f4", {64, 10}},
{"hybrid_analysis_low_kernel.npy", "<f4", {3, 2, 13, 16, 2}},
{"hybrid_synthesis_indices.npy", "<i2", {154, 4}},
{"hybrid_synthesis_values.npy", "<f4", {154}},
{"qmf_synthesis_basis.npy", "<f8", {64, 4, 128}},
{"qmf_synthesis_taps.npy", "<f8", {64, 10, 4}},
};
std::vector<std::uint8_t> raw;
std::vector<std::uint8_t> ordered;
for (const Request& request : requests) {
const std::string name = request.member;
if (!archive.read_member(name, &raw, &error)) {
return Status::fail(JOC_ERR_HRTF_FORMAT, stage::kRender,
"kernel table member " + name + ": " + error);
}
io::NpyArray array;
if (!io::parse_npy(raw.data(), raw.size(), &array, &error)) {
return Status::fail(JOC_ERR_HRTF_FORMAT, stage::kRender,
"kernel table member " + name + ": " + error);
}
if (array.descr != request.shape_text || !io::npy_shape_is(array, request.shape)) {
return Status::fail(JOC_ERR_HRTF_FORMAT, stage::kRender,
"kernel table member " + name + " has an unexpected dtype/shape");
}
// Logical C order: required because the reused kernel indexes the hybrid
// synthesis table row-major while the shipped member is Fortran-order.
if (!io::npy_to_c_order(array, &ordered, &error)) {
return Status::fail(JOC_ERR_HRTF_FORMAT, stage::kRender,
"kernel table member " + name + ": " + error);
}
const std::size_t count = array.element_count();
if (std::strcmp(request.member, "qmf_analysis_coefficients.npy") == 0) {
std::vector<float> values(count);
std::memcpy(values.data(), ordered.data(), count * sizeof(float));
out->qmf_analysis.assign(values.begin(), values.end());
} else if (std::strcmp(request.member, "hybrid_analysis_low_kernel.npy") == 0) {
std::vector<float> values(count);
std::memcpy(values.data(), ordered.data(), count * sizeof(float));
out->hybrid_low.assign(values.begin(), values.end());
} else if (std::strcmp(request.member, "hybrid_synthesis_indices.npy") == 0) {
out->hybrid_indices.resize(count);
std::memcpy(out->hybrid_indices.data(), ordered.data(), count * sizeof(std::int16_t));
} else if (std::strcmp(request.member, "hybrid_synthesis_values.npy") == 0) {
std::vector<float> values(count);
std::memcpy(values.data(), ordered.data(), count * sizeof(float));
out->hybrid_values.assign(values.begin(), values.end());
} else if (std::strcmp(request.member, "qmf_synthesis_basis.npy") == 0) {
std::memcpy(out->qmf_basis.empty() ? (out->qmf_basis.resize(count), out->qmf_basis.data())
: out->qmf_basis.data(),
ordered.data(), count * sizeof(double));
out->qmf_basis.resize(count);
} else {
out->qmf_taps.resize(count);
std::memcpy(out->qmf_taps.data(), ordered.data(), count * sizeof(double));
}
}
out->hybrid_count = static_cast<std::uint32_t>(out->hybrid_values.size());
if (out->hybrid_count == 0u) {
return Status::fail(JOC_ERR_HRTF_FORMAT, stage::kRender,
"kernel tables contain no hybrid synthesis entries");
}
return Status::success();
}
} // namespace joc::hrtf
+66
View File
@@ -0,0 +1,66 @@
#pragma once
#include <cstdint>
#include <string>
#include <vector>
#include "foundation/status.h"
namespace joc::hrtf {
inline constexpr int kShTerms = 36;
inline constexpr int kEars = 2;
inline constexpr int kHybridBands = 77;
inline constexpr int kFormatVersion = 1;
inline constexpr const char* kMagic = "JOC-HRTF-CACHE";
inline constexpr const char* kCacheSchema = "joc-compiled-hrtf-v1";
struct Field {
std::vector<double> coefficients;
std::vector<double> delay_coefficients;
std::vector<double> delay_bounds;
std::vector<double> band_centers_hz;
double measurement_radius_m = 1.0;
long long order = 5;
std::string source_sha256;
std::string cache_key;
std::string payload_sha256;
std::string delay_source;
std::string compiler_version;
std::string phase_policy_version;
std::string sh_convention;
std::string filterbank_json;
std::string metadata_json;
// Compile-side metadata, needed to write the cache back out unchanged.
std::string source_display_name;
std::string fit_report_json;
double projection_ridge = 0.0;
double spherical_harmonic_ridge = 0.0;
};
Status load_jochrtf(const std::string& path, Field* out);
// Binaural filterbank kernels, as the reused kernel expects them (C order, the
// exact dtypes of the ABI parameters).
struct Kernels {
std::vector<double> qmf_analysis;
std::vector<double> hybrid_low;
std::vector<std::int16_t> hybrid_indices;
std::vector<double> hybrid_values;
std::vector<double> qmf_basis;
std::vector<double> qmf_taps;
std::uint32_t hybrid_count = 0;
};
// Loads a kernel-table archive. The file path is an override for verification;
// the shipped tables are embedded (see builtin_kernels) so no data file is needed.
// The Fortran-order index member is transposed into C order on purpose: the reused
// kernel indexes the hybrid synthesis table row-major.
Status load_kernels(const std::string& npz_path, Kernels* out);
// The public filterbank tables compiled into the library (identical values to the
// archive the file loader accepts; the unit test checks their hashes).
const Kernels& builtin_kernels();
} // namespace joc::hrtf
File diff suppressed because it is too large Load Diff
+610
View File
@@ -0,0 +1,610 @@
#include "hrtf/public_filterbank.h"
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <memory>
#include <string>
#include <utility>
#include <vector>
#include "foundation/fft.h"
#include "simd/simd.h"
namespace joc::hrtf {
namespace {
constexpr double kPi = 3.14159265358979323846;
constexpr std::size_t kQmfLength = dsp::kQmfFftSize;
constexpr int kQmfTaps = 10;
constexpr int kSynthesisRank = 4;
constexpr int kSynthesisTaps = 10;
// ----------------------------------------------------------- filterbank -----
// One shared forward plan for the 128-point QMF transform. The analysis bank runs
// it 2 * slots * channels times per chunk, so the twiddle recurrence is built once
// instead of being re-derived inside every butterfly.
const dsp::FftPlan& qmf_fft_plan() {
static const dsp::FftPlan plan(dsp::kQmfFftSize, false);
return plan;
}
// Public 64-band complex QMF analysis (public_filterbank.QmfAnalysis).
class QmfAnalysis {
public:
static_assert(static_cast<std::size_t>(kQmfBands) == simd::kQmfAnalysisBands,
"the dispatched accumulate is written for this band count");
QmfAnalysis(const Kernels& kernels, std::size_t channels)
: channels_(channels), coefficients_(kernels.qmf_analysis) {
history_.assign(9u * channels_ * kQmfBands, 0.0);
// The polyphase MAC consumes one coefficient per band, so the shipped
// [band][tap] layout makes its inner loop a stride-10 gather. Transposing
// once here turns that into a contiguous AXPY. The coefficient values and
// the accumulation order are untouched, so the sums are bit-identical.
coefficients_by_lag_.resize(static_cast<std::size_t>(kQmfTaps) * kQmfBands);
for (int band = 0; band < kQmfBands; ++band) {
for (int tap = 0; tap < kQmfTaps; ++tap) {
coefficients_by_lag_[static_cast<std::size_t>(tap) * kQmfBands +
static_cast<std::size_t>(band)] =
coefficients_[static_cast<std::size_t>(band) * kQmfTaps +
static_cast<std::size_t>(tap)];
}
}
premultiply_.resize(kQmfBands);
post_.resize(kQmfBands);
even_post_.resize(kQmfBands);
for (int band = 0; band < kQmfBands; ++band) {
const double phase = static_cast<double>(band);
premultiply_[static_cast<std::size_t>(band)] =
std::polar(1.0, -kPi * phase / 128.0);
post_[static_cast<std::size_t>(band)] =
std::polar(1.0, -3.0 * (phase + 0.5) * kPi / 128.0);
even_post_[static_cast<std::size_t>(band)] =
Complex(0.0, band % 2 == 0 ? 1.0 : -1.0);
}
}
void reset() { std::fill(history_.begin(), history_.end(), 0.0); }
// samples: [slots*64, channels]; output: [slots, channels, 64] complex.
void process(const std::vector<double>& samples, std::size_t slots,
std::vector<Complex>* output) {
const std::size_t joined_slots = 9u + slots;
const std::size_t history_size = 9u * channels_ * kQmfBands;
const std::size_t joined_size = joined_slots * channels_ * kQmfBands;
// The joined window is filled completely -- the history lands in its first
// 9 * channels * 64 entries and the new samples in the rest -- so it is a
// reusable scratch buffer rather than a fresh zero-filled allocation. The
// history tail is taken by index instead of from end(), because the buffer may
// be longer than the window this call uses.
if (joined_.size() < joined_size) {
joined_.resize(joined_size);
}
std::copy(history_.begin(), history_.end(), joined_.begin());
std::copy(samples.begin(), samples.begin() + static_cast<std::ptrdiff_t>(slots * channels_ * kQmfBands),
joined_.begin() + static_cast<std::ptrdiff_t>(history_size));
// The two polyphase accumulators are read before they are written, so their
// zero fill is load-bearing and stays; only the per-call allocation goes.
const std::size_t accumulator_size = slots * channels_ * kQmfBands;
if (even_.size() < accumulator_size) {
even_.resize(accumulator_size);
}
if (odd_.size() < accumulator_size) {
odd_.resize(accumulator_size);
}
std::fill(even_.begin(), even_.begin() + static_cast<std::ptrdiff_t>(accumulator_size), 0.0);
std::fill(odd_.begin(), odd_.begin() + static_cast<std::ptrdiff_t>(accumulator_size), 0.0);
// The ten lags are ten accumulate passes over the same 64 bands with one
// shared coefficient row; the bands are independent accumulations of a
// single product each, so they are what the dispatched kernel puts in its
// lanes, and every band keeps the caller's own multiply-then-add.
//
// Slots are processed in blocks, with the lag loop inside: one lag pass
// touches every source row once, so running the ten passes over the whole
// chunk re-reads the joined window ten times -- at 1536 slots that is
// hundreds of megabytes per chunk and the loop ends up bound by memory, not
// by arithmetic. A block's ten lag passes instead slide over a window of
// (block + 9) rows that stays in the second-level cache. Lags still run in
// ascending order inside a block, which is the order each output's sum is
// formed in, so nothing about the arithmetic changes.
constexpr std::size_t kSlotBlock = 32;
for (std::size_t first = 0u; first < slots; first += kSlotBlock) {
const std::size_t block = std::min(kSlotBlock, slots - first);
for (int lag = 0; lag < kQmfTaps; ++lag) {
std::vector<double>& target = (lag % 2 == 0) ? even_ : odd_;
const double* row =
coefficients_by_lag_.data() + static_cast<std::size_t>(lag) * kQmfBands;
const std::size_t source_slot = 9u - static_cast<std::size_t>(lag) + first;
simd::qmf_analysis_taps(
target.data() + first * channels_ * kQmfBands,
joined_.data() + source_slot * channels_ * kQmfBands, row,
block * channels_);
}
}
std::copy(joined_.begin() + static_cast<std::ptrdiff_t>(joined_size - history_size),
joined_.begin() + static_cast<std::ptrdiff_t>(joined_size), history_.begin());
// Every output element is assigned below, so the size is all that has to be
// established; a resize of an already correctly sized buffer touches nothing.
output->resize(slots * channels_ * kQmfBands);
std::array<Complex, dsp::kQmfFftSize> even_spectrum{};
std::array<Complex, dsp::kQmfFftSize> odd_spectrum{};
for (std::size_t slot = 0u; slot < slots; ++slot) {
for (std::size_t channel = 0u; channel < channels_; ++channel) {
const double* even_values = even_.data() + (slot * channels_ + channel) * kQmfBands;
const double* odd_values = odd_.data() + (slot * channels_ + channel) * kQmfBands;
transform(even_values, &even_spectrum);
transform(odd_values, &odd_spectrum);
Complex* destination =
output->data() + (slot * channels_ + channel) * kQmfBands;
for (int band = 0; band < kQmfBands; ++band) {
destination[band] = odd_spectrum[static_cast<std::size_t>(band)] +
even_spectrum[static_cast<std::size_t>(band)] *
even_post_[static_cast<std::size_t>(band)];
}
}
}
}
private:
void transform(const double* values, std::array<Complex, dsp::kQmfFftSize>* spectrum) {
for (int band = 0; band < kQmfBands; ++band) {
(*spectrum)[static_cast<std::size_t>(band)] =
Complex(values[band], 0.0) * premultiply_[static_cast<std::size_t>(band)];
}
for (int index = kQmfBands; index < dsp::kQmfFftSize; ++index) {
(*spectrum)[static_cast<std::size_t>(index)] = Complex(0.0, 0.0);
}
dsp::fft_radix2(spectrum, qmf_fft_plan());
for (int band = 0; band < kQmfBands; ++band) {
(*spectrum)[static_cast<std::size_t>(band)] *= post_[static_cast<std::size_t>(band)];
}
}
std::size_t channels_;
std::vector<double> coefficients_; // [64][10]
std::vector<double> coefficients_by_lag_; // [10][64], the same values transposed
std::vector<double> history_; // [9][channels][64]
std::vector<double> joined_; // scratch, [9 + slots][channels][64]
std::vector<double> even_; // scratch, [slots][channels][64], zeroed per call
std::vector<double> odd_; // scratch, [slots][channels][64], zeroed per call
std::vector<Complex> premultiply_;
std::vector<Complex> post_;
std::vector<Complex> even_post_;
};
// Sparse 64-QMF to 77-hybrid analysis (public_filterbank.HybridAnalysis).
class HybridAnalysis {
public:
HybridAnalysis(const Kernels& kernels, std::size_t channels)
: channels_(channels), low_kernel_(kernels.hybrid_low) {
history_.assign(12u * channels_ * 3u * 2u, 0.0);
high_history_.assign(6u * channels_ * 61u, Complex(0.0, 0.0));
// The dispatched join walks one term at a time and adds its 32 weights to
// 32 outputs, so the shipped [tap][band][component] table is regrouped to
// the term order the caller accumulates in. Same weights, same order.
const std::size_t outputs = simd::kHybridOutputs;
low_by_term_.resize(simd::kHybridTerms * outputs);
for (int lag = 0; lag < 13; ++lag) {
for (int point = 0; point < 3; ++point) {
for (int input = 0; input < 2; ++input) {
const std::size_t term =
(static_cast<std::size_t>(lag) * 3u + static_cast<std::size_t>(point)) * 2u +
static_cast<std::size_t>(input);
const std::size_t source = (static_cast<std::size_t>(point) * 2u +
static_cast<std::size_t>(input)) * 13u +
static_cast<std::size_t>(lag);
for (std::size_t output = 0u; output < outputs; ++output) {
low_by_term_[term * outputs + output] =
low_kernel_[source * outputs + output];
}
}
}
}
low_values_.resize(simd::kHybridJoinBlock * simd::kHybridTerms);
low_out_.resize(simd::kHybridJoinBlock * outputs);
}
void reset() {
std::fill(history_.begin(), history_.end(), 0.0);
std::fill(high_history_.begin(), high_history_.end(), Complex(0.0, 0.0));
}
// qmf: [slots, channels, 64]; output: [slots, channels, 77] complex.
void process(const std::vector<Complex>& qmf, std::size_t slots,
std::vector<Complex>* output) {
const std::size_t joined_slots = 12u + slots;
const std::size_t history_size = 12u * channels_ * 6u;
const std::size_t joined_size = joined_slots * channels_ * 3u * 2u;
// Both the joined window and the pending high-band history are written in full
// before they are read, so they are reused scratch buffers; the history tail is
// taken by index because the buffer can be longer than this call's window.
if (joined_.size() < joined_size) {
joined_.resize(joined_size);
}
std::copy(history_.begin(), history_.end(), joined_.begin());
for (std::size_t slot = 0u; slot < slots; ++slot) {
for (std::size_t channel = 0u; channel < channels_; ++channel) {
const Complex* source = qmf.data() + (slot * channels_ + channel) * kQmfBands;
double* destination =
joined_.data() + ((12u + slot) * channels_ + channel) * 6u;
for (int band = 0; band < 3; ++band) {
destination[static_cast<std::size_t>(band) * 2u] = source[band].real();
destination[static_cast<std::size_t>(band) * 2u + 1u] = source[band].imag();
}
}
}
// The low bands are accumulated in a register block and written straight into
// the output, and the high bands are written by the pass below; between them
// every one of the 77 bands is assigned, so only the size has to be set.
output->resize(slots * channels_ * kHybridBands);
// The thirteen taps are summed in a per-output register block and the low
// bands are written straight into the output. Keeping a separate low plane
// and then copying it into the output re-streams tens of megabytes per chunk
// for nothing, and only the first kHybridLow bands are ever touched. The
// join itself is dispatched (see src/simd/simd.h): the 32 outputs of a
// row are 32 independent accumulations over the same 78 terms, which is what
// shares a vector. Every lane keeps the caller's term order -- lag, then
// point, then input -- and its two roundings, and skips exactly the terms
// this loop skips. Rows are staged in blocks so the gathered values do not
// spill out of the first-level cache.
const std::size_t hybrid_rows = slots * channels_;
const std::size_t block = simd::kHybridJoinBlock;
const std::size_t terms = simd::kHybridTerms;
for (std::size_t first = 0u; first < hybrid_rows; first += block) {
const std::size_t count = std::min(block, hybrid_rows - first);
for (std::size_t index = 0u; index < count; ++index) {
const std::size_t row = first + index;
const std::size_t slot = row / channels_;
const std::size_t channel = row % channels_;
double* staged = low_values_.data() + index * terms;
for (int lag = 0; lag < 13; ++lag) {
const std::size_t source_slot = 12u - static_cast<std::size_t>(lag) + slot;
const double* source =
joined_.data() + (source_slot * channels_ + channel) * 6u;
for (int point = 0; point < 3; ++point) {
for (int input = 0; input < 2; ++input) {
staged[(static_cast<std::size_t>(lag) * 3u +
static_cast<std::size_t>(point)) * 2u +
static_cast<std::size_t>(input)] =
source[static_cast<std::size_t>(point) * 2u +
static_cast<std::size_t>(input)];
}
}
}
}
simd::hybrid_low_join(low_values_.data(), low_by_term_.data(),
low_out_.data(), count);
for (std::size_t index = 0u; index < count; ++index) {
Complex* destination = output->data() + (first + index) * kHybridBands;
const double* values = low_out_.data() + index * simd::kHybridOutputs;
for (int band = 0; band < kHybridLow; ++band) {
destination[band] = Complex(values[static_cast<std::size_t>(band) * 2u],
values[static_cast<std::size_t>(band) * 2u + 1u]);
}
}
}
std::copy(joined_.begin() + static_cast<std::ptrdiff_t>(joined_size - history_size),
joined_.begin() + static_cast<std::ptrdiff_t>(joined_size), history_.begin());
// The high bands pass through unchanged but delayed by the six slots of
// history the reference concatenates in front of them. Only the last six
// entries of that concatenation survive into high_history_, so a six-entry
// register replaces the (6 + slots) plane and its full copy. Note the
// output reads the concatenation at index `slot`, not `6 + slot`, so the
// first six output slots come from the history: that offset is part of the
// current output and is preserved verbatim.
for (std::size_t slot = 0u; slot < slots; ++slot) {
for (std::size_t channel = 0u; channel < channels_; ++channel) {
Complex* destination = output->data() +
(slot * channels_ + channel) * kHybridBands + kHybridLow;
if (slot < 6u) {
const Complex* source =
high_history_.data() + (slot * channels_ + channel) * 61u;
for (int band = 0; band < 61; ++band) {
destination[band] = source[band];
}
} else {
const Complex* source =
qmf.data() + ((slot - 6u) * channels_ + channel) * kQmfBands;
for (int band = 3; band < kQmfBands; ++band) {
destination[static_cast<std::size_t>(band - 3)] = source[band];
}
}
}
}
// Every entry of the pending high-band history is written here, so it is a
// reusable scratch buffer; the copy into the live history is kept as it was.
if (next_high_history_.size() < 6u * channels_ * 61u) {
next_high_history_.resize(6u * channels_ * 61u);
}
for (std::size_t entry = 0u; entry < 6u; ++entry) {
const std::size_t combined = slots + entry;
for (std::size_t channel = 0u; channel < channels_; ++channel) {
Complex* destination =
next_high_history_.data() + (entry * channels_ + channel) * 61u;
if (combined < 6u) {
const Complex* source =
high_history_.data() + (combined * channels_ + channel) * 61u;
for (int band = 0; band < 61; ++band) {
destination[band] = source[band];
}
} else {
const Complex* source =
qmf.data() + ((combined - 6u) * channels_ + channel) * kQmfBands;
for (int band = 3; band < kQmfBands; ++band) {
destination[static_cast<std::size_t>(band - 3)] = source[band];
}
}
}
}
std::copy(next_high_history_.begin(), next_high_history_.end(), high_history_.begin());
}
private:
std::size_t channels_;
std::vector<double> low_kernel_; // [3][2][13][16][2]
std::vector<double> low_by_term_; // [78][32], the same weights in the caller's term order
std::vector<double> low_values_; // scratch, [block][78]
std::vector<double> low_out_; // scratch, [block][32]
std::vector<double> history_; // [12][channels][3][2]
std::vector<Complex> high_history_; // [6][channels][61]
std::vector<double> joined_; // scratch, [12 + slots][channels][3][2]
std::vector<Complex> next_high_history_; // scratch, [6][channels][61]
};
// Instantaneous sparse 77-hybrid to 64-QMF synthesis map.
class HybridSynthesis {
public:
explicit HybridSynthesis(const Kernels& kernels) {
const std::size_t rows = kernels.hybrid_indices.size() / 4u;
mapping_.reserve(rows);
for (std::size_t index = 0u; index < rows; ++index) {
Entry entry;
for (int field = 0; field < 4; ++field) {
entry.index[static_cast<std::size_t>(field)] =
kernels.hybrid_indices[index * 4u + static_cast<std::size_t>(field)];
}
entry.gain = kernels.hybrid_values[index];
mapping_.push_back(entry);
}
}
// hybrid: [slots, channels, 77]; output: [slots, channels, 64] complex.
// The sparse map moves a real or imaginary part of one band into a real or
// imaginary part of another, so the two components are accumulated apart.
void process(const std::vector<Complex>& hybrid, std::size_t slots, std::size_t channels,
std::vector<Complex>* output) const {
const std::size_t rows = slots * channels;
std::vector<double> real(rows * kQmfBands, 0.0);
std::vector<double> imaginary(rows * kQmfBands, 0.0);
for (std::size_t slot = 0u; slot < slots; ++slot) {
for (std::size_t channel = 0u; channel < channels; ++channel) {
const std::size_t row = slot * channels + channel;
const Complex* source = hybrid.data() + row * kHybridBands;
for (const Entry& entry : mapping_) {
const double value = entry.index[1] == 0u ? source[entry.index[0]].real()
: source[entry.index[0]].imag();
if (value == 0.0) {
continue;
}
double* destination =
(entry.index[3] == 0u ? real.data() : imaginary.data()) + row * kQmfBands;
destination[entry.index[2]] += value * entry.gain;
}
}
}
// Every output element is assigned from the two accumulators below, so the
// zero fill that `assign` performed was dead; only the size is needed.
output->resize(rows * kQmfBands);
for (std::size_t index = 0u; index < output->size(); ++index) {
(*output)[index] = Complex(real[index], imaginary[index]);
}
}
private:
struct Entry {
std::size_t index[4] = {0u, 0u, 0u, 0u};
double gain = 0.0;
};
std::vector<Entry> mapping_;
};
// Rank-4 64-band synthesis.
class QmfSynthesis {
public:
QmfSynthesis(const Kernels& kernels, std::size_t channels)
: channels_(channels), basis_(kernels.qmf_basis), taps_(kernels.qmf_taps) {
history_.assign(9u * channels_ * kQmfBands * kSynthesisRank, 0.0);
// The dispatched basis kernel reads the four ranks of one (band, tap) as
// one vector, so the shipped [band][rank][tap] table is reordered once
// here. The weights are the same doubles, only their order differs.
const std::size_t bands = static_cast<std::size_t>(kQmfBands);
const std::size_t ranks = static_cast<std::size_t>(kSynthesisRank);
const std::size_t taps = dsp::kQmfFftSize;
basis_by_tap_.resize(bands * taps * ranks);
for (std::size_t band = 0u; band < bands; ++band) {
for (std::size_t tap = 0u; tap < taps; ++tap) {
for (std::size_t rank = 0u; rank < ranks; ++rank) {
basis_by_tap_[(band * taps + tap) * ranks + rank] =
basis_[(band * ranks + rank) * taps + tap];
}
}
}
}
void reset() { std::fill(history_.begin(), history_.end(), 0.0); }
// qmf: [slots, channels, 64]; output: [slots*64, channels] real.
void process(const std::vector<Complex>& qmf, std::size_t slots, std::vector<double>* output) {
const std::size_t rows = slots * channels_;
// [row][band][component] staging for the basis application. Both staging
// planes and the joined window are reusable scratch: every element of each is
// written before it is read, so the buffers are sized once and kept instead of
// being allocated and zero-filled on every call.
const std::size_t flat_size = rows * dsp::kQmfFftSize;
if (flat_.size() < flat_size) {
flat_.resize(flat_size);
}
for (std::size_t row = 0u; row < rows; ++row) {
for (int band = 0; band < kQmfBands; ++band) {
flat_[row * dsp::kQmfFftSize + static_cast<std::size_t>(band) * 2u] =
qmf[row * kQmfBands + static_cast<std::size_t>(band)].real();
flat_[row * dsp::kQmfFftSize + static_cast<std::size_t>(band) * 2u + 1u] =
qmf[row * kQmfBands + static_cast<std::size_t>(band)].imag();
}
}
// The sums are written straight into the joined window: the destination index
// is known up front, the summation order is untouched, and the application
// itself is dispatched -- the four ranks of a band are four independent dot
// products over the same 128 values, so they share a vector while every lane
// keeps the tap order and the two roundings of `sum +=`.
const std::size_t history_size = 9u * channels_ * kQmfBands * kSynthesisRank;
const std::size_t joined_size = history_size + rows * kQmfBands * kSynthesisRank;
if (joined_.size() < joined_size) {
joined_.resize(joined_size);
}
std::copy(history_.begin(), history_.end(), joined_.begin());
simd::qmf_synthesis_basis(flat_.data(), basis_by_tap_.data(),
joined_.data() + history_size, rows);
std::copy(joined_.begin() + static_cast<std::ptrdiff_t>(joined_size - history_size),
joined_.begin() + static_cast<std::ptrdiff_t>(joined_size), history_.begin());
output->assign(rows * kQmfBands, 0.0);
for (int lag = 0; lag < kSynthesisTaps; ++lag) {
for (std::size_t slot = 0u; slot < slots; ++slot) {
const std::size_t source_slot = 9u - static_cast<std::size_t>(lag) + slot;
for (std::size_t channel = 0u; channel < channels_; ++channel) {
const double* source =
joined_.data() +
(source_slot * channels_ + channel) * kQmfBands * kSynthesisRank;
double* destination =
output->data() + (slot * channels_ + channel) * kQmfBands;
for (int band = 0; band < kQmfBands; ++band) {
double sum = 0.0;
for (int rank = 0; rank < kSynthesisRank; ++rank) {
sum += source[static_cast<std::size_t>(band) * kSynthesisRank +
static_cast<std::size_t>(rank)] *
taps_[(static_cast<std::size_t>(band) * kSynthesisTaps +
static_cast<std::size_t>(lag)) * kSynthesisRank +
static_cast<std::size_t>(rank)];
}
destination[band] += sum;
}
}
}
}
}
private:
std::size_t channels_;
std::vector<double> basis_; // [64][4][128]
std::vector<double> basis_by_tap_; // [64][128][4], the same weights transposed
std::vector<double> taps_; // [64][10][4]
std::vector<double> history_; // [9][channels][64][4]
std::vector<double> flat_; // scratch, [rows][128], fully written per call
std::vector<double> joined_; // scratch, [9 + slots][channels][64][4]
};
// public_filterbank.PublicAnalysis77.process: [N, channels] -> [N/64, channels, 77].
void analysis_77(const std::vector<double>& samples, std::size_t slots,
std::size_t channels, QmfAnalysis& qmf,
HybridAnalysis& hybrid_analysis, std::vector<Complex>* hybrid) {
std::vector<double> hops(slots * channels * kQmfHop, 0.0);
for (std::size_t slot = 0u; slot < slots; ++slot) {
for (std::size_t channel = 0u; channel < channels; ++channel) {
for (int index = 0; index < kQmfHop; ++index) {
hops[(slot * channels + channel) * kQmfHop + static_cast<std::size_t>(index)] =
samples[(slot * kQmfHop + static_cast<std::size_t>(index)) * channels + channel];
}
}
}
std::vector<Complex> qmf_bands;
qmf.process(hops, slots, &qmf_bands);
hybrid_analysis.process(qmf_bands, slots, hybrid);
}
// public_filterbank.PublicSynthesis77.process: [slots, channels, 77] -> [slots*64, channels].
void synthesis_77(const std::vector<Complex>& hybrid, std::size_t slots,
std::size_t channels, const HybridSynthesis& synthesis,
QmfSynthesis& qmf, std::vector<double>* time) {
std::vector<Complex> qmf_bands;
synthesis.process(hybrid, slots, channels, &qmf_bands);
std::vector<double> samples;
qmf.process(qmf_bands, slots, &samples);
// The reference transposes (slots, channels, 64) to sample-major output.
time->assign(samples.size(), 0.0);
for (std::size_t slot = 0u; slot < slots; ++slot) {
for (std::size_t channel = 0u; channel < channels; ++channel) {
for (int band = 0; band < kQmfBands; ++band) {
(*time)[(slot * kQmfHop + static_cast<std::size_t>(band)) * channels + channel] =
samples[(slot * channels + channel) * kQmfBands + static_cast<std::size_t>(band)];
}
}
}
}
} // namespace
struct PublicFilterbank::Impl {
Impl(const Kernels& kernels, std::size_t channels)
: channels(channels), qmf(kernels, channels), hybrid_analysis(kernels, channels),
hybrid_synthesis(kernels), qmf_synthesis(kernels, channels) {}
std::size_t channels;
QmfAnalysis qmf;
HybridAnalysis hybrid_analysis;
HybridSynthesis hybrid_synthesis;
QmfSynthesis qmf_synthesis;
};
PublicFilterbank::PublicFilterbank(const Kernels& kernels, std::size_t channels)
: impl_(std::make_unique<Impl>(kernels, channels)) {}
PublicFilterbank::~PublicFilterbank() = default;
void PublicFilterbank::reset() {
impl_->qmf.reset();
impl_->hybrid_analysis.reset();
impl_->qmf_synthesis.reset();
}
void PublicFilterbank::analyze_full_rate(const std::vector<double>& samples, std::size_t slots,
std::vector<Complex>* hybrid) {
analysis_77(samples, slots, impl_->channels, impl_->qmf, impl_->hybrid_analysis, hybrid);
}
void PublicFilterbank::synthesize_full_rate(const std::vector<Complex>& hybrid, std::size_t slots,
std::vector<double>* time) {
synthesis_77(hybrid, slots, impl_->channels, impl_->hybrid_synthesis, impl_->qmf_synthesis,
time);
}
void PublicFilterbank::analyze_qmf(const std::vector<double>& hops, std::size_t slots,
std::vector<Complex>* qmf) {
impl_->qmf.process(hops, slots, qmf);
}
void PublicFilterbank::analyze_hybrid(const std::vector<Complex>& qmf, std::size_t slots,
std::vector<Complex>* hybrid) {
impl_->hybrid_analysis.process(qmf, slots, hybrid);
}
void PublicFilterbank::synthesize_hybrid(const std::vector<Complex>& hybrid, std::size_t slots,
std::vector<Complex>* qmf) {
impl_->hybrid_synthesis.process(hybrid, slots, impl_->channels, qmf);
}
void PublicFilterbank::synthesize_qmf(const std::vector<Complex>& qmf, std::size_t slots,
std::vector<double>* time) {
impl_->qmf_synthesis.process(qmf, slots, time);
}
} // namespace joc::hrtf
+55
View File
@@ -0,0 +1,55 @@
#pragma once
#include <complex>
#include <cstddef>
#include <memory>
#include <vector>
#include "hrtf/jochrtf.h"
// Public 64-QMF / 77-hybrid filterbank, shared by the SOFA field compiler and the
// Rosella renderer (upstream public_filterbank.py and rosella_filterbank.py are
// the same bank). Everything is float64/complex128, as the reference computes it,
// and the stateful half-steps are exposed because Rosella drives them directly.
namespace joc::hrtf {
inline constexpr int kQmfBands = 64;
inline constexpr int kQmfHop = 64;
inline constexpr int kHybridLow = 16;
inline constexpr int kHybridBandCount = 77;
inline constexpr int kLatencySamples = 961;
using Complex = std::complex<double>;
class PublicFilterbank {
public:
PublicFilterbank(const Kernels& kernels, std::size_t channels);
~PublicFilterbank();
PublicFilterbank(const PublicFilterbank&) = delete;
PublicFilterbank& operator=(const PublicFilterbank&) = delete;
void reset();
// Full-rate [slots*64, channels] -> hybrid [slots, channels, 77].
void analyze_full_rate(const std::vector<double>& samples, std::size_t slots,
std::vector<Complex>* hybrid);
// Hybrid [slots, channels, 77] -> full-rate [slots*64, channels].
void synthesize_full_rate(const std::vector<Complex>& hybrid, std::size_t slots,
std::vector<double>* time);
// The stateful half-steps, in the order the reference runs them.
void analyze_qmf(const std::vector<double>& hops, std::size_t slots,
std::vector<Complex>* qmf);
void analyze_hybrid(const std::vector<Complex>& qmf, std::size_t slots,
std::vector<Complex>* hybrid);
void synthesize_hybrid(const std::vector<Complex>& hybrid, std::size_t slots,
std::vector<Complex>* qmf);
void synthesize_qmf(const std::vector<Complex>& qmf, std::size_t slots,
std::vector<double>* time);
private:
struct Impl;
std::unique_ptr<Impl> impl_;
};
} // namespace joc::hrtf
+535
View File
@@ -0,0 +1,535 @@
#include "hrtf/rosella_model.h"
#include <algorithm>
#include <cmath>
#include <cstring>
#include <fstream>
#include <string>
#include <utility>
#include <vector>
#include "foundation/fs_utf8.h"
#include "foundation/mini_json.h"
#include "foundation/sha256.h"
namespace joc::hrtf {
namespace {
// The model's fixed-point lane scale: every stored value is a Q15 integer.
constexpr float kQ15 = 1.0f / 32768.0f;
Status model_fail(joc_error code, const std::string& message) {
return Status::fail(code, stage::kRender, message);
}
float q15(std::int32_t value) { return static_cast<float>(value) * kQ15; }
float q15_exp(std::int32_t value, int exponent) {
return q15(value) * static_cast<float>(std::ldexp(1.0, exponent));
}
std::uint16_t low16(std::int32_t value) {
return static_cast<std::uint16_t>(static_cast<std::uint32_t>(value) & 0xFFFFu);
}
std::string trim(const std::string& text) {
const std::size_t begin = text.find_first_not_of(" \t\r\n");
const std::size_t end = text.find_last_not_of(" \t\r\n");
return begin == std::string::npos ? std::string() : text.substr(begin, end - begin + 1u);
}
// The lane array is read straight out of the JSON text: it is one flat list of
// integers, and building a 15691-node DOM for it would only cost time.
bool parse_int_array(const std::string& raw, std::vector<std::int32_t>* out, std::string* error) {
out->clear();
const char* cursor = raw.c_str();
const char* end = cursor + raw.size();
while (cursor < end && *cursor != '[') {
++cursor;
}
if (cursor == end) {
*error = "rosella_coefficients must be a JSON array";
return false;
}
++cursor;
while (cursor < end) {
while (cursor < end && (*cursor == ' ' || *cursor == '\t' || *cursor == '\r' ||
*cursor == '\n' || *cursor == ',')) {
++cursor;
}
if (cursor >= end) {
break;
}
if (*cursor == ']') {
return true;
}
const bool negative = *cursor == '-';
if (negative) {
++cursor;
}
if (cursor >= end || *cursor < '0' || *cursor > '9') {
*error = "rosella_coefficients contains a non-integer value";
return false;
}
long long value = 0;
while (cursor < end && *cursor >= '0' && *cursor <= '9') {
value = value * 10 + (*cursor - '0');
if (value > (1ll << 40)) {
*error = "rosella_coefficients value is out of range";
return false;
}
++cursor;
}
// A fractional part or an exponent means the value is not an exact integer.
if (cursor < end && (*cursor == '.' || *cursor == 'e' || *cursor == 'E')) {
*error = "rosella_coefficients contains a non-integer value";
return false;
}
if (negative) {
value = -value;
}
if (value < -(1ll << 31) || value > (1ll << 31) - 1) {
*error = "rosella_coefficients value is outside signed int32";
return false;
}
out->push_back(static_cast<std::int32_t>(value));
}
*error = "rosella_coefficients array is truncated";
return false;
}
struct RpHeader {
std::uint16_t stored_checksum = 0;
std::uint16_t computed_checksum = 0;
bool checksum_valid = false;
bool table_a_present = false;
bool table_b_present = false;
bool table_c_present = false;
int table_a_dimension = 0;
int table_a_option = 0;
int table_a_extra = 0;
int table_b_dimension = 0;
int table_b_extra = 0;
int table_b_groups = 0;
int table_c_dimension = 0;
std::size_t active_lanes = 0;
};
Status inspect_rp(const std::vector<std::int32_t>& lanes, RpHeader* out) {
if (lanes.size() < 5u) {
return model_fail(JOC_ERR_HRTF_FORMAT, "Rosella rp must contain whole int32 lanes");
}
if (low16(lanes[0]) != 0x7072u) {
return model_fail(JOC_ERR_HRTF_FORMAT, "bad Rosella rp magic");
}
out->stored_checksum = low16(lanes[1]);
out->table_a_present = low16(lanes[2]) != 0u;
out->table_b_present = low16(lanes[3]) != 0u;
out->table_c_present = low16(lanes[4]) != 0u;
std::size_t index = 5u;
if (out->table_a_present) {
out->table_a_dimension = low16(lanes[index]);
out->table_a_option = low16(lanes[index + 1u]);
out->table_a_extra = low16(lanes[index + 2u]);
index += 5u;
} else {
out->table_a_dimension = 77;
}
if (out->table_b_present) {
if (!out->table_a_present) {
return model_fail(JOC_ERR_HRTF_FORMAT,
"Rosella rp table B cannot be present without table A");
}
out->table_b_dimension = low16(lanes[index]);
out->table_b_extra = low16(lanes[index + 1u]);
out->table_b_groups = low16(lanes[index + 2u]);
index += 3u;
}
if (out->table_c_present) {
out->table_c_dimension = low16(lanes[index]);
index += 1u;
}
const long long payload_words =
static_cast<long long>(index) - 2 +
(out->table_b_present ? (out->table_b_dimension + 380 * out->table_b_groups +
out->table_b_extra + 79)
: 0) +
(out->table_a_present ? (171 * out->table_a_extra + 79 +
2 * (out->table_a_option + 14 * out->table_a_dimension))
: 0) +
11 + (out->table_c_present ? (314 * out->table_c_dimension + 1) : 0);
if (payload_words < 0) {
return model_fail(JOC_ERR_HRTF_FORMAT, "malformed Rosella rp header");
}
out->active_lanes = static_cast<std::size_t>(2 + payload_words);
if (lanes.size() < out->active_lanes) {
return model_fail(JOC_ERR_HRTF_FORMAT, "Rosella rp is truncated");
}
std::uint32_t computed = 0xA569u;
for (std::size_t lane = 2u; lane < out->active_lanes; ++lane) {
computed ^= low16(lanes[lane]);
}
out->computed_checksum = static_cast<std::uint16_t>(computed & 0xFFFFu);
out->checksum_valid = out->computed_checksum == out->stored_checksum;
return Status::success();
}
// _unpack_field: the serialized 154-per-direction field lanes to the padded grid.
void unpack_field(const std::int32_t* serialized, int directions, int exponent,
std::vector<float>* padded) {
padded->assign(static_cast<std::size_t>(160 * directions), 0.0f);
const int stride8 = 8 * directions;
const int stride2 = 2 * directions;
for (int source = 0; source < 154 * directions; ++source) {
const int group4 = (source % stride8) / stride2;
const int destination = (group4 & 3) + 4 * (source % stride2 +
2 * directions * (source / stride8 +
(group4 >> 2)));
(*padded)[static_cast<std::size_t>(destination)] =
q15_exp(serialized[source], exponent);
}
}
// _unpack_table_a_grid: the serialized table-A rows to the padded lane grid.
void unpack_table_a_grid(const std::int32_t* serialized, int dimension, int serialized_rows,
int padded_rows, int lane_group, std::vector<float>* padded) {
padded->assign(static_cast<std::size_t>(padded_rows) * static_cast<std::size_t>(dimension),
0.0f);
const int group_width = lane_group * 4;
for (int source = 0; source < serialized_rows * dimension; ++source) {
const int remainder = source % group_width;
const int destination = (remainder / lane_group) +
4 * (remainder % lane_group +
group_width / 4 * (source / group_width));
(*padded)[static_cast<std::size_t>(destination)] = q15(serialized[source]);
}
}
void unpack_table_a_extra(const std::int32_t* serialized, std::vector<float>* padded) {
padded->assign(160u, 0.0f);
for (int source = 0; source < 154; ++source) {
const int remainder = source & 7;
const int destination = (remainder >> 1) + 4 * ((source & 1) + 2 * (source >> 3));
(*padded)[static_cast<std::size_t>(destination)] = q15(serialized[source]);
}
}
} // namespace
std::string RosellaModel::summary() const {
std::string name = capture.name.empty() ? std::string("unnamed") : capture.name;
return "Rosella personalized_headphone '" + name + "' (" +
(room_model.empty() ? std::string("unknown room") : room_model) + "), " +
std::to_string(table_a_dimension) + " HQMF / 77 hybrid @ " +
std::to_string(sample_rate) + " Hz";
}
Status load_personalized_headphone(const std::string& path, RosellaModel* out) {
if (out == nullptr) {
return model_fail(JOC_ERR_INVALID_ARGUMENT, "null Rosella model destination");
}
if (!fs_utf8::exists(path)) {
return model_fail(JOC_ERR_HRTF_NOT_FOUND, "personalized headphone model not found: " + path);
}
std::ifstream stream = fs_utf8::open_input(path);
if (!stream.good()) {
return model_fail(JOC_ERR_IO, "cannot open " + path);
}
std::string text((std::istreambuf_iterator<char>(stream)), std::istreambuf_iterator<char>());
if (text.empty()) {
return model_fail(JOC_ERR_HRTF_FORMAT, "empty personalized headphone model: " + path);
}
// The checksum is taken over the coefficient lanes, exactly as upstream hashes
// the int32 image of the array.
const std::size_t first = text.find_first_not_of(" \t\r\n");
if (first == std::string::npos || text[first] != '{') {
return model_fail(JOC_ERR_NOT_SUPPORTED,
"raw rp models are not supported; use a .personalized_headphone JSON");
}
std::vector<json::Member> root;
std::string error;
if (!json::parse_object(text, &root, &error)) {
return model_fail(JOC_ERR_HRTF_FORMAT, "invalid personalized headphone JSON: " + error);
}
const json::Member* personalized = json::find(root, "personalized_hrtf");
if (personalized == nullptr) {
return model_fail(JOC_ERR_HRTF_FORMAT, "personalized_hrtf is missing");
}
std::vector<json::Member> inner;
if (!json::parse_object(personalized->raw, &inner, &error)) {
return model_fail(JOC_ERR_HRTF_FORMAT, "invalid personalized_hrtf object: " + error);
}
const json::Member* virtualizer = json::find(inner, "virtualizer_parameters");
if (virtualizer == nullptr) {
return model_fail(JOC_ERR_HRTF_FORMAT, "virtualizer_parameters is missing");
}
std::vector<json::Member> parameters;
if (!json::parse_object(virtualizer->raw, &parameters, &error)) {
return model_fail(JOC_ERR_HRTF_FORMAT, "invalid virtualizer_parameters: " + error);
}
const json::Member* coefficient_member = json::find(parameters, "rosella_coefficients");
if (coefficient_member == nullptr) {
return model_fail(JOC_ERR_HRTF_FORMAT, "rosella_coefficients is missing");
}
RosellaModel model;
model.source_path = path;
if (const json::Member* member = json::find(parameters, "rosella_coefficients_version")) {
json::as_string(*member, &model.coefficient_version);
}
if (const json::Member* member = json::find(parameters, "room_model")) {
json::as_string(*member, &model.room_model);
}
if (const json::Member* capture = json::find(inner, "phrtf_capture_metadata")) {
std::vector<json::Member> fields;
if (json::parse_object(capture->raw, &fields, &error)) {
const std::pair<const char*, std::string*> mapping[] = {
{"capture_submission_date", &model.capture.capture_submission_date},
{"capture_type", &model.capture.capture_type},
{"label", &model.capture.label},
{"name", &model.capture.name},
{"phrtf_algorithm_version", &model.capture.algorithm_version},
{"phrtf_creation_date", &model.capture.creation_date},
{"uuid", &model.capture.uuid},
{"version", &model.capture.version},
};
for (const auto& entry : mapping) {
if (const json::Member* member = json::find(fields, entry.first)) {
json::as_string(*member, entry.second);
}
}
}
}
std::vector<std::int32_t> lanes;
if (!parse_int_array(coefficient_member->raw, &lanes, &error)) {
return model_fail(JOC_ERR_HRTF_FORMAT, error);
}
{
crypto::Sha256 hash;
hash.update(lanes.data(), lanes.size() * sizeof(std::int32_t));
model.coefficient_sha256 = hash.finish_hex();
}
RpHeader header;
Status status = inspect_rp(lanes, &header);
if (!status.ok()) {
return status;
}
if (!header.checksum_valid || header.active_lanes != lanes.size()) {
return model_fail(JOC_ERR_HRTF_FORMAT,
"invalid or non-active Rosella rp coefficient sequence");
}
if (!header.table_a_present || !header.table_b_present || header.table_c_present) {
return model_fail(JOC_ERR_NOT_SUPPORTED,
"the renderer requires table A+B and no table C");
}
if (header.table_a_dimension != 64 || header.table_a_option != 3) {
return model_fail(JOC_ERR_NOT_SUPPORTED,
"the renderer requires the observed 64-channel HQMF layout");
}
if (header.table_b_dimension != 20 || header.table_b_groups != 36) {
return model_fail(JOC_ERR_NOT_SUPPORTED,
"the renderer requires 20 hybrid groups and 36 direction terms");
}
const std::int32_t* values = lanes.data();
const std::size_t total = lanes.size();
std::size_t position = 13u;
const int extra = header.table_a_extra;
model.table_a_dimension = header.table_a_dimension;
model.table_a_option = header.table_a_option;
model.table_a_extra = extra;
model.table_a_header_field = low16(values[8]);
model.table_a_header_25 = low16(values[9]);
model.table_a_control = low16(values[position]);
model.field_exponent = values[position];
position += 1u;
const int option_count = header.table_a_option;
model.table_a_option_ids.resize(static_cast<std::size_t>(option_count));
for (int index = 0; index < option_count; ++index) {
model.table_a_option_ids[static_cast<std::size_t>(index)] =
low16(values[position + static_cast<std::size_t>(index)]);
}
position += static_cast<std::size_t>(option_count);
model.table_a_option_values.resize(static_cast<std::size_t>(option_count));
for (int index = 0; index < option_count; ++index) {
model.table_a_option_values[static_cast<std::size_t>(index)] =
q15(values[position + static_cast<std::size_t>(index)]);
}
position += static_cast<std::size_t>(option_count);
model.table_a_scalar = q15(values[position]);
position += 1u;
const int dimension = header.table_a_dimension;
unpack_table_a_grid(values + position, dimension, 16, 20, 16,
&model.table_a_filter_16x64_padded);
position += static_cast<std::size_t>(16 * dimension);
for (int index = 0; index < 4; ++index) {
model.table_a_four_integers[static_cast<std::size_t>(index)] =
low16(values[position + static_cast<std::size_t>(index)]);
}
position += 4u;
model.table_a_integer = low16(values[position]);
position += 1u;
unpack_table_a_grid(values + position, dimension, 8, 10, 8,
&model.table_a_filter_8x64_padded);
position += static_cast<std::size_t>(8 * dimension);
model.table_a_vector16.resize(16u);
for (int index = 0; index < 16; ++index) {
model.table_a_vector16[static_cast<std::size_t>(index)] =
q15(values[position + static_cast<std::size_t>(index)]);
}
position += 16u;
unpack_table_a_grid(values + position, dimension, 4, 5, 4,
&model.table_a_filter_4x64_padded);
position += static_cast<std::size_t>(4 * dimension);
model.table_a_extra_indices.resize(static_cast<std::size_t>(extra));
for (int index = 0; index < extra; ++index) {
model.table_a_extra_indices[static_cast<std::size_t>(index)] =
low16(values[position + static_cast<std::size_t>(index)]);
}
position += static_cast<std::size_t>(extra);
model.table_a_extra_fields_padded.assign(static_cast<std::size_t>(extra) * 160u, 0.0f);
std::vector<float> unpacked;
for (int index = 0; index < extra; ++index) {
unpack_table_a_extra(values + position, &unpacked);
std::copy(unpacked.begin(), unpacked.end(),
model.table_a_extra_fields_padded.begin() + static_cast<std::ptrdiff_t>(index) * 160);
position += 154u;
}
model.table_a_extra_vectors.assign(static_cast<std::size_t>(extra) * 16u, 0.0f);
for (int index = 0; index < extra; ++index) {
for (int lane = 0; lane < 16; ++lane) {
model.table_a_extra_vectors[static_cast<std::size_t>(index) * 16u +
static_cast<std::size_t>(lane)] =
q15(values[position + static_cast<std::size_t>(lane)]);
}
position += 16u;
}
const std::size_t table_b_start = position;
if (table_b_start != 13u + 1821u + static_cast<std::size_t>(171 * extra)) {
return model_fail(JOC_ERR_HRTF_FORMAT, "Rosella table-A parser lost its place");
}
model.sample_rate = 2 * low16(values[position]);
position += 1u;
model.matrix_exponent = values[position];
position += 1u;
const std::size_t matrix_count = 36u * 36u;
model.matrix_left.resize(matrix_count);
model.matrix_right.resize(matrix_count);
for (std::size_t index = 0; index < matrix_count; ++index) {
model.matrix_left[index] = q15_exp(values[position + index], model.matrix_exponent);
}
position += matrix_count;
for (std::size_t index = 0; index < matrix_count; ++index) {
model.matrix_right[index] = q15_exp(values[position + index], model.matrix_exponent);
}
position += matrix_count;
model.vector_left.resize(36u);
model.vector_right.resize(36u);
for (int index = 0; index < 36; ++index) {
model.vector_left[static_cast<std::size_t>(index)] =
q15_exp(values[position + static_cast<std::size_t>(index)], model.matrix_exponent);
}
position += 36u;
for (int index = 0; index < 36; ++index) {
model.vector_right[static_cast<std::size_t>(index)] =
q15_exp(values[position + static_cast<std::size_t>(index)], model.matrix_exponent);
}
position += 36u;
const std::size_t serialized_count = 154u * 36u;
unpack_field(values + position, 36, model.field_exponent, &model.field_left_padded);
bool odd_zero = true;
for (std::size_t index = 1u; index < serialized_count; index += 2u) {
const float value = q15_exp(values[position + index], model.field_exponent);
if (std::abs(value) > 1.0e-6f) {
odd_zero = false;
break;
}
}
model.field_left_odd_serialized_zero = odd_zero;
position += serialized_count;
unpack_field(values + position, 36, model.field_exponent, &model.field_right_padded);
position += serialized_count;
model.hybrid_flags.resize(20u);
int active_hybrid = 0;
for (int index = 0; index < 20; ++index) {
model.hybrid_flags[static_cast<std::size_t>(index)] =
low16(values[position + static_cast<std::size_t>(index)]);
if (model.hybrid_flags[static_cast<std::size_t>(index)] == 1) {
++active_hybrid;
}
}
position += 20u;
if (active_hybrid != header.table_b_extra) {
return model_fail(JOC_ERR_HRTF_FORMAT, "hybrid value count does not match the header");
}
model.hybrid_values.resize(static_cast<std::size_t>(active_hybrid));
for (int index = 0; index < active_hybrid; ++index) {
model.hybrid_values[static_cast<std::size_t>(index)] =
q15(values[position + static_cast<std::size_t>(index)]);
}
position += static_cast<std::size_t>(active_hybrid);
model.model_scalars.resize(5u);
for (int index = 0; index < 5; ++index) {
model.model_scalars[static_cast<std::size_t>(index)] =
q15(values[position + static_cast<std::size_t>(index)]);
}
position += 5u;
const std::size_t expected_tail =
table_b_start + static_cast<std::size_t>(header.table_b_dimension +
380 * header.table_b_groups +
header.table_b_extra + 79);
if (position != expected_tail) {
return model_fail(JOC_ERR_HRTF_FORMAT, "Rosella table-B parser lost its place");
}
model.header_float_scalars[0] = q15(values[position]);
model.header_float_scalars[1] = q15(values[position + 1u]) * 16.0f;
model.header_integer_fields[0] = values[position + 2u];
model.header_integer_fields[1] = low16(values[position + 3u]);
position += 4u;
for (int profile = 0; profile < 4; ++profile) {
RosellaDistanceProfile parsed;
for (int index = 0; index < 6; ++index) {
parsed.bounds[static_cast<std::size_t>(index)] =
q15(values[position + static_cast<std::size_t>(index)]);
}
position += 6u;
parsed.distance_scale_m =
q15_exp(values[position], values[position + 1u]);
position += 2u;
parsed.inverse_distance_per_m = q15(values[position]);
parsed.axis_scales_internal[0] = q15(values[position + 1u]);
parsed.axis_scales_internal[1] = q15(values[position + 2u]);
parsed.axis_scales_internal[2] = q15(values[position + 3u]);
parsed.minimum_normalized_radius = q15(values[position + 4u]);
position += 5u;
model.profiles[static_cast<std::size_t>(profile)] = parsed;
}
model.profile_tail.resize(8u);
for (int index = 0; index < 8; ++index) {
model.profile_tail[static_cast<std::size_t>(index)] =
q15(values[position + static_cast<std::size_t>(index)]);
}
position += 8u;
for (int index = 0; index < 3; ++index) {
model.post_fields[static_cast<std::size_t>(index)] =
values[position + static_cast<std::size_t>(index)];
}
position += 3u;
if (position != total) {
return model_fail(JOC_ERR_HRTF_FORMAT, "unparsed Rosella coefficient lanes");
}
if (model.sample_rate != 48000) {
return model_fail(JOC_ERR_HRTF_FORMAT,
"Rosella model sample rate must be 48000, got " +
std::to_string(model.sample_rate));
}
*out = std::move(model);
return Status::success();
}
} // namespace joc::hrtf
+88
View File
@@ -0,0 +1,88 @@
#pragma once
#include <array>
#include <cstdint>
#include <string>
#include <vector>
#include "foundation/status.h"
// Parser for the Dolby ".personalized_headphone" model (upstream rosella_model.py).
// The file is JSON whose virtualizer_parameters carry the raw "rp" coefficient
// lanes; everything the renderer needs is unpacked here, in the same float32
// arithmetic the reference uses, because those values are part of the model.
namespace joc::hrtf {
struct RosellaDistanceProfile {
std::array<float, 6> bounds{};
float distance_scale_m = 0.0f;
float inverse_distance_per_m = 0.0f;
std::array<float, 3> axis_scales_internal{};
float minimum_normalized_radius = 0.0f;
};
struct RosellaCaptureMetadata {
std::string capture_submission_date;
std::string capture_type;
std::string label;
std::string name;
std::string algorithm_version;
std::string creation_date;
std::string uuid;
std::string version;
};
struct RosellaModel {
std::string source_path;
std::string coefficient_sha256;
std::string coefficient_version;
std::string room_model;
RosellaCaptureMetadata capture;
int table_a_dimension = 0;
int table_a_option = 0;
int table_a_extra = 0;
int table_a_header_field = 0;
int table_a_header_25 = 0;
int table_a_control = 0;
std::vector<int> table_a_option_ids;
std::vector<float> table_a_option_values;
float table_a_scalar = 0.0f;
std::vector<float> table_a_filter_16x64_padded;
std::array<int, 4> table_a_four_integers{};
int table_a_integer = 0;
std::vector<float> table_a_filter_8x64_padded;
std::vector<float> table_a_vector16;
std::vector<float> table_a_filter_4x64_padded;
std::vector<int> table_a_extra_indices;
std::vector<float> table_a_extra_fields_padded;
std::vector<float> table_a_extra_vectors;
int sample_rate = 0;
int matrix_exponent = 0;
int field_exponent = 0;
std::vector<float> matrix_left;
std::vector<float> matrix_right;
std::vector<float> vector_left;
std::vector<float> vector_right;
std::vector<float> field_left_padded;
std::vector<float> field_right_padded;
bool field_left_odd_serialized_zero = false;
std::vector<int> hybrid_flags;
std::vector<float> hybrid_values;
std::vector<float> model_scalars;
std::array<float, 2> header_float_scalars{};
std::array<int, 2> header_integer_fields{};
std::array<RosellaDistanceProfile, 4> profiles{};
std::vector<float> profile_tail;
std::array<int, 3> post_fields{};
// One line for reports and logs: the capture name and room model are the
// model's own strings, followed by the table layout and sample rate, e.g.
// "Rosella personalized_headphone '<name>' (<room>), <N> HQMF / 77 hybrid @ <rate> Hz".
std::string summary() const;
};
Status load_personalized_headphone(const std::string& path, RosellaModel* out);
} // namespace joc::hrtf
File diff suppressed because it is too large Load Diff
+64
View File
@@ -0,0 +1,64 @@
#pragma once
#include <cstdint>
#include <memory>
#include <string>
#include <vector>
#include "foundation/status.h"
#include "hrtf/rosella_model.h"
#include "oamd/oamd_parser.h"
#include "timeline/position_timeline.h"
// Rosella ".personalized_headphone" binaural renderer (upstream rosella_core.py,
// rosella_direct.py, rosella_room.py and rosella_binaural_renderer.py). It takes
// the same pipeline slot as the SOFA runtime: sixteen object channels per frame in,
// interleaved stereo out, with the OAMD timeline driving the per-block parameters.
namespace joc::hrtf {
// rosella_direct.BINAURAL_PROFILE_NAMES.
enum class RosellaProfile : std::int32_t { Near = 1, Far = 2, Mid = 3 };
struct RosellaRenderOptions {
RosellaProfile profile = RosellaProfile::Mid;
std::int64_t object_delay_samples = 1473;
double tail_seconds = 5.0;
double output_gain = 1.0;
int chunk_frames = 64;
int room_impulse_slots = 4096;
};
class RosellaRuntime {
public:
RosellaRuntime();
~RosellaRuntime();
RosellaRuntime(const RosellaRuntime&) = delete;
RosellaRuntime& operator=(const RosellaRuntime&) = delete;
Status open(const RosellaModel& model, const RosellaRenderOptions& options);
// objects16_planar is channel-major: channel * 1536 + sample.
Status submit_frame(const float* objects16_planar, const oamd::OamdUpdate* update,
std::int64_t frame_index, std::int64_t outer_sample_offset,
std::int64_t object_delay_samples);
// Drains the flush tail: the pending partial chunk plus flush_samples of silence.
Status finish(std::uint32_t flush_samples, std::vector<double>* out);
std::uint32_t finish_capacity(double tail_seconds) const;
Status reset();
const std::vector<double>& output() const;
void take_output(std::vector<double>* out);
std::uint64_t input_samples() const;
std::uint64_t processed_input_samples() const;
std::uint64_t metadata_block_updates() const;
const timeline::OamdPositionTimeline& timeline() const;
private:
struct Impl;
std::unique_ptr<Impl> impl_;
};
} // namespace joc::hrtf
+255
View File
@@ -0,0 +1,255 @@
#include "hrtf/sofa.h"
#include <algorithm>
#include <cmath>
#include <cctype>
#include <cstdio>
#include <fstream>
#include <utility>
#include "foundation/fs_utf8.h"
#include "foundation/sha256.h"
#include "io/hdf5.h"
namespace joc::hrtf {
namespace {
Status sofa_fail(joc_error code, const std::string& message) {
return Status::fail(code, stage::kRender, message);
}
std::string format_number(double value) {
if (std::isfinite(value) && value == std::floor(value) && std::fabs(value) < 1.0e15) {
return std::to_string(static_cast<long long>(value));
}
char buffer[32];
std::snprintf(buffer, sizeof(buffer), "%.6g", value);
return std::string(buffer);
}
// Every array is checked against the element count the convention prescribes, so
// a file whose shape disagrees with its metadata is rejected instead of silently
// producing a shifted impulse response.
Status read_doubles(const io::Hdf5File& file, const std::string& path, std::uint64_t expected,
std::vector<double>* out) {
if (!file.has_dataset(path)) {
return sofa_fail(JOC_ERR_HRTF_FORMAT, "SOFA file has no " + path + " dataset");
}
const Status status = file.read_dataset_double(path, out);
if (!status.ok()) {
return sofa_fail(JOC_ERR_HRTF_FORMAT, "SOFA dataset " + path + ": " + status.message());
}
if (out->size() != expected) {
return sofa_fail(JOC_ERR_HRTF_FORMAT,
"SOFA dataset " + path + " holds " + std::to_string(out->size()) +
" values, expected " + std::to_string(expected));
}
return Status::success();
}
Status read_text(const io::Hdf5File& file, const std::string& name, bool required,
std::string* out) {
io::Hdf5Attribute attribute;
const Status status = file.attribute("", name, &attribute);
if (!status.ok()) {
if (required) {
return sofa_fail(JOC_ERR_HRTF_FORMAT,
"SOFA file has no root attribute " + name + ": " + status.message());
}
return Status::success();
}
*out = attribute.text;
return Status::success();
}
// SHA-256 of the whole file: the compiled-cache key is derived from it, so the
// digest is taken over the exact bytes the parse consumed.
std::string file_digest(const std::string& path) {
std::ifstream stream = fs_utf8::open_input(path);
if (!stream.good()) {
return std::string();
}
crypto::Sha256 hash;
std::vector<char> buffer(1u << 20);
while (stream.good()) {
stream.read(buffer.data(), static_cast<std::streamsize>(buffer.size()));
const std::streamsize count = stream.gcount();
if (count > 0) {
hash.update(buffer.data(), static_cast<std::size_t>(count));
}
}
return hash.finish_hex();
}
// The coordinate declaration of one dataset, when the file carries it.
void read_coordinates(const io::Hdf5File& file, const std::string& dataset,
SofaCoordinate* out) {
io::Hdf5Attribute attribute;
if (file.attribute(dataset, "Type", &attribute).ok()) {
out->type = attribute.text;
}
if (file.attribute(dataset, "Units", &attribute).ok()) {
out->units = attribute.text;
}
}
} // namespace
std::string SofaHrir::summary() const {
return "SOFA " + sofa_conventions + ", " + std::to_string(ir_count) + " IRs x " +
std::to_string(ir_length) + " taps @ " + format_number(sample_rate) + " Hz";
}
Status load_sofa(const std::string& path, SofaHrir* out) {
if (out == nullptr) {
return sofa_fail(JOC_ERR_INVALID_ARGUMENT, "null SOFA destination");
}
if (!fs_utf8::exists(path)) {
return sofa_fail(JOC_ERR_HRTF_NOT_FOUND, "SOFA file not found: " + path);
}
io::Hdf5File file;
Status status = file.open(path);
if (!status.ok()) {
return sofa_fail(JOC_ERR_HRTF_FORMAT, "SOFA file " + path + ": " + status.message());
}
SofaHrir sofa;
sofa.source_path = path;
sofa.source_sha256 = file_digest(path);
for (char& character : sofa.source_sha256) {
character = static_cast<char>(std::toupper(static_cast<unsigned char>(character)));
}
status = read_text(file, "Conventions", true, &sofa.conventions);
if (!status.ok()) {
return status;
}
status = read_text(file, "SOFAConventions", true, &sofa.sofa_conventions);
if (!status.ok()) {
return status;
}
if (sofa.conventions != "SOFA" || sofa.sofa_conventions != "SimpleFreeFieldHRIR") {
return sofa_fail(JOC_ERR_HRTF_UNSUPPORTED_CONVENTION,
"SOFA conventions " + sofa.conventions + "/" + sofa.sofa_conventions +
" are not SimpleFreeFieldHRIR");
}
status = read_text(file, "SOFAConventionsVersion", false, &sofa.convention_version);
if (!status.ok()) {
return status;
}
status = read_text(file, "Version", false, &sofa.version);
if (!status.ok()) {
return status;
}
status = read_text(file, "DataType", false, &sofa.data_type);
if (!status.ok()) {
return status;
}
status = read_text(file, "RoomType", false, &sofa.room_type);
if (!status.ok()) {
return status;
}
status = read_text(file, "Title", false, &sofa.title);
if (!status.ok()) {
return status;
}
status = read_text(file, "DatabaseName", false, &sofa.database_name);
if (!status.ok()) {
return status;
}
status = read_text(file, "ListenerShortName", false, &sofa.listener_short_name);
if (!status.ok()) {
return status;
}
status = read_text(file, "Comment", false, &sofa.comment);
if (!status.ok()) {
return status;
}
io::Hdf5DatasetInfo info;
status = file.dataset_info("Data.IR", &info);
if (!status.ok()) {
return sofa_fail(JOC_ERR_HRTF_FORMAT,
"SOFA file has no usable Data.IR dataset: " + status.message());
}
if (info.shape.size() != 3u || info.shape[1] != 2u || info.shape[0] == 0u ||
info.shape[2] == 0u) {
return sofa_fail(JOC_ERR_HRTF_FORMAT, "SOFA Data.IR is not shaped (M, 2, N)");
}
if (info.shape[0] > 0xFFFFFFFFull || info.shape[2] > 0xFFFFFFFFull) {
return sofa_fail(JOC_ERR_HRTF_FORMAT, "SOFA Data.IR is larger than this reader accepts");
}
sofa.ir_count = static_cast<std::uint32_t>(info.shape[0]);
sofa.ir_length = static_cast<std::uint32_t>(info.shape[2]);
const std::uint64_t taps = static_cast<std::uint64_t>(sofa.ir_count) * 2u * sofa.ir_length;
status = read_doubles(file, "Data.IR", taps, &sofa.ir);
if (!status.ok()) {
return status;
}
std::vector<double> scalar;
status = read_doubles(file, "Data.SamplingRate", 1u, &scalar);
if (!status.ok()) {
return status;
}
sofa.sample_rate = scalar[0];
io::Hdf5Attribute attribute;
if (file.attribute("Data.SamplingRate", "Units", &attribute).ok()) {
sofa.sampling_rate_units = attribute.text;
}
// Data.Delay is optional in the wild; absent means "no delay was measured".
if (file.has_dataset("Data.Delay")) {
std::vector<double> delay;
status = read_doubles(file, "Data.Delay", 2u, &delay);
if (!status.ok()) {
return status;
}
sofa.delay[0] = delay[0];
sofa.delay[1] = delay[1];
}
const std::uint64_t measurements = sofa.ir_count;
status = read_doubles(file, "SourcePosition", measurements * 3u, &sofa.source_position);
if (!status.ok()) {
return status;
}
read_coordinates(file, "SourcePosition", &sofa.source_position_coordinates);
std::vector<double> vector;
status = read_doubles(file, "ListenerPosition", 3u, &vector);
if (!status.ok()) {
return status;
}
std::copy(vector.begin(), vector.end(), sofa.listener_position);
read_coordinates(file, "ListenerPosition", &sofa.listener_position_coordinates);
status = read_doubles(file, "ListenerView", 3u, &vector);
if (!status.ok()) {
return status;
}
std::copy(vector.begin(), vector.end(), sofa.listener_view);
read_coordinates(file, "ListenerView", &sofa.listener_view_coordinates);
status = read_doubles(file, "ListenerUp", 3u, &vector);
if (!status.ok()) {
return status;
}
std::copy(vector.begin(), vector.end(), sofa.listener_up);
read_coordinates(file, "ListenerUp", &sofa.listener_up_coordinates);
status = read_doubles(file, "EmitterPosition", 3u, &vector);
if (!status.ok()) {
return status;
}
std::copy(vector.begin(), vector.end(), sofa.emitter_position);
read_coordinates(file, "EmitterPosition", &sofa.emitter_position_coordinates);
status = read_doubles(file, "ReceiverPosition", 6u, &vector);
if (!status.ok()) {
return status;
}
std::copy(vector.begin(), vector.end(), sofa.receiver_position);
read_coordinates(file, "ReceiverPosition", &sofa.receiver_position_coordinates);
*out = std::move(sofa);
return Status::success();
}
} // namespace joc::hrtf
+61
View File
@@ -0,0 +1,61 @@
#pragma once
#include <cstdint>
#include <string>
#include <vector>
#include "foundation/status.h"
namespace joc::hrtf {
// Coordinate declaration of one SOFA variable: Type ("spherical"/"cartesian") and
// Units. Empty when the file does not declare them (ListenerUp inherits).
struct SofaCoordinate {
std::string type;
std::string units;
};
// SOFA SimpleFreeFieldHRIR as this project consumes it: the impulse responses,
// the measurement geometry and the metadata needed to report what was loaded.
// All angles are degrees, all distances metres, exactly as the file stores them.
struct SofaHrir {
double sample_rate = 0.0;
std::uint32_t ir_count = 0; // M: number of measurements
std::uint32_t ir_length = 0; // N: taps per impulse response
std::vector<double> ir; // C order [M][2][N]
double delay[2] = {0.0, 0.0};
std::vector<double> source_position; // M*3
double listener_position[3] = {0.0, 0.0, 0.0};
double listener_view[3] = {1.0, 0.0, 0.0};
double listener_up[3] = {0.0, 0.0, 1.0};
double emitter_position[3] = {0.0, 0.0, 0.0};
double receiver_position[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
std::string conventions;
std::string sofa_conventions;
std::string convention_version;
std::string version;
std::string data_type;
std::string room_type;
std::string title;
std::string database_name;
std::string listener_short_name;
std::string comment;
std::string sampling_rate_units;
SofaCoordinate source_position_coordinates;
SofaCoordinate listener_position_coordinates;
SofaCoordinate listener_view_coordinates;
SofaCoordinate listener_up_coordinates;
SofaCoordinate emitter_position_coordinates;
SofaCoordinate receiver_position_coordinates;
// Identity of the file itself, needed for the compiled-cache key.
std::string source_path;
std::string source_sha256;
// One line for reports and logs:
// "SOFA SimpleFreeFieldHRIR, <M> IRs x <N> taps @ <rate> Hz".
std::string summary() const;
};
Status load_sofa(const std::string& path, SofaHrir* out);
} // namespace joc::hrtf
+237
View File
@@ -0,0 +1,237 @@
#include "hrtf/sofa_cache.h"
#include <algorithm>
#include <filesystem>
#include <list>
#include <mutex>
#include <utility>
#include <vector>
#include "foundation/fs_utf8.h"
#include "hrtf/sofa.h"
namespace joc::hrtf {
namespace {
namespace fs = std::filesystem;
// Small process-local cache: the reference keeps the last eight compiled fields.
constexpr std::size_t kMemoryCacheEntries = 8;
struct MemoryEntry {
std::string key;
Field field;
};
std::mutex& memory_mutex() {
static std::mutex mutex;
return mutex;
}
std::list<MemoryEntry>& memory_cache() {
static std::list<MemoryEntry> cache;
return cache;
}
bool memory_cache_get(const std::string& key, Field* out) {
std::lock_guard<std::mutex> lock(memory_mutex());
std::list<MemoryEntry>& cache = memory_cache();
for (auto entry = cache.begin(); entry != cache.end(); ++entry) {
if (entry->key == key) {
*out = entry->field;
cache.splice(cache.begin(), cache, entry);
return true;
}
}
return false;
}
void memory_cache_put(const std::string& key, const Field& field) {
std::lock_guard<std::mutex> lock(memory_mutex());
std::list<MemoryEntry>& cache = memory_cache();
for (auto entry = cache.begin(); entry != cache.end(); ++entry) {
if (entry->key == key) {
entry->field = field;
cache.splice(cache.begin(), cache, entry);
return;
}
}
cache.push_front(MemoryEntry{key, field});
while (cache.size() > kMemoryCacheEntries) {
cache.pop_back();
}
}
Status cache_fail(joc_error code, const std::string& message) {
return Status::fail(code, stage::kRender, message);
}
std::string upper(std::string text) {
std::transform(text.begin(), text.end(), text.begin(), [](unsigned char value) {
return static_cast<char>(std::toupper(value));
});
return text;
}
} // namespace
Status parse_cache_policy(const std::string& text, CachePolicy* out) {
if (out == nullptr) {
return cache_fail(JOC_ERR_INVALID_ARGUMENT, "null cache policy");
}
if (text == "none") {
*out = CachePolicy::None;
return Status::success();
}
if (text == "memory") {
*out = CachePolicy::Memory;
return Status::success();
}
if (text == "disk") {
*out = CachePolicy::Disk;
return Status::success();
}
return cache_fail(JOC_ERR_INVALID_CONFIG, "cache_policy must be none, memory, or disk");
}
Status validate_jochrtf(const std::string& path, const std::string& source_sha256,
const std::string& cache_key, Field* out) {
Field field;
const Status status = load_jochrtf(path, &field);
if (!status.ok()) {
return status;
}
if (!source_sha256.empty() && field.source_sha256 != upper(source_sha256)) {
return cache_fail(JOC_ERR_HRTF_HASH, "compiled HRTF source hash mismatch");
}
if (!cache_key.empty() && field.cache_key != upper(cache_key)) {
return cache_fail(JOC_ERR_HRTF_HASH, "compiled HRTF configuration hash mismatch");
}
if (out != nullptr) {
*out = std::move(field);
}
return Status::success();
}
Status save_jochrtf_atomic(const Field& field, const std::string& path) {
if (path.empty()) {
return cache_fail(JOC_ERR_INVALID_ARGUMENT, "empty compiled HRTF cache path");
}
const fs::path target = fs_utf8::to_path(path);
std::error_code error;
if (target.has_parent_path()) {
fs::create_directories(target.parent_path(), error);
if (error) {
return cache_fail(JOC_ERR_OUTPUT_OPEN,
"cannot create " + fs_utf8::from_path(target.parent_path()));
}
}
const std::string temporary = path + ".tmp";
Status status = write_jochrtf(field, temporary);
if (!status.ok()) {
return status;
}
// The rename is what makes a half-written cache impossible to observe.
fs::rename(fs_utf8::to_path(temporary), target, error);
if (error) {
fs::remove(fs_utf8::to_path(temporary), error);
return cache_fail(JOC_ERR_OUTPUT_WRITE, "cannot replace " + path);
}
return Status::success();
}
Status load_or_compile_sofa_field(const SofaFieldRequest& request, Field* out,
std::string* cache_path) {
if (out == nullptr) {
return cache_fail(JOC_ERR_INVALID_ARGUMENT, "null compiled HRTF destination");
}
if (request.sofa_path.empty()) {
return cache_fail(JOC_ERR_INVALID_ARGUMENT, "no SOFA path for the compiled HRTF field");
}
if (request.policy == CachePolicy::Disk && request.cache_dir.empty()) {
return cache_fail(JOC_ERR_INVALID_CONFIG, "the disk cache policy needs a cache directory");
}
SofaHrir sofa;
Status status = load_sofa(request.sofa_path, &sofa);
if (!status.ok()) {
return status;
}
CanonicalHrtf canonical;
status = canonicalize_sofa(sofa, &canonical);
if (!status.ok()) {
return status;
}
// The key depends on the shell that the radius selects, exactly as upstream.
double actual_radius = request.options.shell_radius_m;
(void)canonical_shell_indices(canonical, request.options.shell_radius_m, &actual_radius);
const std::string key = compiled_hrtf_cache_key(
canonical.source_sha256, canonical.sample_rate_hz, actual_radius, request.options.order,
request.options.projection_ridge, request.options.sh_ridge);
if (cache_path != nullptr) {
cache_path->clear();
}
std::string target;
if (request.policy == CachePolicy::Disk) {
target = request.cache_dir;
if (!target.empty() && target.back() != '/' && target.back() != '\\') {
target += "/";
}
target += cache_file_name(canonical.source_path.empty()
? std::string()
: canonical.source_path,
key);
if (fs_utf8::exists(target)) {
Field cached_field;
const Status cached =
validate_jochrtf(target, canonical.source_sha256, key, &cached_field);
if (cached.ok()) {
memory_cache_put(key, cached_field);
*out = std::move(cached_field);
if (cache_path != nullptr) {
*cache_path = target;
}
return Status::success();
}
}
}
Field field;
if (request.policy != CachePolicy::None && memory_cache_get(key, &field)) {
// A memory hit still materialises the disk cache the caller asked for.
if (request.policy == CachePolicy::Disk) {
status = save_jochrtf_atomic(field, target);
if (!status.ok()) {
return status;
}
if (cache_path != nullptr) {
*cache_path = target;
}
}
*out = std::move(field);
return Status::success();
}
status = compile_canonical_field(canonical, request.options, &field);
if (!status.ok()) {
return status;
}
if (field.cache_key != key) {
return cache_fail(JOC_ERR_INTERNAL, "internal compiled HRTF cache-key mismatch");
}
if (request.policy == CachePolicy::Disk) {
status = save_jochrtf_atomic(field, target);
if (!status.ok()) {
return status;
}
if (cache_path != nullptr) {
*cache_path = target;
}
}
if (request.policy != CachePolicy::None) {
memory_cache_put(key, field);
}
*out = std::move(field);
return Status::success();
}
} // namespace joc::hrtf
+38
View File
@@ -0,0 +1,38 @@
#pragma once
#include <string>
#include "foundation/status.h"
#include "hrtf/jochrtf.h"
#include "hrtf/sofa_field.h"
// Compiled-field cache: the .jochrtf is an internal artifact, so the caller only
// names the SOFA file and the policy. "memory" keeps the compiled field in this
// process, "disk" additionally reuses (and writes) <cache_dir>/<name>.<key>.jochrtf.
namespace joc::hrtf {
enum class CachePolicy { None, Memory, Disk };
struct SofaFieldRequest {
std::string sofa_path;
CompileOptions options;
CachePolicy policy = CachePolicy::Memory;
std::string cache_dir; // required for the disk policy
};
// Parses "none"/"memory"/"disk"; anything else is rejected.
Status parse_cache_policy(const std::string& text, CachePolicy* out);
// Returns the compiled field, reusing a valid cache when the policy allows it.
// `cache_path` (optional) receives the cache file that was read or written.
Status load_or_compile_sofa_field(const SofaFieldRequest& request, Field* out,
std::string* cache_path);
// Writes the field to `path` through a temporary file and an atomic rename.
Status save_jochrtf_atomic(const Field& field, const std::string& path);
// Verifies that a cache file belongs to `source_sha256` and `cache_key`.
Status validate_jochrtf(const std::string& path, const std::string& source_sha256,
const std::string& cache_key, Field* out);
} // namespace joc::hrtf
File diff suppressed because it is too large Load Diff
+103
View File
@@ -0,0 +1,103 @@
#pragma once
#include <cstdint>
#include <string>
#include <vector>
#include "foundation/status.h"
#include "hrtf/jochrtf.h"
#include "hrtf/sofa.h"
// SOFA SimpleFreeFieldHRIR -> compiled directional field, ported from the
// reference chain (sofa_canonical.py + sofa_hrtf_field.py + the public
// filterbank): the measurement shell is selected, one delay representation is
// separated, the FIRs are projected onto the 64-QMF/77-hybrid filterbank and the
// result is fitted with fifth-order ACN/N3D real spherical harmonics.
namespace joc::hrtf {
inline constexpr int kFieldOrder = 5;
inline constexpr int kFieldTerms = 36;
inline constexpr double kFieldSampleRateHz = 48000.0;
inline constexpr double kDefaultShellRadiusM = 1.0;
inline constexpr double kDefaultProjectionRidge = 1.0e-3;
inline constexpr double kDefaultSphericalHarmonicRidge = 1.0e-5;
inline constexpr const char* kCompilerVersion = "joc-sofa-compiler-v1";
inline constexpr const char* kPhasePolicyVersion = "sofa-delay-exactly-once-v1";
inline constexpr const char* kShConvention = "ACN/N3D real";
inline constexpr const char* kFilterbankTableVersion = "joc-public-64qmf-77hybrid-v1";
// SHA-256 of the standard filterbank archive the embedded tables came from. It
// participates in the cache key, so it is part of the file-format contract.
inline constexpr const char* kFilterbankArchiveSha256 =
"C05BEF4D26E96ECBD4694E2572F05DA400255C777BA5047300B9D3B1F81081CD";
struct CompileOptions {
double shell_radius_m = kDefaultShellRadiusM;
int order = kFieldOrder;
double projection_ridge = kDefaultProjectionRidge;
double sh_ridge = kDefaultSphericalHarmonicRidge;
};
// Canonical HRIR set: Data.IR and Data.Delay stay separate, the listener frame is
// applied to the source positions and the ears are ordered left/right.
struct CanonicalHrtf {
std::string source_path;
std::string source_sha256;
std::string convention;
std::string convention_version;
std::string processing_label;
double sample_rate_hz = 0.0;
std::uint32_t measurements = 0;
std::uint32_t taps = 0;
int left_receiver_index = 0;
int right_receiver_index = 1;
std::vector<double> source_position_cartesian_m; // [M,3] listener-local
std::vector<double> unit_directions; // [M,3]
std::vector<double> measurement_radius_m; // [M]
std::vector<double> hrir; // [M,2,N] canonical L/R
std::vector<double> delay_samples; // [M,2], not applied
};
// Port of load_simple_free_field_hrir(): strict SimpleFreeFieldHRIR import.
Status canonicalize_sofa(const SofaHrir& sofa, CanonicalHrtf* out);
// Compiles the canonical set into the runtime field (port of SofaHrtfField.fit).
Status compile_sofa_field(const SofaHrir& sofa, const CompileOptions& options, Field* out);
// The measurements on the shell nearest to radius_m; actual_radius_m receives the
// mean radius of that shell (upstream CanonicalHrtf.shell_indices).
std::vector<std::size_t> canonical_shell_indices(const CanonicalHrtf& canonical, double radius_m,
double* actual_radius_m);
// Compiles an already canonicalized set (used by tests and the cache layer).
Status compile_canonical_field(const CanonicalHrtf& canonical, const CompileOptions& options,
Field* out);
// Configuration hash that names the cache file (upstream compiled_hrtf_cache_key).
std::string compiled_hrtf_cache_key(const std::string& source_sha256, double sample_rate_hz,
double shell_radius_m, int order, double projection_ridge,
double sh_ridge);
// Payload hash over the four arrays (upstream _payload_sha256).
std::string field_payload_sha256(const Field& field);
// "<stem>.<first 20 key digits>.jochrtf", the upstream cache file name.
std::string cache_file_name(const std::string& display_name, const std::string& cache_key);
// Serializes the field as a .jochrtf cache the upstream loader also accepts.
Status write_jochrtf(const Field& field, const std::string& path);
// The analysis/gain/synthesis dictionary the projection solves against (dev check).
std::vector<double> hybrid_gain_synthesis_dictionary_for_check(std::size_t sample_count);
// Shell directions and their spherical Voronoi weights (dev check).
void shell_directions_and_weights_for_check(const SofaHrir& sofa, double radius_m,
std::vector<double>* directions,
std::vector<double>* weights);
// PublicAnalysis77 on a unit impulse, interleaved complex (dev check).
std::vector<double> analysis_impulse_for_check(std::size_t total_samples);
// The 77 hybrid-band centre frequencies at 48 kHz.
const std::vector<double>& hybrid_band_center_frequencies_hz();
} // namespace joc::hrtf