Numerically extrapolated discrete layer-peeling algorithm for synthesis of nonuniform fiber Bragg gratings

Youngchol Choi; Joohwan Chun; Jinho Bae

doi:10.1364/OE.19.008254

1. Introduction

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 [1] 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 [2] and [3], 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) [3], although the TDLPA shows degraded performance for a strong uniform grating [9] 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 [16]. 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 [3], 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 [20]. 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 [21] 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.

2. ETDLPA

The synthesis of nonuniform FBG is the one-dimensinoal inverse scattering problem represented by the coupled mode equations in the space-frequency domain [2]

\frac{\partial}{\partial z} = [\begin{array}{l} U (z, δ) \\ V (z, δ) \end{array}] = [\begin{array}{l} i δ \\ q (z) \end{array} \begin{array}{r} q^{*} (z) \\ - i δ \end{array}] [\begin{array}{l} U (z, δ) \\ V (z, δ) \end{array}],

where U(z,δ) is the forward propagating wave, V(z,δ) is the backward propagating wave, q(z) is the coupling coefficient at the grating position z, and δ is the detuning factor. The goal of the synthesis of nonuniform FBG is to find q(z) given the desired spectrum described by the scattering data U(0,δ) and V(0,δ).

This synthesis problem can be solved by the TDLPA based on the discretization of Eq. (1) [3]. 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,

\hat{q} (z, h) = q (z) + \sum_{l = 1}^{\infty} a_{l} h^{l},

where q̂(z,h) is the coupling coefficient computed by the TDLPA with the layer thickness h and a _l is independent with h. A detailed derivation of 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,

\hat{q} (z, 2 h) = q (z) + \sum_{l = 1}^{\infty} a_{l} 2^{l} h^{l} .

Multiplying Eq. (2) by 2 and subtracting from Eq. (3), we get

2 \hat{q} (z, h) - \hat{q} (z, 2 h) = q (z) - 2 a_{2} h^{2} - 6 a_{3} h^{3} + \dots .

Equation (4) is only the first stage of the extrapolation procedure of the ETDLPA, in which the first degree term a ₁ h of Eq. (2) is removed in Eq. (4), and thus, the accuracy order is increased from O(h) to O(h ²). 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 ²) is obtained with the layer thickness of h ². If the number of layers with the layer thickness h is N, the number of layers with the layer thickness h ² 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 ²) to achieve the accuracy O(h ²) 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 ⁺¹). From Eq. (2), the result of the TDLPA with the layer thickness mh can be written by

\hat{q} (z, m h) = q (z) + \sum_{l = 1}^{\infty} a_{l} m^{l} h^{l}, 1 \leq m \leq P + 1.

Because the aim of the P-stage extrapolation procedure is to eliminate a ₁ h, a ₂ h ² , ··· ,a_Ph^P terms in Eq. (2) by combining q̂(z,h),q̂(z,2h),··· ,q̂(z,(P + 1)h), we get

