## Abstract

Characterizing the transmission matrix (TM) of a multimode fiber (MMF) benefits many fiber-based applications and allows in-depth studies on the physical properties. For example, by modulating the incident field, the knowledge of the TM allows one to synthesize any optical field at the distill end of the MMF. However, the extraction of optical fields usually requires holographic measurements with interferometry, which complicates the system design and introduces additional noise. In this work, we developed an efficient method to retrieve the TM of the MMF in a referenceless optical system. With pure intensity measurements, this method uses the extended Kalman filter (EKF) to recursively search for the optimum solution. To facilitate the computational process, a modified speckle-correlation scatter matrix (MSSM) is constructed as a low-fidelity initial estimation. This method, termed EKF-MSSM, only requires 4*N* intensity measurements to precisely solve for *N* unknown complex variables in the TM. Experimentally, we successfully retrieved the TM of the MMF with high precision, which allows optical focusing with the enhancement (>70%) close to the theoretical value. We anticipate that this method will serve as a useful tool for studying physical properties of the MMFs and potentially open new possibilities in a variety of applications in fiber optics.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

When coherent monochromatic light is directed into a multimode fiber (MMF), a certain number of spatial modes are usually excited simultaneously. These modes propagate along the MMF with different wave vectors, exhibiting complicated interference patterns. Besides, inevitable bending and twisting of the MMF also causes strong mode couplings, which further complicates output patterns. As a result, speckle patterns, which are randomly distributed bright and dark spots, are observed at the distil end of the MMF. This phenomenon, termed mode crosstalk, significantly restricts the bandwidth of the MMF and therefore prohibiting it from being used in endoscopy for image delivery.

Wavefront shaping technique emerges as a promising tool to compensate this side effect in MMFs. By modulating the incident light with a special wavefront, it has been shown that arbitrarily designed patterns, such as a bright optical focus, can be synthesized at the distil end of the MMF. There are in general three different approaches to characterize the special wavefront: optical time reversal/optical phase conjugation [1–13], feedback-based wavefront shaping [14–21], and the transmission matrix (TM) measurement [22–27]. Among all three techniques, phase conjugation determines the special wavefront in a single measurement. However, it requires a rather complicated geometry with optical interferometry and a precise alignment between a camera and a spatial light modulator (SLM) [28–31]. Feedback-based wavefront shaping, on the other hand, is referenceless and has the simplest realization. By iteratively updating the modulated wavefront using a greedy optimization algorithm, an optical focus at any desired position can be progressively formed. However, constant information exchange between the SLM and the computer makes this approach less efficient. Alternatively, we can directly measure the TM of the MMF, which enables synthesizing any optical field at the distill end of the MMF. Moreover, the TM also provides rich information, allowing one to study physical properties, including absorption, reflection, and transmission [32–38]. However, the conventional approach to measure optical field usually relies on holographic methods, leading to complicated optical systems that are sensitive to ambient noises.

Common-path interferometry is a widely used approach to minimize the disturbance from ambient environment. By separating part of the active area of the SLM as a static reference, it has been shown that the TM of a thick scattering medium in free space can be characterized [22,39]. This method, however, sacrifices the information carried by the static reference, which is not ideal for characterizing the TM of MMFs with finite spatial modes. Therefore, a method that completely abandons a static reference is highly desirable. Without the reference part, retrieving the TM from pure intensity measurements is essentially solving for the optimum solution of a set of nonlinear equations. To achieve this goal, prVBEM method was demonstrated to extract the TM of a scattering medium in a referenceless optical system using a digital micro-mirror device (DMD) [40,41]. Although being computationally fast, this algorithm sometimes converges to a local optimum solution when the number of measurements is less than 6*N* for solving *N* unknown complex variables, resulting in inaccurate TM. This algorithm was further modified to be implemented in a binary TM format to reduce computational complexity [42]. Semi-definite programming (SDP) was used as a phase retrieval method to compute the TM with an SLM that continuously modulate phases from 0 to 2π [43,44]. Notably, SDP retrieves the TM with sufficient precision, but it requires about *N*ln*N* measurements to solve for *N* unknown complex variables. This nonlinear relationship makes SDP difficult to handle a large TM. Moreover, it was found that the cvx package implemented in previous works caused huge computational burden, which further limits its practical usage [43,44]. Recently, neural networks were also applied to handle similar problems [45–47]. Despite the demonstrated performance, this approach usually does not provide the exact values of the retrieved TM, but instead retrieve the substitute of the TM that is application oriented, such as delivering images.

