Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion ALGORITHMS.md
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ Grouped by family (subdirectory under `lib/src/onehz/`). File paths are relative
### `respiration/`
| Function | File | Method | Citation |
|---|---|---|---|
| `rsaRespRate` | `respiration/resp_rate.dart` | respiratory sinus arrhythmia — HF spectral peak of the RR series | — |
| `rsaRespRate` | `respiration/resp_rate.dart` | respiratory sinus arrhythmia — HF spectral peak of the RR series: Lomb–Scargle on native beat times in 300 s Welch sub-windows (≥ 80 % covered by clean beats), median across them with an agreement gate; beat-rate Nyquist from the median NN | Welch 1967; Lomb 1976; Scargle 1982; Press & Rybicki 1989; DeBoer, Karemaker & Strackee 1984 |
| `riivRespRate` / `fuseRespRate` | `respiration/resp_rate.dart` | respiration-induced intensity variation, fused with the RSA estimate | Pimentel et al. (multi-grid RIIV fusion) |
| `cvhrApnea` / `cvhrApneaScreen` | `respiration/cvhr_apnea.dart` | cyclic-variation-in-HR apnea screening | — |
| `relativeOdi` | `respiration/relative_odi.dart` | ratio-of-ratios relative desaturation index — **never an absolute SpO2 claim** | — |
Expand Down
162 changes: 124 additions & 38 deletions lib/src/onehz/respiration/resp_rate.dart
Original file line number Diff line number Diff line change
Expand Up @@ -18,10 +18,12 @@
// * RIIV rides the 1 Hz green ADC, whose Nyquist caps any rate at 0.5 Hz =
// 30 br/min. We refuse to report a peak at/above that ceiling (aliasing).
// * RSA rides the TACHOGRAM, which is sampled once per BEAT, not once per
// second — so its ceiling is the beat-rate Nyquist (0.5/meanNN), ~37 br/min
// second — so its ceiling is the beat-rate Nyquist (0.5/NN), ~37 br/min
// at HR 75 but only ~25 br/min at HR 50. [rsaRespRate] computes that
// ceiling per window, caps it at [respHiHz], and withholds the whole window
// when that ceiling falls below the classic HF band (the alias would be
// ceiling from the MEDIAN beat interval (whole input and per sub-window —
// never from a time span, which a sensor gap stretches), caps it at
// [respHiHz], and withholds the whole window when that ceiling falls below
// the classic HF band (the alias would be
// indistinguishable from a genuine slower rate). A rate the window cannot
// resolve is ABSENT, never a spurious in-band peak dressed up as a normal
// number.
Expand Down Expand Up @@ -57,13 +59,17 @@ const double rsaHiHz = 0.40;