\sum_{m = 1}^{P + 1} ω_{m} \hat{q} (z, m h) = q (z) + \sum_{l = P + 1}^{\infty} a_{k}^{'} h^{k},

where ω_m is the weight to be determined and

a_{k}^{'}

is independent with h. Inserting Eq. (5) into the left-hand side of Eq. (6), we get

\begin{array}{l} \sum_{m = 1}^{P + 1} ω_{m} \hat{q} (z, m h) & = & \sum_{m = 1}^{P + 1} ω_{m} (q (z) + \sum_{l = 1}^{\infty} a_{l} m^{l} h^{l}) \\ = & q (z) \sum_{m = 1}^{P + 1} ω_{m} + \sum_{l = 1}^{P} h^{l} a_{l} \sum_{m = 1}^{P + 1} ω_{m} m^{l} + \sum_{l = P + 1}^{\infty} h^{l} a_{l} \sum_{m = 1}^{P + 1} ω_{m} m^{l} . \end{array}

Comparing Eq. (6) and Eq. (7), ω_m can be obtained by solving the following system of linear equations,

A_{P} ω_{P} = e_{1},

where

A_{P} = [\begin{array}{l} 1 \\ 1 \\ 1 \\ ⋮ \\ 1 \end{array} \begin{matrix} \begin{array}{l} 1 \\ 2 \end{array} & \begin{array}{l} 1 \\ 3 \end{array} & \begin{array}{l} \dots \\ \dots \end{array} & \begin{array}{l} 1 \\ (P + 1) \end{array} \\ 2^{2} & 3^{2} & \dots & {(P + 1)}^{2} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 2^{P} & 3^{P} & \dots & {(P + 1)}^{P} \end{matrix}], ω_{P} = [\begin{matrix} \begin{array}{l} ω_{1} \\ ω_{2} \end{array} \\ ω_{3} \\ ⋮ \\ ω_{P + 1} \end{matrix}], e_{1} = [\begin{array}{l} 1 \\ 0 \\ 0 \\ ⋮ \\ 0 \end{array}] .

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

ω_{2}^{T} = [ω_{1} ω_{2} ω_{3}] = [3 - 3 1]

The reconstructed grating positions of q̂(z,h,0), q̂(z,2h,0), and q̂(z,3h,0) are given by

\hat{q} (z, h, 0) : 0, h, 2 h, 3 h, 4 h, 5 h, 6 h, 7 h, 8 h, 9 h, 10 h, 11 h, 12 h, 13 h, \dots

\hat{q} (z, 2 h, 0) : 0, 2 h, 4 h, 6 h, 8 h, 10 h, 12 h, 14 h, \dots

\hat{q} (z, 3 h, 0) : 0, 3 h, 6 h, 9 h, 12 h, 15 h, \dots .

The functions q̂(z,2h,0) and q̂(z,3h,0) do not have values at the grating position of h and we cannot get the extrapolated value at this grating position using q̂(z,h,0), q̂(z,2h,0), and q̂(z,3h,0). For grating positions except for the multiples of 6h, we can not perform the extrapolation procedure because the coupling coefficient at these grating positions is not available for at least one layer thickness. Therefore, we can perform the extrapolation process only at the grating positions 6lh,l = 0,1, ···, denoted by the the bold face in Eqs. (11)–(13), using q̂(z,h,0), q̂(z,2h,0), and q̂(z,3h,0) as

q (z) |_{z = (6 l) h} = \sum_{m = 1}^{3} ω_{m} \hat{q} (z, m h, 0) |_{z = (6 l) h}, l = 0, 1, \dots .

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

q (z) |_{z = (6 l + 1) h} = \sum_{m = 1}^{3} ω_{m} \hat{q} (z, m h, h) |_{z = (6 l + 1) h}, l = 0, 1, \dots .

Note that the reconstructed grating positions of q̂(z,h,h), q̂(z,2h,h), and q̂(z,3h,h) are given by

\hat{q} (z, h, h) : h, 2 h, 3 h, 4 h, 5 h, 6 h, 7 h, 8 h, 8 h, 9 h, 10 h, 11 h, 12 h, 13 h, 14 h, \dots

\hat{q} (z, 2 h, h) : h, 3 h, 5 h, 7 h, 9 h, 11 h, 13 h, 15 h, \dots

\hat{q} (z, 3 h, h) : h, 4 h, 7 h, 10 h, 13 h, 16 h, \dots

In general, the grating positions of q̂(z,h,sh),q̂(z,2h,sh),··· ,q̂(z,(P + 1)h,sh) coincide only at z = (s + lL)h,l = 0,1,···, where L is the least common multiple of 1,2,··· ,P + 1, i.e., L = LCM(1,2,··· ,P + 1). Therefore, one execution of the extrapolation procedure in the ET-DLPA can provide the coupling coefficient at equally Lh-spaced grating positions. For the reconstruction of the coupling coefficient at all grating positions, the extrapolation procedure of the ETDLPA should be performed L-times, with the different starting positions spanning from 0 to (L – 1)h. The P-stage extrapolation using q̂(z,mh,sh),m = 1,2,··· ,P + 1 is given by

q (z) |_{z = (L l + s) h} = \sum_{m = 1}^{P + 1} ω_{m} \hat{q} (z, m h, s h) |_{z = (L l + s) h}, l = 0, 1, \dots,

with the accuracy O(h^P ⁺¹). Now ETDLPA(P, h) can be summarized as follows.

Algorithm 1. ETDLPA( P , h )

Set P and h ;

Compute L = LCM(1,2,··· ,P + 1) ;

Construct Eq. (8) and compute ω _N by solving Eq. (8);

For s = 0,1,··· ,L – 1

For m = 1,2,··· ,P + 1

Compute q̂(z,mh,sh);

End

Compute q(z)|_z ₌₍ _Ll ₊ _s ₎ _h, l = 0, 1, ··· using Eq. (19)

End

The TDLPA with N-layers involves N ² + N multiplications. The P-stage ETDLPA with N-layers involves

\sum_{m = 1}^{P + 1} m (\frac{N^{2}}{m^{2}} + \frac{N}{m}) + L (P + 1) N = N^{2} \sum_{m = 1}^{P + 1} \frac{1}{m} + (L + 1) (P + 1) N

multiplications. In the left-hand side of Eq. (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 ²-term of Eq. (20) with that of the TDLPA, the number of the involved multiplications of the ETDLPA is increased by the factor of

\sum_{m = 1}^{P + 1} \frac{1}{m}

. 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 ² to achieve the accuracy O(h ²), and then the number of layers of the TDLPA is increased to 10³ N. Therefore, we need 10⁶ N ² + 10³ N multiplications to get accuracy O(h ²) using the TDLPA. However, the accuracy O(h ²) can be obtained through ETDLPA(1, h), and the involved multiplications of ETDLPA(1, h) are just

\frac{3}{2} N^{2} + 2 N

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

\frac{11}{6} N^{2} + 3 N

and the accuracy is improved to O(h ³). But TDLPA(h ³) requires 10¹² N ² + 10⁶ N multiplications to get the accuracy O(h ³).

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(δ),

g (t) = \sqrt{A} \frac{sin (2 π B t)}{π t}, G (δ) = \sqrt{A} rect (\frac{δ}{B})

with 0 < A ≤ 1, where

rect (x) = {\begin{matrix} 1, & for | x | \leq 1 \\ 0, & for | x | > 1 \end{matrix},

B is the bandwidth. In this example, the grating length is 3.3 cm, A = 0.9, B = 1 nm, and h = 0.005 cm. The impulse response was time-shifted and apodized by a Hann function as done in [2, 3].

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 [1].

Fig. 1 Comparison of the ETDLPA and TDLPA for synthesis of a flat-top bandpass filter: (a) reflectivity, (b) enlarged plot of (a) marked by ellipsoid, (c) coupling coefficient, (d) enlarged plot of (a) marked by ellipsoid.

Download Full Size | PDF

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 [1], 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

\frac{\sum_{z} | q (z) - \hat{q} (z) |}{\sum_{z} | q (z) |},

where q(z) is the exact coupling coefficient and q̂(z) is the coupling coefficient computed by the TDLPA or the ETDLPA.

We consider the third order Butterworth filter with the transfer function,

G (s) = \frac{r_{0}}{\prod_{l = 1}^{3} (s - α_{l})}, α_{1} = exp [i \frac{π}{2} (1 + \frac{2 l - 1}{3})] .

The impulse response g(t) can be computed by the inverse Laplace transform of G(s). The length of the grating is 6 m and

r_{0} = \sqrt{0.3}

. We compare the normalized error and the involved multiplications of the ETDLPA and the TDLPA for synthesis of the third order Butterworth filter. The normalized error and the involved multiplications of the ETDLPA and TDLPA for h = 0.1 are summarized in Table 1. TDLPA(h ²) and ETDLPA(1, h) have the same accuracy O(h ²) and are grouped in Table 1. Similarly, TDLPA(h ³) and ETDLPA(2, h) are grouped in Table 1 because those have the same accuracy O(h ³). The validity of Eq. (2) is verified indirectly in this table by observing the trend of the normalized error of TDLPA(h), TDLPA(h ²), and TDLPA(h ³). 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 ²) 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(h^P ⁺¹) 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 ²) and TDLPA(h ³), respectively. This differences are resulted from the different representation of the remained error terms. For example, the dominant error term of TDLPA(h ²) and ETDLPA(1, h) are a ₁ h ² and −2a ₂ h ², 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 ²) by the factor of 1/64, in spite of the fact that ETDLPA(1, h) and TDLPA(h ²) 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 ³) 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 [1].

