Intensity-only-measurement mode decomposition in few-mode fibers

Egor Manuylovich; Egor Manuylovich; Aleksandr Donodin; Sergei Turitsyn; Sergei Turitsyn

doi:10.1364/OE.437907

1. Introduction

Recent renewed interest in multi-mode fibers (MMFs) and, in particular, in few-mode fibers (FMFs) is dominated by their potential in telecommunications, however, their features are equally important for a wide range of applications. In contrast to conventional single-mode fiber (SMF), FMFs and MMFs provide larger mode areas allowing nonlinearity suppression and damage threshold improvement for high-power laser systems and amplifiers [1,2]. Such fibers open up new ways of imaging [3], including starlight [4], sensing [5], and fundamental studies of novel nonlinear dynamics [6] including optical beam self-cleaning [7], spatiotemporal solitons [8] and ultrafast characterization methods [9,10]. Moreover, an unrelenting capacity increase of the contemporary telecommunication systems dictate using novel tools and principles of signal transmission to seize this demand [11]. Spatial mode-division multiplexing is one of the popular approaches to improve the capacity of optical links especially by using FMFs that can significantly increase communication capacity compared to SMF without a drastic complexity of signal processing [12]. However, all of these applications require beam characterization not only in terms of beam divergence ($M^2$ parameter), but also in terms of waveguide eigenmodes amplitudes and phases. This characterization is also known as mode decomposition (MD) problem. Recovery of a signal from the measured intensity is an attractive approach with various applications ranging from imaging to telecommunications. However, direct detection-based MD methods are usually computationally-intensive and/or capable of recovering few modes only.

Some of the existing MD methods require a reference beam. Digital holography [13,14] and multi-lane light conversion [15] methods are both in this method group. The main drawback of these methods is the requirement of a coherent laser source in the receiver part of the system that limits potential applicability. Those methods that do not require a reference beam can be divided into iterative and non-iterative algorithms. Iterative methods like Gerchberg–Saxton technique [16], line-search [17], and stochastic parallel gradient descent [18] show high accuracy and speed. However, these methods are vitally sensitive to initial conditions and can be looped at a local minimum.

On the other hand, non-iterative MD methods include machine learning methods [19–24] and fractional Fourier system [25]. Machine learning MD algorithms outperform iterative methods in the decomposition speed with a price of high-performance computer requirements, long training times, and poor performance in cases with more than five spatial modes. Recently, a novel non-iterative method, based on dividing inherently non-linear MD problem into a cumbersome linear part and a simple non-linear part, was proposed [26]. This algorithm was based on near-field (NF) MD and allowed to decompose up to 100,000 multi-mode distributions per second for three-, five-, and eight-mode fibers; along with maximum 27-mode decomposition in noiseless case.

Here we propose a significant improvement of the algorithm presented in [26] in terms of the decomposition accuracy by using far-field (FF) and near-far-field (NFF) MD techniques. In addition to the accuracy boost, both techniques allowed a boost in the maximum number of modes decomposed to up to 46 modes for FF method and up to 49 modes for NFF method. In addition, all three possible MD technique variants based on the NF, FF, and NFF decomposition techniques are compared in terms of their accuracy, number of decomposed modes and computational complexity. Moreover, we experimentally verify proposed FF method using a 5-mode fiber and a 21 dB signal-to-noise ratio. We expect the developed technique to find many applications from ultra-fast fiber imaging to laser beam characterization.

2. Methods

To solve the MD problem, it is required to determine the weighting factors $C_k = A_kexp(i\phi _k)$ at which a given intensity distribution is observed at the fiber end facet (near-field) or at the infinite distance from the fiber end facet (far-field). The inference of amplitudes and phases based on the intensity distribution is a non-linear problem that can be solved with different approaches. In the previous article, we presented a fast MD method in few-mode fibers, which is based on dividing of the problem into a cumbersome linear part and a simple non-linear part [26]. The presented method makes it possible to determine the coefficients based on the intensity distribution in the near field:

(1)$$P_{NF}(x,y) =\left\langle\left|E(x,y)\right|^2\right\rangle =\left\langle\left|\sum_{k}C_k\Psi_k(x,y)\right|^2\right\rangle,$$

where $P_{NF}(x,y)$ is the intensity distribution in the near field, and $\Psi _k$ is the near field distribution of $k^{th}$ eigenmode of the fiber.

However, the presented method makes it possible to perform MD only up to the complex conjugation of the weight coefficients, since it is based on the analysis of the intensity distribution in the near field. The near-field (NF) technique requires the use of a rather bulky 4-f system (Fig. 1(a)) to magnify and project the NF distribution onto the image sensor. In this paper, we present a development of the method that allows either the use of a much more compact experimental FF image capture system (Fig. 1(b)) or the use of both near and far-field (NFF) methods (Fig. 1(c)) to determine the true values of the mode coefficients, rather than only the complex conjugations.

Fig. 1. Experimental setups for capturing a) near-field intensity distribution, b) far-field distribution and c) both near and far-field intensity distributions.

Download Full Size | PDF

