Periodic nonlinear Fourier transform for fiber-optic communications , Part I : Theory and numerical methods

In this work, we introduce the periodic nonlinear Fourier transform (PNFT) method as an alternative and efficacious tool for compensation of the nonlinear transmission effects in optical fiber links. In the Part I, we introduce the algorithmic platform of the technique, describing in details the direct and inverse PNFT operations, also known as the inverse scattering transform for periodic (in time variable) nonlinear Schrödinger equation (NLSE). We pay a special attention to explaining the potential advantages of the PNFT-based processing over the previously studied nonlinear Fourier transform (NFT) based methods. Further, we elucidate the issue of the numerical PNFT computation: we compare the performance of four known numerical methods applicable for the calculation of nonlinear spectral data (the direct PNFT), in particular, taking the main spectrum (utilized further in Part II for the modulation and transmission) associated with some simple example waveforms as the quality indicator for each method. We show that the Ablowitz-Ladik discretization approach for the direct PNFT provides the best performance in terms of the accuracy and computational time consumption. © 2016 Optical Society of America OCIS codes: (060.1660) Coherent communications; (060.2330) fiber optics communications; (070.4340) Nonlinear optical signal processing; (000.4430) Numerical approximation and analysis. References and links 1. G. P. Agrawal, Nonlinear Fiber Optics, 5th ed. (Academic, 2013). 2. M. Cvijetic and I. B. Djordjevic, Advanced Optical Communication Systems and Networks (Artech House, 2013). 3. R. J. Essiambre, R. W. Tkach, and R. Ryf, “Fiber nonlinearity and capacity: single mode and multimode fibers,” in Optical Fiber Telecommunications, Vol. VIB: Systems and Networks, 6th ed., Ch. 1, edts. I. Kaminow, T. Li, and A.E. Willner, pp. 1–43 (Academic, 2013). 4. P. J. Winzer, "Scaling optical fiber networks: challenges and solutions," Opt. Photon. News 26, 28–35 (2015). 5. E. Ip and J. M. Kahn “Compensation of dispersion and nonlinear impairments using digital backpropagation,” J. Lightwave Technol. 26, 3416–3425 (2008). 6. X. Chen, X. Liu, S. Chandrasekhar, B. Zhu, and R. W. Tkach, “Experimental demonstration of fiber nonlinearity mitigation using digital phase conjugation,” in Technical Digest of Optical Fiber Communication Conference and Exposition and the National Fiber Optic Engineers Conference (OFC/NFOEC), paper OTh3C.1 (2012). 7. S. L. Jansen, D. van den Borne, B. Spinnler, S. Calabro, H. Suche, P. M. Krummrich, W. Sohler, G. -D. Khoe, and H. de Waardt, "Optical phase conjugation for ultra long-haul phase-shift-keyed transmission," J. Lightwave Technol. 24, 54–64 (2006). 8. I. Phillips, M. Tan, M. F. Stephens, M. McCarthy, E. Giacoumidis, S. Sygletos, P. Rosa, S. Fabbri, S. T. Le, T. Kanesan, S. K. Turitsyn, N. J. Doran, P. Harper, and A. D. Ellis, “Exceeding the nonlinear-shannon limit using Raman laser based amplification and optical phase conjugation,” in Optical Fiber Communication Conference and Exhibition (OFC), Los Angeles, USA, paper M3C.1 (2014). 9. A. R. C. X. Liu, P. J. Winzer, R. W. Tkach, and S. Chandrasekhar, “Phase-conjugated twin waves for communication beyond the Kerr nonlinearity limit,” Nat. Photonics 7, 560–568 (2013). 10. D. Rafique, “Fiber nonlinearity compensation: commercial applications and complexity analysis,” J. Lightwave Technol. 34, 544–553 (2016). 11. V. E. Zakharov and A. B. Shabat, “Exact theory of two-dimensional self-focusing and one-dimensional selfmodulation of waves in nonlinear media,” Soviet Physics-JETP 34, 62–69 (1972). 12. M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, part I: mathematical tools,” IEEE Trans. Inf. Theory 60, 4312–4328 (2014). 13. M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, part II: numerical methods,” IEEE Trans. Inf. Theory 60, 4329–4345 (2014). 14. M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, part III: spectrum modulation,” IEEE Trans. Inf. Theory 60, 4346–4369 (2014). 15. J. E. Prilepsky and S. K. Turitsyn, “Eigenvalue communications in nonlinear fiber channels,” in Odyssey of Light in Nonlinear Optical Fibers: Theory and Applications, Ch. 18, edts. K. Porsezian and R. Ganapathy, pp. 459–490 (CRC, 2015). 16. A. Maruta, Y. Matsuda, H. Terauchi, and A. Toyota, “Digital coherent technology-based eigenvalue modulated optical fiber transmission system,” in Odyssey of Light in Nonlinear Optical Fibers: Theory and Applications, Ch. 19, edts. K. Porsezian and R. Ganapathy, pp. 491–506 (CRC, 2015). 17. A. Hasegawa and T. Nyu, “Eigenvalue communication,” J. Lightwave Technol. 11, 395–399 (1993). 18. S. T. Le, J. E. Prilepsky, and S. K. Turitsyn, “Nonlinear inverse synthesis technique for optical links with lumped amplification,” Opt. Express 23, 8317–8328 (2015). 19. A. Hasegawa and Y. Kodama, "Guiding-center soliton," Phys. Rev. Lett. 66, 161 (1991). 20. E. G. Turitsyna and S. K. Turitsyn, “Digital signal processing based on inverse scattering transform.” Opt. Lett. 38, 4186–4188 (2013). 21. J. E. Prilepsky, S. A. Derevyanko, and S. K. Turitsyn, “Nonlinear spectral management: linearization of the lossless fiber channel,” Opt. Express 21, 24344–24367 (2013). 22. J. E. Prilepsky, S. A. Derevyanko, K. J. Blow, I. Gabitov, and S. K. Turitsyn, “Nonlinear Inverse Synthesis and Eigenvalue Division Multiplexing in Optical Fiber Channels.” Phys. Rev. Lett. 113, 013901 (2014). 23. S. T. Le, J. E. Prilepsky, and S. K. Turitsyn, “Nonlinear inverse synthesis for high spectral efficiency transmission in optical fibers,” Opt. Express 22, 26720–26741 (2014). 24. H. Bülow, “Experimental demonstration of optical signal detection using nonlinear fourier transform.” J. Lightwave Technol. 33, 1433–1439 (2015). 25. V. Aref, H. BÃijlow, K. Schuh, and W. Idler, "Experimental demonstration of nonlinear frequency division multiplexed transmission," in European Conference on Optical Communication (ECOC), Valencia, Spain (2015). 26. H. Bülow, V. Aref, and W. Idler, "Transmission of waveforms determined by 7 eigenvalues with PSK-modulated spectral amplitudes," preprint arXiv:1605.08069 (2016). 27. M. I. Yousefi and X. Yangzhang, "Linear and nonlinear frequency-division multiplexing," preprint arXiv:1603.04389v2 (2016). 28. Z. Dong, S. Hari, T. Gui, K. Zhong, M. I. Yousefi, C. Lu, P. A. Wai, F. R. Kschischang, and A. Lau, "Nonlinear frequency division multiplexed transmissions based on NFT," IEEE Photon. Technol. Lett. 27, 1621–1623 (2015). 29. A. Geisler and C. G. Schaeffer, "Implementation of eigenvalue multiplex transmission with a real Fiber Link using the discrete nonlinear Fourier spectrum," in Proceedings of Photonic Networks (2016), Leipzig, Germany, pp. 1-6. 30. A. Maruta, “Eigenvalue modulated optical transmission system (invited),” in The 20th OptoElectronics and Communications Conference (OECC), Shanghai, China, Paper JThA.21 (2015). 31. S. Turitsyn, “Nonlinear Fourier transform for optical communications,” in Frontiers in Optics 2015, OSA Technical Digest (online), paper FM1E.6 (2015). 32. S. Hari, F. Kschischang, and M. Yousefi, “Multi–eigenvalue communication via the nonlinear Fourier transform,” 27th Biennial Symposium on Communications (QBSC), pp. 92–95 (2014). 33. S. Wahls, S. T. Le, J. E. Prilepsky, H. V. Poor, and S. K. Turitsyn, “Digital backpropagation in the nonlinear Fourier domain,” in Proceedings of IEEE 16th International Workshop in Signal Processing Advances in Wireless Communications (SPAWC), Stockholm, Sweden, pp. 445-449 (2015). 34. S. T. Le, J. E. Prilepsky, M. Kamalian, P. Rosa, M. Tan, J. D. Ania-Castanon, P. Harper, and S. K. Turitsyn “Modified nonlinear inverse synthesis for optical links with distributed Raman amplification,“ in European Conference on Optical Communication (ECOC), Valencia, Spain, paper Tu.1.1.3 (2015). 35. A. Osborne, Nonlinear Ocean Waves and the Inverse Scattering Transform, 1st ed. (Academic, 2010). 36. O. R. Its and V. P. Kotlyarov. “Explicit formulas for the solutions of a nonlinear Schrödinger equation,” Doklady Akad. Nauk Ukrainian SSR, ser. A, 10, 965–968 (1976). 37. E. R. Tracy and H. H. Chen, “Nonlinear self-modulation: an exactly solvable model,” Phys. Rev. A 37, 815–839 (1988). 38. C. F. Elliott, Handbook of Digital Signal Processing: Engineering Applications (Academic, 2013). 39. Q. Zhang, T. H. Chan, and A. Grant, “Spatially periodic signals for fiber channels,” in IEEE International Symposium on Information Theory (ISIT), Honolulu, HI, USA, pp. 2804–2808 (2014). 40. S. J. Savory, “Digital filters for coherent optical receivers,” Opt. Express 16, 804–817 (2008). 41. G. P. Agrawal, Fibre-Optic Communication Systems, 4th ed. (Wiley, 2010). 42. P. Poggiolini A. Carena, V. Curri, and F. Forghieri, "Evaluation of the computational effort for chromatic dispersion compensation in coherent optical PM-OFDM and PM-QAM systems," Opt. Express 17, 1385–1403 (2009). 43. W. Shieh, W. Chen and R.S. Tucker, “Polarisation mode dispersion mitigation in coherent optical orthogonal frequency division multiplexed systems,“ Electron. Lett. 42, 996–997 (2006). 44. Y. Ma and M. J. Ablowitz, “The periodic cubic Schrödinger equation,” Stud. Appl. Math. 65, 113–158 (1981). 45. S. Wahls and H. V. Poor, "Fast numerical nonlinear Fourier transforms," IEEE Trans. Inf. Theory 61, 6957–6974 (2015). 46. A. Bobenko and C. Klein, eds. Computational approach to Riemann surfaces, No. 2013, (Springer Science and Business Media, 2011). 47. B. Deconinck, M. S. Patterson, and C. Swierczewski, “Computing the Riemann constant vector,” preprint, available online at https://depts.washington.edu/bdecon/papers/pdfs/rcv.pdf (2015). 48. J. Frauendiener and C. Klein, “Computational approach to hyperelliptic Riemann surfaces,” Lett. Math. Phys. 105, 379–400 (2015). 49. B. Deconinck, H. Heil, A. Bobenko, M. Van


