1. Introduction

The applicability of plastic optical fibers (POFs) in data transmission networks, as an affordable and robust alternative to glass optical fibers, is limited by signal attenuation and bandwidth [1, 2]. Besides launching conditions and propagation lengths, fibers’ frequency responses are affected by external influences (e.g. fiber bends, mode scramblers) and intrinsic fiber properties (attenuation, mode coupling and modal dispersion). Provided that attenuation and mode coupling properties are known, the frequency responses of step-index (SI) multimode (MM) fibers can be predicted by Gloge’s time-dependent power flow equation [3]. Several methods for the estimation of mode coupling [4–6] and mode-dependent attenuation [7, 8] have been proposed. Whereas mode-dependent attenuation has been measured directly, mode coupling has been estimated from farfield output patterns assuming modal independence of the coupling rate. Mateo et al. demonstrated in a more general approach that taking into account modal dependence of the coupling rate enables more accurate computation of farfield output patterns [9]. By comparison of both models Savović and Djordjevich found small deviations in farfield patterns and argued that it is reasonable to assume constant mode coupling [10].

What the previously published methods have in common is that the mode coupling rates are concluded from farfield output patterns. In contrast, we propose a method that utilizes SI-POF impulse responses for estimation of mode coupling and attenuation instead. We achieve this by optimization of these fiber characteristics so that simulated impulse responses deviate minimally from reference impulse responses of the fiber under test. In order to obtain the simulated impulse responses we apply a matrix-based finite-difference method [11] to solve Gloge’s time-dependent power flow equation [3]. By alternation between fiber lengths we iteratively optimize both, mode coupling and attenuation, by means of downhill simplex algorithms [12]. This method enables estimation of SI-POF characteristics and prediction of impulse responses, frequency responses, fiber bandwidths and coupling lengths. Since no spatial information concerning the field emitted by the fiber is needed, the method reveals fiber properties solely from measurements of impulse responses.

In this paper we apply our proposed method to the case of modal dependence of both, mode coupling and attenuation. However, alternative fiber property models can also be applied. This paper is organized as follows: at first, we briefly discuss Gloge’s power flow equation and recapitulate Mateo’s matrix method to solve it. Then, we explain the downhill simplex optimization scheme as well as our algorithm in general and its specific design for the simulations presented in this paper. Afterwards, we validate our approach against reference impulse responses computed from independently measured and validated fiber characteristics and discuss its applicability for impulse and frequency response prediction as well as for fiber bandwidths and coupling lengths determination.

2. Gloge’s power flow equation and matrix approach

For description of the spatial evolution of light propagating in a SI-POF, Gloge’s power flow equation can be applied [13]. Assuming that the characteristics of neighbouring modes are approximately equal, their discrete values are replaced by the continuous propagation angle θ in Gloge’s approach. As a consequence, the entire modal field can be treated as a continuum and mode coupling, i.e. energy transfer among the propagating modes caused by intrinsic perturbations, can be described by diffusion d(θ). Likewise, mode-dependent attenuation varies with θ and is given by a(θ). Introducing the temporal dimension to Gloge’s equation, the time-dependent power flow equation

For discretization of Eq. (2) a first-order forward difference for the z derivative and first- and second-order central differences for the first- and second-order angular derivatives are used. Thus, Eq. (2) becomes

The transformation matrix M(ω) = (A(ω) + D)^m carries all space-time information of the fiber contributing to power flow of a light pulse in this fiber. Once M is determined, fast computation of fiber optical output in space and time can be accomplished for any known initial spatio-temporal pulse characteristic, especially when matrix-based computation software is used. From the final spatio-temporal field distribution p(L,θ,ω) after fiber length L, the global frequency response h(L,ω) can be calculated by integration over all propagation angles

3. Proposed method for estimation of modal fiber characteristics

3.1. General description of the proposed method

Due to the tight relationship of attenuation and diffusion, the adaption of a simulated to a reference impulse response by optimization of these properties does not necessarily yield the correct, or respectively best possible, solutions for fiber attenuation and diffusion. The sensitivity of an algorithm to such erroneous convergence might be reduced by simultaneous instead of alternating optimization of attenuation and diffusion properties. Instead of a simultaneous optimization, we tackle this problem by alternating both the parameter to be optimized and the fiber length in our approach. Since variations of the reference impulse responses of two different fiber lengths are solely caused by the attenuation and diffusion characteristics of the reference fiber, the proposed algorithm prevents erroneous optimization of these properties by comparison of reference and simulated impulse responses at the two fiber lengths.

