1. INTRODUCTION

Light traveling through a heterogeneous medium is repeatedly scattered and its wavefront is distorted to form a speckle pattern. Similarly, light coupled into a multimode fiber (MMF) propagates in discrete modes that couple to each other and build up phase delays with respect to each other, forming a speckle pattern at the output side. Both these processes are mathematically linear, and so the relationship between the proximal fields and the distal fields can be described by a transmission matrix [1]. The transmission matrix $\textbf{T}$ relates a vector ${E_{\rm{in}}}$ describing the complex fields in the input plane to a vector ${E_{\rm{out}}}$ describing the complex fields in the output plane:

Relating the output wavefront to the input wavefront allows the distortion imposed by the scattering medium to be corrected. This approach has opened up applications such as imaging through a scattering medium [2], and endoscopy through hair-thin fiber endoscopes [3,4]. The transmission matrix can also be used for projection rather than imaging, allowing a diffraction-limited spot [5] or a more complex image [6] to be formed through a scattering medium. By combining imaging and projection, optical tools such as a confocal microscope can even be implemented through an optical fiber [7]. Transmission matrices can also be used to identify eigenchannels that allow maximum energy to be transmitted through the scattering medium [8] or to find wavefronts that optimally probe specific measurement parameters [9].

Measurement noise is inevitable in any experiment but is only rarely considered in the context of transmission matrix measurements. Wang et al. [10] consider static Gaussian noise superimposed on the output fields, corresponding to time-invariant imperfections in the measurement system. Mastiani and Vellekoop [11] propose a measurement technique that maximizes the amount of signal detected, ensuring the signal to noise ratio (SNR) is maximized.

This paper presents several insights that allow the transmission matrix and inverse transmission matrix of a scattering medium or optical waveguide to be robustly characterized in the presence of time-varying measurement noise. In Section 2, a typical interferometric measurement system for obtaining the transmission matrix of an optical fiber is presented. Next, sources of measurement noise in such a setup are discussed in Section 3 and quantified in Section 4. Interferometric measurements are then used to construct the transmission matrix, as simulated in Section 5. It is shown that an appropriate choice of measurement basis allows the effects of measurement noise on the matrix construction to be mitigated. Finally, in Section 6, it is demonstrated that knowledge of the measurement noise floor of the experimental system allows an optimal amount of regularization to be applied when calculating the inverse transmission matrix.

2. EXPERIMENTAL SETUP

The experimental setup considered is very similar to that presented in Ref. [12], except that a single lens is used to couple light into the MMF. The modified experimental system is illustrated in Fig. 1. In brief, a beam is modulated using a phase-modulating spatial light modulator (SLM) before being projected onto the proximal facet of an MMF. An interferometer at the output side is used to map the complex fields exiting the MMF. The SLM used is a Jasper JD8714 nematic liquid crystal on silicon device with a resolution of ${4096} \times {2400}$ and 3.74 µm pixels capable of 256 discrete levels of phase modulation. The scattering medium used is a 2 m optical patch cord (Thorlabs M42L02) with a 50 µm diameter core and a numerical aperture of 0.22 supporting approximately 745 modes in a single polarization. A Thorlabs DC3240M camera was used, triggered off the SLM clock cycle. It is proposed that this experimental setup is broadly representative of most off-axis transmission matrix characterization setups.

 

Fig. 1. Schematic of experimental setup. BS, beam splitter; PBS, polarizing beam splitter; SLM, spatial light modulator; SMF, single mode fiber; MMF, multimode fiber; CAM, camera; POL, polarizer; $O$, object beam; $R$, reference beam; $T$, transmission matrix.

Download Full Size | PDF

The fiber input facet is positioned at the focal plane of the proximal 4 mm lens so that the light pattern projected onto the input fiber facet is the Fourier transform of the light field in the SLM plane (to within a multiplicative phase term). A phase ramp is applied to the displayed hologram and the SLM is tilted, such that the modulated first order is incident on the fiber core while the unmodulated zeroth order misses the fiber core. The SLM pixels are grouped into macropixels. Macropixels can be turned “on” by displaying a phase grating that directs light towards the fiber facet, and “off” by displaying a checkerboard pattern that directs light away from the fiber facet.

