calibration: add private api for motion blur model

lumicks/pylake/force_calibration/calibration_models.py

@@ -9,6 +9,8 @@
 from .detail.power_models import (
     g_diode,
     alias_spectrum,
+    motion_blur_peak,
+    motion_blur_spectrum,
     sphere_friction_coefficient,
     passive_power_spectrum_model,
     theoretical_driving_power_lorentzian,
@@ -414,8 +416,23 @@ def _calculate_power_spectral_density(self, f, fc, diffusion_constant, *filter_p
     def __call__(self, f, fc, diffusion_constant, *filter_params) -> np.ndarray:
         return self._calculate_power_spectral_density(f, fc, diffusion_constant, *filter_params)
 
+    def _motion_blur(self, acquisition_time):
+        """Include effects of motion blur into the model
+
+        Parameters
+        ----------
+        acquisition_time : float
+            Acquisition time in seconds
+        """
+        new_model = copy(self)
+        new_model._calculate_power_spectral_density = motion_blur_spectrum(
+            new_model._calculate_power_spectral_density, acquisition_time
+        )
+
+        return new_model
+
     def _alias_model(self, sample_rate, num_aliases):
-        """Include effects of aliasing into the model
+        """Include effects of aliasing into the model.
 
         Parameters
         ----------
@@ -768,6 +785,24 @@ def _theoretical_driving_power(self, f_corner):
         integrated over the frequency bin corresponding to the driving input."""
         return self._theoretical_driving_power_model(f_corner)
 
+    def _motion_blur(self, acquisition_time):
+        """Include effects of motion blur into the model
+
+        Parameters
+        ----------
+        acquisition_time : float
+            Acquisition time in seconds
+        """
+        new_model = copy(self)
+        new_model._calculate_power_spectral_density = motion_blur_spectrum(
+            self._calculate_power_spectral_density, acquisition_time
+        )
+        new_model._theoretical_driving_power_model = motion_blur_peak(
+            self._theoretical_driving_power_model, self.driving_frequency, acquisition_time
+        )
+
+        return new_model
+
     def calibration_results(
         self,
         fc,

lumicks/pylake/force_calibration/detail/power_models.py

@@ -159,6 +159,8 @@ def theoretical_driving_power_lorentzian(fc, driving_frequency, driving_amplitud
     power spectrum. It corresponds to the driven power spectrum minus the thermal power spectrum
     integrated over the frequency bin corresponding to the driving input.
 
+    Parameters
+    ----------
     fc : float
         Corner frequency [Hz]
     driving_frequency : float
@@ -169,6 +171,54 @@ def theoretical_driving_power_lorentzian(fc, driving_frequency, driving_amplitud
     return driving_amplitude**2 / (2 * (1 + (fc / driving_frequency) ** 2))
 
 
+def motion_blur_peak(peak, driving_frequency, acquisition_time):
+    """Take into account motion blur on the driving peak.
+
+    Parameters
+    ----------
+    peak : callable
+        Function which takes a corner frequency and produces an estimated peak power.
+    driving_frequency : float
+        Driving frequency.
+    acquisition_time : float
+        Acquisition time in seconds.
+
+    References
+    ----------
+    .. [1] Wong, W. P., & Halvorsen, K. (2006). The effect of integration time on fluctuation
+       measurements: calibrating an optical trap in the presence of motion blur. Optics express,
+       14(25), 12517-12531.
+    """
+
+    def blurred(fc, *args, **kwargs):
+        return peak(fc, *args, **kwargs) * np.sinc(driving_frequency * acquisition_time) ** 2
+
+    return blurred
+
+
+def motion_blur_spectrum(psd, acquisition_time):
+    """Take into account motion blur on the spectrum.
+
+    Parameters
+    ----------
+    psd : callable
+        Function which takes a numpy array of frequencies and returns a power spectral density.
+    acquisition_time : float
+        Acquisition time in seconds.
+
+    References
+    ----------
+    .. [1] Wong, W. P., & Halvorsen, K. (2006). The effect of integration time on fluctuation
+       measurements: calibrating an optical trap in the presence of motion blur. Optics express,
+       14(25), 12517-12531.
+    """
+
+    def blurred(freq, *args, **kwargs):
+        return psd(freq, *args, **kwargs) * np.sinc(freq * acquisition_time) ** 2
+
+    return blurred
+
+
 def alias_spectrum(psd, sample_rate, num_alias=10):
     """
     Produce an aliased version of the input PSD function.

lumicks/pylake/force_calibration/tests/data/simulate_ideal.py

@@ -1,6 +1,8 @@
 import numpy as np
 import scipy
 
+from lumicks.pylake.detail.utilities import downsample
+
 
 def calculate_active_term(time, driving_sinusoid, stiffness, gamma):
     """Simulate stage and bead position from simulation parameters.
@@ -124,8 +126,11 @@ def simulate_calibration_data(
         Temperature [C]
     pos_response_um_volt : float
         Response [um/V], also denoted in papers as Rd
-    anti_aliasing : bool
+    anti_aliasing : str, optional
         Should we anti-alias the data?
+        "fir": anti-alias with FIR filter
+        "integrate": just downsample by integration
+        None: just decimate by subsampling
     oversampling : int
         Oversampling factor. Relatively high oversampling ratios are required to reject aliasing
         effectively.
@@ -175,11 +180,15 @@ def simulate_calibration_data(
         positions = apply_diode_filter(positions, *diode, oversampled_dt)
 
     # Decimate the signal
-    positions = (
-        scipy.signal.decimate(positions, oversampling, ftype="fir")
-        if anti_aliasing
-        else positions[::oversampling]
-    )
+    if anti_aliasing == "fir":
+        positions = scipy.signal.decimate(positions, oversampling, ftype="fir")
+    elif anti_aliasing == "integrate":
+        positions = downsample(positions, oversampling, reduce=np.mean)
+    else:
+        if anti_aliasing:
+            raise RuntimeError("Chose invalid anti aliasing method")
+        positions = positions[::oversampling]
+
     time = time[::oversampling]
 
     if driving_sinusoid:

lumicks/pylake/force_calibration/tests/test_camera_calibration.py

@@ -0,0 +1,66 @@
+import numpy as np
+import pytest
+
+from lumicks.pylake.force_calibration.calibration_models import (
+    ActiveCalibrationModel,
+    PassiveCalibrationModel,
+)
+from lumicks.pylake.force_calibration.power_spectrum_calibration import (
+    fit_power_spectrum,
+    calculate_power_spectrum,
+)
+
+from .data.simulate_ideal import simulate_calibration_data
+
+
+# @pytest.mark.slow
+def test_camera_calibration():
+    np.random.seed(909)
+    f_drive = 16.8
+    sample_rate = 500
+    shared_params = {
+        "bead_diameter": 1.01,
+        "viscosity": 1.002e-3,
+        "temperature": 20,
+    }
+    params = {
+        **shared_params,
+        "duration": 600,
+        "sample_rate": sample_rate,
+        "stiffness": 0.05,
+        "pos_response_um_volt": 0.618,
+        "driving_sinusoid": (500, f_drive),
+        "diode": None,
+    }
+
+    oversampling = 128
+    pos, nano = simulate_calibration_data(
+        **params, anti_aliasing="integrate", oversampling=oversampling
+    )
+    passive_model = PassiveCalibrationModel(**shared_params, fast_sensor=True)
+
+    active_model = ActiveCalibrationModel(
+        force_voltage_data=pos,
+        driving_data=nano,
+        **shared_params,
+        sample_rate=sample_rate,
+        driving_frequency_guess=f_drive,
+        fast_sensor=True,
+    )
+
+    # Consider aliasing and motion blur
+    models = (
+        model._motion_blur(1.0 / sample_rate)._alias_model(sample_rate, 20)
+        for model in (passive_model, active_model)
+    )
+
+    distance_responses = [0.6234771740528027, 0.623431094357452]
+    stiffnesses = [0.04950015364112712, 0.04950747132648579]
+    force_responses = [30.86221590734949, 30.864497027941212]
+
+    ps = calculate_power_spectrum(pos, sample_rate, fit_range=(25, 1000), num_points_per_block=200)
+    for model, rd, kappa, rf in zip(models, distance_responses, stiffnesses, force_responses):
+        calibration = fit_power_spectrum(ps, model)
+        np.testing.assert_allclose(calibration["Rd"].value, rd)
+        np.testing.assert_allclose(calibration["Rf"].value, rf)
+        np.testing.assert_allclose(calibration["kappa"].value, kappa)