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-03 17:00:23 -06:00
parent 047852e9f6
commit 583b9234fd
4 changed files with 127 additions and 128 deletions
Download Patch File Download Diff File Expand all files Collapse all files

223
Sources/EcgSynKit/EcgSynKit.swift
View File

@@ -3,31 +3,28 @@ import ComplexModule
import Foundation import Foundation
import OdeInt import OdeInt
import PFFFT import PFFFT
import PFFFTLib
import RealModule
public struct Parameters { public struct TimeParameters {
/// The ratio of power between low and high frequencies. /// The number of beats to simulate.
let lfhfRatio: Double = 0.5 let numBeats: Int = 12
/// The ECT amplitude in mV. /// The ECG sampling frequency in Hz.
let amplitude: Double = 1.4 let sfEcg: Int = 256
/// The internal sampling frequency in Hz.
let sfInternal: Int = 512
/// The mean heart rate in beats per minute.
let hrMean: Double = 70.0
/// The standard deviation of the heart rate.
let hrStd: Double = 1.0
/// RNG seed value. /// RNG seed value.
let seed: UInt64 = 111 let seed: UInt64 = 2
}
/// 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]
public struct RRParameters {
/// Mayer wave frequency in Hz. /// Mayer wave frequency in Hz.
let flo = 0.1 let flo = 0.1
@@ -39,23 +36,26 @@ public struct Parameters {
/// fhi standard deviation. /// fhi standard deviation.
let fhistd = 0.01 let fhistd = 0.01
/// The ratio of power between low and high frequencies.
let lfhfRatio: Double = 0.5
} }
public struct TimeParameters { public struct Parameters {
/// The number of beats to simulate. /// The ECG amplitude in mV.
let numBeats: Int = 8 let range: (Double, Double) = (-0.4, 1.4)
/// The ECG sampling frequency in Hz. /// Amplitude of the noise.
let sfEcg: Int = 256 let aNoise: Double = 0.0
/// The internal sampling frequency in Hz. /// The angle of each attractor (P, Q, R, S, T) around the limit cycle, in radians.
var sfInternal: Int = 512 let theta: [Double] = [-70, -15, 0, 15, 100].map { $0 * .pi / 180 }
/// The mean heart rate in beats per minute. /// The position of attractors (P, Q, R, S, T) above or below the z=0 plane.
let hrMean: Double = 75.0 let a: [Double] = [1.2, -5, 30, -7.5, 0.75]
/// The standard deviation of the heart rate. /// Widths of the attractors (P, Q, R, S, T).
let hrStd: Double = 1.0 let b: [Double] = [0.25, 0.1, 0.1, 0.1, 0.4]
} }
func stdev(_ data: [Double]) -> Double { func stdev(_ data: [Double]) -> Double {
@@ -64,31 +64,85 @@ func stdev(_ data: [Double]) -> Double {
return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1)) return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1))
} }
public struct RRProcess: ~Copyable { public struct RRSeries<T: BinaryFloatingPoint> {
public let timeParameters: TimeParameters
public let rrParamaters: RRParameters
let rng: RandomNumberGenerator
struct Segment {
let end: T
let value: T
}
let segments: [Segment]
let count: Int
public init(timeParameters: TimeParameters, rrParamaters: RRParameters, rng: RandomNumberGenerator, signal: [T]) {
self.timeParameters = timeParameters
self.rrParamaters = rrParamaters
self.rng = rng
let dt = 1.0 / T(timeParameters.sfInternal)
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 / dt).rounded(.toNearestOrEven)) + 1
}
}
segments = rrn
count = signal.count
}
public func valueAt(_ t: T) -> T {
let index = min(segments.partitioningIndex { t < $0.end }, segments.endIndex - 1)
return segments[index].value
}
}
public struct RRGenerator: ~Copyable {
let nrr: Int let nrr: Int
let sfInternal: Int
let fft: FFT<Double> let fft: FFT<Double>
let spectrum: Buffer<Complex<Double>> let spectrum: Buffer<Complex<Double>>
let signal: Buffer<Double> let signal: Buffer<Double>
var rng: RandomNumberGenerator
// mean and standard deviation of RR intervals // mean and standard deviation of RR intervals
let rrMean: Double let rrMean: Double
let rrStd: Double let rrStd: Double
let timeParameters: TimeParameters
public init(params: TimeParameters) { public init(params: TimeParameters) {
sfInternal = params.sfInternal typealias FFT = PFFFT.FFT<Double>
let sfInternal = params.sfInternal
rrMean = 60.0 / params.hrMean rrMean = 60.0 / params.hrMean
rrStd = 60.0 * params.hrStd / (params.hrMean * params.hrMean) rrStd = 60.0 * params.hrStd / (params.hrMean * params.hrMean)
nrr = Int(pow(2.0, ceil(log2(Double(params.numBeats * sfInternal) * rrMean)))) nrr = FFT.nearestValidSize(params.numBeats * sfInternal * Int(rrMean.rounded(.up)), higher: true)
fft = try! FFT(n: nrr)
fft = try! FFT<Double>(n: nrr)
spectrum = fft.makeSpectrumBuffer(extra: 1) spectrum = fft.makeSpectrumBuffer(extra: 1)
signal = fft.makeSignalBuffer() signal = fft.makeSignalBuffer()
timeParameters = params
rng = Xoshiro256Plus(seed: params.seed)
} }
public func generate(params: Parameters) -> [Double] { public mutating func generateSeries(params: RRParameters) -> RRSeries<Double> {
let rr = generateSignal(params: params)
return RRSeries(timeParameters: timeParameters, rrParamaters: params, rng: rng, signal: rr)
}
public mutating func generateSignal(params: RRParameters) -> [Double] {
let w1 = 2.0 * .pi * params.flo let w1 = 2.0 * .pi * params.flo
let w2 = 2.0 * .pi * params.fhi let w2 = 2.0 * .pi * params.fhi
let c1 = 2.0 * .pi * params.flostd let c1 = 2.0 * .pi * params.flostd
@@ -97,12 +151,10 @@ public struct RRProcess: ~Copyable {
let sig2 = 1.0 let sig2 = 1.0
let sig1 = params.lfhfRatio let sig1 = params.lfhfRatio
let sf = Double(sfInternal) let sf = Double(timeParameters.sfInternal)
let dw = (sf / Double(nrr)) * 2.0 * .pi let dw = (sf / Double(nrr)) * 2.0 * .pi
var rng = Xoshiro256Plus(seed: params.seed)
spectrum.mapInPlaceSwapLast { i in spectrum.mapInPlaceSwapLast { i in
let w = dw * Double(i) let w = dw * Double(i)
@@ -131,70 +183,28 @@ public struct RRProcess: ~Copyable {
} }
} }
public struct Generator: ~Copyable { public struct Generator {
let rrp: RRProcess
let hrFact: Double public static func generate(params: Parameters, rrSeries: RRSeries<Double>) -> [Double] {
let hrFactSqrt: Double var rng = rrSeries.rng
let sfInternal = rrSeries.timeParameters.sfInternal
let dt = 1.0 / Double(sfInternal)
let dt: Double let hrFact = 60.0 / rrSeries.timeParameters.hrMean
let hrFactSqrt = sqrt(hrFact)
let timeParams: TimeParameters
public init(params: TimeParameters) {
rrp = RRProcess(params: params)
// stretching factors for intervals based on Bazett's formula
hrFact = sqrt(params.hrMean / 60)
hrFactSqrt = sqrt(hrFact)
// init time scales
dt = 1.0 / Double(params.sfInternal)
timeParams = params
}
public func compute(params: Parameters) -> [Double] {
let ai = params.a let ai = params.a
let bi = params.b.map { $0 * hrFact } let bi = params.b.map { $0 * hrFact }
// adjust extrema parameters for mean heart rate // adjust extrema parameters for mean heart rate
let ti = zip([hrFactSqrt, hrFact, 1, hrFact, hrFactSqrt], params.theta).map(*) let ti = zip([hrFactSqrt, hrFact, 1, hrFact, hrFactSqrt], params.theta).map(*)
let sfInternal = timeParams.sfInternal let fhi = rrSeries.rrParamaters.fhi
let rr = rrp.generate(params: params) let nt = rrSeries.count
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) 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 } 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 result = SIMD3<Double>.integrate(over: ts, y0: x0, tol: 1e-6) { x, t in
let ta = atan2(x[1], x[0]) let ta = atan2(x[1], x[0])
@@ -202,54 +212,35 @@ public struct Generator: ~Copyable {
let r0 = 1.0 let r0 = 1.0
let a0 = 1.0 - sqrt(x[0] * x[0] + x[1] * x[1]) / r0 let a0 = 1.0 - sqrt(x[0] * x[0] + x[1] * x[1]) / r0
let ip = Int(floor(t * Double(sfInternal))) let w0 = 2 * .pi / rrSeries.valueAt(t)
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) let zbase = 0.005 * sin(2 * .pi * fhi * t)
var dxdt = SIMD3<Double>() var dxdt = SIMD3<Double>(a0 * x[0] - w0 * x[1], a0 * x[1] + w0 * x[0], 0.0)
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 { for i in 0 ..< ti.count {
let dt = remainder(ta - ti[i], 2 * .pi) 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] += -ai[i] * dt * exp(-0.5 * (dt * dt) / (bi[i] * bi[i]))
} }
dxdt[2] += -1.0 * (x[2] - zbase) dxdt[2] += -1.0 * (x[2] - zbase)
return dxdt return dxdt
} }
print("mip: \(mip)")
// downsample to ECG sampling frequency // downsample to ECG sampling frequency
let qstep = sfInternal / timeParams.sfEcg let qstep = sfInternal / rrSeries.timeParameters.sfEcg
var zresult = stride(from: 0, to: nt, by: qstep).map { result[$0][2] } var zresult = stride(from: 0, to: nt, by: qstep).map { result[$0][2] }
let (zmin, zmax) = zresult.minAndMax()! let (zmin, zmax) = zresult.minAndMax()!
let zrange = zmax - zmin let zrange = zmax - zmin
var rng = Xoshiro256Plus(seed: params.seed + 1)
// Scale signal between -0.4 and 1.2 mV // Scale signal between -0.4 and 1.2 mV
// add uniformly distributed measurement noise // add uniformly distributed measurement noise
for i in 0 ..< zresult.count { for i in 0 ..< zresult.count {
zresult[i] = 1.6 * (zresult[i] - zmin) / zrange - 0.4 zresult[i] = (params.range.1 - params.range.0) * (zresult[i] - zmin) / zrange - params.range.0
zresult[i] += params.aNoise * (2.0 * rng.nextDouble() - 1.0) 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 return zresult
} }
} }

9
Sources/EcgSynKit/RandomNumberGenerator.swift Normal file
View File

@@ -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
}
}