The complex field distribution of the light exiting the MMF is inferred from the captured camera image using off-axis holography. Off-axis holography is a widely established technique, and excellent reviews can be found in Refs. [13,14]. Phase drift between consecutive interferometric measurements is eliminated using the technique described in Ref. [12].

The transmission matrix is defined as relating SLM macropixels to the complex fields captured by the interferometer. This approach captures various effects, including lens and SLM aberrations. It is approximately equivalent to the transmission matrix relating the SLM plane to the distal fiber facet, which is direclty imaged onto the interferometer camera. It is, however, not equivalent to the transmission matrix relating the proximal fiber facet to the distal fiber facet, as it is the Fourier transform of the SLM that is projected onto the proximal fiber facet.

In the SLM input plane, a region of ${1250} \times {1250}\;{\rm pixels}$ is grouped into ${32} \times {32}$ macropixels. In the output plane, a ${448} \times {448}\;{\rm pixel}$ region of the camera sensor is captured. As part of the off-axis holography calculations, the image is down-sampled in the Fourier domain, giving a smaller ${60} \times {60}$ representation. This yields the ${3600} \times {1024}$ transmission matrix considered throughout this paper.

 

Fig. 2. Standard deviation of (a) magnitude, (b) real component, and (c) imaginary component of object beam noise. Measured in the experimental system by displaying a static hologram and capturing a sequence of 250 complex fields using the off-axis interferometer.

Download Full Size | PDF

3. DOMINANT SOURCES OF EXPERIMENTAL MEASUREMENT NOISE

Sources of measurement noise in interferometers have previously been discussed [15–19], and a brief overview of the anticipated dominant sources is listed here.

Some of the noise sources in this experiment are specific to camera measurements, and others are specific to interferometers. In some cases, the noise contribution is expected to be small from the outset (camera readout noise, camera dark noise, and camera quantization noise). In other cases (optical fiber stability, temperature stability, vibrations, laser wavelength stability, and air flow), it has been possible to reduce the contribution by improving the mechanical stability and shielding of the setup. Unfortunately, there are remaining noise sources such as laser power fluctuations and camera shot noise that cannot be removed and are therefore expected to dominate. It is noted that these dominant sources of measurement noise are expected to be time varying, which will be used to quantify their contributions.

4. EFFECT OF MEASUREMENT NOISE ON AN OFF-AXIS INTERFEROMETRY MEASUREMENT

Several authors have considered the effects of noise on off-axis interferometry measurements. For example, Charrière et al. consider quantization noise [20], finding that 6-bit pixel depths are sufficient for phase imaging. Pandey and Hennelley [21] consider quantization noise in greater depth, deriving equations for the consequent uncertainties on the object beam’s real component, imaginary component, magnitude, and phase. Charrière et al. [20] highlight that an off-axis interferometer is a form of heterodyne detector, and as such, spatially uncorrelated noise can be filtered out and the SNR significantly enhanced. This can be exploited to reach the fundamental shot-noise limit, as was later experimentally demonstrated by others [22–25].

Here, a sequence of measurements has been taken to characterize the noise of the experimental system. First, a static diffraction-limited spot was focused onto the fiber input facet to create a speckle pattern at the fiber output, and the off-axis interferometer was used to capture 250 time-sequential complex field distributions. Next, the time-averaged object beam amplitude is calculated for each pixel, allowing the standard deviation ${\sigma _{{\rm err}}}$ of the magnitude [Fig. 2(a)], the real component [Fig. 2(b)], and the imaginary component [Fig. 2(c)] of the object beam $O$ to be determined as a function of the object beam amplitude.