Far-field measurements are more convenient both because of the compactness of the optical image capturing system and the possibility to easily scale the image resolution. In the case of NF intensity distribution capturing, an increase in the resolution requires an increase in the magnification of the 4-f lens system, which leads to the proportional scaling of the experimental scheme. Thus, to obtain a magnification factor of 67 for the NF capturing scheme, we used a 4-f system with focal lengths of 4.5 mm and 300 mm, which led to an experimental setup overall length of more than 600 mm. To obtain a higher magnification and, accordingly, a larger image size, it is necessary to further increase the focal length of a long-focus lens, which leads to an almost inapplicable length of the experimental setup. In the case of the FF capturing, the magnification coefficient of 67 can be obtained with the length of the free-space beam path of 5 mm, and to obtain a magnification factor of 200, only a length of 15 mm is required. To solve the problem of the FF MD, it is required, as in [26], to calculate the matrix of pairwise products of the fiber eigenmodes. First of all, we calculate the FF distribution for each of the fiber modes and their pairwise products. The FF distribution $\tilde{\Phi }(\xi,\eta )$ can be calculated as the 2D Fourier transform of the NF distribution:

(2)$$\tilde{\Phi}(\xi,\eta)=\mathcal{F}(\Psi(x,y))=\frac{1}{2\pi} \iint\Psi(x,y)\exp({-}i\xi x-i\eta y)dxdy,$$

We can rewrite the Eq. (1) for the FF distribution:

(3)$$P(x,y) =\left\langle\left|\sum_{k}C_k\tilde{\Phi}_k(x,y)\right|^2\right\rangle,$$

where $P(x,y)$ is the FF intensity distribution.

In FF case, either a purely real distribution (for centrally symmetric modes) or a purely imaginary (for centrally antisymmetric modes) is observed in the FF. When calculating the matrix of pairwise products of eigenmodes, we "turn" the phase of the purely imaginary FF modes by an angle of $\pi /2$ to bring them to purely real distributions:

(4)$$\begin{cases} \Phi(\xi,\eta)=\tilde{\Phi}(\xi,\eta),\hspace{62pt}if~~ \tilde{\Phi}(\xi,\eta)=\tilde{\Phi}(-\xi,-\eta), \\ \Phi(\xi,\eta)=\tilde{\Phi}(\xi,\eta)\cdot \exp(i\pi/2),\hspace{15pt}if~~ \tilde{\Phi}(\xi,\eta)={-}\tilde{\Phi}(-\xi,-\eta). \\ \end{cases}$$

In what follows, we will take this rotation of $\pi /2$ into account when calculating the phases of the eigenmodes. We also need to calculate the matrix of pairwise products of eigenmodes for the FF distribution $\mathbf {T}_{FF}$:

(5)$$\tiny \mathbf{T}_{FF}=\begin{pmatrix} \Phi_1^{(1)}\Phi_1^{(1)} \cdots \Phi_N^{(1)}\Phi_N^{(1)} &2\Phi_1^{(1)}\Phi_2^{(1)} \cdots 2\Phi_1^{(1)}\Phi_N^{(1)} &2\Phi_2^{(1)}\Phi_3^{(1)} \cdots 2\Phi_2^{(1)}\Phi_N^{(1)} \cdots 2\Phi_{N-1}^{(1)}\Phi_N^{(1)} \\ \;\;\;\vdots\qquad\vdots\;\;\; &\;\;\;\vdots\qquad\vdots\;\;\; &\;\;\;\vdots\qquad\vdots\qquad\vdots\;\;\; \\ \Phi_1^{(m)}\Phi_1^{(m)}\cdots\Phi_N^{(m)}\Phi_N^{(m)} &2\Phi_1^{(m)}\Phi_2^{(m)} \cdots 2\Phi_1^{(m)}\Phi_N^{(m)} &2\Phi_2^{(m)}\Phi_3^{(m)} \cdots 2\Phi_2^{(m)}\Phi_N^{(m)} \cdots 2\Phi_{N-1}^{(1)}\Phi_N^{(1)} \\ \;\;\;\vdots\qquad\vdots\;\;\; &\;\;\;\vdots\qquad\vdots\;\;\; &\;\;\;\vdots\qquad\vdots\qquad\vdots\;\;\; \\ \Phi_1^{(M^2)}\Phi_1^{(M^2)}\cdots\Phi_N^{(M^2)}\Phi_N^{(M^2)} &2\Phi_1^{(M^2)}\Phi_2^{(M^2)} \cdots 2\Phi_1^{(M^2)}\Phi_N^{(M^2)} &2\Phi_2^{(M^2)}\Phi_3^{(M^2)} \cdots 2\Phi_2^{(M^2)}\Phi_N^{(M^2)} \cdots 2\Phi_{N-1}^{(M^2)}\Phi_N^{(M^2)} \\ \end{pmatrix}$$

Here the upper index ($m$) corresponds to the $m^{th}$ pixel of the image of eigenmode. Figure 2 shows the distributions of the $\mathbf {T}_{FF}$ matrix components for 5-mode fiber. Each image corresponds to one column of matrix $\mathbf {T}_{FF}$. Each image is a point-wise product of the corresponding eigenmodes, and the column of the $\mathbf {T}_{FF}$ matrix is obtained by rearranging these $M$ by $M$ distributions into $M^2$ by 1 vector. The same rearranging is then applied to the intensity distribution to obtain vector $P$.