Introduction
Ever-increasing demand on capacity of communication systems due to the immensely rapid development of versatile on-line services (such as cloud computing, on-demand HD video streaming, on-line games, Internet of Things, etc.) motivates the active search of new approaches capable to increase the data transmission rates in optical networks.This particularly refers to the fiber-optic communication systems, currently carrying an overwhelming fraction of the global data traffic [1][2][3][4].Whilst the "fifth generation" of optical transmission systems, operating with coherent detection, advanced multilevel modulation, and digital signal processing (DSP), has led to the possibility of channel rates exceeding 100 Gb/sec [2][3][4], the current approaches are mostly based on thoroughly-studied techniques borrowed from linear medium, like, e.g., wireless systems.However, when addressing the fiber-optic transmission, one has to deal simultaneously with four different physical effects.These are: signal attenuation along the propagation, signal dispersion, noise (arising due to the amplifier spontaneous emission like the inevitable concomitant of amplification), and the inherent nonlinearity of the fiber medium itself.The availability of efficient methods for the fiber loss compensation, like the lumped erbium-doped amplifiers (EDFA) or Raman distributed amplification (RDA) schemes [1], allows one to somehow simplify the problem and deal only with the signal corruptions induced by the three remaining effects.Moreover, insofar as the mitigation of chromatic and polarization mode dispersion (linear effects) has been already successfully deployed in contemporary transmission systems [1,2], noise and nonlinearity are currently considered to be the major factors limiting the performance of optical networks.A number of various nonlinearity compensation methods, including digital back-propagation (DBP) [5], digital [6] and optical [7,8] phase conjugations (spectral inversion), and phase-conjugated twin waves [9], have been proposed (to mention just a few important advances, see the reviews in e.g.[3,4,10]).However, we note that within the aforementioned techniques are mainly focused on compensating the detrimental effects of linear and nonlinear distortions while within the nonlinear Fourier transform (NFT, originally called the inverse scattering transform, IST [11]) based methods nonlinearity serves as an inherent element of the processing method; see the recent reviews [12][13][14][15][16].In fact, application of this method to optical communications originated from the paper by Hasegawa and Nyu "Eigenvalue communication" [17], where the (relatively simple) modulation of the so-called nonlinear spectrum (or, rather of its solitonic part) was first proposed.
For this proof-of-concept study we deal with the lossless single-polarization nonlinear Schrödinger equation (NLSE) perturbed by the additive white noise as a model for our signal evolution down the optical fiber (see discussion on the applicability of the lossless model, e.g. in [18,19]).The NLSE (together with its generalizations) is a principal master model governing the evolution of the slow-varying optical field envelope along a single-mode fiber [1,3].The pure NLSE (in the absence of losses and noise) belongs to the special class of exactly solvable equations [11] (the so-called integrable equations), the solution of which can be found by means of the NFT method, also widely termed as IST technique in physical and mathematical literature.The NLSE solution by NFT is similar to the application of linear Fourier transform (FT) for finding the solution of a linear partial differential equation given the initial conditions: it converts both the chromatic dispersion and nonlinearity into the linear propagation of uncoupled "nonlinear modes" evolving inside the so-called nonlinear spectral domain in the linear manner [11,15,17].To find the nonlinear spectrum (NS), serving the role of FT spectrum for nonlinear integrable equations, one needs to transform the signal waveform (e.g. the signal at the transmitter) into the nonlinear spectral domain, by performing the direct NFT; for the NLSE the latter is called the Zakharov-Shabat system (ZSS) or scattering problem.For the infinite time interval and vanishing (return-to-zero) signals at the time-domain "boundaries", the NS can generally contain the soliton eigenvalues (non-dispersive part of the solution) and continuous part (describing the dispersive nonlinear waves).In general, it is now well recognized that the NFT method gives us an opportunity to introduce new non-trivial concepts for data transmission.The details of the NFT in applications to optical communications can be found in a great number of recent works [12-16, 18, 20-34].
Despite the numerous recent results and a good potential of the NFT-based methods, there are still many problems associated with applying the traditional infinite-line NFT (attributed to vanishing or return-to-zero signals) to the optical transmission.One of the challenges for standard NFT is the large required processing window.The transmission is supposed to be in the burst-mode [18,23], such that the initial signal is appended with the "wings" usually exceeding the dispersion induced memory.Other challenges include indirect encoding [32] and non-trivial optical signal generation [22,23,32], alongside with some problems associated with the simultaneous processing of the continuous and discrete parts of the NS [33].For the nonlinear inverse synthesis (NIS) idea [18,22,23,34], where the continuous spectrum is used for the modulation, one cannot completely control the time domain occupation of the synthesized signal which time extent usually exceeds that of the initial modulated signal (see Fig. 2 in [22] and the discussion in [23]).The common cause of these problems is the decaying constraint for the solutions of NLSE, which are used for the NFT-based transmission purposes.However, at least some of the afore-mentioned difficulties of the ordinary NFT can be rectified by using the periodic NFT-based processing.
As it will be explained further and in Part II, periodic NFT (PNFT) technique can be more fit to conventional signal processing frameworks (say, built on the inherently periodical FFT processing) and also allows one to cope better with other problems attributed to the NFT-based processing, like a high peak-to-average power ratio (PAPR) and the absence of control over the time domain signal occupation.In this work we introduce PNFT for optical communications and use the new NS, attributed specifically to periodic (in time) solutions of NLSE.Further in Part II we will consider a communication system based on the characteristics of the NS attributed to the PNFT.
The PNFT has been so far mostly studied and developed with regard to the analysis of water (ocean) waves dynamics [35].The mathematical aspects of the periodic IST theory have also been the subject of intensive studying, see e.g.[36,37].In [36] a general structure for the so-called finite-gap solutions of the periodic NLSE is proposed, while in [37] the authors seek for periodic solutions in form of a perturbation to a base solution and derive the analytical expression for the case of perturbed plane wave.The monograph by Osborne [35] follows the lead of afore-mentioned works in a comprehensive review of PNFT properties and applications and equips it with various numerical algorithms to perform the PNFT operations.In this Part I we itemize the mathematical aspects of the PNFT and provide a base for further application of the PNFT theory to fiber-optic communications in the following Part II.Also, in Part I we address the available numerical methods for the direct PNFT computation.
This paper is organized as follows.In Section 2 we compare PNFT and ordinary NFT (for signals vanishing at the boundary, that in what follows we will call a burst-mode NFT) in terms of their suitability for the communication purposes and underline the advantages of the PNFT.After that, in Section 3 we discuss the basic concept of periodic signal processing.Then, in Section 4 we introduce the explicit form of PNFT operations and classify the parts of the associated NS.Using these data, in the next Section 5 we analyze the performance of four numerical methods for finding the periodic NS of two example waveform where the analytical results are available.For completeness, in this section we also briefly mention the methods applicable for the calculation of the inverse PNFT operation though it will not be used further in our Part II.At the end we link our results to the following paper by scrutinizing the communication characteristics of the PNFT-based systems.The Appendix provides some explicit data for the NS calculation of a periodic rectangle pulse train (we note that the latter has not been given in the literature, up to our knowledge).

