1. Introduction

Recent developments in astrophotonics [1, 2] demonstrate a clear need for fast and robust IS algorithms for the design of ultrabroadband, strong FBG-based filters. For these filters, typically iterative (so-called “mixed” scattering schemes (MSS) [3, 4]) are used. These schemes require extreme accuracy and stability of both inverse scattering (IS) and direct scattering (DS) algorithms, which are used for every iteration step. Although few issues with DS methods are reported, it is well known that essentially all IS schemes have convergence limitations, particularly if used to design extremely strong gratings. In this paper, we discuss which of the IS layer peeling methods (LPM) are best suited to the gratings described in Ref. [1, 2]. Here, the goal is to achieve FBGs with extremely large bandwidths (up to 262144 grid points needed), closely located narrow notches and ultrastrong suppression notches. The successful devices designed by us constitute the most complex optical filters ever attempted.

The layer peeling method historically was the first efficient IS procedure adopted for FBG design. Different versions of the method were reported by [5, 6, 7]. The method has a clear physical interpretation of the reflected signal as a superposition of impulse responses from different uniform layers or point reflectors placed along the grating. Each thin layer has small reflectivity and can be taken into account within the so-called first Born approximation. Because of high efficiency (of order N ² operations [6, 7], where N is the number of discretisation points in spectral and spatial domains), this method is now used.

A major issue relating to any IS algorithm is the potential exponential decay of accuracy along the grating because of error accumulation during the reconstruction process [8]. Typically, however, the error accumulation is non-exponential especially for chirped gratings [8]. It is now known [9, 10] that Feced’s method demonstrates less error accumulation in comparison to Skaar’s algorithm [7] and especially Poladian’s approach [6]. Note, however, that as presented in the original work [5], Feced’s method requires ln(N)*N ² operations to run and thus is substantially slower compared to Skaar’s approach. Modifications to Feced’s method suggested in [10] increase the speed of Feced’s LPM up to the speed of Skaar’s method without compromising error accumulation levels.

The other approach to LPM has been implemented in the time domain [7, 8, 11, 12]. The method claims better accuracy than both Feced’s and Skaar’s LPM, but the examples presented in these works, although highly informative, do not fully demonstrate the robustness limits of the algorithm. In fact, we were able to reproduce most of these examples with the spectral-domain LPM IS algorithm of [5, 10]. In this paper, we examine the most efficient version of time domain-based LPM method in detail and compare its speed and robustness with the frequency based LPM approaches of [5, 7].

Encouraging results at high reflectance are demonstrated by a combination of layer peeling and iteration. In the integral layer peeling method proposed by Rosenthal and Horowitz [9], the grating is divided into thin layers, but these layers are not assumed to have a uniform profile. The profile of each layer is found by an iterative solution of Gel’fand-Levitan-Marchenko (GLM) equations. The required number of operations in combined methods is claimed to be of the order of N ². In reality, such efficiency is only possible for really simple FBGs with smooth variation of parameters. For complex multi-channel gratings requiring fast variation of amplitude and phase profiles the actual speed is estimated to be m * N ² were m is the number of required iterations (typically m > 1) resulting in substantial speed reduction, compared to Skaar’s method [7] or the modified Feced’s [5, 10] method.

We note that the authors of [9] mention the possibility for each of their peeled layers to have a single point only, thus eliminating the need for iterations. The method of [9] thus becomes very similar to the frequency-domain LPM IS version of [7] which, in addition, includes a regularisation technique to improve stability (cf. Eq. (6–8) of [7] with Eq. (18–19) of [9]). However, this modification slows the algorithm down in comparison to the original algorithm of [7].

Another recent attempt at a numerical solution of the GLM equations was made in Ref. [13]. The GLM equations were reduced to a matrix form. The inversion of the corresponding matrix and reconstruction of the grating was achieved by a bordering procedure and Cholesky decomposition. This approach takes of order N ³ operations which makes it much too expensive for efficient grating designs. Very recently, the same authors were able to significantly increase the efficiency of their original proposal by discovering and exploiting diagonal symmetry of the GLM matrix [14]. The new method, called Toeplitz Inner Bordering, is claimed to be highly efficient and requires only N ² operations.