Fig. 2. Components of matrix $\mathbf {T}_{FF}$, representing pairwise products of far-field of eigenmodes of 5-mode fiber. The upper indices $^s$ and $^c$ stand for two different orthogonal degenerate sub-modes of non-centrally-symmetrical modes.

Download Full Size | PDF

The dependency of matrix $\mathbf {T}_{FF}$ condition number on the number of fiber eigenmodes is shown in Fig. 3.

Fig. 3. Condition number of matrix $\mathbf {T}_{FF}$ calculated for a various number of modes with M=128.

Download Full Size | PDF

In comparison to the $\mathbf {T}$ matrix for the near field distribution (see [26]) there is no sharp increment in the region around 10 modes. This is due to the fact that unlike the near-field matrix $\mathbf {T}$ the columns of the $\mathbf {T}_{FF}$ matrix which are composed of $LP_{11}$ and $LP_{12}$ modes do not match exactly (see Fig. 4). Such a lack of symmetry positively affects the problem numerical solution, as the $\mathbf {T}_{FF}$ is well-conditioned. The phase recovery method is generally similar to the technique described in [26]. We first rewrite the Eq. (3) in the linear form using Eq. (5) and solve the linear problem:

(6)$$\mathbf{T}_{FF}\cdot z=P,$$

Fig. 4. Non-coinciding components of matrix $\mathbf {T}_{FF}$, related to modes $LP_{11}$ and $LP_{12}$ and their discrepancy.

Download Full Size | PDF

Here $P$ is the intensity distribution in the FF rearranged into a $M^2$ by 1 vector, $\mathbf {T}_{FF}$ is the matrix of pairwise products of the far-field of eigenmodes.

For combinations of eigenmodes including either only centrally symmetric or only centrally antisymmetric modes, we obtain:

(7)$$z_n= \frac{C_kC^*_j+C_jC^*_k}{2}\hspace{20pt}k,j=1..N,\hspace{20pt}n=1..N(N+1)/2.$$

or $z_n=A_kA_jcos(\phi _k-\phi _j)$.

For combinations of eigenmodes including one centrally antisymmetric mode, we get:

(8)$$z_n= \frac{C_kC^*_j-C_jC^*_k}{2i}\hspace{20pt}k,j=1..N,\hspace{20pt}n=1..N(N+1)/2.$$

or $z_n=A_kA_jsin(\phi _k-\phi _j)$, since we have multiplied the central antisymmetric modes by $exp(i\pi /2)$. Including the above mentioned relations, one can represent the vector $z$ as a matrix $\mathbf {Z}$. To do this, we numerate the vector $z$ along the main diagonal first, and then along the columns of the lower triangular matrix $\mathbf {Z}$:

(9)$$\mathbf{Z}=\begin{pmatrix} z_1 & 0 & \cdots & &0 \\ z_{N+1} & z_2 & 0 & &\\ z_{N+2} & z_{2N} & z_3 & \ddots & \vdots\\ \vdots & z_{2N+1} & \vdots & \ddots & 0\\ z_{2N-1} & \cdots & &\cdots &z_N\\ \end{pmatrix}=\begin{pmatrix} Z_{1,1} & 0 & \cdots & &0 \\ Z_{2,1} & Z_{2,2} & 0 & &\\ Z_{3,1} & Z_{3,2} & Z_{3,3} & \ddots & \vdots\\ \vdots & \vdots & \vdots & \ddots & 0\\ Z_{N,1} & Z_{N,2} & &\cdots &Z_{N,N}\\ \end{pmatrix}.$$

And, for example, for 5-mode fiber we can obtain the matrix $\mathbf {Z}$ as following:

(10)$$\mathbf{Z}=\begin{pmatrix} A_1^2 & 0 & 0 & 0 & 0 \\ A_1A_2sin(\phi_2) & A_2^2 & 0 & 0 & 0\\ A_1A_3sin(\phi_3) & A_2A_3cos(\phi_2-\phi_3) & A_3^2 & 0 & 0\\ A_1A_4cos(\phi_4) & A_2A_4sin(\phi_2-\phi_4) & A_3A_4sin(\phi_3-\phi_4) & A_4^2 & 0\\ A_1A_5cos(\phi_5) & A_2A_5sin(\phi_2-\phi_5) & A_3A_5sin(\phi_3-\phi_5) & A_4A_5cos(\phi_4-\phi_5) & A_5^2\\ \end{pmatrix}.$$

We can also obtain the general relations between the amplitudes, phases and the matrix eigenvalues as following:

(11)$$\begin{aligned}& A_k = \sqrt{Z_{k,k}},\hspace{28pt} k=1..N, \\ &cos(\phi_k)=\frac{Z_{k,1}}{A_1A_k},\hspace{10pt} if~ \tilde{\Phi_1}(\xi,\eta)\tilde{\Phi_k}(\xi,\eta)=\tilde{\Phi_1}(-\xi,-\eta)\tilde{\Phi_k}(-\xi,-\eta), \\ &sin(\phi_k)=\frac{Z_{k,1}}{A_1A_k},\hspace{12pt} if~ \tilde{\Phi_1}(\xi,\eta)\tilde{\Phi_k}(\xi,\eta)={-}\tilde{\Phi_1}(-\xi,-\eta)\tilde{\Phi_k}(-\xi,-\eta).\end{aligned}$$

