A Novel Optical-Based Methodology for Improving Nonlinear Fourier Transform

Abstract

The increasing demand for bandwidth and long-haul transmission has led to new methods of signal processing and transmission in optical fiber communication systems. The nonlinear Fourier transform is one of the most recent methods proposed, and is able to represent an integrable nonlinear Schrödinger equation (NLSE) channel in terms of its continuous and discrete spectrum, to overcome the limitation of the bandwidth imposed by the Kerr effect on silica fibers. In this paper, we will propose and investigate the Boffetta-Osburne method for the direct nonlinear Fourier implementation, and the Gel’fand-Levitan-Marchenko equation for the inverse nonlinear Fourier, as only the continuous part of the nonlinear spectrum will be used to encode information. A novel methodology is proposed to improve their numerical implementation with respect to the NLSE, and we analyze in detail how the improved algorithm can be used in a real optical system, by investigating three different modulation schemes. We report increased performance transmission and consistency in the numerical results when the proposed methodology is applied. Our results show that b-modulation will increase the Q-factor by 2 dB with respect to the other two modulations. The improvement results with our proposed methodology suggest that b-modulation applied only to a continuous part of the nonlinear spectrum is a very effective method for maximizing both the transmission bandwidth and distance in optical fiber communication systems.

1. Introduction

The global data traffic is exponentially increasing due to large bandwidth requests coming from different services, such as high-definition interactive video, high-resolution editing video, and distributed systems [1]. The global internet bandwidth has increased by 28% in 2022, reaching approximately 1 Pbps [2]. Optical fiber is one of the most widely used transmission mediums when compared to other backbone networks. As opposed to various communication channels, optical fiber communication systems cannot be modeled as a linear channel under Gaussian perturbations, since the transmission in this medium is indeed nonlinear and mainly affected by three physical effects: (1) Kerr nonlinearity, (2) additive white Gaussian noise, and (3) chromatic dispersion. These effects are considered as undesirable, and digital signal processing is used to mitigate them at the receiver side [3,4].

Nowadays, a completely different approach is applied where nonlinearity is treated as a property of nonlinear channels. This can help to develop new transmission techniques. Nonlinear Fourier Transform (NFT) is one of the most recent methods studied for transmitting information in optical fiber communication systems in order to overpass the current state-of-the-art of multi-user transmission in wavelength division multiplexing (WDM) systems. From a mathematical point of view, the NFT is able to represent an integrable nonlinear dispersive partial differential equation (PDE) in terms of its continuous and discrete spectrum [5]. As the wave propagates over the channel under different nonlinear effects, the impact of the channel on its spectral components will be simply given by independent linear equations, as will be shown later. This technique is very effective in multi-user channels where the information can be encoded on the nonlinear part of the spectrum. This results in reduced inter-channel interference, improved signal to noise ratio, and better spectral efficiency [6,7,8,9]. On the other hand, the implementation of this technique requires physical changes in the current deployed infrastructure [10].

The NLSE is most widely solved using the numerical method, and hence, a larger number of numerical solutions have been proposed for it. The main reason behind the ability for the NLSE ability to be resolved numerically is that it does not have direct solutions known in analytical form, except for some specific cases. Taha et al. [5,11] made a comparative analysis of the numerical methods and the finite difference schemes used to solve the problem associated with NLSE. The analysis included the Ablowitz-Ladik scheme developed using notions related to the Inverse Scattering Transform (IST), nowadays called NFT in the telecommunications system community.

The logic behind NFT is quite simple. If the information signal would have to propagate through a highly nonlinear optical fiber governed by the NLSE equation, then its waveform at a desired distance could be associated with a linear operator that we need to find. This linear operator is an infinite dimensional matrix with a spectrum that consists of continuous and discrete eigenvalues. It turns out that these eigenvalues are invariant of the NLSE during the propagation through the fiber. The goal is to find a modulation scheme where we can modulate these eigenvalues, thus overcoming the nonlinear effect limitation. The first idea based on this simplified observation was published around 30 years ago, where eigenvalue communication was proposed to substitute for the on–off key solution for transmission in optical fiber communication systems [12]. The same idea was re-proposed and experimentally verified later in [13,14,15]. However, many researchers observed that there is no need to have an unchanged shape in the optical pulse in the NFT domain, as the channel effect will be a linear one. Moreover, this effect can be easily equalized at the receiver side, or pre-compensated at the transmitter side.

Several published papers exploit different degrees of freedom of the NFT and the Inverse NFT (INFT), using both discrete and continuous eigenvalues from the nonlinear spectrum. This results in a data transmission over the optical fiber similar to orthogonal frequency division multiplexing (OFDM) [16,17]. Hence, multiple users can access the optical channel.

Other researchers focused only on the discrete part of the spectrum and proved experimentally the reliability of the NFT and INFT [18,19,20,21,22,23,24]. The sole usage of the continuous part was studied in several papers [25,26,27,28,29]. Compared to using only the discrete part, this method resulted in increased performance in terms of spectral efficiency, transmission over long distances, and lower complexity at the receiver side. Dual polarization experimental transmission based on NFT was also reported in [30,31]. Full spectrum modulation using both continuous and discrete nonlinear parts of the spectrum was reported in [32], but this was limited only to back-to-back transmission.

Several numerical algorithms are proposed in the literature to solve the direct NFT and INFT [33]. To compute the direct NFT, we can generalize the algorithms into two main methods: (1) the Ablowitz-Ladik method [34,35] and (2) the Boffetta-Osborne Method [36]. To compute the INFT, we can generalize the algorithms into two main methods: the Darboux transformation method [37] and the Gelfand-Levitan-Marchenko (GLM) method [38]. Fast non linear Fourier transform (FNFT) is also analyzed in several papers, and is published online as a tool to reduce the implementation complexity of the NFT [39,40]. The purpose of FNFT is to develop an algorithm that computes a fast product between polynomials that are used to find the continuous spectrum. The NFT and INFT are relatively new methods, where the associated numerical implementations have not yet been fully addressed, especially its application in the field of optical fiber telecommunications.

Analytical solutions of the NLS have also been recently proposed, in [41,42,43,44,45]. Unfortunately, they cannot be directly applied for optical communication, and numerical solutions with lower complexity are still the only alternative to solve NLS.

In this paper, we propose a methodology to improve both the direct NFT and INFT numerical implementations in order to obtain only the continuous part of the nonlinear spectrum. The reason behind choosing this approach as a preliminary results shows a potential deployment in the future of the optical fiber transmission, reported in [9]. Compared to existing wavelength division multiplexing (WDM) systems, the continuous spectrum has a higher spectral efficiency and reduced complexity at the receiver side.

For the direct NFT numerical implementation, we will investigate in detail the Boffetta-Osborne method, as it cannot be simplified using polynomials, in order to make use of the fast NFT (FNFT), as in [46,47]. An explanation on how this method is improved for NFT, and how it will affect the INFT implementation without performance degradation, will be given below.

For the INFT numerical implementation, we will solve and investigate numerically the GLM integral equations. These equations can represent a complexity that is comparable to FNFT. We have identified the parameters’ values for maximizing the performance, and have used them to compare three different types of modulations that are proposed recently.

In the main strength of our work, we present the complex mathematical tools that are involved in a simplified form, and minimize errors coming from multiple numerical approximations. Being accessible to the NFT research community, this work will help to simulate the optical fiber communication system with high accuracy, avoiding expensive and time consuming experimental setups in laboratories. Compared to the FNFT, the algorithm has higher complexity, but if carefully designed, it can give better performance.

