The discrete layer-peeling algorithm (DLPA) requires to discretize the continuous medium into discrete reflectors to synthesize nonuniform fiber Bragg gratings (FBG), and the discretization step of this discrete model should be sufficiently small for synthesis with high accuracy. However, the discretization step cannot be made arbitrarily small to decrease the discretization error, because the number of multiplications needed with the DLPA is proportional to the inverse square of the layer thickness. We propose a numerically extrapolated time domain DLPA (ETDLPA) to resolve this tradeoff between the numerical accuracy and the computational complexity. The accuracy of the proposed ETDLPA is higher than the conventional time domain DLPA (TDLPA) by an order of magnitude or more, with little computational overhead. To be specific, the computational efficiency of the ETDLPA is achieved through numerical extrapolation, and each addition of the extrapolation depth improves the order of accuracy by one. Therefore, the ETDLPA provides us with computationally more efficient and accurate methodology for the nonuniform FBG synthesis than the TDLPA.
© 2011 OSA
For certain applications of fiber Bragg gratings (FBG), it is important to design the grating profile that produces the desired spectrum. This fiber Bragg grating (FBG) synthesis problem can be solved by one of the inverse scattering techniques, for example, the Gel’fand-Levitan-Marchenko (GLM) method  or the layer peeling approaches [2–4]. The GLM method can synthesize the FBG having a rational spectrum, but cannot guarantee the accuracy of the result for a general spectrum shape, while the approaches based on numerical optimization schemes require intensive computational resources [5–7].
On the other hand, the layer peeling method identifies the medium layer-by-layer using the causality condition of the counter-propagating waves. Since the pioneering works of  and , the discrete layer-peeling algorithm (DLPA) has been regarded as one of the most efficient methods for the synthesis of nonuniform FBG [8–15]. In the usual FBG synthesis problem, the time domain DLPA (TDLPA) has better accuracy than the frequency domain DLPA (FDLPA) , although the TDLPA shows degraded performance for a strong uniform grating  because the TDLPA is sensitive to the roundoff errors and the error accumulation of the TDLPA is maximized when the sign of the coupling profile is not changed . The disadvantage of the TDLPA for a strong uniform grating can be resolved using the integral layer-peeling method [8, 9]. The accuracy of the FDLPA can be made comparable to the TDLPA by increasing the spectral resolution , which increases the computational complexity of the FDLPA. Furthermore, the TDLPA not only has lower computational complexity than the FDLPA, but are appropriate for the parallel implementation because the inner product is not involved in the TDLPA [17–19].
The DLPA computes the coupling coefficient using the equispaced discrete model approximation of the piecewise uniform model for synthesis of nonuniform FBG [2, 3], and the accuracy and the efficiency of the piecewise uniform model for FBG is verified both experimentally and numerically . Similarly to the usual numerical methods, the numerical discretization error of the DLPA is inevitable and the discretization step, i.e., the layer thickness of the DLPA, should be sufficiently small to get high accuracy. However, if the layer thickness of the DLPA is decreased to improve the accuracy, the total number of layers of the DLPA becomes larger. The computational complexity of the DLPA increases proportionally to the inverse square of the layer thickness because the involved multiplication count of the DLPA increases proportionally to the square of the total number of layers. Therefore, the accuracy of the previously known DLPA can be improved only with the greatly increased computational cost.
To overcome this tradeoff between the accuracy and the computational complexity of the DLPA for synthesis of nonuniform FBG, we propose a numerically extrapolated TDLPA (ETDLPA). The proposed ETDLPA employs a numerical extrapolation scheme  to improve the accuracy of the conventional TDLPA without increasing the computational complexity significantly. Similarly to the conventional numerical extrapolation scheme, the ETDLPA combines the results of the TDLPA hierarchically starting from the coarsest discretization, and the order of accuracy of the TDLPA increases as much as that of the extrapolation depth.
However, the extrapolation procedure of the ETDLPA differs from the conventional numerical extrapolation scheme in two ways. First, the extrapolation procedure of the ETDLPA is applied to the vector quantity, unlike the conventional extrapolation method which operates on the scalar quantity. The output of the ETDLPA is a set of reflection coefficients that is a vector quantity. Second, the computation procedure of the ETDLPA is not memoryless because the reflection coefficient at the current position is dependent on the previous reflection coefficients, which is not the case with the conventional extrapolation scheme. In the ETDLPA, the extrapolated reflection coefficients for the whole grating positions cannot be obtained by a single execution of the conventional extrapolation procedure because of these numerically distinctive features of the ETDLPA. In this paper, the overall procedure of the ETDLPA is completed by the repetition, with suitably shifted scattering data, of the basic extrapolation procedure as many times as determined by the involved layer thickness values.
In spite of the increased accuracy order, the computational complexity of the ETDLPA, which employs practically meaningful extrapolation depth, is usually less than three times of that of the DLPA.
The synthesis of nonuniform FBG is the one-dimensinoal inverse scattering problem represented by the coupled mode equations in the space-frequency domain 
This synthesis problem can be solved by the TDLPA based on the discretization of Eq. (1) . Assume that nonuniform FBG is discretized as discrete reflectors with layer thickness h. The error formula of the TDLPA with the layer thickness h is given by the power series of h,Eq. (2) is given in appendix A. It is clear from Eq. (2) that the accuracy of the TDLPA is O(h). The key feature of the ETDLPA is to eliminate the remained error terms in Eq. (2) by combining the results of the TDLPA for various layer thickness values through a numerical extrapolation procedure. For example, the result of the TDLPA with the layer thickness 2h can be written from Eq. (2) as follows, Equation (4) is only the first stage of the extrapolation procedure of the ETDLPA, in which the first degree term a 1 h of Eq. (2) is removed in Eq. (4), and thus, the accuracy order is increased from O(h) to O(h 2). This observation shows that additional information for the coupling coefficient computed by the TDLPA with under-sampled scattering data can improve the accuracy of the TDLPA. Note that the scattering data with the layer thickness 2h is just two times under-sampled data of the scattering data with the layer thickness h. In the TDLPA, the accuracy O(h 2) is obtained with the layer thickness of h 2. If the number of layers with the layer thickness h is N, the number of layers with the layer thickness h 2 is N/h. It is easy to become aware that the computational complexity of Eq. (4) can be remarkably lower than that of q̂(z,h 2) to achieve the accuracy O(h 2) if h is sufficiently small. Note that the value of a l does not need for the extrapolation. The computational efficiency of the ETDLPA will be detailed at the end of this section.
Now we generalize the concept of Eq. (4) to get higher accuracy. The P-stage extrapolation procedure of the ETDLPA can be formulated using q̂(z,h),q̂(z,2h),··· ,q̂(z,(P+1)h) to obtain the accuracy O(h P +1). From Eq. (2), the result of the TDLPA with the layer thickness mh can be written byEq. (2) by combining q̂(z,h),q̂(z,2h),··· ,q̂(z,(P + 1)h), we get Eq. (5) into the left-hand side of Eq. (6), we get
For the successful extrapolation process of the ETDLPA, the coupling coefficient should be available for all layer thickness values involved during the extrapolation procedure. Therefore, the grating position of the ETDLPA should be carefully treated because the TDLPA with different layer thickness produces the coupling coefficient at different grating positions. We introduce the notation q̂(z,mh,sh) to denote the coupling coefficient computed by the TDLPA with the layer thickness mh and the starting position sh. The reconstructed grating positions of q̂(z,mh,sh) are (s + lm)h,l = 0,1, ···. For convenience, we denote the P-stage ETDLPA with the layer thickness h and the TDLPA with the layer thickness h by ETDLPA(P, h) and TDLPA(h), respectively. As an example to explain the grating position problem, let us consider ETDLPA(2, h). ETDLPA(2, h) combines TDLPA(h), TDLPA(2h) and TDLPA(3h) according to Eq. (6). The solution of Eq. (8) for P = 2 is given by
The reconstructed grating positions of q̂(z,h,0), q̂(z,2h,0), and q̂(z,3h,0) are given byEqs. (11)–(13), using q̂(z,h,0), q̂(z,2h,0), and q̂(z,3h,0) as
To get the extrapolated coupling coefficient at the grating positions (6l + 1)h,l = 0,1,2,···, we need to move the starting position from 0 to h. Similarly to the previous discussion about the starting position 0 to compute the extrapolated coupling coefficients at the grating positions 6lh,l = 0,1,···, we can get the extrapolated value at the grating positions (6l + 1)h,l = 0,1,2,··· using q̂(z,h,h), q̂(z,2h,h), and q̂(z,3h,h) as
Algorithm 1. ETDLPA( P , h )
Set P and h ;
Compute L = LCM(1,2,··· ,P + 1) ;
For s = 0,1,··· ,L – 1
For m = 1,2,··· ,P + 1
Compute q(z)|z =( Ll + s ) h, l = 0, 1, ··· using Eq. (19)
The TDLPA with N-layers involves N 2 + N multiplications. The P-stage ETDLPA with N-layers involvesEq. (20), the first term and the second term are the number of the multiplications to compute q̂(z,mh,sh) and Eq. (19), respectively. If we compare the dominant N 2-term of Eq. (20) with that of the TDLPA, the number of the involved multiplications of the ETDLPA is increased by the factor of . This factor is very small for the practical P. The increment of the computational cost of the ETDLPA compared to the TDLPA is negligible because the involved layer thickness values during the extrapolation procedure are larger than the given layer thickness h. For example to illustrate the computationally efficient feature of the proposed ETDLPA, let us consider h = 10–3. The layer thickness of the TDLPA should be h 2 to achieve the accuracy O(h 2), and then the number of layers of the TDLPA is increased to 103 N. Therefore, we need 106 N 2 + 103 N multiplications to get accuracy O(h 2) using the TDLPA. However, the accuracy O(h 2) can be obtained through ETDLPA(1, h), and the involved multiplications of ETDLPA(1, h) are just which is nearly same as that of TDLPA(h). The larger extrapolation stage better reveals the excellence of the ETDLPA. The involved multiplications of ETDLPA(2, h) are and the accuracy is improved to O(h 3). But TDLPA(h 3) requires 1012 N 2 + 106 N multiplications to get the accuracy O(h 3).
3. Numerical examples
To demonstrate the performance of the proposed ETDLPA, we compare the ETDLPA with the TDLPA using two examples. First, for a flat-top bandpass filter, we compare the synthesis results of the ETDLPA and TDLPA to verify the accuracy and the computational efficiency of the proposed ETDLPA. Then, we detail an another example to show more concretely how the ETDLPA is able to improve the accuracy and the numerical aspects of the conventinoal TDLPA.
3.1. Nonuniform FBG synthesis example: flat-top bandpass filter
For nonuniform FBG synthesis example, we consider a flat-top bandpass filter described by the impulse response g(t) and the spectral response G(δ),
The synthesized coupling coefficient and the reflectivity are shown in Fig. 1. In this figure, we compare ETDLPA(2, h) with TDLPA(h). The result of TDLPA(h/1000) is also plotted to show the enhanced accuracy of the ETDLPA in Fig. 1. For a fair comparison, the coupling coefficient of TDLPA(h/1000) is under-sampled 1,000-times and the reflectivity is plotted by using this under-sampled coupling coefficient. The part of Figs. 1(a) and 1(c) marked by ellipsoid enlarged to demonstrate clearly the improved accuracy of the ETDLPA at Figs. 1(b) and 1(d), respectively. The target reflectivity shown in Figs. 1(a) and 1(b) is plotted from the time-shifted and apodized impulse response. ETDLPA(2, h) shows more accurate synthesis result than TDLPA(h) as seen from Fig. 1(b). Also, the reflectivity of ETDLPA(2, h) is even closer to the target reflectivity than TDLPA(h/1000), but the involved multiplications of ETDLPA(2, h) are only 2 × 10–6 of that of TDLPA(h/1000). Therefore, the ETDLPA can synthesize nonuniform FBG in a computationally efficient manner with high accuracy. Note that the computational complexity of ETDLPA(2, h) is just two times of that of TDLPA(h). In Figs. 1(c) and 1(d), we compare the accuracy of the synthesized coupling coefficient by using the result of the TDLPA with smaller layer thickness, although we can not know the exact coupling coefficient for the target reflectivity. It can be seen in Fig. 1(d) that the coupling coefficient reconstructed by ETDLPA(2, h) is nearly same as that of TDLPA(h/1,000). But it shows a gap up to 0.2 in the coupling coefficient between TDLPA(h/1,000) and TDLPA(h). The shape of the coupling coefficient is similar to the impulse response as mentioned in .
Therefore, the proposed ETDLPA provides highly accurate synthesis results with dramatically reduced computational cost.
3.2. Third order Butterworth filter
We verify the accuracy and the computational efficiency of the ETDLPA more quantitatively using the identification problem of the coupling coefficient. The coupling coefficient of the Butterworth filter can be computed exactly by the Gelfand-Levitan-Marchenko method , and therefore, the Butterworth filter is a good example to illustrate the accuracy of the ETDLPA. The measure of the accuracy is the normalized error defined by
We consider the third order Butterworth filter with the transfer function,Table 1. TDLPA(h 2) and ETDLPA(1, h) have the same accuracy O(h 2) and are grouped in Table 1. Similarly, TDLPA(h 3) and ETDLPA(2, h) are grouped in Table 1 because those have the same accuracy O(h 3). The validity of Eq. (2) is verified indirectly in this table by observing the trend of the normalized error of TDLPA(h), TDLPA(h 2), and TDLPA(h 3). The normalized error of the TDLPA is reduced exactly by the decreasing factor of the layer thickness, which corresponds to Eq. (2). For example, the layer thickness and the normalized error of TDLPA(h 2) are 1/10 of those of TDLPA(h), respectively. Also, it is observed in Table 1 that the accuracy order of the ETDLPA is properly increased according to Eq. (6). The normalized error of ETDLPA(1, h) is reduced to 1.1 × 10–2 from 7.1 × 10–2 which is the normalized error of DLPA(h), and the normalized error of ETDLPA(2, h) is reduced to 3.2 × 10–3 from 1.1 × 10–2 which is the normalized error of ETDLPA(1, h). Although we can conclude in this table that ETDLPA(P, h) and TDLPA(hP +1) have the same accuracy, it can be seen that the normalized error of ETDLPA(1, h) and ETDLPA(2, h) are slightly larger than that of TDLPA(h 2) and TDLPA(h 3), respectively. This differences are resulted from the different representation of the remained error terms. For example, the dominant error term of TDLPA(h 2) and ETDLPA(1, h) are a 1 h 2 and −2a 2 h 2, respectively. The computational efficiency of the ETDLPA is cleared by comparing the involved multiplications of the ETDLPA and the TDLPA in Table 1. The multiplications of ETDLPA(1, h) are less than those of TDLPA(h 2) by the factor of 1/64, in spite of the fact that ETDLPA(1, h) and TDLPA(h 2) have the same accuracy. This computational efficiency of the ETDLPA is more increased for higher P as seen by comparing the multiplications of ETDLPA(2, h) and TDLPA(h 3) in Table 1. The coupling coefficient computed by ETDLPA(2, h) and TDLPA(h) are plotted in Fig. 2 for the illustrative purpose. We can see the large gap between the TDLPA(h) and the exact value, but the result of ETDLPA(2, h) can not be discriminated from the exact value with the naked eye in this figure. Also, we can see, again, the similar shape between the computed coupling coefficients and the impulse response as mentioned in .
In Table 2, the layer thickness is decreased from 0.1 to 0.01 and the simulation conditions except for the layer thickness are same as Table 1. In this table, the trend of the normalized error and the multiplications are similar to Table 1, and the accuracy and the computational efficiency of the ETDLPA can be shown more clearly. The computational advantage of the ETDLPA is more increased for the smaller layer thickness. The running time of ETDLPA(2, h) using the latest PC is less than one second, but the running time of TDLPA(h 3) is greater than 168 hours.
We have proposed the ETDLPA that synthesizes nonuniform FBG more accurately than the conventional TDLPA, under comparable computational complexity. The accuracy and the computational efficiency of the ETDLPA have been studied theoritically and verified using numerical examples, concluding that the proposed ETDLPA resolves the tradeoff between the computational complexity and the accuracy of the conventional TDLPA. The ETDLPA is applicable to various optical devices that operate under the one-dimensional inverse scattering principle.
Appendix A Error formula of the TDLPA
We discuss briefly the discretization process of Eq. (1) and derive Eq. (2). It is well known that the continuous problem represented by Eq. (1) can be discretized as the equispaced discrete model, and the TDLPA computes the coupling coefficient using this discrete model. However, for the self-contained derivation of Eq. (2), we start from the fundamental matrix model for a uniform grating to get the discrete model, as done in [2, 3]. The solution of Eq. (1) for the uniform grating of length h is given by [2, 3],Eq. (1), the approximation of F is given by 
From the Taylor series expansion, we getEq. (27) is O(h 2). The matrix R can be rewritten into the familiar hyperbolic rotation matrix , Eq. (34) as follows, Eq. (35) does not degrade the accuracy of Eq. (27) because the accuracy order of Eq. (35) is higher than that of Eq. (27).
The final form of the discretization of Eq. (1) is given byEq. (36) with the layer thickness h and the exact wave variables U and V, and the accuracy of this discretization is O(h 2).Eq. (37), we can compute the discretized wave variables ûh and v̂h with the accuracy O(h 2), and thus, we can write as follows, Eq. (36) and Eq. (37). From the causality of the wave propagation, the reflection coefficient can be computed from Eq. (37) as follows [17, 18] Eq. (38) and Eq. (39) into Eq. (40), it is easy to see that we can rewrite Eq. (40) as follows,
Furthermore, we investigate the relation between and q(z) in Eq. (41). If the arrival time of the right-propagating wave u at position z is t, the left-propagating wave v cannot exist at position z + h for a time lower than t due to the causality of the wave propagation. Therefore, we getEq. (1), we know that Eq. (42) and Eq. (45) into Eq. (43), Eq. (43) can be rewritten as Eq. (46), it is interesting to note that means the reflection coefficient because it is the ratio of the reflected wave v(z,t) to the incident wave u(z,t), and thus hq(z) approximates the reflection coefficient with accuracy O(h 2). Inserting Eq. (47) into Eq. (41), we get Eq. (48) by h, we can obtain Eq. (2). Note that the explicit representation of αl,βl,γl and ξl are not detailed because these values do not need for a numerical extrapolation. It is clear to see that Eq. (2) also holds for the piecewise uniform model by the concatenation of the transfer matrices. The TDLPA is consisted of Eq. (37) and Eq. (40). The accuracy of F̂ is O(h 3) if the discrete reflector is placed at the center of the layer as described in . However, although the accuracy order of F̂ is increased from O(h 2) to O(h 3), it cannot increase the accuracy order of Eq. (2) because the accuracy of Eq. (46) is O(h 2). Moreover, Eq. (37) is appropriate for the general identification problem of the local reflectivity in spite of the lower accuracy order by one.
Note that the TDLPA is appropriate for the parallel computation based on the generator matrix representation . The TDLPA needs only O(N) computing time units with N-processors [17–19], where N is the total number of layers. See the Chapter 10 of  for the detailed procedure of the parallel implementation of the TDLPA.
This work was supported by the Ministry of Land, Transport and Maritime Affairs of Korea under grant 20043001 and D10813810H380000110.
References and links
1. G. H. Song and S. Y. Shin, “Design of corrugated waveguide filters by the Gel’fand-Levitan-Marchenko inverse-scattering method,” J. Opt. Soc. Am. A 2, 1905–1915 (1985). [CrossRef]
2. R. Feced, M. N. Zervas, and M. A. Muriel, “An efficient inverse scattering algorithm for the design of nonuniform fiber bragg gratings,” IEEE J. Quantum Electron . 35, 1105–1115 (1999). [CrossRef]
3. J. Skaar, L. Wang, and T. Erdogan, “On the synthesis of fiber bragg gratings by layer peeling,” IEEE J. Quantum Electron . 37, 165–173 (2001). [CrossRef]
4. J. Bae, J. Chun, and T. Kailath, “The Schur algorithm applied to the design of optical multi-mirror structures,” Numer. Linear Algebra Appl. 12, 283–292 (2005). [CrossRef]
5. J. Skaar and K. M. Risvik, “A genetic algorithm for the inverse problem in synthesis of fiber gratings,” J. Light-wave Technol. 16, 1928–1932 (1998). [CrossRef]
6. J. Bae and J. Chun, “Numerical optimization approach for designing bandpass filters using fiber Bragg gratings,” Opt. Eng. 42, 23–29 (2003). [CrossRef]
8. A. Rosenthal and M. Horowitz, “Inverse scattering algorithm for reconstructing strongly reflecting fiber bragg gratings,” IEEE J. Quantum Electron . 39, 1018–1026 (2003). [CrossRef]
9. J. Skaar, “Inverse scattering for one-dimensional periodic optical structures and application to design and characterization,” presented in Bragg Gratings, Photosensitivity, and Poling in Glass Waveguides, Sydney, Australia, 4–8 July 2005.
10. H. Li, T. Kumagai, and K. Ogusu, “Advanced design of a multichannel fiber Bragg grating based on a layer-peeling method,” J. Opt. Soc. Am. B 21, 1929–1938 (2004). [CrossRef]
11. M. Li, J. Hayashi, and H. Li, “Advanced design of a complex fiber Bragg grating for a multichannel asymmetrical triangular filter,” J. Opt. Soc. Am. B 26228–234 (2009). [CrossRef]
12. A. Buryak, J. Bland-Hawthorn, and V. Steblina, “Comparison of Inverse Scattering Algorithms for Designing Ultrabroadband Fibre Bragg Gratings,” Opt. Express 17, 1995–2004 (2009). [CrossRef] [PubMed]
13. Y. Gong, A. Lin, X. Hu, L. Wang, and X. Liu, “Optimal method of designing triangular-spectrum fiber Bragg gratings with low index modulation and chirp-free structure,” J. Opt. Soc. Am. B 26, 1042–1048 (2009). [CrossRef]
15. M. H. Asghari and J. Azana, “On the Design of Efficient and Accurate Arbitrary-Order Temporal Optical Integrators Using Fiber Bragg Gratings,” J. Lightwave Technol. 27, 3888–3895 (2009). [CrossRef]
16. A. M. Bruckstein, I. Koltracht, and T. Kailath, “Inverse scattering with noisy data,” SIAM J. Sci. Stat. Comput. 7, 1331–1349 (1986). [CrossRef]
17. T. Kailath, “Signal processing applications of some moment problems,” in Proceedings of Symposia in Applied Mathematics (American Mathmatical Society, 1976), pp. 71–109.
18. A. M. Bruckstein and T. Kailath, “Inverse scattering for discrete transmission-line models,” SIAM Review 29, 359–389 (1987). [CrossRef]
19. S. Haykin, Mordern Filters (Macmillan Publishing Company, 1990).
20. T. Erdogan, “Fiber grating spectra,” J. Lightwave Technol. 15, 1277–1294 (1997). [CrossRef]
21. A. Sidi, Practical Extrapolation Methods (Cambridge University Press, 2003). [CrossRef]