Information Transmission using the Nonlinear Fourier Transform, Part II:Numerical Methods 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, e-mails: {mansoor,frank}@comm.utoronto.ca.

# Information Transmission using the Nonlinear Fourier Transform, Part II: Numerical Methods

## Abstract

In this paper, numerical methods are suggested to compute the discrete and the continuous spectrum of a signal with respect to the Zakharov-Shabat 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 layer-peeling 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.

Optical fiber communication, forward nonlinear Fourier transform, Zakharov-Shabat spectral problem, numerical methods, operator eigenproblem.
2

## 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 frequency-division multiplexing (NFDM), a scheme that uses the nonlinear Fourier transform for data communication over integrable channels. NFDM extends traditional orthogonal frequency-division 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 input-output relations in the NFT domain are (see [Part I])

 ^Y(λ) = H(λ)^X(λ)+^Z, ~Y(λj) = H(λj)~X(λj)+~Z,

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 time-domain 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 real-time 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 Zakharov-Shabat 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 layer-peeling method, Ablowitz-Ladik 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 Zakharov-Shabat system, which is a Lax operator for the nonlinear Schrödinger equation.

For later use, we recall that the slowly-varying complex envelope of a narrow-band small-amplitude signal propagating in a dispersive weakly-nonlinear 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:

 jqz = qtt+2|q|2q. (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

 L=j⎛⎝∂∂t−q(t,z)−q∗(t,z)−∂∂t⎞⎠. (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

 vt=(−jλq(t)−q∗(t)jλ)v. (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 two-dimensional 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

 v1(T2,λ)=(01)ejλT2.

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

 v2(t,λ)=a(λ)~v1(t,λ)+b(λ)v1(t,λ). (4)

Following Assumption 1(b), a particular choice for is made by solving the system

 vt=(−jλq(t)−q∗(t)jλ)v,v(T1,λ)=(10)e−jλT1, (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

 a(λ) = v1(T2)ejλT2, (6a) b(λ) = v2(T2)e−jλT2. (7a)

The NFT of a signal consists of a continuous spectral function defined on the real axis

 ^q(λ)=b(λ)a(λ),λ∈R, (8)

and a discrete spectral function defined on the upper half complex plane

 ~q(λj) = b(λj)da(λ)/dλ|λ=λj,j=1,⋯,N,

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:

1. Methods which estimate the continuous spectrum by directly integrating the Zakharov-Shabat system; see Section III.

2. 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 Zakharov-Shabat 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 Zakharov-Shabat 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 .

### Iii-a Forward and Central Discretizations

The most obvious method to attempt to solve (5) is the first-order 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

 P[k]\lx@stackrelΔ=(−jλq[k]−q∗[k]jλ). (9)

Integrating both sides of (5) from to and assuming that the right hand side is constant over this interval, we get

 v[k+1] = Af[k]v[k],k=0,…,N, v[0] = (10)e−jλT1, (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 one-step Euler method does not produce satisfactory results for affordable small step sizes .

One can improve upon the basic Euler method by considering the central-difference iteration [3],

 v[k+1] = v[k−1]+2ϵP[k]v[k].

This makes the discretization second-order, 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).

### Iii-B Fourth-Order Runge-Kutta Method

One can also employ higher-order integration schemes such as the Runge-Kutta methods. Improved results are obtained using the fourth-order Runge-Kutta method [4, 5]. However it takes significant time to estimate the spectrum using such higher-order numerical methods in real-time. 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.

### Iii-C Layer-Peeling Method

In Section IV. C of [Part I], we have calculated the nonlinear spectra of a rectangular pulse. One can approximate as a piece-wise constant signal and use the layer-peeling 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 layer-peeling method read

 (a[k+1],b[k+1]) = (a[k],b[k])∘(x[k],y[k]), (11) (a[0],b[0]) = (1,0),

where the operation is defined as in [1]

 a[k+1]=a[k]x[k]−b[k]¯y[k], b[k+1]=a[k]y[k]+b[k]¯x[k],

in which

 x[k] = (cos(Dϵ)−jλDsin(Dϵ))ejλ(t[k]−t[k−1]), y[k] = −q∗[k]Dsin(Dϵ)e−jλ(t[k]+t[k−1]), (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 layer-peeling 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 (layer-peeling), 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 layer-peeling method gives remarkably accurate results in estimating the nonlinear Fourier transform.

### Iii-D Crank-Nicolson Method

In the Crank-Nicolson method, the derivative of the evolution parameter is approximated by a finite-difference approximation, e.g., forward discretization, and other functions are discretized by taking their average over the end points of the discretization interval:

 v[k+1]−v[k]ϵ=12(P[k]v[k]+P[k+1]v[k+1]),

where is defined in (9). This implicit iteration can be made explicit

 v[k+1]=Acnv[k],k=0,⋯,N,

with the initial condition (10), where

 Acn=(I−ϵ2P[k+1])−1(I+ϵ2P[k]).

As we will see, this simple scheme too gives good results in estimating the nonlinear spectrum.

Ablowitz-Ladik (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 Ablowitz-Ladik discretization of the NLS equation for solving the Zakharov-Shabat eigenproblem in the spectral domain.

Discretization sometimes breaks symmetries, making the discrete version of an integrable equation no longer integrable. A consequence of symmetry-breaking 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 Korteweg-de 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

 Aal[k]=(zQ[k]−Q∗[k]z−1). (13)

 v[k+1]=Aal[k]v[k],v[0]=(10)e−jλT1. (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

 L=(Z−Q[k]−Q∗[k−1]α[k−1]Z−1), (15)

where , and is the shift operator, i.e., , . To the first order in , and (15) can be simplified to

 L=(Z−Q[k]−Q∗[k]Z−1). (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]

 jdq[k]dz = q[k+1]−2q[k]+q[k−1]ϵ2 (17) +|q[k]|2(q[k+1]+q[k−1]).

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 Ablowitz-Ladik 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 continuous-time Zakharov-Shabat system in the spectral domain. This is a non-finite-difference discretization, capable of dealing with oscillations in the Zakharov-Shabat system, which greatly enhances the accuracy of the one-step finite-difference methods.

Following Tao and Thiele’s approach [6] and [4], we can also normalize the matrix

 v[k+1]=1√1+|Q[k]|2(zQ[k]−Q[k]∗z−1)v[k]. (18)

The scale factor does not change the spectrum significantly, since it is canceled out in the ratios and , and also its effects are second-order in . However, numerically, normalization may help in reducing the numerical error. In subsequent sections, we refer to (14) as the Ablowitz-Ladik method (AL1) and to (18) as the modified Ablowitz-Ladik 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 two-dimensional 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.

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

2. 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.

### Iv-a 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)

 a(λ) = v1[N]ejλt[N].

Taking the derivative with respect to , we obtain

 da(λ)dλ = ((v1[N])′+jt[N]v1[N])ejλt[N].

We can update the derivative information along with . In methods of Section III, the transformation of eigenvectors from to can be generally represented as

 v[k+1]=A[k]v[k],

for some suitable one-step update matrix , which varies from method to method. Differentiating with respect to and augmenting with , we get the iterations

 v[k+1] = A[k]v[k], (19a) v′[k+1] = A′[k]v[k]+A[k]v′[k], (20a)

with initial conditions

 v[0]=(10)e−jλt[0],v′[0]=(−jt[0]0)e−jλt[0].

The derivative matrix depends on the method used.

For the forward discretization scheme:

 A′f=M1=(−j00j)ϵ.

For the Crank-Nicolson scheme:

 A′cn=ϵ2(I−ϵ2P[k+1])−1Σ3(I+Acn[k]),

where .

 A′al=(−jz00jz−1)ϵ.

The desired coefficients are obtained at as follows

 a(λ) = v1[N]ejλt[N], b(λ) = v2[N]e−jλt[N], a′(λ) = (v′1[N]+jt[N]v1[N])ejλt[N].

Similarly, the layer-peeling iteration can be augmented to update as well:

 a′[k+1] = a′[k]x[k]+a[k]x′[k]−(b′[k]¯y[k]+b[k]¯y′[k]), b′[k+1] = a′[k]y[k]+a[k]y′[k]+b′[k]¯x[k]+b[k]¯x′[k], a′[0] = b′[0]=0,

where

 x′[k] = jϵ(1−λ2D2)(cos(Dϵ)−sin(Dϵ)Dϵ)ejλϵ, y′[k] = −q∗[k]{λϵD2cos(Dϵ) −(λD3+jt[k]+t[k−1]D)sin(Dϵ)}e−jλ(t[k]+t[k−1]).

The expressions for and are similar, with all ’s replaced with and replaced with .

With the derivative information being available, the Newton-Raphson method is a good scheme to search for the location of the (discrete) eigenvalues. The iteration for the complex-valued Newton-Raphson scheme is

 λk+1=λk−αka(λk)a′(λk), (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 time-domain conserved quantity as the sum of discrete and continuous spectral terms:

 E(k)=E(k)disc+E(k)cont,k=0,1,⋯,

where

 E(k)disc = 4kN∑i=1I(λki), E(k)cont =

and where are time-domain conserved quantities (functionals of the signal). The first few conserved quantities are energy

 E(1)=∞∫−∞|q(t,z)|2dt,

momentum

 E(2)=12j∞∫−∞q(t)dq∗(t)dtdt,

and Hamiltonian

 E(3)=1(2j)2∞∫−∞(|q(t)|4−∣∣∣dq(t)dt∣∣∣2)dt.

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

 Eerror=(E(1)−E(1)cont,E(2)−E(2)cont,E(3)−E(3)cont), (22)

gives an estimate on the number of remaining eigenvalues. When a new eigenvalue is found, this error is updated as

 Eerror:=Eerror−(4Iλ,2Iλ2,43Iλ3).

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.

### Iv-B Discrete Spectrum as a Matrix Eigenvalue Problem

The methods mentioned in Section IV-A 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 small-sized 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 search-based methods. The matrix-based 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.

#### Central-Difference Eigenproblem

The matrix eigenvalue problem can be formulated in the time domain or the frequency domain [5]. Consider the Zakharov-Shabat system in the form

 j⎛⎝∂∂t−q−q∗−∂∂t⎞⎠v=λv. (23)

In the time domain, one can replace the time derivative by the central finite-difference matrix

 D=12ϵ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝010⋯−1−101⋯0⋯00−101100−10⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

and expand (23) as

 j(D−diag(q[k])−diag(q∗[k])−D)v[k]=λv[k]. (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 non-Hermitian matrices. As a result, the eigenvalues of a non-Hermitian matrix (corresponding to a non-self-adjoint operator) are markedly difficult to calculate [12, 13]. Running a general-purpose eigenvalue calculation routine on (24) may therefore not be the most efficient way to get eigenvalues (). One can study the general-purpose numerical eigenvalue algorithms in [12] and tailor them to the discrete Zakharov-Shabat spectral problems in this paper. Note that matrix is (almost) tridiagonal. Thus the eigenproblem for (24), or equivalently discretization of (30) gives a second-order recursive equation . One can then obtain polynomial recursively in operations and proceed to find its roots. Alternatively, can be calculated recursively.

We can also rewrite the Ablowitz-Ladik discretization as a matrix eigenvalue problem. Using the operator (15), we obtain

 v1[k+1]−Q[k]v2[k] = zv1[k], (25a) −Q∗[k−1]v1[k]+α[k−1]v2[k−1] = zv2[k], (26a)

which consequently takes the form

 (U1−diag(Q[k])−diag(Q∗[k−1])UT2)v=zv,

in which

 U1 = ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝010⋯0001⋯0⋯0000110000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, U2 = ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0α[0]0⋯000α[1]⋯0⋯0000α[N−1]α[N]0000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

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 Ablowitz-Ladik 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.

#### Spectral Method

In the frequency domain, one can approximate derivatives with the help of the Fourier transform. Let us assume that

 v(t)=M2∑k=−M2(αkβk)ej2kπtT,q(t)=M2∑k=−M2γkej2kπtT,

where and is the maximum number of frequencies in all variables. Then the Zakharov-Shabat system is

 −αk2πkT−jM2∑m=−M2γk−mβm = λαk, −jM2∑m=−M2γ∗−k+mαm+βk2πkT = λβk.

Thus we obtain

 (ΩΓ−ΓH−Ω)(αβ)=λ(αβ),

where , , and is a Toeplitz matrix with the first row

 −j(γ0γ−1⋯γ−M20⋯0M2),

and the first column

 −j(γ0γ1⋯γM20⋯0M2)T.

The point spectrum is thus found by looking at the eigenvalues of the matrix

 A=(ΩΓ−ΓH−Ω).

## V Running Time, Convergence and Stability of the Numerical Methods

The numerical methods discussed in this paper are first-order 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

 u1 = v1ejλt, u2 = v2e−jλt,

and transform (5) to

 ut=(0q(t)e2jλt−q∗(t)e−2jλt0)u,u(T1,λ)=(10). (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:

 (a[k+1]b[k+1]) = (1Q[k]z−2k−Q∗[k]z2k1)(a[k]b[k]), (28) (a[0]b[0]) = (10).

The nonlinear Fourier coefficients are obtained as and . The discrete nonlinear Fourier transform mentioned in [6] is thus the forward discretization of the normalized Zakharov-Shabat 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 Zakharov-Shabat 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 first-order iteration . The eigenvalues of the matrix in this method are

 s1,2=1±jϵ|q[k]|.

It follows that the forward discretization of (27) gives rise to eigenvalues outside of the unit disk, . As a result, first-order 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 right-hand 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 finite-difference discretizations has been observed in [14] for the sine-Gordon 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 Ablowitz-Ladik discretization of the NLS equation (in time) has the desirable property that chaos and numerical instability, which are sometimes present in finite-difference discretizations of the NLS equation, do not appear at all. Though these results are for the original time-domain 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 Ablowitz-Ladik discretization of the Zakharov-Shabat 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 Zakharov-Shabat 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]

 ⎧⎨⎩<