Characterization of an imaging multimode optical fiber using a digital micro-mirror device based single-beam system

: This work demonstrates experimental approaches to characterize a single multimode fiber imaging system without a reference beam. Spatial light modulation is performed with a digital micro-mirror device that enables high-speed binary amplitude modulation. Intensity-only images are recorded by the camera and processed by a Bayesian inference based algorithm to retrieve the phase of the output optical field as well as the transmission matrix of the fiber. The calculated transmission matrix is validated by three standards: prediction accuracy, transmission imaging, and focus generation. Also, it is found that information on mode count and eigenchannels can be extracted from the transmission matrix by singular value decomposition. This paves the way for a more compact and cheaper single multimode fiber imaging system for many demanding imaging tasks.


Introduction
Recently, single optical multimode fibers (MMF) have become a promising candidate for endoscopic imaging due to the great number of spatial modes travelling within a small crosssection down to tens of micrometers [1]. Typically, waves propagating in MMFs suffer from mode coupling and mode dispersion and generate complex speckle patterns, which imposes a fundamental challenge for a single MMF in imaging applications. Although the intermode dispersion and the mode coupling in an MMF scramble the transmitted image and lead to noise-like output speckle patterns, the light propagation inside a length of MMF can be predicted and the original image can be reconstructed. With appropriate characterization, MMFs can work similarly to conventional refractive imaging components such as optical lenses.
Over the past few years, wavefront shaping technology has been investigated extensively to characterize MMFs so that a single MMF can be used to transmit images. Light propagation in a linear medium is a linear process and can, therefore, be described by its transmission matrix (TM) [2,3], which is calculated from a set of system input-output field pairs. In [4], Choi et al. used a scanning mirror to modulate the input fields and obtained the output fields with off-axis holography. All the data was processed to construct the TM, leading to imaging with diffraction-limited resolution. However, this method necessitates the recording of input fields, which is time-consuming and lowers experiment operability. Alternatively, by offering phase-only modulations, liquid crystal based spatial light modulators (LC-SLM) are more widely used in wavefront shaping. Papadopoulos et al. demonstrated focus generation [5] and point scanning based fluorescence imaging [6] with digital phase conjugation using a LC-SLM. Loterie and associates utilized an LC-SLM to measure the TM of the MMF and further achieved high resolution imaging with a digital confocal [7] and optical confocal setup [8]. The importance of monitoring system stability and phase drift was stressed in their work, as phase correction is essential for TM effectiveness. Digital micromirror devices (DMD) have been used as cheap and high-speed alternatives to LC-SLM. The on board rotation of micro-mirrors can realize binary amplitude modulation, which can also be used to measure the fiber TM and to perform point-scanning [9]. The previously mentioned approaches all require a reference beam, either co-propagating in the fiber [2,3] or in an external form [4][5][6][7][8][9]. However, extracting the full information on the optical field by phase-shifted or off-axis holography comes at the cost of a complicated setup and suffers from instability issues arising from the holographic interference process for phase retrieval.
Since mode coupling is highly dependent on the change of fiber shape and movements, flexible single MMF imaging is a challenge. Normally, MMF imaging systems have limited flexibility, and it has been demonstrated that image transmission for graded-index MMF [10] and high-NA MMF [11] suffers less from fiber deformation. To improve this, Farahi et al. proposed a dynamic bending compensation system with the help of another single mode fiber and a virtual beacon source [12], but the MMF deformation is carefully limited in their setup and sufficient pre-acquired data is needed to perform the re-focusing. Moreover, recent developments in theoretical models enable the prediction of the fiber output after deformation [13], which requires an extremely accurate model to match the experimental results. Since the error accumulates with the image transmission distance, this method is only validated for short fibers (no more than 300 mm). Another related difficulty is that accurate knowledge is required of the fiber shape to model the propagation, which is impractical in real-world applications. To deal with truly free deformation, Caravaca-Aguirre et al. demonstrated a real-time re-focusing system based on feedback and wavefront control [14]. They used in-line interference to extract the TM by dividing each input field into a reference part and a controlled part. This approach does not need any knowledge of fiber deformation and can obtain the TM in tens of milliseconds with the help of a FPGA. But the instability problem due to interference still exists, as in [2,3].
To date, the phase retrieval technique has been applied for scattering materials [15] and multicore fibers [16]. Image transmission through an optical fiber has also been designed based on phase restoration [17]. It has also been demonstrated that one row of the TM can be retrieved in a reference-less setup with a convex optimization algorithm [18]. In this paper, we demonstrate the acquisition of the complete TM of a MMF using a DMD and a phase retrieval algorithm [15] without a reference beam, thereby leading to a simplified and robust optical setup. A sufficiently large number of output fields can be measured quickly thanks to the high speed of DMD. The algorithm itself runs in a parallel format such that distributed computations can be used to accelerate the data acquisition. The calculated TM is then verified using the following three approaches. Firstly, predicted outputs for random inputs are compared with corresponding actual outputs for prediction accuracy estimation. Secondly, imaging in the transmission mode is demonstrated with pre-defined patterns. Lastly, a set of foci is generated on the distal fiber end with the information in TM as a focus-scanning process. Using the TM, we demonstrate that one can obtain accurate mode count and eigenchannels for optimized energy delivery by performing singular value decomposition on the TM.

