Initial commit

Sources/EcgSynKit/EcgSynKit.swift

@@ -0,0 +1,186 @@
+import Algorithms
+import Foundation
+import KissFFT
+
+struct Parameters {
+    /// The number of beats to simulate.
+    let numBeats: Int = 256
+
+    /// The ECG sampling frequency in Hz.
+    let sfEcg: Int = 256
+
+    /// The internal sampling frequency in Hz.
+    var sfInternal: Int = 512
+
+    /// 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
+
+    /// The ratio of power between low and high frequencies.
+    let lfhfRatio: Double = 0.5
+
+    /// The ECT amplitude in mV.
+    let amplitude: Double = 1.4
+
+    /// RNG seed value.
+    let seed: Int = 0
+
+    /// Amplitude of the noise.
+    let aNoise: Double = 0.0
+
+    /// The angle of each attractor (P, Q, R, S, T) around the limit cycle, in radians.
+    let theta: [Double] = [-60, -15, 0, 15, 90].map { $0 * .pi / 180 }
+
+    /// Widths of the attractors (P, Q, R, S, T).
+    let a: [Double] = [1.2, -5, 30, -7.5, 0.75]
+
+    /// The position of attractors (P, Q, R, S, T) above or below the z=0 plane.
+    let b: [Double] = [0.25, 0.1, 0.1, 0.1, 0.4]
+
+    /// 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
+}
+
+/// Compute fft or inverse (based on sign) of vector data, where [i] is real and [i+1] is imaginary, i starts
+/// from 0.
+func fft(data: inout [Double], isign: Int) {
+    let isign = Double(isign)
+    let n = data.count
+
+    var j = 0
+    for i in stride(from: 0, to: n, by: 2) {
+        if j > i {
+            data.swapAt(j, i)
+            data.swapAt(j + 1, i + 1)
+        }
+        var m = n / 2
+        while m >= 2, j >= m {
+            j -= m
+            m /= 2
+        }
+        j += m
+    }
+
+    var mmax = 2
+    while mmax < n {
+        let istep = 2 * mmax
+        let theta = isign * (2.0 * .pi / Double(mmax))
+        let wtemp = sin(0.5 * theta)
+        let wpr = -2.0 * wtemp * wtemp
+        let wpi = sin(theta)
+        var wr = 1.0
+        var wi = 0.0
+
+        for m in stride(from: 0, to: mmax, by: 2) {
+            for i in stride(from: m, to: n, by: istep) {
+                let j = i + mmax
+                let tempr = wr * data[j] - wi * data[j + 1]
+                let tempi = wr * data[j + 1] + wi * data[j]
+                data[j] = data[i] - tempr
+                data[j + 1] = data[i + 1] - tempi
+                data[i] += tempr
+                data[i + 1] += tempi
+            }
+            let wrtemp = wr
+            wr = wr * wpr - wi * wpi + wr
+            wi = wi * wpr + wrtemp * wpi + wi
+        }
+        mmax = istep
+    }
+}
+
+func rrprocess(params: Parameters, nrr: Int) -> [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 rrmean = 60.0 / params.hrMean
+    let rrstd = 60.0 * params.hrStd / (params.hrMean * params.hrMean)
+
+    let sf = Double(params.sfInternal)
+    let df = sf / Double(nrr)
+
+    var swc = (0 ..< nrr / 2 + 1).map {
+        let w = df * Double($0) * 2.0 * .pi
+        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 = (sf / 2.0) * sqrt(hw)
+
+        let ph = 2.0 * .pi * Double.random(in: 0.0 ..< 1.0)
+        return kiss_fft_cpx(r: sw * cos(ph), i: sw * sin(ph))
+    }
+    swc[0].i = 0.0
+    swc[nrr / 2].i = 0.0
+
+    let fft = kiss_fftr_alloc(Int32(nrr), 1, nil, nil)
+
+    let outptr = UnsafeMutablePointer<Double>.allocate(capacity: nrr)
+    kiss_fftri(fft, swc, outptr)
+
+    var rr = (0 ..< nrr).map { outptr[$0] * (1.0 / Double(nrr)) }
+    outptr.deallocate()
+
+    let xstd = stdev(rr)
+    let ratio = rrstd / xstd
+
+    for i in 0 ..< nrr {
+        rr[i] = rr[i] * ratio + rrmean
+    }
+
+    return rr
+}
+
+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))
+}
+
+func compute(params: Parameters) {
+    // adjust extrema parameters for mean heart rate
+    let hrFact = sqrt(params.hrMean / 60)
+    let hrFact2 = sqrt(hrFact)
+
+    let bi = params.b.map { $0 * hrFact }
+
+    /// XXX: Discrepancy here between Java/C and Matlab, the former uses 1.0 for ti[4] adjustment.
+    let ti = zip([hrFact2, hrFact, 1, hrFact, hrFact2], params.theta).map(*)
+
+    let ai = params.a
+
+    let x0 = SIMD3<Double>(1.0, 0.0, 0.04) // XXX: Convert to init from vector3d
+    let rseed = params.seed
+
+    // calculate time scales
+    let h = 1.0 / Double(params.sfInternal)
+    let tstep = 1.0 / Double(params.sfEcg)
+
+    // calculate length of the RR time series
+    let rrmean = (60.0 / params.hrMean)
+
+    let numRr = Int(pow(2.0, ceil(log2(Double(params.numBeats * params.sfInternal) * rrmean))))
+
+    var rr = [Double](repeating: 0.0, count: numRr)
+
+    // TODO: check sfInternal is integer multple of sfEcg
+
+    // define frequency parameters for rr process
+}