Contribution to the Improvement of the Correlation Filter Method for Modal Analysis with a Spatial Light Modulator

Modal decomposition of light is essential to study its propagation properties in waveguides and photonic devices. Modal analysis can be carried out by implementing a computer-generated hologram acting as a match filter in a spatial light modulator. In this work, a series of aspects to be taken into account in order to get the most out of this method are presented, aiming to provide useful operational procedures. First of all, a method for filter size adjustment based on the standard fiber LP-mode symmetry is presented. The influence of the mode normalization in the complex amplitude encoding-inherent noise is then investigated. Finally, a robust method to measure the phase difference between modes is proposed. These procedures are tested by wavefront reconstruction in a conventional few-mode fiber.


Introduction
Characterization of optical fields by means of modal decomposition (MD) is key for the analysis, design, and optimization of multimode waveguide-based devices [1][2][3]. Interest in these techniques is increasing in parallel with the growing research on applications based on multimode structures such as high-power large-mode-area fiber lasers [4,5] and space division multiplexing techniques, aimed at widening the channel capacity of optical communication systems [6][7][8]. Modal decomposition facilitates the study of modal competition [9,10], modal instabilities limiting the maximum power emission in fiber lasers [11,12], nonlinear optical effects where the interaction of many transversal modes is involved [13][14][15], and the design of optical devices based on optical power transfer [16,17] and selective mode excitation [18,19]. This analysis has proven to be a key tool to study the transmission properties of waveguides regarding its modal behavior due to the physical insight that can be obtained [20].
In general, MD methods can be classified as numerical or experimental. On one hand, several numerical MD algorithms based on measured optical intensity distributions have been reported [21]. Iterative methods such as the GS algorithm [22], the SPGD algorithm [23], and the genetic algorithm [24] show a high accuracy but are highly sensitive to the initial values due to the local minima problem and, in some cases, the necessary iterative process can be computationally intensive [25]. The newly emerged neural network methods [26][27][28] outperform the iterative methods in decomposition speed without having the local minima problems and showing a highly accurate performance, but require highperformance computers, a large amount of memory, and a long time for training the neural network. Recently, a non-iterative method based on a sophisticated analytical model has been published [29], achieving a very fast MD performance and without any optimization

Basics of the Correlation Filter Method
The notation to be employed is established in Figure 1, which shows a basic scheme of the setup necessary for modal decomposition using the CFM. Lenses L 1 and L 2 , positioned following a 4f configuration and act as a beam expansor, where the fiber output end is placed at the focus of lens L 1 ; the distance between lenses is the sum of their focal distances and the CGH is located at the focal plane of lens L 2 . Through this beam expansor, the fiber output electric field U(ε, µ) is imaged at the CGH plane, U(x, y). T(x, y) is the CGH transmittance implemented in the SLM, W 0 (x, y) is the electric field after the SLM, and W f (x, y) is the field at the CCD camera position, which is placed at the focal plane of lens L 3 .
Micromachines 2022, 13, x FOR PEER REVIEW 3 of 14 nally, in Section 7 we test the experimental performance of the modal-analysis procedure by a wavefront reconstruction, which leads to the conclusions in Section 8.

