docs: rework force calibration documentation

docs/theory/force_calibration/active.rst

@@ -1,3 +1,5 @@
+.. _active_calibration_theory:
+
 Active Calibration
 ------------------
 

docs/theory/force_calibration/diode.rst

@@ -0,0 +1,66 @@
+.. _diode_theory:
+
+Position sensitive detector
+---------------------------
+
+The previous section introduced the origin of the frequency spectrum of a bead in an optical trap.
+In reality, our measurement is affected by two processes:
+
+1. The motion of the bead in the trap.
+2. The response of the detector to the incident light.
+
+.. image:: figures/diode_filtering.png
+  :nbattach:
+
+This second factor depends on the type of measurement device being used.
+Typical position sensitive detectors are made of silicon.
+Such a detector has a very high bandwidth for visible light (in the MHz range).
+Unfortunately, the bandwidth is markedly reduced for the near infra-red light of the trapping laser :cite:`berg2003unintended,berg2006power`.
+This makes it less sensitive to changes in signal at high frequencies.
+
+Why is the bandwidth limited?
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The high bandwidth of visible light detection of a silicon photodiode is achieved when incoming photons are absorbed in the so-called depletion layer of the diode.
+Unfortunately, silicon has an increased transparency at the near infra-red wavelength of the trapping laser.
+The result of this is that light penetrates deeper into the substrate of the diode, where it generates charge carriers in a different region of the diode.
+These charge carriers then have to diffuse back to the depletion layer, which takes time.
+As a result, a fraction of the signal is detected with a much lower bandwidth.
+
+.. image:: figures/diode.png
+  :nbattach:
+
+This effect is often referred to as the parasitic filtering effect and is frequently modelled as a first order lowpass filter.
+This model is characterized by two numbers whose values depend on the incident laser power :cite:`berg2003unintended`:
+
+- A frequency `f_diode`, given in Hertz.
+- A unit-less relaxation factor `alpha` which reflects the fraction of light that is transmitted instantaneously.
+
+.. _high_corner_freq:
+
+High corner frequencies
+^^^^^^^^^^^^^^^^^^^^^^^
+
+In literature, the diode parameters are frequently estimated simultaneously with the calibration data :cite:`berg2003unintended,hansen2006tweezercalib,berg2006power,tolic2006calibration,tolic2004matlab,berg2004power`.
+Unfortunately, this can cause issues when calibrating at high powers.
+
+Recall that the physical spectrum is characterized by a corner frequency `fc`, and diffusion constant `D`.
+The corner frequency depends on the laser power and bead size (smaller beads resulting in higher corner frequencies).
+The parasitic filtering effect also has a corner frequency (`f_diode`) and depends on the incident intensity :cite:`berg2003unintended`.
+
+When these two frequencies get close, they cannot be estimated from the power spectrum reliably anymore.
+The reason for this is that the effects that these parameters have on the power spectrum becomes very similar.
+When working with small beads or at high laser powers, it is important to verify that the corner frequency `fc` does not approach the frequency of the filtering effect `f_diode`.
+
+Sometimes, the parameters of this diode have been characterized independently.
+In that case, the arguments `fixed_diode` and `fixed_alpha` can be passed to :func:`~lumicks.pylake.calibrate_force()` to fix these parameters to their predetermined values resolving this issue.
+For more information on how to achieve this with Pylake, please refer to the :ref:`diode calibration tutorial<diode_tutorial>`.
+
+Mathematical background
+^^^^^^^^^^^^^^^^^^^^^^^
+
+In literature, it is frequently modelled up to good accuracy with a first order approximation :cite:`berg2003unintended,tolic2006calibration,berg2006power`.
+
+.. math::
+
+    g(f, f_\mathrm{diode}, \alpha) = \alpha^2 + \frac{1 - \alpha ^ 2}{1 + (f / f_\mathrm{diode})^2} \tag{$-$}

docs/theory/force_calibration/fitting.rst

@@ -1,23 +1,38 @@
 Fitting a power spectrum
 ------------------------
 
-In the previous section, the physical origin of the power spectrum was introduced.
-However, there are some practical aspects to consider.
-So far, we have only considered the expectation value of the power spectrum.
-In reality, power spectral values follow a distribution.
+.. only:: html
 
