Initial commit

Sources/ECGSynKit/ECGSyn.swift

@@ -0,0 +1,105 @@
+import Foundation
+import OdeInt
+import RealModule
+
+public struct ECGSyn {
+    public struct Attractor {
+        /// Angle of attractor in radians
+        public let θ: Double
+        /// Position of extremum above or below the z=0 plane.
+        public let a: Double
+        /// Width of the attractor.
+        public let b: Double
+        /// Angle rate factor adjustment `θ * pow(hrMean / 60.0, θrf)`
+        public let θrf: Double
+
+        public init(θ: Double, a: Double, b: Double, θrf: Double = 0.0) {
+            self.θ = θ
+            self.a = a
+            self.b = b
+            self.θrf = θrf
+        }
+
+        public init(deg: Double, a: Double, b: Double, θrf: Double = 0.0) {
+            self.init(θ: deg * .pi / 180, a: a, b: b, θrf: θrf)
+        }
+
+        static func make(deg: Double, _ a: Double, _ b: Double, _ θrf: Double = 0.0) -> Attractor {
+            Attractor(deg: deg, a: a, b: b, θrf: θrf)
+        }
+    }
+
+    public struct Parameters {
+        /// The ECG amplitude in mV.
+        public let range: (Double, Double) = (-0.4, 1.2)
+
+        /// Amplitude of the noise.
+        public let noiseAmplitude: Double = 0.0
+
+        /// Descriptors of the extrema/attractors for the dynamical model.
+        public let attractors: [Attractor] = [
+            .make(deg: -70, 1.2, 0.25, 0.25),
+            .make(deg: -15, -5.0, 0.1, 0.5),
+            .make(deg: 0, 30, 0.1),
+            .make(deg: 15, -7.5, 0.1, 0.5),
+            .make(deg: 100, 0.75, 0.4, 0.25),
+        ]
+    }
+
+    public static func generate(params: Parameters, rrSeries: ECGSynRRSeries<Double>) -> [Double] {
+        var rng = rrSeries.rng
+        let srInternal = rrSeries.timeParameters.srInternal
+
+        let hrSec = rrSeries.timeParameters.hrMean / 60.0
+        let hrFact = sqrt(hrSec)
+
+        // adjust extrema parameters for mean heart rate
+        let ti = params.attractors.map { $0.θ * pow(hrSec, $0.θrf) }
+        let ai = params.attractors.map { $0.a }
+        let bi = params.attractors.map { $0.b * hrFact }
+
+        let fhi = rrSeries.rrParamaters.fhi
+
+        let nt = rrSeries.count
+
+        let dt = 1.0 / Double(srInternal)
+        let ts = (0 ..< nt).map { Double($0) * dt }
+        let x0 = SIMD3<Double>(1.0, 0.0, 0.04)
+
+        let result = SIMD3<Double>.integrate(over: ts, y0: x0, tol: 1e-6) { x, t in
+            let ta = atan2(x[1], x[0])
+
+            let r0 = 1.0
+            let a0 = 1.0 - sqrt(x[0] * x[0] + x[1] * x[1]) / r0
+
+            let w0 = 2 * .pi / rrSeries.valueAt(t)
+
+            let zbase = 0.005 * sin(2 * .pi * fhi * t)
+
+            var dxdt = SIMD3<Double>(a0 * x[0] - w0 * x[1], a0 * x[1] + w0 * x[0], 0.0)
+
+            for i in 0 ..< ti.count {
+                let dt = remainder(ta - ti[i], 2 * .pi)
+
+                dxdt[2] += -ai[i] * dt * exp(-0.5 * (dt * dt) / (bi[i] * bi[i]))
+            }
+            dxdt[2] += -1.0 * (x[2] - zbase)
+
+            return dxdt
+        }
+
+        // extract z and downsample to output sampling frequency
+        var zresult = stride(from: 0, to: nt, by: rrSeries.timeParameters.decimateFactor).map { result[$0][2] }
+
+        let (zmin, zmax) = zresult.minAndMax()!
+        let zrange = zmax - zmin
+
+        // Scale signal between -0.4 and 1.2 mV
+        // add uniformly distributed measurement noise
+        for i in 0 ..< zresult.count {
+            zresult[i] = (params.range.1 - params.range.0) * (zresult[i] - zmin) / zrange + params.range.0
+            zresult[i] += params.noiseAmplitude * (2.0 * rng.nextDouble() - 1.0)
+        }
+        return zresult
+    }
+}

