diff --git a/lisainstrument/fixed_shift_numpy.py b/lisainstrument/fixed_shift_numpy.py
index 399ab162fecd2a4c1b5d1d5b37f0d69f63610a3d..966f41228acacc5897bc89e7f954d17295b613d8 100644
--- a/lisainstrument/fixed_shift_numpy.py
+++ b/lisainstrument/fixed_shift_numpy.py
@@ -381,3 +381,60 @@ def make_fixed_shift_lagrange_numpy(
     """
     fac = FixedShiftLagrange.factory(order)
     return FixedShiftNumpy(left_bound, right_bound, fac)
+
+
+class AdaptiveShiftNumpy:  # pylint: disable=too-few-public-methods
+    """Time shifting that accepts fixed or dynamic shifts as well as constant data
+
+    Instances  act as interpolator function, which delegates to interpolation
+    methods for fixed or dynamic time shifting. The specialised interpolation
+    functions are provided during construction.
+
+    In addition, instances store a fixed sample rate that is assumed for any
+    interpolated data, allowing time shifts to be given in time units.
+    """
+
+    def __init__(
+        self,
+        delay_const: Callable[[np.ndarray, float], np.ndarray],
+        delay_dynamic: Callable[[np.ndarray, np.ndarray], np.ndarray],
+        fsample: float,
+    ):
+        """Construct from fixed and dynamic interpolator functions and sample rate
+
+        Arguments:
+            delay_const: Interpolator function with same interface as FixedShiftNumpy
+            delay_dynamic: Interpolator function with same interface as DynamicShiftNumpy
+            fsample: Sample rate [Hz]
+        """
+        self._delay_const = delay_const
+        self._delay_dynamic = delay_dynamic
+        self._fsample = float(fsample)
+
+    def __call__(
+        self, x: np.ndarray | float, shift_time: np.ndarray | float
+    ) -> np.ndarray | float:
+        """Apply time shift to sequence, accepting scalars to represent constant arrays.
+
+        The shift is given in time units, and the data is assumed to be sampled
+        with rate fsample given in the constructor.
+
+        Both data and time shift can be scalar or 1D numpy arrays. Scalars
+        are interpreted as constant arrays. In case of scalar data, the same
+        scalar is returned. In case of scalar shift, a more efficient algorithm
+        is used, which should yield identical results as for a const shift array.
+
+        Args:
+            x: the data to be shifted, as scalar or 1D array
+            shift_time: time shift [s], as scalar or 1D array
+
+        Returns:
+            The shifted data
+        """
+
+        if isinstance(x, np.ndarray):
+            shift_samples = shift_time * self._fsample
+            if isinstance(shift_samples, np.ndarray):
+                return self._delay_dynamic(x, shift_samples)
+            return self._delay_const(x, float(shift_samples))
+        return float(x)
diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index e096cacc75afcb910fc995461819bc805666f958..301cc8336562e216325c81bc462f68577a3c5db1 100755
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -34,6 +34,7 @@ from .dynamic_delay_dsp import make_dynamic_shift_dsp_dask
 from .dynamic_delay_numpy import ShiftBC
 from .fixed_shift_dask import make_fixed_shift_lagrange_dask
 from .fixed_shift_dsp import make_fixed_shift_dsp_dask
+from .fixed_shift_numpy import AdaptiveShiftNumpy
 from .noises import generate_subseed
 
 logger = logging.getLogger(__name__)
@@ -639,34 +640,6 @@ class Instrument:
         """
         return numpyfy_dask_multi(op, chunks=self.chunking_size)
 
-    def apply_shift(self, x, shift_time):
-        """Apply a time shift using interpolators set up in init_delays()
-
-        The shift is given in time units, and the data is assumed to be sampled
-        with rate self.physics_fs.
-
-        Both data and time shift can be scalar or 1D numpy arrays. Scalars
-        are interpreted as constant arrays. In case of scalar data, the same
-        scalar is returned. In case of scalar shift, a more efficient algorithm
-        is used, which should yield identical results as for a const shift array.
-
-        The properties of the interpolation algorithms are set in init_delays(),
-        see there for details.
-
-        Args:
-            x: the data to be shifted, scalar or 1D array
-            shift_time: the shift in time units, scalar or 1D array
-
-        Returns:
-            The shifted data
-        """
-        if np.isscalar(x):
-            return x
-        shift_samples = shift_time * self.physics_fs
-        if np.isscalar(shift_time):
-            return self._delay_const(x, float(shift_samples))
-        return self._delay_dynamic(x, shift_samples)
-
     def init_delays(self, interpolation):
         """Initialize or design the interpolation functions for the delays
 
@@ -702,18 +675,19 @@ class Instrument:
                 self._delay_const = self._delay_dynamic
             case ("lagrange", int(order), float(min_delay_time), float(max_delay_time)):
                 # ~ self.interpolation_order = order
-                op_dyn = make_dynamic_shift_lagrange_dask(
+                op_dyn_ = make_dynamic_shift_lagrange_dask(
                     order,
                     min_delay_time * self.physics_fs,
                     max_delay_time * self.physics_fs,
                     ShiftBC.ZEROPAD,
                     ShiftBC.ZEROPAD,
                 )
-                op_fix = make_fixed_shift_lagrange_dask(
+                op_fix_ = make_fixed_shift_lagrange_dask(
                     ShiftBC.ZEROPAD, ShiftBC.ZEROPAD, order
                 )
-                self._delay_dynamic = self.numpyfy_dask_multi(op_dyn)
-                self._delay_const = self.numpyfy_dask_multi(op_fix)
+                op_dyn = self.numpyfy_dask_multi(op_dyn_)
+                op_fix = self.numpyfy_dask_multi(op_fix_)
+                self.apply_shift = AdaptiveShiftNumpy(op_fix, op_dyn, self.physics_fs)
             case (
                 "lagrange_dsp",
                 int(order),
@@ -721,18 +695,19 @@ class Instrument:
                 float(max_delay_time),
             ):
                 # ~ self.interpolation_order = order
-                op_dyn = make_dynamic_shift_dsp_dask(
+                op_dyn_ = make_dynamic_shift_dsp_dask(
                     order,
                     min_delay_time * self.physics_fs,
                     max_delay_time * self.physics_fs,
                     ShiftBC.ZEROPAD,
                     ShiftBC.ZEROPAD,
                 )
-                op_fix = make_fixed_shift_dsp_dask(
+                op_fix_ = make_fixed_shift_dsp_dask(
                     ShiftBC.ZEROPAD, ShiftBC.ZEROPAD, order
                 )
-                self._delay_dynamic = self.numpyfy_dask_multi(op_dyn)
-                self._delay_const = self.numpyfy_dask_multi(op_fix)
+                op_dyn = self.numpyfy_dask_multi(op_dyn_)
+                op_fix = self.numpyfy_dask_multi(op_fix_)
+                self.apply_shift = AdaptiveShiftNumpy(op_fix, op_dyn, self.physics_fs)
             case ("lagrange", int(order)):
                 msg = (
                     "interpolation parameter tuples need 4 entries since daskification"