-The real and imaginary part of the frequency spectrum are normally distributed.
+    :nbexport:`Download this page as a Jupyter notebook <self>`
+
+In the previous sections, the physical origin of the power spectrum was introduced.
+However, there are some additional practical aspects to consider.
+
+Spectral down-sampling
+^^^^^^^^^^^^^^^^^^^^^^
+
+So far, we have treated the power spectrum as a simple smooth curve where each frequency corresponds to a given amplitude.
+What we actually have been looking at is the expected value of the power spectrum.
+Basically, if you were to average many power spectra, their average would eventually converge to this curve.
+
+In reality, each data point in a power spectrum follows a distribution.
+The real and imaginary part of the complex spectrum are normally distributed.
+This spectrum is squared to obtain our power spectrum.
 As a consequence, the squared magnitude of the power spectrum is exponentially distributed.
+
 This has two consequences:
 
 - Fitting the power spectral values directly using a simple least squares fitting routine, we would
   get very biased estimates. These estimates would overestimate the plateau and corner frequency,
   resulting in overestimated trap stiffness and force response and an underestimated distance response.
 - The signal to noise ratio is poor (equal to one :cite:`norrelykke2010power`).
 
+.. _blocking_windowing:
+
 A commonly used method for dealing with this involves data averaging, which trades resolution for an
-improved signal to noise ratio. In addition, by virtue of the central limit theorem, data averaging
-leads to a more symmetric data distribution (more amenable to standard least-squares fitting procedures).
+improved signal to noise ratio.
+By virtue of the central limit theorem, as we average more data, the distribution of the data points
+becomes more and more Gaussian and therefore more amenable to standard least-squares fitting procedures.
 
 There are two ways to perform such averaging:
 
@@ -27,23 +42,133 @@ There are two ways to perform such averaging:
   spectral domain by averaging adjacent bins according to :cite:`berg2004power`. This procedure is
   referred to as *blocking*.
 
