diff --git a/audio_field/program_reference.py b/audio_field/program_reference.py
new file mode 100644
index 0000000..f36b858
--- /dev/null
+++ b/audio_field/program_reference.py
@@ -0,0 +1,105 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+from scipy import signal
+
+
+@dataclass(frozen=True)
+class ReferenceSuppressionReport:
+    channel: int
+    input_rms: float
+    predicted_rms: float
+    output_rms: float
+    reduction_db: float
+    phase_mapping: tuple[tuple[float, float, float], ...]
+
+
+def suppress_program_reference(
+    field: np.ndarray,
+    reference: np.ndarray,
+    sample_rate: int,
+    channels: list[int] | tuple[int, ...],
+    *,
+    nperseg: int = 2048,
+    noverlap: int = 1536,
+    regularization: float = 1e-6,
+    subtraction_strength: float = 0.85,
+) -> tuple[np.ndarray, list[ReferenceSuppressionReport]]:
+    x = np.asarray(field, dtype=np.float32)
+    ref = np.asarray(reference, dtype=np.float32).reshape(-1)
+    if x.ndim != 2:
+        raise ValueError("field must be frames-by-channels")
+    ref = pad_or_trim(ref, len(x))
+    cleaned = x.copy()
+    reports = []
+    for channel in channels:
+        ch = int(channel)
+        if ch < 0 or ch >= x.shape[1]:
+            continue
+        transfer, freqs = estimate_transfer(reference=ref, observed=x[:, ch], sample_rate=sample_rate, nperseg=nperseg, noverlap=noverlap, regularization=regularization)
+        predicted = synthesize_reference(reference=ref, transfer=transfer, sample_rate=sample_rate, nperseg=nperseg, noverlap=noverlap)
+        predicted = pad_or_trim(predicted, len(x))
+        before = x[:, ch]
+        after = before - float(subtraction_strength) * predicted
+        cleaned[:, ch] = after.astype(np.float32)
+        before_rms = rms(before)
+        after_rms = rms(after)
+        reports.append(
+            ReferenceSuppressionReport(
+                channel=ch,
+                input_rms=before_rms,
+                predicted_rms=rms(predicted),
+                output_rms=after_rms,
+                reduction_db=20.0 * np.log10(max(after_rms, 1e-12) / max(before_rms, 1e-12)),
+                phase_mapping=phase_mapping(freqs, transfer),
+            )
+        )
+    return cleaned.astype(np.float32), reports
+
+
+def estimate_transfer(*, reference: np.ndarray, observed: np.ndarray, sample_rate: int, nperseg: int, noverlap: int, regularization: float):
+    freqs, _, ref_stft = signal.stft(reference, fs=sample_rate, nperseg=nperseg, noverlap=noverlap, boundary="zeros", padded=True)
+    _, _, obs_stft = signal.stft(pad_or_trim(observed, len(reference)), fs=sample_rate, nperseg=nperseg, noverlap=noverlap, boundary="zeros", padded=True)
+    numerator = np.sum(obs_stft * np.conj(ref_stft), axis=1)
+    denominator = np.sum(np.abs(ref_stft) ** 2, axis=1) + float(regularization)
+    transfer = numerator / denominator
+    return transfer.astype(np.complex64), freqs.astype(np.float32)
+
+
+def synthesize_reference(*, reference: np.ndarray, transfer: np.ndarray, sample_rate: int, nperseg: int, noverlap: int):
+    _, _, ref_stft = signal.stft(reference, fs=sample_rate, nperseg=nperseg, noverlap=noverlap, boundary="zeros", padded=True)
+    predicted_stft = ref_stft * transfer[:, None]
+    _, predicted = signal.istft(predicted_stft, fs=sample_rate, nperseg=nperseg, noverlap=noverlap, input_onesided=True, boundary=True)
+    return predicted.astype(np.float32)
+
+
+def phase_mapping(freqs: np.ndarray, transfer: np.ndarray, *, min_hz: float = 80.0, max_hz: float = 16000.0, bins: int = 24):
+    valid = (freqs >= min_hz) & (freqs <= max_hz) & (np.abs(transfer) > 1e-9)
+    if not np.any(valid):
+        return tuple()
+    selected_freqs = np.geomspace(max(min_hz, float(freqs[valid][0])), min(max_hz, float(freqs[valid][-1])), bins)
+    rows = []
+    for frequency in selected_freqs:
+        index = int(np.argmin(np.abs(freqs - frequency)))
+        magnitude = float(abs(transfer[index]))
+        phase = float(np.angle(transfer[index]))
+        rows.append((float(freqs[index]), phase, magnitude))
+    return tuple(rows)
+
+
+def pad_or_trim(samples: np.ndarray, length: int) -> np.ndarray:
+    x = np.asarray(samples, dtype=np.float32).reshape(-1)
+    if len(x) == length:
+        return x
+    if len(x) > length:
+        return x[:length]
+    out = np.zeros(length, dtype=np.float32)
+    out[: len(x)] = x
+    return out
+
+
+def rms(samples: np.ndarray) -> float:
+    x = np.asarray(samples, dtype=np.float32)
+    return float(np.sqrt(np.mean(x * x))) if x.size else 0.0
diff --git a/docs/audio-field.md b/docs/audio-field.md
index b99bb63..821c019 100644
--- a/docs/audio-field.md
+++ b/docs/audio-field.md
@@ -174,6 +174,14 @@ Suppress room/transient witness energy before FOA encoding:
 
 This first suppression pass treats the Focusrite dialogue anchors as wanted direct sources and the camera mics as room/transient witnesses. It attenuates witness-heavy broadband transients, lightly subtracts witness room energy from anchors, and writes a JSON report. It is deliberately a stream-side cleanup stage, not physical cancellation in the room.
 