To determine the sign of the phase when taking the arcsine or arccosine, we use the second column of the matrix $\mathbf {Z}$ . It should be noted that we use only the first $3N-3$ components of the vector $z$ for calculating mode weights. For a number of modes $N>3$, it means that we do not need to calculate the product of full matrix $\mathbf {T}^{-1}$ with the intensity vector $P$. We need to calculate only the product of the first $3N-3$ rows of this matrix with the intensity vector. Thus, based on the FF intensity distribution, we calculate the vector $z$, using the Eq. (6) by multiplying the pseudo inverse matrix of pairwise products of eigenmodes by the vector of intensity. Then using this vector $z$ we calculate the amplitudes and phases of the modes using Eqs. (9) and (11).

The analysis of only the NF or only the FF does not allow one to unambiguously reconstruct the phases of the modes, since there is phase transformation that leaves the intensity distribution unchanged. In the case of the near field, this is a complex conjugation of the weights (changing the signs of all phases to opposite ones): $\phi _k^{new}=-\phi _k$ , and for the FF, this is the replacement of the phase sign with the opposite one for centrally symmetric modes (for example, $LP_{21}$ or $LP_{02}$): $\phi _k^{new}=-\phi _k$ and transformation $\phi _k^{new}=\pi -\phi _k$ for non-centrally antisymmetric (for example $LP_{11}$ or $LP_{31}$):

(12)$$\begin{cases} \phi_k^{new}={-}\phi_k,\hspace{24pt}if~~ \tilde{\Phi}(\xi,\eta)=\tilde{\Phi}(-\xi,-\eta), \\ \phi_k^{new}=\pi-\phi_k,\hspace{15pt}if~~ \tilde{\Phi}(\xi,\eta)={-}\tilde{\Phi}(-\xi,-\eta). \\ \end{cases}$$

The phase transformations, which conserve the intensity distribution, are different for the NF and FF. This means that knowing both NF and FF distributions allow to calculate the phases exactly, and not only up to complex conjugation (as in the case of the near field) or transformation (12), as in the case of the FF. Therefore, if both FF and NF intensity distributions are measured, then one can decompose mode weights up to a constant phase shift, not only to an intensity-conserving transformation. Figure 5 shows how NF and FF intensity distributions look like for different sets of weights that show either the same NF intensity distribution or the same FF intensity distribution.

Fig. 5. NF and FF intensity distributions before and after weight transformation that conserves either NF or FF intensity distribution.

Download Full Size | PDF

A complete MD can be approached either by calculating the cosine value of the second mode by using NF decomposition and then applying FF MD or it can be done by combining NF and FF matrices in one system of equations. The first approach is just a simple extension of the presented FF method and has the same accuracy as the FF MD method. The second approach requires complex-value calculations which increases the complexity of the method, leading to higher time consumption and accuracy. To do this we simply combine both NF and FF systems of equations and just make a change of variables as simple as:

(13)$$z_n=C_kC^*_j\hspace{20pt}k,j=1..N,\hspace{20pt}n=1..N^2.$$

In this case, we cannot combine $C_kC_j^*$ and $C_k^*C_j$ as there is no symmetry in matrix $\mathbf {T}$ anymore. In the NFF case we use complex-valued pseudo inverse $\mathbf {T^{-1}}$ matrix to calculate the complex-valued z vector. Mode amplitudes and phases are then calculated similarly to (11), but cosine values are proportional to real parts of the corresponding values of matrix $\mathbf {Z}$ and sine values are proportional to imaginary parts of the corresponding values of matrix $\mathbf {Z}$.

3. Results

The MD based on the FF analysis shows better results in the noiseless case compared to the NF case in terms of the decomposition accuracy. In addition, the method works up to 46 modes in the noiseless decomposition (Fig. 6) leading to a significant increase in comparison to 27 mode decomposition for the NF case [26].

Fig. 6. Noiseless MD for 46-mode fiber. True FF intensity (a), reconstructed FF intensity (b), their discrepancy (c), and true and reconstructed weights (d).

Download Full Size | PDF

For 48-mode fiber noiseless decomposition is performed with some errors, but the most of the weights are still determined with a high accuracy (Fig. 7). So, in principle, the proposed method can be applied for the FF MD in up to 46-mode fiber. The effect of the noise presence on the accuracy of the proposed MD algorithm has also been studied using following error metrics:

(14)$$\epsilon_A=\frac{||A^{recov}-A^{true}||}{||A^{true}||},\hspace{10pt} \epsilon_{\phi}=\frac{max(\phi_i^{recov}-\phi_i^{true})}{2\pi},\hspace{10pt} \epsilon_C=\frac{||C^{recov}-C^{true}||}{||C^{true}||}.$$

Fig. 7. Noiseless MD in 48-mode fiber. True FF intensity (a), reconstructed FF intensity (b), their discrepancy (c), and true and reconstructed weights (d).

Download Full Size | PDF