Basics of the Correlation Filter Method
The notation to be employed is established in Figure 1, which shows a basic scheme of the setup necessary for modal decomposition using the CFM. Lenses L1 and L2, positioned following a 4f configuration and act as a beam expansor, where the fiber output end is placed at the focus of lens L1; the distance between lenses is the sum of their focal distances and the CGH is located at the focal plane of lens L2. Through this beam expansor, the fiber output electric field U(ε,µ) is imaged at the CGH plane, U x,y . T x,y is the CGH transmittance implemented in the SLM, W 0 x,y is the electric field after the SLM, and W f x,y is the field at the CCD camera position, which is placed at the focal plane of lens L3. Figure 1. Scheme of the modal decomposition setup. L1, L2, and L3 are lenses. SLM is the spatial light modulator where the CGH is implemented. U(ε,µ) is the waveguide output electric field, U x,y is the electric field at the CGH plane after passing through the beam expansor in a 4f configuration, T x,y is the SLM transmittance, W 0 x,y is the electric field after the SLM, and W f x,y is the electric field at the CCD camera plane.
The multimode waveguide output field to be analyzed constitutes the input CGH field U x,y = c n ϕ n x,y N n=1 .
(1) where x and y are orthogonal coordinates in the CGH plane, N is the total number of propagation modes allowed by the waveguide under study [47], ϕ n x,y is the nth-mode normalized distribution of the transversal electric field, and c n is its complex expansion coefficient, which can be expressed in terms of a modal weight, |c n |, and a relative phase term, φ n , c n = |c n |e iφ n , where i is the imaginary unit and where all mode relative phases are expressed with respect to the 0-mode phase φ 0 = 0 . The electric fields W 0 x,y and U x,y are related by: The Fourier transform of the output distribution, W f (u,v), is obtained at the L3 lens focal plane. Specifically, W f and W 0 are related by means of [48]: Figure 1. Scheme of the modal decomposition setup. L 1 , L 2 , and L 3 are lenses. SLM is the spatial light modulator where the CGH is implemented. U(ε, µ) is the waveguide output electric field, U(x, y) is the electric field at the CGH plane after passing through the beam expansor in a 4f configuration, T(x, y) is the SLM transmittance, W 0 (x, y) is the electric field after the SLM, and W f (x, y) is the electric field at the CCD camera plane.
The multimode waveguide output field to be analyzed constitutes the input CGH field where x and y are orthogonal coordinates in the CGH plane, N is the total number of propagation modes allowed by the waveguide under study [47], φ n (x, y) is the nth-mode normalized distribution of the transversal electric field, and c n is its complex expansion coefficient, which can be expressed in terms of a modal weight, |c n |, and a relative phase term, ϕ n , c n = |c n |e iϕ n , where i is the imaginary unit and where all mode relative phases are expressed with respect to the 0-mode phase (ϕ 0 = 0). The electric fields W 0 (x, y) and U(x, y) are related by: The Fourier transform of the output distribution, W f (u, v), is obtained at the L 3 lens focal plane. Specifically, W f and W 0 are related by means of [48]: where F denotes Fourier transform, λ is the input-light wavelength, f is the lens focal distance, and the coordinates at the focal plane (u, v) are related to the Fourier-transform frequency-space coordinates f x and f y by means of u = λff x , v = λff y . In order to measure the amplitude of mode p, a filter with the following transmittance should be employed [36]: where the asterisk denotes the complex conjugate. After some operations and taking into account orthogonality between mode field distributions [36]: where δ pn is the Kronecker delta. If a CCD detector is placed on axis in the focal plane of L 3 , the use of a filter, as given by Equation (5), allows one to obtain information regarding the p-mode amplitude |c p | at the specific focal plane coordinates (u = 0, v = 0), the intensity at this point being proportional to the modal weight: A simultaneous procedure based on angular multiplexing can be carried out. It consists of employing a filter whose transmittance is a superposition of different T p filters, associating to each p-mode a different wave vector k x,p , k y,p , so that the total transmittance function is: By a convenient choice of these wave vectors, information regarding the different mode modal weights appears at separated-enough points in the focal plane, u p = λfk x,p and v p = λfk y,p being the coordinates of each one of them.
For some applications, it is enough to obtain the modal weights, |c p | 2 . However, sometimes it is required to determine the phase difference between modes. According to [36], in order to measure it, it is necessary to use two transmittance functions, defined as and so that the p-mode relative phase, ϕ p , can be determined with respect to the 0-mode one. In these cases, the intensity at the CCD plane is proportional to the cosine and sine, respectively, of the phase difference [36]. This allows one to obtain an unambiguous solution for the phase differences.

Double-Phase Method for Complex Amplitude Encoding
When implementing the CGH in a phase SLM, all the transmittance functions defined in the previous subsection require the encoding of a complex amplitude in an only-phase device. Among the different proposed procedures for pixel grouping, the double-phase CGHs have been chosen [45] due to their simple implementation and maximum resolution, as in this technique macropixels are composed of only two pixels. This method is based on the principle that any complex value inside the unit circle can be resolved into the sum of two complex values with unit modules.
Consider a CGH implemented in a commercial SLM that only admits phase modulation of its pixels, and a combination of them in macropixels that may offer an approximate complex modulation. Specifically, as explained in [45], a CGH filter with an approximate transmittance to the theoretical one can be achieved by means of grouping pairs of pixels into a macropixel. Considering that each couple of consecutive pixels constitute the macropixel with indexes (m, n), if the desired complex transmittance of the macropixel (m, n) is T mn = |T mn |e iα mn , the following respective phase shifts are assigned to the two consecutive pixels [45]: where |T mn | ≤ 1, 0 ≤ α mn ≤ 2π and ∆ mn = cos −1 (|T mn |). (11) Denoting by Q the transmittance function of the phase-only device where the doublephase CGH are implemented, its Fourier transform can be expressed as: in such a way that [45] ∼ Contribution ∼ Q S , referred to as signal, is the only one that should be expected if the double-phase solution provided an exact T(x, y) transmittance. Nevertheless, another contribution to ∼ Q also appears, due to the double-phase implementation: ∼ Q N , referred to as noise.