Fig. 2 Coupling coefficient computed by ETDLPA(2, h) and TDLPA(h) for the third order Butterworth filter when the layer thickness h is 0.1.

Download Full Size | PDF

Table 1. Comparison of the ETDLPA and TDLPA for Synthesis of the Third Order Butterworth Filter: Normalized Error and Involved Multiplications, h=0.1

View Table | View all tables in this article

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 ³) is greater than 168 hours.

Table 2. Comparison of the ETDLPA and TDLPA for Synthesis of the Third Order Butterworth Filter: Normalized Error and Involved Multiplications, h=0.01

View Table | View all tables in this article

4. Conclusions

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],

[\begin{array}{l} U (z + h, δ) \\ V (z + h, δ) \end{array}] = F [\begin{array}{l} U (z, δ) \\ V (z, δ) \end{array}],

where the transfer matrix F is

F = [\begin{matrix} cosh (γ h) + i \frac{δ}{γ} sinh (γ h) & \frac{q^{*}}{γ} sinh (γ h) \\ \frac{q}{γ} sinh (γ h) & cosh (γ h) - i \frac{δ}{γ} sinh (γ h) \end{matrix}]

and γ ² = |q|² – δ ². To discretize Eq. (1), the approximation of F is given by [2]

F \approx \hat{F} = T R,

