Compensating for phase drifts in holographic measurements

Digital holography allows measuring the complex amplitude of a given wavefront. We presented in detail the off-axis holography approach. However, it requires a separate reference arm. Due to air flow, vibrations, or other perturbations, the optical path length difference between the two arms can fluctuate in time, even in controlled lab experiments on a good optical table. This means that the phase of the measured wavefront is estimated up to a global phase that can randomly change over time. This is very detrimental for transmission matrix measurements as the relative phase between each column has to be precisely estimated. This is particularly true when the measurement time can take few minutes or more when using a liquid crystal spatial light modulator that has a limited frame rate. In [R. Mouthaan et al., Appl. Opt. (2022)], the authors propose a simple yet robust way to compensate for phase fluctuations, even when the phase changes completely between two frames.

Introduction

When using reference-based holographic techniques (typically phase-shifting or off-axis holography) to estimate the phase, even with a rather stable lab setup, it is common to observe a phase drift between the two arms of the system. If not compensated for, it introduces errors in the relative global phase between successive measurements. For many applications, it is not an issue, but when measuring the transmission matrix of a given medium, this introduces random relative phases between pairs of columns of the matrix. It is very detrimental when trying, for instance, to use the matrix to reconstruct an image.

One solution consists in sending, at a regular interval, a known pattern. If the phase between the signal and the reference arm has shifted, we will measure the same complex pattern multiplied by a global phase. We can then estimate this phase drift and correct the other measurements. This is ok as long as the phase fluctuations are slow compared to the frame rate of your measurement system, i.e. if the phase can be considered stable between two consecutive measurements. If one has a slow setup or a quite unstable one, that can be an issue.

Principle

The solution suggested by Mouthaan and his collaborators is the following. First, we chose an input reference field \(\psi_{ref}^{in}\), and measure its complex output field \(\psi_{ref}^{out}\). Then, we want to send a set of input wavefronts and measure the resulting output complex fields (to estimate the transmission matrix). For a given input wavefront \(\psi_{0}^{in}\), the output wavefront is given by the transmission matrix of the system that we want to measure : \(\psi_{0}^{out} = \mathbf{T}\psi_{0}^{in}\). Because of the potential phase drift, we only have access to

\begin{equation}
\hat{\psi}_{0}^{out} = e^{j\phi_0}\psi_{0}^{out} \, ,
\end{equation}

where \(\phi_0\) is an unkonwn phase. The idea is, for each input \(\psi_{0}^{in}\), to also send the sum of this wavefront and the reference one, i.e. \(\psi_{ref}^{in}+\psi_{0}^{in}\). The resulting measured output wavefront is given by:

\begin{equation}
\hat{\psi}_{s}^{out} = e^{j\phi_s}\left(\psi_{ref}^{out}+\psi_{0}^{out}\right) \,,
\end{equation}

with \(\phi_0\) is another unkonwn phase.

To estimate the phase drift  \(\phi_0\), we compute the absolute value of the complex correlation between \(\hat{\psi}_{s}^{out}\) and \(\hat{\psi}_{0}^{out} e^{j\alpha} + \hat{\psi}_{ref}^{out}\), where \(\alpha\) is a free parameter:

\begin{equation}
C(\alpha) =
\frac{
\left|
\left\langle \hat{\psi}_{0}^{out} e^{j\alpha} + \hat{\psi}_{ref}^{out}| \hat{\psi}_{s}^{out} \right\rangle
\right|
}{
\left\langle \hat{\psi}_{0}^{out} e^{j\alpha} + \psi_{ref}^{out}| \hat{\psi}_{0}^{out} e^{j\alpha} + \hat{\psi}_{ref}^{out} \right\rangle
\left\langle \hat{\psi}_{s}^{out}| \hat{\psi}_{s}^{out} \right\rangle
}\, ,
\end{equation}

with \(\left\langle A| B \right\rangle = \sum_k A_k B_k^*\) the scalar product between two vectors.

Let's look at the upper part of this equation, we can rewrite it:

\begin{equation}
\left|\left\langle \hat{\psi}_{0}^{out} e^{j\alpha} + \psi_{ref}^{out}| \hat{\psi}_{s}^{out} \right\rangle \right|
=
\left|\left\langle \hat{\psi}_{0}^{out} e^{j\alpha} + \psi_{ref}^{out}| e^{j\phi_s}\left(\psi_{ref}^{out}+\psi_{0}^{out}\right)\right\rangle \right|
=
\left| \left\langle \psi_{0}^{out} e^{j\left(\alpha+\phi_0\right)} + \psi_{ref}^{out}| \psi_{ref}^{out}+\psi_{0}^{out}\right\rangle \right|
\end{equation}

We see that the correlation is maximal, i.e. \(C(\alpha) = 1\), when \(\alpha = -\phi_0\). We then need to numerically find the value \(\alpha_{max}\) that maximizes \(C(\alpha)\). We can then retrieve the corrected ouptut field \(\psi_{corr}^{out}\) using:

\begin{equation}
\psi_{corr}^{out} = \hat{\psi}_{0}^{out} e^{j\alpha_{max}} \,.
\end{equation}

Conclusion

The proposed technique allows compensating for phase drifts, even fast ones. Indeed, while the phase still need to be stable during one image acquisition, it does not need to be stable between two image acquisitions, as it is the case with other techniques. The drawback is that we need to double the number of measurements, as wee need to measure the response to the wavefront plus the reference in addition to just the wavefront.

Python Implementation

To find \(\alpha_{max}\), the easiest solution is to test different values of \(\alpha\) between \(0\) and \(2\pi\) and keep the best value. One needs to test between 10 and 20 values to have a reasonable compensation of the phase drift, which can slow down the process when we have a lot of images to process. I prefer using the minimize() function from the scipy package. It takes about 6-7 iterations to converge.

Here is a function to retrieve the corrected wavefront from the measurement E (\(\hat{\psi}_{0}^{out}\)), Eref (\(\psi_{ref}^{out}\)) and Esum (\(\hat{\psi}_{sum}^{out}\)).

from scipy.optimize import minimize
from numpy import corrcoef

complex_corr = lambda A,B : corrcoef(A,B)[0,1]

def get_corr_phase(E, Eref, Esum):
    E0 = E
    E1 = (Eref+np.exp(1j*phi)*E).ravel().transpose()
    E2 = Esum.ravel()
    cost = lambda phi: -np.abs(complex_corr(E1,E2))
    res = minimize(cost, x0 = 0., method='Powell', options = {'disp' : False})
    best_phase = res['x'][0]
    return np.exp(1j*best_phase)*E0