From f2f8b35839accd3c626ff6a5dfb15e19c216ccf0 Mon Sep 17 00:00:00 2001 From: "John K. Luebs" Date: Sun, 10 Nov 2024 18:43:50 -0600 Subject: [PATCH] save --- ecgsyn.cpp | 200 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) diff --git a/ecgsyn.cpp b/ecgsyn.cpp index 6a2074b..71dc360 100644 --- a/ecgsyn.cpp +++ b/ecgsyn.cpp @@ -1,8 +1,208 @@ +#include #include +#include +#include + +#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 inline T stdev(std::span 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 class RRSeries { + TimeParameters time_params_; + RRParameters rr_params_; + XoshiroCpp::Xoshiro256Plus rng_; + std::size_t size_; + + struct Segment { + T end; + T value; + }; + std::vector segments_; + +public: + RRSeries(TimeParameters time_params, RRParameters rr_params, + XoshiroCpp::Xoshiro256Plus rng, std::span 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(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::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 class RRGenerator { + TimeParameters time_params_; + T rr_mean_; + T rr_std_; + std::size_t nrr_; + XoshiroCpp::Xoshiro256Plus rng_; + pffft::Fft fft_; + pffft::AlignedVector signal_; + pffft::AlignedVector> 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::nearestTransformSize(time_params.num_beats * + time_params.sr_internal * + std::ceil(rr_mean_))), + fft_(pffft::Fft(nrr_)), signal_(fft_.valueVector()), + spectrum_( + pffft::AlignedVector>(fft_.getSpectrumSize() + 1)) { + if (!fft_.isValid()) { + throw FFTException(); + } + } + + std::vector generateSignal(const RRParameters ¶ms) { + 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 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 rr_gen{TimeParameters{}}; + const auto rr = rr_gen.generateSignal(RRParameters{}); + + RRSeries rr_series{TimeParameters{}, RRParameters{}, + XoshiroCpp::Xoshiro256Plus{}, rr}; + + std::printf("%f", rr_series(0.0)); + // XoshiroCpp::DoubleFromBits(); std::printf("hello world\n"); return 69; }