diff --git a/bilby/gw/detector/calibration.py b/bilby/gw/detector/calibration.py
index cf4fc7d9392e35f369edb31d2379b5430b65a058..5b5d4b1ec8539174a9fdb39f8b69e7e5d4d2dc60 100644
--- a/bilby/gw/detector/calibration.py
+++ b/bilby/gw/detector/calibration.py
@@ -94,6 +94,34 @@ class CubicSpline(Recalibrate):
         return self.__class__.__name__ + '(prefix=\'{}\', minimum_frequency={}, maximum_frequency={}, n_points={})'\
             .format(self.prefix, self.minimum_frequency, self.maximum_frequency, self.n_points)
 
+    @property
+    def delta_amplitude_interp(self):
+        """ Return the updated interpolant or generate a new one if required """
+        amplitude_parameters = [self.params['amplitude_{}'.format(ii)]
+                                for ii in range(self.n_points)]
+
+        if hasattr(self, '_delta_amplitude_interp'):
+            self._delta_amplitude_interp.y = amplitude_parameters
+        else:
+            self._delta_amplitude_interp = interp1d(
+                self.log_spline_points, amplitude_parameters, kind='cubic',
+                bounds_error=False, fill_value=0)
+        return self._delta_amplitude_interp
+
+    @property
+    def delta_phase_interp(self):
+        """ Return the updated interpolant or generate a new one if required """
+        phase_parameters = [
+            self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
+
+        if hasattr(self, '_delta_phase_interp'):
+            self._delta_phase_interp.y = phase_parameters
+        else:
+            self._delta_phase_interp = interp1d(
+                self.log_spline_points, phase_parameters, kind='cubic',
+                bounds_error=False, fill_value=0)
+        return self._delta_phase_interp
+
     def get_calibration_factor(self, frequency_array, **params):
         """Apply calibration model
 
@@ -113,17 +141,10 @@ class CubicSpline(Recalibrate):
             The factor to multiply the strain by.
         """
         self.set_calibration_parameters(**params)
-        amplitude_parameters = [self.params['amplitude_{}'.format(ii)]
-                                for ii in range(self.n_points)]
-        delta_amplitude = interp1d(
-            self.log_spline_points, amplitude_parameters, kind='cubic',
-            bounds_error=False, fill_value=0)(np.log10(frequency_array))
+        log10_frequency_array = np.log10(frequency_array)
 
-        phase_parameters = [
-            self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
-        delta_phase = interp1d(
-            self.log_spline_points, phase_parameters, kind='cubic',
-            bounds_error=False, fill_value=0)(np.log10(frequency_array))
+        delta_amplitude = self.delta_amplitude_interp(log10_frequency_array)
+        delta_phase = self.delta_phase_interp(log10_frequency_array)
 
         calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)
 
diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index ab78a15b634d5de9403ea3c5132016b098717419..0e1d0df0ca39ae93ed7dc47e7ac5b2a89173e81f 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -865,6 +865,18 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
         self.frequency_nodes_quadratic = \
             waveform_generator.waveform_arguments['frequency_nodes_quadratic']
 
+    def get_d_inner_h_interp(self, d_inner_h_tc_array, indices):
+        """ Return the updated interpolant or generate a new one if required """
+        if hasattr(self, '_d_inner_h_interp'):
+            self._d_inner_h_interp.x = self.weights['time_samples'][indices]
+            self._d_inner_h_interp.y = d_inner_h_tc_array
+        else:
+            self._d_inner_h_interp = interp1d(
+                self.weights['time_samples'][indices],
+                d_inner_h_tc_array, kind='cubic', assume_sorted=True)
+
+        return self._d_inner_h_interp
+
     def calculate_snrs(self, waveform_polarizations, interferometer):
         """
         Compute the snrs for ROQ
@@ -916,9 +928,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             'i,ji->j', np.conjugate(h_plus_linear + h_cross_linear),
             self.weights[interferometer.name + '_linear'][indices])
 
-        d_inner_h = interp1d(
-            self.weights['time_samples'][indices],
-            d_inner_h_tc_array, kind='cubic', assume_sorted=True)(ifo_time)
+        d_inner_h = self.get_d_inner_h_interp(d_inner_h_tc_array, indices)(ifo_time)
 
         optimal_snr_squared = \
             np.vdot(np.abs(h_plus_quadratic + h_cross_quadratic)**2,