Information Transmission using the Nonlinear Fourier Transform,
Part II:
Numerical Methods
^{†}^{†}thanks: Submitted for publication on April 3,
2012; accepted December 14, 2012. The authors are with the Edward S. Rogers Sr. Dept. of
Electrical and Computer Engineering, University of Toronto, Toronto,
ON M5S 3G4, Canada, emails: {mansoor,frank}@comm.utoronto.ca.
Abstract
In this paper, numerical methods are suggested to compute the discrete and the continuous spectrum of a signal with respect to the ZakharovShabat system, a Lax operator underlying numerous integrable communication channels including the nonlinear Schrödinger channel, modeling pulse propagation in optical fibers. These methods are subsequently tested and their ability to estimate the spectrum are compared against each other. These methods are used to compute the spectrum of various signals commonly used in the optical fiber communications. It is found that the layerpeeling and the spectral methods are suitable schemes to estimate the nonlinear spectra with good accuracy. To illustrate the structure of the spectrum, the locus of the eigenvalues is determined under amplitude and phase modulation in a number of examples. It is observed that in some cases, as signal parameters vary, eigenvalues collide and change their course of motion. The real axis is typically the place from which new eigenvalues originate or are absorbed into after traveling a trajectory in the complex plane.
I Introduction
The nonlinear Fourier transform (NFT) of a signal is a pair of functions: the continuous spectrum , , and the discrete spectrum , , . The NFT arises in the study of integrable waveform channels as defined in Part I [1]. In such channels, signals propagate (in a potentially complicated manner) according to a given integrable evolution equation, whereas the nonlinear Fourier transform of the signal propagates according to a (simple) multiplication operator.
In [Part I], we proposed nonlinear frequencydivision multiplexing (NFDM), a scheme that uses the nonlinear Fourier transform for data communication over integrable channels. NFDM extends traditional orthogonal frequencydivision multiplexing (OFDM) to channels generatable by a Lax pair. An example is the optical fiber channel, where signal propagation is modeled by the (integrable) nonlinear Schrödinger (NLS) equation. In general, the channel inputoutput relations in the NFT domain are (see [Part I])
where and are continuous and discrete spectra at the input of the channel, and are spectra at the output of the channel, and and represent noise. The channel filter for the NLS equation is given by .
NFDM is able to deal directly with nonlinearity and dispersion, without the need for additional compensation at the transmitter or receiver. In this scheme, information is encoded in the nonlinear spectrum at the channel input, and the corresponding timedomain signal is transmitted. At the receiver, the NFT of the received signal is computed, and the resulting spectra and are subsequently used to recover the transmitted information.
Similar to the ordinary Fourier transform, while the NFT can be computed analytically in a few cases, in general, numerical methods are required. Such methods must be robust, reliable and fast enough to be implemented in realtime at the receiver. In this paper, we suggest and evaluate the performance of a number of numerical algorithms for computing the forward NFT of a given signal. Using these algorithms, we then perform extensive numerical simulations to understand the behavior of the nonlinear spectrum for various pulse shapes and parameters commonly used in data communications.
We are aware of no published work presenting the NFT of various signals numerically, for many pulse shapes and parameters. Such work is necessary to clarify the structure of the nonlinear spectrum and help in its understanding. In part, this has been due to the fact that the NFT has largely remained a theoretical artifice, and practical implementation of the NFT as an applied tool has not yet been pursued in engineering.
We review the relevant literature in Section II and suggest new schemes for the numerical evaluation of the forward NFT. Although these methods are general and work for the AKNS system [2], for the purpose of illustration, we specialize the AKNS system to the ZakharovShabat system. All these methods are put to test in cases where analytical formulae exist and are compared with one another in Section VI. Only some of these methods will be chosen for the subsequent numerical simulations; these are the layerpeeling method, AblowitzLadik integrable discretization, and the spectral matrix eigenvalue scheme. These methods are used in the next sections to numerically compute the nonlinear Fourier transform of a variety of practical pulse shapes encountered in the data communications.
Ii The Nonlinear Fourier Transform
Details of the nonlinear Fourier transform can be found in [Part I]. Here we briefly recall a few essential ingredients required in the numerical computation of the forward transform. As noted earlier, we illustrate numerical methods in the context of the ZakharovShabat system, which is a Lax operator for the nonlinear Schrödinger equation.
For later use, we recall that the slowlyvarying complex envelope of a narrowband smallamplitude signal propagating in a dispersive weaklynonlinear medium, such as an optical fiber, satisfies the cubic nonlinear Schrödinger equation. By proper scaling, the equation can be normalized to the following dimensionless form in dimensions:
(1) 
Here denotes retarded time, and is distance.
The NFT for an integrable evolution equation starts by finding a Lax pair of operators and such that the evolution equation arises as the compatibility condition . For the NLS equation, we may take operator as
(2) 
(The corresponding operator can be found in [Part I].)
The NFT is defined via the spectral analysis of the operator, given in this paper by (2). The spectrum of is found by solving the eigenproblem , where is an eigenvalue of and is its associated eigenvector. It can be shown that the operator in (2) has the isospectral flow property, i.e., its spectrum is invariant even as evolves according to the NLS equation.
The eigenproblem can be simplified to
(3) 
Note that the dependence of is suppressed in (3) (and throughout this paper), as this variable comes into play only in the propagation of the signal, not in the definition and computation of the NFT.
Assumption 1.
Throughout this paper we assume that (a) , and (b) is supported in the finite interval . ∎
The set of eigenvectors associated with eigenvalue in (3) is a twodimensional subspace of the continuously differentiable functions. We define the adjoint of a vector as . If is an element of , then is an element of . It can be shown that any pair of eigenvectors and forms a basis for [Part I]. Using Assumption 1(b), we can select an eigenvector to be a solution of (3) with the boundary condition
The basis eigenvectors and are called canonical eigenvectors.
Having identified a basis for the subspace , we can project any other eigenvector on this basis according to
(4) 
Following Assumption 1(b), a particular choice for is made by solving the system
(5) 
in which we dropped the superscript in for convenience. By solving (5) in the interval for a given and obtaining , the nonlinear Fourier coefficients and can be obtained by considering (4) at . The resulting coefficients obtained in this manner are
(6a)  
(7a) 
The NFT of a signal consists of a continuous spectral function defined on the real axis
(8) 
and a discrete spectral function defined on the upper half complex plane
where are eigenvalues and correspond to the (isolated) zeros of in , i.e., .
From the discussions made, in order to compute the nonlinear spectrum of , the system of differential equations (5) needs to be solved in the interval . Except for special cases, (5) needs to be solved numerically.
Numerical methods for the calculation of the forward nonlinear Fourier transform are divided into two classes in this article:

Methods which estimate the continuous spectrum by directly integrating the ZakharovShabat system; see Section III.

Methods which find the (discrete) eigenvalues. Two approaches are suggested in this paper for this purpose. Similar to the continuous spectrum estimation, we can integrate the ZakharovShabat system numerically and obtain . To find zeros of , the scheme is often supplemented with a search method to locate eigenvalues in the upper half complex plane. One can also discretize and rewrite the ZakharovShabat eigenproblem in the interval as a (large) matrix eigenvalue problem; see Section IV.
We begin by discussing methods which estimate the continuous spectrum.
Iii Numerical Methods for Computing the Continuous Spectrum
In this section, we assume that is given and provide algorithms for calculating the nonlinear Fourier coefficients and . The continuous spectral function is then easily computed from (8). This process can be repeated to compute the spectral amplitudes for any desired finite set of continuous frequencies .
Iiia Forward and Central Discretizations
The most obvious method to attempt to solve (5) is the firstorder Euler method or one of its variations [3].
Recall that the signal is supported in the finite time interval , and partition this interval uniformly according to the mesh with size , i.e., with . Let and let
(9) 
Integrating both sides of (5) from to and assuming that the right hand side is constant over this interval, we get
(10) 
where and is the identity matrix. Equation (10) is iterated from to to find . The resulting vector is subsequently substituted in (6a)–(7a) to obtain and .
We have implemented the Euler method for the calculation of the nonlinear Fourier transform of a number of pulse shapes. Unfortunately, the onestep Euler method does not produce satisfactory results for affordable small step sizes .
One can improve upon the basic Euler method by considering the centraldifference iteration [3],
This makes the discretization secondorder, i.e., the error is of order . Here an additional initial condition is required too, which can be obtained, e.g., by performing one step of the regular forward difference (10).
IiiB FourthOrder RungeKutta Method
One can also employ higherorder integration schemes such as the RungeKutta methods. Improved results are obtained using the fourthorder RungeKutta method [4, 5]. However it takes significant time to estimate the spectrum using such higherorder numerical methods in realtime. Since the method, with its typical parameters, is quite slow and does not outperform some of the schemes suggested in the following sections, we do not elaborate on this method here; see [3] for details. However, for comparison purposes we will include this scheme in our numerical simulations given in Section VI.
IiiC LayerPeeling Method
In Section IV. C of [Part I], we have calculated the nonlinear spectra of a rectangular pulse. One can approximate as a piecewise constant signal and use the layerpeeling property of the nonlinear Fourier transform to estimate the spectrum of any given signal. Let and be the nonlinear Fourier coefficients of in the interval , and and coefficients in the small (rectangular) region . The iterations of the layerpeeling method read
(11)  
where the operation is defined as in [1]
in which
(12) 
and , , . The desired coefficients are obtained as and . Note that the exponential factors in and enter (12) in a telescopic manner. As a result, for the numerical implementation, it is faster to drop these factors and just scale the resulting and coefficients by and , respectively. This, however, reduces the accuracy as it involves the product of large and small numbers.
We are motivated by [6] in which the layerpeeling identity (11) is mentioned as a property of the nonlinear Fourier transform. An equivalent presentation of this method is given in [7] as well.
Note further that a different numerical method, but with the same name (layerpeeling), exists in geophysics and fiber Bragg design [8]; however this method is not directly related to the forward NFT problem considered here.
We shall see in Section VI that the layerpeeling method gives remarkably accurate results in estimating the nonlinear Fourier transform.
IiiD CrankNicolson Method
In the CrankNicolson method, the derivative of the evolution parameter is approximated by a finitedifference approximation, e.g., forward discretization, and other functions are discretized by taking their average over the end points of the discretization interval:
where is defined in (9). This implicit iteration can be made explicit
with the initial condition (10), where
As we will see, this simple scheme too gives good results in estimating the nonlinear spectrum.
IiiE The AblowitzLadik Discretization
AblowitzLadik (AL) discretization is an integrable discretization of the NLS equation in time domain [9]. In this section, we suggest using the Lax pairs of the AblowitzLadik discretization of the NLS equation for solving the ZakharovShabat eigenproblem in the spectral domain.
Discretization sometimes breaks symmetries, making the discrete version of an integrable equation no longer integrable. A consequence of symmetrybreaking is that quantities which are conserved in the continuous model may no longer be invariant in the discrete model. A completely integrable Hamiltonian system with an infinite number of conserved quantities might have a discretized version with no, or few, conserved quantities. The discrete equation therefore does not quite mimic the essential features of the original equation if the step size is not small enough.
However, for some integrable equations, discretizations exist which are themselves completely integrable Hamiltonian systems, i.e., they possess an infinite number of conserved quantities and are linearizable by a Lax pair, and therefore are solvable by the nonlinear Fourier transform. Such developments exist for the NLS and Kortewegde Vries (KdV) equations.
For the NLS equation, the integrable discrete version was introduced by Ablowitz and Ladik [9]. To illustrate the general idea, let us replace for small with in the forward discretization method (10) (the opposite of what is usually done in practice). Let represent the discrete eigenvalue, , and
(13) 
The AblowitzLadik iteration is
(14) 
Under the transformation, the upper half complex plane in domain is mapped to the exterior of the unit circle in the domain. The continuous spectrum therefore lies on the unit circle , while the discrete spectrum lies outside of the unit circle .
One can rewrite the equation (13) in the eigenvalue form , with the following operator
(15) 
where , and is the shift operator, i.e., , . To the first order in , and (15) can be simplified to
(16) 
Given the operator (16), one can consider the operator of the continuous NLS equation and modify its elements such that the compatibility equation represents a discretized version of the NLS equation. It is not hard to verify that after doing so we are led to an operator resulting in the following discrete integrable NLS equation [9, 10]
(17)  
Here the space derivative remains intact and the signal is discretized in time, in such a way that the nonlinearity is somewhat averaged among three time samples. In the continuum limit , (17) approaches the continuous NLS equation and its merits lie in the fact that it is integrable for any , not just in the limit . For example, soliton signals can be observed in this model for any . The equation has its own infinite number constants of motion, approaching integrals of motion in the continuum limit. The operator which leads to (17), and the details of the nonlinear Fourier transform for (17) can be found in [9, 10].
We conclude that the AblowitzLadik discretization can be used not only as a means to discretize the NLS equation in the time domain [11], but also as a means to solve the continuoustime ZakharovShabat system in the spectral domain. This is a nonfinitedifference discretization, capable of dealing with oscillations in the ZakharovShabat system, which greatly enhances the accuracy of the onestep finitedifference methods.
Following Tao and Thiele’s approach [6] and [4], we can also normalize the matrix
(18) 
The scale factor does not change the spectrum significantly, since it is canceled out in the ratios and , and also its effects are secondorder in . However, numerically, normalization may help in reducing the numerical error. In subsequent sections, we refer to (14) as the AblowitzLadik method (AL1) and to (18) as the modified AblowitzLadik method (AL2).
Iv Methods for Calculating the Discrete Spectrum
In order to compute the discrete spectrum, the zeros of in the upper half complex plane must be found. One way to visualize this is to assume a twodimensional mesh in and determine at all mesh points. Discrete eigenvalues are then easily identified by looking at the graph of ; in many cases they correspond to deep and narrow “wells” corresponding to the zeros of the magnitude of .
As noted earlier, two types of methods are suggested to calculate the point spectrum.