-We use the blocking method for spectral averaging, since this allows us to reject noise peaks at high
-resolution prior to averaging. Note however, that the error incurred by this blocking procedure depends
-on :math:`n_b`, the number of points per block, :math:`\Delta f`, the spectral resolution and inversely
-on the corner frequency :cite:`berg2004power`.
+.. image:: figures/blocking.gif
+  :nbattach:
+
+Blocking
+""""""""
+
+Pylake uses the blocking method for spectral averaging, since this allows us to reject noise peaks
+at high resolution prior to averaging (more on this later).
+Note however, that the error incurred by this blocking procedure depends on :math:`n_b`, the number
+of points per block, :math:`\Delta f`, the spectral resolution and inversely on the corner
+frequency :cite:`berg2004power`.
+
+.. math::
+
+    \bar{f} &= \frac{1}{n_b} \sum_{f \in block} f\\
+    \bar{P}_{meas} &= \frac{1}{n_b} \sum_{f \in block} P_{meas}(f)
+
+Setting the number of points per block too low results in a bias from insufficient averaging :cite:`berg2004power`.
+Insufficient averaging would result in an overestimation of the force response :math:`R_f` and an
+underestimation of the distance response :math:`R_d`.
+In practice, one should use a high number of points per block (:math:`n_b \gg 100`),
+unless a very low corner frequency precludes this. In such cases, it is preferable to increase the
+measurement time.
+
+Bias correction
+"""""""""""""""
+
+When sufficient blocking has taken place and noise peaks have been excluded prior to blocking,
+the spectral data points are approximately Gaussian distributed with standard deviation:
+
+.. math::
+
+    \sigma(\bar{f}) = \frac{P(\bar{f})}{\sqrt{n_b}}
+
+This means that regular weighted least squares (WLS) can take place.
+To ensure unbiased estimates in WLS, the data and squared weights must be uncorrelated.
+However, we know that there is a known correlation between these which results
+in a known bias in the estimate for the diffusion constant that can be corrected for after
+fitting :cite:`norrelykke2010power`:
+
+.. math::
+
+    D_{corrected} = D_{wls} \frac{n_b}{n_b + 1}
+
+.. _noise_floor:
+
+Noise floor
+^^^^^^^^^^^
+
+When operating at very low powers (and by extension corner frequencies), a noise floor may be visible at high frequencies.
+It is important to ensure that the upper limit of the fitting range does *not* include the noise floor as it is not taken into account in the calibration model.
+In the dataset below, we can see the effect of the noise floor::
+
+    lk.download_from_doi("10.5281/zenodo.7729823", "test_data")
+    f = lk.File("test_data/noise_floor.h5")
+
+    force_channel = f.force1x
+    reference_calibration = force_channel.calibration[0]
+    pars = {
+        "force_voltage_data": force_channel.data / reference_calibration.force_sensitivity,
+        "sample_rate": force_channel.sample_rate,
+        "bead_diameter": 4.34,
+        "temperature": 25,
+        "hydrodynamically_correct": True,
+        "num_points_per_block": 150,
+        "fit_range": [100, 20000],
+    }
+
+    plt.figure(figsize=(7, 5))
+    plt.subplot(2, 2, 1)
+    calibration = lk.calibrate_force(**pars)
+    calibration.plot()
+    plt.title(f"Stiffness: {calibration.stiffness:.2f}")
+    plt.subplot(2, 2, 3)
+    calibration.plot_spectrum_residual()
+
+    plt.subplot(2, 2, 2)
+    calibration = lk.calibrate_force(**pars | {"fit_range": [100, 3000]})
+    calibration.plot()
+    plt.title(f"Stiffness: {calibration.stiffness:.2f}")
+    plt.subplot(2, 2, 4)
+    calibration.plot_spectrum_residual()
+    plt.tight_layout()
+
+.. image:: figures/noise_floor_free_diode.png
+
+Note that when we have a diode calibration, excluding the noise floor becomes even more important::
+
+    diode_calibration = reference_calibration.diode_calibration
+    pars = pars | diode_calibration(f["Diagnostics"]["Trap power 1"])
+
+    plt.figure(figsize=(7, 5))
+    plt.subplot(2, 2, 1)
+    calibration = lk.calibrate_force(**pars)
+    calibration.plot()
+    plt.title(f"Stiffness: {calibration.stiffness:.2f}")
+    plt.subplot(2, 2, 3)
+    calibration.plot_spectrum_residual()
+
+    plt.subplot(2, 2, 2)
+    calibration = lk.calibrate_force(**pars | {"fit_range": [100, 3000]})
+    calibration.plot()
+    plt.title(f"Stiffness: {calibration.stiffness:.2f}")
+    plt.subplot(2, 2, 4)
+    calibration.plot_spectrum_residual()
+    plt.tight_layout()
+
+.. image:: figures/noise_floor_fixed_diode.png
+
+The reason the effect of the noise floor on the calibration parameters is more pronounced is because with the fixed diode model, the model is less free to adjust the diode model to mitigate its impact.
+As a result, the model uses the corner frequency in an attempt to capture the shape of the noise floor (strongly biasing the result).
+
+Noise peaks
+^^^^^^^^^^^
 
-Setting the number of points per block too low would result in a bias from insufficient averaging
-:cite:`berg2004power`. Insufficient averaging would result in an overestimation of the force response
-(:math:`R_f`) and an underestimation of the distance response (:math:`R_d`). In practice, one should
-use a high number of points per block (:math:`n_b \gg 100`), unless a very low corner frequency precludes this.
-In such cases, it is preferable to increase the measurement time.
+Optical tweezers are precision instruments.
+Despite careful determination and elimination of noise sources, it is not always possible to exclude all potential sources of noise.
+One downside of weighted least squares estimation, is that it is very sensitive to outliers.
+It is therefore important to either exclude noise peaks from the data prior to fitting or use :ref:`robust fitting<robust_fitting>`.
+Noise peaks are always excluded prior to blocking to minimize data loss.
 
 .. _goodness_of_fit:
 
 Goodness of fit
 ---------------
 
-Based on the model and noise assumptions, we can calculate a goodness of fit criterion.
+When working with the Gaussian error model, we can calculate a goodness of fit criterion.
 When sufficient blocking has taken place, the sum of squared residuals that is being minimized during the fitting procedure is distributed according to a chi-squared distribution characterized by :math:`N_{\mathit{dof}} = N_{\mathit{data}} - N_{\mathit{free}}` degrees of freedom.
 Here :math:`N_{\mathit{data}}` corresponds to the number of data points we fitted (after blocking) and :math:`N_{\mathit{free}}` corresponds to the number of parameters we fitted.
 We can use the value we obtain to determine how unusual the fit error we obtained is.

docs/theory/force_calibration/force_calibration.rst

@@ -8,9 +8,8 @@ Introduction
 Why is force calibration necessary?
 -----------------------------------
 
-Optical tweezers typically measure forces and displacements by detecting deflections of a trapping
-laser by a trapped bead. These deflections are measured at the back focal plane of the beam using
-position sensitive detectors (PSDs).
+Optical tweezers typically measure forces and displacements by detecting deflections of a trapping laser by a trapped bead.
+These deflections are measured at the back focal plane of the beam using position sensitive detectors (PSDs).
 
 .. image:: figures/back_focal.png
   :nbattach:
@@ -30,15 +29,57 @@ and
 
 Where :math:`V` is the position-dependent voltage signal from the PSD and :math:`R_d` and :math:`R_f`
 are the displacement and force sensitivity proportionality constants, respectively.
-Force calibration refers to computing these conversion factors.
+
+Force calibration refers to computing the calibration factors needed to convert from raw voltages to actual forces and displacements.
+The values we wish to calculate are:
+
+- Trap stiffness :math:`\kappa`, which reflects how strongly a bead is held by a trap.
+- Force response :math:`R_d`, the proportionality constant between voltage and force.
+- Distance response :math:`R_f`, the proportionality constant between voltage and distance.
 
 Several methods exist to calibrate optical traps based on sensor signals.
 In this section, we will provide an overview of the physical background of power spectral calibration.
 
-Why does the power spectrum look the way it does?
--------------------------------------------------
+How can we calibrate?
+---------------------
+
+We will start by building some intuition about the physical processes.
+Consider a small bead suspended in fluid (no optical trapping taking place).
+This bead moves around due to the random collisions of molecules against the bead.
+If we observe this motion over time, then we can see that it is a lot like the bead is taking aimless steps.
+This unimpeded movement is often called free diffusion or Brownian motion.
+When we look at the motion of the bead over longer time periods, then each random step contributes only a little to the overall movement.
+Because these steps are random, they don't perfectly cancel out in the short term, and lead to a gradual drift.
+If there is no optical trap keeping the bead in place, then the bead slowly drifts off from its starting position.
+
+One way to analyze this motion is to make a power spectrum plot of the bead's position.
+This plot shows how much motion there is at different frequencies of movement.
+Lower frequencies correspond to longer time intervals.
+These were the time intervals we associated with the broader, slow movements of the bead.
+If we think about a bead moving due to random collisions, then we can expect that the bead will move more in these longer time intervals.
+This is why in the power spectrum of free diffusion, we see a lot more energy concentrated at these low frequencies, while the rapid jiggles contribute far less.
+
+.. image:: figures/sim_trap_opt.gif
+  :nbattach:
+
+Now we introduce the optical trap.
+The trap pulls the bead back to the laser focus.
+The strength of this pull depends on how far the bead is from the focus (like a spring).
+If we consider again the spectrum of free diffusion, then we can see why the trap will limit the motion at low frequencies more.
+It is because those motions are larger and will experience a stronger pull of the trap.
+That means they will be damped, because the furthest extents are now limited by the trap keeping the bead in place.
+Since the trap limits the maximum amplitude of motion, this manifests itself as sharp reduction of amplitudes above a certain threshold.
+This is why we typically see a plateau at low frequencies in the power spectrum of a trapped bead.
+
+Important takeaways
+-------------------
+
+- The spectrum of bead motion in a trap can be characterized by a diffusion constant and corner frequency.
+- At low frequencies the trapping force dominates, limiting the amplitudes, while at high frequencies the drag on the bead does.
+
+Mathematical background
+-----------------------
 
-Consider a small bead freely diffusing in a medium (no optical trapping taking place).
 Neglecting hydrodynamic and inertial effects (more on this later), we obtain the following equation of motion:
 
 .. math::