jkl1337/EcgSynKit
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update

Browse Source
This commit is contained in:
John K. Luebs
2024-11-01 23:53:57 -05:00
parent 5f3764ef63
commit 047852e9f6
28 changed files with 191 additions and 5156 deletions
Download Patch File Download Diff File Expand all files Collapse all files

228
Sources/EcgSynKit/EcgSynKit.swift
View File

@@ -1,26 +1,12 @@
import Algorithms
import ComplexModule
import Foundation
import OdeInt
import PFFFT
import PFFFTLib
import RealModule
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
public struct Parameters {
/// The ratio of power between low and high frequencies.
let lfhfRatio: Double = 0.5
@@ -28,7 +14,7 @@ struct Parameters {
let amplitude: Double = 1.4
/// RNG seed value.
let seed: Int = 0
let seed: UInt64 = 111
/// Amplitude of the noise.
let aNoise: Double = 0.0
@@ -55,8 +41,21 @@ struct Parameters {
let fhistd = 0.01
}
struct EcgDerive {
let rrpc: [Double]
public struct TimeParameters {
/// The number of beats to simulate.
let numBeats: Int = 8
/// 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 = 75.0
/// The standard deviation of the heart rate.
let hrStd: Double = 1.0
}
func stdev(_ data: [Double]) -> Double {
@@ -65,21 +64,31 @@ func stdev(_ data: [Double]) -> Double {
return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1))
}
struct RRProcess: ~Copyable {
public struct RRProcess: ~Copyable {
let nrr: Int
let sfInternal: Int
let fft: FFT<Double>
let spectrum: Buffer<Complex<Double>>
let signal: Buffer<Double>
let fft: FFT<Double>
init(nrr: Int) {
self.nrr = nrr
// mean and standard deviation of RR intervals
let rrMean: Double
let rrStd: Double
public init(params: TimeParameters) {
sfInternal = params.sfInternal
rrMean = 60.0 / params.hrMean
rrStd = 60.0 * params.hrStd / (params.hrMean * params.hrMean)
nrr = Int(pow(2.0, ceil(log2(Double(params.numBeats * sfInternal) * rrMean))))
fft = try! FFT<Double>(n: nrr)
spectrum = fft.makeSpectrumBuffer(extra: 1)
signal = fft.makeSignalBuffer()
}
func generate(params: Parameters) -> [Double] {
public func generate(params: Parameters) -> [Double] {
let w1 = 2.0 * .pi * params.flo
let w2 = 2.0 * .pi * params.fhi
let c1 = 2.0 * .pi * params.flostd
@@ -88,74 +97,159 @@ struct RRProcess: ~Copyable {
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(sfInternal)
let sf = Double(params.sfInternal)
let df = sf / Double(nrr)
let dw = (sf / Double(nrr)) * 2.0 * .pi
spectrum.withUnsafeMutableBufferPointer { swc in
for i in 0 ..< nrr / 2 + 1 {
let w = df * Double(i) * 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)
var rng = Xoshiro256Plus(seed: params.seed)
let sw = (sf / 2.0) * sqrt(hw)
let ph = 2.0 * .pi * Double.random(in: 0.0 ..< 1.0)
spectrum.mapInPlaceSwapLast { i in
let w = dw * Double(i)
swc[i].real = sw * cos(ph)
swc[i].imaginary = sw * sin(ph)
}
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)
// pack Nyquist frequency real to imaginary of DC
swc[0].imaginary = swc[nrr / 2].real
let sw = (sf / 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.withUnsafeMutableBufferPointer { outptr in
outptr.map { $0 * 1.0 / Double(nrr) }
}
var rr = signal.map { $0 * 1.0 / Double(nrr) }
let xstd = stdev(rr)
let ratio = rrstd / xstd
let ratio = rrStd / xstd
for i in 0 ..< nrr {
rr[i] = rr[i] * ratio + rrmean
rr[i] = rr[i] * ratio + rrMean
}
return rr
}
}
// func compute(params: Parameters) {
// // adjust extrema parameters for mean heart rate
// let hrFact = sqrt(params.hrMean / 60)
// let hrFact2 = sqrt(hrFact)
public struct Generator: ~Copyable {
let rrp: RRProcess
// let bi = params.b.map { $0 * hrFact }
let hrFact: Double
let hrFactSqrt: Double
// /// 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 dt: Double
// let ai = params.a
// let x0 = SIMD3<Double>(1.0, 0.0, 0.04) // XXX: Convert to init from vector3d
// let rseed = params.seed
let timeParams: TimeParameters
// // calculate time scales
// let h = 1.0 / Double(params.sfInternal)
// let tstep = 1.0 / Double(params.sfEcg)
public init(params: TimeParameters) {
rrp = RRProcess(params: params)
// // calculate length of the RR time series
// let rrmean = (60.0 / params.hrMean)
// stretching factors for intervals based on Bazett's formula
hrFact = sqrt(params.hrMean / 60)
hrFactSqrt = sqrt(hrFact)
// let numRr = Int(pow(2.0, ceil(log2(Double(params.numBeats * params.sfInternal) * rrmean))))
// init time scales
dt = 1.0 / Double(params.sfInternal)
// var rr = [Double](repeating: 0.0, count: numRr)
timeParams = params
}
// // TODO: check sfInternal is integer multple of sfEcg
public func compute(params: Parameters) -> [Double] {
let ai = params.a
let bi = params.b.map { $0 * hrFact }
// adjust extrema parameters for mean heart rate
let ti = zip([hrFactSqrt, hrFact, 1, hrFact, hrFactSqrt], params.theta).map(*)
// // define frequency parameters for rr process
// }
let sfInternal = timeParams.sfInternal
let rr = rrp.generate(params: params)
let fhi = params.fhi
let dt = self.dt
// generate piecewise RR
let rrpc = [Double](unsafeUninitializedCapacity: rr.count * 2) { rrpc, count in
var tecg = 0.0
var i = 0
var j = 0
while i < rr.count {
tecg += rr[i]
j = Int((tecg / dt).rounded(.toNearestOrEven))
for k in i ... j {
rrpc[k] = rr[i]
}
i = j + 1
}
count = i
}
print("rrpc: \(rrpc.count)")
let nt = rr.count
let x0 = SIMD3<Double>(1.0, 0.0, 0.04)
var mip = 0
var mt2 = 0.0
let ts = (0 ..< nt).map { Double($0) * dt }
print("ts: \(ts.count) \(ts.last!)")
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 ip = Int(floor(t * Double(sfInternal)))
mip = max(ip, mip)
mt2 = max(t, mt2)
//print("ip: \(ip) mip: \(mip) mt2: \(mt2)")
let w0 = 2 * .pi / rrpc[ip]
let zbase = 0.005 * sin(2 * .pi * fhi * t)
var dxdt = SIMD3<Double>()
dxdt[0] = a0 * x[0] - w0 * x[1]
dxdt[1] = a0 * x[1] + w0 * x[0]
dxdt[2] = 0.0
for i in 0 ..< ti.count {
let dt = remainder(ta - ti[i], 2 * .pi)
let dt² = dt * dt
dxdt[2] += -ai[i] * dt * exp(-0.5 * dt² / (bi[i] * bi[i]))
}
dxdt[2] += -1.0 * (x[2] - zbase)
return dxdt
}
print("mip: \(mip)")
// downsample to ECG sampling frequency
let qstep = sfInternal / timeParams.sfEcg
var zresult = stride(from: 0, to: nt, by: qstep).map { result[$0][2] }
let (zmin, zmax) = zresult.minAndMax()!
let zrange = zmax - zmin
var rng = Xoshiro256Plus(seed: params.seed + 1)
// Scale signal between -0.4 and 1.2 mV
// add uniformly distributed measurement noise
for i in 0 ..< zresult.count {
zresult[i] = 1.6 * (zresult[i] - zmin) / zrange - 0.4
zresult[i] += params.aNoise * (2.0 * rng.nextDouble() - 1.0)
}
// write zresult to text/CSV file
let filename = "ecg.csv"
try! zresult[0 ..< 800].lazy.map { String($0) }.joined(separator: "\n").write(to: URL(fileURLWithPath: filename), atomically: true, encoding: .utf8)
return zresult
}
}