Optional Accelerate/vDSP FFT and gain-ramped mixing on macOS

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")

include/ear/dsp/gain_interpolator.hpp

@@ -4,6 +4,10 @@
 #include <utility>
 #include <vector>
 
+#ifdef EAR_HAS_VDSP
+#include <Accelerate/Accelerate.h>
+#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<std::vector<std::vector<float>>> {
+#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<float> &ramp_scratch(size_t n) {
+        static thread_local std::vector<float> 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<float> &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

include/ear/fft.hpp

@@ -66,4 +66,12 @@ namespace ear {
   template <typename Real>
   FFTImpl<Real> &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 <typename Real>
+  FFTImpl<Real> &get_fft_vdsp();
+#endif
+
 }  // namespace ear

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 $<TARGET_PROPERTY:xsimd,INTERFACE_INCLUDE_DIRECTORIES>)
 # Only use include directory of Eigen3 target
 # as otherwise Eigen3 needs to be installed and found when linking to a

src/fft_vdsp.cpp

@@ -0,0 +1,127 @@
+#include <Accelerate/Accelerate.h>
+#include <complex>
+#include <cstddef>
+#include <memory>
+#include <vector>
+#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<float> in_real;
+    std::vector<float> in_imag;
+    std::vector<float> out_real;
+    std::vector<float> out_imag;
+  };
+
+  class FFTPlanVdsp : public FFTPlan<float> {
+   public:
+    using Complex = typename FFTPlan<float>::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<WorkBufVdsp &>(workbuf);
+      size_t n_half = n_fft / 2;
+
+      DSPSplitComplex split_in = {wb.in_real.data(), wb.in_imag.data()};
+      vDSP_ctoz(reinterpret_cast<const DSPComplex *>(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<WorkBufVdsp &>(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<DSPComplex *>(output), 2,
+                n_half);
+    }
+
+    std::unique_ptr<FFTWorkBuf> alloc_workbuf() const override {
+      return std::unique_ptr<FFTWorkBuf>(new WorkBufVdsp(n_fft / 2));
+    }
+
+   private:
+    size_t n_fft;
+    vDSP_DFT_Setup setup_forward;
+    vDSP_DFT_Setup setup_reverse;
+  };
+
+  class FFTVdsp : public FFTImpl<float> {
+   public:
+    std::shared_ptr<FFTPlan<float>> 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<vDSP_Length>(n_fft), vDSP_DFT_FORWARD);
+      vDSP_DFT_Setup setup_reverse =
+          setup_forward ? vDSP_DFT_zrop_CreateSetup(
+                              setup_forward, static_cast<vDSP_Length>(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<float>().plan(n_fft);
+      }
+      return std::make_shared<FFTPlanVdsp>(n_fft, setup_forward, setup_reverse);
+    }
+  };
+
+  template <>
+  FFTImpl<float> EAR_EXPORT &get_fft_vdsp<float>() {
+    static FFTVdsp fft;
+    return fft;
+  }
+
+}  // namespace ear

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 $<TARGET_FILE:block_convolver_vdsp_tests>
+    WORKING_DIRECTORY "${PROJECT_BINARY_DIR}/src"
+  )
+endif()
 add_ear_test("bs2051_tests")
 add_ear_test("conversion_tests")
 add_ear_test("delay_buffer_tests")

tests/block_convolver_tests.cpp

@@ -12,7 +12,11 @@
 using namespace ear;
 using namespace ear::dsp::block_convolver;
 
+#ifdef BLOCK_CONVOLVER_TEST_VDSP
+FFTImpl<float> &fft = get_fft_vdsp<float>();
+#else
 FFTImpl<float> &fft = get_fft_kiss<float>();
+#endif
 
 /// A test of the block convolver.
 ///

tests/fft_vdsp_tests.cpp

@@ -0,0 +1,57 @@
+#include <catch2/catch.hpp>
+#include <complex>
+#include <random>
+#include <vector>
+#include "ear/fft.hpp"
+
+using namespace ear;
+
+namespace {
+  std::vector<float> random_signal(size_t n) {
+    std::mt19937 rng(1234);
+    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
+    std::vector<float> 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<float>().plan(n_fft);
+    auto plan_vdsp = get_fft_vdsp<float>().plan(n_fft);
+    auto wb_kiss = plan_kiss->alloc_workbuf();
+    auto wb_vdsp = plan_vdsp->alloc_workbuf();
+
+    std::vector<float> input = random_signal(n_fft);
+    std::vector<float> input_copy = input;
+
+    std::vector<std::complex<float>> spec_kiss(n_fft / 2 + 1);
+    std::vector<std::complex<float>> 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<float>(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<float> 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));
+    }
+  }
+}