Sources/ECGSynKit/ECGSynKit.swift

@@ -0,0 +1,48 @@
+import Algorithms
+import ComplexModule
+import RealModule
+import Foundation
+import PFFFT
+
+public struct TimeParameters {
+    /// The number of beats to simulate.
+    let numBeats: Int = 12
+
+    /// The internal sampling frequency in Hz.
+    let srInternal: Int = 512
+
+    /// Output decimation factor
+    let decimateFactor: Int = 2
+
+    /// The mean heart rate in beats per minute.
+    let hrMean: Double = 60.0
+
+    /// The standard deviation of the heart rate.
+    let hrStd: Double = 1.0
+
+    /// RNG seed value.
+    let seed: UInt64 = 8
+}
+
+public struct RRParameters {
+    /// Mayer wave frequency in Hz.
+    let flo = 0.1
+
+    /// flo standard deviation.
+    let flostd = 0.01
+
+    /// Respiratory rate frequency in Hz.
+    let fhi = 0.25
+
+    /// fhi standard deviation.
+    let fhistd = 0.01
+
+    /// The ratio of power between low and high frequencies.
+    let lfhfRatio: Double = 0.5
+}
+
+func stdev(_ data: [Double]) -> Double {
+    let n = Double(data.count)
+    let mean = data.reduce(0.0, +) / n
+    return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1))
+}

Sources/ECGSynKit/ECGSynRRGenerator.swift

@@ -0,0 +1,81 @@
+import ComplexModule
+import Foundation
+import PFFFT
+import RealModule
+
+public struct ECGSynRRGenerator: ~Copyable {
+    let nrr: Int
+
+    let fft: FFT<Double>
+    let spectrum: Buffer<Complex<Double>>
+    let signal: Buffer<Double>
+
+    var rng: RandomNumberGenerator
+
+    // mean and standard deviation of RR intervals
+    let rrMean: Double
+    let rrStd: Double
+
+    let timeParameters: TimeParameters
+
+    public init(params: TimeParameters) {
+        typealias FFT = PFFFT.FFT<Double>
+
+        let sr = params.srInternal
+        rrMean = 60.0 / params.hrMean
+        rrStd = 60.0 * params.hrStd / (params.hrMean * params.hrMean)
+
+        nrr = FFT.nearestValidSize(params.numBeats * sr * Int(rrMean.rounded(.up)), higher: true)
+        fft = try! FFT(n: nrr)
+        spectrum = fft.makeSpectrumBuffer(extra: 1)
+        signal = fft.makeSignalBuffer()
+
+        timeParameters = params
+        rng = Xoshiro256Plus(seed: params.seed)
+    }
+
+    public mutating func generateSeries(params: RRParameters) -> ECGSynRRSeries<Double> {
+        let rr = generateSignal(params: params)
+        return ECGSynRRSeries(timeParameters: timeParameters, rrParamaters: params, rng: rng, signal: rr)
+    }
+
+    public mutating func generateSignal(params: RRParameters) -> [Double] {
+        let w1 = 2.0 * .pi * params.flo
+        let w2 = 2.0 * .pi * params.fhi
+        let c1 = 2.0 * .pi * params.flostd
+        let c2 = 2.0 * .pi * params.fhistd
+
+        let sig2 = 1.0
+        let sig1 = params.lfhfRatio
+
+        let sr = Double(timeParameters.srInternal)
+
+        let dw = (sr / Double(nrr)) * 2.0 * .pi
+
+        spectrum.mapInPlaceSwapLast { i in
+            let w = dw * Double(i)
+
+            let dw1 = w - w1
+            let dw2 = w - w2
+            let hw = sig1 * exp(-dw1 * dw1 / (2.0 * c1 * c1)) / sqrt(2.0 * .pi * c1 * c1)
+                + sig2 * exp(-dw2 * dw2 / (2.0 * c2 * c2)) / sqrt(2.0 * .pi * c2 * c2)
+
+            let sw = (sr / 2.0) * sqrt(hw)
+            let ph = 2.0 * .pi * rng.nextDouble()
+
+            return Complex(length: sw, phase: ph)
+        }
+
+        fft.inverse(spectrum: spectrum, signal: signal)
+
+        var rr = signal.map { $0 * 1.0 / Double(nrr) }
+
+        let xstd = stdev(rr)
+        let ratio = rrStd / xstd
+
+        for i in 0 ..< nrr {
+            rr[i] = rr[i] * ratio + rrMean
+        }
+        return rr
+    }
+}

