Analytical modeling of a single channel nonlinear fiber optic system based on QPSK

Shiva Kumar; Sina Naderi Shahi; Dong Yang

doi:10.1364/OE.20.027740

1. Introduction

In a highly dispersive fiber, a signal pulse broadens significantly and thereby, it interacts nonlinearly with a large number of pulses in its neighborhood. This nonlinear interaction leads to ghost or echo pulses which is known as intra-channel four-wave mixing (IFWM) [1–3]. The propagation impairments due to IFWM in direct detection systems are analyzed in Refs. [4, 5]. Recently, the modeling of nonlinear distortion in coherent fiber optic systems has drawn significant interest [6–12]. In Ref. [6], propagation impairment due to four wave mixing (FWM) in coherent orthogonal frequency-division multiplexing (OFDM) systems is studied and symmetries are used in the conventional FWM model. In Ref. [7], an analytical expression for the probability density function (PDF) of the IFWM impairments is derived for coherent fiber optic systems based on phase-shift keying (PSK). In Ref. [8], an analytic expression for the nonlinear threshold is found by assuming that signal pulses in each symbol slot are delta functions. In Ref. [9], an analytical expression for the power spectral density (PSD) of the nonlinear interference in a WDM system is developed. To evaluate the PSD, it is necessary to carry out a triple numerical integration and the computational cost scales as M³ where M is the number of samples in frequency domain. In this paper we have developed an analytical expression for the PSD of intrachannel nonlinear distortion. With analytic simplifications, we found that the computational cost scales as ∼ N²M/8 where N is the total number of significant neighboring signal pulses. Typically, N is smaller than M leading to significant reduction in computational time. However, the direct comparison between these two approaches is not appropriate since Ref. [9] primarily focused on interchannel impairments whereas this paper deals with intrachannel impairments only. In Ref. [11], a general first order perturbation theory of a multichannel optical transmission system is developed and stationary phase approximation is done to evaluate the cross-phase modulation fluctuations. In Ref. [12], it is proposed to apply a large predispersion to an optical signal before the fiber transmission and stationary phase approximation is employed to approximate the solution of nonlinear Schrodinger equation (NLSE) in the limit of very strong initial predispersion. In this paper, the stationary phase approach is used to approximate the Fourier transform of the echo pulse in a single channel so that the computational cost of the PSD calculations can be reduced. The stationary phase approximation translates convolutions into simple multiplications leading to a simple closed form expression for the spectrum of the echo pulse. As a result, the spectrum of the echo pulse is found to be proportional to the product of the signal pulse spectra shifted by the amounts proportional to the temporal positions of the signal pulses. Finally, the PSD of the nonlinear distortion is added to the PSD of the amplified spontaneous emission (ASE) and the integration of the PSD over the receiver bandwidth leads to the total variance which is used to calculate the bit error ratio (BER). In this paper, we consider only the case of single polarization. It is straightforward to extend the approach for the case of two polarizations.

This paper is organized as follows. Analytical expressions for the PSD of the nonlinear distortions are derived in section 2. Stationary phase approximation to calculate the spectrum of the echo pulse due to IFWM is also discussed in section 2. In section 3, the analytical expressions for the variance of the nonlinear distortion and BER are validated using numerical simulations. Finally, in section 4, the contributions of this work are summarized.

2. Mathematical derivation of power spectral density

Let the fiber input be

u (t, 0) = \sqrt{P} \sum_{n = - N / 2}^{N / 2} a_{n} p (t - n T_{s}, 0),

where T_s is the symbol interval, p(t, 0) is the pulse shape function at z = 0, and

a_{n} = \frac{x_{n} + i y_{n}}{\sqrt{2}},

x_n and y_n are real random variables that take values ±1 with equal probability, respectively. The evolution of the optical field envelope is governed by the NLSE in the lossless form [7]

i \frac{\partial u}{\partial z} - \frac{β_{2}}{2} \frac{\partial^{2} u}{\partial t^{2}} + γ a^{2} (z) {| u |}^{2} u = 0,

where β₂ is the dispersion coefficient, γ is the nonlinear coefficient, a²(z) = exp(−αz) between amplifiers, and α is fiber loss coefficient. a²(z) is the loss/gain profile that includes the step amplification at the beginning of each span. Using the perturbation technique [13], the field envelope can be expanded as

u = u_{0} + γ u_{1} (t, z) + \dots .

Here u₀ represents the 0^th order solution which satisfies

i \frac{\partial u_{0}}{\partial z} - \frac{β_{2}}{2} \frac{\partial^{2} u_{0}}{\partial t^{2}} = 0 .

The first order correction u₁ is

i \frac{\partial u_{1}}{\partial z} - \frac{β_{2}}{2} \frac{\partial^{2} u_{1}}{\partial t^{2}} = - a^{2} (z) {| u_{0} |}^{2} u_{0} .

Taking the Fourier transform of Eq. (6), we find

i \frac{d {\tilde{u}}_{1}}{d z} - \frac{β_{2}}{2} {(2 π f)}^{2} {\tilde{u}}_{1} = - a^{2} (z) \tilde{b} (f, z),

where

\tilde{b} (f, z) = ℱ [{| u_{0} |}^{2} u_{0}],

{\tilde{u}}_{1} (f, z) = ℱ [u_{1} (t, z)],

ℱ denotes the Fourier transform. Assuming the perfect dispersion compensation at the receiver and with ũ₁(f, 0) = 0, Eq. (7) is solved to yield,

{\tilde{u}}_{1} (f, L_{tot}) = i \int_{0}^{L_{tot}} a^{2} (z) \tilde{b} (f, z) exp [- i β_{2} {(2 π f)}^{2} z / 2] d z,