The amplitude error $\epsilon _A$ represents the normalized amplitude error, the phase error $\epsilon _{\phi }$ shows the maximum phase error across phases of all modes normalized by $2\pi$, and the total error $\epsilon _C$ represents the normalized complex-valued weights error. The latter error includes both amplitude and phase errors. To calculate the accuracy of the presented algorithm we performed multiple decompositions for random sets of mode weights. For FF MD we chose amplitudes from a continuous uniform distribution: $A_{1..N}^{true}~U(0,1)$. To avoid the ambiguity in the MD we chose the phase of the second mode from a continuous uniform distribution $\phi _2^{true}~U(-\pi /2,\pi /2)$ and all other phase coefficients $\phi _{3..N}^{true}~U(0,2\pi )$. Then we calculated the true FF intensity corresponding to each set of weights. We added Gaussian noise, and used this noisy intensity distribution to recover mode weights using our algorithm. Based on the recovered weights we calculated error metrics. Figure 8 shows examples of noisy MD using our algorithm.

Fig. 8. Examples of MD with SNR = 10 dB.

Download Full Size | PDF

We consider a total error $\epsilon _C$ of less than 0.1 as acceptable. So later we refer to this value when we calculate the SNR value necessary to achieve such an accuracy. However, even if $\epsilon _C>0.1$ but less than 0.2-0.3, the recovered weights lay close to their true values. Figure 9 shows how these errors depend on the signal to noise ratio and the number of modes. Each point is averaged over 5000 decompositions. Decomposition is performed by processing 100x100 FF intensity distributions.

Fig. 9. Amplitude, phase, and net errors depending on the noise level in the input image and the number of modes for FF only decomposition. Errors are in dB and color-coded, level lines corresponding to a certain level of error are shown and labeled in the figures. Error metrics are calculated for an image size of 100x100 pixels.

Download Full Size | PDF

The method can be applicable for the MD with an accuracy of 0.1 in 3-mode fiber when SNR is better than 15 dB, in 5-mode fiber for the SNR>22.5 dB and in the 6-mode fiber when the SNR is better than 32 dB. For the 8-mode fiber only amplitudes can be determined with an accuracy of 0.1 when the SNR is better than 25 dB, and for the 10-mode fiber the SNR better than 33 dB is required to determine amplitudes with the same accuracy. This results in a superior performance compared to our previously published MD algorithm based on the decomposition of the NF intensity distribution.

To calculate the accuracy of the NFF method we chose amplitudes from $A_{1..N}^{true}~U(0,1)$ and phases from $\phi _{3..N}^{true}~U(0,2\pi )$ as this method allows decomposing mode weights without any ambiguity. In a noiseless case, the NFF method can perform MD for 49-mode fibers (Fig. 10). This method is also better in the presence of noise. Figure 11 shows how decomposition errors depend on the SNR and the number of modes. This graph was obtained with the same metrics as theFig. 9.

Fig. 10. Noiseless MD for 49-mode fiber. True FF intensity (a), reconstructed FF intensity (b), their discrepancy (c), and true and reconstructed weights (d).

Download Full Size | PDF

Fig. 11. Amplitude, phase, and net errors for NFF decomposition depending on the noise level in the input image and the number of modes. Error metrics are calculated for an image size of 100x100 pixels.

Download Full Size | PDF

This method can be applicable for the MD with an accuracy of 0.1 for 3-mode fiber when the SNR is better than 10 dB, in 5-mode fiber for the SNR>19 dB, and for the 6-mode fiber when the SNR is better than 31 dB. For 8-mode fiber only amplitudes can be determined with an accuracy of 0.1 when the SNR is better than 25 dB, and for the 10-mode fiber the SNR better than 32 dB is required to determine amplitudes with the same accuracy. Thus, the desired accuracy of 0.1 is achieved with 3 dB lower SNR for 3- and 5-mode fibers when using NFF intensity technique than when using only FF intensity distribution.

4. Computational complexity

One of the universal metrics of the computational complexity of algorithms is the number of multiplications used. In this section, we calculate the number of multiplications for our algorithms. The MD method based on the FF analysis consists of the multiplication of pseudo inverse of the matrix of pairwise products of eigenmodes with the vector of intensity. For intensity image of size $M$ by $M$ and the number of modes $N$, one needs to calculate the product of the first $3N-3$ rows of the size of matrix $\mathbf {T}_{FF}^{-1}$ with the intensity vector. So, the linear part of our algorithm requires $Q_{lin}=3M^2(N-1)$ multiplications.

The nonlinear part consists of calculating square roots, arc-cosine, and cosine functions which are not taken into account as these operations can be performed by using lookup tables. The only operation that requires multiplication is

(15)$$\begin{cases} cos(\phi_k)=\frac{Z_{k,1}}{A_1A_k},\hspace{24pt}if~~ \tilde{\Phi}(\xi,\eta)=\tilde{\Phi}(-\xi,-\eta), \\ sin(\phi_k)=\frac{Z_{k,1}}{A_1A_k},\hspace{25pt}if~~ \tilde{\Phi}(\xi,\eta)={-}\tilde{\Phi}(-\xi,-\eta). \end{cases}$$

for each of the N-1 phase coefficients. Here we have 2 multiplications and 2 divisions, so for the nonlinear part of our algorithm we have $Q_{nln}(N,M)=4(N-1)$ multiplications. The total number of multiplications for the FF method is, therefore:

(16)$$Q_{FF}(N,M)=Q_{lin}+Q_{nln}=3M^2(N-1)+4(N-1)=(3M^2+4)(N-1).$$

As $3M^2>>4$ in every practical case, $Q_{FF}\approx 3M^2(N-1)$ a good approximation.