Experimental setup
The diagram of our experimental setup is shown in Fig. 1. A DMD with 1024 × 768 micromirrors (Discovery 1100, Texas Instruments) driven by an interface board (ALP-1, ViALUX) was used to modulate the 632.8-nm laser beam (25-LHP-925-230, Melles Griot), which was direction-adjusted by two mirrors (M1 and M2) and expanded by two lenses (L1 and L2). The DMD modulated light field was coupled into the 50-μm-core MMF (1.5-m long and bent into a circle of diameter ~200-mm, FG050UGA, Thorlabs) with a tube lens TL1 (AC254-200-A-ML, Thorlabs) and a microscope objective OL1 (Nikon CFI Plan Achro 40X, 0.65 NA, 0.56 mm WD, tube lens focal length: 200 mm). Corresponding MMF outputs, magnified by another microscope objective OL2 (Olympus 20X Plan Achro, 0.4 NA, 1.2 mm WD, tube lens focal length: 180 mm) and the tube lens TL2 (AC254-100-A-ML, Thorlabs), were passed through a polarizer (LPVISE100-A, Thorlabs) and recorded by the camera (C1140-22CU, Hamamatsu). In our setup, we set 2 × 2 micromirrors as a single macro pixel in the modulation to enhance the difference between the "ON" and the "OFF" pixels. We also used input patterns consisting of N = 1296 macro pixels and captured output patterns with M = 9216 pixels on camera. The number N is limited by the maximum number of modes supported in the MMF, as sufficient degrees of freedom are required to describe the change of input patterns. Also, if there are too many macro pixels, the beam diffracted by the DMD will be too wide to be collected and resized by the tube and objective lenses. The size of the fiber core and magnification of the objective lens constrains the input pattern size and in addition, the computation time will increase with the size of the input and output patterns. The number M is determined by the parameters of the MMF used in our experiment particularly the core diameter and by the magnification of the objective lens. In order to retrieve phase information successfully using this algorithm, M/N ≥ 2 is needed.
An intuitive method to generate input fields for calibration is to turn the different macro pixels "On" sequentially and record the generated output, however, the diffracted light is very weak and results in poor signal to noise ratio (SNR). Instead, we generated P = 8000 random binary patterned input fields, each with half of the pixels "ON" and half "OFF". The response of a single macro pixel was therefore measured P/2 times over the cycle, improving the SNR / 2 P times compared to the single macro-pixel approach. Although the DMD speed can be up to 8 kHz, it was operated at a maximum rate of 250 Hz due to the limit of camera. In addition to the camera restriction, the speed is further limited by the data rate in the DMD USB cable if the modulation data is bigger than the DMD on-board memory, which in our case only allowed 20 successive images to be displayed. Therefore most of the calibration time of around 25 minutes was used in loading data into DMD memory. With the latest DMD models such as Discovery 4100 from Texas Instruments, the calibration should be possible in one minute. Also note that system stability is critical and this was monitored and shown to be 99% stable during the whole experiment by calculating correlation coefficients of outputs for the same input over time. The input and output data was processed with the prVBEM (phase retrieval Variational Bayes Expectation Maximization) algorithm which was used to calculate the TM for scattering material [15].
The reconstruction of a complex field with only the magnitude information is the socalled phase retrieval problem. With an observation y and a M × N complex-measurement matrix D, the target is to recover the signal x, such that Translating problem (1) into a Bayesian framework, the missing phase information of observations and measurement noise are introduced as new variables. Therefore, each complex value in y is given as where u ∈ {1…M}, θ u ∈ [0,2π] denotes missing conjugate phase of y u and n u stands for a zero-mean Gaussian noise with standard deviation σ n . Moreover, it is also supposed that Vol. 26, No. 14 | 9 Jul 2018 | OPTICS EXPRESS 18439 where CN is a centered isotropic normal distribution. With these assumptions, the missing phases in observations are taken into account in the model as marginalizing on θ u results in a distribution of y u , up to the moduli of y u and Based on the previous model, the retrieval of the complex signal x can be given as the solution of the marginalized Maximum A Posterior (MAP) estimation problem as follows: Due to the marginalization on the uncovered variable θ, it is intractable to calculate the p(x|y) directly. However, it is possible to use sub-optimal methods to solve Eq. (7). Here, the solution by the mean-field approximation is introduced. The mean-field approaches are to approximate the posterior joint distribution ( , | ) p x y θ by a distribution ˆ( , ) and can be addressed efficiently by VB-EM (Variational Bayes Expectation Maximization) algorithm. The VB-EM algorithm works in an iterative manner and successively updates the factors of the MF approximation. Considering the model and the MF approximation described previously, the following update Equation can be obtained where CN denotes circularly symmetric and complex normal distribution and , . In the previous Equations, I 0 (resp. I 1 ) denotes order-0 (resp. 1) the modified Bessel function of the first kind, (.) H stands for the conjugate transpose and (.)* denotes scalar conjugate. In terms of problem Eq. (8), the approximation of p(x|y) can be given as leading to the estimation of row x and ˆi x given by: In short, this algorithm jointly infers one row of the output field phase and one row of the TM iteratively in a Bayesian way. The iteration is stopped if the change to both is below the threshold (i.e., in a stable state), and the steady values are returned as the best estimates. This algorithm is also suitable for parallel computation since estimations of different TM rows are independent. In our implementation, high performance computing facilities (a group of Intel Xeon E5645 and AMD Opteron 6234) were used to perform parallel calculations, requiring 1 to 3 minutes to calculate one TM row, depending on the algorithm convergence time. The acceleration comes from parallelized jobs running concurrently on different nodes, although this approach was not extended to be fully parallel -with a theoretical calculation time of a few minutes -due to the practical limitation of jobs submitted by other users. Typically, the number of parallel jobs (number of CPU cores acquired) was between 100 and 200 requiring approximately 4 hours to calculate the TM, a large improvement over the serial case (calculated at 300 hours).