To overcome some of the limitations that previous methods have, in this work, we developed an efficient method to retrieve the TM of the MMF. With pure intensity measurements, this method uses the extended Kalman filter (EKF) [48] to recursively search for the optimum solution. To facilitate the computational process, a modified speckle-correlation scatter matrix (MSSM) [49] is constructed as a low-fidelity initial estimation. This method, termed EKF-MSSM, only requires 4*N* intensity measurements to solve for *N* unknown complex variables in the TM. This linear relationship makes our proposed method more efficient to the SDP that requires *N*ln*N* measurements. Using numerical simulations, we examined the performance of the EKF-MSSM by retrieving the TM for a wide range of *N* from 100 to 1200. Statistically, we found that the EKF-MSSM always converges to the optimum solution with a probability higher than 96.34%. Further experimental validation was performed in a referenceless optical system. For *N* = 520, we successfully retrieved different rows of the TM of an MMF and realized optical focusing at different positions through the MMF. After loading conjugated phase values of the retrieved rows of the TM to the SLM, high-fidelity optical foci with enhancement ranging from 286 to 362, which account for 70% to 89% of respective theoretical values, have been constantly obtained.

## 2. Principle

#### 2.1 Establishing the EKF solver

Suppose we have *M* nonlinear equations with *N* unknown variables, these equations can be expressed as ${{y}^{({l})}} = {{F}^{({l} )}}({{x}_1},{{x}_2}, \cdots ,{{x}_{{N}}}),\; \;{\rm l} = 1, \cdots ,{M}$. We establish the EKF solver to search for the optimum solution using a recursive approach. We start by constructing the difference equation

*k*denotes the index of iterations. ${\textbf{X}_k} = {({x_{1,k}},{x_{2,k}}, \cdots {x_{n,k}})^{\textbf{T}}}$ is the true solution at the

*k*-th iteration, and the operator ${\textbf{T}}$ denotes matrix transpose. Ideally, the true solution ${\textbf{X}_k}$ should be a constant vector that is invariant to

*k*. However, process noise at the

*k*-th iteration ${\textbf{w}_k}$ inevitably causes small variations on ${\textbf{X}_k}$. Then, we construct the measurement equation as where

*k*-th iteration, which needs to be updated at every iteration. All variables with a “hat” denote the estimated values. ${\textbf{v}_k}$ represents measurement noise at the

*k*-th iteration. The ultimate goal of the iteration process is to make $\Vert{\hat{\textbf{Y}}_k} - {\hat{\textbf{C}}_k}{\hat{\textbf{X}}_k}\Vert_2$ approaches to zero. We note here that through the first-order Taylor expansion adopted in the measurement equation, any to be solved nonlinear equations can be linearized.

With the framework of the EKF introduced, detailed implementations of the iteration process are described in a flowchart shown in Fig. 1. In the first step, we generate the initial estimation ${\hat{\textbf{X}}_0}$ and initialize three matrices ${\textbf{P}_0}$, **R** and **Q**. Here, ${\textbf{P}_0}$ is the estimate error covariance for ${\hat{\textbf{X}}_0} - {\textbf{X}_0}$. A maximal number of *m* = 1000 iterations is set if the convergence test fails. Since we do not have any prior knowledge about ${\textbf{X}_0}$, ${\textbf{P}_0}$ is set to be an identity matrix **I**_{N}_{×}* _{N}*.

**R**and

**Q**are the covariance of $\textbf{w}$ and $\textbf{v}$, respectively, which are set to be

**R**= 10

^{−4}

**I**

_{M}_{×}

*,*

_{M}**Q**= 10

^{−4}

**I**

*in the beginning. The values for $\textbf{P}$ will be updated at each iteration. Iteration process will be executed based on the four equations described in the second box. These four equations have been commonly used in the EKF for kinetic tracking, model estimation, and animation. Rigorous mathematical derivation of these equations is neglected here but can be found in Ref. [48]. This iteration process is repeated until the correlation coefficient between ${\hat{\textbf{X}}_{k - 3}} + {\hat{\textbf{X}}_{k - 2}}$ and ${\hat{\textbf{X}}_{k - 1}} + {\hat{\textbf{X}}_k}$ reaches 99.9999% or*

_{N×N}*k = m*. If the termination criterion is satisfied, the one that minimizes $\Vert{\hat{\textbf{Y}}_k} - {\hat{\textbf{C}}_k}{\hat{\textbf{X}}_k}\Vert_2$ is output as the final estimation of $\hat{\textbf{X}}$.