Modal Analysis Setup
The required experimental setup is depicted in Figure 2. An SMF28 optical fiber supporting six propagation LP-modes is illuminated by a He-Ne laser (λ = 632.8 nm). The first section of optical fiber is repeatedly curved in order to eliminate modes other than the fundamental one, LP 01 . We use this fiber to illuminate a second optical fiber with the same characteristics. This way, when both fibers are placed closely one after another over the same longitudinal axis, the LP 01 mode is the only one that propagates along the second section of optical fiber. It can be used to carry out the system alignment. By sideways displacing the first fiber with respect to the second one, different transversal distributions are generated, due to different combinations of excited modes, which is useful for testing the MD performance. Figure 2A,B show two examples: the first one with both fiber sections totally aligned, used to calibrate the setup, and the second one at an arbitrary position, employed to test the proposed technique's performance.
imate complex modulation. Specifically, as explained in [45], a CGH filter with an approximate transmittance to the theoretical one can be achieved by means of grouping pairs of pixels into a macropixel. Considering that each couple of consecutive pixels constitute the macropixel with indexes (m, n), if the desired complex transmittance of the macropixel (m, n) is T mn = |T mn |e iα mn , the following respective phase shifts are assigned to the two consecutive pixels [45]: α mn (1) = α mn − Δ mn and α mn (2) = α mn +Δ mn , (10) where |T mn | ≤ 1, 0 ≤ α mn ≤ 2π and Δ mn = cos −1 (|T mn |) .
Denoting by Q the transmittance function of the phase-only device where the double-phase CGH are implemented, its Fourier transform can be expressed as: in such a way that [45] Q S = ℱ T x, y .
Contribution Q S , referred to as signal, is the only one that should be expected if the double-phase solution provided an exact T x,y transmittance. Nevertheless, another contribution to Q also appears, due to the double-phase implementation: Q N , referred to as noise.

Modal Analysis Setup
The required experimental setup is depicted in Figure 2. An SMF28 optical fiber supporting six propagation LP-modes is illuminated by a He-Ne laser (λ = 632.8 nm). The first section of optical fiber is repeatedly curved in order to eliminate modes other than the fundamental one, LP . We use this fiber to illuminate a second optical fiber with the same characteristics. This way, when both fibers are placed closely one after another over the same longitudinal axis, the LP 01 mode is the only one that propagates along the second section of optical fiber. It can be used to carry out the system alignment. By sideways displacing the first fiber with respect to the second one, different transversal distributions are generated, due to different combinations of excited modes, which is useful for testing the MD performance. Figure 2A   In order to increase the effective resolution, the fiber output beam size has been magnified at the SLM display using a microscope objective (DIN × 40) and a lens (f 1 = 50 cm) following a 4f configuration (theoretical magnification 125). The 4f-lens combination images the fiber output field distribution through the beam splitter (BS) both in the CCD 1 camera and in the phase-only SLM operating in reflection mode (PLUTO VIS Phase Only Spatial Light Modulator from Holoeye [49]). As it only operates in the horizontal polarization, we select this state of polarization from the input beam with a linear polarizer (LP). In order to obtain a complete MD, the same procedure should be performed for the vertical polarization. In order to do that, the horizontal LP should be replaced by a vertical one, followed by a half-wave plate. The reflected field from the SLM, where the CGH is implemented, is Fourier transformed by lens L 2 (f 2 = 15 cm), and this signal is detected by the CCD 2 camera, placed at the L 2 lens focal plane, where the modal weights are obtained by measuring the intensity in the specified coordinates according to Equation (6). In order to avoid interferences, The CCD 2 and the SLM are covered when measuring intensity profiles with the CCD 1 camera. The SLM was previously calibrated following the procedure described in [50].

