diff --git a/docs/changes/devel/42.feature.rst b/docs/changes/devel/42.feature.rst
new file mode 100644
index 00000000..c93a6a1e
--- /dev/null
+++ b/docs/changes/devel/42.feature.rst
@@ -0,0 +1,2 @@
+**DSS**:
+- ``DSS`` now estimates covariances for :class:`mne.Evoked` inputs with :func:`mne.compute_covariance` instead of falling back to the NumPy path (Issue #39).
diff --git a/mne_denoise/dss/linear.py b/mne_denoise/dss/linear.py
index 9cf7e5e0..c96ba51b 100644
--- a/mne_denoise/dss/linear.py
+++ b/mne_denoise/dss/linear.py
@@ -16,6 +16,7 @@
 from __future__ import annotations
 
 import logging
+import warnings
 from collections.abc import Callable
 
 import numpy as np
@@ -753,23 +754,35 @@ def _fit_mne(
 
         biased_data = self._apply_bias(data)
 
-        if isinstance(inst, BaseEpochs):
-            biased_data = np.transpose(biased_data, (2, 0, 1))
-
         if isinstance(inst, BaseRaw):
             kws.setdefault("tstep", 2.0)
             baseline_cov = mne.compute_raw_covariance(inst, method=method, **kws)
             biased_inst = mne.io.RawArray(biased_data, inst.info, verbose=False)
             biased_cov = mne.compute_raw_covariance(biased_inst, method=method, **kws)
 
-        elif isinstance(inst, BaseEpochs):
-            baseline_cov = mne.compute_covariance(inst, method=method, **kws)
-            biased_inst = mne.EpochsArray(biased_data, inst.info, verbose=False)
-            biased_cov = mne.compute_covariance(biased_inst, method=method, **kws)
-
-        else:  # Evoked - use numpy path since MNE doesn't support Evoked covariance
-            self._fit_numpy(data, weights=weights)
-            return
+        else:
+            if isinstance(inst, BaseEpochs):
+                biased_data = np.transpose(biased_data, (2, 0, 1))
+                biased_inst = mne.EpochsArray(biased_data, inst.info, verbose=False)
+            else:  # Evoked
+                biased_inst = mne.EvokedArray(
+                    biased_data,
+                    inst.info,
+                    tmin=float(inst.times[0]),
+                    verbose=False,
+                )
+            with warnings.catch_warnings():
+                if isinstance(inst, Evoked):
+                    # DSS intentionally uses the raw second moment here.
+                    warnings.filterwarnings(
+                        "ignore",
+                        message="Evoked is not baseline corrected",
+                        category=RuntimeWarning,
+                    )
+                baseline_cov = mne.compute_covariance(inst=inst, method=method, **kws)
+                biased_cov = mne.compute_covariance(
+                    inst=biased_inst, method=method, **kws
+                )
 
         # Extract data from MNE covariances
         self.filters_, self.patterns_, self.eigenvalues_ = compute_dss(
diff --git a/pyproject.toml b/pyproject.toml
index 36f8e806..7c5ccd6e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -45,7 +45,7 @@ classifiers = [
   "Topic :: Scientific/Engineering :: Medical Science Apps.",
 ]
 dependencies = [
-  "mne>=1.9",
+  "mne>=1.13",
   "numpy>=1.22",
   "scipy>=1.8",
   "matplotlib>=3.5",
@@ -78,7 +78,7 @@ dev = [
   "towncrier",
 ]
 docs = [
-  "mne>=1.9",
+  "mne>=1.13",
   "numpydoc>=1.6.0",
   "pydata-sphinx-theme>=0.15.2",
   "sphinx>=7.2.0",
diff --git a/tests/test_linear_dss.py b/tests/test_linear_dss.py
index cd28deaa..63842e20 100644
--- a/tests/test_linear_dss.py
+++ b/tests/test_linear_dss.py
@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import warnings
+from unittest.mock import patch
 
 import mne
 import numpy as np
@@ -817,7 +818,20 @@ def test_dss_mne_evoked_extracts_known_signal():
     bias = BandpassBias(freq_band=(8, 12), sfreq=sfreq)
 
     dss = DSS(bias=bias, n_components=1, normalize_input=False)
-    sources = dss.fit_transform(evoked)
+    with (
+        patch.object(
+            mne, "compute_covariance", wraps=mne.compute_covariance
+        ) as compute_covariance,
+        warnings.catch_warnings(),
+    ):
+        warnings.filterwarnings(
+            "error",
+            message="Evoked is not baseline corrected",
+            category=RuntimeWarning,
+        )
+        sources = dss.fit_transform(evoked)
+
+    assert compute_covariance.call_count == 2
 
     # Top source should correlate with signal
     correlation = np.abs(np.corrcoef(sources[0], signal)[0, 1])