commit 6ca539f334b26762d59fce5294ddcec48e0d88b7
parent 35d5c0c5c1718f076f541c96be5f40e7dd00c04c
Author: d.levin256@gmail.com <d.levin256@gmail.com>
Date: Mon, 12 Nov 2018 16:24:33 +0300
Merge branch 'dev'
Diffstat:
14 files changed, 743 insertions(+), 23 deletions(-)
diff --git a/include/kfr/base/pointer.hpp b/include/kfr/base/pointer.hpp
@@ -208,7 +208,8 @@ CMT_INLINE expression_pointer<T> to_pointer(E&& expr)
{
static_assert(is_input_expression<E>::value, "E must be an expression");
std::shared_ptr<expression_resource> ptr = make_resource(std::move(expr));
- return expression_pointer<T>(ptr->instance(), internal::make_expression_vtable<T, E>(), std::move(ptr));
+ void* instance = ptr->instance();
+ return expression_pointer<T>(instance, internal::make_expression_vtable<T, E>(), std::move(ptr));
}
template <typename T, size_t key>
diff --git a/include/kfr/base/simd_clang.hpp b/include/kfr/base/simd_clang.hpp
@@ -96,7 +96,7 @@ CMT_INLINE void simd_write(T* dest, const simd_type<T, N>& value)
}
template <typename T, size_t N>
-struct vec : public vec_t<T, N>
+struct alignas(alignof(internal::simd_type<T, N>)) vec : public vec_t<T, N>
{
static_assert(is_simd_type<T>::value || !compound_type_traits<T>::is_scalar, "Invalid vector type");
@@ -222,12 +222,12 @@ struct vec : public vec_t<T, N>
// shuffle
template <size_t... indices>
- constexpr vec<value_type, sizeof...(indices)> shuffle(csizes_t<indices...>) const noexcept
+ vec<value_type, sizeof...(indices)> shuffle(csizes_t<indices...>) const noexcept
{
return __builtin_shufflevector(simd, simd, (indices >= N ? -1 : int(indices))...);
}
template <size_t... indices>
- constexpr vec<value_type, sizeof...(indices)> shuffle(const vec& y, csizes_t<indices...>) const noexcept
+ vec<value_type, sizeof...(indices)> shuffle(const vec& y, csizes_t<indices...>) const noexcept
{
return __builtin_shufflevector(simd, y.simd, (indices >= N * 2 ? -1 : int(indices))...);
}
diff --git a/include/kfr/base/types.hpp b/include/kfr/base/types.hpp
@@ -184,7 +184,7 @@ inline datatype operator|(datatype x, datatype y)
inline datatype operator&(datatype x, datatype y)
{
using type = underlying_type<datatype>;
- return static_cast<datatype>(static_cast<type>(x) | static_cast<type>(y));
+ return static_cast<datatype>(static_cast<type>(x) & static_cast<type>(y));
}
template <typename T>
diff --git a/include/kfr/cometa/array.hpp b/include/kfr/cometa/array.hpp
@@ -63,6 +63,10 @@ public:
constexpr array_ref(std::vector<T, Ts...>& vec) noexcept : m_data(vec.data()), m_size(vec.size())
{
}
+ constexpr array_ref(const std::initializer_list<T>& vec) noexcept
+ : m_data(vec.begin()), m_size(vec.size())
+ {
+ }
template <typename InputIter>
constexpr array_ref(InputIter first, InputIter last) noexcept
: m_data(std::addressof(*first)), m_size(std::distance(first, last))
@@ -122,4 +126,4 @@ inline array_ref<const T> make_array_ref(const std::vector<T>& cont)
{
return array_ref<const T>(cont.data(), cont.size());
}
-}
+} // namespace cometa
diff --git a/include/kfr/dsp.hpp b/include/kfr/dsp.hpp
@@ -28,6 +28,7 @@
#include "dsp/biquad_design.hpp"
#include "dsp/dcremove.hpp"
#include "dsp/delay.hpp"
+#include "dsp/ebu.hpp"
#include "dsp/fir.hpp"
#include "dsp/fir_design.hpp"
#include "dsp/fracdelay.hpp"
diff --git a/include/kfr/dsp/biquad.hpp b/include/kfr/dsp/biquad.hpp
@@ -307,7 +307,7 @@ class biquad_filter : public expression_filter<T>
{
public:
biquad_filter(const biquad_params<T>* bq, size_t count)
- : expression_filter<T>(to_pointer(biquad(bq, count, placeholder<T>())))
+ : expression_filter<T>(biquad(bq, count, placeholder<T>()))
{
}
diff --git a/include/kfr/dsp/ebu.hpp b/include/kfr/dsp/ebu.hpp
@@ -0,0 +1,332 @@
+#pragma once
+
+#include <vector>
+
+#include "../base.hpp"
+#include "../testo/assert.hpp"
+#include "biquad.hpp"
+#include "biquad_design.hpp"
+#include "speaker.hpp"
+#include "units.hpp"
+
+#pragma clang diagnostic push
+#if __has_warning("-Winaccessible-base")
+#pragma clang diagnostic ignored "-Winaccessible-base"
+#endif
+
+namespace kfr
+{
+
+template <typename T>
+KFR_SINTRIN T energy_to_loudness(T energy)
+{
+ return T(10) * log10(energy) - T(0.691);
+}
+
+template <typename T>
+KFR_SINTRIN T loudness_to_energy(T loudness)
+{
+ return exp10((loudness + T(0.691)) * T(0.1));
+}
+
+template <typename T>
+struct integrated_vec : public univector<T>
+{
+private:
+ void compute() const
+ {
+ const T z_total = mean(*this);
+ T relative_gate = energy_to_loudness(z_total) - 10;
+
+ T z = 0;
+ size_t num = 0;
+ for (T v : *this)
+ {
+ T lk = energy_to_loudness(v);
+ if (lk >= relative_gate)
+ {
+ z += v;
+ num++;
+ }
+ }
+ z /= num;
+ if (num >= 1)
+ {
+ m_integrated = energy_to_loudness(z);
+ }
+ else
+ {
+ m_integrated = -c_infinity<T>;
+ }
+
+ m_integrated_cached = true;
+ }
+
+public:
+ integrated_vec() : m_integrated(-c_infinity<T>), m_integrated_cached(false) {}
+ void push(T mean_square)
+ {
+ T lk = energy_to_loudness(mean_square);
+ if (lk >= T(-70.))
+ {
+ this->push_back(mean_square);
+ m_integrated_cached = false;
+ }
+ }
+ void reset()
+ {
+ m_integrated_cached = false;
+ this->clear();
+ }
+ T get() const
+ {
+ if (!m_integrated_cached)
+ {
+ compute();
+ }
+ return m_integrated;
+ }
+
+private:
+ mutable bool m_integrated_cached;
+ mutable T m_integrated;
+};
+
+template <typename T>
+struct lra_vec : public univector<T>
+{
+private:
+ void compute() const
+ {
+ m_range_high = -70.0;
+ m_range_low = -70.0;
+ static const T PRC_LOW = 0.10;
+ static const T PRC_HIGH = 0.95;
+
+ const T z_total = mean(*this);
+ const T relative_gate = energy_to_loudness(z_total) - 20;
+
+ if (this->size() < 2)
+ return;
+
+ const size_t start_index =
+ std::upper_bound(this->begin(), this->end(), loudness_to_energy(relative_gate)) - this->begin();
+ if (this->size() - start_index < 2)
+ return;
+ const size_t end_index = this->size() - 1;
+
+ const size_t low_index =
+ static_cast<size_t>(std::llround(start_index + (end_index - start_index) * PRC_LOW));
+ const size_t high_index =
+ static_cast<size_t>(std::llround(start_index + (end_index - start_index) * PRC_HIGH));
+ m_range_low = energy_to_loudness((*this)[low_index]);
+ m_range_high = energy_to_loudness((*this)[high_index]);
+
+ m_lra_cached = true;
+ }
+
+public:
+ lra_vec() : m_range_low(-70), m_range_high(-70), m_lra_cached(false) {}
+ void push(T mean_square)
+ {
+ const T lk = energy_to_loudness(mean_square);
+ if (lk >= -70)
+ {
+ auto it = std::upper_bound(this->begin(), this->end(), mean_square);
+ this->insert(it, mean_square);
+ m_lra_cached = false;
+ }
+ }
+ void reset()
+ {
+ m_lra_cached = false;
+ this->clear();
+ }
+ void get(T& low, T& high) const
+ {
+ if (!m_lra_cached)
+ compute();
+ low = m_range_low;
+ high = m_range_high;
+ }
+
+private:
+ mutable bool m_lra_cached;
+ mutable T m_range_low;
+ mutable T m_range_high;
+};
+
+template <typename T>
+KFR_SINTRIN expression_pointer<T> make_kfilter(int samplerate)
+{
+ const biquad_params<T> bq[] = {
+ biquad_highshelf(T(1681.81 / samplerate), T(+4.0)),
+ biquad_highpass(T(38.1106678246655 / samplerate), T(0.5)).normalized_all()
+ };
+ return to_pointer(biquad(bq, placeholder<T>()));
+}
+
+template <typename T>
+struct ebu_r128;
+
+template <typename T>
+struct ebu_channel
+{
+public:
+ friend struct ebu_r128<T>;
+ ebu_channel(int sample_rate, Speaker speaker, int packet_size_factor = 1, T input_gain = 1)
+ : m_sample_rate(sample_rate), m_speaker(speaker), m_input_gain(input_gain),
+ m_packet_size(sample_rate / 10 / packet_size_factor), m_kfilter(make_kfilter<T>(sample_rate)),
+ m_short_sum_of_squares(3000 / 100 * packet_size_factor),
+ m_momentary_sum_of_squares(400 / 100 * packet_size_factor), m_output_energy_gain(1.0),
+ m_buffer_cursor(0), m_short_sum_of_squares_cursor(0), m_momentary_sum_of_squares_cursor(0)
+ {
+ switch (speaker)
+ {
+ case Speaker::Lfe:
+ case Speaker::Lfe2:
+ m_output_energy_gain = 0.0;
+ break;
+ case Speaker::LeftSurround:
+ case Speaker::RightSurround:
+ m_output_energy_gain = dB_to_power(+1.5);
+ break;
+ default:
+ break;
+ }
+ reset();
+ }
+
+ void reset()
+ {
+ std::fill(m_short_sum_of_squares.begin(), m_short_sum_of_squares.end(), 0);
+ std::fill(m_momentary_sum_of_squares.begin(), m_momentary_sum_of_squares.end(), 0);
+ }
+
+ void process_packet(const T* src)
+ {
+ substitute(m_kfilter, to_pointer(make_univector(src, m_packet_size) * m_input_gain));
+ const T filtered_sum_of_squares = sumsqr(truncate(m_kfilter, m_packet_size));
+
+ m_short_sum_of_squares.ringbuf_write(m_short_sum_of_squares_cursor, filtered_sum_of_squares);
+ m_momentary_sum_of_squares.ringbuf_write(m_momentary_sum_of_squares_cursor, filtered_sum_of_squares);
+ }
+ Speaker get_speaker() const { return m_speaker; }
+
+private:
+ const Speaker m_speaker;
+ const T m_input_gain;
+ const int m_sample_rate;
+ const size_t m_packet_size;
+ expression_pointer<T> m_kfilter;
+ T m_output_energy_gain;
+ univector<T> m_buffer;
+ univector<T> m_short_sum_of_squares;
+ univector<T> m_momentary_sum_of_squares;
+ size_t m_buffer_cursor;
+ size_t m_short_sum_of_squares_cursor;
+ size_t m_momentary_sum_of_squares_cursor;
+};
+
+template <typename T>
+struct ebu_r128
+{
+public:
+ // Correct values for packet_size_factor: 1 (10Hz refresh rate), 2 (20Hz), 3 (30Hz)
+ ebu_r128(int sample_rate, const std::vector<Speaker>& channels, int packet_size_factor = 1)
+ : m_sample_rate(sample_rate), m_running(true), m_need_reset(false),
+ m_packet_size(sample_rate / 10 / packet_size_factor)
+ {
+ for (Speaker sp : channels)
+ {
+ m_channels.emplace_back(sample_rate, sp, packet_size_factor, 1);
+ }
+ }
+
+ int sample_rate() const { return m_sample_rate; }
+
+ size_t packet_size() const { return m_packet_size; }
+
+ void get_values(T& loudness_momentary, T& loudness_short, T& loudness_intergrated, T& loudness_range_low,
+ T& loudness_range_high)
+ {
+ T sum_of_mean_square_momentary = 0;
+ T sum_of_mean_square_short = 0;
+ for (size_t ch = 0; ch < m_channels.size(); ch++)
+ {
+ sum_of_mean_square_momentary += mean(m_channels[ch].m_momentary_sum_of_squares) / m_packet_size *
+ m_channels[ch].m_output_energy_gain;
+ sum_of_mean_square_short += mean(m_channels[ch].m_short_sum_of_squares) / m_packet_size *
+ m_channels[ch].m_output_energy_gain;
+ }
+ loudness_momentary = energy_to_loudness(sum_of_mean_square_momentary);
+ loudness_short = energy_to_loudness(sum_of_mean_square_short);
+ loudness_intergrated = m_integrated_buffer.get();
+ m_lra_buffer.get(loudness_range_low, loudness_range_high);
+ }
+
+ const ebu_channel<T>& operator[](size_t index) const { return m_channels[index]; }
+ size_t count() const { return m_channels.size(); }
+
+ void process_packet(const std::initializer_list<univector_dyn<T>>& source)
+ {
+ process_packet<tag_dynamic_vector>(source);
+ }
+ void process_packet(const std::initializer_list<univector_ref<T>>& source)
+ {
+ process_packet<tag_array_ref>(source);
+ }
+
+ template <size_t Tag>
+ void process_packet(const std::vector<univector<T, Tag>>& source)
+ {
+ T momentary = 0;
+ T shortterm = 0;
+ for (size_t ch = 0; ch < m_channels.size(); ch++)
+ {
+ TESTO_ASSERT(source[ch].size() == m_packet_size);
+ ebu_channel<T>& chan = m_channels[ch];
+ chan.process_packet(source[ch].data());
+ if (m_running)
+ {
+ momentary += mean(m_channels[ch].m_momentary_sum_of_squares) / m_packet_size *
+ m_channels[ch].m_output_energy_gain;
+ shortterm += mean(m_channels[ch].m_short_sum_of_squares) / m_packet_size *
+ m_channels[ch].m_output_energy_gain;
+ }
+ }
+
+ if (m_need_reset)
+ {
+ m_need_reset = false;
+ for (size_t ch = 0; ch < m_channels.size(); ch++)
+ {
+ m_channels[ch].reset();
+ }
+ m_integrated_buffer.reset();
+ m_lra_buffer.reset();
+ }
+ if (m_running)
+ {
+ m_integrated_buffer.push(momentary);
+ m_lra_buffer.push(shortterm);
+ }
+ }
+
+ void start() { m_running = true; }
+ void stop() { m_running = false; }
+ void reset() { m_need_reset = true; }
+
+private:
+ int m_sample_rate;
+ bool m_running;
+ bool m_need_reset;
+ size_t m_packet_size;
+ std::vector<ebu_channel<T>> m_channels;
+ integrated_vec<T> m_integrated_buffer;
+ lra_vec<T> m_lra_buffer;
+};
+
+} // namespace kfr
+
+#pragma clang diagnostic pop
diff --git a/include/kfr/dsp/weighting.hpp b/include/kfr/dsp/weighting.hpp
@@ -43,7 +43,7 @@ KFR_SINTRIN T weight_a_unnorm(T f)
}
template <typename T>
-constexpr static T weight_a_gain = reciprocal(weight_a_unnorm(T(1000.0)));
+const static T weight_a_gain = reciprocal(weight_a_unnorm(T(1000.0)));
template <typename T>
KFR_SINTRIN T aweighting(T f)
@@ -62,7 +62,7 @@ KFR_SINTRIN T weight_b_unnorm(T f)
}
template <typename T>
-constexpr static T weight_b_gain = reciprocal(weight_b_unnorm(T(1000.0)));
+const static T weight_b_gain = reciprocal(weight_b_unnorm(T(1000.0)));
template <typename T>
KFR_SINTRIN T bweighting(T f)
@@ -81,7 +81,7 @@ KFR_SINTRIN T weight_c_unnorm(T f)
}
template <typename T>
-constexpr static T weight_c_gain = reciprocal(weight_c_unnorm(T(1000.0)));
+const static T weight_c_gain = reciprocal(weight_c_unnorm(T(1000.0)));
template <typename T>
KFR_SINTRIN T cweighting(T f)
diff --git a/include/kfr/dsp/window.hpp b/include/kfr/dsp/window.hpp
@@ -127,7 +127,7 @@ struct expression_rectangular : input_expression
template <size_t N>
CMT_INLINE vec<T, N> operator()(cinput_t, size_t index, vec_t<T, N>) const
{
- using TI = utype<T>;
+ using TI = utype<T>;
const vec<TI, N> i = enumerate(vec<TI, N>()) + cast<TI>(index);
return select(i < cast<TI>(m_size), T(1), T(0));
}
@@ -369,11 +369,11 @@ struct expression_flattop : input_expression
CMT_INLINE vec<T, N> operator()(cinput_t cinput, size_t index, vec_t<T, N> y) const
{
const vec<T, N> n = linspace(cinput, index, y) * c_pi<T, 2>;
- constexpr T a0 = 1;
- constexpr T a1 = 1.93;
- constexpr T a2 = 1.29;
- constexpr T a3 = 0.388;
- constexpr T a4 = 0.028;
+ constexpr T a0 = 1;
+ constexpr T a1 = 1.93;
+ constexpr T a2 = 1.29;
+ constexpr T a3 = 0.388;
+ constexpr T a4 = 0.028;
return a0 - a1 * cos(n) + a2 * cos(2 * n) - a3 * cos(3 * n) + a4 * cos(4 * n);
}
size_t size() const { return m_size; }
@@ -399,6 +399,7 @@ struct expression_gaussian : input_expression
}
size_t size() const { return m_size; }
+
private:
window_linspace_m1_1_trunc<T> linspace;
T alpha;
@@ -452,7 +453,7 @@ KFR_WINDOW_BY_TYPE(flattop)
KFR_WINDOW_BY_TYPE(gaussian)
KFR_WINDOW_BY_TYPE(lanczos)
#undef KFR_WINDOW_BY_TYPE
-}
+} // namespace internal
/**
* @brief Returns template expression that generates Rrectangular window of length @c size
@@ -609,11 +610,11 @@ CMT_NOINLINE expression_pointer<T> window(size_t size, window_type type, identit
window_type::bohman, window_type::blackman, window_type::blackman_harris, window_type::kaiser,
window_type::flattop, window_type::gaussian, window_type::lanczos>(),
type,
- [=](auto win) {
+ [size, win_param, symmetry](auto win) {
constexpr window_type window = val_of(decltype(win)());
- return to_pointer<T>(
+ return to_pointer(
typename internal::window_by_type<window>::template type<T>(size, win_param, symmetry));
},
fn::returns<expression_pointer<T>>());
}
-}
+} // namespace kfr
diff --git a/include/kfr/testo/assert.hpp b/include/kfr/testo/assert.hpp
@@ -21,8 +21,11 @@ inline void assertion_failed(const std::string& string, const char* file, int li
errorln("Assertion failed at ", file, ":", line);
errorln(string);
errorln();
+ std::fflush(stderr);
}
+bool check_assertion(...);
+
template <typename Op, typename L, typename R>
bool check_assertion(const comparison<Op, L, R>& comparison, const char* expr, const char* file, int line)
{
@@ -44,7 +47,7 @@ bool check_assertion(const half_comparison<L>& comparison, const char* expr, con
{
assertion_failed(as_string(padleft(22, expr), " | ", comparison.left), file, line);
}
- return false;
+ return result;
}
#if defined(TESTO_ASSERTION_ON) || !(defined(NDEBUG) || defined(TESTO_ASSERTION_OFF))
@@ -67,4 +70,35 @@ bool check_assertion(const half_comparison<L>& comparison, const char* expr, con
#ifndef TESTO_NO_SHORT_MACROS
#define ASSERT TESTO_ASSERT
#endif
+
+template <typename OutType, typename InType>
+inline OutType safe_cast(const InType& val)
+{
+ static_assert(std::is_integral<InType>::value && std::is_integral<OutType>::value,
+ "safe_cast is for numeric types only");
+ if (std::is_signed<InType>::value && std::is_signed<OutType>::value) // S->S
+ {
+ ASSERT(val >= std::numeric_limits<OutType>::min());
+ ASSERT(val <= std::numeric_limits<OutType>::max());
+ }
+ else if (!std::is_signed<InType>::value && !std::is_signed<OutType>::value) // U->U
+ {
+ ASSERT(val <= std::numeric_limits<OutType>::max());
+ }
+ else if (std::is_signed<InType>::value && !std::is_signed<OutType>::value) // S->U
+ {
+ ASSERT(val >= 0);
+ ASSERT(val <= std::numeric_limits<OutType>::max());
+ // val will be converted to an unsigned number for the above comparison.
+ // it's safe because we've already checked that it is positive
+ }
+ else // U->S
+ {
+ ASSERT(val <= std::numeric_limits<OutType>::max());
+ // std::numeric_limits<OutType>::max() will be converted to an unsigned number for the above
+ // comparison. it's also safe here
+ }
+ return static_cast<OutType>(val);
}
+
+} // namespace testo
diff --git a/sources.cmake b/sources.cmake
@@ -84,6 +84,7 @@ set(
${PROJECT_SOURCE_DIR}/include/kfr/dsp/biquad_design.hpp
${PROJECT_SOURCE_DIR}/include/kfr/dsp/dcremove.hpp
${PROJECT_SOURCE_DIR}/include/kfr/dsp/delay.hpp
+ ${PROJECT_SOURCE_DIR}/include/kfr/dsp/ebu.hpp
${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir.hpp
${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir_design.hpp
${PROJECT_SOURCE_DIR}/include/kfr/dsp/fracdelay.hpp
@@ -104,7 +105,7 @@ set(
${PROJECT_SOURCE_DIR}/include/kfr/io/file.hpp
${PROJECT_SOURCE_DIR}/include/kfr/io/python_plot.hpp
${PROJECT_SOURCE_DIR}/include/kfr/io/tostring.hpp
- ${PROJECT_SOURCE_DIR}/include/kfr/testo/testo.hpp
${PROJECT_SOURCE_DIR}/include/kfr/testo/assert.hpp
${PROJECT_SOURCE_DIR}/include/kfr/testo/comparison.hpp
+ ${PROJECT_SOURCE_DIR}/include/kfr/testo/testo.hpp
)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -54,6 +54,8 @@ add_executable(complex_test complex_test.cpp ${KFR_SRC} ${TEST_SRC})
add_executable(base_test base_test.cpp ${KFR_SRC} ${TEST_SRC})
add_executable(expression_test expression_test.cpp ${KFR_SRC} ${TEST_SRC})
+add_executable(ebu_test ebu_test.cpp ${KFR_SRC} ${TEST_SRC})
+
if (ARM)
find_program(EMULATOR "qemu-arm")
else ()
diff --git a/tests/dsp_test.cpp b/tests/dsp_test.cpp
@@ -10,8 +10,230 @@
#include <kfr/dsp.hpp>
#include <kfr/io.hpp>
+#include <numeric>
+
using namespace kfr;
+struct TestFragment
+{
+ float gain; // dB
+ float duration; // seconds
+ float frequency; // Hz
+};
+
+struct TestFragmentMultichannel
+{
+ float gain_L_R; // dB
+ float gain_C; // dB
+ float gain_Ls_Rs; // dB
+ float duration; // seconds
+ float frequency; // Hz
+};
+
+template <typename T>
+static void ebu_test_stereo(int sample_rate, const std::initializer_list<TestFragment>& fragments, T refM,
+ T refS, T refI, T refLRA)
+{
+ ebu_r128<T> loudness(sample_rate, { Speaker::Left, Speaker::Right });
+
+ size_t total_length = 0;
+ for (const TestFragment& f : fragments)
+ {
+ total_length += static_cast<size_t>(f.duration * sample_rate);
+ }
+
+ univector<T> left_right(total_length);
+ size_t pos = 0;
+ for (const TestFragment& f : fragments)
+ {
+ const size_t len = static_cast<size_t>(f.duration * sample_rate);
+ left_right.slice(pos, len) = dB_to_amp(f.gain) * sinenorm(phasor<float>(f.frequency, sample_rate));
+ pos += len;
+ }
+
+ for (size_t i = 0; i < total_length / loudness.packet_size(); i++)
+ {
+ loudness.process_packet({ left_right.slice(i * loudness.packet_size(), loudness.packet_size()),
+ left_right.slice(i * loudness.packet_size(), loudness.packet_size()) });
+ }
+ T M, S, I, RL, RH;
+ loudness.get_values(M, S, I, RL, RH);
+ if (!std::isnan(refM))
+ CHECK(std::abs(M - refM) < 0.05f);
+ if (!std::isnan(refS))
+ CHECK(std::abs(S - refS) < 0.05f);
+ if (!std::isnan(refI))
+ CHECK(std::abs(I - refI) < 0.05f);
+ if (!std::isnan(refLRA))
+ CHECK(std::abs((RH - RL) - refLRA) < 0.05f);
+}
+
+template <typename T>
+static void ebu_test_multichannel(int sample_rate,
+ const std::initializer_list<TestFragmentMultichannel>& fragments, T refM,
+ T refS, T refI, T refLRA)
+{
+ ebu_r128<T> loudness(sample_rate, { Speaker::Left, Speaker::Right, Speaker::Center, Speaker::LeftSurround,
+ Speaker::RightSurround });
+
+ size_t total_length = 0;
+ for (const TestFragmentMultichannel& f : fragments)
+ {
+ total_length += static_cast<size_t>(f.duration * sample_rate);
+ }
+
+ univector<T> left_right(total_length);
+ univector<T> center(total_length);
+ univector<T> surround(total_length);
+ size_t pos = 0;
+ for (const TestFragmentMultichannel& f : fragments)
+ {
+ const size_t len = static_cast<size_t>(f.duration * sample_rate);
+ left_right.slice(pos, len) =
+ dB_to_amp(f.gain_L_R) * sinenorm(phasor<float>(f.frequency, sample_rate));
+ center.slice(pos, len) = dB_to_amp(f.gain_C) * sinenorm(phasor<float>(f.frequency, sample_rate));
+ surround.slice(pos, len) =
+ dB_to_amp(f.gain_Ls_Rs) * sinenorm(phasor<float>(f.frequency, sample_rate));
+ pos += len;
+ }
+
+ for (size_t i = 0; i < total_length / loudness.packet_size(); i++)
+ {
+ loudness.process_packet({ left_right.slice(i * loudness.packet_size(), loudness.packet_size()),
+ left_right.slice(i * loudness.packet_size(), loudness.packet_size()),
+ center.slice(i * loudness.packet_size(), loudness.packet_size()),
+ surround.slice(i * loudness.packet_size(), loudness.packet_size()),
+ surround.slice(i * loudness.packet_size(), loudness.packet_size()) });
+ }
+ T M, S, I, RL, RH;
+ loudness.get_values(M, S, I, RL, RH);
+ if (!std::isnan(refM))
+ CHECK(std::abs(M - refM) < 0.05f);
+ if (!std::isnan(refS))
+ CHECK(std::abs(S - refS) < 0.05f);
+ if (!std::isnan(refI))
+ CHECK(std::abs(I - refI) < 0.05f);
+ if (!std::isnan(refLRA))
+ CHECK(std::abs((RH - RL) - refLRA) < 0.05f);
+}
+
+TEST(ebu_stereo_1_and_2)
+{
+ testo::matrix(named("type") = ctypes_t<float, double>{},
+ named("sample_rate") = std::vector<int>{ 44100, 48000 }, [](auto type, int sample_rate) {
+ using T = type_of<decltype(type)>;
+
+ ebu_test_stereo<T>(sample_rate, { { -23.f, 20.f, 1000.f } }, -23.f, -23.f, -23.f, NAN);
+ ebu_test_stereo<T>(sample_rate, { { -33.f, 20.f, 1000.f } }, -33.f, -33.f, -33.f, NAN);
+ });
+}
+
+TEST(ebu_stereo_3_4_and_5)
+{
+ testo::matrix(named("type") = ctypes_t<float, double>{},
+ named("sample_rate") = std::vector<int>{ 44100, 48000 }, [](auto type, int sample_rate) {
+ using T = type_of<decltype(type)>;
+
+ ebu_test_stereo<T>(
+ sample_rate,
+ { { -36.f, 10.f, 1000.f }, { -23.f, 60.f, 1000.f }, { -36.f, 10.f, 1000.f } }, NAN,
+ NAN, -23.f, NAN);
+ ebu_test_stereo<T>(sample_rate,
+ { { -72.f, 10.f, 1000.f },
+ { -36.f, 10.f, 1000.f },
+ { -23.f, 60.f, 1000.f },
+ { -36.f, 10.f, 1000.f },
+ { -72.f, 10.f, 1000.f } },
+ NAN, NAN, -23.f, NAN);
+ });
+}
+
+TEST(ebu_multichannel_6)
+{
+ testo::matrix(named("type") = ctypes_t<float, double>{},
+ named("sample_rate") = std::vector<int>{ 44100, 48000 }, [](auto type, int sample_rate) {
+ using T = type_of<decltype(type)>;
+
+ ebu_test_multichannel<T>(sample_rate, { { -28.f, -24.f, -30.f, 20.f, 1000.f } }, NAN,
+ NAN, -23.f, NAN);
+ });
+}
+
+TEST(ebu_stereo_9)
+{
+ testo::matrix(named("type") = ctypes_t<float, double>{},
+ named("sample_rate") = std::vector<int>{ 44100, 48000 }, [](auto type, int sample_rate) {
+ using T = type_of<decltype(type)>;
+
+ ebu_test_stereo<T>(sample_rate,
+ { { -20.f, 1.34f, 1000.f },
+ { -30.f, 1.66f, 1000.f },
+ { -20.f, 1.34f, 1000.f },
+ { -30.f, 1.66f, 1000.f },
+ { -20.f, 1.34f, 1000.f },
+ { -30.f, 1.66f, 1000.f },
+ { -20.f, 1.34f, 1000.f },
+ { -30.f, 1.66f, 1000.f },
+ { -20.f, 1.34f, 1000.f },
+ { -30.f, 1.66f, 1000.f } },
+ NAN, -23.f, NAN, NAN);
+ });
+}
+
+TEST(ebu_stereo_12)
+{
+ testo::matrix(named("type") = ctypes_t<float, double>{},
+ named("sample_rate") = std::vector<int>{ 44100, 48000 }, [](auto type, int sample_rate) {
+ using T = type_of<decltype(type)>;
+
+ ebu_test_stereo<T>(
+ sample_rate,
+ { { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f },
+ { -30.f, 0.22f, 1000.f }, { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f },
+ { -20.f, 0.18f, 1000.f }, { -30.f, 0.22f, 1000.f } },
+ -23.f, NAN, NAN, NAN);
+ });
+}
+
+TEST(ebu_lra_1_2_3_and_4)
+{
+ testo::matrix(named("type") = ctypes_t<float, double>{},
+ named("sample_rate") = std::vector<int>{ 44100, 48000 }, [](auto type, int sample_rate) {
+ using T = type_of<decltype(type)>;
+
+ ebu_test_stereo<T>(sample_rate, { { -20.f, 20.f, 1000.f }, { -30.f, 20.f, 1000.f } },
+ NAN, NAN, NAN, 10.f);
+
+ ebu_test_stereo<T>(sample_rate, { { -20.f, 20.f, 1000.f }, { -15.f, 20.f, 1000.f } },
+ NAN, NAN, NAN, 5.f);
+
+ ebu_test_stereo<T>(sample_rate, { { -40.f, 20.f, 1000.f }, { -20.f, 20.f, 1000.f } },
+ NAN, NAN, NAN, 20.f);
+
+ ebu_test_stereo<T>(sample_rate,
+ { { -50.f, 20.f, 1000.f },
+ { -35.f, 20.f, 1000.f },
+ { -20.f, 20.f, 1000.f },
+ { -35.f, 20.f, 1000.f },
+ { -50.f, 20.f, 1000.f } },
+ NAN, NAN, NAN, 15.f);
+ });
+}
+
TEST(delay)
{
const univector<float, 33> v1 = counter() + 100;
@@ -58,7 +280,7 @@ TEST(mixdown_stereo)
TEST(phasor)
{
- constexpr fbase sr = 44100.0;
+ constexpr fbase sr = 44100.0;
univector<fbase, 100> v1 = sinenorm(phasor(15000, sr));
univector<fbase, 100> v2 = sin(constants<fbase>::pi_s(2) * counter(0, 15000 / sr));
CHECK(rms(v1 - v2) < 1.e-5);
diff --git a/tests/ebu_test.cpp b/tests/ebu_test.cpp
@@ -0,0 +1,122 @@
+/**
+ * KFR (http://kfrlib.com)
+ * Copyright (C) 2016 D Levin
+ * See LICENSE.txt for details
+ */
+
+#include <kfr/testo/testo.hpp>
+
+#include <kfr/base.hpp>
+#include <kfr/dsp.hpp>
+#include <kfr/io.hpp>
+
+using namespace kfr;
+
+int main(int argc, char** argv)
+{
+ if (argc < 3)
+ {
+ println("Usage: ebu_test INPUT_IN_F32_RAW_FORMAT CHANNEL_NUMBER");
+ return 1;
+ }
+
+ // Prepare
+ FILE* f = fopen(argv[1], "rb");
+ const int channel_number = atoi(argv[2]);
+ if (channel_number < 1 || channel_number > 6)
+ {
+ println("Incorrect number of channels");
+ return 1;
+ }
+ fseek(f, 0, SEEK_END);
+ uintmax_t size = ftell(f);
+ fseek(f, 0, SEEK_SET);
+ if (size % (sizeof(float) * channel_number))
+ {
+ println("Incorrect file size");
+ return 1;
+ }
+
+ // Read file
+ const size_t length = size / (sizeof(float) * channel_number);
+ univector<float> interleaved(size / sizeof(float));
+ size_t read_len = fread(interleaved.data(), 1, size, f);
+ if (read_len != size)
+ {
+ println("Can't read file");
+ return 1;
+ }
+
+ // Deinterleave
+ univector<univector<float>> data(channel_number, univector<float>(length));
+ for (size_t ch = 0; ch < channel_number; ++ch)
+ {
+ for (size_t i = 0; i < length; ++i)
+ {
+ data[ch][i] = interleaved[i * channel_number + ch];
+ }
+ }
+
+ std::vector<Speaker> speakers;
+ switch (channel_number)
+ {
+ case 1:
+ speakers = { Speaker::Mono };
+ break;
+ case 2:
+ speakers = { Speaker::Left, Speaker::Right };
+ break;
+ case 3:
+ speakers = { Speaker::Left, Speaker::Right, Speaker::Center };
+ break;
+ case 4:
+ speakers = { Speaker::Left, Speaker::Right, Speaker::LeftSurround, Speaker::RightSurround };
+ break;
+ case 5:
+ speakers = { Speaker::Left, Speaker::Right, Speaker::Center, Speaker::LeftSurround,
+ Speaker::RightSurround };
+ break;
+ case 6:
+ speakers = { Speaker::Left, Speaker::Right, Speaker::Center,
+ Speaker::LeftSurround, Speaker::RightSurround, Speaker::Lfe };
+ break;
+ }
+
+ ebu_r128<float> loudness(48000, speakers);
+
+ float M, S, I, RL, RH;
+ float maxM = -HUGE_VALF, maxS = -HUGE_VALF;
+ for (size_t i = 0; i < length / loudness.packet_size(); i++)
+ {
+ std::vector<univector_ref<float>> channels;
+ for (size_t ch = 0; ch < channel_number; ++ch)
+ {
+ channels.push_back(data[ch].slice(i * loudness.packet_size(), loudness.packet_size()));
+ }
+ loudness.process_packet(channels);
+ loudness.get_values(M, S, I, RL, RH);
+ maxM = std::max(maxM, M);
+ maxS = std::max(maxS, S);
+ }
+
+ {
+ // For file-based measurements, the signal should be followed by at least 1.5 s of silence
+ std::vector<univector_dyn<float>> channels(channel_number,
+ univector_dyn<float>(loudness.packet_size()));
+ for (size_t i = 0; i < 15; ++i)
+ loudness.process_packet(channels);
+ float dummyM, dummyS, dummyI;
+ loudness.get_values(dummyM, dummyS, dummyI, RL, RH);
+ }
+
+ println(argv[1]);
+ println("M = ", M);
+ println("S = ", S);
+ println("I = ", I);
+ println("LRA = ", RH - RL);
+ println("maxM = ", maxM);
+ println("maxS = ", maxS);
+ println();
+
+ return 0;
+}