Table 1. A summary of popular IS algorithms including references and information about the required number of operations M needed for one IS cycle. Only the fastest methods M ≤ N ² are potentially suitable for the iterative mixed scattering scheme as in Ref. [2].

View Table | View all tables in this article

Table 1 provides a brief summary of currently available IS algorithms. Here, we limit our attention only to LPM IS methods that require no more than N ² operations for implementation. Thus we consider only algorithms II, III and V from Table 1.

2. Model and analysis

The coupled-mode system of equations describing light propagation in FBGs is

where E_f and E_b are the amplitudes of the forward and backward propagating fields, respectively, δ is the normalized frequency detuning from the central Bragg reflection frequency, z is a local distance along FBG, and q(z) = κ(z)exp[ϕ(z)], where κ = ∣q∣ and ϕ = Arg(q), is a spatial profile of the FBG coupling coefficient. For a reciprocal and lossless (see, e.g., [15] for definitions) FBG of length L we may find complex reflection and transmission coefficients from a transfer matrix T̂, which relates field values at the grating ends

and is given by

where r ₁ is the reflection measured at z = 0, r ₂ is the reflection measured at z = L, and t is the complex transmission coefficient (which is independent of the light propagation direction for reciprocal FBGs). Conventional reflection R_i, group delay τ_i, and chromatic dispersion D ⁽ⁱ⁾ (i = 1,2 depending on the grating side) are related to the complex reflection coefficients as R_i = ∣r_i∲², τ_i = - ∂Arg(r_i)/∂ω) and D ⁽ⁱ⁾ = -(2πc/λ ₀ ²)∂ ²Arg(r_i)/∂ω ², where ω is the light wave frequency, λ ₀ is the FBG central wavelength and c is the speed of light in vacuum. Similar relations exist for transmission T, group delay in transmission τ_t and the complex transmission coefficient t. The calculation of the reflection data r(τ) from the grating profile q(z) is usually referred to as the direct scattering transform (DS), whereas the inverse operation of recovering q(z) profile from the spectral data r(δ) is called the inverse scattering transform (IS). In this paper, we mainly concentrate on the properties of IS algorithms and will comment on DS algorithms only briefly.

In addition to comparing speeds of different IS algorithms we will also analyse their relative robustness, i.e. comparative ability to handle “tough” designs: designs where channel strengths are very high (T < - 100dB) and where spectral channels are located in close proximity to each other. It is typically argued that poor robustness is the result of error accumulation which has three primary sources:

It is still has to be seen which error is most important. We should note, however, that all IS methods are based on finite, rather than infinite, intervals which leaves type 1 errors outside the scope of our analysis. In addition, only algorithm VII allows the FBG envelope segments to have a trapezoidal shape (with a version based on rectangular segments also available), while for methods II, III and V, the trapezoidal versions are either hard or impossible to implement.

Finally, algorithms II and III require a substantial amount of sine and cosine calculations, whereas individual steps of algorithm V are free of them. Thus, the difference in performance of algorithms II and III versus algorithm V will be a very good proxy for the importance of the error of type 3, whereas the difference in performance of algorithms II, III and V versus the trapezoid version of the algorithm VII will be a proxy for the importance of the type 2 error. It may be even better simply to compare the performance of a rectangular segment-based version of algorithm VII with its trapezoid segment-based counterpart.

Before going into the technical details of our analysis, we should briefly described implementation details of algorithms II, III and V. Implementation of algorithm II closely follows the discrete layer peeling (DLP) version of Ref. [7]. Note that Ref. [7] also outlines the corresponding continuous layer peeling (CLP) IS algorithm, which was originally suggested in Ref. [6]. However, CLP version was found to be slower than DLP for the same required accuracy of computation [7] - a result confirmed by our own numerical modelling as well. Below we only talk about DLP implementation of [7].

