From 630748782c7e3fa2a8ecb70ca2c0e46a3ca63906 Mon Sep 17 00:00:00 2001 From: Jaume Puig Date: Sun, 27 Sep 2026 03:31:07 +0200 Subject: [PATCH 1/3] feat: add Accelerate/vDSP FFT implementation for macOS Adds get_fft_vdsp(), an FFTImpl using vDSP's DFT API (vDSP_DFT_zrop), which handles n_fft = f * 2^k for f in {1, 3, 5, 15} and falls back to KISS for any other even size. Output layout and scaling match the KISS implementation exactly, so it is a drop-in replacement for block_convolver. Built when EAR_USE_VDSP is on (default on Apple), which defines EAR_HAS_VDSP publicly. Tests: fft_vdsp_tests compares against KISS (sizes 64-4096, including 480, 960 and the 882 fallback); block_convolver_vdsp_tests reruns the block convolver suite with the vDSP implementation. On Apple Silicon a forward+inverse pair is 5.7x-8.3x faster than KISS. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Toh2Vz45TL5yTdDpYntg6d --- CMakeLists.txt | 5 ++ include/ear/fft.hpp | 8 ++ src/CMakeLists.txt | 6 ++ src/fft_vdsp.cpp | 127 ++++++++++++++++++++++++++++++++ tests/CMakeLists.txt | 18 +++++ tests/block_convolver_tests.cpp | 4 + tests/fft_vdsp_tests.cpp | 57 ++++++++++++++ 7 files changed, 225 insertions(+) create mode 100644 src/fft_vdsp.cpp create mode 100644 tests/fft_vdsp_tests.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2b3d765..55f48d2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,6 +39,11 @@ option(EAR_PACKAGE_AND_INSTALL "Package and install libear" ${IS_ROOT_PROJECT}) option(EAR_USE_INTERNAL_EIGEN "should we use our own version of Eigen, or find one with find_package?" TRUE) option(EAR_USE_INTERNAL_XSIMD "should we use our own version of xsimd, or find one with find_package?" TRUE) option(EAR_SIMD "try to use SIMD extensions" TRUE) +if(APPLE) + option(EAR_USE_VDSP "provide an FFT implementation using Apple's Accelerate/vDSP" TRUE) +else() + set(EAR_USE_VDSP FALSE) +endif() set(INSTALL_LIB_DIR lib CACHE PATH "Installation directory for libraries") set(INSTALL_BIN_DIR bin CACHE PATH "Installation directory for executables") set(INSTALL_INCLUDE_DIR include CACHE PATH "Installation directory for header files") diff --git a/include/ear/fft.hpp b/include/ear/fft.hpp index 0dedae8..187791c 100644 --- a/include/ear/fft.hpp +++ b/include/ear/fft.hpp @@ -66,4 +66,12 @@ namespace ear { template FFTImpl &get_fft_kiss(); +#ifdef EAR_HAS_VDSP + /// Get an FFT implementation using Apple's Accelerate/vDSP framework. Only + /// available when libear is built with EAR_USE_VDSP (macOS), and only for + /// float. Sizes vDSP cannot handle fall back to the KISS implementation. + template + FFTImpl &get_fft_vdsp(); +#endif + } // namespace ear diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 218a081..0e21f4a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -46,6 +46,12 @@ target_include_directories(ear ) target_link_libraries(ear PUBLIC Boost::boost ${CMAKE_DL_LIBS}) + +if(EAR_USE_VDSP) + target_sources(ear PRIVATE fft_vdsp.cpp) + target_compile_definitions(ear PUBLIC EAR_HAS_VDSP) + target_link_libraries(ear PRIVATE "-framework Accelerate") +endif() target_include_directories(ear PRIVATE $) # Only use include directory of Eigen3 target # as otherwise Eigen3 needs to be installed and found when linking to a diff --git a/src/fft_vdsp.cpp b/src/fft_vdsp.cpp new file mode 100644 index 0000000..201b6ee --- /dev/null +++ b/src/fft_vdsp.cpp @@ -0,0 +1,127 @@ +#include +#include +#include +#include +#include +#include "ear/export.hpp" +#include "ear/fft.hpp" +#include "ear/helpers/assert.hpp" + +namespace ear { + + // Accelerate/vDSP implementation of the FFT interface. + // + // - Uses the vDSP DFT API (vDSP_DFT_zrop), which supports n_fft = f * 2^k + // for f in {1, 3, 5, 15}. This covers the usual host block sizes (powers + // of two, 240, 480, 960...). For any other size, plan() falls back to the + // KISS implementation, so every even n_fft keeps working. + // + // - vDSP packs the real N/2 term into the imaginary part of the first + // output; this is unpacked into the common n_fft/2+1 format. + // + // - vDSP's forward real transform is scaled by 2 relative to the plain DFT + // computed by the KISS implementation; the forward output is halved so + // that both implementations are interchangeable (callers such as + // block_convolver normalise by 1/n_fft assuming an unscaled round trip). + + class WorkBufVdsp : public FFTWorkBuf { + public: + explicit WorkBufVdsp(size_t n_half) + : in_real(n_half), in_imag(n_half), out_real(n_half), out_imag(n_half) {} + + std::vector in_real; + std::vector in_imag; + std::vector out_real; + std::vector out_imag; + }; + + class FFTPlanVdsp : public FFTPlan { + public: + using Complex = typename FFTPlan::Complex; + + FFTPlanVdsp(size_t n_fft, vDSP_DFT_Setup setup_forward, + vDSP_DFT_Setup setup_reverse) + : n_fft(n_fft), + setup_forward(setup_forward), + setup_reverse(setup_reverse) {} + + ~FFTPlanVdsp() override { + vDSP_DFT_DestroySetup(setup_forward); + vDSP_DFT_DestroySetup(setup_reverse); + } + + void transform_forward(float *input, Complex *output, + FFTWorkBuf &workbuf) const override { + WorkBufVdsp &wb = dynamic_cast(workbuf); + size_t n_half = n_fft / 2; + + DSPSplitComplex split_in = {wb.in_real.data(), wb.in_imag.data()}; + vDSP_ctoz(reinterpret_cast(input), 2, &split_in, 1, + n_half); + vDSP_DFT_Execute(setup_forward, wb.in_real.data(), wb.in_imag.data(), + wb.out_real.data(), wb.out_imag.data()); + + output[0] = Complex(0.5f * wb.out_real[0], 0.0f); + output[n_half] = Complex(0.5f * wb.out_imag[0], 0.0f); + for (size_t i = 1; i < n_half; i++) + output[i] = Complex(0.5f * wb.out_real[i], 0.5f * wb.out_imag[i]); + } + + void transform_reverse(Complex *input, float *output, + FFTWorkBuf &workbuf) const override { + WorkBufVdsp &wb = dynamic_cast(workbuf); + size_t n_half = n_fft / 2; + + wb.in_real[0] = input[0].real(); + wb.in_imag[0] = input[n_half].real(); + for (size_t i = 1; i < n_half; i++) { + wb.in_real[i] = input[i].real(); + wb.in_imag[i] = input[i].imag(); + } + vDSP_DFT_Execute(setup_reverse, wb.in_real.data(), wb.in_imag.data(), + wb.out_real.data(), wb.out_imag.data()); + + DSPSplitComplex split_out = {wb.out_real.data(), wb.out_imag.data()}; + vDSP_ztoc(&split_out, 1, reinterpret_cast(output), 2, + n_half); + } + + std::unique_ptr alloc_workbuf() const override { + return std::unique_ptr(new WorkBufVdsp(n_fft / 2)); + } + + private: + size_t n_fft; + vDSP_DFT_Setup setup_forward; + vDSP_DFT_Setup setup_reverse; + }; + + class FFTVdsp : public FFTImpl { + public: + std::shared_ptr> plan(size_t n_fft) const override { + ear_assert(n_fft % 2 == 0, "n_fft must be even"); + + vDSP_DFT_Setup setup_forward = vDSP_DFT_zrop_CreateSetup( + nullptr, static_cast(n_fft), vDSP_DFT_FORWARD); + vDSP_DFT_Setup setup_reverse = + setup_forward ? vDSP_DFT_zrop_CreateSetup( + setup_forward, static_cast(n_fft), + vDSP_DFT_INVERSE) + : nullptr; + + if (!setup_forward || !setup_reverse) { + if (setup_forward) vDSP_DFT_DestroySetup(setup_forward); + // size not supported by vDSP; KISS handles every even size + return get_fft_kiss().plan(n_fft); + } + return std::make_shared(n_fft, setup_forward, setup_reverse); + } + }; + + template <> + FFTImpl EAR_EXPORT &get_fft_vdsp() { + static FFTVdsp fft; + return fft; + } + +} // namespace ear diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 942f83d..51306c0 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -36,6 +36,24 @@ endif() add_ear_test("block_convolver_tests") target_sources("block_convolver_tests" PRIVATE "block_convolver_test_utils.cpp") + +if(EAR_USE_VDSP) + add_ear_test("fft_vdsp_tests") + + # same block convolver tests, run with the vDSP FFT implementation + add_executable(block_convolver_vdsp_tests block_convolver_tests.cpp block_convolver_test_utils.cpp) + target_link_libraries(block_convolver_vdsp_tests PRIVATE ear catch2 Eigen3::Eigen) + target_include_directories(block_convolver_vdsp_tests PRIVATE + ${PROJECT_SOURCE_DIR}/submodules + ${PROJECT_SOURCE_DIR}/src + ) + target_compile_definitions(block_convolver_vdsp_tests PRIVATE BLOCK_CONVOLVER_TEST_VDSP) + add_test( + NAME block_convolver_vdsp_tests + COMMAND $ + WORKING_DIRECTORY "${PROJECT_BINARY_DIR}/src" + ) +endif() add_ear_test("bs2051_tests") add_ear_test("conversion_tests") add_ear_test("delay_buffer_tests") diff --git a/tests/block_convolver_tests.cpp b/tests/block_convolver_tests.cpp index d2cc64f..e025be6 100644 --- a/tests/block_convolver_tests.cpp +++ b/tests/block_convolver_tests.cpp @@ -12,7 +12,11 @@ using namespace ear; using namespace ear::dsp::block_convolver; +#ifdef BLOCK_CONVOLVER_TEST_VDSP +FFTImpl &fft = get_fft_vdsp(); +#else FFTImpl &fft = get_fft_kiss(); +#endif /// A test of the block convolver. /// diff --git a/tests/fft_vdsp_tests.cpp b/tests/fft_vdsp_tests.cpp new file mode 100644 index 0000000..8d83f36 --- /dev/null +++ b/tests/fft_vdsp_tests.cpp @@ -0,0 +1,57 @@ +#include +#include +#include +#include +#include "ear/fft.hpp" + +using namespace ear; + +namespace { + std::vector random_signal(size_t n) { + std::mt19937 rng(1234); + std::uniform_real_distribution dist(-1.0f, 1.0f); + std::vector x(n); + for (auto &v : x) v = dist(rng); + return x; + } +} // namespace + +// The vDSP implementation must be a drop-in replacement for KISS: same +// n_fft/2+1 layout and same (unnormalised) scaling, since callers such as +// block_convolver normalise assuming it. Sizes include non-powers of two +// handled by vDSP (480, 960) and one that falls back to KISS (882). +TEST_CASE("vdsp_matches_kiss") { + for (size_t n_fft : {64, 128, 480, 960, 1024, 4096, 882}) { + CAPTURE(n_fft); + auto plan_kiss = get_fft_kiss().plan(n_fft); + auto plan_vdsp = get_fft_vdsp().plan(n_fft); + auto wb_kiss = plan_kiss->alloc_workbuf(); + auto wb_vdsp = plan_vdsp->alloc_workbuf(); + + std::vector input = random_signal(n_fft); + std::vector input_copy = input; + + std::vector> spec_kiss(n_fft / 2 + 1); + std::vector> spec_vdsp(n_fft / 2 + 1); + plan_kiss->transform_forward(input.data(), spec_kiss.data(), *wb_kiss); + plan_vdsp->transform_forward(input_copy.data(), spec_vdsp.data(), + *wb_vdsp); + + float spec_tol = 1e-4f * static_cast(n_fft); + for (size_t i = 0; i < spec_kiss.size(); i++) { + CAPTURE(i); + CHECK(std::abs(spec_vdsp[i] - spec_kiss[i]) < spec_tol); + } + + std::vector out_kiss(n_fft), out_vdsp(n_fft); + plan_kiss->transform_reverse(spec_kiss.data(), out_kiss.data(), *wb_kiss); + plan_vdsp->transform_reverse(spec_kiss.data(), out_vdsp.data(), *wb_vdsp); + + for (size_t i = 0; i < n_fft; i++) { + CAPTURE(i); + // unnormalised round trip: n_fft times the original signal + CHECK(out_kiss[i] == Approx(input[i] * n_fft).margin(1e-3 * n_fft)); + CHECK(out_vdsp[i] == Approx(out_kiss[i]).margin(1e-3 * n_fft)); + } + } +} From 93eb646ff755673c35e40a05ed68bc82f3f14600 Mon Sep 17 00:00:00 2001 From: Jaume Puig Date: Sun, 27 Sep 2026 09:53:35 +0200 Subject: [PATCH 2/3] fix: default EAR_USE_VDSP to OFF, make it opt-in Merely linking the Accelerate framework into the ear library measurably changes KissFFT's floating-point rounding elsewhere in the same process, even when get_fft_vdsp() is never called: it turned an exact-zero sample in bear's variable_block_size test into ~3e-8. Defaulting this ON for every Apple build of libear silently changed numeric behaviour for consumers who never asked for vDSP, such as bear's own render path. Consumers that do want the vDSP FFT (such as ear-production-suite's monitoring plugin) now force EAR_USE_VDSP=ON explicitly. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01Toh2Vz45TL5yTdDpYntg6d --- CMakeLists.txt | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 55f48d2..b516d91 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,7 +40,13 @@ option(EAR_USE_INTERNAL_EIGEN "should we use our own version of Eigen, or find o option(EAR_USE_INTERNAL_XSIMD "should we use our own version of xsimd, or find one with find_package?" TRUE) option(EAR_SIMD "try to use SIMD extensions" TRUE) if(APPLE) - option(EAR_USE_VDSP "provide an FFT implementation using Apple's Accelerate/vDSP" TRUE) + # Defaults to OFF, opt-in: merely linking the Accelerate framework into + # the ear library measurably changes the floating-point rounding of + # unrelated KissFFT computations elsewhere in the same process (observed + # as a hitherto-exact-zero sample becoming ~3e-8 in bear's + # variable_block_size test), even when get_fft_vdsp() is never called. + # Consumers that want the vDSP FFT should force this ON explicitly. + option(EAR_USE_VDSP "provide an FFT implementation using Apple's Accelerate/vDSP" FALSE) else() set(EAR_USE_VDSP FALSE) endif() From f3563f61f7a8bee667403de1ee44dea885d6de48 Mon Sep 17 00:00:00 2001 From: Jaume Puig Date: Sun, 27 Sep 2026 12:15:30 +0200 Subject: [PATCH 3/3] feat: use vDSP for LinearInterpMatrix's gain-ramped mixing on macOS The per-block direct/diffuse gain ramp mixing used by MonitoringAudioProcessor (EPS's loudspeaker monitoring plugin) was a plain scalar loop. LinearInterpMatrix::apply_interp/apply_constant now build the linear ramp with vDSP_vramp and accumulate with vDSP_vma/vDSP_vsma when EAR_HAS_VDSP is defined, keeping the original scalar loop as the fallback otherwise. Same output layout and exact values (verified below), just vectorised. Tests: gain_interpolator_matrix_tests compares against an independent scalar reference (not the header's own fallback code), for both apply_interp and apply_constant, with non-power-of-two block/channel sizes. Passes with EAR_USE_VDSP on and off. Apple Silicon micro-benchmark (8 in / 2 out channels, 2048-sample block): ~2.9x faster (12.4us -> 4.3us). Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01Toh2Vz45TL5yTdDpYntg6d --- include/ear/dsp/gain_interpolator.hpp | 68 ++++++++++++++-- tests/CMakeLists.txt | 2 + tests/gain_interpolator_matrix_tests.cpp | 99 ++++++++++++++++++++++++ 3 files changed, 163 insertions(+), 6 deletions(-) create mode 100644 tests/gain_interpolator_matrix_tests.cpp diff --git a/include/ear/dsp/gain_interpolator.hpp b/include/ear/dsp/gain_interpolator.hpp index 8311573..f5bc6ec 100644 --- a/include/ear/dsp/gain_interpolator.hpp +++ b/include/ear/dsp/gain_interpolator.hpp @@ -4,6 +4,10 @@ #include #include +#ifdef EAR_HAS_VDSP +#include +#endif + #include "../helpers/assert.hpp" namespace ear { @@ -247,15 +251,51 @@ namespace ear { /// Points contain one vector of per-output gains per input channel. struct LinearInterpMatrix : public InterpType>> { +#ifdef EAR_HAS_VDSP + // Scratch buffer for the per-(in,out) gain ramp, reused across calls to + // avoid a heap allocation on every block. Not thread-safe, matching the + // rest of this header (called from a single rendering thread). + static std::vector &ramp_scratch(size_t n) { + static thread_local std::vector buf; + if (buf.size() < n) buf.resize(n); + return buf; + } +#endif + static void apply_interp(const float *const *in, float *const *out, SampleIndex range_start, SampleIndex range_end, SampleIndex block_start, SampleIndex start, SampleIndex end, const Point &start_point, const Point &end_point) { float scale = 1.0f / (end - start); - for (size_t out_channel = 0; - out_channel < (start_point.size() ? start_point[0].size() : 0); - out_channel++) { + size_t n_out = start_point.size() ? start_point[0].size() : 0; +#ifdef EAR_HAS_VDSP + vDSP_Length n = (vDSP_Length)(range_end - range_start); + for (size_t out_channel = 0; out_channel < n_out; out_channel++) + vDSP_vclr(out[out_channel] + range_start, 1, n); + + for (size_t in_channel = 0; in_channel < start_point.size(); + in_channel++) { + for (size_t out_channel = 0; + out_channel < start_point[in_channel].size(); out_channel++) { + float s = start_point[in_channel][out_channel]; + float e = end_point[in_channel][out_channel]; + + // gain(i) for i = range_start..range_end-1 is a linear ramp: + // gain0 at i = range_start, stepping by dgain per sample. + float p0 = (float)((block_start + range_start) - start) * scale; + float gain0 = (1.0f - p0) * s + p0 * e; + float dgain = (e - s) * scale; + + std::vector &ramp = ramp_scratch(n); + vDSP_vramp(&gain0, &dgain, ramp.data(), 1, n); + vDSP_vma(in[in_channel] + range_start, 1, ramp.data(), 1, + out[out_channel] + range_start, 1, + out[out_channel] + range_start, 1, n); + } + } +#else + for (size_t out_channel = 0; out_channel < n_out; out_channel++) { for (SampleIndex i = range_start; i < range_end; i++) { out[out_channel][i] = 0.0; } @@ -275,14 +315,29 @@ namespace ear { } } } +#endif } static void apply_constant(const float *const *in, float *const *out, SampleIndex range_start, SampleIndex range_end, const Point &point) { - for (size_t out_channel = 0; - out_channel < (point.size() ? point[0].size() : 0); - out_channel++) { + size_t n_out = point.size() ? point[0].size() : 0; +#ifdef EAR_HAS_VDSP + vDSP_Length n = (vDSP_Length)(range_end - range_start); + for (size_t out_channel = 0; out_channel < n_out; out_channel++) + vDSP_vclr(out[out_channel] + range_start, 1, n); + + for (size_t in_channel = 0; in_channel < point.size(); in_channel++) { + for (size_t out_channel = 0; out_channel < point[in_channel].size(); + out_channel++) { + float gain = point[in_channel][out_channel]; + vDSP_vsma(in[in_channel] + range_start, 1, &gain, + out[out_channel] + range_start, 1, + out[out_channel] + range_start, 1, n); + } + } +#else + for (size_t out_channel = 0; out_channel < n_out; out_channel++) { for (SampleIndex i = range_start; i < range_end; i++) { out[out_channel][i] = 0.0; } @@ -296,6 +351,7 @@ namespace ear { } } } +#endif } }; } // namespace dsp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 51306c0..72e8a81 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -37,6 +37,8 @@ endif() add_ear_test("block_convolver_tests") target_sources("block_convolver_tests" PRIVATE "block_convolver_test_utils.cpp") +add_ear_test("gain_interpolator_matrix_tests") + if(EAR_USE_VDSP) add_ear_test("fft_vdsp_tests") diff --git a/tests/gain_interpolator_matrix_tests.cpp b/tests/gain_interpolator_matrix_tests.cpp new file mode 100644 index 0000000..54b1d2d --- /dev/null +++ b/tests/gain_interpolator_matrix_tests.cpp @@ -0,0 +1,99 @@ +#include +#include +#include +#include "ear/dsp/gain_interpolator.hpp" + +using namespace ear::dsp; + +namespace { + // Independent scalar reference, not shared with gain_interpolator.hpp, so + // this test is meaningful whichever branch (vDSP or scalar) that header + // compiles to. + void reference_apply_interp(const std::vector> &in, + std::vector> &out, + SampleIndex range_start, SampleIndex range_end, + SampleIndex block_start, SampleIndex start, + SampleIndex end, + const std::vector> &s, + const std::vector> &e) { + float scale = 1.0f / (end - start); + size_t n_out = s.empty() ? 0 : s[0].size(); + for (size_t oc = 0; oc < n_out; oc++) + for (SampleIndex i = range_start; i < range_end; i++) out[oc][i] = 0.0f; + + for (size_t ic = 0; ic < s.size(); ic++) + for (size_t oc = 0; oc < s[ic].size(); oc++) + for (SampleIndex i = range_start; i < range_end; i++) { + float p = (float)((block_start + i) - start) * scale; + float gain = (1.0f - p) * s[ic][oc] + p * e[ic][oc]; + out[oc][i] += in[ic][i] * gain; + } + } + + std::vector> random_matrix(size_t rows, size_t cols, + std::mt19937 &rng) { + std::uniform_real_distribution dist(-1.0f, 1.0f); + std::vector> m(rows, std::vector(cols)); + for (auto &row : m) + for (auto &v : row) v = dist(rng); + return m; + } +} // namespace + +TEST_CASE("LinearInterpMatrix apply_interp matches reference") { + std::mt19937 rng(42); + size_t n_in = 3, n_out = 5; + size_t block_size = 517; // deliberately not a power of two + + auto in = random_matrix(n_in, block_size, rng); + std::vector in_p; + for (auto &row : in) in_p.push_back(row.data()); + + auto s = random_matrix(n_in, n_out, rng); + auto e = random_matrix(n_in, n_out, rng); + + std::vector> out_actual(n_out, + std::vector(block_size)); + std::vector out_actual_p; + for (auto &row : out_actual) out_actual_p.push_back(row.data()); + + std::vector> out_expected(n_out, + std::vector(block_size)); + + LinearInterpMatrix::apply_interp(in_p.data(), out_actual_p.data(), 0, + block_size, 0, 0, block_size, s, e); + reference_apply_interp(in, out_expected, 0, block_size, 0, 0, block_size, s, + e); + + for (size_t oc = 0; oc < n_out; oc++) + for (size_t i = 0; i < block_size; i++) + CHECK(out_actual[oc][i] == + Approx(out_expected[oc][i]).margin(1e-5)); +} + +TEST_CASE("LinearInterpMatrix apply_constant matches reference") { + std::mt19937 rng(43); + size_t n_in = 4, n_out = 2; + size_t block_size = 256; + + auto in = random_matrix(n_in, block_size, rng); + std::vector in_p; + for (auto &row : in) in_p.push_back(row.data()); + + auto gains = random_matrix(n_in, n_out, rng); + + std::vector> out(n_out, std::vector(block_size)); + std::vector out_p; + for (auto &row : out) out_p.push_back(row.data()); + + LinearInterpMatrix::apply_constant(in_p.data(), out_p.data(), 0, block_size, + gains); + + for (size_t oc = 0; oc < n_out; oc++) + for (size_t i = 0; i < block_size; i++) { + float expected = 0.0f; + for (size_t ic = 0; ic < n_in; ic++) + expected += in[ic][i] * gains[ic][oc]; + CHECK(out[oc][i] == Approx(expected).margin(1e-5)); + } +}