where

T = [\begin{matrix} exp (i δ h) & 0 \\ 0 & exp (- i δ h) \end{matrix}], R = [\begin{matrix} cosh (| q | h) & \frac{q^{*}}{| q |} sinh (| q | h) \\ \frac{q}{| q |} sinh (| q | h) & cosh (| q | h) \end{matrix}] .

From the Taylor series expansion, we get

\begin{array}{l} F_{11} & = & cosh (γ h) + i \frac{δ}{γ} sinh (γ h) = (1 + \frac{γ^{2} h^{2}}{2} + O (h^{4})) + i \frac{δ}{γ} (γ h + O (h^{3})) \\ = & 1 + i δ h + \frac{γ^{2} h^{2}}{2} + O (h^{3}) \end{array}

\begin{array}{l} {\hat{F}}_{11} & = & exp (i δ h) cosh (| q | h) = (1 + i δ h - \frac{δ^{2} h^{2}}{2} + O (h^{3})) (1 + \frac{{| q |}^{2} h^{2}}{2} + O (h^{4})) \\ = & 1 + i δ h + \frac{γ^{2} h^{2}}{2} + O (h^{3}), \end{array}

F_{21} = \frac{q}{γ} sinh (γ h) = \frac{q}{γ} (γ h + O (h^{3})) = q h + O (h^{3}),

\begin{array}{l} {\hat{F}}_{21} & = & exp (- i δ h) \frac{q}{| q |} sinh (| q | h) = (1 + O (h)) (\frac{q}{| q |} (| q | h + O (h^{3}))) \\ = & q h + O (h^{2}), \end{array}

where A_ij is the (i, j) entry of the matrix A. Therefore, the accuracy of Eq. (27) is O(h ²). The matrix R can be rewritten into the familiar hyperbolic rotation matrix [3],

R = {(1 - {| κ (z) |}^{2})}^{- 1 / 2} [\begin{matrix} 1 & - κ^{*} (z) \\ - κ (z) & 1 \end{matrix}],

where

κ (z) = - tanh (h | q (z) | \frac{q^{*} (z)}{| q (z) |} .

This approximation means that the coupled mode equations can be discretized as the equispaced discrete reflectors. The matrix R explains the interaction of the waves using the reflection coefficient κ(z) at the grating position z, and the matrix T describes the propagation of the waves in the homogeneous interval (z,z + h). The reflection coefficient κ(z) can be approximated using the Taylor series expansion of the hyperbolic tangent function in Eq. (34) as follows,

κ (z) = - h q^{*} (z) + O (h^{3}) .

Note that the accuracy of 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 by

[\begin{array}{l} {\hat{U}}_{h} (z + h, δ) \\ {\hat{V}}_{h} (z + h, δ) \end{array}] = [\begin{matrix} exp (i δ h) & 0 \\ 0 & exp (- i δ h) \end{matrix}] {(1 - {| h q (z) |}^{2})}^{- 1 / 2} [\begin{matrix} 1 & h q (z) \\ h q^{*} (z) & 1 \end{matrix}] [\begin{matrix} U (z, δ) \\ V (z, δ) \end{matrix}],

where Û_h and V̂_h are used to discriminate between the approximated wave variables computed by Eq. (36) with the layer thickness h and the exact wave variables U and V, and the accuracy of this discretization is O(h ²).

The space-time representation of Eq. (36) is obtained by the inverse Fourier transforms of Eq. (36) as follows [3],

[\begin{matrix} {\hat{u}}_{h} (z + h, t + h) \\ {\hat{v}}_{h} (z + h, t - h) \end{matrix}] = {(1 - {| h q (z) |}^{2})}^{- 1 / 2} [\begin{matrix} 1 & h q (z) \\ h q^{*} (z) & 1 \end{matrix}] [\begin{matrix} u (z, t) \\ v (z, t) \end{matrix}],

where u, v, û_h and v̂_h are the inverse fourier transform of U, V, Û_h and V̂_h, respectively, and the accuracy of this discretization is also O(h ²). Using Eq. (37), we can compute the discretized wave variables û_h and v̂_h with the accuracy O(h ²), and thus, we can write as follows,

{\hat{u}}_{h} (z + h, t + h) = u (z + h, t + h) + \sum_{l = 2}^{\infty} α_{l} h^{l},

{\hat{v}}_{h} (z + h, t + h) = v (z + h, t + h) + \sum_{l = 2}^{\infty} β_{l} h^{l},

where α_l and β_l are independent with h. These power series representations follow from the Taylor series expansion used for the derivation of 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]

