diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 00000000..1e0bdf57
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,32 @@
+name: Tests
+
+on:
+  pull_request:
+  push:
+    branches: [main]
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install package and test dependencies
+        run: python -m pip install -e ".[test]"
+
+      - name: Run tests
+        run: python -m pytest tests/ -v
diff --git a/pyproject.toml b/pyproject.toml
index 0d0bda0f..0c4941a5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ authors = [
 ]
 description = "Python port of the MUST toolbox for ultrasound signal processing and generation of simulated images"
 readme = "README.md"
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 classifiers = [
     "Programming Language :: Python :: 3",
     "Operating System :: OS Independent",
@@ -24,6 +24,9 @@ dependencies = [
 ]
 dynamic = ["license", "version"]
 
+[project.optional-dependencies]
+test = ["pytest"]
+
 [project.urls]
 Homepage = "https://www.biomecardio.com/MUST"
 Repository = "https://github.com/creatis-ULTIM/PyMUST"
diff --git a/src/pymust/txdelay3.py b/src/pymust/txdelay3.py
index 712abeaf..50b72b00 100644
--- a/src/pymust/txdelay3.py
+++ b/src/pymust/txdelay3.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 from . import utils
 import numpy as np
 import scipy.optimize
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 00000000..8ade035b
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,12 @@
+import numpy as np
+import pytest
+
+
+@pytest.fixture(params=["L11-5v", "L12-3v", "C5-2v", "P4-2v"])
+def probe_name(request):
+    return request.param
+
+
+@pytest.fixture
+def rng():
+    return np.random.default_rng(0)
diff --git a/tests/data/pfield.npz b/tests/data/pfield.npz
new file mode 100644
index 00000000..40f15ded
Binary files /dev/null and b/tests/data/pfield.npz differ
diff --git a/tests/data/simus.npz b/tests/data/simus.npz
new file mode 100644
index 00000000..d4c22e2f
Binary files /dev/null and b/tests/data/simus.npz differ
diff --git a/tests/generate_golden_data.py b/tests/generate_golden_data.py
new file mode 100644
index 00000000..c2087bed
--- /dev/null
+++ b/tests/generate_golden_data.py
@@ -0,0 +1,23 @@
+"""Regenerate the golden reference .npz files used by test_regression.py.
+
+Run this deliberately after a change that intentionally alters simus/pfield
+numerics, then review the diff of the resulting .npz files:
+
+    python tests/generate_golden_data.py
+"""
+import numpy as np
+
+from golden_scenarios import DATA_DIR, SCENARIOS
+
+
+def main():
+    DATA_DIR.mkdir(exist_ok=True)
+    for name, build in SCENARIOS.items():
+        outputs = build()
+        path = DATA_DIR / f"{name}.npz"
+        np.savez(path, **outputs)
+        print(f"wrote {path}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/golden_scenarios.py b/tests/golden_scenarios.py
new file mode 100644
index 00000000..0aa150f4
--- /dev/null
+++ b/tests/golden_scenarios.py
@@ -0,0 +1,52 @@
+"""Fixed input scenarios shared by the golden-data generator and the regression
+tests, so both sides always compute from the exact same inputs.
+"""
+import numpy as np
+import pymust
+
+DATA_DIR = __import__("pathlib").Path(__file__).parent / "data"
+
+
+PLANE_WAVE_TILTS = [-0.3, -0.1, 0.0, 0.1, 0.3]  # radians, several steering directions
+
+
+def pfield_scenario():
+    param = pymust.getparam("L11-5v")
+    param.Nelements = 16
+
+    x, z = pymust.impolgrid(10, 0.04, np.pi / 4, param)
+    y = np.zeros_like(x)
+
+    outputs = {}
+    for tilt in PLANE_WAVE_TILTS:
+        delays = pymust.txdelayPlane(param, tilt)
+        rp, _, _ = pymust.pfield(x, y, z, delays, param)
+        outputs[f"rms_pressure_plane_tilt_{tilt:+.1f}"] = rp
+
+    focused_delays = pymust.txdelayFocused(param, 0.005, 0.03)
+    rp_focused, _, _ = pymust.pfield(x, y, z, focused_delays, param)
+    outputs["rms_pressure_focused"] = rp_focused
+
+    return outputs
+
+
+def simus_scenario():
+    param = pymust.getparam("L11-5v")
+    param.Nelements = 8
+
+    x = np.array([0.0, 0.3e-2, -0.3e-2])
+    y = np.array([0.0, 0.0, 0.0])
+    z = np.array([0.2e-2, 0.4e-2, 0.6e-2])
+    rc = np.array([1.0, 0.7, 1.3])
+    delays = pymust.txdelayPlane(param, 0.05)
+
+    rf, _ = pymust.simus(x, y, z, rc, delays, param)
+    iq = pymust.rf2iq(rf, param)
+    img = pymust.bmode(iq)
+    return {"rf": rf, "iq": iq, "bmode": img}
+
+
+SCENARIOS = {
+    "pfield": pfield_scenario,
+    "simus": simus_scenario,
+}
diff --git a/tests/test_getparam.py b/tests/test_getparam.py
new file mode 100644
index 00000000..c01d68ae
--- /dev/null
+++ b/tests/test_getparam.py
@@ -0,0 +1,18 @@
+import numpy as np
+import pymust
+
+
+def test_getparam_returns_sane_probe_parameters(probe_name):
+    param = pymust.getparam(probe_name)
+
+    assert param.Nelements > 0
+    assert param.fc > 0
+    assert param.pitch > 0
+    assert param.bandwidth > 0
+
+
+def test_getparam_unknown_probe_raises():
+    import pytest
+
+    with pytest.raises(Exception):
+        pymust.getparam("not-a-real-probe")
diff --git a/tests/test_imports.py b/tests/test_imports.py
new file mode 100644
index 00000000..c6dc81b8
--- /dev/null
+++ b/tests/test_imports.py
@@ -0,0 +1,38 @@
+import pymust
+
+
+EXPECTED_FUNCTIONS = [
+    "bmode",
+    "dasmtx",
+    "dasmtx3",
+    "getparam",
+    "impolgrid",
+    "iq2doppler",
+    "getNyquistVelocity",
+    "pfield",
+    "pfield3",
+    "rf2iq",
+    "simus",
+    "simus3",
+    "tgc",
+    "txdelay",
+    "txdelayCircular",
+    "txdelayPlane",
+    "txdelayFocused",
+    "txdelay3",
+    "txdelay3Plane",
+    "txdelay3Diverging",
+    "txdelay3Focused",
+    "getDopplerColorMap",
+    "genscat",
+    "mkmovie",
+    "getpulse",
+    "smoothn",
+    "sptrack",
+]
+
+
+def test_all_public_functions_are_importable():
+    for name in EXPECTED_FUNCTIONS:
+        assert hasattr(pymust, name), f"pymust.{name} is missing"
+        assert callable(getattr(pymust, name))
diff --git a/tests/test_regression.py b/tests/test_regression.py
new file mode 100644
index 00000000..6ba29f53
--- /dev/null
+++ b/tests/test_regression.py
@@ -0,0 +1,87 @@
+"""Golden-master regression tests.
+
+These compare simus/pfield output, computed fresh on every test run, against
+reference arrays saved in tests/data/. simus and pfield are deterministic
+within one environment (no randomness - see tests/golden_scenarios.py), but
+NOT bit-reproducible across environments: pfield accumulates results in
+single precision (complex64), and different BLAS/LAPACK backends (e.g.
+OpenBLAS vs Apple's Accelerate) round matrix reductions differently. simus
+then runs a sharp threshold ("zero out samples below -100dB relative to the
+peak", see the RelThresh/tanh step in simus.py) on top of that, which turns
+tiny backend-dependent noise into visibly different values for the handful
+of samples that sit right at the threshold. This was caught in practice by
+comparing a conda-env-generated reference against a plain pip install: exact
+comparison failed on ~36% of RF samples even though every value was
+numerically negligible (all well under 1e-4 of the peak amplitude).
+
+So a sample's relative tolerance depends on how loud it is compared to the
+signal's own peak (in dB, 20*log10(|value|/peak)):
+  - at or above NOISE_FLOOR_DB: tight rtol - this is "real" signal, and a
+    regression here (wrong scaling, wrong timing/shape, wrong physics)
+    should still fail loudly.
+  - below NOISE_FLOOR_DB: loose rtol - down in the noise floor, relative
+    comparison is meaningless (dividing by a near-zero reference blows up
+    the ratio) and this is exactly where backend rounding causes samples to
+    snap across simus's threshold. A small absolute tolerance (scaled to
+    the peak) still catches a sample that shouldn't be near-zero at all.
+
+If you intentionally change the numerics of simus/pfield (e.g. a bug fix that
+changes the output), regenerate the references and review the diff:
+
+    python tests/generate_golden_data.py
+"""
+import numpy as np
+import pytest
+
+from golden_scenarios import DATA_DIR, SCENARIOS
+
+NOISE_FLOOR_DB = -30.0
+TIGHT_RTOL = 1e-4
+LOOSE_RTOL = 1.0
+ATOL_FRACTION_OF_PEAK = 1e-5
+
+
+def _load_reference(name):
+    path = DATA_DIR / f"{name}.npz"
+    if not path.exists():
+        pytest.skip(f"no golden reference at {path}; run generate_golden_data.py")
+    with np.load(path) as data:
+        return {key: data[key] for key in data.files}
+
+
+def _assert_matches_reference(fresh_value, ref_value, label):
+    fresh_value = np.asarray(fresh_value, dtype=np.complex128 if np.iscomplexobj(fresh_value) else np.float64)
+    ref_value = np.asarray(ref_value, dtype=np.complex128 if np.iscomplexobj(ref_value) else np.float64)
+    assert fresh_value.shape == ref_value.shape, label
+
+    peak = np.max(np.abs(ref_value))
+    atol = max(peak * ATOL_FRACTION_OF_PEAK, 1e-9)
+
+    with np.errstate(divide="ignore"):
+        ref_db = 20 * np.log10(np.abs(ref_value) / peak)
+    rtol = np.where(ref_db < NOISE_FLOOR_DB, LOOSE_RTOL, TIGHT_RTOL)
+
+    diff = np.abs(fresh_value - ref_value)
+    allowed = atol + rtol * np.abs(ref_value)
+    bad = diff > allowed
+
+    if bad.any():
+        worst = np.unravel_index(np.argmax(diff - allowed), diff.shape)
+        raise AssertionError(
+            f"{label}: {bad.sum()}/{bad.size} elements exceed tolerance "
+            f"(no longer matches the golden reference); worst at {worst}: "
+            f"fresh={fresh_value[worst]!r} ref={ref_value[worst]!r} "
+            f"diff={diff[worst]:.3g} allowed={allowed[worst]:.3g} "
+            f"({ref_db[worst]:.1f} dB relative to peak)"
+        )
+
+
+@pytest.mark.parametrize("scenario_name", sorted(SCENARIOS))
+def test_matches_golden_reference(scenario_name):
+    reference = _load_reference(scenario_name)
+    fresh = SCENARIOS[scenario_name]()
+
+    assert set(fresh) == set(reference), "scenario outputs changed - regenerate golden data"
+
+    for key, fresh_value in fresh.items():
+        _assert_matches_reference(fresh_value, reference[key], f"{scenario_name}.{key}")
diff --git a/tests/test_simus_pipeline.py b/tests/test_simus_pipeline.py
new file mode 100644
index 00000000..e45adcf5
--- /dev/null
+++ b/tests/test_simus_pipeline.py
@@ -0,0 +1,28 @@
+import numpy as np
+import pymust
+
+
+def test_simus_rf2iq_bmode_pipeline():
+    param = pymust.getparam("L11-5v")
+    param.Nelements = 8
+
+    x = np.array([0.0, 0.0])
+    y = np.array([0.0, 0.0])
+    z = np.array([0.2e-2, 0.5e-2])
+    rc = np.array([1.0, 1.0])
+    tx_delays = np.zeros(param.Nelements).reshape((1, -1))
+
+    rf, _ = pymust.simus(x, y, z, rc, tx_delays, param)
+    assert rf.ndim == 2
+    assert rf.shape[1] == param.Nelements
+    assert np.isfinite(rf).all()
+    assert np.any(rf != 0)
+
+    iq = pymust.rf2iq(rf, param)
+    assert iq.shape == rf.shape
+    assert np.iscomplexobj(iq)
+    assert np.isfinite(iq).all()
+
+    img = pymust.bmode(iq)
+    assert img.shape == rf.shape
+    assert img.dtype == np.uint8
diff --git a/tests/test_txdelay.py b/tests/test_txdelay.py
new file mode 100644
index 00000000..46aaab12
--- /dev/null
+++ b/tests/test_txdelay.py
@@ -0,0 +1,23 @@
+import numpy as np
+import pymust
+
+
+def test_txdelay_plane_and_focused(probe_name):
+    param = pymust.getparam(probe_name)
+
+    plane = pymust.txdelayPlane(param, 0.1)
+    focused = pymust.txdelayFocused(param, 0, 0.03)
+
+    for delays in (plane, focused):
+        assert delays.shape[-1] == param.Nelements
+        assert np.isfinite(delays).all()
+
+
+def test_txdelay_circular_wave_on_linear_array():
+    # txdelayCircular is only defined for linear (non-curved) arrays.
+    param = pymust.getparam("L11-5v")
+
+    circular = pymust.txdelayCircular(param, 0.1, np.pi / 3)
+
+    assert circular.shape[-1] == param.Nelements
+    assert np.isfinite(circular).all()
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 00000000..acd7a76c
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,40 @@
+import numpy as np
+import pymust
+
+
+def test_genscat_generates_scatterers_within_roi():
+    roi_width, roi_depth = 0.02, 0.02
+    x, y, z, rc = pymust.genscat(np.array([roi_width, roi_depth]), 0.001)
+
+    assert x.shape == y.shape == z.shape == rc.shape
+    assert x.size > 0
+    # The ROI is centered on x=0 with its top edge at z=0, per genscat's docstring.
+    assert np.abs(x).max() <= roi_width / 2 + 1e-9
+    assert (z >= -1e-9).all()
+    assert z.max() <= roi_depth + 1e-9
+
+
+def test_smoothn_reduces_noise(rng):
+    t = np.linspace(0, 10, 200)
+    clean = np.sin(t)
+    noisy = clean + 0.3 * rng.standard_normal(t.size)
+
+    smoothed, _, _ = pymust.smoothn(noisy)
+
+    assert smoothed.shape == noisy.shape
+    assert np.isfinite(smoothed).all()
+    error_before = np.mean((noisy - clean) ** 2)
+    error_after = np.mean((smoothed - clean) ** 2)
+    assert error_after < error_before
+
+
+def test_impolgrid_returns_grid_coordinates():
+    param = pymust.getparam("L11-5v")
+
+    grid = pymust.impolgrid(50, 0.05, np.pi / 3, param)
+
+    assert len(grid) == 2
+    x, z = grid
+    assert x.shape == z.shape
+    assert np.isfinite(x).all()
+    assert np.isfinite(z).all()