The rest of the paper is organized as follows: in Section 2, we analyze in detail the direct NFT, and explain the numerical procedure to recover the continuous scattering coefficient. In Section 3, we analyze the procedures to be followed, starting from the INFT theory until an accurate numerical simulation of Marchenko system is reached. In Section 4, we validate our numerical algorithms, and in Section 5, we compare the performance of three different types of modulations proposed in the literature. The conclusion is given in the last section.

2. Scattering Coefficient for Nonlinear Fourier Transform

The evolution of the slowly varying complex-valued envelope of the scalar electric field $y(z,t)$ propagating in a single mode optical fiber, neglecting the effects of polarization and attenuation, but taking into account the chromatic dispersion and nonlinear refraction (Kerr effect), is the following:

$$\frac{\partial y(z,t)}{\partial z}=i\frac{{\beta }_2}{2}\frac{{\partial }^{2}y(z,t)}{\partial t^2}-i\gamma {\left|y(z,t)\right|}^{2}y(z,t)$$

(1)

where z is the space coordinate, t is the time reference frame co-moving with the group velocity of the signal envelope, ${\beta }_2\in \mathbb{R}$ chromatic dispersion, and $\gamma >0$ nonlinear parameter. In particular, ${\beta }_2$ represents the dispersion of the group velocity, which causes an enlargement of the impulse during propagation; parameter $\gamma$ models the nonlinear refraction, which is the dependence of the refractive index on the optical field intensity.

To simplify Equation (1), we will first multiply it by the imaginary number i, and then we will make the following substitution: $t\to t_{0}t$, $z\to -\sigma z_{0}z$, $y\to \sqrt{P_0}u$, with $P_0=\frac{\left|{\beta }_2\right|}{\left(\gamma T_0^2\right)}$ in power units, where $T_0$ is a free normalization parameter in units of time, $\sigma =-sgn\left({\beta }_2\right)$, a constant number that represents the nature of the signal inside the optical fiber; the following parameters are all positive quantities: $P_0$, $t_0$, $z_0$ > 0. By selecting all these parameters to have the following relation: $\frac{z_0\left|{\beta }_2\right|}{2t_0^2}=1,z_0\gamma P_0=2$, for example by choosing $P_0=1,z_0=\frac{2}{\gamma },t_0=\sqrt{\frac{\left|{\beta }_2\right|}{\gamma }}$, then we will obtain a simplified form of Equation (1) as the following:

$$i\frac{\partial u}{\partial z}+\frac{{\partial }^{2}u}{\partial t^2}+2\sigma {\left|u\right|}^{2}u=0$$

(2)

when $\sigma =-1$, then the solution envelope of Equation (2) will be solitons. These are usually called deep water waves. In the case of $\sigma =+1$, then the soliton cannot exist and the solution is called shallow water as it is similar to the Korteweg-de Vries (KdV) equation. In this work, we are interested in a deep water waves solution with $\sigma =-1$.

Next, we will explain how the NFT method works. Let us suppose that we have a differential equation that describes the behavior of $u(t,z)$ in (2), and we want to solve it using NFT. First, we have to find two linear operators $\mathcal{L}$ and A that are dependent on $\left(t,z,\left\{u(t,z)\right\}\right)$ (and possibly the derivatives of $u(t,z)$) that will describe a scattering problem of the following system:

$$\left\{\begin{array}{c}\mathcal{L}w=\lambda w\\ w_z=Aw\end{array}\right.$$

(3)

where $\lambda$ is the eigenvalue, w is the eigenfunction, and the subscript denotes the derivative.

The operators $\mathcal{L}$ and A will have to be chosen according to a certain criteria as explained below: we look for an adequate evolutionary equation over the spatial dimension for w, assuming that the Schwarz Theorem applies; that is, $w_{tz}=w_{zt}$. We are going to take the derivative of the first equation in (3) with respect to space, supposing that: ${\lambda }_z=0$. After some simple manipulation, we can find the following relation: $\left|{\mathcal{L}}_z+\left(\mathcal{L}A-A\mathcal{L}\right)\right|w=0$. If we want this to be valid for every w, then we have to impose the following Lax equation [48]:

$${\mathcal{L}}_z+\left(\mathcal{L}A-A\mathcal{L}\right)=0$$

(4)

This equation contains a condition of compatibility for the operators, which in this case are called Lax pairs and can be possibly solved using the PDE for u. To use the NFT, we need to associate Equation (2) with some Lax pairs. Suppose we have the following Lax pairs [49]:

$$\mathcal{L}=\left[\begin{array}{cc}i{\partial }_t& -iu(t,z)\\ -i\sigma u^{*}\left(t,z\right)& -i{\partial }_t\end{array}\right]$$

(5)

$$A=\left[\begin{array}{cc}2i{\partial }_t^2+i\sigma {\left|u\right|}^2& -2iu{\partial }_t-i{u}_t\\ -2i\sigma u^{*}{\partial }_t-i\sigma u_t^{*}& -2i{\partial }_t^2-i\sigma {\left|u\right|}^2\end{array}\right]$$

(6)

The first differential operator in (5) gives us the associated eigenvalue problem to the NLSE, while the second one is the temporal evolution of the eigenfunction, as explained in the following theorem.

Theorem 1.

The Lax equation given in Equation (4) is equivalent to the NLSE equation in Equation (2) on the hypothesis that the following conditions are true: $\mathcal{L}\Phi =\lambda \Phi$, ${\Phi }_z=A\Phi$, ${\lambda }_z=0$, where $\Phi =\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right]$ is an autofunction of two components.

Proof of Theorem 1.

We first apply the operator A to $\Phi$ as

$${\Phi }_z=A\Phi =\left(\genfrac{}{}{0pt}{}{2i{\varphi }_{1,tt}+i\sigma {\left|u\right|}^2{\varphi }_1-2iu{\varphi }_{2,t}-i{u}_t{\varphi }_2}{-2i\sigma u^{*}{\varphi }_{1,t}-i\sigma u_t^{*}{\varphi }_1-2i{\partial }_t^2{\varphi }_{2,tt}-i\sigma {\left|u\right|}^2{\varphi }_2}\right)$$

(7)

and after this, we apply $\mathcal{L}$ to obtain two components: $\mathcal{L}A{\Phi }_1$ and $\mathcal{L}A{\Phi }_2$. We perform the same thing, but this time we first apply $\mathcal{L}$ and then A to obtain the first component $A\mathcal{L}{\Phi }_1$. Finally, we can find: $\left(\mathcal{L}A-A\mathcal{L}\right){\Phi }_1=-u_{tt}{\varphi }_2-2\sigma {\left|u\right|}^{2}u{\varphi }_2$ and $\left\{{\mathcal{L}}_z+\left(\mathcal{L}A-A\mathcal{L}\right)\right\}{\Phi }_1=-\left(i{u}_z+u_{tt}+2\sigma {\left|u\right|}^{2}u\right){\varphi }_2=0$, which is the NLSE in Equation (2). □

After choosing the Lax pair, we need to analyze the direct NFT associated with the following equation:

$$\mathcal{L}\Phi =\lambda \Phi$$

(8)