To compare the computational complexity with the best CNN-based MD algorithm presented before [24] we need to calculate the complexity of the VGG-16-like CNN that was used for MD. The adaptation of the VGG-16 convolutional neural network that was used for the MD problem consists of 13 convolutional layers (CL) and 2 fully connected (FC) layers. The number of convolutional filters and sizes of layers is shown in Table 1.

Table 1. Structure of a VGG-16-like CNN used for MD in [24]

View Table

The VGG-16-like CNN consists of $N_c = 13$ convolutional layers. Each convolutional filter is of 3x3 pixels size. So, the number of multiplications per 1 convolution is $n = 9$. Padding: same was used on each convolutional layer. So, the number of convolutions per layer is the number of pixels on this layer. For a layer of size $M$ by $M$ pixels, it is $M^2$ convolutions.

The total number of multiplications needed during propagation through $N_c=13$ convolutional layers is

(17)$$Q_{CL}=\sum_{i=1}^{N_c} nM_i^2C_i,$$

where $M_i$ is the size of $i$-th layer and $C_i$ is the number of convolutional filters applied at $i$-th layer. The total number of multiplications per 1 FC layer is the product of the number of nodes in the previous layer and the number of nodes in the current layer. So, for the first FC layer with $F_1$ elements, it is defined as

(18)$$Q^{(1)}_{FC}=M_{13}^2C_{13}\cdot F_1.$$

The second FC layer is of size $F_2=2N-1$, so the total number of multiplications for the second layer is

(19)$$Q^{(2)}_{FC}=F_1\cdot F_2=F_1(2N-1).$$

Then to recover phases from the cosine values the CNN-based method requires the calculation of intensity distributions for each set of possible phase values and then compare them to the true intensity distribution. For each of the $N-1$ cosine values, there are 2 possible phase values that correspond to this cosine value, so the total $2^{N-1}$ intensity calculations are needed. Each intensity calculation consists of $(2N-1)M^2$ multiplications as each of $M^2$ pixels for each of $N$ modes must be multiplied by its amplitude ($N$ multiplications) and phase ($N-1$ multiplications) values. So, the total number of multiplications for deriving phase coefficients from their cosine values is

(20)$$Q_{\phi}=2^{N-1}\cdot(2N-1)M^2.$$

This leads to the total number of multiplications per 1 MD using CNN:

(21)$$\begin{aligned}Q_{CNN}&=Q_{CL}+Q_{FC}+Q_{\phi}= \\ &=\sum^{13}_{i=1}9M^2_iC_i+M^2_{13}C_{13}\cdot F_1+F_1\cdot F_2+2^{N-1}\cdot (2N-1)M^2 = \\ &=2430\cdot M^2+2048\cdot M^2 + 1024\cdot(2N-1)+2^{N-1}\cdot(2N-1)M^2 = \\ &=M^2\cdot(4478+2^{N-1}\cdot(2N-1))+1-24\cdot(2N-1). \end{aligned}$$

An approximation of the above equation can be given by

(22)$$Q_{CNN}\approx M^2\cdot(4478+2^{N-1}\cdot(2N-1)).$$

The NFF MD method requires multiplications of complex values. Each multiplication of complex-valued numbers requires 4 multiplications of real-value numbers. Moreover, the number of pixels is doubled as we analyze both NF and FF intensity distributions. Therefore, for the complete MD algorithm the number of multiplications is

(23)$$Q_{NFF}(N,M)\approx24M^2(N-1).$$

Figure 12 shows how the performance gain depends on the number of modes. The NFF MD algorithm requires 8 times more multiplications than the FF only method, but it benefits from both higher accuracy and ability to recover phase coefficients unambiguously (not only up to transformation (12)) which is the main advantage in comparison to CNN-based MD algorithm that does not allow to completely decompose phases [24].

Fig. 12. The ratio between the numbers of multiplications needed for a single-mode decomposition using the CNN algorithm and the proposed algorithm.

Download Full Size | PDF

Thus, the proposed techniques require orders of magnitude less multiplications than the best CNN-based MD algorithm presented to date.

5. Experimental verification

We experimentally tested the FF algorithm for a 5-mode fiber. We assembled an experimental setup shown in Fig. 1(b). We used a piece of SM2000 fiber and a laser with a central wavelength of 1060 nm. This fiber supports 5 modes at the laser operating wavelength. The corresponding V number of the fiber at the specified wavelength is 3.91. A lens with a focal length of 4.5 mm was used to convert NF distribution at the fiber facet to the FF distribution on the camera. We used a polarizer beyond the lens to select only one polarization. The output intensity distribution after the lens was captured by a CCD camera Ximea mq013rg-e2. The digital decomposition speed was much higher than the speed of the camera (60 frames per second) and to test the accuracy and performance of the algorithm without limitations of camera speed, a series of images was recorded and then processed. The complex mode structure was achieved by bending and twisting the FMF.

During the pre-processing each picture was cropped to size of 400x400 and the average background intensity was subtracted. The background intensity was measured as a captured intensity distribution far outside the main beam region where the intensity of the main beam is guaranteed to be less than one thousandth of the maximum intensity of the main beam. The noise intensity was measured as a variance of the background intensity. The captured intensity distribution, the background intensity distribution and the histogram of the background are shown in Fig. 13.