Channel model
We first introduce PNFT for a simplified NLSE model, with the chromatic dispersion coefficient β 2 < 0 and the Kerr nonlinearity coefficient γ (see e.g.[1][2][3] for details): where q = q(t, z) is the slow-varying envelope of the signal, z stands for distance along the fiber and t denotes the retarded time in the frame co-moving with the envelope [1].This equation describes the interplay of linear dispersion, responsible for signal broadening in time domain, and Kerr nonlinearity, causing the signal spectrum broadening.The dispersion brings about two major effects, namely, the inter-symbol interference (ISI) due to the spreading of the initial symbols, and it simultaneously distorts the waveform of our initial signal.

ISI mitigation
In a noiseless system the received signal at the end of the fiber is given by Here x m is the m th input symbol, q 0 (t) is the pulse shape, T is the time interval between two consecutive symbols, h 0 (t) is the impulse response of fiber, and the sign " * " implies convolution.The received signal y m (t) is typically dispersed over the interval larger that T due to the dispersion-induced channel memory, and therefore the ISI takes place.In a linear regime (when signal power is small) in a coherent detection system, one can remove the chromatic dispersion by filtering the signal using a filter whose impulse response is governed by the transfer function, h(t) [40,41]: where L is the fiber length, and H 0 (ω) and H(ω) are the FT of h 0 (t) and h(t), respectively.Given the whole signal y(t), in the absence of noise this filtering removes the ISI completely (Fig. 1).To implement this all-pass filter as a time-domain equalizer one has to truncate the impulse response to a finite duration and then execute the filter as a FIR filter.For a large number of samples and a large filter length, it is more efficient to implement the filter in the frequency domain [5].Although, in either form, one has to filter the whole stream of consecutive symbols for recovering the signal by considering a finite channel memory T d , it is possible to merely process the extended window around the target symbol.This approach is called overlap-and-save, where the input signal is divided into some overlapping segments which, after having been filtered, are truncated and cascaded to construct the continuous filtered signal (Fig. 1).The procedure of removing ISI in overlap-and-save processing is as follows [38].After windowing the input signal with rectangular window W (t) of unit height centred at the n-th symbol, having a total duration T w = T + 2T d , and assuming that the channel memory is less than the initial signal duration T , we have where the only important property of y r n−1 (t) and y r n+1 (t), which are defined through Eq. ( 4) as the non-interfering part of (n − 1)-th and (n + 1)-th symbol, respectively, is that the closest non zero part of them is at least T d far from two ends of xn (t); the latter quantity is defined by the expression below.The next step is to filter ŷ(t) to eliminate the signal broadening: where x r n−1 (t) and x r n+1 (t) are the results of filtering y r n−1 (t) and y r n+1 (t), respectively.The last stage in recovering the signal xn (t) = x n q 0 (t − nT ) is to select the middle part of x(t) (with the length T ) while discarding the rest.Note that, since h(t) is a filter with the same memory as h 0 (t), its signal broadening does not lead y r n−1 (t) and y r n+1 (t) to expand to the middle part.Thus, when we pass them through h(t), the result remains in the discarded part of the signal and does not participate in xn (t).
This filtering to combat the ISI comes with a delay and costs of implementing a filter [40,42], and the alternative is to insert a guard interval between the adjacent information-bearing signals to prevent the ISI at the expense of losing some effective data rate (a burst mode).The latter is especially effective for high spectral efficiency (SE) modulation formats such as orthogonal frequency division multiplexing (OFDM) [42,43].We note that the chromatic dispersion (CD) compensation in OFDM by using the cyclic prefix is supposed to eliminate the linear distortions of the signal, whereas in NFT formalism we merely need to mitigate the ISI to collect all the samples of the symbol.To be more precise, the linear distortion causes two major effects: the widening of the pulse (leading to ISI) and the change of the signal waveform.To recover the transmitted signal in conventional optical communication system one needs to restore the original state of the pulse for the sake of getting back the encoded data.This is rendered via removing the linear distortion using the techniques such as the one mentioned earlier.On the other hand, within NFT formalism, since the evolution of the NS is linear, one can ignore the linear filtering at the receiver to compensate this linear distortion.However, to perform the NFT one has to get back the broadened pulse without any interference, i.e. to mitigate the ISI (see Refs. [12,22,23,32] for the details of the NFT-based methods).The simple filtering is not capable to restore the signal in the presence of nonlinearity, despite the fact that for the anomalous dispersion the overall dispersion gets effectively reduced due to the nonlinearity action [41].The back-propagation (BP) technique was proposed to compensate both linear and nonlinear distortions by means of solving the NLSE in a backward direction [5].However, the processing complexity of BP is considerable because the processing window has to be kept larger than the pulse duration plus twice the channel memory (Fig. 2).Moreover, the BP technique is sensitive to high nonlinearity and/or high powers, when the degradation due to the nonlinear interaction between the signal and noise is essential [5].From the other side, by using cyclic extension instead of empty intervals, the processing window can be made as small as the length of the data-bearing part of the signal.This is due to the fact that after removing the cyclic prefix, the data-bearing part of the signal is not affected by the ISI due to CD and thus, still contains the complete (though, distorted) information about the transmitted signal.As a result, the computational complexity can be significantly reduced due to the reduced processing widow, see Fig. 3 for the visualization of the processing steps.At the same time it is worth mentioning that in soliton-based communication there is no any dispersive-type signal broadening so no special provision is required to avoid ISI [25,26].