#### 2.2 Retrieving the TM using the EKF solver

The TM of the MMF is assumed to be a square matrix with a dimension of *N*×*N* for simplicity. Moreover, the uncorrelated transmission coefficients (UTC) model is employed in this work, in which each element ${t}({{p},{q}} )= {a}({{p},{q}} )+ {ib}({{p},{q}} )$ satisfies a circular Gaussian distribution [50]. Since most commercialized liquid crystal based SLMs only support phase-only modulation and in order to mimic realistic experimental environment, each element of the input field is assumed to have the same amplitude but random phase values ${\left[ {\begin{array}{ccc} {{{e}^{{i}{{\theta }_1}}},}&{ \cdots ,}&{{{e}^{{i}{{\theta }_{{N}}}}}} \end{array}} \right]^{\textbf{T}}}$. After passing through the MMF, the *p*-th element in the output field is

*p*-th row of the TM.

From now on, we restrict ourselves to solve for the *p*-th row (*p* is a dummy variable) of the TM. With this premise, the notation of $a({p,q} )$, $b({p,q} )$, and $E(p )$ can be simplified to $a(q ),$ $b(q )$, and *E*, respectively. To estimate the values for $a(q )$ and $b(q )$, *M* different input fields are generated and their corresponding intensity measurements at the output end of the MMF are measured. In this condition, a set of *M* nonlinear equations is established to solve for 2*N* (including both real and imaginary parts) unknown variables:

*l*denotes the index of measurements/equations, ${y^{(l )}}$ is the intensity at the

*p*-th pixel when the

*l-*th training phase pattern is applied. To use the EKF solver, we substitute Eq. (5) into Eqs. (1)–(3) and get

*M*= 4

*N*intensity measurements are sufficient to retrieve one row of the elements in the TM (

*N*elements, 2

*N*unknown variables). In other words, the number of rows in the TM times 4

*N*becomes the required number of intensity measurements to retrieve the complete TM. The number of input trials can still be kept at 4

*N*if parallel measurements are performed by detector arrays.

#### 2.3 Constructing the MSSM as the initial estimation

A good choice of the initial estimation ${\hat{\textbf{X}}_0}$ can effectively facilitate the computational process. Inspired by the SSM theory [49], we define a **Z** matrix

*m*-th column of a

**P**matrix. Different from the SSM that constructs the

**P**matrix using a known TM [49], we construct the

**P**matrix from the known input fields, as follows

**Z**matrix, the complex eigenvector that has the maximum eigenvalue is chosen as the initial estimation ${\hat{\textbf{X}}_0}$. In the following simulation, we will show that such choice can greatly facilitate the computational process.

## 3. Numerical results

We start by discussing how many intensity measurements/equations are required for the EKF-MSSM method to accurately retrieve one row of the TM (2*N* unknown variables). From now on, the discussions will focus on retrieving only one row of the TM for simplicity, however, we emphasize that the EKF-MSSM can work for retrieving the entire TM. In general, for linear equations, 2*N* independent measurements/equations can provide enough information to estimate the values for 2*N* independent variables. However, all the equations here have quadratic form, indicating more intensity measurements/equations are needed. Following the strategy in Ref. [40], we also define $\gamma = M/N$. In general, a small $\gamma $ may lead to a fast convergence to a local optimal solution. A large $\gamma $, on the other hand, increases the probability to converge to the global optimal solution, but it may also introduce redundant information and consume unnecessary computational resources. To find an appropriate $\gamma $ for the EKF-MSSM, we quantified the performance of the proposed method under different values of $\gamma $, using the quality of the optical focus formed with the estimated values of the TM. In practice, we simply sent in the conjugation of the elements in the TM for focusing. The more accurate the retrieved TM is, the better quality the focus has. Following the convention of wavefront shaping, the quality of the focus is quantified by the enhancement, which is defined as the ratio between the intensity of focus after optimization *I*_{foc} and the initial ensemble averaged intensity before optimization ${\langle{I}_{\textrm{avg}}\rangle}$ [14]_{:}

*N*is the number of degrees of freedom in the system, and the factor

*α*∈[0,1] depends on the types of modulation. For the phase-only modulation used in this work,

*α*is

*π*/4.

