This commit is contained in:
2024-11-03 17:48:59 -06:00
parent 583b9234fd
commit 3b1a207b3f
10 changed files with 262 additions and 255 deletions

View File

@@ -4,12 +4,12 @@
import PackageDescription import PackageDescription
let package = Package( let package = Package(
name: "EcgSynKit", name: "ECGSynKit",
products: [ products: [
// Products define the executables and libraries a package produces, making them visible to other packages. // Products define the executables and libraries a package produces, making them visible to other packages.
.library( .library(
name: "EcgSynKit", name: "ECGSynKit",
targets: ["EcgSynKit"] targets: ["ECGSynKit"]
), ),
], ],
dependencies: [ dependencies: [
@@ -19,7 +19,7 @@ let package = Package(
], ],
targets: [ targets: [
.target( .target(
name: "EcgSynKit", name: "ECGSynKit",
dependencies: [ dependencies: [
.product(name: "PFFFT", package: "SwiftPFFFT"), .product(name: "PFFFT", package: "SwiftPFFFT"),
.product(name: "Algorithms", package: "swift-algorithms"), .product(name: "Algorithms", package: "swift-algorithms"),
@@ -27,8 +27,8 @@ let package = Package(
] ]
), ),
.testTarget( .testTarget(
name: "EcgSynKitTests", name: "ECGSynKitTests",
dependencies: ["EcgSynKit"] dependencies: ["ECGSynKit"]
), ),
] ]
) )

View File

@@ -0,0 +1,80 @@
import Foundation
import OdeInt
import RealModule
public enum ECGSyn {
public struct Parameters {
/// The ECG amplitude in mV.
let range: (Double, Double) = (-0.4, 1.4)
/// 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] = [-70, -15, 0, 15, 100].map { $0 * .pi / 180 }
/// The position of attractors (P, Q, R, S, T) above or below the z=0 plane.
let a: [Double] = [1.2, -5, 30, -7.5, 0.75]
/// Widths of the attractors (P, Q, R, S, T).
let b: [Double] = [0.25, 0.1, 0.1, 0.1, 0.4]
}
public static func generate(params: Parameters, rrSeries: RRSeries<Double>) -> [Double] {
var rng = rrSeries.rng
let srInternal = rrSeries.timeParameters.srInternal
let dt = 1.0 / Double(srInternal)
// adjust extrema parameters for mean heart rate
let hrFact = 60.0 / rrSeries.timeParameters.hrMean
let hrFactSqrt = sqrt(hrFact)
let ai = params.a
let bi = params.b.map { $0 * hrFact }
let ti = zip([hrFactSqrt, hrFact, 1, hrFact, hrFactSqrt], params.theta).map(*)
let fhi = rrSeries.rrParamaters.fhi
let nt = rrSeries.count
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 ECG sampling frequency
let qstep = srInternal / rrSeries.timeParameters.srEcg
var zresult = stride(from: 0, to: nt, by: qstep).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.aNoise * (2.0 * rng.nextDouble() - 1.0)
}
return zresult
}
}

View File

@@ -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 ECG sampling frequency in Hz.
let srEcg: Int = 256
/// The internal sampling frequency in Hz.
let srInternal: 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.
let seed: UInt64 = 10
}
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))
}

View File

@@ -0,0 +1,81 @@
import ComplexModule
import Foundation
import PFFFT
import RealModule
public struct RRGenerator: ~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) -> 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 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
}
}

View File

@@ -0,0 +1,44 @@
import Foundation
import RealModule
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 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
}
public func valueAt(_ t: T) -> T {
let index = min(segments.partitioningIndex { t < $0.end }, segments.endIndex - 1)
return segments[index].value
}
}

View File

@@ -1,246 +0,0 @@
import Algorithms
import ComplexModule
import Foundation
import OdeInt
import PFFFT
public struct TimeParameters {
/// The number of beats to simulate.
let numBeats: Int = 12
/// The ECG sampling frequency in Hz.
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.
let seed: UInt64 = 2
}
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
}
public struct Parameters {
/// The ECG amplitude in mV.
let range: (Double, Double) = (-0.4, 1.4)
/// 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] = [-70, -15, 0, 15, 100].map { $0 * .pi / 180 }
/// The position of attractors (P, Q, R, S, T) above or below the z=0 plane.
let a: [Double] = [1.2, -5, 30, -7.5, 0.75]
/// Widths of the attractors (P, Q, R, S, T).
let b: [Double] = [0.25, 0.1, 0.1, 0.1, 0.4]
}
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))
}
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 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 sfInternal = params.sfInternal
rrMean = 60.0 / params.hrMean
rrStd = 60.0 * params.hrStd / (params.hrMean * params.hrMean)
nrr = FFT.nearestValidSize(params.numBeats * sfInternal * 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) -> 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 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 sf = Double(timeParameters.sfInternal)
let dw = (sf / 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 = (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.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
}
}
public struct Generator {
public static func generate(params: Parameters, rrSeries: RRSeries<Double>) -> [Double] {
var rng = rrSeries.rng
let sfInternal = rrSeries.timeParameters.sfInternal
let dt = 1.0 / Double(sfInternal)
let hrFact = 60.0 / rrSeries.timeParameters.hrMean
let hrFactSqrt = sqrt(hrFact)
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(*)
let fhi = rrSeries.rrParamaters.fhi
let nt = rrSeries.count
let x0 = SIMD3<Double>(1.0, 0.0, 0.04)
let ts = (0 ..< nt).map { Double($0) * dt }
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
}
// downsample to ECG sampling frequency
let qstep = sfInternal / rrSeries.timeParameters.sfEcg
var zresult = stride(from: 0, to: nt, by: qstep).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.aNoise * (2.0 * rng.nextDouble() - 1.0)
}
return zresult
}
}

View File

@@ -1,5 +1,5 @@
import Testing import Testing
@testable import EcgSynKit @testable import ECGSynKit
import PFFFT import PFFFT
import ComplexModule import ComplexModule
import Foundation import Foundation
@@ -10,8 +10,8 @@ import Foundation
var rrg = RRGenerator(params: timeParameters) var rrg = RRGenerator(params: timeParameters)
let parameters = Parameters() let parameters = ECGSyn.Parameters()
let ecg = Generator.generate(params: parameters, rrSeries: rrg.generateSeries(params: rrParameters)) let ecg = ECGSyn.generate(params: parameters, rrSeries: rrg.generateSeries(params: rrParameters))
// write ecg to file // write ecg to file
let url = URL(fileURLWithPath: "ecg.txt") let url = URL(fileURLWithPath: "ecg.txt")
let ecgString = ecg.map { String($0) }.joined(separator: "\n") let ecgString = ecg.map { String($0) }.joined(separator: "\n")