Fast Inverse Nonlinear Fourier Transform

TL;DR

This paper introduces two efficient algorithms for the inverse nonlinear Fourier transform based on the non-Hermitian Zakharov-Shabat problem, significantly improving computational speed and accuracy for practical signal processing applications.

Contribution

The paper presents novel fast inverse NFT algorithms with $O(KN+N ext{log}^2N)$ complexity, combining a new layer-peeling scheme and Darboux transformation for practical use.

Findings

Algorithms achieve $O(N ext{log}^2N)$ complexity.

Convergence rate of $O(N^{-2})$ demonstrated.

Numerical tests confirm high performance.

Abstract

This paper considers the non-Hermitian Zakharov-Shabat (ZS) scattering problem which forms the basis for defining the SU$(2)$-nonlinear Fourier transform (NFT). The theoretical underpinnings of this generalization of the conventional Fourier transform is quite well established in the Ablowitz-Kaup-Newell-Segur (AKNS) formalism; however, efficient numerical algorithms that could be employed in practical applications are still unavailable. In this paper, we present two fast inverse NFT algorithms with $O(KN+N\log^2N)$ complexity and a convergence rate of $O(N^{-2})$ where $N$ is the number of samples of the signal and $K$ is the number of eigenvalues. These algorithms are realized using a new fast layer-peeling (LP) scheme ($O(N\log^2N)$) together with a new fast Darboux transformation (FDT) algorithm ($O(KN+N\log^2N)$) previously developed by the author. The proposed fast inverse NFT algorithm proceeds in two steps: The first step involves computing the radiative part of the potential using the fast LP scheme for which the input is synthesized under the assumption that the radiative potential is nonlinearly bandlimited, i.e., the continuous spectrum has a compact support and the discrete spectrum is empty. The second step involves addition of bound states using the FDT algorithm. Finally, the performance of these algorithms is demonstrated through exhaustive numerical tests.

Tables

Table 1. Table 1: Implicit Adams Method

Method	$𝜷$	Order of Convergence
IA1	$(\frac{1}{2},\frac{1}{2})$	$2$
IA2	$(-\frac{1}{12},\frac{8}{12},\frac{5}{12})$	$3$
IA3	$(\frac{1}{24},-\frac{5}{24},\frac{19}{24},\frac{9}{24})$	$4$

Equations

$${\mathbf{v}}_t=-i\zeta {\sigma }_3\mathbf{v}+U\mathbf{v},$$

$${\mathfrak{S}}_K=\{({\zeta }_k,b_k)\in {\mathbb{C}}^2|\mathrm{Im}{\zeta }_k>0,k=1,2,\dots ,K\}.$$

$$a_S(\zeta )=\prod _{k=1}^K\left(\frac{\zeta -{\zeta }_k}{\zeta -{\zeta }_k^{\ast }}\right),$$

$$i{\partial }_{x}q={\partial }_t^{2}q+2|q{|}^{2}q,(t,x)\in \mathbb{R}\times {\mathbb{R}}_{+},$$

$$M_{n+1}(z^2)=\frac{z^{-1}}{{\Theta }_{n+1}}\left(\begin{array}{cc}1+z^2{Q}_{n+1}{R}_n& z^2{Q}_{n+1}+Q_n\\ R_{n+1}+z^2{R}_n& R_{n+1}{Q}_n+z^2\end{array}\right).$$

$$\begin{array}{rl}{\mathbf{v}}_{n+1}& =\frac{1}{{\Theta }_{n+1}}\left(\begin{array}{cc}1& Q_{n+1}\\ R_{n+1}& 1\end{array}\right)\left(\begin{array}{cc}{z}^{-1}& 0\\ 0& z\end{array}\right)\left(\begin{array}{cc}1& Q_n\\ R_n& 1\end{array}\right){\mathbf{v}}_n,\end{array}$$

$${\mathbf{P}}_n(z)=\left(\begin{array}{c}{P}_1^{(n)}(z)\\ P_2^{(n)}(z)\end{array}\right)=\sum _{k=0}^n{\mathbf{P}}_k^{(n)}{z}^k=\sum _{k=0}^n\left(\begin{array}{c}{P}_{1,k}^{(n)}\\ P_{2,k}^{(n)}\end{array}\right)z^k.$$

$${\mathbf{P}}_{n+1}(z^2)=M_{n+1}(z^2){\mathbf{P}}_n(z^2).$$

$$P^{(n+1)}_{1,0}=\Theta^{-1}_{n+1}\prod_{k=1}^{n}\biggl{(}\frac{1+Q_{k}R_{k}}{1-Q_{k}R_{k}}\biggl{)}=\Theta^{-1}_{n+1}\prod_{k=1}^{n}\biggl{(}\frac{2-\Theta_{k}}{\Theta_{k}}\biggl{)},$$

$$R_{n+1}=\frac{P_{2,0}^{(n+1)}}{P_{1,0}^{(n+1)}},R_n=\frac{\chi }{1+\sqrt{1+|\chi {|}^2}},$$

$$\chi =\frac{[P_{2,1}^{(n+1)}-R_{n+1}{P}_{1,1}^{(n+1)}]}{[P_{1,0}^{(n+1)}-Q_{n+1}{P}_{2,0}^{(n+1)}]}.$$

$$\rho (\xi )=\sum _{k\in \mathbb{Z}}{\rho }_k{e}^{\frac{ik\pi \xi }{\Lambda }}.$$

$$\breve{\rho }(\xi )=\rho (\xi )z^{2n}=\sum _{k=0}^{\infty }{\breve{\rho }}_k{z}^{2k}+z^{2n}{\mathcal{R}}_n(z^2).$$

$$\breve{p}(\tau )=F^{-1}[\rho ](\tau )=\frac{1}{2\pi }{\int }_{-\Lambda }^{\Lambda }\breve{\rho }(\xi )e^{-i\xi \tau }d\xi .$$

$$P_1^{(N)}(z^2)=1,P_2^{(N)}(z^2)=\sum _{k=0}^{N-1}{\breve{\rho }}_k{z}^{2k}.$$

$$\Delta \xi =\frac{\pi }{2Mh}.$$

$$2hk\breve{p}(2hk)$$

$$=\frac{1}{2M}\sum _{j=-M}^{M-1}\breve{\rho }({\xi }_j)e^{-i\frac{2\pi jk}{2M}},$$

$${\mathcal{E}}_{+}(T)={\int }_T^{\infty }|q(t){|}^{2}dt\le \frac{2{\mathcal{I}}_2^2(T)}{[1-{\mathcal{I}}_1^2(T)]},$$

$${\mathcal{I}}_m(T)={\left[{\int }_{2T}^{\infty }|p(-\tau ){|}^{m}d\tau \right]}^{1/m}$$