Figure 2(a) shows the achieved enhancement as a function of $\gamma $, when the TM is retrieved by the EKF, the MSSM, and the EKF-MSSM, respectively. To generate this plot, *N* is fixed at 100, and 500 independent runs were averaged to minimize statistical error. The purple line denotes the theoretical value of about 78. The light blue line denotes the performance of the EKF when the initial estimation ${\hat{\textbf{X}}_0}$ is randomly chosen. As we can see from the figure that the enhancement is quite low for small $\gamma $, but it gradually increases with increasing $\gamma $. At $\gamma = \textrm{7}$, the enhancement reaches 98.85% of the theoretical value, which implies a high probability to converge to a global optimal solution. This observation indicates that the EKF can easily trapped into local optimum solution when $\gamma $ is not large enough. This problem can be potentially alleviated by trying different initial estimation, at the cost of computational time. Interestingly, the MSSM, shown as the black line in Fig. 2(a), can also generate optical foci with relatively low enhancement. Even at $\gamma $ = 10, the MSSM can only generate optical focus with the enhancement reaching 70.3% of the theoretical value. Surprisingly, if we choose the MSSM as the initial low-fidelity estimation ${\hat{\textbf{X}}_0}$ for the EKF, the performance of the EKF can be greatly improved. This method, termed EKF-MSSM [dark blue line in Fig. 2(a)], can generate optical focus with the enhancement reaches 99.53% of the theoretical value at $\gamma = 4$. Moreover, we check the performance of the EKF-MSSM for different *N* ranging from 100 to 1200, and the results are shown in Fig. 2(b). $\gamma = 4$ is fixed and 100 independent runs were performed for averaging in each case. Specifically, at *N* = 100, 250, 400, 650, 900, and 1200, the enhancement can reach 99.25%, 98.30%, 99.07%, 99.15%, 100% and 96.34% of the theoretical value, respectively, indicating that the EKF-MSSM can always retrieve the TM with high precision. Therefore, we conclude that *M* = 4*N* ($\gamma = 4$) is sufficient to solve for one row of the TM (2*N* unknown real variables) by using the EKF-MSSM. As a comparison, the MSSM, denoted by the black line, always serves as a low-fidelity initial estimation ${\hat{\textbf{X}}_0}$. Figures 2(c) and 2(d) further show the histogram of the enhancement and the number of iterations cost for the EKF and the EKF-MSSM. These data correspond to $\gamma = 4$ case in Fig. 2(a). As we can see from the figure that, the EKF-MSSM costs a lot less iterations and is more likely to converge to the global optimum solution, compared to the EKF with an arbitrarily chosen initial estimation.

We further show the computational time consumed by the EKF-MSSM in Fig. 3(a). We note that unless otherwise specified, all the simulation results were obtained with a computer having an i5-8600k 3.6 GHz CPU and a 16-GB RAM. No independent graphics card assistance was added. 100 independent runs were used for averaging. For example, when *N* = 100 and $\gamma = 4$, the EKF-MSSM only takes about 0.16 second to retrieve one row of the TM. Although the computational time dramatically increases as *N* becomes large, a personal computer with 16-GB RAM can still handle the computational burden with *N* up to a few thousands, which is enough for most MMFs. We speculate that the nonlinear increase in the computational time is possibly due to the shortage of available physical computational memory, which is essential in matrix operation. In contrast, the SDP [43,44] requires about $\gamma = \textrm{ln}N$ to recover the TM at huge computational burden. Even replaced with $\gamma = 4$, the SDP takes about 28s seconds at *N* = 100 [estimated using our computer, shown as a red star in Fig. 3(a)]. When *N* goes beyond 150, the SDP stops working as the embedded cvx package crashes due to the shortage of computational memory. The prVBEM with binary-intensity measurements is substantially faster, which requires 1 to 3 minutes to retrieve one row of TM at *N* = 1296 and $\gamma = 6$ [estimated in Ref. [41], shown as a black star in Fig. 3(a)]. Although a direct observation indicates the prVBEM is slightly efficient than the EKF-MSSM, we note that the value of the prVBEM was obtained using high-performance computing facilities. The existence of measurement noise will certainly deteriorate the performance of the EKF-MSSM. We also exam the performance of the EKF-MSSM in noisy conditions when *N* = 100, $\gamma = 4$. As with any method, the existence of noise can prohibit the EKF-MSSM from converging to the global optimum solution. Figure 3(b) plots the normalized enhancement (with respect to the theoretical value) as a function of the signal-to-noise ratio (SNR) in blue. The computational time is not largely affected by the existence of the noise. When the SNR is larger than 10, we can see that the normalized enhancement can maintain at a relatively high level about 0.88. Beyond this level, the normalized enhancement gradually decreases, but can still be around 0.25 when the SNR is as low as 1. As a comparison, the performance of the SDP under the same condition was also plotted in red. When the SNR is larger than 5, the EKF-MSSM slightly outperforms the SDP.

