# Fast Nonlinear Fourier Transform using Chebyshev Polynomials

We explore the class of exponential integrators known as exponential time differencing (ETD) method in this letter to design low complexity nonlinear Fourier transform (NFT) algorithms that compute discrete approximations of the scattering coefficients in terms of the Chebyshev polynomials for real values of the spectral parameter. In particular, we discuss ETD Runge-Kutta methods which yield algorithms with complexity O(N^2N) (where N is the number of samples of the signal) and an order of convergence that matches the underlying one-step method.

03/26/2019

### Efficient Nonlinear Fourier Transform Algorithms of Order Four on Equispaced Grid

We explore two classes of exponential integrators in this letter to desi...
09/09/2019

### A Fast Chebyshev Spectral Method for Nonlinear Fourier Transform

In this letter, we present a fast and well-conditioned spectral method b...
12/03/2018

### Fast Nonlinear Fourier Transform Algorithms Using Higher Order Exponential Integrators

The nonlinear Fourier transform (NFT) has recently gained significant at...
09/16/2019

### A Chebyshev Spectral Method for Nonlinear Fourier Transform: Norming Constants

In this paper, we present a Chebyshev based spectral method for the comp...
02/01/2019

### Multi-layered Cepstrum for Instantaneous Frequency Estimation

We propose the multi-layered cepstrum (MLC) method to estimate multiple ...
12/14/2020

### Is FFT Fast Enough for Beyond-5G Communications?

In this work, we consider the complexity and throughput limits of the Fa...
04/01/2020

### From Fourier to Koopman: Spectral Methods for Long-term Time Series Prediction

We propose spectral methods for long-term forecasting of temporal signal...

## I Introduction

This letter addresses the signal processing aspect of a nonlinear Fourier transform (NFT) based optical fiber communication system which has emerged as one of the potential ways of mitigating nonlinear signal distortions at higher signal powers [1, 2]. The idea consists in encoding information in the nonlinear Fourier spectrum of the signal which is synthesized using the inverse NFT. The information can be decoded using the direct NFT which forms the primary focus of this letter. Let us note that inverse NFT algorithm of extremely high order of convergence for computing the radiative potential have been proposed [3], which coupled with a higher order convergent direct NFT via the Darboux transformation procedure, could yield a full inverse NFT with higher order of convergence.

Recently, integrating factor (IF) based exponential integrators have been successfully exploited for the purpose of developing fast direct as well as inverse NFT algorithms [4, 5, 6, 7]. In this letter, we introduce another class exponential integrators which are known as exponential time differencing (ETD) method, a name originally used in computational electrodynamics, for the purpose of designing fast NFT algorithms. This particular class of algorithms rely on FFT based fast polynomial arithmetic in the Chebyshev basis as opposed to the monomial basis considered in the aforementioned works. While ETD schemes do not have any of the drawbacks associated with the IF schemes [8], the ETD schemes offer an additional advantage that there are no restrictions on the kind of the grid employed for sampling the potential. All ETD schemes within the Runge-Kutta as well as linear multistep methods can be sped up. In this letter, however, we have chosen the Runge-Kutta methods for our exposition.

We begin our discussion with a brief review of the scattering theory closely following the formalism presented in [4]. The NFT of any signal is defined via the Zakharov-Shabat (ZS) scattering problem (henceforth referred to as the ZS problem) which can be stated as follows: For and , where . The potential is defined by and with (). Here, is known as the spectral parameter and is the complex-valued signal. The solution of the ZS problem consists in finding the so called scattering coefficients which are defined through special solutions known as the Jost solutions which are linearly independent solutions of the ZS problem such that they have a plane-wave like behavior at or . The Jost solution of the second kind, denoted by , has the asymptotic behavior as . The asymptotic behavior as determines the scattering coefficients and for . In this letter, we primarily focus on the continuous spectrum, also referred to as the reflection coefficient, which is defined by for .

## Ii The Numerical Scheme

The exponential time differencing Runge-Kutta (ETD–RK) methods with internal stages can be characterized much in the same way as the standard Runge-Kutta method adapted to the ZS problem: Given the step-size , an ordered set of nodes defining the abscissas and the matrix-valued weight functions , , the internal stages (for ) of the ETD–RK method together with the final update read as

 (1)

where we have used the convention . The linear system stated above can be written as

 (I2s×2s−A(ζ;h)Dn02s×2−B(ζ;h)⊺Dnσ0)(Υnvn+1)=(E(ζ;h)vne−iζhσ3vn), (2)

where , , and, , and are defined as

 (3)

The existence of a solution of the linear system (2) depends on the determinant which can be shown to be nonzero for sufficiently small . For explicit RK methods, . Introducing the transfer matrix , the numerical scheme can be stated as

 vn+1=e−iζh[Δn+1]−1Mn+1(ζ;h)vn. (4)

Let us remark that collocation method [9] with Legendre–Gauss–Lobatto nodes is one of the simplest methods of constructing implicit ETD–RK methods. The weight functions in this case can be stated as follows: For , we have and