Benefits of the PNFT over the ordinary NFT
In digital signal processing, even when we work with aperiodic signals, we basically assume periodicity when using the discrete FT as the most widely used computation tool.Thus, considering periodic signals can be reckoned as a "natural choice" (see, e.g., Chapter 1 of [35]).Together with this there are some additional benefits in the PNFT utilization (as opposed to the ordinary NFT) that are listed below.
• As explained in the previous subsection, the use of periodic signals with cyclic extension reduces the processing window at the receiver side excluding the "wings" having the double dispersion memory extent, 2T d .Specifically, since the channel memory depends on the data rate [42], the amount of two guard intervals difference between periodic and vanishing signals can be considerable.
• Within the PNFT, due to its controllable constant time duration, we can have a projection from the modulated signal in the NS domain into the space-time domain with the predefined characteristics.Mapping data on the NS could be done by using the inverse PNFT, similarly to the NIS idea [22,23]: At the transmitter side one synthesizes a signal from some modulated (encoded) NS profile and then launches it through the fiber.Since for the ordinary NFT the signal characteristics for different modulations and power levels (burst energies) inside the NS domain usually have different and weakly controllable time-and frequency-domain occupation, this data mapping is complicated and not straightforward [22,23,32] (note that in the latter Ref.only nearly 2,500 signals from over 2.3 millions of possible realizations were taken to reach the 3.14 bit/s/Hz SE).Moreover, in the PNFT for some signals constructed from special kind of constellation, the bandwidth of signals have small variation, which helps to make the projection simpler.Hence, the PNFT encoding schemes can be built by means of the encoders of currently used communication systems and, due to the similarities between the PNFT-generated signals and OFDM ideology in special cases, specifically with the encoders of the OFDM systems (see Part II and [37]).
• The periodic solutions of NLSE can be expressed in terms of Riemann theta functions, see Section 5; the latter can be regarded as a direct generalization of the linear FT with timevarying coefficients with some existing efficient fast algorithms for its computation [35].
• Another benefit of using the periodic signals over the vanishing ones is that by using the PNFT-based processing we deal with a continuous stream of signals without sudden (perhaps, very wide) gaps in power distribution, so that the PAPR of PNFT-based transmission methods is considerably lower than that of the signals used in ordinary NFT-based transmission methods.