Another interesting question for the EKF-MSSM is what will happen if the phase values used for generating input trials are discrete with finite numbers. To answer this question, we investigate the influence on the enhancement by using both the EKF-MSSM and the SDP through numerical simulations. Figure 4(a) plots the normalized enhancement as a function of the number of discrete phase values *d*. For simplicity, the phase values are equally distributed. Numerical results show that both the EKF-MSSM and the SDP cannot accurately retrieve the TM when *d* = 2, leading to reduced enhancement. Such a condition cannot be improved even if we enlarge the $\gamma $ to 10 (not shown in the figure). When *d* becomes larger than or equal to 2, both the EKF-MSSM and the SDP retrieve the TM with high accuracy. As a result, full-phase modulation-based algorithms do not have strict restrictions on the phase levels provided by the SLM. Figure 4(b) further shows the computational time cost for each case, and it can be observed that the computational time is not largely affected by the choice of *d*. Interestingly, *d* = 2 takes the longest time among all cases, possibly because insufficient information needs more iterations to converge. Besides considering SLMs with finite number of phase values, we also considered the possibility of using a DMD with “on” and “off” states. Unfortunately, however, the EKF-MSSM cannot be directly extended as Eqs. (5), (8), and (9) become ambiguous.

## 4. Experimental results

Figure 5(a) illustrates the experimental procedure. To retrieve the TM that has a dimension of *N*${\times} $*N*, 4*N* input trials were sent into the MMF and their corresponding 4*N* output speckle patterns were measured. The phase values of the 4*N* input trails are randomly chosen within a range from 0 to 2π. The detailed experimental setup is shown in Fig. 5(b). A continuous wave laser (MDL-C-642-30mW) that operates at 642 nm was used as the light source. A half-wave plate and a polarizing beam splitter were used to adjust the optical power of the system. The beam was then expanded by a pair of lenses L1 and L2 to cover the active area of a phase-only SLM (PLUTO2-NIR, Holoeye, 1920${\times} $1080, 8 µm/pixel). After being modulated and reflected by the SLM, the light was directed by a beam splitter to a 1-meter long MMF. The step-index MMF is silica based (refractive index ∼ 1.48), and has a core diameter of 50${\pm} 2.5$ µm and a numerical aperture of 0.200${\pm} 0.01$ (FC/PC-FC/PC-50/125-900µm-1m, Shenzhen Optics-forest Inc), which supports about 1200 modes at 642nm. Two objective lenses were used to couple light into and to collect the output light from the MMF. An 8-bit CCD (Point Grey Research Grasshopper3 GS3-U3-32S4C) captured the output speckle pattern with RGB 1024 × 768 Mode1. Both the CCD and the SLM were controlled by a computer. Without intended interventions, the system is quite stable and the speckle correlation can be kept above 90% even after 1 hour. This time scale is much longer than the typical operation time of algorithm, which guarantees the accuracy of the retrieved TM.

In experiments, we equally subdivided the active area of the SLM into $26 \times 20 = 520$ independent segments. The pixels within the same segment were modulated with the same phase. A sequence of $4 \times 520 = 2080$ input trials (*N* = 520), with random phase values ranging from 0 to 2*π*, were sent into the MMF. Figure 6(a) shows a typical speckle pattern observed by the CCD, which exhibits randomly distributed bright and dark spots. The corresponding 2080 intensity measurements were performed by the CCD. As a typical example, the intensity of a target pixel on the CCD, labeled as pixel 1, was plotted in Fig. 6(b). These intensity measurements were then used to retrieve one row of the TM by the EKF-MSSM. As the first demonstration, the retrieved TM enables optical focusing through the MMF at the target pixel (pixel 1 or P1), and the CCD captured image is shown in Fig. 6(c). We note that the intensity of the formed focus, denoted by a black dot encircled in red in Fig. 6(b), is much higher than the averaged intensity of the previous measurements [red solid line in Fig. 6(b)]. The enhancement of the focus $\eta \approx 286$, which is about 70% of the theoretical value. As a comparison, Fig. 6(d) shows the CCD captured focus at pixel 1 (P1’), when the TM was retrieved by the MSSM. In this case, weak speckle background can still be seen, and the enhancement $\eta \approx 76$, which is only 18% of the theoretical value.