Fig. 13. a) Measured intensity, b) Measured intensity far outside of the main beam region or background and c) Background intensity histogram

Download Full Size | PDF

The calculated SNR value was of approximately 21 dB. We believe that two main sources of noise are present in our case: the light scattered by optical elements and the temperature noise of the camera. To find the center of the beam profile we captured 100 images and guessed some initial values of coordinates of the center as well as the scaling factor of the image. Then we calculated eigenmodes corresponding to these guessed parameters, then performed MD, then reconstructed intensity patterns from these recovered mode weights and finally calculated the mean correlation between measured and reconstructed intensity distributions. Then we numerically optimized the 3 free parameters (coordinates of the center and scaling factor) to maximize mean correlation across all these 100 images. The numerical optimization was performed by Nelder-Mead (simplex) method. The optimized parameters were then used for processing all other experimentally captured intensity distributions. Then the MD algorithm was applied to each intensity distribution. The acquired weights were used to reconstruct the intensity using the known FF distributions and finally the accuracy of MD was analyzed by calculating correlation between the measured and reconstructed intensity distributions. Figure 14 shows examples of captured far field distributions, reconstructed intensities, discrepancies, and corresponding correlations.

Fig. 14. Examples of the experimental MD if 5-mode fiber.

Download Full Size | PDF

Figure 15 shows the correlation for the whole series of images captured and processed by FF algorithm. The mean decomposition time was 52 µs. The correlation is above 0.95 for the whole series of captured intensity distributions.

Fig. 15. Correlation between captured and reconstructed intensity distributions for a series of images of FF distribution in 5-mode fiber.

Download Full Size | PDF

6. Conclusion

We developed an intensity-only MD technique based on FF and NFF analysis which is superior compared to the fastest previously published MD algorithm in terms of both number of modes in noiseless case and noise sustainability for MD in the presence of noise. Compared to the best previous results, the number of modes in noiseless case increased from 27 to 46 for FF analysis only and 49 in case of using both FF and NF (NFF). Importantly, we demonstrate a substantial improvement of performance in the presence of noise: the SNR needed for MD with a full error of 0.1 decreased from 20 dB to 15 dB (FF) or 10 dB (NFF) for 3-mode fiber; from 30 dB to 22.5 dB (FF) or 19 dB (NFF) for 5-mode fiber. The NFF also benefits from the ability to completely decompose phase weights of higher modes relative to the fundamental mode. And another advantage of the FF technique is the use of simplified optical scheme compared to NF and NFF methods, that does not only decreases the number of required optical components but also increases the setup compactness what is crucial for the in field applications. The rigorous analysis of computational complexity is provided, and it has shown that the presented algorithm requires orders of magnitude less multiplications than a comparable machine learning-based MD technique. The developed progress in intensity-only MD is another step towards fully passive phase-sensitive multi-mode receivers.

Funding

Horizon 2020 Framework Programme (814276); Engineering and Physical Sciences Research Council (EP/R035342/1).

Acknowledgments

EM and SKT acknowledge the support of EPSRC project TRANSNET. AD acknowledges funding by the EU Horizon 2020 Marie Skłodowska-Curie Action ETN project WON (grant agreement 814276).

Disclosures

The authors declare no conflicts of interest.

Data availability

The data that supports the findings of this study is available from the corresponding author upon reasonable request.

References

1. Y. Wang, Y. Tang, S. Yan, and J. Xu, “High-power mode-locked 2 μm multimode fiber laser,” Laser Phys. Lett. 15(8), 085101 (2018). [CrossRef]

2. N. Bai, E. Ip, T. Wang, and G. Li, “Multimode fiber amplifier with tunable modal gain using a reconfigurable multimode pump,” Opt. Express 19(17), 16601–16611 (2011). [CrossRef]

3. Y. Qiu, J. Xu, K. K. Wong, and K. K. Tsia, “Exploiting few mode-fibers for optical time-stretch confocal microscopy in the short near-infrared window,” Opt. Express 20(22), 24115–24123 (2012). [CrossRef]

4. M. Diab, A. N. Dinkelaker, J. Davenport, K. Madhav, and M. M. Roth, “Starlight coupling through atmospheric turbulence into few-mode fibres and photonic lanterns in the presence of partial adaptive optics correction,” Mon. Not. R. Astron. Soc. 501(2), 1557–1567 (2020). [CrossRef]

5. A. Li, Y. Wang, Q. Hu, and W. Shieh, “Few-mode fiber based optical sensors,” Opt. Express 23(2), 1139–1150 (2015). [CrossRef]

6. L. G. Wright, D. N. Christodoulides, and F. W. Wise, “Controllable spatiotemporal nonlinear effects in multimode fibres,” Nat. Photonics 9(5), 306–310 (2015). [CrossRef]

7. K. Krupa, A. Tonello, B. M. Shalaby, M. Fabert, A. Barthélémy, G. Millot, S. Wabnitz, and V. Couderc, “Spatial beam self-cleaning in multimode fibres,” Nat. Photonics 11(4), 237–241 (2017). [CrossRef]

8. W. H. Renninger and F. W. Wise, “Optical solitons in graded-index multimode fibres,” Nat. Commun. 4(1), 1719 (2013). [CrossRef]