h \hat{q} (z + h, h) = \frac{{\hat{v}}_{h} (z + h, z + h)}{{\hat{u}}_{h} (z + h, z + h)},

where q̂(z + h,h) is the coupling coefficient reconstructed by the TDLPA with the layer thickness h at the position z + h. Inserting Eq. (38) and Eq. (39) into Eq. (40), it is easy to see that we can rewrite Eq. (40) as follows,

h \hat{q} (z + h, h) = \frac{v (z + h, t + h)}{u (z + h, t + h)} + \sum_{l = 2}^{\infty} γ_{l} h^{l},

where γ_l is independent with h.

Furthermore, we investigate the relation between $\frac{v (z + h, t + h)}{u (z + h, t + h)}$ 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 get

v (z + h, t - h) = 0.

The Taylor series for v(z + h,t – h) at (z,t) is

v (z + h, t - h) = v (z, t) + h v^{(1)} (z, t) + \sum_{l = 2}^{\infty} \frac{h^{l}}{l!} v^{(l)} (z, t),

where

v^{(n)} (z, t) = {(\frac{\partial}{\partial z} - \frac{\partial}{\partial t})}^{n} v (z, t) .

From the space-time representation of Eq. (1), we know that

v^{(1)} (z, t) = (\frac{\partial}{\partial z} - \frac{\partial}{\partial t}) v (z, t) = - q (z) u (z, t) .

Inserting Eq. (42) and Eq. (45) into Eq. (43), Eq. (43) can be rewritten as

h q (z) = \frac{v (z, t)}{u (z, t)} + \sum_{l = 2}^{\infty} \frac{h^{l}}{l!} \frac{v^{(1)} (z, t)}{u (z, t)} .

Therefore,

h q (z + h) = \frac{v (z + h, t + h)}{u (z + h, t + h)} + \sum_{l = 2}^{\infty} \frac{h^{l}}{l!} \frac{v^{(l)} (z + h, t + h)}{u (z + h, t + h)} .

In Eq. (46), it is interesting to note that

\frac{v (z, t)}{u (z, t)}

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 ²). Inserting Eq. (47) into Eq. (41), we get

h \hat{q} (z + h, h) = h q (z + h) + \sum_{l = 2}^{\infty} ξ_{l} h^{l},

where ξ_l is independent with h. Dividing the both sides of 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 ³) if the discrete reflector is placed at the center of the layer as described in [2]. However, although the accuracy order of F̂ is increased from O(h ²) to O(h ³), it cannot increase the accuracy order of Eq. (2) because the accuracy of Eq. (46) is O(h ²). 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 [17]. 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 [19] for the detailed procedure of the parallel implementation of the TDLPA.

Acknowledgments

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]

7. K. Aksnes and J. Skaar, “Design of short fiber Bragg gratings by use of optimization,” Appl. Opt. 43, 2226–2230 (2004). [CrossRef] [PubMed]

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]

14. M. Li and J. Yao, “All-fiber temporal photonic fractional Hilbert transformer based on a directly designed fiber Bragg grating,” Opt. Lett. 35, 223–225 (2010). [CrossRef] [PubMed]

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]

Numerically extrapolated discrete layer-peeling algorithm for synthesis of nonuniform fiber Bragg gratings

Abstract

1. Introduction

2. ETDLPA

3. Numerical examples

3.1. Nonuniform FBG synthesis example: flat-top bandpass filter

3.2. Third order Butterworth filter

4. Conclusions

Appendix A Error formula of the TDLPA

Acknowledgments

References and links

Cited By

Figures (2)

Tables (2)

Equations (48)

Optics Express

Synthesis Method	Normalized Error	Multiplications	Accuracy

TDLPA(h)	7.1 × 10^–2	3.6 × 10³	O(h)

TDLPA(h ²)	7.1 × 10^–3	3.6 × 10⁵	O(h ²)
ETDLPA(1, h)	1.1 × 10^–2	5.6 × 10³	O(h ²)

TDLPA(h ³)	7.1 × 10^–4	3.6 × 10⁷	O(h ³)
ETDLPA(2, h)	3.2 × 10^–3	7.7 × 10³	O(h ³)