With the same input trails, we further show focusing light through the MMF at different positions. With the TM retrieved by the EKF-MSSM, Figs. 7(a) and 7(b) show the CCD captured images of the optical foci formed at two different positions, i.e., labeled as pixel 2 and pixel 3. The enhancement of these two foci reach about 76% (P2) and 89% (P3) of the theoretical value. As a comparison, when the TM was retrieved using the MSSM with the same set of data, the formed optical foci, shown in Figs. 7(c) and 7(d), exhibit much poorer contrast. Quantitatively, the enhancement of these two foci only reach about 17% (P2’) and 13% (P3’) of the theoretical value. The above results clearly demonstrate the superior performance of the EKF-MSSM in retrieving the TM in a referenceless optical system. As a final remark, we recall that the existence of measurement noise, contributed from inaccurate phase values of the SLM, laser intensity fluctuations, electronic noises of the CCD, and time-varying characteristics of the MMF can downgrade the performance of the EKF-MSSM. These realistic conditions suggest that further development on the algorithm should continuously focus on the robustness, efficiency, and accuracy of the algorithm. Nonetheless, we still retrieved the TM of the MMF with high precision, which allows optical focusing with the enhancement close (>70%) to the theoretical value.

## 5. Conclusion

In summary, we developed a non-holographic method, termed the EKF-MSSM, to retrieve the TM of the MMF. Using numerical approaches, we demonstrated that 4*N* intensity measurements are enough for the EKF-MSSM to accurately retrieve the TM. The feasibility of the EKF-MSSM was further demonstrated experimentally. For *N* = 520, we built a referenceless optical system, successfully retrieved the TM, and demonstrated optical focusing through the MMF at different pixels of the CCD. These optical foci exhibit good quality, exhibiting the enhancement larger than 70% of the theoretical value. Due to the superior performance and relatively simple framework of the EKF-MSSM, this method has the potential to become a powerful tool in biomedical imaging by efficiently calibrating MMF-based endoscopes.

## Funding

National Natural Science Foundation of China (61605252).

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **M. Cui and C. Yang, “Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation,” Opt. Express **18**(4), 3444–3455 (2010). [CrossRef]

**2. **Y. M. Wang, B. Judkewitz, C. A. DiMarzio, and C. Yang, “Deep-tissue focal fluorescence imaging with digitally time-reversed ultrasound-encoded light,” Nat. Commun. **3**(1), 928 (2012). [CrossRef]

**3. **K. Si, R. Fiolka, and M. Cui, “Fluorescence imaging beyond the ballistic regime by ultrasound-pulse-guided digital phase conjugation,” Nat. Photonics **6**(10), 657–661 (2012). [CrossRef]

**4. **T. R. Hillman, T. Yamauchi, W. Choi, R. R. Dasari, M. S. Feld, Y. Park, and Z. Yaqoob, “Digital optical phase conjugation for delivering two-dimensional images through turbid media,” Sci. Rep. **3**(1), 1909 (2013). [CrossRef]

**5. **E. H. Zhou, H. Ruan, C. Yang, and B. Judkewitz, “Focusing on moving targets through scattering samples,” Optica **1**(4), 227–232 (2014). [CrossRef]

**6. **D. Wang, E. H. Zhou, J. Brake, H. Ruan, M. Jang, and C. Yang, “Focusing through dynamic tissue with millisecond digital optical phase conjugation,” Optica **2**(8), 728–735 (2015). [CrossRef]

**7. **H. Ruan, M. Jang, and C. Yang, “Optical focusing inside scattering media with time-reversed ultrasound microbubble encoded light,” Nat. Commun. **6**(1), 8968 (2015). [CrossRef]

**8. **B. Judkewitz, R. Horstmeyer, I. M. Vellekoop, I. N. Papadopoulos, and C. Yang, “Translation correlations in anisotropically scattering media,” Nat. Phys. **11**(8), 684–689 (2015). [CrossRef]

**9. **E. E. Morales-Delgado, S. Farahi, I. N. Papadopoulos, D. Psaltis, and C. Moser, “Delivery of focused short pulses through a multimode fiber,” Opt. Express **23**(7), 9109–9120 (2015). [CrossRef]

**10. **Y. Liu, C. Ma, Y. Shen, J. Shi, and L. V. Wang, “Focusing light inside dynamic scattering media with millisecond digital optical phase conjugation,” Optica **4**(2), 280–288 (2017). [CrossRef]

**11. **Y. Shen, Y. Liu, C. Ma, and L. V. Wang, “Sub-Nyquist sampling boosts targeted light transport through opaque scattering media,” Optica **4**(1), 97–102 (2017). [CrossRef]