NFT for vanishing signals
In this section, first we recall some basics of the ordinary NFT for the sake of our further introduction of its periodic counterpart and the comparison of both methods.We provide just some basic elements of NFT; more details can be found in [12-15, 21-24, 30].We write down the normalized unperturbed NLSE as where the variables are normalized according to In ( 7) T s is an arbitrary characteristic time scale (a free parameter), and Z s = 2T 2 s /|β 2 | is the associated spatial scale.

Direct NFT operations
As the NLSE ( 6) is integrable, it can be recast as a compatibility conditions for two linear equations, the so-called Lax pair [11].The first equation from this pair corresponds to the direct NFT operation, i.e. to the ZSS.The ZSS defines the decomposition of our profile in time domain, say at z = 0, q(t, 0) = q(t), into the set of nonlinear spectral data, or, in other words, projects q(t) onto the NS.We note that for the ordinary NFT we explicitly assume that |q(t, 0)| → 0 as |t| → ∞, and the decay of q(t) is stronger than any power of t (for the processing one usually truncates q(t) to a finite duration).The second equation of the Lax pair defines the evolution of the NS with z; the evolution of scattering data appears to be trivial and linear [11].
The ZSS is given by the set of two linear equations where q(t) enters as an effective potential: Here φ 1 and φ 2 are the auxiliary eigenfunctions and λ is the (generally complex) spectral parameter.Solutions to this equation depend on the chosen conditions in t-variable.In ordinary NFT, the eigenfunctions tend to the exponential type solution [φ 1 , φ 2 ] T ∼ [e −iλt , e iλt ] T for large |t|, when |q(t)| → 0. For real λ we can fix two special solutions of ZSS (8) as The two sets of linearly independent solution corresponding to Ψ i , where the "boundary condition" for each component Ψ k i is defined by Eq. ( 9), are: Now one defines the scattering coefficients (also called Jost coefficients or NFT amplitudes), a(λ ) and b(λ ), forming the core of the direct NFT, by expressing Ψ 2 through Ψ 1 as t → +∞: The NS associated with the arbitrary vanishing profile q(t) generally consists of two parts.The reflection coefficient defined for real λ , r(λ ), is the continuous part of NS referring to the dispersive nonlinear waves The discrete spectrum consists of complex eigenvalues λ n lying in the upper complex halfplane of λ , C + , defined by the condition a(λ n ) = 0 (i.e.λ n corresponds to the pole of r(λ )).Each λ n corresponds to a soliton component of our solution (we assume that each λ n is nondegenerate).Each soliton is also characterised by the second complex-valued parameter, the norming constants γ n , defined as γ n = b(λ n )/a (λ n ), where the derivative is taken with respect to λ .The set Σ = {r(λ ), λ n , γ n |λ ∈ R, λ n ∈ C + } forms the NS for a vanishing signal.The evolution of NS components (modes) along the fiber is uncoupled and linear: while λ n remain constant, λ n (z) = λ n (z 0 ) for any z and z 0 .