+Suppress known speaker/program bleed using output loopback as ground truth:
+
+```powershell
+.\.venv\Scripts\python.exe .\scripts\audio_field.py suppress-reference --profile .\config\audio-field.json --input .\calibration\runs\<run-folder>\field-aligned.wav --reference .\calibration\runs\<run-folder>\ground_truth_loopback.wav --output .\calibration\runs\<run-folder>\field-program-suppressed.wav
+```
+
+This is the music-as-both-enemy-and-teacher path. The reference program audio is transformed into a time/frequency transfer estimate per mic channel, reconstructed as predicted bleed, subtracted from the aligned field, and exported as phase/frequency mapping evidence. A future live loop should publish the known program stem as its own OBS/Aquarium channel while using the same ground truth to keep learning the room.
+
 The `play-record` and `record-field` commands remain for a `shared-input-device` profile only. They deliberately refuse the distributed rig until an alignment stage emits one coherent six-channel WAV.
 
 ## Invariants
diff --git a/scripts/audio_field.py b/scripts/audio_field.py
index ceecfea..03a70b8 100644
--- a/scripts/audio_field.py
+++ b/scripts/audio_field.py
@@ -47,6 +47,19 @@ RoomSuppressionConfig = _room_suppression.RoomSuppressionConfig
 suppress_room_field = _room_suppression.suppress_room_field
 
 