Correlation-Filter Size adjustment
The MD performance is based on the overlap in the SLM between the incident light beam and the implemented CGH. Inaccuracies in the position or size of the encoded CGH with respect to the incident light can result in the detection of erroneous modes or in the incorrect determination of the modal weights. Notwithstanding its importance, these issues have not been studied in depth, except for in a recently published tutorial [41]. Their proposed approach to transversally align the CGH position consists of displacing the filter position in such a way that, for a mismatched overlap, the on-axis null appears centered regarding the input light beam. If the first and second optical fiber sections are aligned so as to achieve single LP 01 mode propagation in the second one, there should be zero on-axis intensity when implementing any other mode into the match filter. We can benefit from LP e 11 and LP o 11 mode symmetry in order to center the CGH vertically and horizontally, respectively. The simulated transversal intensity distributions at the CCD 2 plane, when the LP e 11 and LP e 11 modes are implemented in the SLM, are shown in Figure 3a,b, respectively, for different cases of transversal centering. As one can see, when the filters are centered with respect to the incident field, there is zero intensity in the optical axis point (marked with a red cross). As the filters decenter, the intensity at the measurement position increases.
In order to increase the effective resolution, the fiber output beam size has been magnified at the SLM display using a microscope objective (DIN × 40) and a lens (f1 = 50 cm) following a 4f configuration (theoretical magnification 125). The 4f-lens combination images the fiber output field distribution through the beam splitter (BS) both in the CCD1 camera and in the phase-only SLM operating in reflection mode (PLUTO VIS Phase Only Spatial Light Modulator from Holoeye [49]). As it only operates in the horizontal polarization, we select this state of polarization from the input beam with a linear polarizer (LP). In order to obtain a complete MD, the same procedure should be performed for the vertical polarization. In order to do that, the horizontal LP should be replaced by a vertical one, followed by a half-wave plate. The reflected field from the SLM, where the CGH is implemented, is Fourier transformed by lens L2 (f2 = 15 cm), and this signal is detected by the CCD2 camera, placed at the L2 lens focal plane, where the modal weights are obtained by measuring the intensity in the specified coordinates according to Equation (6). In order to avoid interferences, The CCD2 and the SLM are covered when measuring intensity profiles with the CCD1 camera. The SLM was previously calibrated following the procedure described in [50].

Correlation-Filter Size adjustment
The MD performance is based on the overlap in the SLM between the incident light beam and the implemented CGH. Inaccuracies in the position or size of the encoded CGH with respect to the incident light can result in the detection of erroneous modes or in the incorrect determination of the modal weights. Notwithstanding its importance, these issues have not been studied in depth, except for in a recently published tutorial [41]. Their proposed approach to transversally align the CGH position consists of displacing the filter position in such a way that, for a mismatched overlap, the on-axis null appears centered regarding the input light beam. If the first and second optical fiber sections are aligned so as to achieve single LP mode propagation in the second one, there should be zero on-axis intensity when implementing any other mode into the match filter. We can benefit from LP 11 e and LP 11 o mode symmetry in order to center the CGH vertically and horizontally, respectively. The simulated transversal intensity distributions at the CCD2 plane, when the LP 11 e and LP 11 e modes are implemented in the SLM, are shown in Figure 3a,b, respectively, for different cases of transversal centering. As one can see, when the filters are centered with respect to the incident field, there is zero intensity in the optical axis point (marked with a red cross). As the filters decenter, the intensity at the measurement position increases. We propose that this same idea can be used to adjust the filter size. Despite it being straightforward to calculate the theoretical 4f-magnification, this value can be slightly modified in the experimental implementation, mainly because it is not easy to place the We propose that this same idea can be used to adjust the filter size. Despite it being straightforward to calculate the theoretical 4f-magnification, this value can be slightly modified in the experimental implementation, mainly because it is not easy to place the optical fiber output end at the exact focal length distance from the microscope principal plane. In Ref. [51], a method for filter-size adjustment is proposed, for which it is required to obtain some parameters, such as the beam quality factor and the second moment. Here, we present a different approach, based on the LP 02 mode symmetry. If we implement the LP 02 mode into the CGH and illuminate it with the LP 01 mode, only when both their sizes are matched will there be zero intensity at the optical axis, as shown in the simulations from Figure 4.
plane. In Ref. [51], a method for filter-size adjustment is proposed, for which it is required to obtain some parameters, such as the beam quality factor and the second moment. Here, we present a different approach, based on the LP 02 mode symmetry. If we implement the LP 02 mode into the CGH and illuminate it with the LP 01 mode, only when both their sizes are matched will there be zero intensity at the optical axis, as shown in the simulations from Figure 4. We can use this behavior to find the correct scale just by minimizing the optical axis intensity as a function of the theoretical magnification value. This is shown experimentally in Figure 5. Overlaying the graph, three measured transversal intensity distributions are shown for a smaller, adjusted, and larger size (from left to right) of the filter with respect to the incident light. The measured intensity distributions follow a similar behavior than the simulated ones in Figure 4. According to the minimum intensity, the magnification value is 141, which has a 13% relative error with respect to the theoretical one. This error can be explained due to the inherent uncertainties of the optical system positioning. This method can be used as long as the set of modes has the required symmetry, which is common in cylindrical-like structures.