where L_tot is the total transmission distance. The solution of Eq. (5) with the initial condition given by Eq. (1) is

u_{0} (t, z) = \sqrt{P} \sum_{n = - N / 2}^{N / 2} a_{n} p (t - n T_{s}, z),

where

p (t, z) = ℱ^{- 1} [\tilde{p} (f, z)],

\tilde{p} (f, z) = \tilde{p} (f, 0) exp [i β_{2} {(2 π f)}^{2} z / 2] .

From Eq. (8), we have

\begin{array}{l} \tilde{b} (f, z) & = & ℱ [u_{0} u_{0} u_{0}^{*}] \\ = & {\tilde{u}}_{0} (f, z) * {\tilde{u}}_{0} (f, z) * {\tilde{u}}_{0}^{*} (- f, z), \end{array}

where “*” denotes convolution. Using Eq. (11) in Eq. (14), we find

\begin{array}{l} \tilde{b} (f, z) & = & P^{3 / 2} \sum_{l = - N / 2}^{N / 2} \sum_{m = - N / 2}^{N / 2} \sum_{n = - N / 2}^{N / 2} a_{l} a_{m} a_{n}^{*} {[\tilde{p} (f, z) exp (i 2 π f l T_{s})] * \\ [\tilde{p} (f, z) exp (i 2 π f m T_{s})] * [{\tilde{p}}^{*} (- f, z) exp (i 2 π f n T_{s})]} \\ = & P^{3 / 2} \sum_{l} \sum_{m} \sum_{n} a_{l} a_{m} a_{n}^{*} {\tilde{X}}_{l, m, n} (f, z), \end{array}

where

{\tilde{X}}_{l, m, n} (f, z) = [\tilde{p} (f, z) exp (i 2 π f l T_{s})] * [\tilde{p} (f, z) exp (i 2 π f m T_{s})] * [{\tilde{p}}^{*} (- f, z) exp (i 2 π f n T_{s})] .

The summation in Eq. (15) is assumed to be from −N/2 to N/2. Substituting Eq. (15) into Eq. (10), we find

{\tilde{u}}_{1} (f, L_{tot}) = i P^{3 / 2} \sum_{l} \sum_{m} \sum_{n} a_{l} a_{m} a_{n}^{*} {\tilde{Y}}_{l, m, n} (f),

{\tilde{Y}}_{l, m, n} (f) = \int_{0}^{L_{tot}} a^{2} (z) exp (- i β_{2} {(2 π f)}^{2} z / 2) {\tilde{X}}_{l, m, n} (f, z) d z .

The distortion due to fiber nonlinearity is δũ_NL = γu₁. The PSD of the nonlinear distortion is defined as

ρ_{N L} (f) = lim_{N \to \infty} \frac{1}{(N + 1) T_{s}} E {{| δ {\tilde{u}}_{N L} (f) |}^{2}},

where E{} denotes the ensemble average. Using Eq. (17), Eq. (19) can be written as