Sources/ECGSynKit/ECGSynRRSeries.swift

@@ -0,0 +1,46 @@
+import Foundation
+import RealModule
+import Algorithms
+
+public struct ECGSynRRSeries<T: BinaryFloatingPoint> {
+    let timeParameters: TimeParameters
+    let rrParamaters: RRParameters
+    let rng: RandomNumberGenerator
+    let count: Int
+
+    struct Segment {
+        let end: T
+        let value: T
+    }
+    let segments: [Segment]
+
+    public init(timeParameters: TimeParameters, rrParamaters: RRParameters, rng: RandomNumberGenerator, signal: [T]) {
+        self.timeParameters = timeParameters
+        self.rrParamaters = rrParamaters
+        self.rng = rng
+
+        let sr = T(timeParameters.srInternal)
+
+        var rrn = [Segment]()
+        // generate piecewise RR time series
+        do {
+            var tecg = T.zero
+            var i = 0
+            while i < signal.count {
+                tecg += signal[i]
+                rrn.append(Segment(end: tecg, value: signal[i]))
+                i = Int((tecg * sr).rounded(.toNearestOrEven)) + 1
+            }
+        }
+
+        segments = rrn
+        count = signal.count
+    }
+
+    @inline(__always)
+    public func valueAt(_ t: T) -> T {
+        let index = min(segments.partitioningIndex { t < $0.end }, segments.endIndex - 1)
+        return segments[index].value
+    }
+
+}

Sources/ECGSynKit/RandomNumberGenerator.swift

@@ -0,0 +1,9 @@
+extension RandomNumberGenerator {
+    mutating func nextDouble() -> Double {
+        Double(next() >> 11) * 0x1.0p-53
+    }
+
+    mutating func nextFloat() -> Float {
+        Float(next() >> 40) * 0x1.0p-24
+    }
+}

Sources/ECGSynKit/SplitMix64.swift

@@ -0,0 +1,21 @@
+
+struct SplitMix64 : RandomNumberGenerator {
+    public typealias State = UInt64
+    public private(set) var state: State
+
+    init(state: UInt64) {
+        self.state = state
+    }
+
+    public mutating func next() -> UInt64 {
+        state &+= 0x9E3779B97F4A7C15
+        var z = state
+        z = (z ^ (z >> 30)) &* 0xBF58476D1CE4E5B9
+        z = (z ^ (z >> 27)) &* 0x94D049BB133111EB
+        return z ^ (z >> 31)
+    }
+
+    public mutating func nextDouble() -> Double {
+        Double(next() >> 11) * 0x1.0p-53
+    }
+}

Sources/ECGSynKit/Xoshiro256Plus.swift

@@ -0,0 +1,59 @@
+import Foundation
+
+public enum Xoshiro256: Equatable {
+    public typealias State = (UInt64, UInt64, UInt64, UInt64)
+
+    internal static var invalidState: State { (0, 0, 0, 0) }
+
+    internal static func isValid(state: State) -> Bool {
+        state != invalidState
+    }
+}
+
+@inlinable
+@inline(__always)
+internal func rotl(_ x: UInt64, _ k: UInt64) -> UInt64 {
+    (x << k) | (x >> (64 &- k))
+}
+
+struct Xoshiro256Plus: RandomNumberGenerator {
+    public typealias State = Xoshiro256.State
+
+    private var state: State
+
+    public init() {
+        var generator = SystemRandomNumberGenerator()
+        self.init(seed: generator.next())
+    }
+
+    public init(seed: UInt64) {
+        var generator = SplitMix64(state: seed)
+        var state = Xoshiro256.invalidState
+
+        repeat {
+            state = (generator.next(), generator.next(), generator.next(), generator.next())
+        } while !Xoshiro256.isValid(state: state)
+
+        self.init(state: state)
+    }
+
+    public init(state: State) {
+        precondition(Xoshiro256.isValid(state: state), "The state must not be zero")
+        self.state = state
+    }
+
+    public mutating func next() -> UInt64 {
+        let result = state.0 &+ state.3
+        let t = state.1 << 17
+
+        state.2 ^= state.0
+        state.3 ^= state.1
+        state.1 ^= state.2
+        state.0 ^= state.3
+
+        state.2 ^= t
+        state.3 = rotl(state.3, 45)
+
+        return result
+    }
+}