For the NFT method to work, we need to impose that at infinity, our solution u will decay to zero, $u\to 0$ for $\left|t\right|\to +\infty$. Equation (8) will be approximated with the following result: $i\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]{\Phi }_t=\lambda \Phi ,$ with ${\varphi }_{1,t}=-i\lambda {\varphi }_1$ and ${\varphi }_{2,t}=+i\lambda {\varphi }_2$, whose solutions are simple to calculate.

To form the basis of the solutions space of Equation (8), we can choose from the following Jost functions, defined as limit [48]:

$$\varphi \left(t,\lambda \right)\sim \left(\genfrac{}{}{0pt}{}{1}{0}\right)e^{-i\lambda t},\overline{\varphi }\left(t,\lambda \right)\sim \left(\genfrac{}{}{0pt}{}{0}{1}\right)e^{i\lambda t},fort\to -\infty$$

(9)

$$\phi \left(t,\lambda \right)\sim \left(\genfrac{}{}{0pt}{}{0}{1}\right)e^{i\lambda t},\overline{\phi }\left(t,\lambda \right)\sim \left(\genfrac{}{}{0pt}{}{1}{0}\right)e^{-i\lambda t},fort\to +\infty$$

(10)

The overbar sign means that two autofunctions are in straight relation with each other. For example, the existence of $\varphi \left(t,{\lambda }^{*}\right)$ will imply also that of $\overline{\varphi }\left(t,\lambda \right)$, given that if the first function is a solution to $\mathcal{L}\Phi =\lambda \Phi$, then the second will also be a solution that is linearly independent from the first one:

$$\overline{\varphi }\left(t,\lambda \right)=\left[\begin{array}{cc}0& -\sigma \\ 1& 0\end{array}\right]{\varphi }^{*}\left(t,\lambda \right)$$

(11)

Proposition 1.

The four solutions ϕ, $\overline{\varphi }$, φ, $\overline{\phi }$ are independent in pairs.

Proof of Proposition 1.

The Wronskian of two autofunctions can be easily found to be 1 as: $W\left(\varphi ,\overline{\varphi }\right)={lim}_{x\to \infty }W\left(\varphi ,\overline{\varphi }\right)=e^{-i\lambda t}{e}^{i\lambda t}=1$. □

The equation $\mathcal{L}\Phi =\lambda \Phi$ should have only two linear independent solutions, and so $\varphi$ and $\overline{\varphi }$ should be linearly dependent with respect to the pair $\phi$ and $\overline{\phi }$, as follows:

$$\varphi \left(t,\lambda \right)=b\left(\lambda \right)\phi \left(t,\lambda \right)+a\left(\lambda \right)\overline{\phi }\left(t,\lambda \right)$$

(12)

$$\overline{\varphi }\left(t,\lambda \right)=\overline{a}\left(\lambda \right)\phi \left(t,\lambda \right)+\overline{b}\left(\lambda \right)\overline{\phi }\left(t,\lambda \right)$$

(13)

This is true when $\lambda \in \mathbb{R}$, as the four functions are limited, and in that case, these equations will define the relative four scattering coefficients data denoted by $a\left(\lambda \right)$, $b\left(\lambda \right)$, $\overline{a}\left(\lambda \right)$, $\overline{b}\left(\lambda \right)$. We note here that these scattering coefficients are time invariant, a property that will be used for transmitting data over the optical fiber.

To find the spatial evolution of the scattering data, under the hypothesis that: $\mathcal{L}\Phi =\lambda \Phi$ holds for every z, the operator A is such that:

$${\Phi }_z=A\Phi =\left(\genfrac{}{}{0pt}{}{-2i{\lambda }^2+i\sigma {\left|u\right|}^{2}2u\lambda +i{u}_t}{-2i\sigma u^{*}\lambda +i\sigma u_t^{*}2i{\lambda }^2-i\sigma {\left|u\right|}^2}\right)\Phi$$

(14)

which is the matrix representation that we will use to find the spatial evolution of the scattering data. The new Jost functions will be defined now as:

$$\Phi \left(t,z;\lambda \right)=e^{-2i\lambda z}\varphi \left(t,z;\lambda \right),\overline{\Phi }\left(t,z;\lambda \right)=e^{+2i\lambda z}\overline{\varphi }\left(t,z;\lambda \right)$$

(15)

$$\Psi \left(t,z;\lambda \right)=e^{+2i\lambda z}\phi \left(x,t;\lambda \right),\overline{\Psi }\left(t,z;\lambda \right)=e^{-2i\lambda z}\overline{\phi }\left(t,z;\lambda \right)$$

(16)

They are valid for every z, they will solve equation Equation (8), and they will satisfy both Equations (9) and (10).

Let’s take the spatial derivative of Equation (12) with the new Jost function, and after some manipulation, we obtain the following relations:

$$a_z(\lambda ,z)=0,b_z\left(\lambda ,z\right)=4i{\lambda }^{2}b\left(\lambda ,z\right)$$

(17)

$$\overline{a_z}(\lambda ,z)=0,\overline{b_z}\left(\lambda ,z\right)=-4i{\lambda }^2\overline{b}\left(\lambda ,z\right)$$

(18)

We can write the spatial evolution of the continuous scattering data, starting from the initial ones as follows:

$$a\left(\lambda ;z\right)=a\left(\lambda ;0\right),\overline{a}\left(\lambda ;z\right)=\overline{a}\left(\lambda ;0\right)$$

(19)

$$b\left(\lambda ;z\right)=b\left(\lambda ;0\right)e^{+4i{\lambda }^{2}z},\overline{b}\left(\lambda ;z\right)=\overline{b}\left(\lambda ;0\right)e^{-4i{\lambda }^{2}z}$$

(20)

When the condition $a\left(\lambda \right)=0$ is true, then we can define the discrete spectrum of the Lax operator, which consists of a finite number of discrete eigenvalues. The discrete scattering coefficient corresponding to these discrete eigenvalues can be computed by using Equations (12) and (13). In this work, we will focus only on the continuous parameter $\lambda$ and its associated continuous scattering coefficient. Altogether, they will define our nonlinear spectrum of the signal $u(t,z)$.

Numerical Method for the Continuous Scattering Coefficient

We have developed a numerical algorithm to find the continuous scattering coefficient of the Lax equation: $\mathcal{L}\Phi =\lambda \Phi ,$ with the initial condition signal $u(t,0)$. We followed the idea introduced in the article by Boffetta and Osborne [36].

The algorithm will use a discretization of $u(t,0)$ in $2M+1$ samples as: $u_n=u\left(t_n,0\right)=u\left(n\Delta t,0\right)=u(n\frac{T}{M},0)$, where $2T$ is the total width of the function u centered at $t=0$, such that $u\left(t,0\right)=0$ for $\left|t\right|\ge T$ and $M>0$.

The steps followed to find out the numerical algorithm for the continuous scattering coefficient when $\lambda \in \mathbb{R}$ are given below:

Step 1. We need to solve numerically the Lax equation: $\mathcal{L}\Phi =\lambda \Phi$, and find an equation that provides us with a relationship between $\varphi \left(t,\lambda \right)$ and $\varphi \left(-t,\lambda \right)$ for $\left|t\right|\ge T$.

Instead of solving it directly, we will first assume that $u\left(t,0\right)=u_n$ is constant in the interval $\left[t_n-\frac{\Delta t}{2},t_n+\frac{\Delta t}{2}\right]$. After this, we are going to associate to it the equation, ${\Phi }_t=B_n\Phi$, which is equivalent to $\mathcal{L}\Phi =\lambda \Phi$. $B_n$ is an approximation of $\mathcal{L}$ with constant step. We are going to solve numerically the following Cauchy problem, in an iterative way:

$$\left\{\begin{array}{c}{\Phi }_t=B_n\Phi \\ \Phi \left(t_n-\frac{\Delta t}{2},\lambda \right)={\Phi }_n\end{array}\right.,t\in (t_n-\frac{\Delta t}{2},t_n+\frac{\Delta t}{2})$$

(21)

with

$$B_n\left(\lambda \right)=\left[\begin{array}{cc}-i\lambda & u_n\\ -\sigma u_n^{*}& i\lambda \end{array}\right],\forall t\in \left(t_n-\frac{\Delta t}{2},t_n+\frac{\Delta t}{2}\right)$$

(22)

in order to find the new starting point $\Phi \left(t_{n+1}-\frac{\Delta t}{2},\lambda \right)={\Phi }_{n+1}$, having already previously found ${\Phi }_n$. Solving equation Equation (21) until $t_{n+1}$ will give us the following iteration:

$$\Phi \left(t_{n+1}-\frac{\Delta t}{2}\right)=\Phi \left(t_n+\frac{\Delta t}{2}\right)=U\left(u_n\right)\Phi \left(t_n-\frac{\Delta t}{2}\right)$$

(23)

where $U\left(u_n\right)$ represents the effect of the $\mathcal{L}$ operator on interval $(t_n-\frac{\Delta t}{2},t_n+\frac{\Delta t}{2}]$. In particular, we have that:

$$U\left(u_n\right)=\exp (\Delta x{B}_n)$$

(24)

An explicit expression of $U\left(u_n\right)$ can be calculated for our case through the definition of the exponential of a matrix:

$$exp\left(\Delta t{B}_n\right)=\sum _{k=0}^{\infty }\frac{\Delta t^k}{k!}{B}_n^k$$

(25)

and taking into account Taylor’s series of the cosh and sinh functions, for the matrix taken into consideration. The obtained result is:

$$U\left(u_n\right)=\left[\begin{array}{cc}cosh\left(q_n\Delta t\right)-\frac{i\lambda }{q_n}sinh(q_n\Delta t)& \frac{u_n}{q_n}sinh(q_n\Delta t)\\ -\frac{\sigma \overline{u_n}}{q_n}sinh(q_n\Delta t)& cosh(q_n\Delta t)+\frac{i\lambda }{q_n}sinh(q_n\Delta t)\end{array}\right]$$

(26)

with $q_n^2=-\sigma {\left|u_n\right|}^2-{\lambda }^2$ being a complex number. If we introduce a new vector as:

$$\Xi \left(t,\lambda \right)=\left[\begin{array}{c}\Phi \\ {\Phi }^{{}^{\prime }}\end{array}\right]$$

(27)

with ${\Phi }^{\prime }=\frac{\partial \Phi }{\partial \lambda }$, Equation (23) will be transformed into:

$$\Xi \left(t_{n+1}-\frac{\Delta t}{2}\right)=T\left(u_n\right)\Xi \left(t_n-\frac{\Delta t}{2}\right)$$

(28)

with

$$T\left(u_n\right)=\left[\begin{array}{cc}U\left(u_n\right)& 0\\ U^{{}^{\prime }}\left(u_n\right)& U\left(u_n\right)\end{array}\right]$$

(29)

We take $t_{-M}$ as the starting point and iterate it until we find the following relation for a generic $t_n$:

$$\Xi \left(t_{n+1}+\frac{\Delta t}{2}\right)=\prod _{i=n}^{-M}T\left(u_i\right)\Xi \left(t_{-M}-\frac{\Delta t}{2}\right)$$

(30)

having the relation:

$$\prod _{i=M}^{-M}T\left(u_i\right)\Xi \left(t_{-M}\right)=\left[\begin{array}{cc}\Sigma & 0\\ {\Sigma }^{{}^{\prime }}& \Sigma \end{array}\right]\Xi \left(t_{-M}-\frac{\Delta t}{2}\right)=S\Xi \left(t_{-M}-\frac{\Delta t}{2}\right)$$

(31)

and define S the scattering matrix as:

$$S=\left[\begin{array}{cc}\Sigma & 0\\ {\Sigma }^{{}^{\prime }}& \Sigma \end{array}\right],\Sigma \left(\lambda \right)=\prod _{i=M}^{-M}U\left(u_i\right)$$

(32)

In Equation (31), we found a relation that connects $\varphi \left(t,\lambda \right)$ and $\varphi \left(-t,\lambda \right)$ for $\left|t\right|\ge T$. In particular, we performed it for the case when $t=T+\frac{\Delta t}{2}$. Therefore, we can move to the second step of our algorithm, starting from this value of t.

Step 2. Assume that the conditions (9) and (10) are satisfied for $\left|t\right|\ge T$, given that $u\left(t,0\right)=0$ for $\left|t\right|\ge T$. The operator $\mathcal{L}$ (depending on $u\left(t,0\right))$ will assume its asymptotic configuration as:

$${\mathcal{L}}_{\infty }=i\left[\begin{array}{cc}{\partial }_t& 0\\ 0& -{\partial }_t\end{array}\right]$$

(33)

and so for every $t>T$, the following relation should hold:

$$\varphi \left(-t,\lambda \right)=\left(\genfrac{}{}{0pt}{}{1}{0}\right)e^{i\lambda t}$$

(34)

For $t=T+\frac{\Delta t}{2}$, we can rewrite Equation (34) as:

$$\varphi \left(-T-\frac{\Delta t}{2},\lambda \right)=\left(\genfrac{}{}{0pt}{}{1}{0}\right)e^{i\lambda (T+\frac{\Delta t}{2})}$$

(35)

Step 3. Given that: $\varphi \left(t,\lambda \right)=b\left(\lambda \right)\phi \left(t,\lambda \right)+a\left(\lambda \right)\overline{\phi }\left(t,\lambda \right)$, we can substitute Equation (10) into it, to find the following relation for the same value of t:

$$\varphi \left(T+\frac{\Delta t}{2},\lambda \right)=b\left(\lambda \right)\left(\genfrac{}{}{0pt}{}{0}{1}\right)e^{i\lambda \left(T+\frac{\Delta t}{2}\right)}+a\left(\lambda \right)\left(\genfrac{}{}{0pt}{}{1}{0}\right)e^{-i\lambda \left(T+\frac{\Delta t}{2}\right)}$$

(36)

Use Equations (27), (31), (35) and (36) to find an expression for $\varphi \left(t,\lambda \right)$. Using simple algebraic steps, we can have a numerical approximation of both the continuous scattering data $a\left(\lambda \right)$ and $b\left(\lambda \right)$ as:

$$a\left(\lambda \right)=S_{11}{e}^{2i\lambda \left(T+\frac{\Delta t}{2}\right)}={\left(\prod _{i=M}^{-M}U\left(u_i\right)\right)}_{11}{e}^{2i\lambda \left(T+\frac{\Delta t}{2}\right)}$$

(37)

$$b\left(\lambda \right)=S_{21}={\left(\prod _{i=M}^{-M}U\left(u_i\right)\right)}_{21}$$

(38)

A summary of the main steps involved in the calculation of continuous scattering data is presented below:

• Determine a good approximation of T so that we can say that $u\left(t,0\right)=0$ for $\left|t\right|\ge T$;
• Descretize the interval $(-T,T)$ in $2M+1$ samples $t_n=n\Delta t=n\frac{T}{M}$, with $n=-M,-M+1,\dots M-1,M;$
• Find ${\prod }_{i=M}^{-M}U\left(u_i,\lambda \right)$ for the necessary $\lambda$;
• Use Equations (37) and (38) to find the continuous scattering data for the chosen $\lambda$;

3. The Inverse Nonlinear Fourier Transform Method

In following section, we will use the scattering data found in the previous section to find its relationship with the solution $u(t,z)$ of the NLSE equation. We will first give without proof the following theorem that will be an important step towards the Gel’fand-Levitan-Marchenko equation [49].

Theorem 2.

The solutions$\phi \left(t,\lambda \right)$ and $\overline{\phi }\left(t,\lambda \right)$ of Lax Equation (8) also satisfying (10), will be represented as follows:

$$\phi \left(t,\lambda \right)=\left(\genfrac{}{}{0pt}{}{0}{1}\right)e^{i\lambda t}+{\int }_t^{\infty }K(t,s)e^{i\lambda t}ds$$

(39)

$$\overline{\phi }\left(t,\lambda \right)=\left(\genfrac{}{}{0pt}{}{1}{0}\right)e^{-i\lambda t}+{\int }_t^{\infty }\overline{K(t,s)}{e}^{-i\lambda t}ds$$

(40)

with $K(t,s)$ and $\overline{K(t,s)}$ being independent from λ, such that:

$$K\left(t,s\right)=-\frac{1}{2}u\left(t\right)$$

(41)

$$\overline{K(t,s)}=\frac{1}{2}u{\left(t\right)}^{*}$$

(42)

In order to obtain the Gel’fand-Levitan-Marchenko equation, we have to consider that $u\left(t\right)$ should decay fast enough in order to extend Equation (12) to the upper part of the complex plane in such a way as to satisfy the following relation of the Jost functions and continuous data scattering:

$$\frac{1}{a\left(\lambda \right)}\varphi \left(t,\lambda \right)=\frac{b\left(\lambda \right)}{a\left(\lambda \right)}\phi \left(t,\lambda \right)+\overline{\phi }\left(t,\lambda \right)$$

(43)

If we substitute (39) and (40) in (43) and we also define the function $F\left(t\right)=\frac{1}{2\pi }{\int }_{}^{}\frac{b\left(\lambda \right)}{a\left(\lambda \right)}{e}^{i\lambda t}d\lambda$, we can easily recover the Gel’fand-Levitan-Marchenko equation given as:

$$K\left(t,y\right)+F^{*}\left(t+y\right)-\int {\int }_t^{\infty }K\left(t,r\right)F\left(r+s\right)F^{*}\left(s+y\right)drds=0$$

(44)

The important case studied in this paper is when the discrete scattering data do not exist, $\left(a\left(\lambda \right)\ne 0,\mathfrak{I}\left(\lambda \right)\ge 0\right)$, and $F\left(t\right)$ is defined as the inverse Fourier transform of $\frac{b\left(\lambda \right)}{a\left(\lambda \right)}$. The ratio $\frac{b\left(\lambda \right)}{a\left(\lambda \right)}$ is also called the reflection coefficient or spectral amplitude of the nonlinear spectrum. Fast Fourier transform (FFT) can be used to compute it.

Making the following variable change in (44), $\beta =y-t$ will transform the Gel’fand-Levitan-Marchenko equation into two equations called Marchenko systems, which will highly simplify our numerical simulation. If we define: $B_2\left(t,0;z\right)=-2u(t,z)$ and $B_1\left(t,\beta \right)={\int }_0^{\infty }{B}_2\left(\beta ,r\right)F\left(r+\beta +2t\right)dr$ as our two auxiliary functions, we can substitute Equation (44) equivalently with the following Marchenko system:

$$\left\{\begin{array}{c}{B}_1\left(t,\beta \right)={\int }_0^{\infty }{B}_2\left(t,r\right)F\left(r+\beta +2t\right)dr\\ B_2\left(t,\beta \right)={\int }_0^{\infty }{B}_2\left(t,r\right)F^{*}(r+\beta +2t)dr-F^{*}(2t+\beta )\end{array}\right.$$

(45)

Equation (45) should be now solved numerically. Observe that in this way the quadrature formulas that we will use for the approximation of the two integrals that will be applied on the same interval and on the same variable. Note that to apply such quadrature formulas, it will be necessary to approximate the domain of function F, as explained in the next subsection.

3.1. Procedure for the Numerical Simulation of $F(t;z)$

The purpose of this section is to find a procedure to numerically simulate Equation (45). We need first to define the domain of $F(t;z)=\frac{1}{2\pi }{\int }_{}^{}\frac{b(\lambda ;z)}{a(\lambda ;z)}{e}^{i\lambda t}d\lambda$. FFT will be used in this regard, together with some sampling constraints to satisfy the sampling theorem. This can be easily recognized inside the integral with some minor changes.

Let us suppose that the domain of $\frac{b(\lambda ;z)}{a(\lambda ;z)}$ will be in a range between $[{\lambda }_{Max},{\lambda }_{Min}]$ with $N\in \mathbb{N}$ total discrete points where we have to apply FFT. We denote these discrete points as: ${\lambda }_k=k\Delta \lambda =\frac{k({\lambda }_{Max}-{\lambda }_{Min})}{N}$, with $k=-\frac{N}{2},\dots ,\frac{N}{2}-1$ and $\Delta \lambda$ frequency resolution. In a strict relation with these points, we will have the non-zero time domain of $F(t;z)$ where we will define the sampling points as: $T_n=n\Delta T$ with $n=-\frac{N}{2},\dots ,\frac{N}{2}-$ and $\Delta T$ time resolution. The number of points for both domains $F(t;z)$ and $\frac{b(\lambda ;z)}{a(\lambda ;z)}$ will be the same in this work. Our scope is to find the relationship between ${\lambda }_k$ and $T_n$ for both the time and frequency domains.

Starting with the smallest frequency associated with the interval $[{\lambda }_{Max},{\lambda }_{Min}]$ that will be related to $\Delta T$ as: $\Delta T=\frac{1}{{\lambda }_{Max}-{\lambda }_{Min}}$, we also have $N\Delta T=\frac{N}{{\lambda }_{Max}-{\lambda }_{Min}}$ and $\Delta T\Delta \lambda =\frac{1}{N}$.

In order to solve Equation (44), we developed the following procedure:

• Approximate the values ${\lambda }_{Max},{\lambda }_{Min}$, such that $\frac{b(\lambda ;z)}{a(\lambda ;z)}=0$ if $\lambda <{\lambda }_{Min}$ and $\lambda >{\lambda }_{Max}$ to estimate the domain for $\frac{b(\lambda ;z)}{a(\lambda ;z)}$.
• Control if $2\Delta T<\frac{2\pi }{{\lambda }_{Max}-{\lambda }_{Min}}$.

  (a)

  If True, define first $p=\frac{2\pi }{{\lambda }_{Max}-{\lambda }_{Min}}-2\Delta T$ and perform the following change in the variables ${\lambda }_{Min}$ and ${\lambda }_{Max}$:

  ${\lambda }_{Min}={\lambda }_{Min}-\frac{p}{2}$ and ${\lambda }_{Max}={\lambda }_{Max}-\frac{p}{2}$.

  (b)

  If False, divide by 2 the step $\Delta t$ until point (a) is satisfied.

• Select the resolution in frequency based on the total number of points as:

  $$\Delta \lambda =\frac{{\lambda }_{Max}-{\lambda }_{Min}}{N}$$

  (46)
• Select the resolution in time as:

  $$\Delta T=\frac{1}{{\lambda }_{Max}-{\lambda }_{Min}}$$

  (47)
• Find the spatial evolution of the spectral amplitude by sampling $\frac{b(\lambda ;z=0)}{a(\lambda ;z=0)}{e}^{4i{\lambda }^{2}z}$ on the grid constructed for $\frac{b(\lambda ;z=0)}{a(\lambda ;z=0)}$, and apply the inverse FFT in order to find the vector $F(t;z)$ on the grid with amplitude:

  $$N\Delta T=\frac{N}{{\lambda }_{Max}-{\lambda }_{Min}}$$

  (48)
• Select a user-defined variable called tolerance, $\epsilon$, and check if the function $F(t;z)$ is less than it at the edge only, in order to better approximate its domain. We will differentiate between two cases:

  (a)

  If $max\left[P\left(\pm \frac{N}{2}\Delta T\right)\right]<\epsilon$, then the approximated support and the corresponding grid are chosen such that we can apply the quadrature formulas approximating the integrals to be computed in the Marchenko system of Equation (3.5).

  (b)

  Otherwise, we multiply by 2 the total number of points and restart from the second step of this numerical algorithm.

Following this procedure, we will be able to properly represent the signal $u\left(t,z\right)$ and the integral of the function $F(t;z)$.

3.2. The Numerical Method Algorithm for the Marchenko System

Given the above considerations, the points on which we will decide to approximate the function $u\left(t,z\right)$ will be equally spaced at a distance $\Delta T$. The points where we will find $u\left(t,z\right)$ are collected in the following vector: $T=\left[t_1=a,\dots ,t_n=a+\left(n-1\right)\Delta T,\dots ,t_N=b\right],n=1,\dots ,N$, where $a<0$ and $b>0$ are extreme points of this vector.

To solve the Marchenko system in Equation (45), we are going to use the Nyström method, approximating the integrals by means of a quadrature formula to derive a discrete system to be solved through an appropriate method. In particular, we want to find $B_2\left(t,\beta \right)$ for $\beta =0$. The system in (45) with $\beta =0$ will be:

$$\left\{\begin{array}{c}{B}_1\left(t_n,0\right)={\int }_0^{\infty }{B}_2\left(t,r\right)F\left(r+2t_n\right)dr\\ B_2\left(t_n,0\right)={\int }_0^{\infty }{B}_1\left(t,r\right)F^{*}(r+2t_n)dr-F^{*}\left(2t\right)\end{array}\right.$$

(49)

Equation (49) should be solved numerically. We observe that inside the integral, we need to evaluate $B_2\left(t,r\right)$ and $B_1\left(t,r\right)$ for values of higher than zero. As a result, we need to solve for the values of $\beta \ne 0$. In particular, we are going to take $\beta$ and r from the same domain $D={\left[2m\Delta t\right]}_{m=0}^l$ where $l=ceil\left[\frac{N}{2}-n\right]$ is the first value of m, such that we are outside the domain of F. The reason for this choice is that such a domain should take into account the fact that we know the support of F, which is evaluated in both integrals as $r+2t_n=r+2n\Delta t$. The sampling step of variable r where we will integrate is: $\Delta \beta =2\Delta t$. Now, we can finally approximate both integrals in Equation (49) using the Simpson rule to find the system of equations:

$$\left\{\begin{array}{c}{B}_1\left(t_n,r_j\right)=\sum _{m=0}^l{d}_m{B}_2\left(t_n,r_m\right)F\left(r_{n+j+m}\right)\\ B_2\left(t_n,r_j\right)=\sum _{m=0}^l{d}_m{B}_1\left(t_n,r_m\right)F^{*}\left(r_{n+j+m}\right)-F^{*}\left(r_{n+j}\right)\end{array}\right.$$

(50)

where $r_{n+j+m}=\left(k+j+l\right)\Delta \beta$, $d_m=\frac{\Delta \beta }{3}\{1,4,2,4,2..\}$, $j=0,\dots ,l$. The system in Equation (50) can be represented using a matrix notation such as:

$$\left|\begin{array}{cc}I& -\mathcal{F}D\\ {\mathcal{F}}^{*}D& I\end{array}\right|\left|\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{B}}_2\end{array}\right|=\left|\begin{array}{c}0\\ {\mathcal{F}}^{*}\end{array}\right|$$

(51)

where $\mathcal{F}={\left(F\left(r_{n+j+m}\right)\right)}_{m,j=0}^l$, ${\mathcal{F}}^{*}={\left(F\left(r_{n+j}\right)\right)}_{j=0}^l$, $D=diag{\left\{d_m\right\}}_{m=0}^l$, are known variables. The unknown variables are:

$${\mathbf{B}}_1={\left(B_1\left(t_n,r_m\right)\right)}_{m=0}^l,{\mathbf{B}}_2={\left(B_2\left(t_n,r_m\right)\right)}_{m=0}^l$$

(52)

In Equation (51), we are interested in the first component of ${\mathbf{B}}_2$. After finding ${\mathbf{B}}_1=\mathcal{F}D{\mathbf{B}}_2$ from the first part of the system, we substitute it in the second part, and then solve the following system:

$$\left({\mathcal{F}}^{*}D\mathcal{F}D+I\right){\mathbf{B}}_2={\mathcal{F}}^{*}$$

(53)

Once we have solved it, we are able to find out the signal propagated at a desired distance z:

$$u\left(t,z\right)={\mathbf{B}}_2\left(r_0\right)$$

(54)

In summary, to find out the solution $u\left(t,z\right)$ of Equation (2) through the implementation of the INFT (in this case, when $a\left(\lambda \right)\ne 0$ for every $\lambda \in \mathbb{C}$), we need to follow the steps below:

• Based on the above algorithm, we first need to find the scattering data $a\left(\lambda \right)$ and $b\left(\lambda \right)$, and their ratio $\frac{b\left(\lambda \right)}{a\left(\lambda \right)}$.
• We need to find the spatial evolution of continuous spectral amplitudes $\frac{b\left(\lambda \right)}{a\left(\lambda \right)}$ by multiplying them with $e^{4i{\lambda }^{2}z}$ (step 2 should be applied only in case we will transmit the signal in an optical fiber described by Equation (2)).
• Find $F\left(t,z\right)=\frac{1}{2\pi }{\int }_{}^{}\frac{b\left(\lambda \right)}{a\left(\lambda \right)}{e}^{i\lambda t}d\lambda$ using IFFT, as explained above.
• Solve the system in (53) to find out $u\left(t,z\right)$ given in (54).

4. Validation of the NFT and INFT Algorithm Implementations

To verify our algorithm implementation for both the direct and inverse nonlinear Fourier methods, we have run multiple tests on soliton pulses taken as the initial condition with a well-known analytically closed form solution. The initial condition of our simulation is the hyperbolic secant function shown also in Figure 1:

$$u\left(t,0\right)=e^{it}sech\left(t\right)$$

(55)

We can theoretically find the closed form solution of the scattering data, as in [50]:

$$a\left(\lambda \right)={\varphi }_1(t,\lambda )e^{i\lambda },b\left(\lambda \right)={\varphi }_2(t,\lambda )e^{i\lambda }$$

(56)

$$\varphi =\left[\begin{array}{c}-\frac{\left(i/2\right)tanh\left(t\right)+\lambda +1/2}{i/2+\lambda +1/2}{e}^{-i\lambda t}\\ -\frac{i/2}{\left(i/2+\lambda +1/2\right)cosh\left(t\right)}{e}^{-i(\lambda +1)t}\end{array}\right]$$

(57)

In Figure 2 and Figure 3, we show the continuous scattering data obtained from Equation (55) using the NFT algorithm described in the previous sections.

The difference between the analytical continuous scattering, $a\left(\lambda \right)$, and the numerical one, $\widehat{a}\left(\lambda \right)$, is measured in terms of the normalized means squared errors (NMSEs), defined as:

$$NMSE=\frac{{\int }_{-\infty }^{+\infty }{\left|a\left(\lambda \right)-\widehat{a}\left(\lambda \right)\right|}^{2}d\lambda }{{\int }_{-\infty }^{+\infty }{\left|a\left(\lambda \right)\right|}^{2}d\lambda }$$

(58)

We compare the numerical results obtained with the proposed modified algorithm similiar to [36] and theoretical one given in Equations (56) and (57), for different numbers of discrete points $M,$ as shown in Figure 4. We have assumed a truncated soliton pulse in the range $t\in \left[-6,6\right]$, and a spectral parameter $\lambda \in \left[-10,10\right]$. The NMSE is used to measure the differences for both scattering coefficients $a\left(\lambda \right)$ and $b\left(\lambda \right)$ in order to understand how much the numerical simulations will deviate from the analytical one as a function of M. If $M=128$, it means that $u\left(t,0\right)$ is discretized in $2\ast M+1=257$ points. The results are highly correlated to M, as can be observed from Figure 4. If M will increase, the NMSE will also increase more and more between the numerical and theoretical results.

For further verification of our numerical algorithm procedure for the calculation of scattering data, we have checked if the following condition is satisfied (which we know that it should be always true):

$${\left|a\left(\lambda \right)\right|}^2+{\left|b\left(\lambda \right)\right|}^2=1$$

(59)

For a fixed $M=128$, we have obtained an error that is approximately ${10}^{-10}$, that is: ${\left|a\left(\lambda \right)\right|}^2+{\left|b\left(\lambda \right)\right|}^2-1\approx {10}^{-10}$; consequently we can conclude that the proposed methodology has been well implemented and that it is working.

To verify for $F(t;0)$, we recall that $u\left(t,0\right)$ can be found from (41) and (44), and assuming $z=0$ and $y=t$, we have that: $-2K\left(t,0\right)=u\left(t\right)$ and $K\left(t,t\right)+F^{*}\left(2t\right)-{\int \int }_t^{\infty }$ $K\left(t,r\right)F\left(r+s\right)F^{*}\left(s+t\right)drds=0$. When $u\left(t,0\right)$ is very small, we can assume that K and F are also small. It follows that the third addendum of Equation (44) can be negligible and we can consider: $\left|u\left(t,0\right)\right|\approx 2\left|F(2t,0)\right|$.

Assuming as an initial condition: $u\left(t,0\right)={10}^{-3}{e}^{-x^2}$, we can simply use the FFT to find out $F\left(t\right)$, taking into account strictly following the algorithm shown in the section above. The domain of $\lambda$ is selected based on the maximum value that we have observed from the ratio of $b\left(\lambda \right)/a\left(\lambda \right)$, which is in the order of ${10}^{-3}$, such that: $\left|b\left(\lambda \right)/a\left(\lambda \right)\right|>{10}^{-6}$ $\forall \lambda \in \left[{\lambda }_{min},{\lambda }_{max}\right]$. In addition, to verify the hypothesis according to which $\left|u\left(t,0\right)\right|\approx 2\left|F(2t,0)\right|$ and knowing that the maximum of $u\left(t,0\right)$ is ${10}^{-3}$, we can assume that the domain of $F\left(t\right)$ is the interval where $F\left(t\right)>{10}^{-6}$. We are going to look for the number of points N on which to calculate $b\left(\lambda \right)/a\left(\lambda \right)$ and $F\left(t\right)$ based on the above observation.

In Figure 5, we show the result of both functions |u(t, 0)| and 2|F(2t, 0)| with N = 128, which looks very similar to each other. Selecting the right number of points N is very important, as the two functions will not look so similar to each other in the case where N is lower than 128.

5. Applications of NFT and INFT on an Optical Communication System

In this section, we will apply our proposed implementation of the NFT/INFT algorithms on a modulated signal used in optical communications. The information is encoded in continuous scattering data at the transmitter side and is decoded on the receiver side using both NFT and INFT numerical methods (as explained in the previous sections). We call this a naive implementation as it is in principle simple, but it has some disadvantages that we are going to note below. The second modulation that we will investigate is similar to the previous one, with the difference of better controlling the pulses in the time and frequency domains. Lastly, we will apply our NFT/INFT algorithm to the recently proposed method called b-modulation, where only $b\left(\lambda \right)$ will be used to encode the information at the transmitter side.

The nonlinear Fourier transform transmission system will exploit its orthogonality property to allow different users to transmit at the same time in optical fiber medium without any interference with each other. This is similar to the traditional OFDM that is used in almost all new wireless and wired technologies. The scenario that we are going to investigate will involve only the continuous part of the nonlinear spectrum, as the discrete part will increase the complexity at the receiver side. The use of a root searching algorithm is needed, and it is difficult to control their number. In a numerical simulation, we will make sure that there is no discrete part of the nonlinear spectrum.

The first modulation where we are going to test our implementation of NFT/INFT algorithm is originally proposed in [16]. Suppose we have the following quadrature phase shift keying (QPSK) modulated signal expressed in the time domain as:

$$s\left(t,0\right)=\sum _k^{}{s}_{k}p(t-kT)$$

(60)

where $s_k\in \left\{1\pm i,-1\pm i\right\}$ complex symbols, $T=0.02\left[ns\right]$ is the pulse duration, and $p\left(t\right)$ is the root raised cosine function with unit energy and roll off factor $\rho =0.2$. Let us also define: $S\left(f\right)=\mathcal{F}\left\{s\left(t\right)\right\}$, the linear Fourier transform of our modulated signal. In the first modulation, the linear spectrum $S\left(f\right)$ is mapped directly to the continuous nonlinear spectrum defined as the ratio between scattering data, called the reflection coefficient: ${\eta }_c\left(\lambda ,0\right)=\frac{b\left(\lambda ,0\right)}{a\left(\lambda ,0\right)}=-\frac{1}{2\pi }S\left(-\frac{\lambda }{\pi }\right)$. In order to compute the numerical sample in the time domain, we use the INFT algorithm as explained in the previous sections. The transmitted signal in the optical fiber will be:

$${\eta }_c\left(t,0\right)=INFT\left({\eta }_c\left(\lambda ,0\right)\right)$$

(61)

At the receiver side, the NFT algorithm is used to obtain:

$${\eta }_c\left(\lambda ,z\right)=NFT\left({\eta }_c\left(t,z\right)\right)$$

(62)

The channel equalization is also applied to obtain the transmitted non linear spectrum as:

$${\eta }_c\left(\lambda ,0\right)={\eta }_c\left(\lambda ,z\right)e^{-4j{\lambda }^{2}z}$$

(63)

The equalization is performed by inverting the effect of the channel, which is simply a phase shift in the nonlinear Fourier transform. An inverse Fourier transform followed by a match filter are used to recover the received symbols:

$$\widehat{s_k}={\int }_{-\infty }^{\infty }s\left(t,z\right)p^{*}\left(t-kT\right)dt$$

(64)

where $\widehat{s_k}$ is an estimation of $s_k$ in case the channel is affected by white Gaussian noise.

In the second modulation that we are going to investigate, the reflection coefficient ${\eta }_c\left(\lambda ,0\right)$ is not obtained directly from mapping the linear Fourier transform, but from the following transformation to ensure that the signals have the same energy in linear and nonlinear space:

$${\eta }_{c,n}\left(\lambda ,0\right)={\left(e^{{\left|S\left(\lambda ,0\right)\right|}^2}-1\right)}^{1/2}{e}^{jarg\left(S\left(\lambda ,0\right)\right)}$$

(65)

The third modulation to test our algorithm is called b-modulation. The data are encoded only on $b\left(\lambda ,0\right)$ as part of the nonlinear spectrum with the following transformation:

$$b\left(\lambda ,0\right)=AS\left(\lambda ,0\right)$$

(66)

where:

$$A=\frac{\sqrt{1-e^{-{\left|S\left(\lambda ,0\right)\right|}^2}}}{\sqrt{{\left|S\left(\lambda ,0\right)\right|}^2}}$$

(67)

The continuous scattering data $a\left(\lambda ,0\right)$ are related to $b\left(\lambda ,0\right)$ as:

$$a\left(\lambda ,0\right)=\sqrt{1-{\left|b\left(\lambda ,0\right)\right|}^2}{e}^{j\mathcal{H}\left(\frac{log\left(1-{\left|b\left(\lambda ,0\right)\right|}^2\right)}{2}\right)}$$

(68)

where $\mathcal{H}(·)$ denotes the Hilbert transform. The reflection coefficient will be denoted as ${\eta }_{c,b}\left(\lambda ,0\right)$. The details of b-modulation are already discussed in several other papers [30,51,52].

We conducted different simulation in MATLAB to transmit the three modulated signals in an optical fiber medium of total length 1600 km. The fiber was divided into 20 spans of 80 km each, where at the end of each span, white noise originating from distributed Raman amplification was added, with a noise figure of 4 dB. The optical fiber dispersion parameter is ${\beta }_2$ = −20.4 ps²/km, and the nonlinear coefficient is $\gamma$ = 1.3 W${}^{-1}$ km${}^{-1}$. At the receiver side, we are going to measure the performance using the Q-factor in dB found directly from the constellation affected by the white noise, as in [53]:

$$Q=10{log}_{10}{\left(\frac{{\sum }_k^{}{\left|s_k-\widehat{s_k}\right|}^2}{{\sum }_k^{}{\left|s_k\right|}^2}\right)}^{-1}$$

(69)

For historical reasons, the performance in optical systems are usually measured in terms of the quality Q-factor, given as: $Q=\sqrt{2}erf{c}^{-1}\left(2BER\right)$, where $erf{c}^{-1}$ is the inverse complementary error function, and $BER$ is the bit error rate. For a fixed $BER={10}^{-3}$, Q-factor is equal to 9.8 dB. Systems with a Q-factor of less than 9.8 dB will fail to work, as the number of errors is very large. We can think of this value as a lower bound for our system to work. The reason for why we are using Equation (69) as a quality parameter is because the number of transmitted symbols in our numerical simulation is only 1024, and the minimum $BER$ that we can achieve is only $4.8\times {10}^{-4}$. To fully investigate the system, even for lower $BER$, we use Equation (69). Because the oversampling factor used in this work is 40, we cannot increase the number of symbols to achieve a lower $BER$, as the computation cost will be very high.

The local oscillator at the receiver side has no frequency drift and has a linewidth of 100 kHz. The total number of symbols transmitted was $2^{10}$. We added $25\%$ of overhead to satisfy Step 1 of the direct NFT algorithm. An oversampling factor of 40 was selected for the simulations. The system performance for the three different modulations applying our algorithm are shown in Figure 6 and Figure 7. In Figure 6, we show the Q-factor versus the average power in the input of the optical fiber for all three modulations. We observe that the b-modulation performs better compared to the other two modulations, which on the other hand, show very similar performances. The maximum value of the Q-factor, Q = 25.9 dB, is achieved when the average power is equal to −1 dBm, for the b-modulation shown in the blue line in Figure 6. Compared to the other two modulations, the Q-factor is increased by 2 dB. If the average power is less than −1 dBm, the performance is deteriorating as the Q-factor is decreasing. The reason behind this is the noisy channel assumed in the system, which is added by the distributed optical Raman amplifiers assumed over the optical links. If the average power is greater than −1 dBm, then the Q-factor is also decreasing due to the nonlinear effects and the numerical errors from many approximation values in our method.

In Figure 7, we show the Q-factor versus transmission distances for all three modulations. We also observe better performance when b-modulation is used, shown in the blue line, and similar performances for the other two modulations. We observe that in the first 400 km of transmission distance, b-modulation has a lower Q-factor with respect to the other two. After 500 km of transmission distance, b-modulation start to shows a better Q-factor for all of the distances that we have simulated.

6. Conclusions

In this work, we have improved the NFT and INFT algorithms by going into the numerical details of the integral equation, proposing a novel methodology for their implementation, validating the results, and applying them on three different modulated signals for an optical communication system.

Only the continuous part of the nonlinear scattering data is investigated, as the spectral efficiency is higher and the complexity implementation is lower compared to the discrete nonlinear spectrum. We have considered the Boffetta-Osburne method for direct NFT, and the numerical solution of the GLM integral equation for INFT implementation. For both methods, we have simplified their complex mathematical derivation and have shown how to apply them step-by-step for an improved numerical implementation, as they are subject to various approximation errors that may drastically affect the final performance results.

We have identified four parameters that affect the performances of our numerical implementations of NFT and INFT during the validation process. They consist of two temporal parameters: T function width and M discrete points, and two spectral parameters: resolution in frequency and the number of points N. The temporal parameter, T, requires the function to be considered zero outside of $2T$, and this will be the source of additional overhead data that are used in a real modulated signal. The total discrete points, $2M+1$, will define how many data points we can transmit at the same time, similar to the OFDM scheme. A lower number will affect the performance of the system. The number of points N was selected to be a high number with the assumption of having a large oversampling factor at the transmitter side. All these parameters should be optimized in order to gain the maximum performance from the nonlinear Fourier transform. We also observed that for large amplitudes, the performance drastically decreased because when truncating the integrals, due to the Gibbs effect, the errors will also increase. In addition, during each iteration, the errors are accumulated by affecting the performance of the system. In Section 3.1, we showed the steps of the algorithm to be used to properly represent the signal in the time and frequency domains, and how we should select the domain of the Marchenko system for INFT implementation.

We then conducted extensive numerical simulations on three different modulation schemes. Our results suggest that b-modulation is an effective method to be used in optical fiber communication systems. It shows better performance with a higher Q-factor at a fixed input average power and distance.

The results presented in this paper show the feasibility of a continuous nonlinear optical-based spectrum with our proposed methodology for NFT/INFT implementation. Despite the fact that this is a significant step forward for numerical simulation, it needs further research targeting of least NFT/INFT complexity and higher accuracy algorithm. For future work, we plan to use only the Boffeta-Osburne method for both direct NFT and INFT using the duality principle, and we compare the results with a system where only the Marchenko system is used. Using only one algorithm without performance degradation will probably make the system less prone to numerical errors.