diff --git a/lib/src/onehz/sleep/cardio_stager.dart b/lib/src/onehz/sleep/cardio_stager.dart
index 27c8cf1..9b25244 100644
--- a/lib/src/onehz/sleep/cardio_stager.dart
+++ b/lib/src/onehz/sleep/cardio_stager.dart
@@ -441,6 +441,12 @@ CardioStagerResult cardioStager(
     for (var e = 0; e < nEpoch; e++)
       if (still(e) && !rk[e].isNaN) rk[e]
   ];
+  // These three samples are fixed for the rest of the stage, so take their
+  // median+MAD once instead of re-deriving them per epoch inside the classify
+  // loop below. Same values, ~2,880 fewer sorts on a full night.
+  final rmssdScale = RobustScale.of(sleepRmssd);
+  final lfhfScale = RobustScale.of(sleepLfhf);
+  final rkScale = RobustScale.of(sleepRk);
   final hrSdSample = [for (final s in hrSd) if (s > 0) s];
   final hrSdMed = median(hrSdSample) ?? double.infinity;
   // RMSSD reference for the deep "not-elevated-HRV" gate, blended toward profile.
@@ -475,17 +481,17 @@ CardioStagerResult cardioStager(
     // large movement) AND an HR floor (HR ≥ the local p25 — REM is not the
     // quiescent cardiac trough). This recovers REM the RMSSD-only rule missed.
     final rmZ = (rmssdMed != null && !rmssd[e].isNaN && sleepRmssd.length >= 4)
-        ? robustZ(rmssd[e], sleepRmssd)
+        ? rmssdScale?.z(rmssd[e])
         : null;
     final rmssdDown = rmZ != null && rmZ < -0.4; // RMSSD notably below sleep base
     // LF/HF elevated (sympathetic shift) vs the night's sleeping LF/HF.
     final lfhfZ = (sleepLfhf.length >= 4 && !lfhf[e].isNaN)
-        ? robustZ(lfhf[e], sleepLfhf)
+        ? lfhfScale?.z(lfhf[e])
         : null;
     final lfhfHigh = lfhfZ != null && lfhfZ > _remLfhfZ;
     // R(k) burst: instantaneous-HR variability elevated vs sleeping R(k).
     final rkZ = (sleepRk.length >= 4 && !rk[e].isNaN)
-        ? robustZ(rk[e], sleepRk)
+        ? rkScale?.z(rk[e])
         : null;
     final rkBurst = rkZ != null && rkZ > _remRkZ;
     final remAutonomic = rmssdDown || lfhfHigh || rkBurst;
diff --git a/lib/src/onehz/util.dart b/lib/src/onehz/util.dart
index 377c7f0..f90c596 100644
--- a/lib/src/onehz/util.dart
+++ b/lib/src/onehz/util.dart
@@ -70,12 +70,30 @@ double? mad(List<double> xs, {bool scaled = true}) {
 /// Iglewicz–Hoaglin modified z-score of [x] against a sample, using
 /// median + MAD. Returns null if MAD is 0 (degenerate / fully-quantized) so
 /// the caller can fall back to a coarser test rather than divide by zero.
-double? robustZ(double x, List<double> sample) {
-  if (sample.length < 2) return null;
-  final m = median(sample)!;
-  final s = mad(sample);
-  if (s == null || s == 0) return null;
-  return (x - m) / s;
+double? robustZ(double x, List<double> sample) => RobustScale.of(sample)?.z(x);
+
+/// The median + MAD of a sample, computed once.
+///
+/// [robustZ] re-derives both on every call, which sorts the sample twice. When
+/// scoring many values against the SAME unchanging sample — the usual case
+/// inside a per-epoch loop — build one of these outside the loop instead.
+/// Results are identical to [robustZ], including its null semantics: absent
+/// for a sample under two values or a MAD of zero.
+class RobustScale {
+  final double med;
+  final double scale;
+
+  const RobustScale(this.med, this.scale);
+
+  static RobustScale? of(List<double> sample) {
+    if (sample.length < 2) return null;
+    final m = median(sample)!;
+    final s = mad(sample);
+    if (s == null || s == 0) return null;
+    return RobustScale(m, s);
+  }
+
+  double z(double x) => (x - med) / scale;
 }
 
 /// Ordinary z-score against a sample (mean+SD). Null if SD undefined or 0.