diff --git a/src/backend/models/api/statistics.py b/src/backend/models/api/statistics.py
index 9bef81a..1bf11f9 100644
--- a/src/backend/models/api/statistics.py
+++ b/src/backend/models/api/statistics.py
@@ -50,7 +50,8 @@ class CapabilityData(BaseModel):
     pp: Optional[float] = None
     ppk: Optional[float] = None
     mean: float
-    std_dev: float
+    std_dev: float  # overall sample std dev (n-1), used for Pp/Ppk
+    std_dev_within: Optional[float] = None  # MRbar/1.128, used for Cp/Cpk
     n: int
     utl: Optional[float] = None
     ltl: Optional[float] = None
diff --git a/src/backend/services/spc_service.py b/src/backend/services/spc_service.py
index b946a70..33e0a7c 100644
--- a/src/backend/services/spc_service.py
+++ b/src/backend/services/spc_service.py
@@ -55,6 +55,27 @@ def compute_summary(
     )
 
 
+def _sigma_within(values: list[float]) -> float | None:
+    """Estimate within-subgroup sigma from the moving range of individuals.
+
+    sigma_within = MRbar / d2, with d2 = 1.128 (moving range of consecutive
+    pairs, i.e. subgroup size 2). `values` must be in chronological order
+    (the order in which measurements arrived) — the moving range is not
+    meaningful on a re-sorted series.
+
+    Returns None when fewer than 2 values are available, or when MRbar == 0
+    (no variation between consecutive points: the index would be undefined).
+    """
+    n = len(values)
+    if n < 2:
+        return None
+    moving_ranges = [abs(values[i] - values[i - 1]) for i in range(1, n)]
+    mr_bar = stats.mean(moving_ranges)
+    if mr_bar == 0:
+        return None
+    return mr_bar / 1.128
+
+
 def compute_capability(
     values: list[float],
     utl: float | None,
@@ -63,12 +84,18 @@ def compute_capability(
 ) -> CapabilityData:
     """Compute capability indices Cp, Cpk, Pp, Ppk.
 
+    `values` must be in chronological order (arrival order), since Cp/Cpk are
+    derived from the moving range between consecutive points:
+
+    - sigma_within = MRbar / 1.128 (MRbar = mean of |x_i - x_(i-1)|)
     - Cp  = (UTL - LTL) / (6 * sigma_within)
-    - Cpk = min((UTL - mean) / (3 * sigma), (mean - LTL) / (3 * sigma))
-    - Pp/Ppk: same formulas using population std dev (same data, no subgrouping)
+    - Cpk = min((UTL - mean) / (3 * sigma_within), (mean - LTL) / (3 * sigma_within))
+    - Pp/Ppk: same formulas using the overall sample std dev (n-1, `stats.stdev`)
+      of the whole series — this is the "process performance" pair and is
+      allowed to be wider than Cp/Cpk when there is drift between subgroups.
 
     Args:
-        values: Measured values.
+        values: Measured values, in chronological order.
         utl: Upper Tolerance Limit (None if not defined).
         ltl: Lower Tolerance Limit (None if not defined).
         nominal: Nominal value (None if not defined).
@@ -81,39 +108,39 @@ def compute_capability(
         return CapabilityData(
             cp=None, cpk=None, pp=None, ppk=None,
             mean=values[0] if n == 1 else 0.0,
-            std_dev=0.0, n=n,
+            std_dev=0.0, std_dev_within=None, n=n,
             utl=utl, ltl=ltl, nominal=nominal,
         )
 
     mean = stats.mean(values)
-    # Population std dev for Pp/Ppk
-    std_dev_pop = stats.pstdev(values)
-    # Sample std dev for Cp/Cpk
+    # Overall sample std dev (n-1) for Pp/Ppk
     std_dev_sample = stats.stdev(values)
+    # Within-subgroup sigma (moving range) for Cp/Cpk
+    sigma_within = _sigma_within(values)
 
     cp = cpk = pp = ppk = None
 
-    if utl is not None and ltl is not None and std_dev_sample > 0:
-        cp = round((utl - ltl) / (6 * std_dev_sample), 4)
+    if sigma_within is not None:
+        if utl is not None and ltl is not None:
+            cp = round((utl - ltl) / (6 * sigma_within), 4)
 
-    if std_dev_sample > 0:
         cpk_values = []
         if utl is not None:
-            cpk_values.append((utl - mean) / (3 * std_dev_sample))
+            cpk_values.append((utl - mean) / (3 * sigma_within))
         if ltl is not None:
-            cpk_values.append((mean - ltl) / (3 * std_dev_sample))
+            cpk_values.append((mean - ltl) / (3 * sigma_within))
         if cpk_values:
             cpk = round(min(cpk_values), 4)
 
-    if utl is not None and ltl is not None and std_dev_pop > 0:
-        pp = round((utl - ltl) / (6 * std_dev_pop), 4)
+    if std_dev_sample > 0:
+        if utl is not None and ltl is not None:
+            pp = round((utl - ltl) / (6 * std_dev_sample), 4)
 
-    if std_dev_pop > 0:
         ppk_values = []
         if utl is not None:
-            ppk_values.append((utl - mean) / (3 * std_dev_pop))
+            ppk_values.append((utl - mean) / (3 * std_dev_sample))
         if ltl is not None:
-            ppk_values.append((mean - ltl) / (3 * std_dev_pop))
+            ppk_values.append((mean - ltl) / (3 * std_dev_sample))
         if ppk_values:
             ppk = round(min(ppk_values), 4)
 
@@ -121,6 +148,7 @@ def compute_capability(
         cp=cp, cpk=cpk, pp=pp, ppk=ppk,
         mean=round(mean, 6),
         std_dev=round(std_dev_sample, 6),
+        std_dev_within=round(sigma_within, 6) if sigma_within is not None else None,
         n=n,
         utl=utl, ltl=ltl, nominal=nominal,
     )
@@ -135,10 +163,12 @@ def compute_control_chart(
     ltl: float | None,
     nominal: float | None,
 ) -> ControlChartData:
-    """Compute control chart data with UCL/LCL and out-of-control detection.
+    """Compute control chart data (individuals chart) with UCL/LCL and
+    out-of-control detection.
 
-    UCL = mean + 3*sigma
-    LCL = mean - 3*sigma
+    UCL = mean + 3*sigma_within = mean + 2.66*MRbar
+    LCL = mean - 3*sigma_within = mean - 2.66*MRbar
+    with sigma_within = MRbar/1.128 (see `_sigma_within`).
     Out-of-control: points outside UCL/LCL.
 
     Args:
@@ -160,10 +190,8 @@ def compute_control_chart(
 
     mean = stats.mean(values)
 
-    if n >= 2:
-        sigma = stats.stdev(values)
-    else:
-        sigma = 0.0
+    sigma_within = _sigma_within(values)
+    sigma = sigma_within if sigma_within is not None else 0.0
 
     ucl = mean + 3 * sigma
     lcl = mean - 3 * sigma
diff --git a/src/backend/tests/test_spc_service.py b/src/backend/tests/test_spc_service.py
new file mode 100644
index 0000000..ca9beee
--- /dev/null
+++ b/src/backend/tests/test_spc_service.py
@@ -0,0 +1,123 @@
+"""Unit tests for spc_service.py: Cp/Cpk (within, moving-range) vs Pp/Ppk (overall).
+
+Expected values are computed in each test from the formulas themselves
+(not copied from the implementation's output), so a regression in the
+formula would also break the test's own expectation derivation... except
+we compute them independently here with plain stdlib code, mirroring the
+spec in gitea:Adriano/TieMeasureFlow#19.
+"""
+import statistics as stats
+from datetime import datetime, timedelta
+
+from src.backend.services.spc_service import compute_capability, compute_control_chart
+
+
+def _expected_sigma_within(values: list[float]) -> float | None:
+    """Independent re-implementation of MRbar/1.128 for test expectations."""
+    if len(values) < 2:
+        return None
+    mr = [abs(values[i] - values[i - 1]) for i in range(1, len(values))]
+    mr_bar = stats.mean(mr)
+    if mr_bar == 0:
+        return None
+    return mr_bar / 1.128
+
+
+# 13 values around 10.00, then 12 values around 10.08: a level shift that
+# a global (overall) std dev smears, but that stays invisible to the
+# consecutive-point moving range within each stable block.
+_SHIFTED_SERIES = (
+    [10.00, 10.01, 9.99, 10.02, 9.98, 10.01, 9.99, 10.00, 10.02, 9.98, 10.01, 9.99, 10.00]
+    + [10.08, 10.09, 10.07, 10.08, 10.09, 10.07, 10.08, 10.09, 10.07, 10.08, 10.09, 10.07]
+)
+
+
+class TestComputeCapability:
+    def test_cp_uses_sigma_within_pp_uses_overall_sample_stdev(self):
+        utl, ltl = 10.2, 9.8
+        result = compute_capability(_SHIFTED_SERIES, utl, ltl, nominal=10.0)
+
+        sample_std = stats.stdev(_SHIFTED_SERIES)
+        sigma_within = _expected_sigma_within(_SHIFTED_SERIES)
+        assert sigma_within is not None
+
+        expected_cp = round((utl - ltl) / (6 * sigma_within), 4)
+        expected_pp = round((utl - ltl) / (6 * sample_std), 4)
+
+        assert result.cp == expected_cp
+        assert result.pp == expected_pp
+        # Cp (within, tight) must differ from Pp (overall, widened by the
+        # level shift) — this is exactly the drift #19 says must surface.
+        assert result.cp != result.pp
+        assert result.cp > result.pp
+        assert result.std_dev == round(sample_std, 6)
+        assert result.std_dev_within == round(sigma_within, 6)
+
+    def test_constant_series_returns_none_indices(self):
+        values = [10.0] * 10
+        result = compute_capability(values, 10.2, 9.8, 10.0)
+        assert result.cp is None
+        assert result.cpk is None
+        assert result.pp is None
+        assert result.ppk is None
+        assert result.std_dev_within is None
+
+    def test_n_less_than_2_returns_none(self):
+        result = compute_capability([10.0], 10.2, 9.8, 10.0)
+        assert result.cp is None
+        assert result.cpk is None
+        assert result.pp is None
+        assert result.ppk is None
+        assert result.n == 1
+
+    def test_one_sided_tolerance_only_cpk(self):
+        values = [10.0, 10.01, 9.99, 10.02, 9.98, 10.01, 9.99, 10.0]
+        result = compute_capability(values, utl=10.2, ltl=None, nominal=10.0)
+        assert result.cp is None  # both UTL and LTL required
+        assert result.pp is None
+        assert result.cpk is not None
+        assert result.ppk is not None
+
+
+class TestComputeControlChart:
+    def test_ucl_lcl_use_sigma_within_not_overall_stdev(self):
+        base = datetime(2026, 1, 1)
+        timestamps = [base + timedelta(minutes=i) for i in range(len(_SHIFTED_SERIES))]
+
+        result = compute_control_chart(
+            _SHIFTED_SERIES, timestamps, 10.2, None, None, 9.8, 10.0
+        )
+
+        mean = stats.mean(_SHIFTED_SERIES)
+        sigma_within = _expected_sigma_within(_SHIFTED_SERIES)
+        expected_ucl = round(mean + 3 * sigma_within, 6)
+        expected_lcl = round(mean - 3 * sigma_within, 6)
+
+        assert result.ucl == expected_ucl
+        assert result.lcl == expected_lcl
+
+        # Points beyond the within-sigma limits must be flagged: the level
+        # shift is detectable this way but was invisible with the overall
+        # sample std dev (limits too wide to ever trip).
+        assert result.out_of_control != []
+        overall_std = stats.stdev(_SHIFTED_SERIES)
+        overall_ucl = mean + 3 * overall_std
+        overall_lcl = mean - 3 * overall_std
+        overall_ooc = [
+            i for i, v in enumerate(_SHIFTED_SERIES) if v > overall_ucl or v < overall_lcl
+        ]
+        assert overall_ooc == []
+
+    def test_constant_series_zero_limits(self):
+        values = [10.0] * 6
+        timestamps = [datetime(2026, 1, 1) + timedelta(minutes=i) for i in range(6)]
+        result = compute_control_chart(values, timestamps, 10.2, None, None, 9.8, 10.0)
+        assert result.ucl == 10.0
+        assert result.lcl == 10.0
+        assert result.out_of_control == []
+
+    def test_single_value_zero_limits(self):
+        timestamps = [datetime(2026, 1, 1)]
+        result = compute_control_chart([10.0], timestamps, 10.2, None, None, 9.8, 10.0)
+        assert result.ucl == 10.0
+        assert result.lcl == 10.0