Algorithm III is based on the ideas of Ref. [5, 10] for which Ref. [5] outlines the general description and Ref. [10] provides speed-enhancing changes. Here we should note one important point: we found that the Feced algorithm is rather unstable if implemented exactly as in [5]. To fix the algorithm convergence, we had to change Eq. (16) of [5] to a slightly different form:

instead of the original

We suspect that this misprint is one of the main reasons why algorithm III has been criticised in the literature (see, for example, a discussion of the Gibbs phenomenon in [12]; we indeed found that Gibbs-like effect is present if Eq. (5) is used, but is eliminated by changing to Eq. (4) in [5]). Paradoxically, it is quite possible that the existence of this misprint has contributed positively to the recent rapid progress in the field of IS algorithm development: difficulties with implementation of the stable version of the spectral domain-based algorithm III as described in [5] may have encouraged a few researchers to develop their own more robust versions of the IS algorithm.

Algorithm V is implemented in line with [8], which is a time-domain realisation of algorithm II. A slightly different implementation of [12] is a time-domain analog of algorithm III. It has a very similar accuracy to algorithm V (described below) and slightly higher complexity resulting in 10% – 15% slower speed. However, such a small difference in algorithm performance may be attributed to a sub-optimal optimisation of our implementation of the method of [12], because we put a relatively larger emphasis on optimising the version of [8]. We describe our implementation of [8] in detail below and compare it to the description on pages 2230–2231 of [8], making small changes to improve efficiency.

Algorithm V. At step 1 of algorithm V, we initialise the forward- and backward-propagating fields at the first layer (j = 1):

where h is the discrete Fourier transform of target reflection spectrum r ₁. This step is essentially identical to step 1 of [8] apart from index shifts by unity so that the first FBG layer now corresponds to j = 1.

At step 2, we compute the reflection coefficient ρ_j = v ^(j) ₁/u ^(j) ₁ (identical to step 2 of [8]), apart from the index shift: j = 1, rather than 0 in [8], for the first layer.

At step 3, we combine propagation step 3 and time-shifting step 4 of [8] and also make some efficiency increasing modifications as follows. First, we calculate u ^(j+1) ₁ = u ^(j) ₁ -v ^(j) ₁ ρ_j ^*. Then we calculate all other required u and v elements in a cycle:

where index k is incremented by 1 from the initial value of k = 1 until the final value of k = N - j inclusive, where j is the index of the currently “peeled” FBG layer. Note that with j increasing, the cycle is getting shorter reaching a single point at j = N - 1. No operations are performed at j = N. Note that the form and the order of Eqs. (7) is ideally suited for the numerical implementation of the cycle and do not require any temporary variable use: one only needs to perform v_k ← v _k+1 - ρ_ju _k+1 and then u _k+1 ← u _k+1 - ρ_j ^* v _k+1 operations.

At step 4, we check the j count. If j < N then we increase j by 1, go to step 2 and continue. If j = N, then we stop. Once all coefficients ρ_j are calculated, complex coupling coefficients q_j (for FBG layers in z intervals of (j - 1)Δz ≤ z < jΔz, where Δz = L/N) can be calculated as q_j = - (1/Δz)arctanh(∣ρ_j∣)ρ_j ^*/∣ρ_j∣, where j = 1,2,3,…N.

In implementing of all of the above-mentioned IS algorithms, we only aim to obtain FBG profiles as an output without simultaneously calculating the corresponding final spectra, which may be slightly different from target spectra for complex designs. This is exactly what is needed for MSS of [3,4], where FBG profiles are obtained, then modified and only then used to obtain new spectral data. To complete an MSS cycle, efficient DS algorithms are also needed but these are outside the scope of the present paper. However, we note that each of the discussed LPM IS algorithms can be straightforwardly transformed into a corresponding DS algorithm by backward propagating E_b(δ,z = L) = 1.0 and E_f(δ,z = L) = 0 fields from z = L to z = 0 in spectral or time domain representation and calculating the reflection r ₁(δ) = E_b(δ,z = 0)/E_f(δ,z = 0).