The overall trend of these data follows a square root trend, indicating that shot noise is dominant. We can infer that other sources of noise have been sufficiently attenuated so as to be negligible. As previously noted, Charrière et al. [20] posit that the off-axis holography calculation filters out spatially uncorrelated noise. While noise introduced by the camera is mostly spatially uncorrelated and will be reduced by the off-axis interferometry calculation, other noise sources such as that introduced by the SLM, air flow, mechanical stability, or wavelength drift do introduce spatially correlated noise and can therefore not be attenuated with this approach. Furthermore, a spatial correlation may also be introduced by optical memory effects in the fiber. A least-squares trendline has been fit to both the real component noise and the imaginary component noise, and a relationship of ${\sigma _{\rm{err}}} = 0.21\sqrt {|O|}$ is obtained in both cases.

 

Fig. 3. (a)–(c) Different sequences of patterns are used to acquire the transmission matrix, with associated colormap shown in (d). (e) shows that measurement noise must be considered when deciding which sequence to use. Simulations have been performed with a constant power incident on the SLM. (a) Macropixel patterns. (b) Random patterns. (c) Hadamard patterns. (d) Color map. (e) Transmission matrix reconstruction error.

Download Full Size | PDF

5. EFFECT OF MEASUREMENT NOISE ON CONSTRUCTION OF THE FORWARD TRANSMISSION MATRIX

The transmission matrix is characterized by displaying a sequence of $o$ input patterns on the SLM and measuring the corresponding output fields. The transmission matrix is then constructed by concatenating the obtained input patterns to form a matrix ${\textbf{E}_{\textbf{in}}}$, concatenating the obtained output fields to form a matrix ${\textbf{E}_{\textbf{out}}}$, and performing the following calculation:

Similarly, the inverse transmission matrix can be constructed using

Here, ${^\dagger}$ denotes the Moore–Penrose pseudoinverse, a generalized inverse that yields a least-squares solution to the system of equations, and can be calculated for non-square matrices. No regularization is applied in this section when calculating the Moore–Penrose pseudoinverse. Instead, optimal regularization when performing this calculation will be considered in the next section.

A variety of pattern sequences can be displayed on the SLM during the data acquisition process, including macropixel patterns, where a single macropixel is translated across the SLM [Fig. 3(a)]; random patterns, where randomly generated phase profiles are displayed on the SLM [Fig. 3(b)]; and Hadamard patterns, which are constructed from the rows or columns of a Hadamard matrix [Fig. 3(c)].

A simple computational model was developed to evaluate the effects of measurement noise on the transmission matrix calculation when different pattern sequences are used. An $m \times n = 3600 \times 1024$ complex transmission matrix is randomly generated and the singular values are modified such that these are one for the first 745 modes, with the remaining singular values following an inverse decay. The input pattern ${\textbf{E}_{\textbf{in}}}$ and the transmission matrix $\textbf{T}$ are used to generate the output fields ${\textbf{E}_{\textbf{out}}} = \textbf{T}{\textbf{E}_{\textbf{in}}}$. Noise following the trends of Fig. 2 is superimposed on the output fields, and Eq. (2) is used to calculate the estimated transmission matrix $\tilde{\textbf{T}}$. This is then compared to the ground truth transmission matrix $\textbf{T}$ using the following phase-sensitive normalized mean squared error (NMSE) error metric, where sums are performed over all matrix elements:

These simulations have been performed for varying numbers of measurements $o$ with results shown in Fig. 3(e).

Additionally, each reconstructed transmission matrix is used to model an imaging process. This is done by calculating the speckle pattern ${E_{\rm speckle}} = \textbf{T}{E_{\rm image}}$ formed at the distal end of the fiber, and then using ${\tilde E_{\rm image}} = {\textbf{T}^\dagger}{E_{\rm speckle}}$ to obtain an estimate of the displayed image. In each case, this is done using the transmission matrix obtained using 1024 input and output patterns. The obtained images are shown in Figs. 4(a)–4(d).

 

Fig. 4. (a) Test image is reconstructed using a transmission matrix unaffected by noise. (c), (d) The same image is reconstructed using transmission matrices affected by noise. Results are shown for transmission matrices constructed using (b) macropixel input patterns, (c) random input patterns, and (d) Hadamard input patterns.

