diff --git a/bayesian_listener/bayesian_listener.py b/bayesian_listener/bayesian_listener.py
index c577abd..afa01d6 100644
--- a/bayesian_listener/bayesian_listener.py
+++ b/bayesian_listener/bayesian_listener.py
@@ -485,40 +485,5 @@ def plot_post(self, posterior, estimations):
                                              np.log(np.finfo(amps.dtype).eps)),
                                              estimations.squeeze())
 
-def test_interp():
-    # doing it from scratch
-    sofa_file = 'data/P0001_FreeFieldCompMinPhase_48kHz.sofa'
-    am = BayesianListener(sofa_file)
-    am.prepare_features(use_cache=False, interpolation='SHMAX')
-
-    # Get target spectral cues to compute color limits
-    side = 0
-    dirs = am.coords.sph()
-    median_idx = np.abs(dirs[:, 0] - 0) < 2
-    elevations = dirs[median_idx, 1]
-    amps_target = am.spectral_cues[median_idx, :, side]
-    sorted_indices = np.argsort(elevations)
-    amps_target = amps_target[sorted_indices, :]
-
-    # Compute color limits from target
-    clim = (np.min(amps_target), np.max(amps_target))
-
-    # Create figure with two subplots
-    fig, axes = plt.subplots(1, 2, figsize=(16, 6))
-
-    # Plot target and template features side by side with same color limits
-    am.plot_cues('- Target features (i.e. original)',
-                 fig=fig,
-                 ax=axes[0],
-                 clim=clim)
-    am.template.plot_cues('- Template features (i.e. interpolated)',
-                          fig=fig, ax=axes[1],
-                          clim=clim,
-                          elev_min=-45)
-
-    plt.tight_layout()
-    plt.show()
-
-
 
 
diff --git a/tests/test_bayesian_listener.py b/tests/test_bayesian_listener.py
index 8e90b21..c99cf4a 100644
--- a/tests/test_bayesian_listener.py
+++ b/tests/test_bayesian_listener.py
@@ -68,4 +68,49 @@ def test_model_single():
 
     # Compare with fixed expected spherical coordinates (azimuth, elevation)
     expected_dir_sph = np.array([[126.58887 ,  -9.108036]])
-    np.testing.assert_allclose(estimated_dir, expected_dir_sph, rtol=1e-2)
\ No newline at end of file
+    np.testing.assert_allclose(estimated_dir, expected_dir_sph, rtol=1e-2)
+
+def test_interp():
+    """Test SHMAX interpolation produces valid template features."""
+    sofa_file = get_sofa_file()
+    am = BayesianListener(sofa_file)
+    am.prepare_features()
+
+    # Verify template was created
+    assert hasattr(am, 'template'), "Template should be created after prepare_features"
+    assert am.template is not None, "Template should not be None"
+
+    # Verify template has spectral cues
+    assert hasattr(am.template, 'spectral_cues'), "Template should have spectral_cues"
+    assert am.template.spectral_cues is not None, "Template spectral_cues should not be None"
+
+    # Verify shapes match expected dimensions
+    # spectral_cues shape should be (n_directions, n_frequencies, n_sides)
+    assert am.spectral_cues.ndim == 3, "Spectral cues should be 3D array"
+    assert am.template.spectral_cues.ndim == 3, "Template spectral cues should be 3D array"
+
+    # Template and original should have same frequency and side dimensions
+    assert am.spectral_cues.shape[1:] == am.template.spectral_cues.shape[1:], \
+        "Template should have same frequency and side dimensions as original"
+
+    # Verify interpolated values are reasonable (finite and within expected range)
+    assert np.all(np.isfinite(am.template.spectral_cues)), \
+        "Template spectral cues should all be finite"
+
+    # Get target spectral cues for comparison
+    side = 0
+    amps_target = am.spectral_cues[260, :, side]
+
+    # Verify template values are in similar range to original
+    template_min = np.min(am.template.spectral_cues[260, :, side])
+    template_max = np.max(am.template.spectral_cues[260, :, side])
+    target_min = np.min(amps_target)
+    target_max = np.max(amps_target)
+
+    # Template should be within a reasonable range of original data
+    # Allow some margin since interpolation might extend slightly
+    margin = 0.2 * (target_max - target_min)
+    assert template_min >= target_min - margin, \
+        "Template min should not be significantly below original min"
+    assert template_max <= target_max + margin, \
+        "Template max should not be significantly above original max"