Evidently, where

 (n−1)!φn(z)=∫10e(1−τ)zτn−1dτ. (6)

The -functions defined above are central to all ETD schemes; therefore, we discuss it in more detail in Sec. II-B. In the following, we use the convention , and in order to represent the ETD–RK methods using the Butcher Tableau. An example of an implicit three-stage ETD–RK method is the Lobatto IIIA method of order , labelled as ETD–IRK, which is given by

 00001212φ3,2−34φ2,2+12φ1,2φ2,2−φ3,212φ3,2−14φ2,214φ3−3φ2+φ14φ2−8φ34φ3−φ24φ3−3φ2+φ14φ2−8φ34φ3−φ2

Noting that , the transfer matrix relation can be written as where

with . Example of an explicit four-stage ETD–RK method of order is given by [8]

 000001212φ1,200012012φ1,200112φ1,2(φ0,2−σ0)0φ1,204φ3−3φ2+φ12φ2−4φ32φ2−4φ34φ3−φ2

which we label as ETD–ERK.

### Ii-a Fast algorithm for Chebyshev nodes

Let the computational domain be and set . Let the number of steps be and . The grid is defined by with where is the number of samples. Consider the Jost solution : and . Let the cumulative products be denoted by and so that the discrete approximation to the scattering coefficients, and , can be worked out to be and .

Now, let us consider the problem of computing the samples of the continuous spectrum at the Chebyshev-Gauss-Lobatto (CGL) nodes given by for . Let us first note that the matrix-valued weight functions can be computed beforehand. The naive approach to forming the cumulative products and by direct multiplication evidently yields a complexity of . The fast algorithm, on the other hand, is based on the idea that each of the transfer matrices and the determinants can be represented using a small number of Chebyshev polynomials. The motivation behind this step is to use fast polynomial arithmetic in Chebyshev basis for computing the discrete scattering coefficients. To this end, define

 Mn→n+m−1=Mn+m−1×Mn+m−2×…×Mn, (7)

where is a power of . The strategy is to evaluate on CGL nodes by direct multiplication of the constituting matrices followed by a discrete Chebyshev transform to obtain the Chebyshev expansion coefficients:

 Mn→n+m−1(ξ;h)=d∑j=0M(j)n→n+m−1Tj(2πξh), (8)

for where denotes the Chebyshev polynomial of degree . Next, we also extend the strategy to the determinants . If the total number of matrices (or determinants) to be multiplied in the Chebyshev basis was originally, the aforementioned step reduces this number to at an initial cost of where is the cost of discrete Chebyshev transform of size .

The next step is to compute the products of the matrices (or determinants) pairwise as in [6] except we use FCT instead of FFT. Let , where is a power of , be the complexity of multiplying two polynomials of degree in the Chebyshev basis. Let be a power of and denote the complexity of multiplying matrices, then

. The overall complexity now depends on how the polynomial products are carried out in the Chebyshev basis which is considered next. Let us note that the Chebyshev transform of a complex vector of size

is equivalent to the DFT of a vector of size ; therefore, the complexity of the fast Chebyshev algorithm is  [10] where denotes the complexity of FFT of size .

In the following, we describe a fast algorithm for multiplying polynomials in the Chebyshev basis which proceeds by defining an equivalent multiplication problem in the monomial basis [10]: Let us consider the polynomials and of degree stated in the Chebyshev basis: and . Let be expressed as . Define , and . Also define and , then it can be shown by direct multiplication and using the property , for all that

 2rj=⎧⎪⎨⎪⎩γ0+2δD,j=0,γj+δD−j+δD+j,j=1,…,D−1,γj,j=D,…,2D. (9)

Let us show that the algorithm discussed above can be implemented using FFTs as opposed to FFTs of size

leading to a complexity estimate of

. The computation of is straightforward yielding a complexity of . The computation of can be carried out by employing no more than FFT of size as follows: Define whose degree is . Putting and , we have which can be easily computed using the quantities involved in computing . Finally, the inverse FFT of the sequence followed by a shift is required in order to determine which costs multiplications yielding an overall complexity of where we have retained only the leading term.

Noting that , the asymptotic complexity of computing the scattering coefficients using the polynomial multiplication algorithm described above works out to be . The original problem of evaluating the continuous spectrum at CGL nodes now reduces to carrying out an inverse discrete Chebyshev transform.

### Ii-B The exponential functions in ETD

In ETD–RK methods, the weights turn out to be linear combinations of the so-called –functions defined by (6). The first few of these functions are , and . These functions have a removable singularity at and it can be easily shown that . The higher-order -functions can be obtained from the recurrence relation which is well-conditioned for . For , the computation of -functions using the expressions above suffers from cancellation errors [11]

; therefore, we resort to their Padé approximants for

(). The –Padé approximant [11] is given by with

 N(ℓ)d=d∑j=0(j∑k=0νj,k)zj=d∑j=0νjzj,D(ℓ)d=d∑j=0δjzj, (10)