/// The highest respiratory frequency an RSA window can actually resolve (Hz).
///
/// The tachogram is sampled once per BEAT, so its Nyquist is `0.5 / meanNN`,
/// not 0.5 Hz. At HR 75 (NN 800 ms) that is 0.625 Hz; at HR 50 it is 0.417 Hz.
/// [respHiHz] caps it because nothing downstream of a 1 Hz record should claim
/// more than 30 br/min.
double rsaCeilingHz(double meanNnSec) {
if (!meanNnSec.isFinite || meanNnSec <= 0) return respHiHz;
final nyq = 0.5 / meanNnSec;
/// The tachogram is sampled once per BEAT, so its Nyquist is `0.5 / NN`, not
/// 0.5 Hz (DeBoer, Karemaker & Strackee 1984). At HR 75 (NN 800 ms) that is
/// 0.625 Hz; at HR 50 it is 0.417 Hz. [respHiHz] caps it because nothing
/// downstream of a 1 Hz record should claim more than 30 br/min.
///
/// [beatIntervalSec] is the BEAT INTERVAL — the median of the NN values, each one a
/// real interval between adjacent beats. Not span/(n−1): missing beats do not
/// lower the Nyquist, they only stretch the span.
double rsaCeilingHz(double beatIntervalSec) {
if (!beatIntervalSec.isFinite || beatIntervalSec <= 0) return respHiHz;
final nyq = 0.5 / beatIntervalSec;
return nyq < respHiHz ? nyq : respHiHz;
}

Expand All @@ -73,12 +79,20 @@ class RespEstimate {
final double? peakHz; // the spectral peak (Hz)
final double? power; // peak power (normalized)
final String source; // 'rsa' | 'riiv'
const RespEstimate(this.brpm, this.peakHz, this.power, this.source);

/// RSA only: the sub-windows the rate rests on, out of all it tried — how
/// much of the input was actually usable.
final int? usableSubwindows;
final int? subwindows;
const RespEstimate(this.brpm, this.peakHz, this.power, this.source,
{this.usableSubwindows, this.subwindows});
Map<String, dynamic> toJson() => {
'brpm': brpm == null ? null : round6(brpm!),
if (peakHz != null) 'peak_hz': round6(peakHz!),
if (power != null) 'power': round6(power!),
'source': source,
if (usableSubwindows != null) 'usable_subwindows': usableSubwindows,
if (subwindows != null) 'subwindows': subwindows,
};
}

Expand All @@ -99,6 +113,12 @@ class RespEstimate {
/// inside it.
const double rsaSegmentSec = 300;

/// Usable sub-windows at which [rsaRespRate]'s confidence stops being
/// discounted for thin evidence: ≈ 1 h of beats, since the 300 s sub-windows
/// step by 150 s — (3600 − 300) / 150 + 1 = 23. Counting them as 5-min bins
/// (12) would call ~33 minutes an hour.
const double _rsaFullWeightSubwindows = 23;

/// RSA respiratory rate from cleaned NN beat times (PRIMARY 24/7 source).
///
/// [nnMs] cleaned NN intervals (ms), [nnTimesMs] their cumulative beat times
Expand Down Expand Up @@ -133,6 +153,15 @@ const double rsaSegmentSec = 300;
/// nights, real consensus 52–85%. Uniform peaks over the ~21 br/min searchable
/// band would put ±2 br/min agreement at 19% by chance, which is what the
/// surrogates show. 50% is ~2.5× chance and sits in the measured gap.
///
/// GAPS. Both Nyquist ceilings (whole input and per sub-window) come from the
/// MEDIAN NN, the beat interval itself. Gaps belong to the completeness test:
/// a sub-window is used only when its kept beats COVER ≥ 80 % of it (Σ NN, not
/// first-to-last span). Lomb–Scargle handles the uneven sampling a dropout
/// leaves (Press & Rybicki 1989); the coverage test discards the stretches too
/// empty to spectrum. Confidence scales with how many sub-windows the rate
/// rests on, reaching full weight at [_rsaFullWeightSubwindows] (≈ 1 h), so a
/// rate read from a few minutes is not "confident".
Metric<RespEstimate> rsaRespRate(
List<double> nnMs,
List<double> nnTimesMs, {
Expand Down Expand Up @@ -170,8 +199,17 @@ Metric<RespEstimate> rsaRespRate(
// SEARCH CEILING. The window's own beat-rate Nyquist, capped at [respHiHz].
// The search used to stop at [rsaHiHz] while the guard below tested against
// [respHiHz], so the guard was unreachable and the 24 br/min cap was silent.
final meanNnSec = spanSec / (nnMs.length - 1);
final hiHz = rsaCeilingHz(meanNnSec);
//
// THE BEAT INTERVAL IS A PROPERTY OF THE BEATS, NOT OF THE CLOCK BETWEEN THE
// FIRST AND THE LAST ONE. `nnTimesMs` (correctRr) is re-anchored to the wall
// clock across sensor dropouts and advances across every rejected beat, so
// span/(n−1) is the beat interval DIVIDED BY COVERAGE: a 60 bpm night missing
// 25 % of its beats read as 45 bpm and the whole night was withheld as an
// alias. Each NN value is one real interval between adjacent beats; their
// median is the beat interval, gap-free by construction, and it is what the
// tachogram's Nyquist (0.5/NN; DeBoer et al. 1984) depends on.
final beatSec = median(nnMs)! / 1000.0;
final hiHz = rsaCeilingHz(beatSec);
// The ceiling must at least cover the classic HF band. Below that, a rate
// inside the band the literature defines would fold back down into the band
// as an alias and be indistinguishable from a genuine slower one — measured:
Expand All @@ -181,9 +219,11 @@ Metric<RespEstimate> rsaRespRate(
return Metric<RespEstimate>.absent(
tier: Tier.high,
inputs_used: inputs,
note: 'beat rate ${round6(60 / meanNnSec)} bpm resolves only to '
'${round6(hiHz * 60)} br/min — below the HF band, any peak could be '
'an alias; rate withheld',
note: 'sleeping heart rate ${round6(60 / beatSec)} bpm (median beat '
'interval ${round6(beatSec)} s) can resolve breathing from beat '
'timing only up to ${round6(hiHz * 60)} br/min, below the top of the '
'normal 9–24 br/min band, so any peak could be an alias; rate '
'withheld',
);
}

Expand All @@ -206,26 +246,40 @@ Metric<RespEstimate> rsaRespRate(
final peakHz = <double>[]; // the same peaks in Hz, index-aligned
final peakPwr = <double>[]; // their spectral power, index-aligned
var atCeiling = 0;
// The ceilings those sub-windows peaked at, br/min — their own, which
// move with heart rate; not the whole input's.
var ceilingLo = double.infinity, ceilingHi = 0.0;
var belowBand = 0;
var thin = 0;
var gappy = 0;
for (var s = tSec.first; s + segSec <= tSec.last + 1e-9; s += segSec / 2) {
final lo = _lowerBound(tSec, s);
final hi = _lowerBound(tSec, s + segSec);
final k = hi - lo;
// A sub-window has to be BOTH beat-dense and time-complete: a dropout in
// A sub-window has to be BOTH time-complete and beat-dense: a dropout in
// the middle leaves few beats spanning the full 5 minutes, and its
// periodogram is a window function, not a spectrum.
if (k < 30 || tSec[hi - 1] - tSec[lo] < segSec * 0.8) {
final segT = tSec.sublist(lo, hi);
final segNn = nnMs.sublist(lo, hi);
// TIME-COMPLETE = the kept beats COVER ≥ 80 % of the sub-window. Σ NN, not
// first-to-last span: two beats near the edges with a 3-minute hole between
// them passed the span test and handed Lomb–Scargle a window function.
// Judged first, so a sub-window a sensor gap emptied is named a gap.
final coveredSec = segNn.fold<double>(0, (a, b) => a + b) / 1000.0;
if (coveredSec < segSec * 0.8) {
gappy++;
continue;
}
if (k < 30) {
thin++;
continue;
}
final segT = tSec.sublist(lo, hi);
final segNn = nnMs.sublist(lo, hi);
final segSpan = segT.last - segT.first;
final segSpan = segT.last - segT.first; // still sets the grid resolution
// Per sub-window Nyquist: heart rate moves through the night, so the
// resolvable ceiling does too. A sub-window whose beat rate cannot cover
// the HF band is dropped for the SAME alias reason as the whole window.
final segHi = rsaCeilingHz(segSpan / (k - 1));
// the HF band is dropped for the SAME alias reason as the whole window —
// from its median beat interval, gap-free, like the whole-input check.
final segHi = rsaCeilingHz(median(segNn)! / 1000.0);
if (segHi < rsaHiHz) {
belowBand++;
continue;
Expand All @@ -246,25 +300,51 @@ Metric<RespEstimate> rsaRespRate(
// rate: reporting the edge would publish the ceiling as a measurement.
if (pk >= segHi - (segHi - rsaLoHz) / (grid - 1)) {
atCeiling++;
ceilingLo = math.min(ceilingLo, segHi * 60);
ceilingHi = math.max(ceilingHi, segHi * 60);
continue;
}
peaks.add(pk * 60.0);
peakHz.add(pk);
peakPwr.add(_powerAt(ls, pk));
}
final dropped = atCeiling + belowBand + thin + gappy;
final total = peaks.length + dropped;
if (peaks.length < 3) {
final dropped = atCeiling + belowBand + thin;
final why = dropped == 0
? 'only ${peaks.length} usable sub-windows'
: (belowBand >= atCeiling && belowBand >= thin
? '$belowBand of ${peaks.length + dropped} sub-windows had a beat '
'rate too low to cover the HF band (any peak could be an alias)'
: (atCeiling >= thin
? '$atCeiling of ${peaks.length + dropped} sub-windows peaked '
'at/above the resolvable ceiling '
'(${round6(hiHz * 60)} br/min)'
: '$thin of ${peaks.length + dropped} sub-windows were too '
'sparse or gappy to spectrum'));
// The largest counter is the reason; ties go to the earlier entry.
final reasons = <(int, String)>[
(
belowBand,
'$belowBand of $total sub-windows had a beat rate too low to cover '
'the HF band (any peak could be an alias)'
),
(
atCeiling,
atCeiling == 0
? ''
: '$atCeiling of $total sub-windows peaked at/above their '
'resolvable ceiling (${round6(ceilingLo)}'
'${ceilingHi > ceilingLo ? '–${round6(ceilingHi)}' : ''} '
'br/min)'
),
(
gappy,
Comment thread
sourcery-ai[bot] marked this conversation as resolved.
'$gappy of $total ${round6(segSec)}s sub-windows had less than 80 % '
'of their time covered by clean beats (sensor gaps or rejected '
'beats)'
),
(
thin,
'$thin of $total sub-windows had fewer than 30 beats or no spectral '
'peak'
),
];
var top = reasons.first;
for (final r in reasons.skip(1)) {
if (r.$1 > top.$1) top = r;
}
final why =
dropped == 0 ? 'only ${peaks.length} usable sub-windows' : top.$2;
return Metric<RespEstimate>.absent(
tier: Tier.high,
inputs_used: inputs,
Expand Down Expand Up @@ -303,10 +383,16 @@ Metric<RespEstimate> rsaRespRate(
final bestPeakHz = peakHz[best];
final bestPower = peakPwr[best];
// Confidence: how much of the window agrees with itself, discounted by
// artifacts. Cap below 1 (PRV ceiling).
final conf = ((1 - artifactFraction) * consensus).clamp(0.2, 0.9);
// artifacts and by how few sub-windows it rests on. Cap below 1 (PRV
// ceiling).
final conf = ((1 - artifactFraction) *
consensus *
(peaks.length / _rsaFullWeightSubwindows).clamp(0.0, 1.0))
.clamp(0.2, 0.9);
return Metric<RespEstimate>(
value: RespEstimate(brpm, bestPeakHz, bestPower, 'rsa'),
value: RespEstimate(brpm, bestPeakHz, bestPower, 'rsa',
usableSubwindows: peaks.length,
subwindows: total),
confidence: conf,
tier: Tier.high,
inputs_used: inputs,
Expand Down
Loading
Loading