Inverse NFT
The INFT gives the inverse mapping of the NS Σ onto the time domain producing the profile q(t) (or q(t, z)).The INFT is defined through the solutions of the Gelfand-Levitan-Marchenko equations for the unknown functions M 1,2 (τ, τ ) [11]: The function R(τ) is defined via the NS Σ and can contain contributions from both solitonic (discrete) and radiation (continuous) spectrum parts, R(τ) = R s (τ) + R rd (τ), where we have assumed that all discrete eigenvalues have a unit multiplicity.Having solved Eqs. ( 14) for M 1,2 (τ, τ ), the sought solution in the space-time domain is recovered as q(t) = −2M 1 (t,t) (or the same expression with the z-dependence of functions).

Definition of the direct PNFT
The NFT approach in the case of periodic signal with some period T 0 , q(t + nT 0 , 0) = q(t, 0) for any integer n, brings about the signal decomposition into the different kind of NS, that now consists of the so-called main spectrum and auxiliary spectrum.A finite-band solution for the NLSE is the one described by the finite number of parameters in main spectrum, e.g. by the set of complex numbers λ 1 , λ 2 , . . .∈ C. Auxiliary spectrum can be defined in two senses, namely, in the sense of Kotlyarov [36] and in the sense of Ma and Ablowitz [44]; each definition is suitable for constructing the periodic finite-gap NLSE solution [45].Here, since in Part II we will use only the main spectrum for the modulation and data transmission, we do not give a lot of details for the auxiliary spectrum properties; see [35][36][37]44].
To identify the main spectrum we again employ the ZSS, Eq. ( 8), but now we assume that the "potential" is a periodic function of t.We choose an arbitrary base point t = t 0 and introduce two independent solutions of ZSS, φ (t,t 0 ; λ ) and φ (t,t 0 ; λ ), which at t = t 0 have the values φ (t 0 ,t 0 ; λ ) = [1, 0] T and φ (t 0 ,t 0 ; λ ) = [0, 1] T .The 2×2 fundamental matrix for Eq. ( 8) is defined as At t = t 0 , Φ(t 0 ,t 0 ; λ ) is equal to the rank 2 identity matrix.The monodromy matrix is then defined as M(t 0 ; λ ) = Φ(t 0 + T 0 ,t 0 ; λ ) which is used to find the NS of the periodic signal.The main spectrum is defined as the solutions (in λ variable) of the equation [35] which corresponds to periodic or anti-periodic eigenfunctions.The main spectrum is the analogue of the discrete spectrum for the ordinary NFT.The important property of the main spectrum is that it remains invariant during the pulse evolution along z-direction.It is this property that will be utilized in Part II for the encoding transmission of data, similar to the invariant soliton spectrum used in the ordinary NFT-based eigenvalue communications idea [17,32].
For the test of the numerical methods we will need the exactly solvable periodic ZSS cases, where the main spectrum is known.As a first exactly-solvable case we can take the plane wave solution, q(t, z) = q 0 e iµt , for which the main spectrum consists of complex points [45]: The main spectrum of a plane wave with q 0 = 5, µ = 3 is shown in the Appendix, Fig. 6 (for the comparison reasons).Another solution where the explicit expression for the main spectrum can be obtained is the periodic sequence of rectangular pulses, see the Appendix for the calculations.Both these periodic waveforms attributed to the exactly-solvable ZSS will be used in the next section for the comparison of different numerical methods.
The auxiliary spectrum of a signal with g−1 modes is defined at any time and space point as Time and space evolution of auxiliary spectrum in governed by the following equations: where σ j ∈ {+1, −1} is a Riemann sheet index.