Before starting the speed and robustness analysis, we adopt a sample FBG spectrum which will be used for testing FBG grating profile reconstructions by different IS algorithms. We opt to use a 7-channel spectrum with channel spacing of Δλ_CS = 0.8nm. Apart from the central wavelength position [taken as λ_i = 1550 + Δλ_CS * (i - 4) nm, where i= 1,2,3… 7], each individual channel has identical spectral characteristics: spectral half-width of Δλ = 0.33nm and second order dispersion value of D ₂ = - 500ps/nm. The reflection profile has the following shape:

where Δλ = 0.33nm, λ_i are channel centres (i = 1,2,3…7), T_min is the minimum transmission value (transmission in channel centres) and a high value of n is chosen: n = 30. This profile is equivalent to a high-order Butterworth profile that is essentially ideal for astronomical filters [16]. It is reasonably hard to reconstruct for T_min → 0 especially for multichannel designs with bandwidth utilisation Q = 2Δλ/Δλ_CS values exceeding 80%. Finally, we dephase spectral channels with respect to each other to reduce the peak value of reconstructed designs (see [4] for details of this technique) with the dephasing angle values of ϕ ₀ ⁽ⁱ⁾ = (0, 0.063300π, 0.768574π, 0.782320π, 0.101907π, 0.729966π, 0).

3. Results

For speed comparisons, we use our sample design with T = -10 dB and construct FBG designs over a 40cm length, although the length where q(z) value is reasonably strong (say, κ = ∣q∣ > 0.001) is less than 20cm. Extra length improves LPM IS convergence. Note that for all three analysed algorithms the complexity of design (e.g. change in channel number or channel strength) has a limited effect on algorithm speed: for the same algorithm and the same number of grid points variations of speed do not exceed 5% for any well-converging design. We used a Windows Vista-based laptop PC with Intel Core 2 Duo T7300 2 GHz CPU and 2 GB RAM. The results of speed analysis are presented in Table 2. Note that an absolute speed value is significantly dependent on detailed software and hardware configuration and our results should only be used for comparison of relative trends.

Table 2 shows that algorithms II and III have similar speeds, whereas algorithm V outperforms both of them by a factor 10 to 20 for different numbers of grid points. We checked that this trend continues for higher number of grid points as well. For example, for 65536 grid points, the run time of algorithm V is 35 seconds versus 680 seconds for algorithm III.

Table 2. The results for comparative speed analysis of algorithms II, III and V. All data represent CPU clock time in seconds for different algorithms and different number of grid points.

View Table | View all tables in this article

Now we proceed to convergence analysis which is critical for MSS applications. We use the same test design, fix the number of grid points to 8192 and increase channel strength for all channels until the algorithm exhibits convergence problems. We vary T_min strength by 1dB increments. For simplicity we limit all convergence problems to two types: (i) T ripple issue, where transmission ripple for any of channels starts to exceed 0.1dB level and (ii) “tail” issue, where the tail of FBG grating design starts to display small artificial ripples (we take κ= 0.001 as an acceptable limit), which rapidly grow with further increase of the FBG strength. See Fig. 1 for an example of the “tail ripple” effect.

 

Fig. 1. An example of a small ripple development in the tail of a strong FBG design (see around z = 30 region).

Download Full Size | PDF

Our main algorithm convergence results can be summarised as follows: algorithm II is the first to fail (due to T ripple development) at T_min = -31dB. Algorithm III fails next (also due to T ripple development) at T_min = -40dB. Algorithm V does not fail until remarkably significant levels of about T_min = -102dB. At these levels a small artificial ripple starts to develop, but no ripple in T or group delay characteristics is present until at least T_min ≈ -120dB. The corresponding critical FBG designs and their spectral characteristics are presented in Figs. 2, 3 and 4.

