diff --git a/CMakeLists.txt b/CMakeLists.txt index 2b3d765..b516d91 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,6 +39,17 @@ 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) + # 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() 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/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/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..72e8a81 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -36,6 +36,26 @@ 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") + + # 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)); + } + } +} 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)); + } +}