Update

ecgsyn.cpp

@@ -3,8 +3,8 @@
 #include <numeric>
 #include <span>
 
+#include "mini-odeint.hpp"
 #include "pffft.hpp"
-
 #include "xoshiro.hpp"
 
 #if defined(__clang__)
@@ -29,7 +29,7 @@ FORCE_INLINE constexpr auto sqr(auto &&x) noexcept(noexcept(x * x))
 }
 
 template <typename T> inline T stdev(std::span<const T> data) {
-  const auto n = T{data.size()};
+  const auto n = static_cast<T>(data.size());
   const auto mean = std::accumulate(data.begin(), data.end(), T{}) / n;
   const auto variance =
       std::accumulate(data.begin(), data.end(), T{},
@@ -38,43 +38,53 @@ template <typename T> inline T stdev(std::span<const T> data) {
   return std::sqrt(variance);
 }
 
+template <typename T>
+constexpr bool is_wasm_bindable_v =
+    std::conjunction_v<std::is_trivially_copyable<T>,
+                       std::is_standard_layout<T>>;
+
 } // namespace
 
 namespace ecgsyns {
 struct TimeParameters {
-  int num_beats = 12;
-  int sr_internal = 512;
-  int decimate_factor = 2;
-  double hr_mean = 60.0;
-  double hr_std = 1.0;
-  std::uint64_t seed = 0;
+  int num_beats{12};
+  int sr_internal{512};
+  int decimate_factor{2};
+  double hr_mean{60.0};
+  double hr_std{1.0};
+  std::uint64_t seed{8};
 };
+static_assert(is_wasm_bindable_v<TimeParameters>);
 
 struct RRParameters {
-  double flo = 0.1;
-  double flostd = 0.01;
-  double fhi = 0.25;
-  double fhistd = 0.01;
-  double lf_hf_ratio = 0.5;
+  double flo{0.1};
+  double flostd{0.01};
+  double fhi{0.25};
+  double fhistd{0.01};
+  double lf_hf_ratio{0.5};
 };
+static_assert(is_wasm_bindable_v<RRParameters>);
 
-template <typename T> class RRSeries {
+template <typename T = double> class RRSeries {
+public:
+  struct Segment {
+    T end;
+    T value;
+  };
+
+private:
   TimeParameters time_params_;
   RRParameters rr_params_;
   XoshiroCpp::Xoshiro256Plus rng_;
   std::size_t size_;
 
-  struct Segment {
-    T end;
-    T value;
-  };
   std::vector<Segment> segments_;
 
 public:
   RRSeries(TimeParameters time_params, RRParameters rr_params,
            XoshiroCpp::Xoshiro256Plus rng, std::span<const T> signal)
-      : time_params_(std::move(time_params)), rr_params_(std::move(rr_params)),
-        rng_(std::move(rng)), size_(signal.size()) {
+      : time_params_{time_params}, rr_params_{rr_params}, rng_{rng},
+        size_{signal.size()} {
 
     const auto sr = static_cast<T>(time_params_.sr_internal);
 
@@ -82,16 +92,18 @@ public:
       T tecg{};
       std::size_t i{};
 
-      while (i < size_) {
+      while (i < signal.size()) {
         tecg += signal[i];
         segments_.emplace_back(tecg, signal[i]);
-        i = static_cast<std::size_t>(std::nearbyint(tecg * sr * T{0.5}) *
-                                     T{2.0}) +
-            1;
+        i = 1 + static_cast<std::size_t>(std::nearbyint(tecg * sr * T{0.5}) *
+                                         T{2.0});
       }
     }
   }
 
+  constexpr std::size_t segments_size() const { return segments_.size(); }
+  constexpr const Segment *segments_data() const { return segments_.data(); }
+
   T operator()(T t) const noexcept {
     auto lower = std::lower_bound(
         segments_.begin(), segments_.end(), t,
@@ -102,30 +114,38 @@ public:
     return lower->value;
   }
 
-  const TimeParameters &timeParams() const noexcept { return time_params_; }
-  const RRParameters &rrParams() const noexcept { return rr_params_; }
-  const XoshiroCpp::Xoshiro256Plus &rng() const noexcept { return rng_; }
-  std::size_t size() const noexcept { return size_; }
+  [[nodiscard]] const TimeParameters &time_params() const noexcept {
+    return time_params_;
+  }
+  [[nodiscard]] const RRParameters &rr_params() const noexcept {
+    return rr_params_;
+  }
+  XoshiroCpp::Xoshiro256Plus &rng() noexcept { return rng_; }
+  [[nodiscard]] std::size_t size() const noexcept { return size_; }
+  [[nodiscard]] std::size_t output_size() const noexcept {
+    return size_ / time_params_.decimate_factor;
+  }
 };
 
 class FFTException : public std::exception {
 public:
-  FFTException() {}
-  virtual const char *what() const throw() { return "FFTException"; }
+  [[nodiscard]] const char *what() const noexcept override {
+    return "FFTException";
+  }
 };
 
-template <typename T> class RRGenerator {
+template <typename T = double> class RRGenerator {
   TimeParameters time_params_;
+  XoshiroCpp::Xoshiro256Plus rng_;
   T rr_mean_;
   T rr_std_;
   std::size_t nrr_;
-  XoshiroCpp::Xoshiro256Plus rng_;
   pffft::Fft<T> fft_;
   pffft::AlignedVector<T> signal_;
   pffft::AlignedVector<std::complex<T>> spectrum_;
 
 public:
-  RRGenerator(const TimeParameters &time_params)
+  explicit RRGenerator(const TimeParameters &time_params)
       : time_params_(time_params), rng_(time_params.seed),
         rr_mean_(60.0 / time_params.hr_mean),
         rr_std_(60.0 * time_params.hr_std / sqr(time_params.hr_mean)),
@@ -140,7 +160,20 @@ public:
     }
   }
 
-  std::vector<T> generateSignal(const RRParameters &params) {
+  constexpr std::size_t nrr() const noexcept { return nrr_; }
+
+  constexpr const TimeParameters &time_params() const noexcept {
+    return time_params_;
+  }
+
+  RRSeries<T> generate(const RRParameters &params) {
+    auto buf = std::make_unique_for_overwrite<T[]>(nrr_);
+    std::span signal{buf.get(), nrr_};
+    generate_signal(params, signal);
+    return {time_params_, params, rng_, signal};
+  }
+
+  void generate_signal(const RRParameters &params, std::span<T> output) {
     const T w1 = T{2.0 * M_PI} * params.flo;
     const T w2 = T{2.0 * M_PI} * params.fhi;
     const T c1 = T{2.0 * M_PI} * params.flostd;
@@ -151,10 +184,10 @@ public:
 
     const T sr = time_params_.sr_internal;
 
-    const T dw = (sr / T{nrr_}) * T{2.0 * M_PI};
+    const T dw = (sr / T(nrr_)) * T{2.0 * M_PI};
 
     for (std::size_t i{}; auto &p : spectrum_) {
-      const T w = dw * T{i};
+      const T w = dw * T(i);
 
       const T dw1 = w - w1;
       const T dw2 = w - w2;
@@ -176,34 +209,272 @@ public:
     fft_.inverse(spectrum_, signal_);
 
     std::ranges::transform(signal_, signal_.begin(),
-                           [this](T x) { return x * T{1.0} / T{nrr_}; });
-    const T xstd = stdev(signal_);
+                           [this](T x) { return x * T(1.0 / nrr_); });
+    const T xstd = stdev<T>(signal_);
     const T ratio = rr_std_ / xstd;
     std::vector<T> result;
     result.reserve(nrr_);
 
-    std::ranges::transform(signal_, std::back_inserter(result),
+    std::ranges::transform(signal_, output.begin(),
                            [ratio, this](T x) { return x * ratio + rr_mean_; });
-    return result;
   }
 };
 
+template <typename T = double> struct Attractor {
+  T theta;
+  T a;
+  T b;
+  T theta_rf;
+
+  constexpr static Attractor make(T degrees, T a, T b, T rf = 0.0) {
+    return {degrees * M_PI / 180.0, a, b, rf};
+  }
+};
+
+template <typename T = double, typename U = std::vector<Attractor<T>>>
+struct Parameters {
+  std::pair<const T, const T> range{-0.4, 1.2};
+  T noise_amplitude{};
+
+  U attractors;
+};
+
+template <typename T> struct Parameters<T, std::vector<Attractor<T>>> {
+  std::pair<T, T> range{-0.4, 1.2};
+  T noise_amplitude{};
+
+  std::vector<Attractor<T>> attractors{
+      Attractor<T>::make(-70.0, 1.2, 0.25, 0.25),
+      Attractor<T>::make(-15.0, -5.0, 0.1, 0.5),
+      Attractor<T>::make(0.0, 30, 0.1),
+      Attractor<T>::make(15.0, -7.5, 0.1, 0.5),
+      Attractor<T>::make(100.0, 0.75, 0.4, 0.25),
+  };
+};
+
+template <typename T, typename U = std::vector<Attractor<T>>>
+std::size_t generate(const Parameters<T, U> &params, RRSeries<T> &rr_series,
+                     std::span<T> zresult) {
+
+  struct ExPoint {
+    T ti;
+    T ai;
+    T bi;
+  };
+  auto &rng = rr_series.rng();
+
+  const auto sr_internal{rr_series.time_params().sr_internal};
+
+  const T hr_sec{rr_series.time_params().hr_mean / 60.0};
+  const T hr_fact{std::sqrt(hr_sec)};
+
+  // adjust extrema parameters for mean heart rate
+  std::vector<ExPoint> ex;
+  ex.reserve(params.attractors.size());
+  for (const auto &a : params.attractors) {
+    ex.emplace_back(a.theta * std::pow(hr_sec, a.theta_rf), a.a, a.b * hr_fact);
+  }
+
+  const T fhi{rr_series.rr_params().fhi};
+  const auto nt{rr_series.size()};
+
+  const T dt{1.0 / static_cast<T>(sr_internal)};
+
+  std::vector<T> ts;
+  ts.reserve(nt);
+  for (std::size_t i{}; i < nt; ++i) {
+    ts.emplace_back(static_cast<T>(i) * dt);
+  }
+
+  using namespace mini_odeint;
+  Vec3<T> x0{1.0, 0.0, 0.04};
+  std::vector<Vec3<T>> ys(nt);
+  explicitRungeKutta<Vec3<T>>(
+      std::span{ys}, std::span<const T>{ts}, x0, 1e-6, [&](Vec3<T> v, T t) {
+        const T ta{std::atan2(v.y, v.x)};
+
+        const T r0{1.0};
+        const T a0{T{1.0} - std::sqrt(sqr(v.x) + sqr(v.y)) / r0};
+
+        const T w0{T{2.0 * M_PI} / rr_series(t)};
+
+        const T zbase{T{0.005} * std::sin(T{2.0 * M_PI} * fhi * t)};
+
+        Vec3<T> f{a0 * v.x - w0 * v.y, a0 * v.y + w0 * v.x, 0.0};
+
+        for (const auto &e : ex) {
+          const T dt{std::remainder(ta - e.ti, T{2.0 * M_PI})};
+
+          f.z += -e.ai * dt * std::exp(T{-0.5} * sqr(dt) / sqr(e.bi));
+        }
+        f.z += T{-1.0} * (v.z + zbase);
+        return f;
+      });
+
+  // extract z and downsample to output rate
+  for (std::size_t i{}, j{}; i < nt && j < zresult.size();
+       i += rr_series.time_params().decimate_factor, ++j) {
+    zresult[j] = ys[i].z;
+  }
+
+  const auto [zmin, zmax] = std::ranges::minmax(zresult);
+  const T zrange = zmax - zmin;
+
+  // Scale signal between -0.4 and 1.2 mV
+  // add uniformly distributed measurement noise
+  std::ranges::transform(zresult, zresult.begin(), [&](T z) {
+    return (params.range.second - params.range.first) * (z - zmin) / zrange +
+           params.range.first +
+           T{2.0 * XoshiroCpp::DoubleFromBits(rng()) - 1.0} *
+               params.noise_amplitude;
+  });
+  return ys.size() / rr_series.time_params().decimate_factor;
+}
+
 } // namespace ecgsyns
 
+namespace {
+
+template <typename T> class WasmSpan {
+public:
+  T *ptr;
+  std::size_t n;
+
+public:
+  using iterator = T *;
+  constexpr std::size_t size() const noexcept { return n; }
+  constexpr iterator begin() const noexcept { return ptr; }
+  constexpr iterator end() const noexcept { return ptr + n; }
+};
+
+class WrapperBase {
+public:
+  virtual ~WrapperBase() = default;
+
+  WrapperBase(const WrapperBase &) = delete;
+  WrapperBase &operator=(const WrapperBase &) = delete;
+
+protected:
+  WrapperBase() = default;
+};
+
+template <typename T> struct Wrapper : public WrapperBase {
+  T value;
+
+  explicit Wrapper(T value) : value{std::move(value)} {}
+};
+
+struct WObject {
+  enum class type : std::uint32_t {
+    kRRGenerator = 1,
+    kRRSeries,
+  } t;
+
+  WObject(type t) : t{t} {}
+};
+
+struct WRRGenerator : WObject {
+  std::size_t nrr;
+  std::size_t output_size;
+
+  WRRGenerator(std::size_t nrr, std::size_t output_size)
+      : WObject{type::kRRGenerator}, nrr{nrr}, output_size{output_size} {}
+};
+
+struct NRRGenerator : WRRGenerator {
+  ecgsyns::RRGenerator<double> obj;
+
+  NRRGenerator(ecgsyns::RRGenerator<double> o)
+      : WRRGenerator{o.nrr(), o.nrr() / o.time_params().decimate_factor},
+        obj{std::move(o)} {}
+};
+
+struct WRRSeries : WObject {
+  std::size_t size;
+  std::size_t output_size;
+  std::size_t nsegments;
+  const ecgsyns::RRSeries<double>::Segment *segments;
+  WRRSeries(std::size_t size, std::size_t output_size, std::size_t nsegments,
+            const ecgsyns::RRSeries<double>::Segment *segments)
+      : WObject{type::kRRSeries}, size{size}, output_size{output_size},
+        nsegments{nsegments}, segments{segments} {}
+};
+
+struct NRRSeries : WRRSeries {
+  ecgsyns::RRSeries<double> obj;
+
+  NRRSeries(ecgsyns::RRSeries<double> o)
+      : WRRSeries{o.size(), o.output_size(), o.segments_size(),
+                  o.segments_data()},
+        obj{std::move(o)} {}
+};
+
+} // namespace
+
+#include <span>
+
 extern "C" {
 
-int ecgsyn() {
+using namespace ecgsyns;
+
+WRRGenerator *rr_gen_new(const TimeParameters *time_params) {
+  return new NRRGenerator{ecgsyns::RRGenerator<double>{*time_params}};
+}
+
+WRRSeries *rr_gen_new_series(WRRGenerator *generator,
+                             const RRParameters *params) {
+  return new NRRSeries{
+      static_cast<NRRGenerator *>(generator)->obj.generate(*params)};
+}
+
+void rr_gen_generate(WRRGenerator *generator, const RRParameters *params,
+                     double *output) {
+  auto &rr = static_cast<NRRGenerator *>(generator)->obj;
+  rr.generate_signal(*params, {output, rr.nrr()});
+}
+
+void destroy_obj(const WObject *o) {
+  switch (o->t) {
+  case WObject::type::kRRGenerator:
+    delete static_cast<const NRRGenerator *>(o);
+    break;
+  case WObject::type::kRRSeries:
+    delete static_cast<const NRRSeries *>(o);
+    break;
+  }
+}
+
+using WasmParameters = Parameters<double, WasmSpan<const Attractor<double>>>;
+static_assert(is_wasm_bindable_v<WasmParameters>);
+
+void ecgsyn(const WasmParameters *params, WRRSeries *rr_series,
+            double *output) {
+
+  auto &rr = static_cast<NRRSeries *>(rr_series)->obj;
+  ecgsyns::generate(*params, rr, {output, rr.output_size()});
+}
+}
+
+#if defined(ECGSYN_HOST_BUILD)
+
+#include <fstream>
+int main() {
   using namespace ecgsyns;
 
-  RRGenerator<float> rr_gen{TimeParameters{}};
-  const auto rr = rr_gen.generateSignal(RRParameters{});
+  TimeParameters time_params{};
+  time_params.num_beats = 4;
+  auto rr = rr_init(&time_params);
 
-  RRSeries<float> rr_series{TimeParameters{}, RRParameters{},
-                            XoshiroCpp::Xoshiro256Plus{}, rr};
+  const RRParameters rr_params{};
+  auto rrs = rr_generate(rr, &rr_params);
 
-  std::printf("%f", rr_series(0.0));
-  // XoshiroCpp::DoubleFromBits();
-  std::printf("hello world\n");
-  return 69;
-}
+  Parameters params;
+  std::vector<double> output(rrs->value.output_size());
+  generate(params, rrs->value, std::span(output));
+  std::ofstream f{"ecg.csv"};
+  for (const auto &x : output) {
+    f << x << '\n';
+  }
+  return 0;
 }
+#endif