**12. **Z. Yu, J. Huangfu, F. Zhao, M. Xia, X. Wu, X. Niu, D. Li, P. Lai, and D. Wang, “Time-reversed magnetically controlled perturbation (TRMCP) optical focusing inside scattering media,” Sci. Rep. **8**(1), 2927 (2018). [CrossRef]

**13. **C. Ma, J. Di, Y. Zhang, P. Li, F. Xiao, K. Liu, X. Bai, and J. Zhao, “Reconstruction of structured laser beams through a multimode fiber based on digital optical phase conjugation,” Opt. Lett. **43**(14), 3333–3336 (2018). [CrossRef]

**14. **I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. **32**(16), 2309–2311 (2007). [CrossRef]

**15. **I. M. Vellekoop and A. P. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Commun. **281**(11), 3071–3080 (2008). [CrossRef]

**16. **O. Katz, E. Small, Y. Bromberg, and Y. Silberberg, “Focusing and compression of ultrashort pulses through scattering media,” Nat. Photonics **5**(6), 372–377 (2011). [CrossRef]

**17. **D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, “Genetic algorithm optimization for focusing through turbid media in noisy environments,” Opt. Express **20**(5), 4840–4849 (2012). [CrossRef]

**18. **P. Lai, L. Wang, J. W. Tay, and L. V. Wang, “Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media,” Nat. Photonics **9**(2), 126–132 (2015). [CrossRef]

**19. **B. Blochet, L. Bourdieu, and S. Gigan, “Focusing light through dynamical samples using fast continuous wavefront optimization,” Opt. Lett. **42**(23), 4994–4997 (2017). [CrossRef]

**20. **Z. Wu, J. Luo, Y. Feng, X. Guo, Y. Shen, and Z. Li, “Controlling 1550-nm light through a multimode fiber using a Hadamard encoding algorithm,” Opt. Express **27**(4), 5570–5580 (2019). [CrossRef]

**21. **Q. Feng, F. Yang, X. Xu, B. Zhang, Y. Ding, and Q. Liu, “Multi-objective optimization genetic algorithm for multi-point light focusing in wavefront shaping,” Opt. Express **27**(25), 36459–36473 (2019). [CrossRef]

**22. **S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media,” Phys. Rev. Lett. **104**(10), 100601 (2010). [CrossRef]

**23. **T. Chaigne, O. Katz, A. C. Boccara, M. Fink, E. Bossy, and S. Gigan, “Controlling light in scattering media non-invasively using the photoacoustic transmission matrix,” Nat. Photonics **8**(1), 58–64 (2014). [CrossRef]

**24. **M. Kim, W. Choi, Y. Choi, C. Yoon, and W. Choi, “Transmission matrix of a scattering medium and its applications in biophotonics,” Opt. Express **23**(10), 12648–12668 (2015). [CrossRef]

**25. **J. Yoon, K. Lee, J. Park, and Y. Park, “Measuring optical transmission matrices by wavefront shaping,” Opt. Express **23**(8), 10158–10167 (2015). [CrossRef]

**26. **W. Xiong, P. Ambichl, Y. Bromberg, B. Redding, S. Rotter, and H. Cao, “Spatiotemporal Control of Light Transmission through a Multimode Fiber with Strong Mode Coupling,” Phys. Rev. Lett. **117**(5), 053901 (2016). [CrossRef]

**27. **J. Xu, H. Ruan, Y. Liu, H. Zhou, and C. Yang, “Focusing light through scattering media by transmission matrix inversion,” Opt. Express **25**(22), 27234–27246 (2017). [CrossRef]

**28. **M. Azimipour, F. Atry, and R. Pashaie, “Calibration of digital optical phase conjugation setups based on orthonormal rectangular polynomials,” Appl. Opt. **55**(11), 2873–2880 (2016). [CrossRef]

**29. **A. Hemphill, Y. Shen, J. Hwang, and L. Wang, “High-speed alignment optimization of digital optical phase conjugation systems based on autocovariance analysis in conjunction with orthonormal rectangular polynomials,” J. Biomed. Opt. **24**(03), 1 (2018). [CrossRef]

**30. **M. Jang, H. Ruan, H. Zhou, B. Judkewitz, and C. Yang, “Method for auto-alignment of digital optical phase conjugation systems based on digital propagation,” Opt. Express **22**(12), 14054–14071 (2014). [CrossRef]

