diff --git a/Sources/EcgSynKit/EcgSynKit.swift b/Sources/EcgSynKit/EcgSynKit.swift index c1834f4..737ab3c 100644 --- a/Sources/EcgSynKit/EcgSynKit.swift +++ b/Sources/EcgSynKit/EcgSynKit.swift @@ -3,31 +3,28 @@ import ComplexModule import Foundation import OdeInt import PFFFT -import PFFFTLib -import RealModule -public struct Parameters { - /// The ratio of power between low and high frequencies. - let lfhfRatio: Double = 0.5 +public struct TimeParameters { + /// The number of beats to simulate. + let numBeats: Int = 12 - /// The ECT amplitude in mV. - let amplitude: Double = 1.4 + /// 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 = 111 - - /// 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] + let seed: UInt64 = 2 +} +public struct RRParameters { /// Mayer wave frequency in Hz. let flo = 0.1 @@ -39,23 +36,26 @@ public struct Parameters { /// fhi standard deviation. let fhistd = 0.01 + + /// The ratio of power between low and high frequencies. + let lfhfRatio: Double = 0.5 } -public struct TimeParameters { - /// The number of beats to simulate. - let numBeats: Int = 8 +public struct Parameters { + /// The ECG amplitude in mV. + let range: (Double, Double) = (-0.4, 1.4) - /// The ECG sampling frequency in Hz. - let sfEcg: Int = 256 + /// Amplitude of the noise. + let aNoise: Double = 0.0 - /// The internal sampling frequency in Hz. - var sfInternal: Int = 512 + /// 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 mean heart rate in beats per minute. - let hrMean: Double = 75.0 + /// 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] - /// The standard deviation of the heart rate. - let hrStd: Double = 1.0 + /// 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 { @@ -64,31 +64,85 @@ func stdev(_ data: [Double]) -> Double { return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1)) } -public struct RRProcess: ~Copyable { +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 sfInternal: 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) { - sfInternal = params.sfInternal + typealias FFT = PFFFT.FFT + + let 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(n: nrr) + 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 func generate(params: Parameters) -> [Double] { + 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 @@ -97,12 +151,10 @@ public struct RRProcess: ~Copyable { let sig2 = 1.0 let sig1 = params.lfhfRatio - let sf = Double(sfInternal) + let sf = Double(timeParameters.sfInternal) let dw = (sf / Double(nrr)) * 2.0 * .pi - var rng = Xoshiro256Plus(seed: params.seed) - spectrum.mapInPlaceSwapLast { i in let w = dw * Double(i) @@ -131,70 +183,28 @@ public struct RRProcess: ~Copyable { } } -public struct Generator: ~Copyable { - let rrp: RRProcess +public struct Generator { - let hrFact: Double - let hrFactSqrt: Double + 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 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 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 sfInternal = timeParams.sfInternal + let fhi = rrSeries.rrParamaters.fhi - 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 nt = rrSeries.count let x0 = SIMD3(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.integrate(over: ts, y0: x0, tol: 1e-6) { x, t in let ta = atan2(x[1], x[0]) @@ -202,54 +212,35 @@ public struct Generator: ~Copyable { 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 w0 = 2 * .pi / rrSeries.valueAt(t) let zbase = 0.005 * sin(2 * .pi * fhi * t) - var dxdt = SIMD3() - - dxdt[0] = a0 * x[0] - w0 * x[1] - dxdt[1] = a0 * x[1] + w0 * x[0] - - dxdt[2] = 0.0 + 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) - 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) return dxdt } - print("mip: \(mip)") - // 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] } 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.range.1 - params.range.0) * (zresult[i] - zmin) / zrange - params.range.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 } } diff --git a/Sources/EcgSynKit/RandomNumberGenerator.swift b/Sources/EcgSynKit/RandomNumberGenerator.swift new file mode 100644 index 0000000..7b68d9d --- /dev/null +++ b/Sources/EcgSynKit/RandomNumberGenerator.swift @@ -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 + } +} diff --git a/Sources/EcgSynKit/Xoshiro256Plus.swift b/Sources/EcgSynKit/Xoshiro256Plus.swift index 50998f2..6c1dfbd 100644 --- a/Sources/EcgSynKit/Xoshiro256Plus.swift +++ b/Sources/EcgSynKit/Xoshiro256Plus.swift @@ -56,8 +56,4 @@ struct Xoshiro256Plus: RandomNumberGenerator { return result } - - public mutating func nextDouble() -> Double { - Float64(next() >> 11) * 0x1.0p-53 - } } diff --git a/Tests/EcgSynKitTests/EcgSynKitTests.swift b/Tests/EcgSynKitTests/EcgSynKitTests.swift index 7f85da9..492210c 100644 --- a/Tests/EcgSynKitTests/EcgSynKitTests.swift +++ b/Tests/EcgSynKitTests/EcgSynKitTests.swift @@ -2,16 +2,19 @@ import Testing @testable import EcgSynKit import PFFFT import ComplexModule +import Foundation @Test func fftTest () { - let p0 = TimeParameters() - let c = Generator(params: p0) - let p = Parameters() - let z = c.compute(params: p) + let timeParameters = TimeParameters() + let rrParameters = RRParameters() - // print("z: \(z)") - print("PFFFT.simdArch: \(FFT>.simdArch)") - //fft(data: &a, isign: -1) - //print("out: \(rr)") + var rrg = RRGenerator(params: timeParameters) + + let parameters = Parameters() + 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) }