The algorithm of the proposed method starts with the reasonable assumption of mode-independent attenuation (see Fig. 1). For this assumption diffusion parameters are optimized so that the simulated impulse response H_sim(L_A, t) at fiber length L_A deviates minimally from the reference impulse response H_ref(L_A, t). In a next step the fiber length is altered and modal dependence of attenuation is introduced. Now, the attenuation parameters are optimized adapting H_sim(L_B, t) to H_ref(L_B, t). In the next iteration, these new attenuation parameters are used for another optimization of diffusion parameters at fiber length L_A, and as a results, it is most probable that better diffusion parameters will be found. The iterations proceed until the deviations of H_sim(L_A,B, t) and H_ref(L_A,B, t) converge and the best possible solution for fiber diffusion and attenuation are found.

 

Fig. 1 Program flow of optimization algorithm

Download Full Size | PDF

3.2. The downhill simplex optimization algorithm

As a measure for conformity of H_sim(L_A,B, t) and H_ref(L_A,B, t) their root-mean-square deviation (RMSD) is computed. The minimization of RMSD is performed by a downhill simplex (DS) algorithm. In general, DS algorithms optimize any objective function with q parameters by repetitive subsitution of the worst out of (q + 1) simplex vertices [12]. The substitutions comply with distinct rules that cause the simplex to move towards better values and to contract to the optimum (local) value in the q-dimensional space. In our approach the q parameters are given by the modal diffusion and attenuation functions, that will be introduced later. The simplex vertices correspond to different sets of these parameters and the objective is to minimize RMSD. Two of such DS algorithms are embedded in one main loop performing optimization at different fiber lengths (see Fig. 1). In addition to the launched spatio-temporal pulse distribution, P(L = 0,θ,t), an initial guess for preliminary constant attenuation, a(θ) = γ, is passed on to the very first DS algorithm. In this DS algorithm the d(θ)-parameters are optimized at fiber length L_A and the best diffusion parameter set is passed on to the second DS algorithm, which optimizes the a(θ)-parameters at fiber length L_B. Afterwards, the best a(θ)-parameters are handed over to the next DS algorithm, substituting the initial guess, a(θ) = γ, and potentially better d(θ)-parameters are computed this time owing to the improved attenuation parameters. As the algorithm advances, the deviations between the simulated and the reference pulse shapes reduce and converge to the best possible solution of fiber diffusion and attenuation.

Concerning the initial parameter sets forming the very first simplex in each DS algorithm, two points have to be emphasized: firstly, these sets are not equivalent to boundary conditions, since the simplex can move towards values beyond the initial vertices. But, secondly, depending on the choice of initial data sets it might happen that the DS algorithms converge to local, instead of global, minima. For reducing the likeliness of this erroneous convergence we vary size and position of the initial simplex at each main loop roundtrip based on previously obtained best diffusion or attenuation parameters, respectively. We have found that this strategy works for all of the data which we analyzed. An alternative solution to the local minima problems may be to run a genetic algorithm prior to our algorithm, which is more suitable for a global scan of the parameter space [17]. A more detailed discussion of the simplex initialization is given below in section 4.

3.3. Assumed modal functions for fiber characteristics

The algorithm, as it has been described so far, can be applied with any kind of diffusion or attenuation function and a generally unlimited number of free parameters, on condition that convergence takes place at physically reasonable values. However, increasing the number of free parameters drastically reduces robustness of the method, when no additional restrictions or more sampling fiber lengths L_C, L_D,... are introduced.

For the results discussed later in this paper we applied following functions for modal diffusion and attenuation:

Besides the physical restriction of a(θ) not to be smaller than zero at any θ, we set d(θ) > 0 and d_q > 0 to ensure positive diffusion, decreasing with propagation angle [19]. θ_c is the inner critical angle of the fiber given by the fiber’s numerical aperture (NA). a₀ depends on a₁ and a₂ as well as on the independently measured fiber attenuation γ and is given by