Inverse PNFT
Now we briefly outline the procedure of how to retrieve the (periodic) profile q(t) starting from the given NS distribution (assuming that we know both two kinds of periodic nonlinear spectral data).There are several methods to construct finite-gap periodic solutions of NLSE using NS from which, perhaps, the most studied one is Darboux-Bäcklund transforms vastly used for burst-mode signal construction [13,45].However, it is more convenient to use the theta-function representation to form the temporal profile of the solution due to the guaranteed periodicity of the resulting signal.One can reconstruct the solution of NLSE by means of squared ZSS eigenfunctions in a way explained in [36,37], which for the finite-gap case leads to the solutions of NLSE in the following form [35,37]: where k 0 and ω 0 are real constants which can be obtained from the NS data, and we have introduced the Riemann theta function, Θ(W|τ), which is the multidimensional generalization of the ordinary FT [35]: Here g is a finite number defining the number of modes for our g-gap solution, m is a point of the (g − 1)-dimensional integer lattice, and W ± = π(kz +ω ω ωt +δ δ δ ± )/2 is a vector calculated from NS.The set {k,ω ω ω,δ δ δ , τ} is called Riemann spectrum and τ is the Riemann (period) matrix [35]; their particular values can be, again, obtained from the NS data.Thus, the inverse PNFT procedure for a finite-gap solution can be reformulated as the problem of finding the Riemann spectrum from the given two types of NS, namely from main and auxiliary spectra.Then one uses the Riemann spectrum obtained for constructing the solution via Eq.( 21), i.e. using the Riemann theta functions, Eq. ( 22).The procedure of finding the Riemann spectrum from NS is as follows: • Since the solutions of Eqs. ( 20) lie on a specific Riemann surface, we can change the variables in such a way that the auxiliary parameters have a trivial dependency on z and t [37,46].So, at this stage one has to construct the Riemann surface and define the appropriate differential equations on it to carry out integrations rendering the solutions of Eqs.(20) on this surface.
• Using the above-mentioned differential equations and by defining proper loops on the surface (paths of integration), one can find the Riemann matrix τ by calculating some loop integrals (for more details see [46]).
• From these loop integrals it is possible to get the simple evolution of the auxiliary spectrum recast in a form containing k, ω ω ω, and δ δ δ .
Although there is still currently a lack of a generic approach for how to deal with the inverse PNFT, there are several software packages and routines allowing one to construct the profile in time domain using the periodic NS data.For the first stage of the inversion, i.e. for one's finding the Riemann spectrum, there are some useful packages and codes embedded in mathematical software like Maple, Sage, and Mathematica [46][47][48].For the second stage, which is to construct the Riemann theta functions (22) using the Riemann spectrum, in addition to the symbolic implementations [49,50], some "hyper-fast" algorithms for the numerical reconstruction of the signal were proposed [35].

Numerical methods to for the main PNFT spectrum calculation
In this section we compare the performance of four numerical methods for the calculation of the main spectrum (used further in Part II for the transmission example), see also Ref. [45].As a quality indicator for each method we use the deviation from the analytical expressions for the main spectra of two waveforms, namely for that corresponding to the plane wave and the periodic train of rectangle pulses.For these waveforms the periodic ZSS, Eq. ( 8), can be solved analytically, see Eq. ( 18) and Appendix.We do here with four methods: the layer peeling method [55], Ablowitz-Ladik integrable discretization method, Crank-Nicolson method and spectral Fourier collocation method [52].The first three methods basically deal with finding the elements of the monodromy matrix of ZSS (8) as a function of spectral parameter λ .Then one has to append the computation with a root-finding routine, such as the Newton-Raphson method.
The last Fourier collocation method recasts the ZSS into an eigenvalue problem in the Fourier space.For convenience reasons we rewrite the ZSS (8) as

Layer peeling algorithm
For this method we discretize q(t) in time at points t n = n ∆t, q n = q(t n ), where the time interval t n+1 − t n = ∆t = T 0 /N, n = 1, 2, • • • N, and T 0 is the signal period.Now we write down the ZSS solution at the end of each segment as Ψ(t n + ∆t) = U(q n )Ψ(t n ), where and k 2 = |q n | 2 + λ 2 is constant inside each interval ∆t [13,35].At the end of one period the monodromy matrix is then approximately given by the expression:

Crank-Nicolson method
Approximating the time derivative in Eq. ( 23) with the central difference as [13,45] leads to the expression for the eigenfunction at the (n + 1)-th step: From this expression we can obtain the approximation of the monodromy matrix as It is worth mentioning that the structure of matrix P(q) allows one to consider the elements of the monodromy matrix as rational polynomial expressions [45], which substantially simplifies computations.It also means that the method can be recast into a fast algorithmic form by using the fast polynomial multiplication [45].

Ablowitz-Ladik discretization method
Using the approximation e ±i∆tλ ≈ 1 ± i∆tλ and introducing ω = e i∆tλ , we arrive at the Ablowitz-Ladik approximation of the ZSS (26), where after each iteration we have [13,45]: In this formula α n = 1 + ∆t 2 |q n | 2 .The resulting monodromy matrix is given by the expression As same as we had for the Crank-Nicolson method, here we can use the fast polynomial multiplication algorithms to simplify the computations and recast the method's implementation into the fast form by using fast polynomial arithmetic [45].Thanks to the structure of Eq. ( 29), the elements of the monodromy matrix have only even (or odd) powers; this doubles the speed of computations in comparison with the Crank-Nicolson method.In contrast to the Crank-Nicolson method where the accuracy depends on the approximation of the derivative operation only, the Ablowitz-Ladik method also suffers from the approximation in the change of variable λ → ω.However, since the order of the associated polynomial for the Ablowitz-Ladik method is half of that for the Crank-Nicolson method, the accuracy of the former can be even better than that of the latter.