TM evaluation
In order to test the calculated TM, we determined the correlation coefficient of the TMpredicted outputs with the recorded outputs for the same inputs using the corr2/corr functions in MATLAB and averaged over around 8000 input patterns. The relationship between the quantity of data and the accuracy of the calculated TM is shown in Fig. 2, where P denotes the number of input patterns and N is the macro pixel count on DMD. It is found that more input fields result in stronger prediction ability. From a Bayesian inference perspective, the posterior probability distribution is updated and narrowed down by feeding the algorithm more and more data. Since the expectation value of the latest posterior probability is returned as the inferred value in the algorithm, the estimated value tends to deviate more from the true value if there is insufficient data. As the amount of data increases, the posterior is sharp enough for the estimated value to be very close to the true value. For our setup, P/N = 6 is sufficient for a successful implementation as the corresponding prediction accuracy is around 98%, demonstrating high-quality TM estimation and prediction as precise as for holography based methods. In the following sections, P/N = 6 was therefore used to ensure good performance with calculations accelerated in a parallel way.
With TM, we were able to reconstruct the displayed DMD patterns based on the recorded camera images, which is a demonstration of image transmission through the MMF. To illustrate this, three pre-defined patterns ('LIV', '3.14' and a smiley face) were intentionally included in the input set and the phases of corresponding outputs were recovered with the algorithm. The input DMD pattern recovery performed here is to verify the calculated TM and the system resolution. The calculated TM can be then used to determine the input DMD patterns for focused spots raster-scanning at the distal end [9]. The inferred DMD patterns E OP were calculated by solving the following linear Equation: , where E OP and E IP stand for the electric fields on the proximal and distal ends of the fiber respectively. During the calibration, the DMD was functioning as an object and so the proximal end was defined as OP while the distal end was defined as IP. The E OP may be calculated by using the inverse of the TM, however, as the TM is not a square matrix (N ≠ M) in most cases, and has no strict inverse. The pinv function in MATLAB could be used to calculate the so-called pseudoinverse as a substitute of the strict inverse, although considering that for our TM N is a few times larger than M, such replacement may cause additional noise. Solving the Eq. (18) directly can avert this problem and is faster, especially for large TM. This step was done in MATLAB and the effectiveness of the TM was checked afterwards.
Due to the noise, the TM is imperfect but can still achieve 98% correlation with 25% less calculation time than the pseudoinverse method. The inferred patterns are compared to the actual input patterns in Fig. 3. The patterns could be reproduced accurately through the MMF, demonstrating that the TM was estimated correctly and capable of undoing the distortion caused by MMF. The line widths in the three patterns were designed to be around 1.4 μm to evaluate the system imaging resolution. Since the adjacent lines, in Figs. 3(d) and 3(e), and small dots in Figs. 3(e) and 3(f), were well separated from the background and imaged clearly, the resolution is close to the Abbe diffraction limit calculated as λ/(2NA) = 1.44 μm. Therefore, our MMF is being fully exploited and working at or near diffraction-limited resolution. To further improve imaging performance, higher-NA MMF may be utilized or a shorter-wavelength laser may be used. Note that although DMD patterns were used here as objects, for general imaging the object or camera should be matched with the modulation device [7,8] as the TM depends on the modulator's specific position. For the best performance, the optical distance between the modulator and the proximal fiber end should be equal to that between object/camera and the proximal fiber end. Additionally, when matching the camera with the modulator, it is also suggested that they have same pixel size and pitch and that transverse displacement is minimized.
In order to further evaluate our system, we utilized the information in the calculated TM to perform focus generation through the MMF. Specifically, the complex conjugate of one TM row was chosen as a necessary input field for generating a focus at the distal fiber end. The position of the focus was controlled by choosing the corresponding row in the TM. To implement this complex input wave with the DMD, we selectively turned on the macro pixels where the real part was positive. By doing so, we were able to create a clean focus through the MMF, such as the one around the fiber end center shown in Figs. 4(c) and 4(d) by fitting the data with a Gaussian function. The result shows that the full width at half maximum (FWHM) of the focus in Fig. 4(a) is about 2.2 μm along the y-axis and 1.9 μm along the xaxis. We attribute the difference to the fact that the MMF is not an absolutely axisymmetric waveguide due to random fiber bending. The theoretical value of FWHM in our case is around 1.44 μm (Abbe Criterion) and the degradation results from the difference between ideal input field and practical input field. Note that input fields are processed to be binary as DMD amplitude modulation is not continuous. Such field transformation will inevitably cause that input to differ from the ideal case.
By choosing different TM rows, we scanned the focus across a 50 by 50 grid within the fiber end and recorded the intensity enhancement factors (EFs), defined as the ratio of the focus peak value to the average value of output generated by a random input. The EF distribution is presented as a heatmap in Fig. 4(b). It is found that the maximum EF = 70 was obtained due to fully exploiting the spatial modes in the MMF and accurate phase retrieval. Further improvement could be obtained by using more DMD macro pixels and higher-NA fibers [8]. It was also observed that the focus quality is not uniform over the fiber end. The DMD modulation is binary, which creates an error between the actual input fields and the ideal ones and this error changes for the different positions creating a variation in the focus quality. To get a higher EF or a more balanced EF distribution, one may want to utilize more macro pixels on the DMD or to use an LC-SLM instead (at the expense of speed) as both approaches will bring stronger modulation and more precise control of the input field.

