save

ecgsyn.cpp

@@ -1,8 +1,208 @@
+#include <cstdint>
 #include <cstdio>
+#include <numeric>
+#include <span>
+
+#include "pffft.hpp"
+
+#include "xoshiro.hpp"
+
+#if defined(__clang__)
+#define FORCE_INLINE [[gnu::always_inline]] [[gnu::gnu_inline]] extern inline
+
+#elif defined(__GNUC__)
+#define FORCE_INLINE [[gnu::always_inline]] inline
+
+#elif defined(_MSC_VER)
+#pragma warning(error : 4714)
+#define FORCE_INLINE __forceinline
+
+#else
+#define FORCE_INLINE forceinline
+#endif
+
+namespace {
+
+FORCE_INLINE constexpr auto sqr(auto &&x) noexcept(noexcept(x * x))
+    -> decltype(x * x) {
+  return x * x;
+}
+
+template <typename T> inline T stdev(std::span<const T> data) {
+  const auto n = T{data.size()};
+  const auto mean = std::accumulate(data.begin(), data.end(), T{}) / n;
+  const auto variance =
+      std::accumulate(data.begin(), data.end(), T{},
+                      [mean](T acc, T x) { return acc + sqr(x - mean); }) /
+      (n - 1);
+  return std::sqrt(variance);
+}
+
+} // 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;
+};
+
+struct RRParameters {
+  double flo = 0.1;
+  double flostd = 0.01;
+  double fhi = 0.25;
+  double fhistd = 0.01;
+  double lf_hf_ratio = 0.5;
+};
+
+template <typename T> class RRSeries {
+  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()) {
+
+    const auto sr = static_cast<T>(time_params_.sr_internal);
+
+    {
+      T tecg{};
+      std::size_t i{};
+
+      while (i < 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;
+      }
+    }
+  }
+
+  T operator()(T t) const noexcept {
+    auto lower = std::lower_bound(
+        segments_.begin(), segments_.end(), t,
+        [](const auto &seg, const auto t) { return seg.end < t; });
+    if (lower == segments_.end()) {
+      lower = std::prev(lower);
+    }
+    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_; }
+};
+
+class FFTException : public std::exception {
+public:
+  FFTException() {}
+  virtual const char *what() const throw() { return "FFTException"; }
+};
+
+template <typename T> class RRGenerator {
+  TimeParameters time_params_;
+  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)
+      : 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)),
+        nrr_(pffft::Fft<T>::nearestTransformSize(time_params.num_beats *
+                                                 time_params.sr_internal *
+                                                 std::ceil(rr_mean_))),
+        fft_(pffft::Fft<T>(nrr_)), signal_(fft_.valueVector()),
+        spectrum_(
+            pffft::AlignedVector<std::complex<T>>(fft_.getSpectrumSize() + 1)) {
+    if (!fft_.isValid()) {
+      throw FFTException();
+    }
+  }
+
+  std::vector<T> generateSignal(const RRParameters &params) {
+    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;
+    const T c2 = T{2.0 * M_PI} * params.fhistd;
+
+    const T sig2 = T{1};
+    const T sig1 = params.lf_hf_ratio;
+
+    const T sr = time_params_.sr_internal;
+
+    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 dw1 = w - w1;
+      const T dw2 = w - w2;
+
+      const T hw = sig1 * std::exp(-sqr(dw1) / (T{2} * sqr(c1))) /
+                       std::sqrt(T{2.0 * M_PI} * sqr(c1)) +
+                   sig2 * std::exp(-sqr(dw2) / (T{2} * sqr(c2))) /
+                       std::sqrt(T{2.0 * M_PI} * sqr(c2));
+
+      const T sw = (sr / T{2}) * std::sqrt(hw);
+
+      const T ph = 2.0 * M_PI * XoshiroCpp::DoubleFromBits(rng_());
+
+      p = std::polar(sw, ph);
+      ++i;
+    }
+    spectrum_.front().imag(spectrum_.back().real());
+
+    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_);
+    const T ratio = rr_std_ / xstd;
+    std::vector<T> result;
+    result.reserve(nrr_);
+
+    std::ranges::transform(signal_, std::back_inserter(result),
+                           [ratio, this](T x) { return x * ratio + rr_mean_; });
+    return result;
+  }
+};
+
+} // namespace ecgsyns
 
 extern "C" {
 
 int ecgsyn() {
+  using namespace ecgsyns;
+
+  RRGenerator<float> rr_gen{TimeParameters{}};
+  const auto rr = rr_gen.generateSignal(RRParameters{});
+
+  RRSeries<float> rr_series{TimeParameters{}, RRParameters{},
+                            XoshiroCpp::Xoshiro256Plus{}, rr};
+
+  std::printf("%f", rr_series(0.0));
+  // XoshiroCpp::DoubleFromBits();
   std::printf("hello world\n");
   return 69;
 }