+def load_program_reference_module():
+    spec = importlib.util.spec_from_file_location("localcast_program_reference", ROOT / "audio_field" / "program_reference.py")
+    module = importlib.util.module_from_spec(spec)
+    assert spec.loader is not None
+    sys.modules[spec.name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+_program_reference = load_program_reference_module()
+suppress_program_reference = _program_reference.suppress_program_reference
+
+
 def utc_stamp() -> str:
     return datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%S.%fZ")
 
@@ -1759,6 +1772,52 @@ def cmd_suppress_room(args):
     print(json.dumps(result, indent=2))
 
 
+def cmd_suppress_reference(args):
+    profile = load_profile(args.profile)
+    sample_rate, data = wavfile.read(args.input)
+    field = as_float_matrix(data)
+    ref_rate, reference = wavfile.read(args.reference)
+    reference = mono(as_float_matrix(reference))
+    reference = resample_if_needed(reference, int(ref_rate), int(sample_rate))
+    channels = [int(value) for value in args.channel] if args.channel else list(range(field.shape[1]))
+    cleaned, reports = suppress_program_reference(
+        field,
+        reference,
+        int(sample_rate),
+        channels,
+        nperseg=int(args.window_size),
+        noverlap=int(args.overlap),
+        regularization=float(args.regularization),
+        subtraction_strength=float(args.subtraction_strength),
+    )
+    output = args.output or args.input.with_name(args.input.stem + "-program-suppressed.wav")
+    wavfile.write(output, int(sample_rate), cleaned.astype(np.float32))
+    result = {
+        "output": str(output),
+        "reference": str(args.reference),
+        "sampleRate": int(sample_rate),
+        "shape": list(cleaned.shape),
+        "channels": [
+            {
+                "channel": report.channel,
+                "inputRms": report.input_rms,
+                "predictedRms": report.predicted_rms,
+                "outputRms": report.output_rms,
+                "reductionDb": report.reduction_db,
+                "phaseMapping": [
+                    {"frequencyHz": frequency, "phaseRadians": phase, "magnitude": magnitude}
+                    for frequency, phase, magnitude in report.phase_mapping
+                ],
+            }
+            for report in reports
+        ],
+        "notes": "Known program/reference bleed is transformed, mapped into each mic channel, subtracted, and exported as phase/frequency evidence for runtime fitting.",
+    }
+    report_path = args.report or output.with_suffix(".program-reference.json")
+    write_json(report_path, result)
+    print(json.dumps(result, indent=2))
+
+
 def infer_anchor_channels(profile):
     anchors = []
     for mic in profile["microphones"]:
@@ -2038,6 +2097,19 @@ def main():
     p.add_argument("--limit-peak", type=float, default=0.98)
     p.set_defaults(func=cmd_suppress_room)
 
+    p = sub.add_parser("suppress-reference", help="Subtract known speaker/program reference bleed from an aligned field and export phase mapping.")
+    p.add_argument("--profile", type=Path, default=ROOT / "config" / "audio-field.example.json")
+    p.add_argument("--input", type=Path, required=True)
+    p.add_argument("--reference", type=Path, required=True)
+    p.add_argument("--output", type=Path)
+    p.add_argument("--report", type=Path)
+    p.add_argument("--channel", type=int, action="append")
+    p.add_argument("--window-size", type=int, default=2048)
+    p.add_argument("--overlap", type=int, default=1536)
+    p.add_argument("--regularization", type=float, default=1e-6)
+    p.add_argument("--subtraction-strength", type=float, default=0.85)
+    p.set_defaults(func=cmd_suppress_reference)
+
     args = parser.parse_args()
     args.func(args)
 
diff --git a/tests/test_program_reference.py b/tests/test_program_reference.py
new file mode 100644
index 0000000..c845902
--- /dev/null
+++ b/tests/test_program_reference.py
@@ -0,0 +1,36 @@
+import unittest
+
+import numpy as np
+
+from audio_field.program_reference import suppress_program_reference
+
+
+class ProgramReferenceTests(unittest.TestCase):
+    def test_reference_bleed_is_suppressed_from_channel(self):
+        sample_rate = 48000
+        t = np.arange(sample_rate, dtype=np.float32) / sample_rate
+        reference = 0.2 * np.sin(2 * np.pi * 440 * t) + 0.1 * np.sin(2 * np.pi * 1400 * t)
+        voice = 0.03 * np.sin(2 * np.pi * 220 * t)
+        field = np.zeros((sample_rate, 2), dtype=np.float32)
+        field[:, 0] = voice + 0.7 * reference
+        field[:, 1] = voice
+
+        cleaned, reports = suppress_program_reference(
+            field,
+            reference,
+            sample_rate,
+            channels=[0],
+            nperseg=2048,
+            noverlap=1536,
+            subtraction_strength=1.0,
+        )
+
+        self.assertEqual((sample_rate, 2), cleaned.shape)
+        self.assertLess(reports[0].output_rms, reports[0].input_rms)
+        self.assertLess(float(np.sqrt(np.mean((cleaned[:, 0] - voice) ** 2))), 0.01)
+        np.testing.assert_allclose(cleaned[:, 1], field[:, 1])
+        self.assertGreater(len(reports[0].phase_mapping), 0)
+
+
+if __name__ == "__main__":
+    unittest.main()