Download Full Size | PDF

If no noise is superimposed on the output fields, then all pattern sequences give the same NMSE trend, decreasing linearly for $o \lt n$ and perfectly predicting the transmission matrix for $o \gt n$. It is noted that if the number of camera pixels had been fewer than the number of SLM pixels ($m \lt n$), this cutoff would have occurred at $o = m$. Figure 4(a) shows that the reconstructed image is perfectly formed.

It can be seen that the use of the macropixel patterns as input fields gives an error of more than 100% up to the $o = n$ point, after which the error starts to decrease gradually. This is because only $\frac{1}{n} = \frac{1}{{1024}}$ of the power is transmitted through the fiber, so the signal is lost in the noise. The NMSE improves only after the $o = n$ point, when input patterns are repeated. Although it may not be practically possible, the SNR can be improved by increasing the power transmitted through the fiber. Increasing the amplitude of the modeled input fields by a factor of $n = 1024$ gives an NMSE curve that exactly matches that of the Hadamard patterns. Figure 4(b) shows that using this transmission matrix can yield images that are recognizable, but are faint and have low contrast.

The sequence of random input patterns gives an NMSE curve that peaks at the $o = n$ point. This is because the random input patterns are not orthogonal, so overfitting occurs. The random input patterns can yield low errors when $o \gg n$. It can be seen from Fig. 4(c) that these transmission matrices yield images that seem noisy.

The sequence of Hadamard patterns gives a linearly decreasing NMSE trend for $o \lt n$, settling to an approximately steady value of ${\rm NMSE} = 3.81\%$ for $o \gt n$. This pattern significantly outperforms the other patterns considered in the presence of measurement noise. This can be attributed to both the fact that as much power as possible is transmitted through the scattering medium, hence maximizing the SNR, and the fact that these input patterns are mathematically orthogonal to each other, hence avoiding overfitting. These input patterns give the best imaging performance [Fig. 4(d)].

The results presented assumed a waveguide supporting 745 modes. Waveguides supporting both more modes and fewer modes were also modeled, and very similar trends were observed.

Inspection of Fig. 3(e) shows that orthogonal input patterns that couple as much power as possible through the scattering medium should be preferred when constructing a transmission matrix. In the case considered here, where the SLM plane is related to the fiber output plane, these requirements are met by the patterns based on the Hadamard matrix. Linear gratings meet the requirement of coupling as much power as possible through the fiber. However, these are not orthogonal, as they are truncated by the finite SLM size, and will therefore lead to overfitting.

Alternatively, a transmission matrix used to relate the input facet of the MMF to the output facet of the MMF might be desired. In our experimental setup, this corresponds to using the Fourier transform of the fields in the SLM plane as the input fields in the transmission matrix calculation. Additionally, physical Fourier transforms as performed by a lens are susceptible to vignetting and aberrations, which may also affect the analysis. The profile of input fields is hence altered, affecting their orthogonality as well as power coupling through the medium. Consequently, different input field sequences become feasible and optimal. For instance, scanning a diffraction-limited spot over the fiber input facet [26,27] satisfies orthogonality and power requirements. A basis consisting of fiber modes [28], Bessel modes [29], or Laguerre–Gauss [30] modes are orthogonal and hence not prone to overfitting. In theory, these can maximize the amount of power coupled through the system, but in practice generating such bases with an SLM cannot be done without inherent losses [31].

 

Fig. 5. (a) Singular values of ${\textbf{E}_{\textbf{out}}}$, measured for a ${60} \times {60}$ region on the interferometer using a sequence of 1024 Hadamard patterns. (b) Singular values of ${\textbf{E}_{\textbf{in}}}$, calculated for 1024 patterns of ${32} \times {32}$ hologram pixels, where different sequences of input patterns have been considered.

Download Full Size | PDF