Spectral Fourier collocation method
For a linear differential equation of the form ẏ = A(t)y, where A(t) is a piecewise continuous periodic matrix with period T 0 (such as Eq. ( 23)), the fundamental matrix (the matrix formed by two independent solutions of the equation) is represented as Φ(t) = Φ(t)e Rt , with R being a constant matrix and Φ(t) = Φ(t + T 0 ) is a periodic matrix.Introducing the new eigenfunctions where µ ∈ [0, 2π T 0 ], one can explicitly rewrite our original problem (23) in the form [53]: Solving this eigenvalue problem for µ = 0 and µ = π/T 0 (corresponding to periodic and antiperiodic eigenfunctions, respectively) gives the main spectrum of the signal.Since our signal is periodic, we can substitute for the periodic functions their Fourier series [52]: where k = 2π T 0 .In this way we obtain for the ZSS the eigenvalue problem: and Q is a matrix with elements q i,i− j = q j for i = 1, • • • , 2N if −N +i ≤ j < i, and zero otherwise.The main spectrum is obtained by solving the eigenvalue problem (34).The accuracy of this method is spectral due to the use of Fourier expansion: For smooth enough functions, the error in approximating the potential by a Fourier series decays exponentially [52].To reconstruct the main spectrum of our periodic signal we just need to use the two values, µ = 0 and µ = π/T 0 , and the corresponding eigenvalues [51].The computational complexity of this method is mainly attributed to finding the eigenvalues of a complex-valued non-Hermitian matrix, which, by using ordinary algorithms, is O(M 3 ), where M is the dimension of the matrix in Eq. ( 34); here it is equal to 4N.Thus, one of the disadvantages of this method is its large computational complexity.18), with q 0 = 5, µ = 3 for Ablowitz-Ladik, layer peeling (on this method's runtime refer to the text), Crank-Nicolson and spectral methods vs. the number of temporal samples.

Comparison of the direct PNFT algorithms in terms of accuracy and runtime for the main spectrum computations
Using two example signals, namely the plane wave and the periodic train of rectangles, with the known man spectrum (see Section 4 and Appendix), we compare the accuracy and runtime of the four methods described in the previous subsection.The method's quality measure is defined as the difference between the ideal (theoretical) and computed main spectrum of these signals.This difference is obtained by averaging the Euclidean distance between the analytically predicted and computed spectrum points.For the plane wave with q 0 =5, µ =3, the error profiles of four numerical methods are depicted in Fig. 4. One can see that the best performance for this case is attributed to the layer peeling and spectral methods.However, considering the runtime (the right panel of Fig. 4) and since the computational complexity of layer-peeling method with a search area of size K and number of samples N is O(NK), the advantage of using Ablowitz-Ladik and Crank-Nicolson methods is evident.From Fig. 4 we can deduce that considering run-time and accuracy, the Ablowitz-Ladik method is the most efficient one.Fig. 5 shows the error of finding the main spectrum by means of four methods, and their normalized runtime, for a rectangular pulse train with period T 0 = 2π, duration T = 2 and amplitude A = 3.From Figs. 4-5 one can conclude that the layer peeling and spectral methods are not convenient in the numerical PNFT computation for the large number of samples due to their significant computational complexity.On the other hand, the Ablowitz-Ladik method seems to be the best numerical routine for its acceptable accuracy and excellent computational load, so in Part II we use the Ablowitz-Ladik method for the decoding of the modulated main spectrum.

Conclusion
In the Part I we introduced the basic idea behind implementation of the PNFT in fiber-optic communication systems.We outlined the potential advantages of the PNFT over the ordinary NFT, in particular, the processing load and signal generation, lower PAPR level, and others.Yet another important advantage of the PNFT is that it allows one to keep the control over the time-domain signal extent, while for the ordinary NFT the time duration of the synthesized signals is not usually known apriori (this refers to both the multisoliton [32] and continuous spectrum based [18,22,23] methods).This means that we can potentially optimize the PNFTbased methods more efficiently in order to achieve a high spectral efficiency of the transmission.First we note that the ZSS.( 8) can be solved exactly for a constant "potential" q(t) = A, as it was used in the layer peeling method, subsection 5.1.We divide the rectangular pulse into three intervals, T 1 : −T 0 /2 < t < −T /2, T 2 : −T /2 < t < +T /2, and T 3 : T /2 < t < T 0 /2.For the intervals T 1 and T 3 the signal is zero and has the constant amplitude A inside T 2 .After some straightforward algebra we arrive at the expression for the monodromy matrix M rec in the form: To find the main spectrum one has to find the zeros of Tr M rec (λ ) = ±2.The main spectrum for the train of rectangular pulses with the individual rectangle duration T = 2, period T 0 = 2π and amplitude A = 3, is shown in Fig. 6 (b).As can be seen in this figure, the main spectrum consists of the set of points with small imaginary parts and of a few purely imaginary ones.In the limit of large period [51], the former points approach the real axis to form the continuous NS part of the ordinary NFT spectrum, while the imaginary points describe two soliton eigenvalues nucleated from a single rectangle with a given extent and duration.

Fig. 1 .
Fig. 1.ISI mitigation in linear regime by using the linear filtering.

Fig. 2 .
Fig. 2. ISI mitigation in nonlinear regime with three neighboring pulses using BP.

Fig. 3 .
Fig. 3. Processing window in the case of a burst-mode signal and periodic signal with cyclic extension.

Fig. 4 .
Fig.4.a) The error, and b) the normalized (to the number of samples) runtime, for finding the main spectrum of the plane wave, Eq. (18), with q 0 = 5, µ = 3 for Ablowitz-Ladik, layer peeling (on this method's runtime refer to the text), Crank-Nicolson and spectral methods vs. the number of temporal samples.