$$\frac{2{\mathcal{I}}_2^2(T)}{[1-{\mathcal{I}}_1^2(T)]}\le ϵ,$$

$$\rho (\xi )=\sum _{k\in \mathbb{Z}}{\rho }_k{z}^{2k},z=e^{i\xi h}.$$

$$\tilde{g}(z^2)=\{\begin{array}{ll}g(z^2)& |z|<1,\\ -g^{\ast }(1/z^{\ast 2})& |z|>1,\end{array}$$

$${\tilde{g}}^{(-)}(z^2)-{\tilde{g}}^{(+)}(z^2)=\log \left[\frac{|\rho (z^2){|}^2}{1+|\rho (z^2){|}^2}\right],|z|=1,$$

$$f(z^2)=\sum _{k\in \mathbb{Z}}{f}_k{z}^{2k},|z|=1.$$

$$\tilde{g}(z^2)=\frac{1}{2\pi i}{\oint }_{|w|=1}\frac{f(w)}{z^2-w}dw.$$

$$g(z^2)=\sum _{k\in {\mathbb{Z}}_{+}\cup \{0\}}{f}_k{z}^{2k},|z|<1.$$

$$\breve{b}(\xi )=b(\xi )z^{2n}=\left[\sum _{k=0}^{\infty }{\breve{\rho }}_k{z}^{2k}+z^{2n}{\mathcal{R}}_n(z^2)\right]\exp [g(z^2)].$$

$${\breve{b}}_N(z^2)={\left\{a_N(z^2)\sum _{k=0}^{N-1}{\breve{\rho }}_k{z}^{2k}\right\}}_N=\sum _{k=0}^{N-1}{\breve{b}}_k{z}^{2k},$$

$$P_1^{(N)}(z^2)=\sum _{k=0}^{N-1}{a}_k{z}^{2k},P_2^{(N)}(z^2)=\sum _{k=0}^{N-1}{\breve{b}}_k{z}^{2k}.$$

$$e_{\text{rel.}}=\Vert q^{(\text{num.})}-q{\Vert }_{{\mathsf{L}}^2}/\Vert q{\Vert }_{{\mathsf{L}}^2},$$

Full text

Notations

The set of non-zero positive real numbers ($\mathbb{R}$) is denoted by $\mathbb{R}_{+}$. Non-zero positive (negative) integers are denoted by $\mathbb{Z}_{+}$ ($\mathbb{Z}_{-}$). For any complex number $\zeta$, $\operatorname{Re}(\zeta)$ and $\operatorname{Im}(\zeta)$ refer to the real and the imaginary parts of $\zeta$, respectively. Its complex conjugate is denoted by $\zeta^{*}$. The upper-half (lower-half) of complex plane, $\mathbb{C}$, is denoted by $\mathbb{C}_{+}$ ($\mathbb{C}_{-}$).

I Introduction

The nonlinear Fourier (NF) spectrum offers a novel way of encoding information in optical pulses where the nonlinear effects are adequately taken into account as opposed to being treated as a source of distortion. This idea has its origin in the work of Hasegawa and Nyu Hasegawa and Nyu (1993) who were the first to propose the use of discrete eigenvalues of the NF spectrum for encoding information. Recent advances in coherent optical communication have made it possible to reconsider this old idea with some extensions and improvements. Extension of this scheme consists in using additional degrees of freedom offered by the NF spectrum such as the norming constants and the continuous spectrum. For an overview of the recent progress in theoretical as well as experimental aspects of various optical communication methodologies that are based on the nonlinear Fourier transform (NFT), we refer the reader to the review article Turitsyn et al. (2017) and the references therein.

In order to realize any NFT-based modulation methodology, it is imperative to have a suitable low-complexity NFT algorithm which forms the primary motivation behind this work. The central idea is to use a fast version of the well-known layer-peeling (LP) algorithm within the framework of an appropriate discretization scheme applied to the Zakharov-Shabat (ZS) problem. This approach has been characterized as the differential approach by Bruckstein et al. Bruckstein et al. (1985); Bruckstein and Kailath (1987) where fast realizations of the LP algorithm which achieves a complexity of $\mathop{\mathscr{O}}\left(N\log^{2}N\right)$ for $N$ samples of the reflection data are also discussed. However, the earliest work on fast LP is that of McClary McClary (1983) which appeared in the geophysics literature. More recently, this method has been adopted by Brenne and Skaar Brenne and Skaar (2003) in the design of grating-assisted codirectional couplers. However, this paper reports a complexity of $\mathop{\mathscr{O}}\left(N^{2}\right)$111In this paper, we do not consider the method of discretization presented in Brenne and Skaar (2003); however, let us briefly mention that on account of the piecewise constant assumption used in this work for the scattering potential, the order of convergence gets artificially restricted to $\mathop{\mathscr{O}}\left(N^{-1}\right)$. This problem has been remedied in Vaibhav (2017) where this discretization scheme is termed as the split-Magnus method.. It is interesting to note that, at the heart of it, all of the aforementioned versions of LP are similar; however, the manner in which the discrete system is obtained seem to vary. In this work, we consider the discrete system obtained as a result of applying (exponential) trapezoidal rule to the ZS problem as discussed in Vaibhav (2017).

The next important idea is to recognize that the Darboux transformation (DT) provides a promising route to the most general inverse NFT algorithm. A fast version of DT (referred to as FDT) is developed in Vaibhav (2017) which is based on the pioneering work of Lubich on convolution quadrature Lubich (1994) and a fast LP algorithm. The schematic of the fast inverse NFT is shown in Fig. 1 where we note that FDT is capable of taking a seed potential $q_{R}(t)$ and augmenting it by introducing the bound states corresponding to $\mathfrak{S}_{K}$ (the discrete spectrum to be introduced in Sec. II and $K$ is the number of bound states or eigenvalues). If $q_{R}(t)$ is the radiative part of $q(t)$, i.e., it is generated from NF spectrum which has an empty discrete spectrum and $\rho_{R}(\xi)$ as the reflection coefficient, then $q(t)$ is the full inverse of the NF spectrum characterized by $\mathfrak{S}_{K}$ and $\rho(\xi)$. The preliminary results of this approach were reported in Vaibhav and Wahls (2017). In this paper, we describe two fast inverse NFT algorithms that exhibit a complexity of $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ and a rate of convergence of $\mathop{\mathscr{O}}\left(N^{-2}\right)$ where $N$ is the number of samples and $K$ is the number of eigenvalues (or bound states).

Finally, we note that the LP algorithm (irrespective of the underlying discrete system) has the reputation of being ill-conditioned or unstable in the presence of noise Bruckstein et al. (1986); Skaar and Feced (2002) in the reflection coefficient. For optical communication, this observation is important but not critical as the reflection coefficient is known exactly at the stage of encoding of information at the transmitter end. A more relevant question here, therefore, is the stability of the algorithm in the presence of round-off errors. We provide exhaustive numerical tests in order to understand the ill-conditioning effects; however, no theoretical results for stability are provided.

This paper is organized as follows: Sec. II discusses the basic theory of scattering. Sec. III introduces the discrete framework for forward/inverse scattering, which admits of the layer-peeling property. This section also introduces a recipe for computing a class of signals dubbed as the nonlinearly bandlimited signals. Finally, the inverse NFT is described in Sec. III.3 and the numerical results are presented in Sec. IV. Sec. V concludes this paper.

II The AKNS System

The NFT of any complex-valued signal $q(t)$ is introduced via the associated Zakharov-Shabat scattering problem Zakharov and Shabat (1972) which can be stated as follows: Let $\zeta\in\mathbb{R}$ and $\bm{v}=(v_{1},v_{2})^{\intercal}\in\mathbb{C}^{2}$, then

[TABLE]

where $\sigma_{3}=\operatorname{diag}(1,-1)$, and, the matrix elements of $U$ are $U_{11}=U_{22}=0$ and $U_{12}=q(t)=-U_{21}^{*}=-r^{*}(t)$. Here, $q(t)$ is identified as the scattering potential. The solution of the scattering problem (1), henceforth referred to as the ZS problem, consists in finding the so called scattering coefficients which are defined through special solutions of (1) known as the Jost solutions. The Jost solutions of the first kind, denoted by $\bm{\psi}(t;\zeta)$, has the asymptotic behavior $\bm{\psi}(t;\zeta)e^{-i\zeta t}\rightarrow(0,1)^{\intercal}$ as $t\rightarrow\infty$. The Jost solutions of the second kind, denoted by $\bm{\phi}(t,\zeta)$, has the asymptotic behavior $\bm{\phi}(t;\zeta)e^{i\zeta t}\rightarrow(1,0)^{\intercal}$ as $t\rightarrow-\infty$. The so-called scattering coefficients, $a(\zeta)$ and $b(\zeta)$, are obtained from the asymptotic behavior $\bm{\phi}(t;\zeta)\rightarrow(a(\zeta)e^{-i\zeta t},b(\zeta)e^{i\zeta t})^{\intercal}$ as $t\rightarrow\infty$. The process of computing these scattering coefficients will be referred to as forward scattering.

In general, the nonlinear Fourier spectrum for the potential $q(t)$ comprises a discrete and a continuous spectrum. The discrete spectrum consists of the so called eigenvalues $\zeta_{k}\in\mathbb{C}_{+}$, such that $a(\zeta_{k})=0$, and, the norming constants $b_{k}$ such that $\bm{\phi}(t;\zeta_{k})=b_{k}\bm{\psi}(t;\zeta_{k})$. Note that $(\zeta_{k},\,b_{k})$ describes a bound state or a solitonic state associated with the potential. For convenience, let the discrete spectrum be denoted by the set

[TABLE]

The continuous spectrum, also referred to as the reflection coefficient, is defined by $\rho(\xi)={b(\xi)}/{a(\xi)}$ for $\xi\in\mathbb{R}$. In preparation for the discussion in the following sections, let us define

[TABLE]

and $\rho_{R}(\xi)=a_{S}(\xi)\rho(\xi)$. The reflection coefficient $\rho_{R}(\xi)$ now corresponds to a purely radiative potential.

Next, let us note that the class of integrable nonlinear evolution problems that can be treated by the methods proposed in this article are those described by the Ablowitz-Kaup-Newell-Segur formalism Ablowitz et al. (1974); Ablowitz and Segur (1981). In optical fiber communication, the propagation of optical field in a loss-less single mode fiber under Kerr-type focusing nonlinearity is governed by the nonlinear Schrödinger equation (NSE) Kodama and Hasegawa (1987); Agrawal (2013) which can be cast into the following standard form

[TABLE]

where $q(t,x)$ is a complex valued function associated with the slowly varying envelope of the electric field, $t$ is the retarded time and $x$ is position along the fiber. If the potential evolves according to (4), then, the scattering data evolves as: $b_{k}(x)=b_{k}e^{-4i\zeta_{k}^{2}x}$ and $\rho(\xi,x)=\rho(\xi)e^{-4i\xi^{2}x}$ ($a(\zeta)$ and, consequently, $\zeta_{k}$ do not evolve). In the rest of the paper, we suppress the dependence on $x$ for the sake brevity.

III Discrete Inverse Scattering

In order to discuss the discretization scheme, we take an equispaced grid defined by $t_{n}=T_{1}+nh,\,\,n=0,1,\ldots,N,$ with $t_{N}=T_{2}$ where $h$ is the grid spacing. Define $\ell_{-},\ell_{+}\in\mathbb{R}$ such that $h\ell_{-}=-T_{1}$, $h\ell_{+}=T_{2}$. Further, let us define $z=e^{i\zeta h}$. For the potential functions supported in $[T_{1},T_{2}]$, we set $Q_{n}=2hq(t_{n})$, $R_{n}=2hr(t_{n})$. In the following, we summarize the discrete framework reported in Vaibhav (2017) which is based on the trapezoidal rule of integration. Setting $\Theta_{n}=1-Q_{n}R_{n}$, the recurrence relation for the Jost solution reads as $\bm{v}_{n+1}=z^{-1}M_{n+1}(z^{2})\bm{v}_{n}$, which is referred to as the discrete scattering problem. Here $M_{n+1}(z^{2})$ is known as the transfer matrix which is given by

[TABLE]

Note that the transfer matrix approach introduced above is analogous to that used to solve wave-propagation problems in dielectric layered-media (Born and Wolf, 1999, Chap. 1). In particular, from the factorization

[TABLE]

it can be inferred that the continuous system in (1) is approximated by two instantaneous scatterers with “free-space” propagation between them in each of the layers as shown in Fig. 2. The error analysis of the discrete system presented above is carried out in Vaibhav (2017) where it is shown that the global order of convergence is $\mathop{\mathscr{O}}\left(h^{2}\right)$ for fixed $\zeta$.

In order to express the discrete approximation to the Jost solutions, let us define the vector-valued polynomial

[TABLE]

The Jost solution $\bm{\phi}$ can be written in the form $\bm{\phi}_{n}=z^{\ell_{-}}z^{-n}\bm{P}_{n}(z^{2})$ with the initial condition given by $\bm{\phi}_{0}=z^{\ell_{-}}(1,0)^{\intercal}$ that translate into $\bm{P}_{0}=(1,0)^{\intercal}$. The recurrence relation for $\bm{P}_{n}(z^{2})$ takes the form

[TABLE]

The discrete system discussed above facilitated the development of a fast forward scattering algorithm in Vaibhav (2017). This relied on the fact that the transfer matrices have polynomial entries–a form that is amenable to FFT-based fast polynomial arithmetic Henrici (1993).

In the following sections, we provide details of the fast inverse NFT algorithm by first developing the methods needed for inversion of the continuous spectrum to compute what can be viewed as a purely radiative potential. The general version of the inverse NFT is then developed using the FDT algorithm presented in Vaibhav (2017).

III.1 The layer-peeling algorithm

Borrowing the terminology from the theory of layered dielectric media (Born and Wolf, 1999, Chap. 1), let the interval $[t_{n},t_{n+1}]$ correspond to the $(n+1)$-th layer which is completely characterized by the transfer matrix $M_{n+1}(z^{2})$ (see Fig. 2). The discrete forward scattering consists in “accumulating” all the layers to form $\bm{P}_{N}(z^{2})$. The problem of recovering the discrete samples of the scattering potential from the discrete scattering coefficients or $\bm{P}_{N}(z^{2})$ is referred to as the discrete inverse scattering which is facilitated by the so-called layer-peeling (LP) algorithm. Starting from the recurrence relation (7), one LP step consists in using $\bm{P}_{n+1}(z^{2})$ to retrieve the samples of the potential needed to compute the transfer matrix $\widetilde{M}_{n+1}(z^{2})=z^{-2}[M_{n+1}(z^{2})]^{-1}$ so that the entire step can be repeated with $\bm{P}_{n}(z^{2})$ until all the samples of the potential are recovered. The mathematical details of this algorithm can be found in Vaibhav (2017). For the sake of reader’s convenience, some of the main results are summarized below.

Assume $Q_{0}=0$. Then the recurrence relation (7) yields

[TABLE]

and $\bm{P}^{(n+1)}_{n+1}=0$. The last relationship follows from the assumption $Q_{0}=0$. For sufficiently small $h$, it is reasonable to assume that $1+Q_{n}R_{n}>0$ so that $P^{(n)}_{1,0}>0$ (it also implies that $|Q_{n}|=|R_{n}|<1$). The layer-peeling step consists in computing the samples of the potential, $R_{n+1}$ and $R_{n}$ (with $Q_{n+1}=-R^{*}_{n+1}$ and $Q_{n}=-R^{*}_{n}$) as follows:

[TABLE]

where

[TABLE]

Note that $P^{(n+1)}_{1,0}\neq 0$ and ${P^{(n+1)}_{1,0}-Q_{n+1}P^{(n+1)}_{2,0}}\neq 0$. As evident from (5), the transfer matrix, $M_{n+1}(z^{2})$, connecting $\bm{P}_{n}(z^{2})$ and $\bm{P}_{n+1}(z^{2})$ is completely determined by these relations.

If the steps of the LP algorithm are carried out sequentially, one ends up with a complexity of $\mathop{\mathscr{O}}\left(N^{2}\right)$. It turns out that a fast implementation of this LP algorithm does exist Vaibhav (2017), which has a complexity of $\mathop{\mathscr{O}}\left(N\log^{2}N\right)$ for the discrete system considered in this article. In the following sections, we describe how to synthesize the input for the LP algorithm in order to compute the radiative part of the scattering potential.

III.2 Nonlinearly bandlimited signals

A signal is said to be nonlinearly bandlimited if it has an empty discrete spectrum and a reflection coefficient $\rho(\xi)$ that is compactly supported in $\mathbb{R}$. This is a direct generalization of the notion of bandlimited signals for conventional Fourier transform. However, nonlinearly bandlimited signals are not bandlimited, in general. Let us consider the reflection coefficient $\rho(\xi)$ as input. Let the support of $\rho(\xi)$ be contained in $[-\Lambda,\Lambda]$ so that its Fourier series representation is

[TABLE]

If $|\rho_{k}|$ is significant only for $k\geq-n$ ($n\in\mathbb{Z}_{+}$), then $\rho(\xi)=\sum_{k=-n}^{\infty}\rho_{k}z^{2k}+\mathcal{R}_{n}(z^{2})$, where $z=\exp(i\pi\xi/2\Lambda)$ and $\mathcal{R}_{n}$ denotes the remainder terms. Putting $h=\pi/2\Lambda$ and $T_{2}=nh\equiv h\ell_{+}$, we have $\exp(2i\xi T_{2})=z^{2n}$ so that

[TABLE]

Now, it follows that $\breve{\rho}_{k}=2h\breve{p}(2hk)$ where

[TABLE]

Let $2\Lambda_{0}$ be the fundamental period and $\Lambda=m\Lambda_{0}$, where $m\in\mathbb{Z}_{+}$; then, $h=\pi/2m\Lambda_{0}\equiv h_{0}/m$; therefore, $h\leq h_{0}$. Now, if we ignore the remainder term and truncate the series after $N$ terms in (11), the input to the fast LP algorithm can be

[TABLE]

This accomplishes the inversion of the reflection coefficient which is assumed to be compactly supported. Let $\xi_{j}=j\Delta\xi$ for $j=-M,\ldots,M-1$, where

[TABLE]

Then the coefficients $\breve{\rho}_{k}$ can be estimated using the Fourier sum

[TABLE]

for $k=0,1,\ldots,N$. The quantity $M$ is chosen to be some multiple of $N$, say, $M=n_{\text{os}}\times N$ where $n_{\text{os}}$ is referred to as the oversampling factor. Therefore, the overall complexity of synthesizing the input for the LP algorithm works out to be $\mathop{\mathscr{O}}\left(N\log N\right)$.

Before we conclude this discussion, let us consider the problem of estimation of $T_{2}$. It is of interest to determine a $T_{2}$ such that the energy in the tail of the scattering potential, which is to be neglected, is below a certain threshold, say, $\epsilon$. Fortunately, there is an interesting result due to Epstein Epstein (2004) that allows us to do exactly that. From the theory of Gelfand-Levitan-Marchenko equations, it can be shown that there exists a time $T$ such that

[TABLE]

assuming $\mathcal{I}_{1}(T)<1$ where

[TABLE]

for $m=1,2$ (see Appendix B for a proof which, in essence, is contained in the work of Epstein Epstein (2004)). Let $T=T(\epsilon)$ be such that

[TABLE]

then $\mathcal{E}_{+}(T)\leq\epsilon$. Consequently, it suffices to choose $T_{2}\geq T(\epsilon)$.

III.2.1 Alternative approach

It is possible to compute the polynomial approximation to the scattering coefficients $a(\xi)$ and $b(\xi)$ using $\rho(\xi)$, which can be then used to synthesize the input to the fast LP algorithm. There is no apparent benefit of this approach compared to the method described above; however, we describe it for the sake of completeness. The first step consists of constructing a polynomial approximation for $a(\zeta)$ in $|z|<1$ where $z=e^{i\zeta h}$ (under the assumption that no bound states are present). To this end, let

[TABLE]

With a slight abuse of notation, let us denote this expansion as $\rho(z^{2})$. Let us note that in this case, $a(\xi)$ is not analytic in $\mathbb{R}$ which means that it is also not analytic on the unit circle $|z|=1$. Here, the relation Ablowitz et al. (1974); Ablowitz and Segur (1981) $|a(\xi)|^{2}+|b(\xi)|^{2}=1$ allows us to set up a Riemann-Hilbert (RH) problem for a sectionally analytic function

[TABLE]

such that the jump condition is given by

[TABLE]

where $\tilde{g}^{(-)}(z^{2})$ and $\tilde{g}^{(+)}(z^{2})$ denotes the boundary values when approaching the unit circle from $|z|<1$ and $|z|>1$, respectively. Let the jump function on the RHS of (18) be denoted by $f(z^{2})$ which can be expanded as a Fourier series

[TABLE]

Now, the solution to the RH problem can be stated using the Cauchy integral (Henrici, 1993, Chap. 14)

[TABLE]

The function $g(z^{2})$ analytic in $|z|<1$ then works out to be

[TABLE]

Finally, $a_{N}(z^{2})=\{\exp[g(z^{2})]\}_{N}$ with $z=e^{i\zeta h}$ where $\{\cdot\}_{N}$ denotes truncation after $N$ terms. The implementation of the procedure laid out above can be carried out using the FFT algorithm, which involves computation of the coefficients $f_{k}$ and the exponentiation in the last step (Henrici, 1993, Chap. 13). Note that, in the computation of $g(z^{2})$, we discarded half of the coefficients; therefore, in the numerical implementation it is necessary to work with at least $2N$ number of samples of $f(z^{2})$ in order to obtain $a_{N}(z^{2})$ which is a polynomial of degree $N-1$.

The next step is to compute the polynomial approximation for $\breve{b}(\xi)$. To this end, consider

[TABLE]

In the following, we will again discard the remainder term. The polynomial approximation for $\breve{b}(\xi)$ reads as

[TABLE]

Now, the input to the fast LP algorithm works out to be

[TABLE]

III.3 Fast inverse NFT

In the previous sections, we restricted ourselves to the case of empty discrete spectrum. In this section, we describe how a fast inverse NFT algorithm can be developed for the general NF spectrum using either the Classical DT (CDT) or the FDT algorithm reported in Vaibhav (2017). Given a reflection coefficient $\rho(\xi),\,\xi\in\mathbb{R},$ and the discrete spectrum $\mathfrak{S}_{K}$, define $a_{S}(\xi)$ as in (3) and $\rho_{R}(\xi)=a_{S}(\xi)\rho(\xi)$. Let $q(t)$ denote the scattering potential corresponding to the aforementioned NF spectrum.

Now, as illustrated in Fig. 1, the inverse NFT can be carried out in the following two steps:

• I.

Generate the signal $q_{R}(t)$ corresponding to the reflection coefficient $\rho_{R}(\xi)$ using the method described in Sec. III.2. This amounts to computing the purely radiative part of the complete potential $q(t)$. The complexity of this step is $\mathop{\mathscr{O}}\left(N\log^{2}N\right)$ if the number of nodes used for the FFT operation involved there is given by $M=n_{\text{os}}N$ where $n_{\text{os}}\ll N$. Here, $n_{\text{os}}$ can be identified as the oversampling factor (typically $\leq 8$).

• II.

Use the signal $q_{R}(t)$ as the seed potential and add bound states described by $\mathfrak{S}_{K}$ using the CDT or the FDT algorithm to obtain $q(t)$. The complexity of this step is $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ when FDT is employed while $\mathop{\mathscr{O}}\left(K^{2}N\right)$ when CDT is employed. Here we also consider the partial-fraction (PF) variant of the FDT algorithm (labeled as FDT-PF), which is shown to offer a small increase in speed Vaibhav (2017).

Finally, let us note that the overall complexity of the inverse NFT is given by $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ when FDT is used and $\mathop{\mathscr{O}}\left(N(K^{2}+\log^{2}N)\right)$ when CDT is used.

IV Numerical Experiments

Let $q^{(\text{num}.)}$ denote the numerically computed potential for a given NF spectrum. If the exact potential $q$ is known, then we quantify the error as

[TABLE]

where the integrals are evaluated numerically using the trapezoidal rule. For the purpose of convergence analysis, only those examples are deemed to be admissible where closed-form solutions are available. However, on account of scarcity of such examples, an exhaustive test for universality of the algorithm cannot be carried out in this manner. To remedy this, we choose a higher-order convergent algorithm for the forward scattering problem and compute the NF spectrum of the potential generated by the fast inverse NFT. The error between the computed NF spectrum and the provided NF spectrum can serve as a good metric to measure the robustness of the algorithm.

For the higher-order scheme, we choose the (exponential) $3$-step implicit Adams method (IA3) Hairer et al. (1993) which has an order of convergence $4$, i.e., $\mathop{\mathscr{O}}\left(N^{-4}\right)$ (see Appendix A for details). Fortunately, this method can also be made fast by the use of FFT-based polynomial arithmetic which allows us to test for large number of samples ($N\in\{2^{10},2^{11},\ldots,2^{20}\}$). Note that this procedure by no means qualifies as the test for total numerical error on account of the fact that the error metric is not the true numerical error. Therefore, the results in this case must be interpreted with caution. Further, for the sake of comparison, we also consider the Töplitz inner bordering (TIB) algorithm for inverse scattering (Belai et al. Belai et al. (2007)) whenever the discrete spectrum is empty. We use the second order convergent version of this algorithm which has also been reported in Frumin et al. Frumin et al. (2015). The latter paper provides an improved understanding of the original algorithm presented in Belai et al. (2007); therefore, we choose to refer to Frumin et al. (2015) in this article whenever we mention the TIB algorithm.

Finally, let us emphasize that the primary objective behind the numerical tests in this section is to verify the trends expected from the theory. The actual values of any defined performance metric observed in the results are merely representative of what can be achieved222The total run-time, for instance, may differ on different computing machines; therefore, we would only be interested in trends as far as the complexity analysis of the algorithms are concerned., and, admittedly, better results can be obtain by appropriately tuning the parameters used in the test. For instance, a good choice of the computational domain helps to maintain a smaller step-size in the numerical discretization and, hence, lowers the numerical error.

IV.1 Secant-hyperbolic potential

Here, we would like to devise tests to confirm the order of convergence and the complexity of computations for the algorithms proposed thus far. To this end, we choose the secant-hyperbolic potential given by $q(t)=(A_{R}+K)\operatorname{sech}t$, which is treated exactly in Satsuma and Yajima (1974). Here $A_{R}\in[0,0.5)$ and $K$ is a positive integer. The discrete spectrum can be stated as

[TABLE]

and the continuous spectrum is given by $\rho=\rho_{R}/a_{S}$ where $a_{S}(\xi)$ is defined by (3) and

[TABLE]

with $b(\xi)=-\sin[(A_{R}+K)\pi]\operatorname{sech}(\pi\xi)$. This test consist in studying the behavior of the fast INFTs for different number of samples ($N$) as well as eigenvalues ($K$). We set $A_{R}=0.4$. The scattering potential is scaled by $\kappa=2(\sum_{k=1}^{K}\operatorname{Im}\zeta_{k})^{1/2}$ and $[-T,\,T]$, $T=30\kappa/\min_{k}(\operatorname{Im}\zeta_{k})$, is taken as the computational domain and we set $N_{\text{th}}=N/8$ for FDT-PF as in Vaibhav (2017).

Let us first consider the case $K=0$ so that $\rho=\rho_{R}$ (setting the convention that $a_{S}=1$ when $K=0$). Note that on account of the exponential decay of $\rho$, it can be assumed to be effectively supported in a bounded domain. Besides the knowledge of the true potential allows us to provide a good estimate of the computational domain. Set $T=\log(2A_{R}/\epsilon)\approx 30$ for $\epsilon=10^{-12}$, then $[-T,\,T]$ can be taken as the computational domain333For the ZS problem, let us note that the error in the initial condition at the left-boundary can be kept below $\epsilon>0$, if $\|q\chi_{(-\infty,T_{1}]}\|_{\mathsf{L}^{1}}\leq\sinh^{-1}\epsilon$ Vaibhav (2017).. The result for $A_{R}=0.4$ is plotted in Fig. 3 which shows that the performance of LP is comparable to that of TIB. Further, each of these algorithms exhibit a second order of convergence (i.e., error vanishing as $\mathop{\mathscr{O}}\left(N^{-2}\right)$). The run-time behavior in Fig. 3 shows that LP-based INFTs have a poly-log complexity per sample as opposed to the $\mathop{\mathscr{O}}\left(N\right)$ complexity per sample exhibited by TIB.

For $K>0$, the results are plotted in Fig. 4 which reveal that the fast INFTs based on FDT (labeled as ‘INFT-B’) and FDT-PF (labeled as ‘INFT-C’) are superior to that based on CDT (labeled as ‘INFT-A’) which becomes unstable with increasing number of eigenvalues. The latter, however, can be useful for a small number eigenvalues. The figure also confirms the second order of convergence of INFT-B/-C which is consistent with the underlying one-step method, namely, the trapezoidal rule. For small number of eigenvalues, INFT-A also exhibits a second order of convergence. Finally, let us observe that, for fixed $N$, INFT-A has a complexity of $\mathop{\mathscr{O}}\left(K^{2}\right)$ and that for INFT-B/-C is $\mathop{\mathscr{O}}\left(K\right)$. While these trends can be confirmed from Fig. 5, let us mention that, with an improved implementation, INFT-B/-C can be made even more competitive to INFT-A in complexity.

IV.2 Nonlinearly bandlimited signals

Let us consider a soliton-free signal whose continuous spectrum is given by

[TABLE]

where $\Xi=|\tau_{s}\xi|-(1-\beta)$ with $\beta\in[0,1]$, and, $A_{\text{rc}}$ and $\tau_{s}$ are positive constants. The nonlinear impulse response (NIR) $h_{\text{rc}}(\tau)$ can be worked out exactly; however, we do not use this information for constructing the input to the fast LP algorithm. Note that $H_{\text{rc}}(\xi)$ and $h_{\text{rc}}(\tau)$ describe the well-known raised-cosine filter in the frequency-domain and the time-domain, respectively.

In order to estimate the computational domain, we use Epstein’s result discussed in Sec. III.2 which consists in finding a time $T$ such that

[TABLE]

assuming $\mathcal{I}_{1}(T)<1$ where

[TABLE]

for $m=1,2$. A crude estimate for $T$ such that $\mathcal{I}^{2}_{2}(T)=\epsilon$ is given by

[TABLE]

which uses the asymptotic form of $h_{\text{rc}}(\tau)$. If $\epsilon\ll 1$, we may assume that the potential is effectively supported444The Epstein’s theorem provides an estimate for the right boundary if the right NIR is used; therefore, strictly speaking, the computational domain must be of the form $(-\infty,T(\epsilon)]$. in $[-T(\epsilon),T(\epsilon)]$ where we set $\epsilon=10^{-9}$. Also, let $\beta=0.5$ and $\tau_{s}=1$ in the following.

For this example, we devise two kinds of tests. For the first kind of tests, we disregard any modulation scheme and carry out the inverse NFT for varying number of samples ($N$) for each of the values of $A_{\text{rc}}\in\{10,\ldots,50\}$. In the second kind of tests, we considers the quadrature-phase-shift-keyed (QPSK) modulation scheme which is described later. Let $\Omega_{h}=[-\pi/2h,\pi/2h]$, then the error is quantified by

[TABLE]

where the integrals are computed from $N$ equispaced samples in $\Omega_{h}$ using the trapezoidal rule. As stated in the beginning, the quantity $\rho^{(\text{num.})}$ is computed using the (exponential) IA3.

The results of the first kind of tests are shown in Fig. 7 where a comparison is made between LP and TIB555The complexity of TIB becomes prohibitive for increasing $N$, therefore, we restrict ourselves to $N\leq 2^{18}$.. From the plots in the top row of Fig. 7, the second order of convergence is readily confirmed for both of these algorithms with LP performing somewhat better than TIB. The plateauing of the error in these plots can be attributed to accumulating numerical errors in the inverse NFT algorithm as well as the implicit Adams method. The behavior of the error with respect to $A_{\text{rc}}$ is shown in the bottom row of Fig. 7 where LP shows better performance than TIB.

Now, for the second kind of tests, we consider the QPSK modulation of the continuous spectrum as follows

[TABLE]

where the index set is $J=\{-N_{\text{sym}}/2,\ldots,N_{\text{sym}}/2-1\}$ and $s_{n}\in\{\pm 1,\pm i\}$ with $N_{\text{sym}}>0$ being an even integer. The estimate for the right boundary works out to be

[TABLE]

however, an estimate for the for the left boundary is not available in a closed form. Here, we take a heuristic approach by setting $T_{1}=-W\times T_{2}$ where $W$ is chosen by trial and error. The scale factor $A_{\text{rc}}$ is chosen such that $A_{\text{eff.}}=10$ where

[TABLE]

It is important to observe here that the signal generated from (32) is highly asymmetric with poor decay behavior as $t\rightarrow-\infty$ (see Fig. 6). The higher values of the quantities $N_{\text{sym}}$ and $A_{\text{eff.}}$, both, tend to worsen this phenomenon. Therefore, this example turns out to be very challenging for the numerical algorithm. In Fig. 8, we provide results of numerical experiments conducted with $N_{\text{sym}}\in\{4,8,\ldots,256\}$ number of symbols where $W=5\log_{2}N_{\text{sym}}$ is used to determine the computational domain. The accuracy of LP and TIB, both, tends to worsen with increasing number of symbols where LP performs slightly better than TIB. Based on these results it is evident that any method of pulse-shaping must take into account the relationship between the signal and its NF spectrum as opposed to directly applying conventional Fourier transform based techniques of pulse-shaping.

IV.2.1 Addition of Bound states

Here, we fix $A_{\text{rc}}=20$ and assume no modulation of the continuous spectrum. The bound states to be added are described by (26). Let us observe that the “augmented” potential has a reflection coefficient which is given by $\rho^{(\text{aug.})}=\rho/a_{S}$. Now the error can be quantified by (31). The potentials are scaled by $\kappa$ as in Sec. IV.1 and the computational domain is chosen such that $-T_{1}=T_{2}=T(\epsilon)\kappa/\min_{k}(\operatorname{Im}\zeta_{k})$.

The results for the continuous spectrum are shown in Fig. 9 where the order of convergence can be confirmed from the plots in the top row. The plots in the bottom row reveal that INFT-A, which is based on CDT, is unstable for increasing number of eigenvalues. On the other hand, the algorithms INFT-B/-C, which are based on FDT/FDT-PF, respectively, seem to perform equally well without showing any signs of instability.

For the discrete spectrum, we assume that the discrete eigenvalues are known exactly, and, then use this information to compute the norming constants using the method discussed in Vaibhav (2017). The error is quantified by

[TABLE]

where $b^{(\text{num.})}_{k}$ is the numerically computed norming constant using IA3. The results are shown in Fig. 10 where the order of convergence turns out to be $\mathop{\mathscr{O}}\left(N^{-1}\right)$ from the plots in the top row. This decrease of order of convergence can be attributed to the use of the true eigenvalues as opposed to the numerically computed one to compute the norming constants. Again, the plots in the bottom row reveal that INFT-A is unstable for increasing number of eigenvalues. On the other hand, the algorithms INFT-B/-C seem to perform equally well while showing no signs of instability.

V Conclusion

To conclude, we have presented two new fast INFT algorithms with $\mathop{\mathscr{O}}\left(KN+N\log^{2}N\right)$ complexity and a convergence rate of $\mathop{\mathscr{O}}\left(N^{-2}\right)$. These algorithms are based on the discrete framework introduced in Vaibhav (2017) for the ZS scattering problem where the well-known one-step method, namely, the trapezoidal rule is employed for the numerical discretization. Further, our algorithm depends on the fast LP and the FDT algorithm presented in Vaibhav (2017). Numerical tests reveal that both variants of the INFT algorithm are capable of dealing with a large number of eigenvalues (within the limitations of the double precision arithmetic) previously unreported. Further, for the cases considered in this article, our algorithms perform better than the TIB algorithm Belai et al. (2007); Frumin et al. (2015) in terms of accuracy while being faster by an order of magnitude. Let us also note that the TIB algorithm has no consequence for the fast inverse NFT in the general case.

Next, let us mention that we have not included simulations of a realistic optical fiber link in order to demonstrate the effectiveness of our algorithms. A thorough testing for various NFT-based modulation schemes for a realistic optical fiber link is beyond the scope of this paper. This omission however does not impact the study of the limitation of the proposed algorithms from a numerical analysis perspective.

Future research on fast INFTs will further focus on the stability properties of the LP algorithm and the DT iterations. Moreover, we would also like to consider other linear multistep methods to obtain a higher-order convergent forward/inverse NFTs. The implicit Adams method used in this paper for the purpose of testing already demonstrates that such possibilities do exist, at least, for the solution of the direct ZS problem.

Appendix A Implicit Adams Method

In order to develop the numerical scheme based on the implicit Adams (IA) method, we begin with the transformation $\tilde{\bm{v}}=e^{i\sigma_{3}\zeta t}\bm{v}$ so that the ZS problem in (1) reads as

[TABLE]

Let the grid $\{t_{n}\}$ be as defined in III, and, set $U_{n}=U(t_{n})$ and $\widetilde{U}_{n}=\widetilde{U}(t_{n})$. The discretization of (36) using the $m$-step IA method ($m\in\{1,2,3\}$) reads as

[TABLE]

where $\bm{\beta}=(\beta_{0},\beta_{1},\ldots,\beta_{m})$ are known constants (Hairer et al., 1993, Chap. III.1) (also summarized in Table 1). Solving for $\tilde{\bm{v}}_{n+m}$, we have

[TABLE]

or, equivalently,

[TABLE]

where $\sigma_{0}=\operatorname{diag}(1,1)$. The individual matrices can be worked out as

[TABLE]

where $Q_{n}=(h\beta_{m})q_{n}$, $R_{n}=(h\beta_{m})r_{n}$, $\Theta_{n}=1-Q_{n}R_{n}$ and

[TABLE]

Also,

[TABLE]

The $m$-step IA methods lead to transfer matrices $\mathcal{M}_{n}\in\mathbb{C}^{2m\times 2m}$ of the form

[TABLE]

where $M^{(s)}_{n+m}(z^{2})\in\mathbb{C}^{2\times 2}$ so that

[TABLE]

where $\bm{w}_{n}=z^{n}\bm{v}_{n}$ and

[TABLE]

Let us consider the Jost solution $\bm{\phi}(t;\zeta)$. We assume that $q_{n}=0$ for $n=-m+1,-m+2,\ldots,0$ so that $\bm{\phi}_{n}=z^{\ell_{-}}z^{-n}(1,0)^{\intercal}$ for $n=-m+1,-m+2,\ldots,0$. The discrete approximation to the Jost solution can be expressed as $\bm{\phi}_{n}=z^{\ell_{-}}z^{-n}\bm{P}_{n}(z^{2})$. The initial condition works out to be

[TABLE]

yielding the recurrence relation

[TABLE]

where $\boldsymbol{\mathcal{P}}_{n}(z^{2})=(\bm{P}_{n}(z^{2}),\bm{P}_{n-1}(z^{2}),\ldots,\bm{P}_{n-m+1}(z^{2}))^{\intercal}\in\mathbb{C}^{2m}$. The discrete approximation to the scattering coefficients is obtained from the scattered field: $\bm{\phi}_{N}=(a_{N}z^{-\ell_{+}},b_{N}z^{\ell_{+}})^{\intercal}$ yields $a_{N}(z^{2})={P}^{(N)}_{1}(z^{2})$ and $b_{N}(z^{2})=(z^{2})^{-\ell_{+}}{P}^{(N)}_{2}(z^{2})$. The quantities $a_{N}$ and $b_{N}$ are referred to as the discrete scattering coefficients uniquely defined for $\operatorname{Re}\zeta\in[-{\pi}/{2h},\,{\pi}/{2h}]$.

Finally, let us mention that, for $\zeta$ varying over a compact domain, the error in the computation of the scattering coefficients can be shown to be $\mathop{\mathscr{O}}\left(N^{-p}\right)$ provided that $q(t)$ is at least $p$-times differentiable (Hairer et al., 1993, Chap. III).

It is evident from the preceding paragraph that the forward scattering step requires forming the cumulative product: ${\mathcal{M}}_{N}(z^{2})\times{\mathcal{M}}_{N-1}(z^{2})\times\ldots\times{\mathcal{M}}_{2}(z^{2})\times{\mathcal{M}}_{1}(z^{2})$. Let $\bar{m}$ denote the nearest base-$2$ number greater than or equal to $(m+1)$, then pairwise multiplication using FFT Henrici (1993) yields the recurrence relation for the complexity $\varpi(n)$ of computing the scattering coefficients with $n$ samples: $\varpi(n)=8m^{3}\nu(\bar{m}n/2)+2\varpi(n/2),\,\,n=2,\,4,\,\ldots,\,N,$ where $\nu(n)=\mathop{\mathscr{O}}\left(n\log n\right)$ is the cost of multiplying two polynomials of degree $n-1$ (ignoring the cost of additions). Solving the recurrence relation yields $\varpi(N)=\mathop{\mathscr{O}}\left(m^{3}N\log^{2}N\right)$.

Finally, the results of the tests for benchmarking are shown Fig. 11.

Appendix B An extension of the theorem of Epstein

In the following, we would like to extend Theorem $4$ of Epstein (2004) to obtain the result (29). Define the nonlinear impulse response

[TABLE]

and assume $p(\tau)\in\mathsf{L}^{1}\cap\mathsf{L}^{2}$. Consider the Jost solutions with prescribed asymptotic behavior as $x\rightarrow\infty$:

[TABLE]

where $\bm{A}$ is independent of $\zeta$. Our starting point for the analysis of the inverse problem would be the Gelfand-Levitan-Marchenko (GLM) integral equations. In the following we fix $t\in\mathbb{R}$ so that the GLM equations for $y\in\Omega_{t}=[t,\infty)$ is given by

[TABLE]

where $f(\tau)=p(-\tau)$. The solution of the GLM equations allows us to recover the scattering potential using $q(t)=-2A_{1}(t,t)$ together with the estimate $\|q\chi_{[t,\infty)}\|^{2}_{2}=-2A_{2}(t,t)$ where $\chi_{\Omega}$ denotes the characteristic function of $\Omega\subset\mathbb{R}$. Define the operator

[TABLE]

whose Hermitian conjugate, denoted by $\mathscr{P}^{\dagger}$, works out to be

[TABLE]

Define $\mathscr{K}=\mathscr{P}^{\dagger}\circ\mathscr{P}$, so that

[TABLE]

where the kernel function $\mathcal{K}(y,x;t)$ is given by

[TABLE]

The GLM equations in (47) can now be stated as

[TABLE]

which is a Fredholm integral equation of the second kind where

[TABLE]

Let $\mathcal{I}_{m}(t)=\|f\chi_{[2t,\infty)}\|_{\mathsf{L}^{m}}$ for $m=1,2,\infty$, then

[TABLE]

and, $\|\Phi_{2}(t,\cdot)\|_{\mathsf{L}^{\infty}(\Omega_{t})}\leq[\mathcal{I}_{2}(t)]^{2}$. If $\mathcal{I}_{1}(t)<1$, then the standard theory of Fredholm equations suggests that the resolvent of the operator $\mathscr{K}$ exists Gripenberg et al. (1990). Under this assumption, certain estimates for $q(t)$ can be easily obtained Epstein (2004): From (52), we have

[TABLE]

which yields

[TABLE]

Given that from here one can only assert that $|A_{j}(t,y)|\leq\|{A}_{j}(t,\cdot)\|_{\mathsf{L}^{\infty}(\Omega_{t})}$ almost everywhere (a.e.), we need to ascertain the continuity of $A_{j}(t,y)$ with respect to $y$ throughout the domain $\Omega_{t}$ or as $y\rightarrow t$ from above. Assume that $f(\tau)$ is continuous, then $\Phi_{j}(t,y)$ is continuous with respect to $y$. It can be seen that the kernel function $\mathcal{K}(y,x;t)$ is also continuous with respect to $y$. Therefore, if the resolvent kernel is continuous (w.r.t. $y$) then the result follows. To this end, consider the Neumann series for the resolvent $\mathscr{R}=\sum_{n\in\mathbb{Z}_{+}}(-1)^{n}\mathscr{K}_{n}$ where $\mathscr{K}_{n}=\mathscr{K}\circ\mathscr{K}_{n-1}$ with $\mathscr{K}_{1}=\mathscr{K}$. For fixed $t$, the partial sums $\sum_{1\leq n\leq N}\|\mathscr{K}_{n}\|_{\mathsf{L}^{\infty}(\Omega_{t})}\leq[1-\mathcal{I}^{2}_{1}(t)]^{-1}$ for all $N<\infty$. Therefore, uniform convergence of the partial sums allows us to conclude the continuity of the limit of the partial sums.

Now using the identities $q(t)=-2A_{1}(t,t)$ and $\|q\chi_{[t,\infty)}\|^{2}_{2}=-2A_{2}(t,t)$, we have

[TABLE]

If $\mathcal{I}_{1}(t)<1$ does not hold for all $t\in\mathbb{R}$, one can find a $T>0$ such that $\mathcal{I}_{1}(t)<1$ holds for $t\in[T,\infty)$. The estimates obtained above would then be valid in $[T,\infty)$.

The second inequality in (55) can be used to choose the computational domain for the inverse NFT. Let us consider the example considered in Sec. IV.2: The nonlinear impulse response works out to be

[TABLE]

Note that $p_{\text{rc}}(-\tau)=p_{\text{rc}}(\tau)$. From the asymptotic form

[TABLE]

it follows that

[TABLE]

If $\mathcal{I}_{1}(T)\ll 1$, then setting $\mathcal{I}^{2}_{2}(T)=\epsilon$ gives

[TABLE]