diff --git a/lumicks/pylake/force_calibration/power_spectrum.py b/lumicks/pylake/force_calibration/power_spectrum.py
index 4934f6685..186a48d96 100644
--- a/lumicks/pylake/force_calibration/power_spectrum.py
+++ b/lumicks/pylake/force_calibration/power_spectrum.py
@@ -193,7 +193,34 @@ def _variance(self):
         return self._apply_transforms(self._raw_variance)
 
     def downsampled_by(self, factor, reduce=np.mean) -> "PowerSpectrum":
-        """Returns a spectrum downsampled by a given factor."""
+        """Returns a spectrum downsampled by a given factor.
+
+        Parameters
+        ----------
+        factor : int
+            Factor to down-sample the spectrum by.
+        reduce : callable
+            (Deprecated) Function to use for down-sampling the data. Only `np.mean` will be
+            supported going forward.
+        """
+        if reduce != np.mean:
+            warnings.warn(
+                DeprecationWarning(
+                    "Providing other reduction functions than `np.mean` is deprecated and will be "
+                    "removed in a future version of Pylake"
+                )
+            )
+            return PowerSpectrum(
+                downsample(self.frequency, factor, reduce),
+                downsample(self.power, factor, reduce),
+                self.sample_rate,
+                self.total_duration,
+                self.unit,
+                num_points_per_block=self.num_points_per_block * factor,
+                total_samples_used=self.total_sampled_used,
+                variance=None,
+            )
+
         ba = copy(self)
         ba._downsampling_factor = ba._downsampling_factor * factor
 
diff --git a/lumicks/pylake/force_calibration/tests/test_power_spectrum.py b/lumicks/pylake/force_calibration/tests/test_power_spectrum.py
index 372ccf115..309d327d8 100644
--- a/lumicks/pylake/force_calibration/tests/test_power_spectrum.py
+++ b/lumicks/pylake/force_calibration/tests/test_power_spectrum.py
@@ -84,6 +84,9 @@ def test_power_spectrum_attrs(frequency, num_data, sample_rate):
     np.testing.assert_allclose(power_spectrum.frequency_bin_width, sample_rate / num_data)
 
 
+@pytest.mark.filterwarnings(
+    "ignore:Providing other reduction functions than `np.mean` is deprecated*"
+)
 @pytest.mark.parametrize(
     "frequency, num_data, sample_rate, num_blocks, f_blocked, p_blocked",
     [
@@ -100,28 +103,35 @@ def test_power_spectrum_blocking(
     frequency, num_data, sample_rate, num_blocks, f_blocked, p_blocked
 ):
     """Functional test whether the results of blocking the power spectrum are correct"""
-    data = np.sin(2.0 * np.pi * frequency / sample_rate * np.arange(num_data)) / np.sqrt(2)
-    power_spectrum = PowerSpectrum.from_data(data, sample_rate)
 
-    downsampling_factor = len(power_spectrum.frequency) // num_blocks
-    blocked = power_spectrum.downsampled_by(downsampling_factor)
-    np.testing.assert_allclose(blocked.frequency, f_blocked)
-    np.testing.assert_allclose(blocked.power, p_blocked, atol=1e-16)
-    np.testing.assert_allclose(blocked.num_samples(), len(f_blocked))
-    np.testing.assert_allclose(len(power_spectrum.power), num_data // 2 + 1)
-    np.testing.assert_allclose(
-        blocked.num_points_per_block, np.floor(len(power_spectrum.power) / num_blocks)
-    )
-    np.testing.assert_allclose(
-        blocked.frequency_bin_width, blocked.frequency[1] - blocked.frequency[0]
-    )
+    def custom_reduce(x, *args, **kwargs):
+        return np.mean(x, *args, **kwargs)
+
+    for reduce in (np.mean, custom_reduce):
+        data = np.sin(2.0 * np.pi * frequency / sample_rate * np.arange(num_data)) / np.sqrt(2)
+        power_spectrum = PowerSpectrum.from_data(data, sample_rate)
+
+        downsampling_factor = len(power_spectrum.frequency) // num_blocks
+        blocked = power_spectrum.downsampled_by(downsampling_factor, reduce=reduce)
+        np.testing.assert_allclose(blocked.frequency, f_blocked)
+        np.testing.assert_allclose(blocked.power, p_blocked, atol=1e-16)
+        np.testing.assert_allclose(blocked.num_samples(), len(f_blocked))
+        np.testing.assert_allclose(len(power_spectrum.power), num_data // 2 + 1)
+        np.testing.assert_allclose(
+            blocked.num_points_per_block, np.floor(len(power_spectrum.power) / num_blocks)
+        )
+        np.testing.assert_allclose(
+            blocked.frequency_bin_width, blocked.frequency[1] - blocked.frequency[0]
+        )
 
-    # Downsample again and make sure the num_points_per_block is correct
-    dual_blocked = blocked.downsampled_by(2)
-    np.testing.assert_allclose(dual_blocked.num_points_per_block, blocked.num_points_per_block * 2)
-    np.testing.assert_allclose(
-        dual_blocked.frequency_bin_width, dual_blocked.frequency[1] - dual_blocked.frequency[0]
-    )
+        # Downsample again and make sure the num_points_per_block is correct
+        dual_blocked = blocked.downsampled_by(2, reduce=reduce)
+        np.testing.assert_allclose(
+            dual_blocked.num_points_per_block, blocked.num_points_per_block * 2
+        )
+        np.testing.assert_allclose(
+            dual_blocked.frequency_bin_width, dual_blocked.frequency[1] - dual_blocked.frequency[0]
+        )
 
 
 @pytest.mark.parametrize(