diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index d499e5fc..16d0582e 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -14,7 +14,7 @@ jobs:
     strategy:
       max-parallel: 4
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
 
     steps:
     - uses: actions/checkout@v4
diff --git a/meegkit/utils/asr.py b/meegkit/utils/asr.py
index cab47ca9..6c5bf371 100755
--- a/meegkit/utils/asr.py
+++ b/meegkit/utils/asr.py
@@ -158,6 +158,11 @@ def fit_eeg_distribution(X, min_clean_fraction=0.25, max_dropout_fraction=0.1,
     # calculate the distribution's standard deviation from alpha and beta
     sig = np.sqrt((alpha ** 2) * gamma(3 / beta) / gamma(1 / beta))
 
+    # Ensure scalar values are returned (extract from arrays if needed)
+    alpha = float(np.asarray(alpha).squeeze())
+    mu = float(np.asarray(mu).squeeze())
+    sig = float(np.asarray(sig).squeeze())
+
     return mu, sig, alpha, beta
 
 
@@ -432,7 +437,8 @@ def numf(h, a, nb):
         toeplitz(impr, np.concatenate((1, np.zeros((1, nb))), axis=None)),
         h.T, rcond=None)[0].T
 
-    return b
+    # Ensure 1D array is returned
+    return np.atleast_1d(b.squeeze())
 
 
 def denf(R, na):
diff --git a/meegkit/utils/sig.py b/meegkit/utils/sig.py
index 4451407d..600fc654 100644
--- a/meegkit/utils/sig.py
+++ b/meegkit/utils/sig.py
@@ -612,6 +612,7 @@ def prony(h, nb, na):
 
     # K-M by N
     a_right = np.linalg.lstsq(-H2, h1, rcond=None)[0]
-    a = np.r_[(np.array([1]), a_right)][None, :]
+    a = np.r_[(np.array([1]), a_right)]
     b = np.dot(np.dot(c, a), H1.T)
-    return b, a
+    # Ensure 1D arrays are returned
+    return np.atleast_1d(b.squeeze()), np.atleast_1d(a.squeeze())
diff --git a/meegkit/utils/stats.py b/meegkit/utils/stats.py
index 830c6571..b5686c47 100644
--- a/meegkit/utils/stats.py
+++ b/meegkit/utils/stats.py
@@ -385,7 +385,7 @@ def snr_spectrum(X, freqs, n_avg=1, n_harm=1, skipbins=1):
 
             # Ratio
             with np.errstate(divide="ignore", invalid="ignore"):
-                SNR[i_bin, i_trial] = np.sqrt(A) / np.sqrt(B)
+                SNR[i_bin, i_trial] = np.sqrt(A) / np.sqrt(np.mean(B))
 
             del A
             del B
diff --git a/pyproject.toml b/pyproject.toml
index b92301dc..4f69a56a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ license = {text = "BSD (3-clause)"}
 dynamic = ["version", "dependencies"]
 description = "M/EEG denoising in Python"
 readme = {file = "README.md", content-type = "text/markdown"}
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 
 [project.urls]
 homepage = "https://nbara.github.io/python-meegkit"
diff --git a/tests/test_asr.py b/tests/test_asr.py
index f417c2e7..91f57502 100644
--- a/tests/test_asr.py
+++ b/tests/test_asr.py
@@ -189,7 +189,7 @@ def test_asr_class(method, reref, show=False):
 
     if reref:
         if method == "riemann":
-            with pytest.raises(ValueError, match="Add regularization"):
+            with pytest.raises(ValueError, match="add regularization"):
                 blah = ASR(method=method, estimator="scm")
                 blah.fit(raw2[:, train_idx])