Although there is no necessity of choosing it that way, we found that the algorithm is more robust for L_A < L_B, especially for underfilled launch and fiber lengths shorter than the coupling length. We attribute this to the fact that, firstly, attenuation tends to increase with propagation angle and that, secondly, the amount of higher order modes propagating at large angles scales with fiber length for underfilled launch and fiber lengths shorter than the coupling length. Thus, the impact of mode-dependent attenuation becomes more pronounced as the fiber is elongated.

4. Validation of approach

4.1. Setup and initialization of the algorithm

In order to validate our approach we apply it to frequently used SI-POFs with known diffusion and attenuation characteristics: ESKA-PREMIER GH4001 (GH) from Mitsubishi and PGU-FB1000 (PGU) from Toray, both with NA = 0.5 and n = 1.49. The diffusion and attenuation characteristics of both fibers have been experimentally determined in [9]. With these fiber properties farfield patterns and frequency responses have been computed and close match to experimental data has been demonstrated [9, 11]. Owing to the direct correlation between frequency and impulse response we argue that computed impulse responses based on these fiber properties must be very close to real impulse responses of the same fiber under similar launching conditions. Therefore, we use the fiber properties presented in [9] for computation of reference impulse responses, H_ref(L_A, t) and H_ref(L_B, t).

As initial field distribution we set

As mentioned above simplex initialization is crucial for avoiding contraction at local minima in the DS algorithms. In Fig. 2 the simplex initializations we applied to all but the first main loop roundtrip are shown. For the attenuation DS algorithm the initial simplex is an equilateral triangle with the best data set of the previous main loop, $S_{f,\left(i-1\right)}^a$, in its centroid. Size and position of the simplex are given by the vertex $S_{1,i}^a$, which is obtained by mirroring $S_{f,\left(i-1\right)}^a$ at the origin of the coordinate system. For the diffusion DS algorithm the parameters of the previously best set are multiplied and/or divided by the variable ε, which starts at a value of 20 and approaches two as the algorithm advances. Thus, four possible simplex vertices surrounding $S_{f,\left(i-1\right)}^d$, the previously best diffusion parameter set, are computed. The points indicated as $S_{2,i}^d$ and $S_{3,i}^d$ as well as one of the $S_{1,i}^d$ points, which is alternated at each main loop iteration, form the initial right angle simplex of the algorithm. The applied simplex initializations serve both: the likeliness of convergence to local minima is reduced by varying the appearance of the initial simplex and inheritance of the previously best parameter sets is ensured.

 

Fig. 2 Simplex initialization for attenuation (a) and diffusion (b) DS algorithms

Download Full Size | PDF

The initial simplex vertices for the first main loop roundtrip are set as $S_{1,1}^d$ (d₀ = 1, d_q = 1), $S_{2,1}^d\left(0,1\right)$ and $S_{3,1}^d\left(1,0\right)$ for the diffusion function DS algorithm and $S_{1,1}^a$ (a₁ = −100, a₂ = −100), $S_{2,1}^a$ (−100, 100) and $S_{3,1}^a$ (100, −100) for the attenuation function DS algorithm. These initial data sets are exceeding reasonable values on purpose to demonstrate robustness of our method. Further, we treat exceeded boundary conditions as infinite RMSD. This ensures substitution of these sets as the algorithm advances and retains the general structure of the DS scheme. In case that none of the simplex vertices lies within the defined parameter space, a new simplex is computed from randomized points within the prior simplex.

Discretized step sizes are Δz = 1 mm and Δθ = 5 mrad. For reduction of computational effort we roughly optimized the fiber characteristics at the first six roundtrips (128 sampling points) and increased accuracy to 2048 sampling points afterwards.

4.2. Simulation results and discussion

As the algorithm advances the modeled fiber characteristics and computed impulse responses improve. As an example, the evolution of PGU fiber characteristics and corresponding impulse responses at 100 m fiber length are shown in Fig. 3. At each iteration step the pulse shape gets closer to the reference shape and close match is reached in the end. Major fitting is accomplished within the first iteration steps, whereas conformity improves gradually in the final iteration steps. Likewise, simulated fiber characteristics approach the reference data. We attribute the relatively small variations in attenuation function compared to the vast variations in diffusion function to the constraint given to the attenuation function by Eq. (11).

 

Fig. 3 Evolution of impulse responses of 100 m PGU fiber (left plot), modal attenuation (top right) and modal diffusion (bottom right) over iteration steps of the proposed method. Simulated data is shown as red lines, reference data from [9] in blue lines.

