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 6N 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 NlnN 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 4N 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 NlnN 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

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⁻⁴ I_M×M, Q = 10⁻⁴ I_N×N 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 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}}$.

 

Fig. 1. A flowchart of the EKF solver.

Download Full Size | PDF

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

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 2N (including both real and imaginary parts) unknown variables:

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

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 (2N 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, 2N independent measurements/equations can provide enough information to estimate the values for 2N 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]_:

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 = 4N ($\gamma = 4$) is sufficient to solve for one row of the TM (2N 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.

 

Fig. 2. Simulation results of the EKF-MSSM. (a) The enhancement as a function of $\gamma $, when the TM is retrieved by the EKF, the MSSM, and the EKF-MSSM, respectively. N = 100 is fixed. Error bar: standard error of 500 independent runs. (b) The enhancement as a function of the degrees of freedom N. $\gamma = 4$ is fixed. Error bar: standard error of 100 independent runs. (c) The histogram of the enhancement, based on the results from 500 independent runs with $\gamma = 4$. (d) The histogram of the number of iterations, based on the results from 500 independent runs with $\gamma = 4$.

Download Full Size | PDF

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.

 

Fig. 3. Quantification of the performance when $\gamma = 4$. (a) The computational time consumed by the EKF-MSSM as a function of the degrees of freedom N. Error bar: standard error of 100 independent runs. Two distinct data points are also shown to represent the performance of the SDP and the prVBEM. (b) Normalized enhancement (with respect to the theoretical value) as a function of the signal-to-noise ratio for the SDP (red) and the EKF-MSSM (blue). N = 100. Error bar: standard error of 200 independent runs.

Download Full Size | PDF

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.

 

Fig. 4. Quantification of the performance of the algorithms when the SLM can only modulate d discrete phases values. $\gamma = 4$, N = 100 are adopted to generate these results. Error bar: standard error of 200 independent runs. (a) Normalized enhancement (with respect to the theoretical value) as a function of d. (b) The computational time as a function of d.

Download Full Size | PDF

4. Experimental results

Figure 5(a) illustrates the experimental procedure. To retrieve the TM that has a dimension of N${\times} $N, 4N input trials were sent into the MMF and their corresponding 4N output speckle patterns were measured. The phase values of the 4N 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.

 

Fig. 5. Schematics of the experiments. (a) The experimental procedure. 4N input trials with random phase maps were sent into the MMF and their corresponding 4N output speckle patterns were measured. (b) The experimental setup. BB1, BB2: beam block; BS: beam splitter; CCD: charge-coupled device; HWP: half-wave plate; L1, L2: lens; M: mirror; MMF: multimode fiber; OBJ1, OBJ2: objective lens; PBS: polarizing beam splitter; SLM: spatial light modulator; P: polarizer.

Download Full Size | PDF

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.

 

Fig. 6. Experimental demonstration of retrieving the TM using the EKF-MSSM. (a) A typical speckle pattern observed by the CCD at the output of the MMF. (b) Distributions of 2080 intensity measurements at pixel 1. The averaged intensity is denoted by the red solid line. After retrieving the TM, the intensity of the focus is denoted by a black dot encircled in red. (c) The CCD captured focus at pixel 1, when the TM was retrieved by the EKF-MSSM. $\eta \approx 286\; ({\textrm{P}1})$. Computational time: about 23 seconds. (d) The CCD captured focus at pixel 1, when the row of TM was retrieved by the MSSM. $\eta \approx 76\; ({\textrm{P}1^{\prime}} )$, Computational time: about 10 seconds.

Download Full Size | PDF

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.

 

Fig. 7. Focusing light through the MMF at different positions with the computed TM. (a)(b) The CCD captured images of the optical foci formed at pixel 2 and pixel 3, when the TM were retrieved by the EKF-MSSM. The enhancement of these two foci reach about 76% (P2) and 89% (P3) of the theoretical value. Computational time: about 21 seconds and 25 seconds, respectively. (c)(d) The CCD captured images of the optical foci formed at pixel 2 and pixel 3, when the TM were retrieved by the MSSM. The enhancement of these two foci reach about 17% (P2’) and 13% (P3’) of the theoretical value. Computational time: both around 10 seconds.

Download Full Size | PDF

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 4N 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.