Double Phase Method Noise Term Analysis
Consider a double-phase CGH with transmittance Q x,y (Equation (18)), implemented in order to emulate an ideal CGH filter with complex transmittance T x,y (Equation (13)). According to Equations (3), (4), and (12), this double-phase filter yields a field distribution in the L2 back focal plane, such as: We can use this behavior to find the correct scale just by minimizing the optical axis intensity as a function of the theoretical magnification value. This is shown experimentally in Figure 5. Overlaying the graph, three measured transversal intensity distributions are shown for a smaller, adjusted, and larger size (from left to right) of the filter with respect to the incident light. The measured intensity distributions follow a similar behavior than the simulated ones in Figure 4.
plane. In Ref. [51], a method for filter-size adjustment is proposed, for which it is required to obtain some parameters, such as the beam quality factor and the second moment. Here, we present a different approach, based on the LP 02 mode symmetry. If we implement the LP 02 mode into the CGH and illuminate it with the LP 01 mode, only when both their sizes are matched will there be zero intensity at the optical axis, as shown in the simulations from Figure 4. We can use this behavior to find the correct scale just by minimizing the optical axis intensity as a function of the theoretical magnification value. This is shown experimentally in Figure 5. Overlaying the graph, three measured transversal intensity distributions are shown for a smaller, adjusted, and larger size (from left to right) of the filter with respect to the incident light. The measured intensity distributions follow a similar behavior than the simulated ones in Figure 4. According to the minimum intensity, the magnification value is 141, which has a 13% relative error with respect to the theoretical one. This error can be explained due to the inherent uncertainties of the optical system positioning. This method can be used as long as the set of modes has the required symmetry, which is common in cylindrical-like structures.

Double Phase Method Noise Term Analysis
Consider a double-phase CGH with transmittance Q x,y (Equation (18)), implemented in order to emulate an ideal CGH filter with complex transmittance T x,y (Equation (13)). According to Equations (3), (4), and (12), this double-phase filter yields a field distribution in the L2 back focal plane, such as: According to the minimum intensity, the magnification value is 141, which has a 13% relative error with respect to the theoretical one. This error can be explained due to the inherent uncertainties of the optical system positioning. This method can be used as long as the set of modes has the required symmetry, which is common in cylindrical-like structures.