Download Full Size | PDF

The characteristics of the GH fiber evolve in similar way. The computed parameters of both fibers after ten main loop iterations are given in Table 1, the corresponding functions are displayed as solid red lines in Fig. 4. The iteration procedure yields reasonably close match of simulation and reference characteristics for both fibers, even so GH characteristics are reproduced to a better degree than PGU characteristics. The least squares fits of Eqs. (9)–(11) to the reference fiber characteristics reveal that the remaining deviations of simulated to reference data are mainly attributed to intrinsic differences of modal functions assumed by Mateo et al. and by us. In fact, RMSDs of the least squares fits to the reference characteristics are only marginally lower (by a factor of about 3 or less) than those of the simulated modal functions. We take the reasonably close match of simulated and reference fiber properties as evidence for functionality of our approach that reconstructs these characteristics solely from iterative adaption of fiber impulse responses. Furthermore, significant improvement of conformity for refined modal functions can be expected from the presented results.

 

Fig. 4 Diffusion (left) and attenuation (right) functions for GH (top row) and PGU fiber (bottom row). Blue lines represent reference data from [9], simulation results are shown as solid red lines. Dashed red lines are least squares fits of Eqs. (9)–(11) to the reference characteristics for θ ≤ θ_c.

Download Full Size | PDF

Table 1. Fiber parameters for modal diffusion and attenuation

View Table

Note, that the simulated modal functions are obviously no perfect imitations of the original functions, which certainly results in deviations of the simulated spatio-temporal responses to the reference data. However, in the subsequent discussion it will be demonstrated that the accuracy of the simulation results is precise enough to predict fiber responses in the temporal domain.

With the results given in Table 1 impulse and frequency responses have been computed and are exemplarily shown for different lengths of the PGU fiber in Fig. 5. Both, the impulse (left) and frequency responses (right), are displayed for fiber lengths from 50 m to 700 m (red lines) and again compared to results obtained from reference data (blue lines) [9].

 

Fig. 5 Impulse (left) and frequency responses (right) at varying lengths of the PGU fiber. Blue lines represent reference data from [9], simulation results are shown as red lines.

Download Full Size | PDF

Except for a slight overestimation of the pulse’s trailing edge, which becomes apparent for large L because of the Fig. scaling, close match of simulation to reference impulse responses for the entire fiber length range is achieved. We attribute this overestimation to differences between the simulated modal functions and the reference data as shown in Fig. 4. For propagation angles larger 4°, attenuation characteristics of the PGU fiber are reconstructed reasonably well. However, attenuation at narrow angles is overrated significantly. In our simulation this deviation is compensated for by the net diffusion, that is, in contrast to the reference net diffusion, continuously directed towards lower angles. Consequently, signal power at higher angles is lowered, which is in turn compensated for by a slightly underestimated attenuation for 4° <θ < 12°.

As apparent from the impulse responses in Fig. 5, this interactive compensation of discrepancies results in wide conformity of impulse responses. However, its limitations are demonstrated by the deviations at the pulse’s trailing edge, that corresponds to large average propagation angles. At these angles, the reference net diffusion is directed towards smaller angles at a rate twice as high as for the simulated net diffusion. Therefore, increased attenuation of the simulated signal would be necessary to compensate for the differences in diffusion. However, within the constraints given by the attenuation function, this cannot be accomplished without a global reduction of impulse response conformity. Consequently, root-mean-square (RMS) pulse durations of the simulated impulse responses are marginally longer than those of the reference impulse responses (Fig. 6(b)).

 

Fig. 6 Bandwidths and RMS pulse durations over lengths of GH (a) and PGU (b) fiber. Reference data from [9] is given by blue symbols, simulation results are given in red (dashed) lines. Black symbols denote measured bandwidths for two samples of each fiber (filled circles and open squares, respectively) from [11]. The arrows indicate which vertical axis the data refers to.

Download Full Size | PDF