**31. **Y.-W. Yu, C.-C. Sun, X.-C. Liu, W.-H. Chen, S.-Y. Chen, Y.-H. Chen, C.-S. Ho, C.-C. Lin, T.-H. Yang, and P.-K. Hsieh, “Continuous amplified digital optical phase conjugator for focusing through thick, heavy scattering medium,” OSA Continuum **2**(3), 703–714 (2019). [CrossRef]

**32. **Y. Chong and A. D. Stone, “Hidden black: Coherent enhancement of absorption in strongly scattering media,” Phys. Rev. Lett. **107**(16), 163901 (2011). [CrossRef]

**33. **M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park, and W. Choi, “Maximal energy transport through disordered media with the implementation of transmission eigenchannels,” Nat. Photonics **6**(9), 581–585 (2012). [CrossRef]

**34. **A. Yamilov, S. Petrenko, R. Sarma, and H. Cao, “Shape dependence of transmission, reflection, and absorption eigenvalue densities in disordered waveguides with dissipation,” Phys. Rev. B **93**(10), 100201 (2016). [CrossRef]

**35. **M. I. Mishchenko, L. Liu, and J. W. Hovenier, “Effects of absorption on multiple scattering by random particulate media: exact results,” Opt. Express **15**(20), 13182–13187 (2007). [CrossRef]

**36. **R. Sarma, A. Yamilov, and H. Cao, “Enhancing light transmission through a disordered waveguide with inhomogeneous scattering and loss,” Appl. Phys. Lett. **110**(2), 021103 (2017). [CrossRef]

**37. **Y. Huang, Y. Shen, C. Min, and G. Veronis, “Switching of the direction of reflectionless light propagation at exceptional points in non-PT-symmetric structures using phase-change materials,” Opt. Express **25**(22), 27283 (2017). [CrossRef]

**38. **Y. Huang, Y. Shen, and G. Veronis, “Non-PT-symmetric two-layer cylindrical waveguide for exceptional-point-enhanced optical devices,” Opt. Express **27**(26), 37494–37507 (2019). [CrossRef]

**39. **D. B. Conkey, A. M. Caravaca-Aguirre, and R. Piestun, “High-speed scattering medium characterization with application to focusing light through turbid media,” Opt. Express **20**(2), 1733–1740 (2012). [CrossRef]

**40. **A. Drémeau, A. Liutkus, D. Martina, O. Katz, C. Schülke, F. Krzakala, S. Gigan, and L. Daudet, “Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques,” Opt. Express **23**(9), 11898–11911 (2015). [CrossRef]

**41. **L. Deng, J. D. Yan, D. S. Elson, and L. Su, “Characterization of an imaging multimode optical fiber using a digital micro-mirror device based single-beam system,” Opt. Express **26**(14), 18436–18447 (2018). [CrossRef]

**42. **T. Zhao, L. Deng, W. Wang, D. S. Elson, and L. Su, “Bayes’ theorem-based binary algorithm for fast reference-less calibration of a multimode fiber,” Opt. Express **26**(16), 20368–20378 (2018). [CrossRef]

**43. **M. N’Gom, M.-B. Lien, N. M. Estakhri, T. B. Norris, E. Michielssen, and R. R. Nadakuditi, “Controlling light transmission through highly scattering media using semi-definite programming as a phase retrieval computation method,” Sci. Rep. **7**(1), 1–9 (2017). [CrossRef]

**44. **M. N’Gom, T. B. Norris, E. Michielssen, and R. R. Nadakuditi, “Mode control in a multimode fiber through acquiring its transmission matrix from a reference-less optical system,” Opt. Lett. **43**(3), 419–422 (2018). [CrossRef]

**45. **A. Turpin, I. Vishniakou, and J. d Seelig, “Light scattering control in transmission and reflection with neural networks,” Opt. Express **26**(23), 30911–30929 (2018). [CrossRef]

**46. **N. Borhani, E. Kakkava, C. Moser, and D. Psaltis, “Learning to see through multimode fibers,” Optica **5**(8), 960–966 (2018). [CrossRef]

**47. **P. Fan, T. Zhao, and L. Su, “Deep learning the high variability and randomness inside multimode fibers,” Opt. Express **27**(15), 20241–20258 (2019). [CrossRef]

**48. **G. Welch and G. Bishop, “An introduction to the Kalman filter,” (1995).

**49. **K. Lee and Y. Park, “Exploiting the speckle-correlation scattering matrix for a compact reference-free holographic image sensor,” Nat. Commun. **7**(1), 13359 (2016). [CrossRef]

**50. **J. W. Goodman, * Speckle phenomena in optics: theory and applications* (Roberts and Company Publishers, 2007).