ρ_{N L} (f) = lim_{N \to \infty} \frac{γ^{2} P^{3}}{(N + 1) T_{s}} \sum_{l} \sum_{m} \sum_{n} \sum_{l^{'}} \sum_{n^{'}} \sum_{m^{'}} E {a_{l} a_{l^{'}}^{*} a_{m} a_{m^{'}}^{*} a_{n}^{*} a_{n^{'}}} {\tilde{Y}}_{l, m, n} (f) {\tilde{Y}}_{l^{'}, m^{'}, n^{'}}^{*} (f) .

The PSD can be divided into two groups. They are (i) non-degenerate intra-channel four-wave mixing (ND-IFWM), and (ii) degenerate intra-channel four-wave mixing (D-IFWM). For constant intensity modulation such as QPSK, it can be shown that self-phase modulation (SPM) and intrachannel cross-phase modulation (IXPM) produce only a deterministic phase shift, which can be removed by the electrical equalizer. So, in this paper, we ignore SPM and IXPM.

2.1. ND-IFWM

Let us first consider the case l ≠ m ≠ n and l′ ≠ m′ ≠ n′. For QPSK signals, we have

E {a_{l} a_{l^{'}}^{*}} = K_{1} δ_{l l^{'}},

E {a_{l} a_{l^{'}}} = 0,

where

\begin{matrix} K_{1} = E {{| a_{l} |}^{2}} = 1, \\ E {a_{l} a_{l^{'}}^{*} a_{m} a_{m^{'}}^{*} a_{n}^{*} a_{n^{'}}} = [δ_{l l^{'}} δ_{m m^{'}} δ_{n n^{'}} + δ_{l m^{'}} δ_{l^{'} m} δ_{n n^{'}}] . \end{matrix}

δ is Kronecker delta function. Using Eqs. (21)–(23), Eq. (20) becomes,

ρ_{N D - I F W M} (f) = lim_{N \to \infty} \frac{2 γ^{2} P^{3}}{(N + 1) T_{s}} \sum_{l} \sum_{m} \sum_{n} {| {\tilde{Y}}_{l, m, n} (f) |}^{2} .

Equation (24) can be rewritten as

\begin{array}{l} ρ_{N D - I F W M} (f) & = & lim_{N \to \infty} \frac{2 γ^{2} P^{3}}{(N + 1) T_{s}} {\underset{l \neq m, l + m - n = - N / 2}{\sum_{l} \sum_{m} \sum_{n}} {| {\tilde{Y}}_{l, m, n} (f) |}^{2} \\ + \underset{l \neq m, l + m - n = - N / 2 + 1}{\sum_{l} \sum_{m} \sum_{n}} {| {\tilde{Y}}_{l, m, n} (f) |}^{2} + \dots + \underset{l \neq m, l + m - n = - N / 2}{\sum_{l} \sum_{m} \sum_{n}} {| {\tilde{Y}}_{l, m, n} (f) |}^{2}} . \end{array}

The signal pulses located at lT_s, mT_s, and nT_s generate a echo pulse at qT_s = (l + m − n)T_s [1–5]. Therefore, q^th term on the right-hand side (RHS) of Eq. (25) represents the nonlinear distortion on the q^th symbol interval. Due to symmetry, the ensemble average of the nonlinear distortion should be the same on each symbol interval. In the other words, each term on the RHS of Eq. (25) should be equal, which yields

\begin{array}{l} ρ_{N D - I F W M} (f) & = & \frac{2 γ^{2} P^{3}}{T_{s}} \underset{l \neq m, l + m - n = 0}{\sum_{l} \sum_{m} \sum_{n}} {| {\tilde{Y}}_{l, m, n} (f) |}^{2} \\ = & \frac{2 γ^{2} P^{3}}{T_{s}} \underset{l \neq m}{\sum_{l} \sum_{m}} {\tilde{Z}}_{l, m} (f), \end{array}

where

{\tilde{Z}}_{l, m} (f) = {| {\tilde{Y}}_{l, m, l + m} (f) |}^{2} .

2.2. D-IFWM

Next, let us consider the case, l = m ≠ n and l′ = m′ ≠ n′. In this case

E {a_{l}^{2} {(a_{l^{'}}^{*})}^{2} a_{n}^{*} a_{n^{'}}} = K_{1} K_{2} {δ_{l l^{'}} δ_{n n^{'}}},

where

K_{2} = E {{| a_{l} |}^{4}} = 1 .

Using Eqs. (28) and (29) in Eq. (20), PSD due to D-IFWM is

ρ_{D - I F W M} = lim_{N \to \infty} \frac{γ^{2} P^{3}}{(N + 1) T_{s}} \sum_{l} \sum_{n} {| {\tilde{Y}}_{l, l, n} |}^{2} .

Proceeding as before, we find

ρ_{D - I F W M} = \frac{γ^{2} P^{3}}{T_{s}} \sum_{l} {\tilde{Z}}_{l, l} (f) .

2.3. Correlation between D-IFWM and ND-IFWM

Consider the case l = m ≠ n and l′ = m′ ≠ n′. Now

\begin{array}{l} E {a_{l} a_{l^{'}}^{*} a_{m} a_{m^{'}}^{*} a_{n}^{*} a_{n^{'}}} & = & E {a_{l}^{2} a_{l^{'}}^{*} . a_{m^{'}}^{*} a_{n}^{*} a_{n^{'}}} \\ = & 0, \end{array}

since

E {a_{l}^{2} a_{l^{'}}^{*}} = {\begin{array}{l} 0, & if l \neq l^{'} \\ E {{| a_{l} |}^{2} a_{l}} = 0, & if l = l^{'} \end{array}

Therefore, from Eq. (20), we have ρ_NL(f) = 0. In other words, there is no correlation between D-IFWM and ND-IFWM.

2.4. Total PSD

Total PSD is given by

ρ_{N L} (f) = ρ_{N D - I F W M} (f) + ρ_{D - I F W M} (f),

where ρ_ND_−IFWM(f), and ρ_D_−IFWM(f) are given by Eqs. (26) and (31), respectively. For Gaussian pulses, X̃_l,m,n(f) can be calculated as (see Appendix A)

{\tilde{X}}_{l, m, l + m} = D exp (A f^{2} + B f + C),

A = - \frac{(ξ^{2} + δ^{2})}{3 ξ + i δ},

B = \frac{i 4 π f T_{s} (l + m) ξ}{3 ξ + i δ},

C = - \frac{2 π^{2} T_{s}^{2} [(l^{2} + m^{2}) ξ - l m (ξ + i δ)]}{(3 ξ + i δ) (ξ - i δ)},

D = \frac{k^{3} π}{\sqrt{(ξ - i δ) (3 ξ + i δ)}},

\begin{array}{l} δ = 2 π^{2} β_{2} z, \\ k = \sqrt{2 π} T_{0}, \\ ξ = 2 π^{2} T_{0}^{2} . \end{array}

and T₀ is 1/e pulse width.

2.5. Stationary phase approximation

For non-Gaussian pulse shapes, the convolution in Eq. (16) cannot be evaluated analytically. Due to the rapidly varying phase of p̃(f, z) in Eq. (13), stationary phase approximation can be employed to approximate X̃_l,m,n(f). Stationary phase method is a standard technique for evaluating the integrals of the form [14]

I = \int_{- \infty}^{\infty} G (x) e^{i y (x)} d x,

where y(x) is a fast-varying function of x over most of the range of integration and G(x) is a slowly varying function. At the rapidly varying regions of y(x), the contribution to the integral is approximately zero as the area under the high frequency sinusoids with its slowly varying envelope G(x) is close to zero. The only significant contributions to the integral occurs in the regions where dy/dx = 0, i.e. at the points where the phase is stationary. At the vicinity of stationary phase point, x₀, y(x) may be written as

y (x) = y (x_{0}) + \frac{1}{2} y^{″} (x_{0}) {(x - x_{0})}^{2} .

Using Eq. (42), Eq. (41) may be approximated as

I \approx G (x_{0}) e^{i y (x_{0})} \int_{- \infty}^{\infty} e^{i y^{″} (x_{0}) {(x - x_{0})}^{2} / 2} d x,

\approx G (x_{0}) e^{i y (x_{0})} \sqrt{\frac{2 π}{i y^{″} (x_{0})}}

Now returning to Eq. (16), it has double convolutions which can be analytically integrated using the stationary phase approximation when the dispersion is sufficiently large. Now Eq. (16) becomes (see Appendix B)

\begin{array}{l} {\tilde{X}}_{l, m, l + m} (f, z) = \frac{π}{| δ |} \tilde{p} (f - \frac{π l T_{s}}{δ}, 0) \tilde{p} (f - \frac{π m T_{s}}{δ}, 0) \\ \tilde{p} (- f + \frac{π (l + m) T_{s}}{δ}, 0) exp [i (δ f^{2} + \frac{2 π^{2} T_{s}^{2} l m}{δ})] . \end{array}

Note that the convolution in Eq. (16) is hard to evaluate numerically unless the pulse shape is Gaussian. But the stationary phase approximation translates the convolutions into simple multiplications as shown in Eq. (45), which can be easily computed. When the Nyquist pulse such as sinc pulse is used, Eq. (45) can be further simplified. A sinc pulse has a rectangular spectrum,

\tilde{p} (f) = {\begin{array}{l} 1, & | f | \leq B_{s} / 2 \\ 0, & otherwise \end{array}

where B_s = 1/T_s. Using Eq. (46), Eq. (45) can be approximated as

{\tilde{X}}_{l, m, l + m} (f, z) = \frac{π}{| δ |} {\tilde{p}}_{l, m}^{'} (f) exp [i (δ f^{2} + \frac{2 π^{2} T_{s}^{2} l m}{δ})],

where

{\tilde{p}}_{l, m}^{'} (f) = {\begin{array}{l} 1, & lt \leq f \leq rt \\ 0, & otherwise \end{array}

lt = \max (l, m, l + m) \frac{π T_{s}}{δ} - \frac{B_{s}}{2},

rt = \min (l, m, l + m) \frac{π T_{s}}{δ} + \frac{B_{s}}{2} .

Property 1:

X_l,m,l+m is invariant under the exchange of l and m.

X_{l, m, l + m} (f) = X_{m, l, l + m} (f) .

This property holds true in general (see Eq. (16)) even without the stationary phase approximation.

Property 2:

Since we have assumed that p(t) is real and symmetric, it follows that p̃(f) is symmetric and from Eq. (45), it is easy to see that there is a mirror symmetry,

X_{l, m, l + m} (f) = X_{- m, - l, - (l + m)} (- f) .

2.6. Variance

Using the Wiener-Khinchin theorem, the variance is obtained as

σ_{N L}^{2} = \int_{- \infty}^{+ \infty} ρ_{N L} (f) H_{rec} (f) d f,

where H_rec(f) is the receiver transfer function. Using Eq. (35), Eq. (53) may be written as

σ_{N L}^{2} = σ_{N D - I F W M}^{2} + σ_{D - I F W M}^{2},

σ_{r}^{2} = \int_{- \infty}^{+ \infty} ρ_{r} (f) H_{rec} (f) d f,

where r = ND − IFWM,D − IFWM, and H_rec is receiver transfer function. From the property 2, it follows that,

ρ_{r} (f) = ρ_{r} (- f), r = N D - I F W M, D - I F W M .

Now Eq. (57) may be written as

σ_{r}^{2} = \int_{- \infty}^{+ \infty} ρ_{r} (f) [H_{rec} (f) + H_{rec} (- f)] d f .

2.7. Computational cost

Figure 1 shows the classification of intrachannel impairments for N = 6 and computational cost associated with SPM, IXPM and IFWM. When Property 1 and Property 2 are not used, the computational cost per frequency calculations are as shown in Table 1.

Fig. 1 Classification of intrachannel nonlinear impairments (N = 6).

Download Full Size | PDF

Table 1. Computational cost of nonlinear impairments per frequency

View Table

When Property 1 is used, the computational cost per frequency for ND-IFWM is N(N −1)/2. If both Property 1 and Property 2 are used, the cost per frequency for ND-IFWM and D-IFWM are N²/4 and N/2, respectively, and total computational cost per frequency (ND-IFWM + D-IFWM) is N²/4+N/2. If there are M samples in the frequency domain, total computational cost is (N²/4+N/2)M. In addition, if |l +m| > N/2, |n| > N/2, and from Eq. (26), it follows that the signal pulse centered at nT_s with |n| > N/2 does not contribute significantly for the formation of the echo pulse at t = 0 and hence, such a triplet may be ignored. With this approximation, total computational cost scales as ∼ N²M/8 for large N. Validation of the stationary phase approximation is carried out in section 3.1.

3. Results and discussions

We carried out the numerical simulations of the fiber optic sytem using the split-step Fourier method in order to test the validity of our analytical model. The following parameters are used: fiber loss α = 0.2 dB/km, fiber nonlinear coefficient γ = 1.1 (W.km)⁻¹, symbol rate = 25 Gbaud, and modulation = QPSK. Gaussian pulses with full width at half maximum (FWHM) of T_FWHM = 20 psec are launched to the fiber to obtain Figs. 2, 3 and 6. Amplifiers spacing is 80 km. Multi-span fiber-optic system is simulated here. The dispersion is uncompensated in each span. At the receiver, the transmission fiber dispersion is fully compensated either optically or electrically. Laser phase noise, polarization effects, and the coherent receiver imperfections are ignored since the primary focus of this paper is to validate our analytical model for the fiber nonlinear impairments. For numerical simulation, a pseudo-random bit sequence (PRBS) of length 2¹⁵ − 1 is used for the calculation of the PSD as well as BER. A Gaussian filter with a full bandwidth of 100 GHz is used as the receiver filter. The significant number of neighbors, N = 20. Equation (45) provides a guideline for choosing the frequency resolution. The minimum frequency shift of the pulse spectrum is πT_s/δ. If the frequency resolution Δf is larger than πT_s/δ, errors occur in the computation of X̃_l,m,l+m(f, z). For a 20-span system and for |β₂| = 21 ps²/km, Δf = πT_s/δ = 0.189 GHz. With the 100 GHz bandwidth of the receiver filter, number of frequency samples, M = 527. Since N is much smaller than M, in this example the computational cost savings would be ∼ O ((527/20)²).

Fig. 2 Analytical and numerical variances vs. peak power for (a) 5-spans, and (b) 20-spans system (β₂ = −21 ps²/km).

Download Full Size | PDF

Fig. 3 Analytical and numerical variances vs. dispersion parameter for (a) 5-spans, and (b) 20-spans system (P_peak = 0 dBm).

Download Full Size | PDF

Figures 2(a) and 2(b) show the analytical and numerical variances as a function of the launch peak power for a 5-spans and 20-spans systems, respectively. To calculate the PSD numerically, we proceed as follows. The SPM and IXPM introduce a constant phase shift which is removed by multiplying the received signal by exp(iθ) where θ is found adaptively. We used adaptive least mean square (LMS) equalizer to compensate the phase shift. We assumed the following parameters for the LMS algorithm: Number of filter taps = 10, number of training sequence = 2¹⁰, number of samples/symbol = 2, and step size = 0.1. The numerical PSD due to the nonlinear distortion is computed by subtracting the optical field envelope at the transmitter from that at the receiver (after dispersion compensation and the phase shift removal) and then taking the Fourier transform of the difference. To account for the bit-pattern variations, numerical simulations were performed 20 times with different bit patterns and the average PSD is computed. Over a range of powers that is of practical interest for QPSK-based system (−6 dBm to 0 dBm), the discrepancy between the analytical model and the numerical model is less than 4% and 12% in 5-span (Fig. 2(a)) and 20-span (Fig. 2(b)) systems, respectively. Figures 3(a) and 3(b) show the analytical and numerical variances versus the accumulated dispersion for 5-spans and 20-spans systems, respectively. As can be seen, there is a good agreement between numerical simulations and analytical results for a 5-span system. For a 20-span system, there is a small discrepancy at large launch powers which is probably due to the truncation of the field up to the first order (see Eq. (4)). For the 5-span system, when the dispersion is small, the variance of the nonlinear distortion is quite small. However, it grows quickly and beyond 7 ps/nm.km, it decays slowly. For the 20-span system, the variance decreases slowly with the transmission fiber dispersion.

3.1. Stationary phase approximation with raised-cosine pulse

The raised-cosine spectrum is commonly used in communication because of its compact spectrum. In this case, p̃(f) is of the form [15],

\tilde{p} (f) = {\begin{array}{l} 1, & | f | \leq \frac{1 - a}{2 T_{s}} \\ \frac{1}{2} [1 - sin (\frac{π T_{s}}{a} (| f | - \frac{1}{2 T_{s}}))], & \frac{1 - a}{2 T_{s}} < | f | \leq \frac{1 + a}{2 T_{s}} \\ 0, & | f | > \frac{1 + a}{2 T_{s}} \end{array}

where a is the roll-off factor, and T_s is the symbol time interval. X̃_l,m,n can not be calculated exactly for this pulse shape. Using the stationary phase approximation (Eq. (45)), total PSD and the variance is calculated. The following parameters are used for raised-cosine pulse. Roll-off factor a = 1 and symbol time interval T_s = 40 psec are assumed. Figures 4(a) and 4(b) show the variance as a function of the launch peak power for a 5-span and 20-span systems, respectively. In the range of −6 dBm to 0 dBm, the discrepancy between the analytical model and the numerical model is less than 7% in Figs. 4(a) and 4(b). Figures 5(a) and 5(b) show the dependence of the variance on the fiber dispersion for a 5-span and 20-span systems, respectively. For a 5-span system, when the dispersion is very low, we see that the stationary phase approximation becomes inaccurate. This inaccuracy is due to the fact that the phase (∝ β₂ f²) does not vary rapidly at low dispersions. However, practical fiber-optic systems use fibers with moderate to large dispersions to suppress nonlinear effects and therefore, stationary phase approximation leads to reasonably accurate results for dispersion parameter range that is of practical interest.

Fig. 4 Validation of the stationary phase approximation. Analytical and numerical variances vs. launch peak power for (a) 5-spans, and (b) 20-spans system. Raised-cosine pulses are used. β₂ = −21 ps²/km.

Download Full Size | PDF

Fig. 5 Validation of the stationary phase approximation. Analytical and numerical variances vs. dispersion parameter for (a) 5-spans, and (b) 20-spans system. Raised-cosine pulses are used. P = 0 dBm.

Download Full Size | PDF

Fig. 6 Analytical and numerical BER vs. average launch power. Gaussian pulses are used. Number of spans=20 and β₂ = −9.1 ps²/km.

Download Full Size | PDF

3.2. BER calculation

In this section, analytical BER is compared with numerical BER. For the analytical BER, total noise variance is calculated first and then SNR is obtained. Ignoring the interplay between the amplifier noise and nonlinearity (so called Gordon-Mollenauer noise), total PSD is given by

ρ_{tot} (f) = ρ_{N L} (f) + ρ_{A S E} (f),

where ρ_ASE (f) is ASE noise PSD due to amplifiers,

ρ_{A S E} = \sum_{j = 1}^{N_{A}} (e^{α L_{a, j}} - 1) h \bar{f} n_{s p},

where N_A is number of amplifiers, L_a,j is amplifiers’ spacing, α is fiber loss, h is Planck constant, f̄ is the mean frequency of the channel, and n_sp is spontaneous noise factor. Using Eq. (59), total noise variance σ_tot is calculated and the probability of error P_e for QPSK is given by [15]

P_{e} = 2 Q (\sqrt{2 S N R}) - Q^{2} (\sqrt{2 S N R}),

where Q is Q-function [15] and SNR is the signal to noise ratio,

S N R = \frac{P_{a v}}{σ_{tot}^{2}},

P_{a v} = \frac{P}{1.88},

for Gaussian pulses with 50% duty cycle, and

σ_{tot}^{2} = \int ρ_{tot} (f) H_{rec} (f) d f .

The Q-factor is defined as

Q (lin) = \sqrt{2} {erfc}^{- 1} (2 P_{e}),

Q (d B Q_{20}) = 20 {log}_{10} Q (lin) .

Figure 6 shows the analytical and numerical BER versus the launch power. We assumed the following parameters: n_sp = 10 dB, and number of spans = 20. We chose a relatively large noise figure intentionally so as to reduce the computational time of Monte Carlo simulations at large launch power. We found that the maximum discrepancy between the analytical and numerical Q-factor is less than 0.6 dBQ₂₀ which is attributed to non-Gaussian distribution of the IFWM pdf [7]. In Eq. (61), it is assumed that the noise is Gaussian distributed. However, in Ref. [7], it is shown that the pdf of intrachannel impairments are asymmetric and non-Gaussian. Ref. [16] has modeled the pdf of the nonlinear interference as Gaussian distribution. When the perturbation includes the Gaussian-distributed ASE and the non-Gaussian-distributed IFWM, accurate evaluation of the BER may be carried out using the Gram-Charlier technique [17], which would be the subject of a future investigation.

4. Conclusion

We have developed analytical expressions for the PSD of the nonlinear distortions due to IFWM. Combining this PSD with that of the ASE, BER is estimated analytically which is found to be in good agreement with numerical simulations. For non-Gaussian pulse shapes, the spectrum of the the echo pulse can not be calculated analytically. When the phase is varying rapidly as in the case of a high dispersion transmission fiber, stationary phase approximation can be employed to calculate the spectrum of the echo pulse and hence, the PSD can be calculated analytically. The stationary phase approximation translates convolutions into simple multiplications leading to a simple analytical expression for the spectrum of the echo pulse. These analytical expressions significantly reduce the computational time to estimate the BER.

5. Appendix A: Gaussian pulse case

For a Gaussian pulse shape, we have

p (t, 0) = exp (- \frac{t^{2}}{2 T_{0}^{2}}),

and

\tilde{p} (f, 0) = k exp (- ξ f^{2}),

where

\begin{array}{l} k = \sqrt{2 π} T_{0}, \\ ξ = 2 π^{2} T_{0}^{2} . \end{array}

Equation (16) can be rewritten as

\begin{array}{l} {\tilde{X}}_{l, m, l + m} (f, z) & = & \int \int {\tilde{p}}_{1} (f_{1}) exp (i 2 π f_{1} l T_{s}) {\tilde{p}}_{2} (f_{2} - f_{1}) exp [i 2 π (f_{2} - f_{1}) m T_{s}] \\ {\tilde{p}}_{3} (f - f_{2}) exp [i 2 π (f - f_{2}) n T_{s}] d f_{1} d f_{2}, \end{array}

n = l + m,

{\tilde{p}}_{1} (f) = \tilde{p} (f, z) = \tilde{p} (f, 0) exp [i {(2 π f)}^{2} β_{2} z / 2],

{\tilde{p}}_{2} (f) = {\tilde{p}}_{1} (f),

{\tilde{p}}_{3} (f) = {\tilde{p}}_{1}^{*} (- f) = k exp [- ξ f^{2} - i (2 π f^{2}) β_{2} z / 2] .

Let

X_{l, m, n} (f, z) = \int M (f_{2}) {\tilde{p}}_{3} (f - f_{2}) exp [i 2 π (f - f_{2}) n T_{s}] d f_{2},

\begin{array}{l} M (f_{2}) & = & \int \tilde{p} (f_{1}, 0) exp [i 2 π f_{1} l T_{s} + i δ f_{1}^{2}] \tilde{p} (f_{2} - f_{1}, 0) \\ exp [i 2 π (f_{2} - f_{1}) m T_{s} + i δ {(f_{2} - f_{1})}^{2}] d f_{1}, \end{array}

where 2π²β₂z = δ. Using Eq. (68) in Eq. (76), we find

\begin{array}{l} M (f_{2}) & = & k^{2} \int exp [- ξ f_{1}^{2} + i δ f_{1}^{2} - ξ {(f_{2} - f_{1})}^{2} + i 2 π f_{1} l T_{s} \\ + i 2 π (f_{2} - f_{1}) m T_{S} + i δ {(f_{2} - f_{1})}^{2}] d f_{1}, \\ = & k^{2} exp [- ξ f_{2}^{2} + i δ f_{2}^{2} + i 2 π f_{2} m T_{s}] \\ \int exp {- 2 [ξ - i δ] f_{1}^{2} + i 2 π f_{1} (l - m) T_{s} + 2 ξ f_{1} f_{2} - i 2 δ f_{1} f_{2}} d f_{1} . \end{array}

Equation (77) may be rewritten as

M (f_{2}) = k^{2} exp [- (ξ - i δ) f_{2}^{2} + i 2 π f_{2} m T_{s}] I (f_{2}),

where

\begin{array}{l} I (f_{2}) & = & \int exp {- 2 (ξ - i δ) f_{1}^{2} + 2 f_{1} [i π (l - m) T_{s} + (ξ - i δ) f_{2}]} d f_{1}, \\ = & exp {\frac{{[(ξ - i δ) f_{2} + i π (l - m) T_{s}]}^{2}}{2 (ξ - i δ)}} \frac{\sqrt{π}}{\sqrt{2 (ξ - i δ)}} . \end{array}

Substituting Eq. (79) in Eq. (78), we find

M (f_{2}) = J exp (- β f_{2}^{2} + i μ f_{2}),

J = \frac{k^{2} \sqrt{π}}{\sqrt{2 (ξ - i δ)}} exp (\frac{- π^{2} {(l - m)}^{2} T_{s}^{2}}{2 (ξ - i δ)}),

β = \frac{ξ - i δ}{2},

μ = π T_{s} (l + m) .

Substituting Eq. (80) in Eq. (75) and noting that

{\tilde{p}}_{3} (f - f_{2}) = k exp [- (ξ + i δ) {(f - f_{2})}^{2}],

Eq. (75) becomes

\begin{array}{l} X_{l, m, l + m} & = & J exp [- (ξ + i δ) f^{2} + i 2 π f n T_{s}] \int exp [- (ξ \\ - i δ) f_{2}^{2} / 2 - (ξ + i δ) f_{2}^{2} + (ξ + i δ) 2 f f_{2}] d f_{2}, \\ = D exp (A f^{2} + B f + C), \end{array}

A = - \frac{(ξ^{2} + δ^{2}) f^{2}}{3 ξ + i δ},

B = \frac{i 4 π f T_{s} (l + m) ξ}{3 ξ + i δ},

C = - \frac{2 π^{2} T_{s}^{2} [(l^{2} + m^{2}) ξ - l m (ξ + i δ)]}{(3 ξ + i δ) (ξ - i δ)},

D = \frac{k^{3} π}{\sqrt{(ξ - i δ) (3 ξ + i δ)}} .

6. Appendix B: Stationary phase approximation

From Eq. (76), we have

\begin{array}{l} M (f_{2}) & = & \int \tilde{p} (f_{1}, 0) \tilde{p} (f_{2} - f_{1}, 0) exp [i 2 π f_{1} l T_{s} + i δ f_{1}^{2} + i δ {(f_{2} - f_{1})}^{2} + i 2 π (f_{2} - f_{1}) m T_{s}] d f_{1} \\ = & exp (i 2 π f_{2} m T_{s} + i δ f_{2}^{2}) I (f_{2}), \end{array}

I (f_{2}) = \int \tilde{p} (f_{1}, 0) \tilde{p} (f_{2} - f_{1}, 0) exp [i 2 δ f_{1}^{2} + i 2 π f_{1} (l - m) T_{s} - 2 i δ f_{1} f_{2}] d f_{1} .

Let

θ (f_{1}) = 2 δ f_{1}^{2} - 2 δ f_{1} f_{2} + f_{1} 2 π (l - m) T_{s},

so that

I (f_{2}) = \int \tilde{p} (f_{1}, 0) \tilde{p} (f_{2} - f_{1}, 0) exp (i θ (f_{1})) d f_{1} .

We assume that p(t) is real and symmetric with respect to origin so that p̃(f) is real. Since θ(f₁) is varying rapidly, the dominant contribution to the integral in Eq. (93) comes when

\frac{d θ}{d f_{1}} = 0,

or

f_{1, opt} = \frac{π (m - l) T_{s}}{2 δ} + \frac{f_{2}}{2} .

Substituting Eq. (95) in Eq. (92), we find

θ (f_{1, opt}) = - \frac{δ}{2} {[f_{2} + \frac{π (m - l) T_{s}}{δ}]}^{2} .

Under the stationary phase approximation, Eq. (93) becomes

I (f_{2}) ≅ \tilde{p} (f_{1, opt}, 0) \tilde{p} (f_{2} - f_{1, opt}, 0) exp [i θ (f_{1, opt})] l_{1},

\begin{array}{l} l_{1} & = & \sqrt{\frac{2 π}{| θ^{″} (f_{1, opt}) |}} exp [sgn (θ^{″} (f_{1, opt})) i π / 4], \\ = & \sqrt{\frac{π}{2 | δ |}} exp (i sgn (δ) π / 4), \end{array}

where ″ denotes differentiating twice. Using Eqs. (95) and (98) in Eq. (97), we find

\begin{array}{l} I (f_{2}) & = & l_{1} \tilde{p} (\frac{π (m - l) T_{s} + f_{2} δ}{2 δ}, 0) \tilde{p} (\frac{δ f_{2} - π (m - l) T_{s}}{2 δ}, 0) \\ exp {- \frac{i δ}{2} [f_{2} + \frac{π (m - l) T_{s}}{δ}]} . \end{array}

Substituting Eq. (99) in Eq. (90), we find

M (f_{2}) = s_{1} (f_{2}) exp [\frac{i δ f_{2}^{2}}{2} - \frac{i π^{2} {(m - l)}^{2} T_{s}^{2}}{2 δ} + i π f_{2} (m + l) T_{s}],

where

s_{1} (f_{2}) = \tilde{p} (\frac{π (m - l) T_{s} + f_{2} δ}{2 δ}, 0) \tilde{p} (\frac{δ f_{2} - π (m - l) T_{s}}{2 δ}, 0) .

Substituting Eq. (100) in Eq. (75) and simplifying, we obtain

{\tilde{X}}_{l, m, l + m} (f) = exp [- \frac{i π^{2} {(m - l)}^{2} T_{s}^{2}}{2 δ} - i δ f_{2}] \int s_{2} (f_{2}) exp (i θ_{2} (f_{2})) d f_{2},

where

s_{2} (f_{2}) = s_{1} (f_{2}) \tilde{p} (f - f_{2}, 0),

θ_{2} (f_{2}) = - \frac{δ f_{2}^{2}}{2} + 2 δ f f_{2} - π T_{s} f_{2} (l + m) .

Since θ₂ is varying rapidly, we use the stationary phase approximation again to evaluate the integral in Eq. (102). Differentiating Eq. (104) with respect to f₂ and setting the result to zero yields

f_{2, opt} = 2 f - \frac{π T_{s} (l + m)}{δ} .

Proceeding as before, Eq. (102) can be simplified as

\begin{array}{r} {\tilde{X}}_{l, m, l + m} (f) = \frac{π}{| δ |} \tilde{p} (f - \frac{π l T_{s}}{δ}, 0) \tilde{p} (f - \frac{π m T_{s}}{δ}, 0) \tilde{p} (- f \\ + \frac{π (l + m) T_{s}}{δ}, 0) exp [i (δ f^{2} + \frac{2 π^{2} T_{s}^{2} l m}{δ})] . \end{array}

References and links

1. I. Shake, H. Takara, K. Mori, S. Kawanishi, and Y. Yamabayashi, “Influence of inter-bit four-wave mixing in optical TDM transmission,” Electron. Lett. 34, 1600–1601 (1998). [CrossRef]

2. R. J. Essiambre, B. Mikkelsen, and G. Raybon, “Intra-channel crossphase modulation and four-wave mixing in high-speed TDM systems,” Electron. Lett. 35, 1576–1578 (1999). [CrossRef]

3. P. V. Mamyshev and N. A. Mamysheva, “Pulse-overlapped dispersion managed data transmission and intrachan-nel four-wave mixing,” Opt. Lett. 24, 1454–1456 (1999). [CrossRef]

4. A. Mecozzi, C. B. Clausen, and M. Shtaif, “Analysis of intrachannel nonlinear effects in highly dispersed optical pulse transmission,” IEEE Photon. Technol. Lett. 12, 292–294 (2000).

5. S. Kumar, “Intrachannel four-wave mixing in dispersion managed RZ systems,” IEEE Photon. Technol. Lett. 13, 800–802 (2001). [CrossRef]

6. M. Nazarathy, J. Khurgin, R. Weidenfeld, Y. Meiman, P. Cho, R. Noe, I. Shpantzer, and V. Karagodsky, “Phased-array cancellation of nonlinear FWM in coherent OFDM dispersive multi-span links,” Opt. Express 16, 15777–15810 (2008). [CrossRef] [PubMed]

7. D. Yang and S. Kumar, “Intra-channel four-wave mixing impairments in dispersion-managed coherent fiber-optic systems based on binary phase-shift keying,” J. Lightw. Technol. 27, 2916–2923 (2009). [CrossRef]

8. A. Bononi, P. Serena, N. Rossi, E. Grellier, and F. Vacondio, “Modeling monlinearity in coherent transmissions with dominant intrachannel-four-wave-mixing,” Opt. Express 20, 7777–7791 (2012). [CrossRef] [PubMed]

9. A. Carena, V. Curri, G. Bosco, P. Poggiolini, and F. Forghieri, “Modeling of the impact of nonlinear propagation effects in uncompensated optical coherent transmission links,” J. Lightw. Technol. 30, 1524–1539 (2012). [CrossRef]

10. P. Poggiolini, A. Carena, V. Curri, G. Bosco, and F. Forghieri, “Analytical modeling of non-Linear propagation in uncompensated optical transmission links,” IEEE Photon. Technol. Lett. 23, 742–744 (2011). [CrossRef]

11. A. Mecozzi and Rene Essiambre, “Nonlinear Shannon limit in pseudolinear coherent systems,” J. Lightw. Technol. 30, 2011–2024 (2012). [CrossRef]

12. S. Turitsyn, M. Sorokina, and S. Derevyanko, “Dispersion-dominated nonlinear fiber-optic channel,” Opt. Lett. 37, 2931–2933 (2012). [CrossRef] [PubMed]

13. S. Kumar and D. Yang, “Second-order theory for self-phase modulation and cross-phase modulation in optical fibers,” J. Lightw. Technol. 23, 2073–2080 (2005). [CrossRef]

14. J. D. Jackson, Classical Electrodynamics (Wiley, 1998).

15. U. Madhow, “Chapter 3: Demodulation,” in Fundamentals of Digital Communication (Cambridge University Press, 2008). [CrossRef]

16. A. Carena, G. Bosco, V. Curri, P. Poggiolini, M. T. Taiba, and F. Forghieri, “Statistical characterization of PM-QPSK signals after propagation in uncompensated fiber links,” in Proceedings ECOC, (2010), pp. 1–3.

17. M. Nazarathy, “Accurate evaluation of bit-error rates of optical communication systems using the Gram-Charlier series,” IEEE Trans. Commun. 54, 295–301 (2006). [CrossRef]

Analytical modeling of a single channel nonlinear fiber optic system based on QPSK

Abstract

1. Introduction

2. Mathematical derivation of power spectral density

2.1. ND-IFWM

2.2. D-IFWM

2.3. Correlation between D-IFWM and ND-IFWM

2.4. Total PSD

2.5. Stationary phase approximation

2.6. Variance

2.7. Computational cost

3. Results and discussions

3.1. Stationary phase approximation with raised-cosine pulse

3.2. BER calculation

4. Conclusion

5. Appendix A: Gaussian pulse case

6. Appendix B: Stationary phase approximation

References and links

Cited By

Figures (6)

Tables (1)

Equations (106)

Optics Express