As well as the impulse responses, the simulated frequency responses match the reference frequency responses very closely. The remaining deviations are especially apparent for short fiber lengths. The reference traces have slightly less amplitudes at low frequencies and significantly more contribution of high frequencies than the simulated traces at these fiber lengths. The differences result from the deviations in the impulse response shapes, which are narrower and steeper for the reference data. As the fiber is elongated, the improvement of impulse response simulation is apparent from the increased conformity of the frequency responses. The left deviation, namely that the reference pulse’s trailing edge falls steeper than the simulated one, is displayed by a slightly higher contribution of the frequency components at the falling edge of the frequency response. The discussed differences in both regimes, for short and for long fibers, result in slightly overestimated, or respectively underestimated, fiber bandwidths as shown in Fig. 6(b). Both simulated and reference bandwidths capture the general tendency of the experimental results from [11] denoted as black symbols. The improved conformity of computed and measured bandwidths for longer fiber lengths is consistent with the finding of Mateo et al. that the signal suffers from strong initial diffusion at the fiber input [11]. As in [11], we also achieve better match of computed to experimental bandwidths for the PGU fiber than for the GH fiber.

Despite the discussed deviations, the simulation results for the PGU SI-POF coincide very well with the reference data. Thus, the results enable estimation of further MM fiber characteristics as e.g. the coupling length L_c indicating equilibrium mode distribution (EMD) at which pulse broadening rate changes from ∝ L to $\propto \sqrt{L}$[3]. From the results displayed in Fig. 6(b) we obtain L_c ≈ 40 m, which is close to recently published coupling lengths of 38 m or 46 m, obtained from numerically computed spatial field distributions [10]. Similarly, slopes of bandwidth-length dependence have been estimated by a two-line fit to the corresponding curves. For the simulation data we obtain 1.30 and 0.58. These results match data from [11] quite well, who obtained 1.3 and 0.8 from experimental data up to L ≈ 100 m and similar launching conditions.

For the GH fiber we obtained L_c ≈ 70 m as well as 1.32 and 0.65 as two-line fit bandwidth slopes. Again, the results are in good agreement to previously published results that are 75 m or 80 m for coupling length L_c[10] and 1.2 for bandwidth slope of GH obtained from experimental data up to L ≈ 100 [11].

Regarding accuracy of our algorithm we notice less deviation of simulation to reference data for the GH fiber compared to the PGU fiber. We attribute this to the notably better reconstruction of the reference fiber characteristics by our algorithm, as apparent from Fig. 4. We argue that the main reason for lower accuracy at PGU modeling is given by the comparably vast variations of the reference attenuation function at narrow propagation angles, which cannot be reconstructed by the modal attenuation function we assumed in Eq. (10). For the GH fiber, the reference attenuation fluctuates about six times less as for the PGU fiber, so more accurate attenuation modeling is possible. In this way, the diffusion function might be modeled more precisely as well, since there is no need to compensate for intrinsic deviations of simulated and reference attenuation functions and an overall better simulation result is obtained.

Comparing the results of simulation and reference demonstrates both the tight relationship between diffusion and attenuation and the successful reconstruction of these fiber properties by the algorithm. Although very different modal dependencies have been assumed for the fiber properties, good agreement of impulse and frequency responses as well as pulse durations and bandwidths is achieved. Furthermore, compliance of simulated results and experimental results from [11] supports the reasonable assumption that the computed reference impulse responses must correspond closely to real impulse responses and, consequently, that applying those is an applicable approach for validation of our method.

5. Conclusion

We have presented a computational method for estimation of mode coupling and attenuation characteristics of step-index plastic optical fibers that can be used for prediction of impulse and frequency responses, coupling lengths and bandwidths. In contrast to previously published methods, we concluded global fiber characteristics solely from fiber impulse responses at two different fiber lengths.

We validated our approach against reference fibers and their characteristics, published by other authors, and demonstrated its applicability for frequency and impulse response prediction. Further, we computed coupling lengths and bandwidth slopes of the fibers under test and achieved good agreement to recently published data for identical fibers. We discussed the intrinsic limitations of our method and demonstrated that our method yields reliable pulse duration and bandwidth results even when significant differences to reference fiber characteristics are present.

Although reliable measurement of impulse responses is challenging, especially for short fiber lengths, our method might significantly reduce experimental effort, since no spatial information concerning the field emitted by the fiber is needed. Measuring temporal fiber responses to impulses or step function signals at two fiber lengths suffices to estimate the fiber’s characteristic properties. Furthermore, the characteristic functions underlying our algorithm can be varied and/or expanded to match other, potentially more sophisticated, modal dependencies. However, it has to be mentioned that the proposed method is computationally more intensive than the determination of fiber properties from farfield patterns.