An additional consideration is ensuring that the entire facet and numerical aperture are sampled at the input and output planes. This corresponds to selecting appropriate lenses as well as suitable holograms. For input patterns based on a Hadamard matrix, the overall size of the hologram as well as the macropixel size must be optimized. Similarly, when using input patterns based on linear gratings, the overall size and range of tilt angles must be appropriately selected. Alternatively, this can be thought of as selecting a basis that samples all modes guided by the waveguide. This can be achieved by using waveguide modes, Bessel modes, or Laguerre–Gauss modes with sizes tailored to the waveguide diameter.

6. EFFECT OF MEASUREMENT NOISE ON THE CONSTRUCTION OF THE INVERSE TRANSMISSION MATRIX

Once a sequence of input fields and corresponding output fields has been measured, the inverse transmission matrix can be calculated using Eq. (4). The problem is not necessarily of full rank, though, and can therefore be sensitive to measurement noise. We consider the singular value decomposition of a matrix $\textbf{M}$, such that $\textbf{M} = \textbf{U}{\boldsymbol \Sigma}{\textbf{V}^\textbf{*}}$, where $\textbf{U}$ and $\textbf{V}$ are unitary rotation matrices and ${\boldsymbol \Sigma}$ is a diagonal matrix with decreasing singular values along the diagonal. The Moore–Penrose inverse of $\textbf{M}$ can be calculated from ${\textbf{M}^\dagger} = \textbf{V}{{\boldsymbol \Sigma}^\dagger}{\textbf{U}^\textbf{*}}$. The Moore–Penrose inverse of ${\boldsymbol \Sigma}$ is calculated by replacing each diagonal element (singular value) with its inverse. This highlights the cause of instabilities when inverting a matrix: small singular values significantly impact the inverse transmission matrix and will be disproportionately affected by noise.

Previously, this has been resolved by using the conjugate transpose transmission matrix instead of the inverse transmission matrix [32], or by using a Tikhonov-regularized inverse transmission matrix [7], which strikes a balance between the inverse transmission matrix and the conjugate transpose transmission matrix [33]. A Tikhonov-regularized inverse matrix is calculated using an inverse ${\boldsymbol \Sigma}$ matrix where the diagonal elements are instead calculated using

Here, $\alpha$ is a user-defined parameter that is not known beforehand. $\alpha$ has been selected in other works by sweeping through a range of values, and inspecting the so-called L-curve [34], or by imaging a target and comparing the results obtained with the known ground truth [7]. An alternative approach based on a truncated singular value decomposition (TSVD) is proposed in this paper. This entails setting any singular values below a given threshold $\tau$ to zero, and discarding the corresponding rows and columns of $\textbf{U}$ and $\textbf{V}$. Again, it is not clear what threshold to use, but it shall be shown that the noise characteristics of the experimental setup can inform the optimal value of $\tau$.

Inspection of Eq. (4) reveals that $\textbf{E}_{\textbf{out}}^\dagger$ can be inverted instead of ${\textbf{T}^\dagger}$, in which case we can exploit the fact that the noise properties of ${\textbf{E}_{\textbf{out}}}$ are known to determine what value of $\tau$ should be used. An appropriate approach might hence consist of discarding all singular values of ${\textbf{E}_{\textbf{out}}}$ below the noise floor. Gavish and Donoho [35] show that this is sub-optimal though, as singular values marginally above the threshold may still be affected by noise. Instead, Gavish and Donoho show that a singular value threshold of ${\tau ^*}(\sigma) = \lambda (\beta)\sqrt n \sigma$ is optimal for an $m \times n$ matrix where $\beta = \frac{m}{n}$, and $\lambda$ is calculated from

While Gavish and Donoho assume additive Gaussian white noise, their approach has been reported to be equally effective when other types of noise are considered [36].