9. Y. Leventoux, G. Granger, K. Krupa, A. Tonello, G. Millot, M. Ferraro, F. Mangini, M. Zitelli, S. Wabnitz, S. Février, and V. Couderc, “3d time-domain beam mapping for studying nonlinear dynamics in multimode optical fibers,” Opt. Lett. 46(1), 66–69 (2021). [CrossRef]

10. H. K. Chandrasekharan, K. Ehrlich, M. G. Tanner, D. M. Haynes, S. Mukherjee, T. A. Birks, and R. R. Thomson, “Observing mode-dependent wavelength-to-time mapping in few-mode fibers using a single-photon detector array,” APL Photonics 5(6), 061303 (2020). [CrossRef]

11. P. J. Winzer, D. T. Neilson, and A. R. Chraplyvy, “Fiber-optic transmission and networking: the previous 20 and the next 20 years,” Opt. Express 26(18), 24190–24239 (2018). [CrossRef]

12. D. Richardson, J. M. Fini, and L. E. Nelson, “Space-division multiplexing in optical fibres,” Nat. Photonics 7(5), 354–362 (2013). [CrossRef]

13. M. Leo, R. Piccolo, C. Distante, P. Memmolo, M. Paturzo, and P. Ferraro, “Multilevel bidimensional empirical mode decomposition: a new speckle reduction method in digital holography,” Opt. Eng. 53(11), 112314 (2014). [CrossRef]

14. M. Lyu, Z. Lin, G. Li, and G. Situ, “Fast modal decomposition for optical fibers using digital holography,” Sci. Rep. 7(1), 6556 (2017). [CrossRef]

15. G. Labroille, N. Barré, O. Pinel, B. Denolle, K. Lenglé, L. Garcia, L. Jaffrès, P. Jian, and J.-F. Morizur, “Characterization and applications of spatial mode multiplexers based on multi-plane light conversion,” Opt. Fiber Technol. 35, 93–99 (2017). [CrossRef]

16. O. Shapira, A. F. Abouraddy, J. D. Joannopoulos, and Y. Fink, “Complete modal decomposition for optical waveguides,” Phys. Rev. Lett. 94(14), 143902 (2005). [CrossRef]

17. R. Brüning, P. Gelszinnis, C. Schulze, D. Flamm, and M. Duparré, “Comparative analysis of numerical methods for the mode analysis of laser beams,” Appl. Opt. 52(32), 7769–7777 (2013). [CrossRef]

18. H. Lü, P. Zhou, X. Wang, and Z. Jiang, “Fast and accurate modal decomposition of multimode fiber based on stochastic parallel gradient descent algorithm,” Appl. Opt. 52(12), 2905–2908 (2013). [CrossRef]

19. A. Liu, T. Lin, H. Han, X. Zhang, Z. Chen, F. Gan, H. Lv, and X. Liu, “Analyzing modal power in multi-mode waveguide via machine learning,” Opt. Express 26(17), 22100–22109 (2018). [CrossRef]

20. L. Hofer, L. Jones, J. Goedert, and R. Dragone, “Hermite–gaussian mode detection via convolution neural networks,” J. Opt. Soc. Am. A 36(6), 936–943 (2019). [CrossRef]

21. H. Hu, Z. Li, Y. Zhu, and L. Sun, “Non-degenerate mode power decomposition in optical fiber without prior knowledge,” Optik 206, 164354 (2020). [CrossRef]

22. S. Rothe, Q. Zhang, N. Koukourakis, and J. W. Czarske, “Intensity-only mode decomposition on multimode fibers using a densely connected convolutional network,” arXiv preprint arXiv:2008.00864 (2020).

23. Y. An, L. Huang, J. Li, J. Leng, L. Yang, and P. Zhou, “Deep learning-based real-time mode decomposition for multimode fibers,” IEEE J. Sel. Top. Quantum Electron. 26(4), 1–6 (2020). [CrossRef]

24. Y. An, L. Huang, J. Li, J. Leng, L. Yang, and P. Zhou, “Learning to decompose the modes in few-mode fibers with deep convolutional neural network,” Opt. Express 27(7), 10127–10137 (2019). [CrossRef]

25. W. Yan, X. Xu, and J. Wang, “Modal decomposition for few mode fibers using the fractional fourier system,” Opt. Express 27(10), 13871–13883 (2019). [CrossRef]

26. E. S. Manuylovich, V. V. Dvoyrin, and S. K. Turitsyn, “Fast mode decomposition in few-mode fibers,” Nat. Commun. 11(1), 5507 (2020). [CrossRef]

Block	1	2	3	4	5	6	7
Number of layers	2	2	3	3	3	1	1
Image size	$M \times M$	$\frac{M}{2} \times \frac{M}{2}$	$\frac{M}{4} \times \frac{M}{4}$	$\frac{M}{8} \times \frac{M}{8}$	$\frac{M}{16} \times \frac{M}{16}$	$1 \times 1$	$1 \times 1$
Number of filters	64	128	256	512	512	1024	2N-1

Intensity-only-measurement mode decomposition in few-mode fibers

Abstract

1. Introduction

2. Methods

3. Results

4. Computational complexity

5. Experimental verification

6. Conclusion

Funding

Acknowledgments

Disclosures

Data availability

References

Data availability

Cited By

Figures (15)

Tables (1)

Equations (23)

Optics Express