Our results quite clearly demonstrate that the time domain LPM IS algorithm substantially outperforms spectral-based LPM algorithms both in speed and robustness producing truly outstanding results and eliminating concerns of [9] about the intrinsic issues with LPM IS algorithms. Since the only major difference between algorithms II and III and algorithm V is the lack of sine and cosine calculations in the latter method, we conclude that the evaluation of trigonometric functions is the primary cause of IS performance problems and should be avoided if possible. This conclusion is also supported by our extra numerical experiment where we changed double precision to single precision in a few key subroutines of algorithms II, III, and V. Interestingly no noticeable deterioration of algorithm V performance is observed, but both algorithms II and III were seriously affected and displayed a large error increase. This result, in turn, may be attributed to “noise-magnifying” properties of trigonometric functions.

 

Fig. 2. The maximum high-quality convergence for the method of [7]: Top value for transmission reduction is limited to -31dB (for accepted T ripple of ±0.1 dB).

Download Full Size | PDF

 

Fig. 3. The maximum high-quality convergence for the method of [5] with modifications of [10]: Top value for transmission reduction is limited to -40dB (for accepted T ripple of ±0.1 dB).

Download Full Size | PDF

 

Fig. 4. The maximum high-quality convergence for the method of [11, 12]: Top value for transmission reduction is -102dB. For lower values of T_min the tail ripple starts to exceed 0.001 level. It is important to note that this tail ripple does note manifest itself in substantial T and group delay ripples for T_min > -120dB.

Download Full Size | PDF

 

Fig. 5. An example of 147-channel FBG filter design with 400nm bandwidth (obtained with algorithm V). Algorithms II and III failed to converge on this design.

Download Full Size | PDF

Similar performance results can be obtained for the corresponding DS algorithms, i.e. modified algorithms II, III and V for the direct propagation of E_b(δ, z) and E_f (δ, z) fields for a given FBG profile q(z). We implemented these corresponding DS algorithms primarily for IS accuracy check purposes: for a -10dB strong FBG sample 7-channel design q(z) concurrent application of DS and IS algorithms results in ∣q_initial(z)- q_final(z)∣ difference of less or equal than 10^-12 for all 0 < z < L, when calculations were performed in double precision. Similar to the IS results, a DS algorithm based on a time-domain approach is at least one order of magnitude faster than its spectral-based counterparts (none of the algorithms have convergence issues). However, spectral-based methods are a lot more flexible allowing spectral calculations for a selected range of frequencies, whereas the time-domain approach always requires a calculation of the complete frequency range. As a result the time-domain DS algorithm, in contrast to spectral-based methods, cannot be parallelised straightforwardly.

Finally, we note that it is preferable, but not essential, to have IS and DS methods of the same algorithm type in a mixed scattering scheme (MSS). For example, although we now typically use matching-type IS and DS in our MSS implementations (both based on algorithm V), sometimes the use of type V IS together with type II or type III DS algorithms is preferred due to the spectral range flexibility and parallelisation reasons described above.

4. Discussion and conclusion

We compared several popular LPM algorithms and analysed their relative speed and robustness. In addition we suggested a few minor efficiency improving modifications of these methods. Unexpectedly we found a previously unreported result that the time domain-based algorithm of [11, 12] (algorithm V) is an order of magnitude faster and significantly more robust than any of analysed spectral domain-based LPM IS algorithms. We attribute the superior speed and robustness to the small amount of trigonometric function evaluations used in the algorithm of [11, 12] compared to spectral domain-based LPM IS algorithms.

As a result we identified the algorithm which is well suited for MSS applications. In fact, only using algorithm V as an element of MSS scheme of [2] we were able to design a strong FBG-based filter covering the bandwidth of 400nm and containing 147 narrow notches with notch T_min ranging from -10 to -40dB without any convergence problems (see Fig. 5). We attribute this success to superior robustness properties of the underlying IS algorithm (algorithm V).

In the future we plan to extend our analysis to include IS algorithms of [9] and [14]. Although our preliminary estimates indicate that these methods are very unlikely to match the speed of the time domain-based LPM IS algorithm V (see §1), they still may have better robustness properties. If this turns out to be the case, then these IS methods may be useful for MSS-based construction of grating designs, where even algorithm V has convergence issues.