TM properties
To the best of our knowledge, TMs have mainly been used to generate specific patterns or to realize imaging. In this section, we analyze other TM information mathematically and extract information on fiber mode count and eigenchannels for optimal energy delivery.
In most practical applications multimode fibers are not straight and deformation-free, and the bending or twisting impacts mode properties such as propagation constant, transmittance as well as mode count, leading to potential problems in information delivery. It is much harder to confirm the mode properties of a practical fiber with deformation compared to an ideal multimode fiber, for which one can easily estimate the mode count N m by 2 0.5 , where F stands for the normalised frequency of the multimode fiber. Equation (1) is obtained by solving the Helmholtz function for step-index fibers and it is assumed that F is large (e.g. 100). Further, the normalized frequency can be calculated by where a denotes fiber core radius, NA stands for its numerical aperture and λ is wavelength. When dealing with deformed fibers, it is impractical to solve the wave Equations and find the mode properties due to the asymmetry, complex boundary conditions and unknown deformation. In the previous section, we demonstrated that the TM of a multimode fiber imaging system could be found in a reference-less setup and can describe how the light field is transformed by the whole system by relating independent modes at the input plane with those at the output plane. Here, we further demonstrate that it is possible to find the actual mode count in a multimode fiber by calculating the rank of its TM. The elements in the TM are the coefficients of linear Equations relating input and output pixels. Therefore the rank of the matrix indicates the number of independent Equations, which is also the number of degrees of freedom in the imaging system. Note that the number of the spatial modes is the count of degrees of freedom of the system. Therefore we can deduce that the rank of the recorded TM [19] is actually the number of modes travelling in the multimode fiber, which can be obtained regardless of the specific deformations and as an alternative to solving the complex wave Equations.
To record the TM of a specific multimode fiber, we used the experimental setup introduced section 2 and Fig. 1. The prVBEM algorithm was utilized to infer the phase information for output images and to calculate the elements of the TM. The singular values of the TM were obtained by singular value decomposition: Here, U and V are unitary matrices while S is a diagonal matrix of non-negative real numbers, which are the singular values of the TM. In theory, the rank of the TM is the count of the positive singular values. Figure 5 shows the squared singular value distribution of the TM of a multimode fiber. It is clear that the singular value drops with increasing mode index. Since the squared singular value indicates mode transmittance, our observation is in line with the fact that higher-order modes experience larger loss, as also found in [19]. Also note that all the singular values are positive in this TM although the values are quite small in the tail, indicating the rank of TM is 1296. However, considering system noise, which can be estimated by calculating the correlation between recorded output patterns and predicted output patterns, the effective rank of the matrix is defined by the specific mode index, N R , such that the area below the curve up to N R is 98% of the total area. This produces a mode count of 840, significantly different from the figure of 755, which is obtained from half of the estimation of from Eq. (19). The reason for halving the estimation is that a polarizer was used in our experimental system and therefore the fiber modes in one polarization direction were filtered out. The difference is due to two reasons: 1) Eq. (2) is more accurate for large normalized frequency, and 2) Eq. (19) is for an ideal straight fiber without bends or twists. We also adjusted an additional translation stage (MT1/M, Thorlabs) attached to the middle of the fiber and observed the change of TM rank. As shown in Fig. 6, fewer fiber modes were supported in the fiber for larger translations, i.e. for more severe bending. Theoretically, energy is more likely to couple into cladding modes or radiation modes in the case of increased bending, explaining the decrease in fiber mode count. The energy coupling efficiency depends on the input field distribution and the fiber mode distribution: higher overlap between the input field and the mode field will result in better coupling efficiency. Also, different modes experience varying transmission loss and this is also related to the specific fiber geometry. Choi's group studied eigenchannels in scattering media to maximize energy delivery efficiency [20], but to the best of our knowledge, similar work has not been done for multimode fibers. Here, we further demonstrate that the column vectors in matrix V are actually eigenchannels in the TM, the first one of which corresponds to the input pattern with maximum energy delivery efficiency.
With a TM, we can first theoretically study the energy transport for each eigenchannel. In order to do this, we virtually injected column vectors of matrix V into the fiber and obtained corresponding outputs with the help of the TM. The normalized energy transmittance is shown in Fig. 5(b) and one can clearly see the descending trend, indicating the first column vector is optimal for energy transport. In addition, the normalized energy transmittance values are exactly the normalized squared singular values, indicating the physical meaning of the singular values of the TM. The results are similar to those in [20] where the decrease ofenergy efficiency with the increase of eigenchannel index is clear. Note that we measured output images for random input patterns in order to construct the TM. Here, we used these input and output data to calculate the energy transport efficiency for random input patterns. In addition, we implemented the best eigenchannel and recorded the corresponding output. The energy delivery efficiency for random patterns was normalized to that of best eigenchannel and the results are presented in Fig. 7. In this experiment, we observe that the energy transport efficiency for the eigenchannel is over 1.5 times higher than that of random patterns. However, this value is less than the theoretical enhancement of around 8 times in our case. The primary reason for this observation is that it is only possible to realize binary amplitude modulation with the DMD.  The DMD-limited and perfect eigenchannels are compared in Fig. 8 where their corresponding output patterns are also presented. Due to the difference between the realistic pattern and the ideal pattern, the outputs are distinct, leading to lower energy delivery efficiency. This drawback significantly limits the ability to generate accurate patterns on the fiber proximal end. Nevertheless, the improvement on the energy transport efficiency with imperfect input patterns demonstrates the existence of eigenchannels. A great improvement may be predicted if using a phase modulator such as LC-SLM or by using a DMD with a lens to offer phase modulation as in [10]. In our experiment, the fiber proximal end is aligned with OL1 in Fig. 1 by maximizing the fiber output power with all DMD macropixels turned on. Since the transmission matrix describes the optical transmission from the DMD to the camera, the transmission properties (e.g. eigenchannels) are affected by both the distance and the angle between OL1 and the fiber proximal end.

Conclusions
In summary, it has been demonstrated that the full TM of a multimode fiber is achievable in a reference-less setup and with phase retrieval techniques. Such a configuration leads to not only a simplified system but also better system stability. The TM is calculated by the phase retrieval algorithm, which can be accelerated by distributed or parallel calculation. The generated TM has also been examined to have a strong predictive ability, to be able to perform diffraction-limited imaging and to be capable of concentrating light to a focus. Furthermore, it has been demonstrated that information on mode count and eigenchannels is available from the TM by singular value decomposition. We believe that the proposed setup is an advance over current approaches and opens a new method for MMF characterization. As a fast and robust implementation, it is a promising candidate for MMF focusing and endoscopic imaging.