where the coefficients can be computed using the recurrence relations

 νj,k+1=−(d−k)(ℓ+j−k)(k+1)(2d+ℓ−k)νj,k,νj,0=1(ℓ+j)!,δj+1=(d−j)(j+1)(2d+ℓ−j)δj,δ0=1. (11)

For , we first scale the variable to , where is such that , and use the –Padé approximants stated above. Once the values of are available, can be computed by a scaling procedure given by [11]

 22kφ2k(2z)=φk(z)φk(z)+2k∑j=k+12φj(z)(2k−j)!,22k+1φ2k+1(2z)=[φk(z)+1k!]φk+1(z)+2k+1∑j=k+22φj(z)(2k+1−j)!.

The first few of these scaling relations are as follows: , and .

### Ii-C Polynomial approximation

For with , it is possible to develop polynomial approximations similar to that proposed in [12] for the functions which converge exponentially. This can serve as an alternative to the rational Padé approximation method presented above. In addition, this furnishes a direct proof of the fact that the approximation (8) is exponentially accurate. Employing the expansion where denotes the Legendre polynomial of degree and denotes the spherical Bessel function of order , we have

 φk(iξ)=∞∑n=0(2n+1)injn(σ)∫10τk−1Ln(1−τσξ)(k−1)!dτ=∞∑n=0(2n+1)injn(σ)Pn(ξσ;k),k≥1. (12)

The first few polynomials work out to be: , and . For , we have the recurrence relations

 P2m=2x(2m−12)2m+kP2m−1−(2m−k−1)2m+kP2m−2+(−1)m(2m−12)Γ(m−12)√π(2m+k)Γ(m+1)Γ(k),P2m+1=2x(2m+12)2m+1+kP2m−(2m−k)2m+1+kP2m−1, (13)

where we have suppressed the dependence on and for the sake of brevity.

Now, let denote the approximation using the first terms of the expansion in (12) and the corresponding error. From the inequalities and , and, setting (i.e. assuming ), we have where which confirms that the rate of convergence is exponential with respect to (see Fig. 1 for a numerical example).

### Ii-D Numerical tests

For the numerical experiments, we employ the chirped secant-hyperbolic potential given by () considered in [7]. We set the computational domain to be and let . Let denote the -domain; then, the error in computing is quantified by where is the numerical approximation and the integrals are computed using the trapezoidal rule. For the purpose of testing, we include the RK methods presented in [7] which are labelled as IF–ERK and IF–IRK. The results are plotted in Fig. 2 which comprises convergence analysis (top row) and the trade-off between accuracy and complexity (bottom row). The results show that the ETD methods tend to be more accurate that the IF methods but their complexity-accuracy trade-off is comparable. A superior multiplication algorithm for polynomials in the Chebyshev basis may, however, dramatically change this situation.

## References

• [1] M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, Part I,” IEEE Trans. Inf. Theory, vol. 60, no. 7, pp. 4312–4369, 2014.
• [2] S. K. Turitsyn, J. E. Prilepsky, S. T. Le, S. Wahls, L. L. Frumin, M. Kamalian, and S. A. Derevyanko, “Nonlinear Fourier transform for optical data processing and transmission: advances and perspectives,” Optica, vol. 4, no. 3, pp. 307–322, Mar 2017.
• [3] V. Vaibhav, “Nonlinearly bandlimited signals,” J. Phys. A: Math. Theor., vol. 52, no. 10, p. 105202, 2019.
• [4] ——, “Fast inverse nonlinear Fourier transformation using exponential one-step methods: Darboux transformation,” Phys. Rev. E, vol. 96, p. 063302, 2017.
• [5] ——, “Fast inverse nonlinear Fourier transform,” Phys. Rev. E, vol. 98, p. 013304, 2018.
• [6] ——, “Higher order convergent fast nonlinear Fourier transform,” IEEE Photonics Technol. Lett., vol. 30, no. 8, pp. 700–703, 2018.
• [7] ——, “Efficient nonlinear Fourier transform algorithms of order four on equispaced grid,” IEEE Photonics Technol. Lett., vol. 31, no. 15, pp. 1269–1272, 2019.
• [8] S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems,” J. Comput. Phys., vol. 176, no. 2, pp. 430–455, 2002.
• [9] M. Hochbruck and A. Ostermann, “Exponential Runge–Kutta methods for parabolic problems,” Appl. Numer. Math., vol. 53, no. 2, pp. 323–339, 2005.
• [10] P. Giorgi, “On polynomial multiplication in Chebyshev basis,” IEEE Trans. Comput., vol. 61, no. 6, pp. 780–789, 2011.
• [11] H. Berland, B. Skaflestad, and W. M. Wright, “EXPINT—A MATLAB package for exponential integrators,” ACM Trans. Math. Softw., vol. 33, no. 1, 2007.
• [12] A. Y. Suhov, “An accurate polynomial approximation of exponential integrators,” J. Sci. Comput., vol. 60, no. 3, pp. 684–698, 2014.