One can use the integrationbased algorithms mentioned in Section III which calculate nonlinear Fourier coefficients, and search for eigenvalues using a root finding method, such as the NewtonRaphson method. Such methods require good initial points and one needs to be careful about convergence [3].

It is also possible to rewrite the spectral problem for an operator as a (large) matrix eigenvalue problem. The point spectrum of the operator can be found in this way too.
Iva Discrete Spectrum via Search Methods
To calculate the discrete spectral amplitude , we require as well. As we will show, information about the derivative of can be updated recursively along with the information about , without resorting to approximate numerical differentiation.
Recall that the nonlinear Fourier coefficient is given by (6a)
Taking the derivative with respect to , we obtain
We can update the derivative information along with . In methods of Section III, the transformation of eigenvectors from to can be generally represented as
for some suitable onestep update matrix , which varies from method to method. Differentiating with respect to and augmenting with , we get the iterations
(19a)  
(20a) 
with initial conditions
The derivative matrix depends on the method used.
For the forward discretization scheme:
For the CrankNicolson scheme:
where .
For the AblowitzLadik method:
The desired coefficients are obtained at as follows
Similarly, the layerpeeling iteration can be augmented to update as well:
where
The expressions for and are similar, with all ’s replaced with and replaced with .
With the derivative information being available, the NewtonRaphson method is a good scheme to search for the location of the (discrete) eigenvalues. The iteration for the complexvalued NewtonRaphson scheme is
(21) 
where is some step size modifier; usually . The iteration stops if is almost stationary, i.e., if for a small . In practice, the quadratic convergence of the scheme is often very fast and occurs in just a few iterations.
In data communications, since noise is usually small, the points in the transmitted constellation can (repeatedly) serve as the initial conditions for (21). In this case, convergence is usually achieved in a couple of iterations. For an unknown signal, random initial conditions are chosen. In either case, one or more sequence of Newton iterations have to be performed for any single eigenvalue.
To make sure that all of the eigenvalues are found, we can check the trace formula for [Part I]. The trace formula is a time frequency identity relating the hierarchy of infinitely many conserved quantities to the spectral components.
In general, the trace formula represents a timedomain conserved quantity as the sum of discrete and continuous spectral terms:
where
and where are timedomain conserved quantities (functionals of the signal). The first few conserved quantities are energy
momentum
and Hamiltonian
For , the trace formula is a kind of Parseval’s identity (Plancherel’s theorem), representing the total energy of the signal in time as the sum of the energy of the discrete and continuous spectral functions. When satisfied, the Parseval’s identity ensures that all of the signal energy has been accounted for.
We first calculate the continuous spectrum and its “energy terms” for a sufficiently fine mesh on the real axis. The norm of the error vector
(22) 
gives an estimate on the number of remaining eigenvalues. When a new eigenvalue is found, this error is updated as
The process is repeated until is less than a small prescribed tolerance value.
To summarize, given the signal , its nonlinear Fourier transform can be computed based on Algorithm 1.
IvB Discrete Spectrum as a Matrix Eigenvalue Problem
The methods mentioned in Section IVA find the discrete spectrum by searching for eigenvalues in the upper half complex plane. Sometimes it is desirable to have all eigenvalues at once, which can be done by solving a matrix eigenvalue problem [5]. These schemes obviously estimate only (discrete) eigenvalues and do not give information on the rest of the spectrum. Since the matrix eigenvalue problem can be solved quickly for smallsized problems, it might take less computational effort to compute the discrete spectrum in this way. In addition, one does not already need the continuous spectrum to estimate the size of the discrete spectrum. On the other hand, for large matrices that arise when a large number of signal samples are used, the matrix eigenproblem (which is not Hermitian) is slow and it may be better to find the discrete spectrum using the searchbased methods. The matrixbased methods also have the disadvantages that they can generate a large number of spurious eigenvalues, and one may not be able to restrict the algorithm for finding eigenvalues of a matrix to a certain region of the complex plane.
Below we rewrite some of the methods mentioned in the Section III as a regular matrix eigenvalue problem, for the computation of the discrete spectrum.
IvB1 CentralDifference Eigenproblem
The matrix eigenvalue problem can be formulated in the time domain or the frequency domain [5]. Consider the ZakharovShabat system in the form
(23) 
In the time domain, one can replace the time derivative by the central finitedifference matrix
and expand (23) as
(24) 
The point spectrum is contained in the eigenvalues of the matrix in the left hand side of (24).
Eigenvalues of a real symmetric or Hermitian matrix can be found relatively efficiently, owing to the existence of a complete orthonormal basis and the stability of the eigenvalues. In this case, a unitary matrix can be designed such that the sequence of similarity transformations , , converges to an almost diagonal (or a triangular) matrix. The similarity transformation can be found using, e.g., the QR factorization or the Householder transformation [12].
Unfortunately, most of the useful statements about computations using Hermitian matrices cannot usually be generalized to nonHermitian matrices. As a result, the eigenvalues of a nonHermitian matrix (corresponding to a nonselfadjoint operator) are markedly difficult to calculate [12, 13]. Running a generalpurpose eigenvalue calculation routine on (24) may therefore not be the most efficient way to get eigenvalues (). One can study the generalpurpose numerical eigenvalue algorithms in [12] and tailor them to the discrete ZakharovShabat spectral problems in this paper. Note that matrix is (almost) tridiagonal. Thus the eigenproblem for (24), or equivalently discretization of (30) gives a secondorder recursive equation . One can then obtain polynomial recursively in operations and proceed to find its roots. Alternatively, can be calculated recursively.
IvB2 AblowitzLadik Eigenproblem
We can also rewrite the AblowitzLadik discretization as a matrix eigenvalue problem. Using the operator (15), we obtain
(25a)  
(26a) 
which consequently takes the form
in which
and , . Note that all shifting operations are cyclic, so that all vector indices remain in the interval .
Similar results can be obtained for the simplified operator (16) instead of (15). This corresponds to the above discretization with and replaced with . Similarly, one can rewrite the normalized AblowitzLadik iteration (18) as a matrix eigenvalue problem. This corresponds to the operator (15) with the top left and bottom right entries replaced, respectively, with and . Accordingly and in (25a) and (26a) are replaced, respectively, with and and in the given expressions.
IvB3 Spectral Method
In the frequency domain, one can approximate derivatives with the help of the Fourier transform. Let us assume that
where and is the maximum number of frequencies in all variables. Then the ZakharovShabat system is
Thus we obtain
where , , and is a Toeplitz matrix with the first row
and the first column
The point spectrum is thus found by looking at the eigenvalues of the matrix
V Running Time, Convergence and Stability of the Numerical Methods
The numerical methods discussed in this paper are firstorder matrix iterations and therefore the running time of all of them to get is multiplications and additions per eigenvalue. This corresponds to a complexity of operations for the calculation of the continuous spectrum on a mesh with eigenvalues. The exact number of operations depends on the details of the implementation and the memory requirement of the method. All iterative methods thus take about the same time asymptotically, albeit with different coefficients.
An important observation is that, while the Fast Fourier Transform (FFT) takes operations to calculate the spectral amplitudes of a vector with length at equispaced frequencies, the complexity of the methods described in this paper to compute the continuous spectrum are . Similarly, it takes operations to calculate the discrete spectrum. In other words, so far we do not exploit the potentially repetitive operations in our computations.
It is evident from (6a) that as , should grow as so that is a finite complex number. The canonical eigenvector thus has an unbounded component as (i.e., ). One can, however, normalize and according to
and transform (5) to
(27) 
The desired coefficients are simply and . Consequently, if one is interested in obtaining eigenvectors in addition to the coefficients and , the discretization of the normalized system (27) has better numerical properties:
(28)  
The nonlinear Fourier coefficients are obtained as and . The discrete nonlinear Fourier transform mentioned in [6] is thus the forward discretization of the normalized ZakharovShabat system (27).
We are interested in the convergence of (or and ) as a function of for fixed values of and . That is to say, we require that the error as (for fixed and ). The (global) error in all methods described in this paper is at least , and therefore all these methods are convergent.
Some of these methods are, however, not stable. This is partially because the ZakharovShabat system can have unbounded solutions for , i.e., as . Errors can potentially be amplified by the system unstable dynamics. One should be cautious about the normalized system (27) as well. For example, forward discretization of the normalized system (28) gives a firstorder iteration . The eigenvalues of the matrix in this method are
It follows that the forward discretization of (27) gives rise to eigenvalues outside of the unit disk, . As a result, firstorder discretization of (27) are also unstable. In cases where , we can consider normalizing the iterations by dividing by (in the case of (28), dividing the righthand side by ). The resulting iteration has eigenvalues inside the unit disk. For the effect is only second order in , however it helps in managing the numerical error if larger values of are chosen.
An issue pertinent to numerical methods is chaos. Chaos and numerical instability of finitedifference discretizations has been observed in [14] for the sineGordon equation, which is also integrable and shares a number of basic properties with the NLS equation. In [15], the authors conclude that the standard discretizations of the cubic nonlinear Schrödinger equation may lead to spurious numerical behavior. This instability is related to the homoclinic orbits of the NLS equation, i.e., it occurs if the initial signal is chosen to be close to the homoclinic orbit of the equation. It disappears only if the step size is made sufficiently small, which can be smaller than what is desired in practice.
It is shown in [15] that the AblowitzLadik discretization of the NLS equation (in time) has the desirable property that chaos and numerical instability, which are sometimes present in finitedifference discretizations of the NLS equation, do not appear at all. Though these results are for the original timedomain equations, the issue can occur in the spectral eigenvalue problem (5) as well, if the signal is close to a certain family of functions (related to and ). Therefore, among discretizations studied, the AblowitzLadik discretization of the ZakharovShabat system is immune to chaos and numerical instability. This is particularly important in the presence of amplifier noise, where chaos can be more problematic.
Vi Testing and Comparing the Numerical Methods
In this section, we test and compare the ability of the suggested numerical schemes to estimate the nonlinear Fourier transform (with respect to the ZakharovShabat system) of various signals. Numerical results are compared against analytical formulae, in a few cases where such expressions exist. Our aim is to compare the speed and the precision of these schemes for various pulse shapes in order to determine which ones are best suited for subsequent simulation studies.
To derive the analytical formulae, recall that the continuous spectral function can be written as , in which satisfies [1]
(29) 
Similarly, one can solve the secondorder differential equation