4
Sources/EcgSynKit/Xoshiro256Plus.swift
View File

@@ -56,8 +56,4 @@ struct Xoshiro256Plus: RandomNumberGenerator {
return result return result
} }
public mutating func nextDouble() -> Double {
Float64(next() >> 11) * 0x1.0p-53
}
} }

19
Tests/EcgSynKitTests/EcgSynKitTests.swift
View File

@@ -2,16 +2,19 @@ import Testing
@testable import EcgSynKit @testable import EcgSynKit
import PFFFT import PFFFT
import ComplexModule import ComplexModule
import Foundation
@Test func fftTest () { @Test func fftTest () {
let p0 = TimeParameters() let timeParameters = TimeParameters()
let c = Generator(params: p0) let rrParameters = RRParameters()
let p = Parameters()
let z = c.compute(params: p)
// print("z: \(z)") var rrg = RRGenerator(params: timeParameters)
print("PFFFT.simdArch: \(FFT<Complex<Float32>>.simdArch)")
//fft(data: &a, isign: -1) let parameters = Parameters()
//print("out: \(rr)") let ecg = Generator.generate(params: parameters, rrSeries: rrg.generateSeries(params: rrParameters))
// write ecg to file
let url = URL(fileURLWithPath: "ecg.txt")
let ecgString = ecg.map { String($0) }.joined(separator: "\n")
try! ecgString.write(to: url, atomically: true, encoding: .utf8)
} }