An experimental example of this is now considered; 1024 Hadamard input patterns are used to characterize the MMF under test. The matrix ${\textbf{E}_{\textbf{out}}}$ is hence of size $3600 \times 1024$, and as this waveguide supports only 745 modes in the single polarization considered, it can be anticipated that, beyond some threshold, the singular values will be dominated by noise. This size matrix gives $\beta \approx 3.5$, and $\lambda \approx 3.36$. From Fig. 2(a), the standard deviation of the noise on the beam is estimated to be at most ${\sigma _{\rm{err}}} = 0.8$, and an optimal threshold ${\tau ^*} \approx 86$ is hence obtained. Singular values of ${\textbf{E}_{\textbf{out}}}$ below this threshold are discarded, as illustrated in Fig. 5(a), and Eq. (4) is used to calculate the inverse transmission matrix.

Next, the amplitude image shown in Fig. 6(a) is displayed on the SLM. The regularized inverse transmission matrix can be used to recover the image from the complex speckle pattern captured by the interferometer. This gives the result shown in Fig. 6(b) where the spade symbol is clearly visible. For comparison, the image obtained using an inverse transmission matrix constructed without regularization is shown in Fig. 6(c), and the image obtained using the conjugate transmission matrix is shown in Fig. 6(d).

 

Fig. 6. (a) Imaging target displayed on the SLM. (b) Image obtained using inverse transmission matrix constructed with regularization. (c) Image obtained using inverse transmission matrix constructed without regularization. (d) Image obtained using conjugate transmission matrix.

Download Full Size | PDF

To illustrate the strengths of this technique, an image with more pixels than the number of modes guided by the waveguide was considered. An improved image will be obtained if the number of modes guided by the waveguide exceeds the number of macropixels imaged. Alternatively, this technique could have been tested by generating a diffraction-limited spot at the fiber output facet; however, this would require only a single row of the transmission matrix to be correct, and would therefore not test the validity of the entire transmission matrix.

To validate this approach to setting the TSVD threshold, $\tau$ has been swept through a range of values, with the obtained image reconstruction NMSE values shown in Fig. 7. It can be seen that the ideal threshold as calculated from the experimental noise in the system, indicated by a dotted line, does indeed correspond to the minimum in the NMSE curve. The image error obtained using a Tikhonov regularized inverse transmission matrix is also shown, and it is interesting to note that it follows a very similar trend, with neither approach significantly out-performing the other.

 

Fig. 7. Image reconstruction error for different amounts of regularization.

Download Full Size | PDF

The dominant computational overhead for the proposed approach is that of calculating the singular value decomposition of the dense matrix of output fields, done in $O(mn\max (m,n))$ [37]. This is potentially significantly more than the $O(mn)$ overhead required for calculating the conjugate transpose, but is equal to the overhead required for calculating a pseudoinverse. This is because most implementations of the pseudoinverse use a singular value decomposition.

In the analysis, ${\textbf{E}_{\textbf{in}}}$ has been assumed to have no uncertainty, with effects such as flicker due to pulse width modulation of the liquid crystal being accounted for in the noise on ${\textbf{E}_{\textbf{out}}}$. As such, the pseudoinverse of ${\textbf{E}_{\textbf{in}}}$ can be directly taken when calculating $\textbf{T}$. Should uncertainties be assigned to ${\textbf{E}_{\textbf{out}}}$, then regularization may again be required. Inspection of Fig. 5(b) reveals that this is the case for non-orthogonal patterns such as the random basis, but that the constant singular values of the Hadamard and macropixel patterns mean that regularization would not be necessary in these cases. This hence illustrates another strength of the Hadamard basis.

7. CONCLUSION

To the authors’ knowledge, this paper is the first to consider the effects of measurement noise on the construction of a scattering medium transmission matrix. The experimental noise for a representative interferometric measurement system has been quantified and used to model how the choice of input patterns affects the quality of the constructed transmission matrix. Superior results were obtained for orthogonal input patterns that couple maximum power through the scattering medium. It was also shown that the construction of the inverse transmission matrix is particularly susceptible to measurement noise, and regularization is required to obtain results with low errors. Importantly, we have demonstrated that the optimal amount of regularization can be estimated from the noise properties of the experimental system, rather than by trial and error, which significantly speeds up the regularization procedure.