diff --git a/notebooks/HippocampalPlaceCellExample.ipynb b/notebooks/HippocampalPlaceCellExample.ipynb
index 04e3ad7d..9789a3ac 100644
--- a/notebooks/HippocampalPlaceCellExample.ipynb
+++ b/notebooks/HippocampalPlaceCellExample.ipynb
@@ -2,20 +2,20 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "id": "4110300d",
+   "id": "39256c2a",
    "metadata": {},
    "source": [
     "<!-- parity-note -->\n",
     "## MATLAB Parity Note\n",
     "- Source MATLAB helpfile: `HippocampalPlaceCellExample.mlx`\n",
     "- Fidelity status: `high_fidelity`\n",
-    "- Remaining justified differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with real figures; the Python port still uses an approximate Zernike-like basis rather than the original MATLAB toolbox implementation.\n"
+    "- Remaining justified differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with the same normalized 10-term Zernike basis used by MATLAB; exact AIC/BIC values and surface styling still vary modestly because the Python GLM solver and plotting backend are not byte-identical to MATLAB.\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8c6412bd",
+   "id": "2b552e70",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -85,7 +85,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "07fa2765",
+   "id": "02b3faec",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -105,7 +105,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "1a5876df",
+   "id": "3e5df761",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -125,7 +125,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "db6754b1",
+   "id": "22988fc3",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -152,7 +152,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "5c98b75c",
+   "id": "2bb2d885",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -179,7 +179,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "e4dee3d3",
+   "id": "34d08f71",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -246,7 +246,7 @@
     "            f\"Cells analyzed: {int(summary['num_cells_fit'])}\",\n",
     "            f\"Mean Gaussian-Zernike ΔAIC: {summary['mean_delta_aic_gaussian_minus_zernike']:.2f}\",\n",
     "            f\"Mean Gaussian-Zernike ΔBIC: {summary['mean_delta_bic_gaussian_minus_zernike']:.2f}\",\n",
-    "            \"Negative values favor the Zernike-like model.\",\n",
+    "            \"Negative values favor the Zernike model.\",\n",
     "        ]\n",
     "    ),\n",
     "    va=\"top\",\n",
diff --git a/notebooks/nSTATPaperExamples.ipynb b/notebooks/nSTATPaperExamples.ipynb
index 39f762b9..093de5b2 100644
--- a/notebooks/nSTATPaperExamples.ipynb
+++ b/notebooks/nSTATPaperExamples.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "id": "f7aeaead",
+   "id": "989e2794",
    "metadata": {},
    "source": [
     "<!-- parity-note -->\n",
@@ -15,7 +15,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "5224c304",
+   "id": "e3ae0b3f",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -75,7 +75,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "7479e214",
+   "id": "ad04d480",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -86,7 +86,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8aad9697",
+   "id": "f5fd7ff0",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -102,7 +102,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "7190f56b",
+   "id": "2b8b5cc6",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -113,7 +113,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "1c1ee772",
+   "id": "36d369c6",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -138,7 +138,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "be2211f0",
+   "id": "780da944",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -158,7 +158,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "33028dc0",
+   "id": "3c220adb",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -176,7 +176,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8afc80fe",
+   "id": "78a302f9",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -194,7 +194,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "52f1ba93",
+   "id": "2b09c5e3",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -205,7 +205,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "6b5263c4",
+   "id": "03f25c28",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -224,7 +224,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "45a8e34e",
+   "id": "03ee9c55",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -240,7 +240,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "6e33b0e3",
+   "id": "09f0d4f1",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -259,7 +259,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "6376eee2",
+   "id": "b32115d9",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -281,7 +281,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "d3e45147",
+   "id": "d209b1e8",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -301,7 +301,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "d859424b",
+   "id": "570fa694",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -321,7 +321,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "cbf7d0c0",
+   "id": "2e5c472a",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -332,7 +332,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "1765fe2c",
+   "id": "ad2e93f8",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -348,7 +348,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "2b42cc95",
+   "id": "03fe133c",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -366,7 +366,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "5806ab16",
+   "id": "2be08d3f",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -377,7 +377,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "5354140f",
+   "id": "4470cf1f",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -394,7 +394,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "cec81566",
+   "id": "2c42e63b",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -410,7 +410,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "f9f22950",
+   "id": "e193a719",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -426,7 +426,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "92b61189",
+   "id": "3718c5d4",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -437,7 +437,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0b46271e",
+   "id": "96e2af7c",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -453,7 +453,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "ab1e1e3b",
+   "id": "ccd871af",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -469,7 +469,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "71966f40",
+   "id": "113c1903",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -483,21 +483,21 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "449c2417",
+   "id": "64654717",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# SECTION 26: Zernike-like place-field mesh\n",
+    "# SECTION 26: Zernike place-field mesh\n",
     "fig = _fig(\"experiment4 zernike mesh\", figsize=(9.0, 6.5))\n",
     "ax = fig.add_subplot(111, projection=\"3d\")\n",
     "ax.plot_surface(exp4[\"mesh\"][\"grid_x\"], exp4[\"mesh\"][\"grid_y\"], exp4[\"mesh\"][\"zernike_field\"], cmap=\"Greens\", linewidth=0.0, antialiased=True)\n",
-    "ax.set_title(\"Zernike-like place-field estimate\")\n"
+    "ax.set_title(\"Zernike place-field estimate\")\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "198fdb9f",
+   "id": "115929b8",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -508,7 +508,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8a6f7766",
+   "id": "26af3ac2",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -526,7 +526,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "5d51d0c1",
+   "id": "22ee7cab",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -537,7 +537,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "cffe4e46",
+   "id": "f0312bf9",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -556,7 +556,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "17a6d341",
+   "id": "216a748d",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -575,7 +575,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a38816f8",
+   "id": "358d0d12",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -586,7 +586,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "24daef12",
+   "id": "be288a10",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -604,7 +604,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "42056ccb",
+   "id": "7e158b30",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -623,7 +623,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a0dcaad0",
+   "id": "2b0b149c",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -641,7 +641,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "61804aaf",
+   "id": "74cb3002",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -664,7 +664,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "3f0a28f1",
+   "id": "a52f8e5c",
    "metadata": {},
    "outputs": [],
    "source": [
diff --git a/nstat/core.py b/nstat/core.py
index 2190da94..31b03234 100644
--- a/nstat/core.py
+++ b/nstat/core.py
@@ -723,12 +723,33 @@ def resample(self, sample_rate: float) -> "SignalObj":
         return copied
 
     def resampleMe(self, newSampleRate: float) -> None:
+        try:
+            from scipy.interpolate import interp1d
+        except Exception as exc:  # pragma: no cover
+            raise ImportError("scipy is required for SignalObj.resampleMe") from exc
+
         rate = float(newSampleRate)
         if rate <= 0:
             raise ValueError("sampleRate must be > 0.")
+        if self.sampleRate == rate:
+            return
+        self.restoreToOriginal()
         dt = 1.0 / rate
         newTime = np.arange(self.time[0], self.time[-1] + 0.5 * dt, dt, dtype=float)
-        newData = np.column_stack([np.interp(newTime, self.time, self.data[:, i]) for i in range(self.dimension)])
+        if self.data.shape[0] > 1:
+            columns = []
+            for index in range(self.dimension):
+                interpolator = interp1d(
+                    self.time,
+                    self.data[:, index],
+                    kind="cubic",
+                    bounds_error=False,
+                    fill_value=0.0,
+                )
+                columns.append(np.asarray(interpolator(newTime), dtype=float))
+            newData = np.column_stack(columns)
+        else:
+            newData = np.asarray(self.data, dtype=float).copy()
         self.time = newTime
         self.data = newData
         self.sampleRate = rate
@@ -737,7 +758,9 @@ def resampleMe(self, newSampleRate: float) -> None:
 
     @property
     def derivative(self) -> "SignalObj":
-        deriv = np.column_stack([np.gradient(self.data[:, i], self.time) for i in range(self.dimension)])
+        deriv = np.zeros_like(self.data, dtype=float)
+        if self.data.shape[0] > 1:
+            deriv[1:, :] = np.diff(self.data, axis=0) * float(self.sampleRate)
         deriv[~np.isfinite(deriv)] = 0.0
         labels = [f"d_{label}" if label else "" for label in self.dataLabels]
         return self._spawn(self.time, deriv, data_labels=labels)
@@ -1260,7 +1283,7 @@ def computeStatistics(self, makePlots: int = 0) -> None:
         self.Lstatistic = np.nan
 
         if isi.size:
-            sigma = float(np.std(isi))
+            sigma = float(np.std(isi, ddof=1)) if isi.size > 1 else 0.0
             mu = float(np.mean(isi))
             if np.isfinite(mu) and mu > 0:
                 r = sigma / mu
diff --git a/nstat/paper_examples.py b/nstat/paper_examples.py
index 3a78ca26..f2bda945 100644
--- a/nstat/paper_examples.py
+++ b/nstat/paper_examples.py
@@ -13,6 +13,7 @@
 from .decoding_algorithms import DecodingAlgorithms
 from .glm import fit_poisson_glm
 from .simulation import simulate_poisson_from_rate
+from .zernike import zernike_basis_from_cartesian
 
 
 Summary = dict[str, float]
@@ -165,25 +166,6 @@ def _spike_indicator_from_times(time: np.ndarray, spike_times: np.ndarray) -> np
     return y
 
 
-def _zernike_like_basis(x: np.ndarray, y: np.ndarray) -> np.ndarray:
-    theta = np.arctan2(y, x)
-    r = np.sqrt(x * x + y * y)
-    return np.column_stack(
-        [
-            np.ones_like(r),
-            r,
-            r**2,
-            np.cos(theta),
-            np.sin(theta),
-            r * np.cos(theta),
-            r * np.sin(theta),
-            r**2 * np.cos(2.0 * theta),
-            r**2 * np.sin(2.0 * theta),
-            r**3,
-        ]
-    )
-
-
 def run_experiment4(data_dir: Path, *, return_payload: bool = False) -> Summary | tuple[Summary, Payload]:
     mat_path = data_dir / "Place Cells" / "PlaceCellDataAnimal1.mat"
     d = loadmat(mat_path, squeeze_me=True, struct_as_record=False)
@@ -197,7 +179,7 @@ def run_experiment4(data_dir: Path, *, return_payload: bool = False) -> Summary
     offset = np.full(time.shape[0], np.log(max(dt, 1e-12)), dtype=float)
 
     x_gauss = np.column_stack([x, y, x * x, y * y, x * y])
-    x_zern = _zernike_like_basis(x, y)
+    x_zern = zernike_basis_from_cartesian(x, y)
 
     n_eval = int(min(8, neurons.shape[0]))
     delta_aic = []
diff --git a/nstat/paper_examples_full.py b/nstat/paper_examples_full.py
index 0cb94845..83652164 100644
--- a/nstat/paper_examples_full.py
+++ b/nstat/paper_examples_full.py
@@ -14,6 +14,7 @@
 from .decoding_algorithms import DecodingAlgorithms
 from .glm import fit_poisson_glm
 from .simulation import simulate_poisson_from_rate
+from .zernike import zernike_basis_from_cartesian
 
 
 def _default_repo_root() -> Path:
@@ -460,23 +461,6 @@ def _spike_indicator(time: np.ndarray, spike_times: np.ndarray) -> np.ndarray:
     return y
 
 
-def _zernike_like_basis(x: np.ndarray, y: np.ndarray) -> np.ndarray:
-    theta = np.arctan2(y, x)
-    r = np.sqrt(x * x + y * y)
-    return np.column_stack([
-        np.ones_like(r),
-        r,
-        r**2,
-        np.cos(theta),
-        np.sin(theta),
-        r * np.cos(theta),
-        r * np.sin(theta),
-        r**2 * np.cos(2.0 * theta),
-        r**2 * np.sin(2.0 * theta),
-        r**3,
-    ])
-
-
 def _load_placecell_dataset(path: Path):
     d = _loadmat_checked(path)
     if d is None:
@@ -484,6 +468,12 @@ def _load_placecell_dataset(path: Path):
         time = np.linspace(0.0, 20.0, 2400, dtype=float)
         x = 0.8 * np.sin(0.6 * time) + 0.2 * np.sin(1.7 * time + 0.5)
         y = 0.7 * np.cos(0.5 * time + 0.3)
+        radius = np.sqrt(x * x + y * y)
+        max_radius = float(np.max(radius)) if radius.size else 0.0
+        if max_radius > 0.98:
+            scale = 0.98 / max_radius
+            x = x * scale
+            y = y * scale
         n_cells = 8
         neurons = []
         for _ in range(n_cells):
@@ -504,13 +494,13 @@ def _evaluate_place_models(x: np.ndarray, y: np.ndarray, time: np.ndarray, neuro
     dt = float(np.median(np.diff(time)))
     offset = np.full(time.shape[0], np.log(max(dt, 1e-12)), dtype=float)
     x_gauss = np.column_stack([x, y, x * x, y * y, x * y])
-    x_zern = _zernike_like_basis(x, y)
+    x_zern = zernike_basis_from_cartesian(x, y)
 
     x_grid = np.linspace(float(np.min(x)), float(np.max(x)), 40, dtype=float)
     y_grid = np.linspace(float(np.min(y)), float(np.max(y)), 40, dtype=float)
     xx, yy = np.meshgrid(x_grid, y_grid)
     grid_gauss = np.column_stack([xx.ravel(), yy.ravel(), xx.ravel() ** 2, yy.ravel() ** 2, xx.ravel() * yy.ravel()])
-    grid_zern = _zernike_like_basis(xx.ravel(), yy.ravel())
+    grid_zern = zernike_basis_from_cartesian(xx.ravel(), yy.ravel(), fill_value=np.nan)
 
     d_aic: list[float] = []
     d_bic: list[float] = []
diff --git a/nstat/release_check.py b/nstat/release_check.py
new file mode 100644
index 00000000..e8848f8d
--- /dev/null
+++ b/nstat/release_check.py
@@ -0,0 +1,96 @@
+from __future__ import annotations
+
+import argparse
+import shutil
+import subprocess
+import sys
+from pathlib import Path
+
+
+def _repo_root() -> Path:
+    return Path(__file__).resolve().parents[1]
+
+
+def default_matlab_repo_root(repo_root: Path | None = None) -> Path:
+    base = _repo_root() if repo_root is None else repo_root.resolve()
+    return base.parent / "nSTAT"
+
+
+def _matlab_quote(value: str) -> str:
+    return value.replace("'", "''")
+
+
+def build_release_gate_commands(
+    repo_root: Path | None = None,
+    *,
+    matlab_repo_root: Path | None = None,
+    skip_matlab: bool = False,
+) -> list[list[str]]:
+    base = _repo_root() if repo_root is None else repo_root.resolve()
+    matlab_root = default_matlab_repo_root(base) if matlab_repo_root is None else matlab_repo_root.resolve()
+
+    commands: list[list[str]] = [
+        [sys.executable, str(base / "tools" / "parity" / "export_matlab_gold_fixtures.py"), "--repo-root", str(base), "--matlab-repo", str(matlab_root)],
+        [
+            sys.executable,
+            "-m",
+            "pytest",
+            "-q",
+            "tests/test_zernike.py",
+            "tests/test_matlab_gold_fixtures.py",
+            "tests/test_signalobj_fidelity.py",
+            "tests/test_nspiketrain_fidelity.py",
+            "tests/test_workflow_fidelity.py",
+            "tests/test_matlab_reference.py",
+            "tests/test_simulink_fidelity_audit.py",
+            "tests/test_parity_report.py",
+        ],
+        [sys.executable, str(base / "tools" / "notebooks" / "run_notebooks.py"), "--group", "parity_core", "--timeout", "900"],
+        [sys.executable, str(base / "tools" / "parity" / "build_report.py")],
+    ]
+    if not skip_matlab:
+        commands.append(
+            [
+                "matlab",
+                "-batch",
+                (
+                    f"cd('{_matlab_quote(str(matlab_root))}'); "
+                    "addpath(pwd); "
+                    "addpath(fullfile(pwd,'helpfiles')); "
+                    "results = runtests('tests/python_port_fidelity'); "
+                    "assertSuccess(results); exit"
+                ),
+            ]
+        )
+    return commands
+
+
+def run_release_gate(
+    repo_root: Path | None = None,
+    *,
+    matlab_repo_root: Path | None = None,
+    skip_matlab: bool = False,
+) -> None:
+    base = _repo_root() if repo_root is None else repo_root.resolve()
+    if not skip_matlab and shutil.which("matlab") is None:
+        raise FileNotFoundError("`matlab` is not on PATH; rerun with --skip-matlab or install MATLAB CLI access")
+    for command in build_release_gate_commands(base, matlab_repo_root=matlab_repo_root, skip_matlab=skip_matlab):
+        subprocess.run(command, cwd=base, check=True)
+
+
+def main(argv: list[str] | None = None) -> int:
+    parser = argparse.ArgumentParser(description="Run the Python plus MATLAB fidelity release gate.")
+    parser.add_argument("--repo-root", type=Path, default=None)
+    parser.add_argument("--matlab-repo", type=Path, default=None)
+    parser.add_argument("--skip-matlab", action="store_true")
+    args = parser.parse_args(argv)
+
+    run_release_gate(args.repo_root, matlab_repo_root=args.matlab_repo, skip_matlab=args.skip_matlab)
+    return 0
+
+
+__all__ = ["build_release_gate_commands", "default_matlab_repo_root", "main", "run_release_gate"]
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv[1:]))
diff --git a/nstat/zernike.py b/nstat/zernike.py
new file mode 100644
index 00000000..24513ee7
--- /dev/null
+++ b/nstat/zernike.py
@@ -0,0 +1,118 @@
+from __future__ import annotations
+
+from math import comb, sqrt
+from typing import Iterable
+
+import numpy as np
+
+
+DEFAULT_ZERNIKE_MODES: tuple[tuple[int, int], ...] = (
+    (0, 0),
+    (1, -1),
+    (1, 1),
+    (2, -2),
+    (2, 0),
+    (2, 2),
+    (3, -3),
+    (3, -1),
+    (3, 1),
+    (3, 3),
+)
+
+
+def _as_mode_vector(values: Iterable[int] | np.ndarray, name: str) -> np.ndarray:
+    array = np.asarray(list(values) if not isinstance(values, np.ndarray) else values, dtype=int).reshape(-1)
+    if array.ndim != 1:
+        raise ValueError(f"{name} must be a vector")
+    return array
+
+
+def _as_float_vector(values: Iterable[float] | np.ndarray, name: str) -> np.ndarray:
+    array = np.asarray(values, dtype=float).reshape(-1)
+    if array.ndim != 1:
+        raise ValueError(f"{name} must be a vector")
+    return array
+
+
+def _radial_polynomial(n: int, m_abs: int, r: np.ndarray) -> np.ndarray:
+    radial = np.zeros(r.shape[0], dtype=float)
+    for s in range((n - m_abs) // 2 + 1):
+        coefficient = ((-1) ** s) * comb(n - s, s) * comb(n - 2 * s, (n - m_abs) // 2 - s)
+        radial += float(coefficient) * np.power(r, n - 2 * s, dtype=float)
+    return radial
+
+
+def zernfun(
+    n: Iterable[int] | np.ndarray,
+    m: Iterable[int] | np.ndarray,
+    r: Iterable[float] | np.ndarray,
+    theta: Iterable[float] | np.ndarray,
+    *,
+    normalized: bool = False,
+) -> np.ndarray:
+    """Evaluate MATLAB-style Zernike functions on the unit disk."""
+
+    n_vec = _as_mode_vector(n, "n")
+    m_vec = _as_mode_vector(m, "m")
+    r_vec = _as_float_vector(r, "r")
+    theta_vec = _as_float_vector(theta, "theta")
+
+    if n_vec.shape != m_vec.shape:
+        raise ValueError("n and m must be the same length")
+    if r_vec.shape != theta_vec.shape:
+        raise ValueError("r and theta must be the same length")
+    if np.any((r_vec < 0.0) | (r_vec > 1.0)):
+        raise ValueError("All radial coordinates must be between 0 and 1")
+
+    out = np.zeros((r_vec.shape[0], n_vec.shape[0]), dtype=float)
+    for column, (n_val, m_val) in enumerate(zip(n_vec.tolist(), m_vec.tolist(), strict=False)):
+        m_abs = abs(int(m_val))
+        if n_val < 0:
+            raise ValueError("All n values must be non-negative")
+        if m_abs > n_val or (n_val - m_abs) % 2 != 0:
+            raise ValueError("Each |m| must be <= n and differ from n by a multiple of 2")
+
+        radial = _radial_polynomial(int(n_val), m_abs, r_vec)
+        if normalized:
+            radial *= sqrt(float(n_val + 1)) if m_val == 0 else sqrt(float(2 * (n_val + 1)))
+
+        if m_val > 0:
+            out[:, column] = radial * np.sin(theta_vec * float(m_val))
+        elif m_val < 0:
+            out[:, column] = radial * np.cos(theta_vec * float(m_val))
+        else:
+            out[:, column] = radial
+    return out
+
+
+def zernike_basis_from_cartesian(
+    x: Iterable[float] | np.ndarray,
+    y: Iterable[float] | np.ndarray,
+    *,
+    modes: tuple[tuple[int, int], ...] = DEFAULT_ZERNIKE_MODES,
+    normalized: bool = True,
+    fill_value: float | None = None,
+) -> np.ndarray:
+    x_vec = _as_float_vector(x, "x")
+    y_vec = _as_float_vector(y, "y")
+    if x_vec.shape != y_vec.shape:
+        raise ValueError("x and y must be the same length")
+
+    theta = np.arctan2(y_vec, x_vec)
+    r = np.sqrt(x_vec * x_vec + y_vec * y_vec)
+    inside = r <= (1.0 + 1e-12)
+    n = np.asarray([mode[0] for mode in modes], dtype=int)
+    m = np.asarray([mode[1] for mode in modes], dtype=int)
+
+    if fill_value is None:
+        if not np.all(inside):
+            raise ValueError("Cartesian coordinates must lie within the unit disk")
+        return zernfun(n, m, np.clip(r, 0.0, 1.0), theta, normalized=normalized)
+
+    out = np.full((x_vec.shape[0], len(modes)), float(fill_value), dtype=float)
+    if np.any(inside):
+        out[inside] = zernfun(n, m, np.clip(r[inside], 0.0, 1.0), theta[inside], normalized=normalized)
+    return out
+
+
+__all__ = ["DEFAULT_ZERNIKE_MODES", "zernfun", "zernike_basis_from_cartesian"]
diff --git a/parity/README.md b/parity/README.md
index 8ceccd58..e0123ae4 100644
--- a/parity/README.md
+++ b/parity/README.md
@@ -14,6 +14,18 @@ Generated report:
 python tools/parity/build_report.py
 ```
 
+Refresh MATLAB-derived fixtures:
+
+```bash
+python tools/parity/export_matlab_gold_fixtures.py
+```
+
+Run the combined Python plus MATLAB release gate:
+
+```bash
+nstat-release-check --matlab-repo ../nSTAT
+```
+
 Current headline status:
 - Public API coverage matches the MATLAB inventory except for the explicitly non-applicable `nstatOpenHelpPage`.
 - Class-fidelity auditing is tracked separately from name-mapping parity in `class_fidelity.yml`, and it remains intentionally stricter and more conservative than the mapping manifest.
diff --git a/parity/class_fidelity.yml b/parity/class_fidelity.yml
index 10db0f1e..14f6e1a2 100644
--- a/parity/class_fidelity.yml
+++ b/parity/class_fidelity.yml
@@ -40,8 +40,9 @@ items:
     unported.
   - Structure serialization is close but not exhaustive for every MATLAB-only field.
   required_remediation:
-  - Add MATLAB-derived fixtures for filter outputs, xcorr/autocorrelation traces, plotting
-    selectors, and the remaining spectral utility methods.
+  - Extend the committed MATLAB-derived fixtures beyond derivative, integral, spline
+    resampling, filtering, and xcorr to cover autocorrelation selectors and the remaining
+    spectral utility methods.
   plotting_report_parity: Core plotting and correlation helpers are implemented; some
     MATLAB-only plot selectors, spectral utilities, and report-style helpers remain
     lighter.
@@ -102,8 +103,9 @@ items:
   - Some MATLAB visual styling and distribution-fit detail in the ISI plotting helpers
     remains lighter than MATLAB.
   required_remediation:
-  - Add MATLAB-derived fixtures for partitionNST, burst/statistics outputs, and ISI
-    plotting traces.
+  - Extend the committed MATLAB-derived fixtures beyond getSigRep, partitionNST, and
+    burst-stat summaries to cover ISI plotting traces and the remaining visualization
+    details.
   plotting_report_parity: Raster, ISI, and burst-oriented plotting helpers now execute
     on the canonical class, though visual detail remains lighter than MATLAB.
 - matlab_name: nstColl
@@ -331,10 +333,14 @@ items:
     in the expected workflow positions.
   known_remaining_differences:
   - Simulink-backed recursive-CIF behavior is represented by a native Python implementation,
-    but it is not yet fixture-matched one-for-one against MATLAB/Simulink outputs.
+    but it is not yet fixture-matched one-for-one against MATLAB/Simulink stochastic
+    trajectories.
   required_remediation:
-  - Add MATLAB-derived fixtures for CIF evaluation and thinning outputs.
-  - Add MATLAB/Simulink comparison fixtures for recursive CIF simulation semantics.
+  - Extend the committed MATLAB-derived fixtures beyond analytic lambda/gradient/Jacobian
+    outputs and the deterministic recursive lambda prefix to cover additional thinning
+    and seeded simulation summaries.
+  - Add MATLAB/Simulink comparison fixtures for recursive CIF simulation trajectories
+    when the random-stream alignment question is resolved.
   plotting_report_parity: Simulation/report plotting is limited; downstream notebooks
     generate figures with helper code rather than a full MATLAB-equivalent CIF report
     API.
diff --git a/parity/notebook_fidelity.yml b/parity/notebook_fidelity.yml
index 46b2c474..e44ee7f7 100644
--- a/parity/notebook_fidelity.yml
+++ b/parity/notebook_fidelity.yml
@@ -155,9 +155,10 @@ items:
   python_notebook: notebooks/HippocampalPlaceCellExample.ipynb
   fidelity_status: high_fidelity
   remaining_differences: The notebook now reproduces the dataset-backed place-cell
-    model-comparison and field-visualization workflow with real figures; the Python
-    port still uses an approximate Zernike-like basis rather than the original MATLAB
-    toolbox implementation.
+    model-comparison and field-visualization workflow with the same normalized 10-term
+    Zernike basis used by MATLAB; exact AIC/BIC values and surface styling still vary
+    modestly because the Python GLM solver and plotting backend are not byte-identical
+    to MATLAB.
   python_sections: 5
   python_expected_figures: 11
   python_uses_figure_tracker: true
diff --git a/pyproject.toml b/pyproject.toml
index 5485d16c..86ea335b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -38,3 +38,4 @@ include = ["nstat*"]
 [project.scripts]
 nstat-paper-examples = "nstat.paper_examples:main"
 nstat-install = "nstat.install:main"
+nstat-release-check = "nstat.release_check:main"
diff --git a/tests/parity/fixtures/matlab_gold/cif_exactness.mat b/tests/parity/fixtures/matlab_gold/cif_exactness.mat
new file mode 100644
index 00000000..a9b37b02
Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/cif_exactness.mat differ
diff --git a/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat b/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat
new file mode 100644
index 00000000..64522508
Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat differ
diff --git a/tests/parity/fixtures/matlab_gold/point_process_exactness.mat b/tests/parity/fixtures/matlab_gold/point_process_exactness.mat
new file mode 100644
index 00000000..840c58d9
Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/point_process_exactness.mat differ
diff --git a/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat b/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat
new file mode 100644
index 00000000..d931cfe1
Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat differ
diff --git a/tests/test_matlab_gold_fixtures.py b/tests/test_matlab_gold_fixtures.py
new file mode 100644
index 00000000..6693a368
--- /dev/null
+++ b/tests/test_matlab_gold_fixtures.py
@@ -0,0 +1,112 @@
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+from scipy.io import loadmat
+
+from nstat import CIF, Covariate, SignalObj, nspikeTrain
+
+
+REPO_ROOT = Path(__file__).resolve().parents[1]
+FIXTURE_ROOT = REPO_ROOT / "tests" / "parity" / "fixtures" / "matlab_gold"
+
+
+def _load_fixture(name: str) -> dict[str, np.ndarray]:
+    return loadmat(FIXTURE_ROOT / name, squeeze_me=True, struct_as_record=False)
+
+
+def _scalar(payload: dict[str, np.ndarray], key: str) -> float:
+    return float(np.asarray(payload[key], dtype=float).reshape(-1)[0])
+
+
+def _vector(payload: dict[str, np.ndarray], key: str) -> np.ndarray:
+    return np.asarray(payload[key], dtype=float).reshape(-1)
+
+
+def test_signalobj_matches_matlab_gold_fixture() -> None:
+    payload = _load_fixture("signalobj_exactness.mat")
+    signal = SignalObj(_vector(payload, "time"), np.asarray(payload["data"], dtype=float), "sig", "time", "s", "u", ["x1", "x2"])
+
+    filtered = signal.filter(_vector(payload, "filter_b"), _vector(payload, "filter_a"))
+    derivative = signal.derivative
+    integral = signal.integral()
+    resampled = signal.resample(_scalar(payload, "resample_rate"))
+    xcorr = signal.getSubSignal(1).xcorr(signal.getSubSignal(2), int(_scalar(payload, "xcorr_maxlag")))
+
+    np.testing.assert_allclose(filtered.data, np.asarray(payload["filtered_data"], dtype=float), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(derivative.data, np.asarray(payload["derivative_data"], dtype=float), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(integral.data, np.asarray(payload["integral_data"], dtype=float), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(resampled.time, _vector(payload, "resampled_time"), rtol=1e-12, atol=1e-12)
+    np.testing.assert_allclose(resampled.data, np.asarray(payload["resampled_data"], dtype=float), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(xcorr.time, _vector(payload, "xcorr_time"), rtol=1e-12, atol=1e-12)
+    np.testing.assert_allclose(xcorr.data.reshape(-1), _vector(payload, "xcorr_data"), rtol=1e-8, atol=1e-10)
+
+
+def test_nspiketrain_matches_matlab_gold_fixture() -> None:
+    payload = _load_fixture("nspiketrain_exactness.mat")
+    nst = nspikeTrain(
+        _vector(payload, "spikeTimes"),
+        "nst",
+        _scalar(payload, "binwidth"),
+        _scalar(payload, "minTime"),
+        _scalar(payload, "maxTime"),
+        "time",
+        "s",
+        "spikes",
+        "spk",
+        0,
+    )
+
+    sig = nst.getSigRep(_scalar(payload, "binwidth"), _scalar(payload, "minTime"), _scalar(payload, "maxTime"))
+    parts = nst.partitionNST([_scalar(payload, "minTime"), 0.2, _scalar(payload, "maxTime")])
+
+    np.testing.assert_allclose(sig.time, _vector(payload, "sig_time"), rtol=1e-12, atol=1e-12)
+    np.testing.assert_allclose(sig.data[:, 0], _vector(payload, "sig_data"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(nst.getISIs(), _vector(payload, "isis"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(float(nst.avgFiringRate), _scalar(payload, "avgFiringRate"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(float(nst.B), _scalar(payload, "B"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(float(nst.An), _scalar(payload, "An"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(float(nst.burstIndex), _scalar(payload, "burstIndex"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(parts.getNST(1).spikeTimes, _vector(payload, "part1_spikes"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(parts.getNST(2).spikeTimes, _vector(payload, "part2_spikes"), rtol=1e-8, atol=1e-10)
+
+
+def test_cif_eval_surface_matches_matlab_gold_fixture() -> None:
+    payload = _load_fixture("cif_exactness.mat")
+    cif = CIF(
+        beta=_vector(payload, "beta"),
+        Xnames=["stim1", "stim2"],
+        stimNames=["stim1", "stim2"],
+        fitType="binomial",
+    )
+
+    stim_val = _vector(payload, "stimVal")
+
+    np.testing.assert_allclose(cif.evalLambdaDelta(stim_val), _scalar(payload, "lambda_delta"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(cif.evalGradient(stim_val).reshape(-1), _vector(payload, "gradient"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(cif.evalGradientLog(stim_val).reshape(-1), _vector(payload, "gradient_log"), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(cif.evalJacobian(stim_val), np.asarray(payload["jacobian"], dtype=float), rtol=1e-8, atol=1e-10)
+    np.testing.assert_allclose(cif.evalJacobianLog(stim_val), np.asarray(payload["jacobian_log"], dtype=float), rtol=1e-8, atol=1e-10)
+
+
+def test_point_process_lambda_trace_matches_matlab_gold_fixture() -> None:
+    payload = _load_fixture("point_process_exactness.mat")
+
+    time = np.arange(0.0, 50.0 + 0.001, 0.001, dtype=float)
+    stim = Covariate(time, np.sin(2 * np.pi * 1.0 * time), "Stimulus", "time", "s", "Voltage", ["sin"])
+    ens = Covariate(time, np.zeros_like(time), "Ensemble", "time", "s", "Spikes", ["n1"])
+    _, lambda_cov = CIF.simulateCIF(
+        -3.0,
+        np.array([-1.0, -2.0, -4.0], dtype=float),
+        np.array([1.0], dtype=float),
+        np.array([0.0], dtype=float),
+        stim,
+        ens,
+        numRealizations=5,
+        simType="binomial",
+        seed=int(_scalar(payload, "seed")),
+        return_lambda=True,
+    )
+
+    np.testing.assert_allclose(lambda_cov.data[: _vector(payload, 'lambda_head').shape[0], 0], _vector(payload, "lambda_head"), rtol=1e-8, atol=1e-10)
diff --git a/tests/test_release_check.py b/tests/test_release_check.py
new file mode 100644
index 00000000..7d91b856
--- /dev/null
+++ b/tests/test_release_check.py
@@ -0,0 +1,24 @@
+from __future__ import annotations
+
+from pathlib import Path
+
+from nstat.release_check import build_release_gate_commands
+
+
+REPO_ROOT = Path(__file__).resolve().parents[1]
+
+
+def test_release_gate_includes_fixture_generation_and_matlab_suite() -> None:
+    commands = build_release_gate_commands(REPO_ROOT, matlab_repo_root=REPO_ROOT.parent / "nSTAT")
+    flattened = [" ".join(command) for command in commands]
+
+    assert any("export_matlab_gold_fixtures.py" in item for item in flattened)
+    assert any("tests/python_port_fidelity" in item for item in flattened)
+    assert any("addpath(fullfile(pwd,'helpfiles'))" in item for item in flattened)
+
+
+def test_release_gate_can_skip_matlab() -> None:
+    commands = build_release_gate_commands(REPO_ROOT, skip_matlab=True)
+    flattened = [" ".join(command) for command in commands]
+
+    assert not any(item.startswith("matlab ") for item in flattened)
diff --git a/tests/test_zernike.py b/tests/test_zernike.py
new file mode 100644
index 00000000..09e6c6bf
--- /dev/null
+++ b/tests/test_zernike.py
@@ -0,0 +1,50 @@
+from __future__ import annotations
+
+import numpy as np
+
+from nstat.zernike import zernfun, zernike_basis_from_cartesian
+
+
+def test_normalized_zernfun_matches_closed_form_first_ten_modes() -> None:
+    x = np.array([0.3, -0.5], dtype=float)
+    y = np.array([0.4, 0.25], dtype=float)
+    r = np.sqrt(x * x + y * y)
+    theta = np.arctan2(y, x)
+
+    basis = zernike_basis_from_cartesian(x, y)
+    expected = np.column_stack(
+        [
+            np.ones_like(r),
+            2.0 * r * np.cos(theta),
+            2.0 * r * np.sin(theta),
+            np.sqrt(6.0) * (r**2) * np.cos(2.0 * theta),
+            np.sqrt(3.0) * (2.0 * r**2 - 1.0),
+            np.sqrt(6.0) * (r**2) * np.sin(2.0 * theta),
+            np.sqrt(8.0) * (r**3) * np.cos(3.0 * theta),
+            np.sqrt(8.0) * (3.0 * r**3 - 2.0 * r) * np.cos(theta),
+            np.sqrt(8.0) * (3.0 * r**3 - 2.0 * r) * np.sin(theta),
+            np.sqrt(8.0) * (r**3) * np.sin(3.0 * theta),
+        ]
+    )
+
+    np.testing.assert_allclose(basis, expected, rtol=1e-12, atol=1e-12)
+
+
+def test_zernfun_matches_matlab_sign_convention_for_positive_and_negative_m() -> None:
+    r = np.array([0.2, 0.7], dtype=float)
+    theta = np.array([0.1, 1.2], dtype=float)
+
+    values = zernfun(np.array([1, 1], dtype=int), np.array([-1, 1], dtype=int), r, theta, normalized=True)
+
+    np.testing.assert_allclose(values[:, 0], 2.0 * r * np.cos(theta), rtol=1e-12, atol=1e-12)
+    np.testing.assert_allclose(values[:, 1], 2.0 * r * np.sin(theta), rtol=1e-12, atol=1e-12)
+
+
+def test_zernike_basis_from_cartesian_can_mask_outside_unit_disk() -> None:
+    x = np.array([0.0, 1.1], dtype=float)
+    y = np.array([0.0, 0.0], dtype=float)
+
+    basis = zernike_basis_from_cartesian(x, y, fill_value=np.nan)
+
+    np.testing.assert_allclose(basis[0, 0], 1.0, rtol=1e-12, atol=1e-12)
+    assert np.isnan(basis[1]).all()
diff --git a/tools/notebooks/build_helpfile_fidelity_notebooks.py b/tools/notebooks/build_helpfile_fidelity_notebooks.py
index da1260c3..0bc53d8a 100644
--- a/tools/notebooks/build_helpfile_fidelity_notebooks.py
+++ b/tools/notebooks/build_helpfile_fidelity_notebooks.py
@@ -554,7 +554,7 @@ def _plot_isi_hist(ax, train, lambda_hz, *, title):
 ## MATLAB Parity Note
 - Source MATLAB helpfile: `HippocampalPlaceCellExample.mlx`
 - Fidelity status: `high_fidelity`
-- Remaining justified differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with real figures; the Python port still uses an approximate Zernike-like basis rather than the original MATLAB toolbox implementation.
+- Remaining justified differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with the same normalized 10-term Zernike basis used by MATLAB; exact AIC/BIC values and surface styling still vary modestly because the Python GLM solver and plotting backend are not byte-identical to MATLAB.
 """
 
 
@@ -753,7 +753,7 @@ def _plot_field_grid(fig, animal_key, field_key, title):
                 f"Cells analyzed: {int(summary['num_cells_fit'])}",
                 f"Mean Gaussian-Zernike ΔAIC: {summary['mean_delta_aic_gaussian_minus_zernike']:.2f}",
                 f"Mean Gaussian-Zernike ΔBIC: {summary['mean_delta_bic_gaussian_minus_zernike']:.2f}",
-                "Negative values favor the Zernike-like model.",
+                "Negative values favor the Zernike model.",
             ]
         ),
         va="top",
diff --git a/tools/notebooks/build_nstat_paper_notebook.py b/tools/notebooks/build_nstat_paper_notebook.py
index 7805adf7..ce5ffeb8 100644
--- a/tools/notebooks/build_nstat_paper_notebook.py
+++ b/tools/notebooks/build_nstat_paper_notebook.py
@@ -335,11 +335,11 @@ def _fig(label: str, *, figsize=(8.5, 4.5)):
     ax.set_title("Gaussian place-field estimate")
     """,
     """
-    # SECTION 26: Zernike-like place-field mesh
+    # SECTION 26: Zernike place-field mesh
     fig = _fig("experiment4 zernike mesh", figsize=(9.0, 6.5))
     ax = fig.add_subplot(111, projection="3d")
     ax.plot_surface(exp4["mesh"]["grid_x"], exp4["mesh"]["grid_y"], exp4["mesh"]["zernike_field"], cmap="Greens", linewidth=0.0, antialiased=True)
-    ax.set_title("Zernike-like place-field estimate")
+    ax.set_title("Zernike place-field estimate")
     """,
     """
     # SECTION 27: Experiment 5
diff --git a/tools/notebooks/parity_notes.yml b/tools/notebooks/parity_notes.yml
index f0601f62..06a43428 100644
--- a/tools/notebooks/parity_notes.yml
+++ b/tools/notebooks/parity_notes.yml
@@ -39,7 +39,7 @@ notes:
     file: notebooks/HippocampalPlaceCellExample.ipynb
     source_matlab: HippocampalPlaceCellExample.mlx
     fidelity_status: high_fidelity
-    remaining_differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with real figures; the Python port still uses an approximate Zernike-like basis rather than the original MATLAB toolbox implementation.
+    remaining_differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with the same normalized 10-term Zernike basis used by MATLAB; exact AIC/BIC values and surface styling still vary modestly because the Python GLM solver and plotting backend are not byte-identical to MATLAB.
   - topic: HybridFilterExample
     file: notebooks/HybridFilterExample.ipynb
     source_matlab: HybridFilterExample.mlx
diff --git a/tools/parity/export_matlab_gold_fixtures.py b/tools/parity/export_matlab_gold_fixtures.py
new file mode 100644
index 00000000..8313279f
--- /dev/null
+++ b/tools/parity/export_matlab_gold_fixtures.py
@@ -0,0 +1,51 @@
+from __future__ import annotations
+
+import argparse
+import shutil
+import subprocess
+import sys
+from pathlib import Path
+
+
+THIS_DIR = Path(__file__).resolve().parent
+REPO_ROOT = THIS_DIR.parents[1]
+HELPER_DIR = THIS_DIR / "matlab"
+
+
+def _matlab_quote(value: str) -> str:
+    return value.replace("'", "''")
+
+
+def default_matlab_repo(repo_root: Path | None = None) -> Path:
+    base = REPO_ROOT if repo_root is None else repo_root.resolve()
+    return base.parent / "nSTAT"
+
+
+def build_matlab_batch(repo_root: Path, matlab_repo: Path) -> str:
+    return (
+        f"addpath('{_matlab_quote(str(HELPER_DIR))}'); "
+        f"export_matlab_gold_fixtures('{_matlab_quote(str(repo_root))}','{_matlab_quote(str(matlab_repo))}'); "
+        "exit"
+    )
+
+
+def main(argv: list[str] | None = None) -> int:
+    parser = argparse.ArgumentParser(description="Generate MATLAB-derived gold fixtures for Python parity tests.")
+    parser.add_argument("--repo-root", type=Path, default=REPO_ROOT)
+    parser.add_argument("--matlab-repo", type=Path, default=None)
+    args = parser.parse_args(argv)
+
+    repo_root = args.repo_root.resolve()
+    matlab_repo = default_matlab_repo(repo_root) if args.matlab_repo is None else args.matlab_repo.resolve()
+    if not matlab_repo.exists():
+        raise FileNotFoundError(f"MATLAB reference repo not found at {matlab_repo}")
+    if shutil.which("matlab") is None:
+        raise FileNotFoundError("`matlab` is not on PATH")
+
+    command = ["matlab", "-batch", build_matlab_batch(repo_root, matlab_repo)]
+    subprocess.run(command, cwd=repo_root, check=True)
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv[1:]))
diff --git a/tools/parity/matlab/export_matlab_gold_fixtures.m b/tools/parity/matlab/export_matlab_gold_fixtures.m
new file mode 100644
index 00000000..ac44f983
--- /dev/null
+++ b/tools/parity/matlab/export_matlab_gold_fixtures.m
@@ -0,0 +1,122 @@
+function export_matlab_gold_fixtures(repoRoot, matlabRepoRoot)
+if nargin < 1 || isempty(repoRoot)
+    error('repoRoot is required');
+end
+if nargin < 2 || isempty(matlabRepoRoot)
+    matlabRepoRoot = fullfile(fileparts(repoRoot), 'nSTAT');
+end
+
+repoRoot = char(repoRoot);
+matlabRepoRoot = char(matlabRepoRoot);
+
+addpath(matlabRepoRoot);
+addpath(fullfile(matlabRepoRoot, 'helpfiles'));
+addpath(genpath(fullfile(matlabRepoRoot, 'libraries')));
+
+fixtureRoot = fullfile(repoRoot, 'tests', 'parity', 'fixtures', 'matlab_gold');
+if ~exist(fixtureRoot, 'dir')
+    mkdir(fixtureRoot);
+end
+
+export_signalobj_fixture(fixtureRoot);
+export_nspiketrain_fixture(fixtureRoot);
+export_cif_fixture(fixtureRoot);
+export_point_process_fixture(fixtureRoot);
+end
+
+function export_signalobj_fixture(fixtureRoot)
+t = (0:0.1:0.4)';
+data = [1 0; 2 1; 4 0; 8 -1; 16 0];
+s = SignalObj(t, data, 'sig', 'time', 's', 'u', {'x1', 'x2'});
+
+filtered = s.filter([0.25 0.5 0.25], 1);
+resampled = s.resample(20);
+derivative = s.derivative;
+integral_sig = s.integral();
+xc = xcorr(s.getSubSignal(1), s.getSubSignal(2), 2);
+
+payload = struct();
+payload.time = s.time;
+payload.data = s.data;
+payload.filter_b = [0.25 0.5 0.25];
+payload.filter_a = 1;
+payload.filtered_data = filtered.data;
+payload.derivative_data = derivative.data;
+payload.integral_data = integral_sig.data;
+payload.resample_rate = 20;
+payload.resampled_time = resampled.time;
+payload.resampled_data = resampled.data;
+payload.xcorr_maxlag = 2;
+payload.xcorr_time = xc.time;
+payload.xcorr_data = xc.data;
+
+save(fullfile(fixtureRoot, 'signalobj_exactness.mat'), '-struct', 'payload');
+end
+
+function export_nspiketrain_fixture(fixtureRoot)
+spikeTimes = [0.05; 0.10; 0.11; 0.30; 0.47];
+binwidth = 0.05;
+nst = nspikeTrain(spikeTimes, 'nst', binwidth, 0.0, 0.5, 'time', 's', 'spikes', 'spk', 0);
+sig = nst.getSigRep(binwidth, 0.0, 0.5);
+parts = nst.partitionNST([0.0 0.2 0.5]);
+
+payload = struct();
+payload.spikeTimes = spikeTimes;
+payload.binwidth = binwidth;
+payload.minTime = 0.0;
+payload.maxTime = 0.5;
+payload.sig_time = sig.time;
+payload.sig_data = sig.data;
+payload.isis = nst.getISIs();
+payload.avgFiringRate = nst.avgFiringRate;
+payload.B = nst.B;
+payload.An = nst.An;
+payload.burstIndex = nst.burstIndex;
+payload.numBursts = nst.numBursts;
+payload.numSpikesPerBurst = nst.numSpikesPerBurst;
+payload.part1_spikes = parts.getNST(1).spikeTimes;
+payload.part2_spikes = parts.getNST(2).spikeTimes;
+
+save(fullfile(fixtureRoot, 'nspiketrain_exactness.mat'), '-struct', 'payload');
+end
+
+function export_cif_fixture(fixtureRoot)
+cif = CIF([0.1 0.5], {'stim1', 'stim2'}, {'stim1', 'stim2'}, 'binomial');
+stimVal = [0.6; -0.2];
+
+payload = struct();
+payload.beta = [0.1 0.5];
+payload.stimVal = stimVal;
+payload.lambda_delta = cif.evalLambdaDelta(stimVal);
+payload.gradient = cif.evalGradient(stimVal);
+payload.gradient_log = cif.evalGradientLog(stimVal);
+payload.jacobian = cif.evalJacobian(stimVal);
+payload.jacobian_log = cif.evalJacobianLog(stimVal);
+
+save(fullfile(fixtureRoot, 'cif_exactness.mat'), '-struct', 'payload');
+end
+
+function export_point_process_fixture(fixtureRoot)
+rng(5);
+Ts = .001;
+t = (0:Ts:50)';
+mu = -3;
+H = tf([-1 -2 -4], [1], Ts, 'Variable', 'z^-1');
+S = tf([1], 1, Ts, 'Variable', 'z^-1');
+E = tf([0], 1, Ts, 'Variable', 'z^-1');
+stim = Covariate(t, sin(2*pi*1*t), 'Stimulus', 'time', 's', 'Voltage', {'sin'});
+ens = Covariate(t, zeros(length(t), 1), 'Ensemble', 'time', 's', 'Spikes', {'n1'});
+[spikeColl, lambda] = CIF.simulateCIF(mu, H, S, E, stim, ens, 5, 'binomial');
+
+spikeCounts = zeros(1, spikeColl.numSpikeTrains);
+for i = 1:spikeColl.numSpikeTrains
+    spikeCounts(i) = length(spikeColl.getNST(i).spikeTimes);
+end
+
+payload = struct();
+payload.seed = 5;
+payload.lambda_head = lambda.data(1:8, 1);
+payload.spike_counts = spikeCounts;
+
+save(fullfile(fixtureRoot, 'point_process_exactness.mat'), '-struct', 'payload');
+end
diff --git a/tools/release/run_fidelity_gate.py b/tools/release/run_fidelity_gate.py
new file mode 100644
index 00000000..7121d0c7
--- /dev/null
+++ b/tools/release/run_fidelity_gate.py
@@ -0,0 +1,16 @@
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+
+THIS_DIR = Path(__file__).resolve().parent
+REPO_ROOT = THIS_DIR.parents[1]
+if str(REPO_ROOT) not in sys.path:
+    sys.path.insert(0, str(REPO_ROOT))
+
+from nstat.release_check import main
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv[1:]))