Double Phase Method Noise Term Analysis
Consider a double-phase CGH with transmittance Q(x, y) (Equation (18)), implemented in order to emulate an ideal CGH filter with complex transmittance T(x, y) (Equation (13)). According to Equations (3), (4), and (12), this double-phase filter yields a field distribution in the L 2 back focal plane, such as: where U = F [U(x, y)] and the ( * ) sign represents convolution product. We consider two different contributions to W f , called here signal and noise terms, respectively: Because Q S = F [T(x, y)], W f,S represents the field distribution that would be obtained if an ideal CGH filter with complex transmittance T(x, y) was employed. Nevertheless, a superimposed noise term W f,N is also present. As stated in Equation (16), this term contains the convolution of two factors. One of them, ∼ Q N , is intrinsic to the filter, while the other, ∼ U, is directly related to the input electric field distribution, which is obviously unpredictable as it is the object of analysis. Therefore, the magnitude and distribution of the noise term is expected to be very different depending on each working condition. For this reason, its impact cannot be analyzed by means of a general mathematical treatment.
A case in which the noise factor strongly affects the MD is shown in Figure 6. Suppose an incident beam with all six LP-modes equally present and a CGH in which all the match filters (one for each mode) are simultaneously multiplexed with a different grating to spatially separate the signals at the Fourier plane. Specifically, the vertexes of a regular hexagon, as shown in Figure 6a-c, show the simulated intensity at the L 2 back focal plane, separated into the signal term W f,S (Equation (15)) and the noise term W f,N (Equation (16)), respectively.
different contributions to W f , called here signal and noise terms, respectively: Because Q S = ℱ T x,y , W f,S represents the field distribution that would be obtained if an ideal CGH filter with complex transmittance T x,y was employed. Nevertheless, a superimposed noise term W f,N is also present. As stated in Equation (16), this term contains the convolution of two factors. One of them, Q N , is intrinsic to the filter, while the other, U , is directly related to the input electric field distribution, which is obviously unpredictable as it is the object of analysis. Therefore, the magnitude and distribution of the noise term is expected to be very different depending on each working condition. For this reason, its impact cannot be analyzed by means of a general mathematical treatment.
A case in which the noise factor strongly affects the MD is shown in Figure 6. Suppose an incident beam with all six LP-modes equally present and a CGH in which all the match filters (one for each mode) are simultaneously multiplexed with a different grating to spatially separate the signals at the Fourier plane. Specifically, the vertexes of a regular hexagon, as shown in Figure 6a-c, show the simulated intensity at the L2 back focal plane, separated into the signal term W f,S (Equation (15)) and the noise term W f,N (Equation (16)), respectively. By comparing the modal weight percentages with and without noise terms, relative errors up to 50% are obtained, as shown in Table 1, thus preventing the correct performance of the MD. The modal weight in percentage and the relative error are computed by |c l | 2 (%) = 100|c l | 2 / ∑ |c n | 2 N n=1 and ε r = 100 |c| Exp 2 -|c| Teo 2 /|c| Teo 2 , respectively, where |c| 2 and |c| 2 are the experimental and theoretical modal weights, the last one being 16.67% as all modes are equally present. Table 1. Simulated modal weights (obtained with the W f function) and relative errors (obtained by comparing them with the ones determined through the W f,S function) following the simultaneous configuration of Figure 6 when all six LP-modes are equally present in the incident light.  By comparing the modal weight percentages with and without noise terms, relative errors up to 50% are obtained, as shown in Table 1, thus preventing the correct performance of the MD. The modal weight in percentage and the relative error are computed by |c l | 2 (%) = 100|c l | 2 / ∑ N n=1 |c n | 2 and ε r = 100(|c| 2 Exp − |c| 2 Teo )/|c| 2 Teo , respectively, where |c| 2 Exp and |c| 2 Teo are the experimental and theoretical modal weights, the last one being 16.67% as all modes are equally present. Table 1. Simulated modal weights (obtained with the W f function) and relative errors (obtained by comparing them with the ones determined through the W f,S function) following the simultaneous configuration of Figure 6 when all six LP-modes are equally present in the incident light. 20. 3 22 In view of the findings, it is essential to reduce the noise term influence on the MD performance. One way to achieve this is to realize the MD sequentially rather than simultaneously. However, it is not possible to obtain all modal weights in real time with this procedure, which can be an important disadvantage depending on the system dynamics. Another approach is to normalize the transmittance function (Equation (7)). Theoretically, any transmittance value between zero and one is valid when implementing the double phase method. Nevertheless, in order to use the whole range of the SLM and reduce the noise impact, the maximum value of the transmittance function amplitude should be as close as possible to one. This affects parameter ∆ mn from Equation (11), which is necessary for the implementation of the double phase method. The normalization has a direct impact on the signal to noise ratio, as shown in Figure 7, where we have used the same configuration than the one in Figure 6 and summed the intensities at the six information positions for both the signal (Equation (15)) and the noise (Equation (16)) terms. procedure, which can be an important disadvantage depending on the system dyn Another approach is to normalize the transmittance function (Equation (7)). Th cally, any transmittance value between zero and one is valid when implementin double phase method. Nevertheless, in order to use the whole range of the SLM a duce the noise impact, the maximum value of the transmittance function amp should be as close as possible to one. This affects parameter Δ mn from Equatio which is necessary for the implementation of the double phase method. The norm tion has a direct impact on the signal to noise ratio, as shown in Figure 7, where w used the same configuration than the one in Figure 6 and summed the intensities six information positions for both the signal (Equation (15)) and the noise (Equatio terms.

Figure 7.
Signal and noise intensity evolution as a function of the maximum value in the tra tance function (Equation (7)) for the configuration case shown in Figure 6.
As shown in Figure 7, by adjusting the maximum value of the transmittanc plitude to one, we improve the signal to noise ratio, while if we keep the set of orthonormality (situation marked with the red arrows), the noise intensity cannot glected with respect to the signal, thus affecting the MD. When all the modes are ured using the same transmittance function (simultaneous MD case), the renormali appears as a constant factor equal for all measured modal weights. However, the malization factors must be considered when measuring different modal weight different transmittance functions in order not to break the mode orthonormality.

Robustness Improvement in Relative Phase Measurements
In theory, the phase difference between modes can be unambiguously deter from two measurements by using Equations (8) and (9), each one having two po solutions. Nevertheless, the mode relative phase measurements can be easily influ by system instabilities (i.e., laser power fluctuations, micro-positioners looseness, e by misalignments due to position uncertainties (i.e., lenses, fiber, CCD, SLM). As sequence, there may be an error in the arccosine and arcsine result that could p correct phase retrieval. For such cases, we propose a different transmittance funct order to measure the relative phases and also to study the phase determination tainty.
Suppose we implement in the SLM the following transmittance function: T x,y = 1 √2 ϕ 0 * x,y +ϕ l * x,y ·e iθ , Figure 7. Signal and noise intensity evolution as a function of the maximum value in the transmittance function (Equation (7)) for the configuration case shown in Figure 6.
As shown in Figure 7, by adjusting the maximum value of the transmittance amplitude to one, we improve the signal to noise ratio, while if we keep the set of modes orthonormality (situation marked with the red arrows), the noise intensity cannot be neglected with respect to the signal, thus affecting the MD. When all the modes are measured using the same transmittance function (simultaneous MD case), the renormalization appears as a constant factor equal for all measured modal weights. However, the renormalization factors must be considered when measuring different modal weights with different transmittance functions in order not to break the mode orthonormality.

Robustness Improvement in Relative Phase Measurements
In theory, the phase difference between modes can be unambiguously determined from two measurements by using Equations (8) and (9), each one having two possible solutions. Nevertheless, the mode relative phase measurements can be easily influenced by system instabilities (i.e., laser power fluctuations, micro-positioners looseness, etc.) or by misalignments due to position uncertainties (i.e., lenses, fiber, CCD, SLM). As a consequence, there may be an error in the arccosine and arcsine result that could prevent correct phase retrieval. For such cases, we propose a different transmittance function in order to measure the relative phases and also to study the phase determination uncertainty.
In order to improve the robustness of the relative phase determination, we perform a series of intensity measurements as a function of θ and fit them to the function in Equation (18).
As an example, Figure 8 shows the so mentioned measurements of the phase difference between the LP e 21 and LP 01 modes, when the input distribution is the one shown in Figure 2b, together with its best fit by least squares. It is clear to see that, despite the fact that the experimental dots follow a similar behavior than the cosine function, there are some significant differences with respect to the fitting function. Both instabilities (i.e., laser fluctuations, fiber disturbances) and experimental inaccuracies (i.e., devices alignments) may be responsible for the lack of agreement between the measurements and the theoretical function.
In order to improve the robustness of the relative phase determination, we perform a series of intensity measurements as a function of θ and fit them to the function in Equation (18).
As an example, Figure 8 shows the so mentioned measurements of the phase difference between the LP 21 e and LP 01 modes, when the input distribution is the one shown in Figure 2b, together with its best fit by least squares. It is clear to see that, despite the fact that the experimental dots follow a similar behavior than the cosine function there are some significant differences with respect to the fitting function. Both instabilities (i.e., laser fluctuations, fiber disturbances) and experimental inaccuracies (i.e., devices alignments) may be responsible for the lack of agreement between the measurements and the theoretical function. Superimposed to the graph in Figure 8, four points have been highlighted. First, two points, θ = 0 and θ = π/2, are the particular cases from Equations (8) and (9). In this case the points show a good agreement with the fitting curve. In fact, following the two-measurements procedure, a phase difference of 0.23π rad is obtained, which is almost the same value as the one determined through the fitting approach: 0.24π rad Based on this result, one could think that the fitting approach does not improve the phase determination. However, the relative phase obtained through the two-measurements procedure may depend heavily on the selected points, as one can see in Figure 8. Suppose we choose the particular cases θ = π and θ = 3π/2, whose experimental values are far from the fitting curve, and compute the phase difference. In this case, the obtained angles do not match each other. The obtained phases that fall in the same quadrant are 0 and 0.39π rad, respectively. It is clear to see that, in general, the function fit allows us to obtain a relative phase less influenced by possible instabilities and uncertainties. In the next section, we observe its impact by performing a wavefront reconstruction.

Wavefront Reconstruction
In order to test the proposed size adjustment technique and phase retrieval procedure, together with the MD method itself, we have reconstructed the transversal intensity distribution shown in Figure 2b. Figure 9 shows the comparison between the experi- Superimposed to the graph in Figure 8, four points have been highlighted. First, two points, θ = 0 and θ = π/2, are the particular cases from Equations (8) and (9). In this case, the points show a good agreement with the fitting curve. In fact, following the two-measurements procedure, a phase difference of 0.23π rad is obtained, which is almost the same value as the one determined through the fitting approach: 0.24π rad. Based on this result, one could think that the fitting approach does not improve the phase determination. However, the relative phase obtained through the two-measurements procedure may depend heavily on the selected points, as one can see in Figure 8. Suppose we choose the particular cases θ = π and θ = 3π/2, whose experimental values are far from the fitting curve, and compute the phase difference. In this case, the obtained angles do not match each other. The obtained phases that fall in the same quadrant are 0 and 0.39π rad, respectively. It is clear to see that, in general, the function fit allows us to obtain a relative phase less influenced by possible instabilities and uncertainties. In the next section, we observe its impact by performing a wavefront reconstruction.

Wavefront Reconstruction
In order to test the proposed size adjustment technique and phase retrieval procedure, together with the MD method itself, we have reconstructed the transversal intensity distribution shown in Figure 2b. Figure 9 shows the comparison between the experimental distribution (a) and the numerically reconstructed ones (b-d). The reconstructed distributions have been obtained by measuring all six LP-mode weights and relative phases. However, in the reconstructions shown in Figure 9b,c, the proposed CGH size adjustment and phase retrieval techniques have not been used, respectively. The reconstruction shown in Figure 9b has been obtained by employing the proposed phase retrieval fitting approach but not the CGH size adjustment. Instead, its size has been computed theoretically with the 4f magnification. On the contrary, for the reconstruction shown in Figure 9c, the CGH size adjustment technique has been used, but the relative phases have been determined through the two-measurements method. Finally, the distribution shown in Figure 9d has been computed using both of the proposed techniques. As one can see, a good agreement is shown between Figure 9a,d distributions, highlighting the importance of both the size adjustment and the phase retrieval. The residual intensity patterns, computed where m and r mean measured and reconstructed, are respectively provided, together with each reconstructed image. distributions have been obtained by measuring all six LP-mode weights and relative phases. However, in the reconstructions shown in Figure 9b,c, the proposed CGH size adjustment and phase retrieval techniques have not been used, respectively. The recon struction shown in Figure 9b has been obtained by employing the proposed phase re trieval fitting approach but not the CGH size adjustment. Instead, its size has been computed theoretically with the 4f magnification. On the contrary, for the reconstruction shown in Figure 9c, the CGH size adjustment technique has been used, but the relative phases have been determined through the two-measurements method. Finally, the dis tribution shown in Figure 9d has been computed using both of the proposed techniques As one can see, a good agreement is shown between Figure 9a,d distributions, high lighting the importance of both the size adjustment and the phase retrieval. The residua intensity patterns, computed where m and r mean measured and reconstructed, are re spectively provided, together with each reconstructed image. By comparing the transversal intensity distributions one can see that the CGH siz adjustment is quite critical and that, although a look alike distribution is obtained with out determining the relative phases with our proposed procedure, a better reconstruction is obtained when employed. To quantitatively compare the intensity distributions from Figure 9, the root-mean-square error (RMSE) has been computed for each of the recon structed images with respect to the experimentally registered one, all four of them pre viously normalized to the unit, obtaining the RMSE values 2.4 × 10 −2 , 1.8 × 10 −2 , and 1.7 × 10 −3 for the Figure 9b-d intensity distributions, respectively. Distribution from Figure 9d gives a RMSE one order of magnitude lower than the ones from Figure 9b,c, showing the good performance of the reconstruction when using the proposed procedures. Table 2 summarizes the modal weights and relative phases obtained by using the proposed techniques, with which distribution from Figure 9d is computed. By comparing the transversal intensity distributions one can see that the CGH size adjustment is quite critical and that, although a look alike distribution is obtained without determining the relative phases with our proposed procedure, a better reconstruction is obtained when employed. To quantitatively compare the intensity distributions from Figure 9, the root-mean-square error (RMSE) has been computed for each of the reconstructed images with respect to the experimentally registered one, all four of them previously normalized to the unit, obtaining the RMSE values 2.4 × 10 −2 , 1.8 × 10 −2 , and 1.7 × 10 −3 for the Figure 9b-d intensity distributions, respectively. Distribution from Figure 9d gives a RMSE one order of magnitude lower than the ones from Figure 9b,c, showing the good performance of the reconstruction when using the proposed procedures. Table 2 summarizes the modal weights and relative phases obtained by using the proposed techniques, with which distribution from Figure 9d is computed.

Conclusions
We have presented a series of practical procedures to perform an accurate modal decomposition of light, and we tested them satisfactorily by wavefront reconstruction in a few-mode fiber, allowing us to verify both the good performance of the modal analysis procedure and the proposed techniques. In order to reduce the SNR when performing a simultaneous MD and exploit the SLM range, the set of modes orthonormality needs to be broken, and thus a scale factor taken into account. It has been shown that it is possible to benefit from the LP-mode symmetry, not only to transversally center the CGH, but also to adjust its size with respect to the incident beam. This procedure is not restricted to the LP-modes, and it can be used as long as the set of modes provides the necessary symmetry, which is usually fulfilled in cylindrical-like waveguides. When setup instabilities make the correct phase difference determination through the classical approach difficult, a possibility to improve phase retrieval robustness consists of implementing a transmittance as a function of an angle parameter and adjust the set of measurements. The good performance of these techniques has been successfully tested by wavefront reconstruction.