diff --git a/Package.swift b/Package.swift index a9b1598..356b407 100644 --- a/Package.swift +++ b/Package.swift @@ -4,12 +4,12 @@ import PackageDescription let package = Package( - name: "EcgSynKit", + name: "ECGSynKit", products: [ // Products define the executables and libraries a package produces, making them visible to other packages. .library( - name: "EcgSynKit", - targets: ["EcgSynKit"] + name: "ECGSynKit", + targets: ["ECGSynKit"] ), ], dependencies: [ @@ -19,7 +19,7 @@ let package = Package( ], targets: [ .target( - name: "EcgSynKit", + name: "ECGSynKit", dependencies: [ .product(name: "PFFFT", package: "SwiftPFFFT"), .product(name: "Algorithms", package: "swift-algorithms"), @@ -27,8 +27,8 @@ let package = Package( ] ), .testTarget( - name: "EcgSynKitTests", - dependencies: ["EcgSynKit"] + name: "ECGSynKitTests", + dependencies: ["ECGSynKit"] ), ] ) diff --git a/Sources/ECGSynKit/ECGSyn.swift b/Sources/ECGSynKit/ECGSyn.swift new file mode 100644 index 0000000..ac8e634 --- /dev/null +++ b/Sources/ECGSynKit/ECGSyn.swift @@ -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] { + 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(1.0, 0.0, 0.04) + + let result = SIMD3.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(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 + } +} diff --git a/Sources/ECGSynKit/ECGSynKit.swift b/Sources/ECGSynKit/ECGSynKit.swift new file mode 100644 index 0000000..038ccd5 --- /dev/null +++ b/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 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)) +} diff --git a/Sources/ECGSynKit/RRGenerator.swift b/Sources/ECGSynKit/RRGenerator.swift new file mode 100644 index 0000000..b3b8d74 --- /dev/null +++ b/Sources/ECGSynKit/RRGenerator.swift @@ -0,0 +1,81 @@ +import ComplexModule +import Foundation +import PFFFT +import RealModule + +public struct RRGenerator: ~Copyable { + let nrr: Int + + let fft: FFT + let spectrum: Buffer> + let signal: Buffer + + 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 + + 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 { + 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 + } +} diff --git a/Sources/ECGSynKit/RRSeries.swift b/Sources/ECGSynKit/RRSeries.swift new file mode 100644 index 0000000..8f0083c --- /dev/null +++ b/Sources/ECGSynKit/RRSeries.swift @@ -0,0 +1,44 @@ +import Foundation +import RealModule + +public struct RRSeries { + 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 + } + +} diff --git a/Sources/EcgSynKit/RandomNumberGenerator.swift b/Sources/ECGSynKit/RandomNumberGenerator.swift similarity index 100% rename from Sources/EcgSynKit/RandomNumberGenerator.swift rename to Sources/ECGSynKit/RandomNumberGenerator.swift diff --git a/Sources/EcgSynKit/SplitMix64.swift b/Sources/ECGSynKit/SplitMix64.swift similarity index 100% rename from Sources/EcgSynKit/SplitMix64.swift rename to Sources/ECGSynKit/SplitMix64.swift diff --git a/Sources/EcgSynKit/Xoshiro256Plus.swift b/Sources/ECGSynKit/Xoshiro256Plus.swift similarity index 100% rename from Sources/EcgSynKit/Xoshiro256Plus.swift rename to Sources/ECGSynKit/Xoshiro256Plus.swift diff --git a/Sources/EcgSynKit/EcgSynKit.swift b/Sources/EcgSynKit/EcgSynKit.swift deleted file mode 100644 index 737ab3c..0000000 --- a/Sources/EcgSynKit/EcgSynKit.swift +++ /dev/null @@ -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 { - 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 - let spectrum: Buffer> - let signal: Buffer - - 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 - - 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 { - 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] { - 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(1.0, 0.0, 0.04) - - let ts = (0 ..< nt).map { Double($0) * dt } - - let result = SIMD3.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(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 - } -} diff --git a/Tests/EcgSynKitTests/EcgSynKitTests.swift b/Tests/ECGSynKitTests/ECGSynKitTests.swift similarity index 71% rename from Tests/EcgSynKitTests/EcgSynKitTests.swift rename to Tests/ECGSynKitTests/ECGSynKitTests.swift index 492210c..f4eb863 100644 --- a/Tests/EcgSynKitTests/EcgSynKitTests.swift +++ b/Tests/ECGSynKitTests/ECGSynKitTests.swift @@ -1,5 +1,5 @@ import Testing -@testable import EcgSynKit +@testable import ECGSynKit import PFFFT import ComplexModule import Foundation @@ -10,8 +10,8 @@ import Foundation var rrg = RRGenerator(params: timeParameters) - let parameters = Parameters() - let ecg = Generator.generate(params: parameters, rrSeries: rrg.generateSeries(params: rrParameters)) + let parameters = ECGSyn.Parameters() + let ecg = ECGSyn.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")