ecgsyn.js/ecgsyn.cpp

#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;
}
}