diff --git a/lib/onehz.dart b/lib/onehz.dart index 00c7a8f..bc5c7f9 100644 --- a/lib/onehz.dart +++ b/lib/onehz.dart @@ -20,6 +20,7 @@ export 'src/onehz/foundations/ppg_sqi.dart'; export 'src/onehz/foundations/baseline.dart'; export 'src/onehz/foundations/ewma_baselines.dart'; export 'src/onehz/foundations/fusion.dart'; +export 'src/onehz/foundations/gen5_ppg_hr.dart'; // Tier-1 clinical. export 'src/onehz/clinical/hrv_time.dart'; diff --git a/lib/src/onehz/foundations/gen5_ppg_hr.dart b/lib/src/onehz/foundations/gen5_ppg_hr.dart new file mode 100644 index 0000000..f55858d --- /dev/null +++ b/lib/src/onehz/foundations/gen5_ppg_hr.dart @@ -0,0 +1,102 @@ +/// Derive BPM from gen5 v26 PPG waveform bursts (24 Hz int16 samples). +/// +/// Pure, isolate-safe. Absent / degenerate input returns null — never a +/// fabricated BPM. +library; + +import 'dart:math' as math; + +/// Minimum concatenated PPG length for resting-HR ACF (~10 s @ 24 Hz). +/// Four 1 s bursts cannot resolve 40–55 bpm (needs ≥4 beats in-window). +const int kGen5PpgHrMinSamples = 240; + +/// Derive BPM from gen5 v26 PPG (24 Hz int16 samples). +/// [samples] may be one or more concatenated 24-sample bursts. +/// Returns null if too short, low variance, or no clear ACF peak in 25–230 bpm. +int? deriveHrFromGen5PpgWaveform(List samples, {double sampleHz = 24.0}) { + if (samples.length < kGen5PpgHrMinSamples || sampleHz <= 0) return null; + + final n = samples.length; + final xs = List.filled(n, 0); + var sum = 0.0; + for (var i = 0; i < n; i++) { + xs[i] = samples[i].toDouble(); + sum += xs[i]; + } + final mean = sum / n; + var energy = 0.0; + for (var i = 0; i < n; i++) { + xs[i] -= mean; + energy += xs[i] * xs[i]; + } + // Population stddev; flatline / near-constant ADC → abstain. + final std = math.sqrt(energy / n); + if (std < 1.0 || energy <= 0) return null; + + // Lag ↔ BPM: bpm = 60 * hz / lag. Search 25–230 bpm inclusive. + final minLag = math.max(1, (60.0 * sampleHz / 230.0).ceil()); + final maxLag = math.min(n - 1, (60.0 * sampleHz / 25.0).floor()); + if (minLag > maxLag) return null; + + // Normalized autocorrelation (lag-0 energy in the denominator). + final acf = List.filled(maxLag + 1, 0); + for (var lag = minLag; lag <= maxLag; lag++) { + var s = 0.0; + final lim = n - lag; + for (var i = 0; i < lim; i++) { + s += xs[i] * xs[i + lag]; + } + acf[lag] = s / energy; + } + + // Best INTERIOR local peak, with prominence against immediate neighbours. + // + // The search deliberately starts at minLag + 1 and stops at maxLag - 1. A + // periodicity claim needs a real turning point — a lag whose ACF is higher + // than the lag on EITHER side — and the boundaries of the search range + // cannot supply that evidence, because one of their neighbours lies outside + // the range and was never computed. + // + // Treating a boundary as a peak is not a harmless edge case, it is the whole + // fabrication mode. With `left = -infinity` at lag == minLag the `c <= left` + // rejection can never fire, and the prominence fallback substitutes `c - 1`, + // which makes prominence 1.0 — maximal. So ANY smoothly decaying ACF (a + // monotonic decay is exactly what you get from baseline wander, a DC drift, + // or a motion artifact — signals with no cardiac content whatsoever) was + // accepted at lag == minLag and reported as a confident 206 bpm, the top of + // the search range. Measured on this implementation before the change: + // + // baseline wander (pure sine drift, no heartbeat) -> 206 + // linear ramp (pure DC drift) -> 206 + // single step -> 206 + // + // 206 bpm from a resting wrist is not a plausible reading, and this function + // documents that it never fabricates a BPM. Requiring an interior peak makes + // all three abstain, which is the honest answer. + var bestLag = -1; + var bestScore = double.negativeInfinity; + const minAcf = 0.15; + // Need at least one lag with a computed neighbour on both sides. + if (maxLag - minLag < 2) return null; + for (var lag = minLag + 1; lag <= maxLag - 1; lag++) { + final c = acf[lag]; + if (c < minAcf) continue; + final left = acf[lag - 1]; + final right = acf[lag + 1]; + // Strict on both sides: a genuine turning point, not a shoulder. + if (c <= left || c <= right) continue; + final prominence = c - math.max(left, right); + if (prominence <= 0) continue; + // Prefer higher ACF; break ties toward the physiologically mid-range lag. + final score = c + 0.01 * prominence; + if (score > bestScore) { + bestScore = score; + bestLag = lag; + } + } + if (bestLag < 0) return null; + + final hr = (60.0 * sampleHz / bestLag).round(); + if (hr < 25 || hr > 230) return null; + return hr; +} diff --git a/test/gen5_ppg_hr_test.dart b/test/gen5_ppg_hr_test.dart new file mode 100644 index 0000000..758d3b4 --- /dev/null +++ b/test/gen5_ppg_hr_test.dart @@ -0,0 +1,165 @@ +import 'dart:math' as math; + +import 'package:openstrap_analytics/onehz.dart'; +import 'package:test/test.dart'; + +List sinePpg({ + required double bpm, + required int n, + double sampleHz = 24.0, + double amplitude = 1000, + double offset = 500, + double noiseStd = 0, + int? seed, +}) { + final rng = seed != null ? math.Random(seed) : null; + final f = bpm / 60.0; + return List.generate(n, (i) { + final t = i / sampleHz; + var v = offset + amplitude * math.sin(2 * math.pi * f * t); + if (noiseStd > 0 && rng != null) { + // Box–Muller-ish: uniform noise is enough for a regression guard. + v += (rng.nextDouble() * 2 - 1) * noiseStd; + } + return v.round(); + }); +} + +void main() { + group('deriveHrFromGen5PpgWaveform', () { + test('synthetic 72 bpm sine (10 s window) recovers ~72', () { + final samples = sinePpg(bpm: 72, n: kGen5PpgHrMinSamples); + final hr = deriveHrFromGen5PpgWaveform(samples); + expect(hr, isNotNull); + expect(hr!, closeTo(72, 3)); + }); + + test('synthetic 120 bpm sine (10 s window) recovers ~120', () { + final samples = sinePpg(bpm: 120, n: kGen5PpgHrMinSamples); + final hr = deriveHrFromGen5PpgWaveform(samples); + expect(hr, isNotNull); + expect(hr!, closeTo(120, 3)); + }); + + test('low HR 45 bpm needs a longer window (12 s)', () { + final samples = sinePpg(bpm: 45, n: 288); + final hr = deriveHrFromGen5PpgWaveform(samples); + expect(hr, isNotNull); + expect(hr!, closeTo(45, 3)); + }); + + test('noisy 72 bpm sine still recovers within tolerance', () { + final samples = sinePpg( + bpm: 72, + n: kGen5PpgHrMinSamples, + noiseStd: 80, + seed: 42, + ); + final hr = deriveHrFromGen5PpgWaveform(samples); + expect(hr, isNotNull); + expect(hr!, closeTo(72, 5)); + }); + + test('flatline abstains', () { + expect( + deriveHrFromGen5PpgWaveform(List.filled(kGen5PpgHrMinSamples, 42)), + isNull, + ); + }); + + test('too-short abstains (under 10 s)', () { + expect(deriveHrFromGen5PpgWaveform(sinePpg(bpm: 72, n: 239)), isNull); + expect(deriveHrFromGen5PpgWaveform(sinePpg(bpm: 72, n: 23)), isNull); + }); + }); + + // FABRICATION REGRESSION — "never a fabricated BPM" has to hold for signals + // with no cardiac content, not only for short or flat ones. + // + // The peak search used to consider the BOUNDARIES of its own lag range. At + // `lag == minLag` the left neighbour was `-infinity`, so the `c <= left` + // rejection could never fire, and the prominence fallback substituted + // `c - 1`, making prominence 1.0 — maximal. Any smoothly DECAYING + // autocorrelation was therefore accepted at the shortest lag and reported as + // a confident 206 bpm, the top of the 25–230 search range. + // + // A monotonically decaying ACF is exactly what baseline wander, a DC drift + // or a motion artifact produce, so the failure mode was not exotic: the + // quieter and smoother the input, the more certain the fabricated + // tachycardia. The fix requires an INTERIOR local maximum — a real turning + // point with computed neighbours on both sides. + group('never fabricates a BPM from non-cardiac input', () { + List wander() => [ + for (var i = 0; i < 480; i++) + (500 * math.sin(2 * math.pi * i / 400)).round() + ]; + List ramp() => [for (var i = 0; i < 480; i++) i * 3]; + List step() => [for (var i = 0; i < 480; i++) i < 240 ? 0 : 1000]; + List decay() => + [for (var i = 0; i < 480; i++) (1000 * math.exp(-i / 150.0)).round()]; + + test('pure baseline wander (slow sine drift, no heartbeat) abstains', () { + expect( + deriveHrFromGen5PpgWaveform(wander()), + isNull, + reason: 'used to return a confident 206 bpm', + ); + }); + + test('a pure DC ramp abstains', () { + expect(deriveHrFromGen5PpgWaveform(ramp()), isNull); + }); + + test('a single step edge abstains', () { + expect(deriveHrFromGen5PpgWaveform(step()), isNull); + }); + + test('no non-cardiac input may report the range boundary (206 bpm)', () { + final inputs = >{ + 'wander': wander(), + 'ramp': ramp(), + 'step': step(), + 'decay': decay(), + }; + inputs.forEach((name, wave) { + expect( + deriveHrFromGen5PpgWaveform(wave), + anyOf(isNull, isNot(206)), + reason: '$name pinned to the top of the search range', + ); + }); + }); + + test('white noise almost never resolves, and never at the boundary', () { + final rng = math.Random(42); + var hits = 0; + for (var t = 0; t < 400; t++) { + final noise = [ + for (var i = 0; i < 480; i++) rng.nextInt(2001) - 1000 + ]; + final hr = deriveHrFromGen5PpgWaveform(noise); + if (hr != null) { + hits++; + expect(hr, isNot(206)); + } + } + expect(hits / 400, lessThan(0.05)); + }); + }); + + group('real cardiac signals still resolve after the boundary fix', () { + test('clean sines across the physiological range stay accurate', () { + // Integer lags at 24 Hz quantize the answer, and one lag step is worth + // more at high BPM — hence a proportional tolerance rather than a fixed one. + for (final bpm in [45, 60, 75, 100, 140]) { + final hr = deriveHrFromGen5PpgWaveform(sinePpg(bpm: bpm * 1.0, n: 480)); + expect(hr, isNotNull, reason: '$bpm bpm must resolve'); + expect( + (hr! - bpm).abs(), + lessThanOrEqualTo((bpm * 0.04).ceil()), + reason: '$bpm bpm resolved as $hr', + ); + } + }); + }); +}