Nonlinear Self-Interference Cancellation with Learnable Orthonormal Polynomials for Full-Duplex Wireless Systems

The contributions of this paper are summarized as follows:

• •

  We first show the existence of a set of orthonormal nonlinear basis functions which relate the transmit signal to received SI. The key idea to constructing the proposed basis functions is to generalize the IH polynomials according to the (possibly time-varying) moments of the transmitted signal.

• •

  Then, we propose a computationally efficient algorithm that constructs the orthonormal nonlinear basis functions. By harnessing the recursive structure in identifying the coefficients of IH polynomials, we show that the construction of $p$th-order orthonormal polynomial is possible with a computational complexity of $\mathcal{O}(p^{2})$.

• •

  Cascading the constructed orthonormal basis functions with the LMS algorithm, we present a SIC algorithm for non-stationary input data, which we refer to as the AOP-LMS algorithm. In the moment learning phase of this algorithm, the moments of transmitted data symbols are estimated and then the set of orthonormal basis functions is systematically constructed. Then, in the filter coefficient learning phase, our algorithm uses the LMS algorithm to refine and adapt these coefficients.

• •

  We further accelerate the proposed approach with a look-up table containing precomputed orthonormal basis functions for various signal distributions (e.g, $M$-QAM). In practice, such a look-up table can be constructed a priori based on the known set of modulation and coding schemes (MCSs) employed by a transceiver and can then be referenced in real time as the MCS changes.

• •

  Using numerical simulation, we show that our proposed SIC algorithm provides significant improvement in both mean squared error (MSE) and convergence speed compared to the state-of-the-art algorithms, including the HP-RLS algorithm, under non-stationary input distributions.

• •

  We also verify the effectiveness of the proposed SIC algorithm using a real-time wireless testbed. These experimental results further confirm that our proposed approach—and its ability to adapt to changes in modulation—translates from theory to implementation.

This paper is organized as follows. Section II defines the system model of the FD wireless communication system and motivates the need for pair-wise orthogonality with respect to the basis of nonlinear LMS algorithms. Section III presents a method to construct orthonormal polynomials using partial moments of the input signal. Section IV provides orthonormal polynomials for various input signal distributions. Section V introduces the proposed digital SIC algorithms based on the formulations laid forth in Sections III-IV. Section VI and Section VII show results of the proposed SIC algorithms through numerical simulation and real-world experiments. Section VIII concludes the paper.

The main task of digital SIC is to find the best approximation of the nonlinear function $f(\cdot)$ in an online manner—a problem akin to those appearing in contexts of adaptive filter theory as system identification [24]. To develop a model for this nonlinear function approximation, we harness established models for nonlinear PAs and linear filters.

We begin by considering Saleh’s PA model, a time-honored wideband PA model with memory effects initially introduced in [34, 35]. Under such a model, the output of a PA with memory $M$, whose input is $x[n]$, is modeled as

where $\beta\geq 0$ and $\gamma\geq 0$ capture the transition sharpness and small signal gain of the PA, respectively. Here, $h^{\sf PA}[m]$ is the $m$th coefficient of the PA impulse response. We assume that $\beta,\gamma,$ and $\{h^{\sf PA}[m]\}$ are initially unknown to the system and may vary with time (e.g., due to temperature).

Comparatively, Saleh’s model introduces more significant AM/PM distortion than most off-the-shelf solid-state PAs [34, 36, 37]. Consequently, employing this model will provide flexibility and robustness in capturing (potentially other) sources of nonlinearity beyond solely AM/AM distortion.

Using the truncated power series expansion with maximum degree $P$ (odd integer), we obtain an approximation of the nonlinear PA transfer function (2) in a form of the linear combination of $\frac{P+1}{2}$ HPs $|x[n]|^{2p}x[n]$ for $p\in\left\{0,1,\ldots,\frac{P-1}{2}\right\}$:

Without loss of generality, we generalize the $p$th order nonlinear basis function as a linear combination of the HPs with maximum degree $p$ (odd integer), namely

where $c_{p,k}\in\mathbb{C}$ is the $k$th coefficient of the $p$th basis function and ${\bf c}_{p}=[c_{p,0},\ldots,c_{p,p-1}]$. This overparameterized representation of the basis function provides additional degrees of freedom to capture PA nonlinearity. For instance, when $c_{p,k}=0$ for $k\in\left\{0,1,\ldots,p-2\right\}$ and $c_{p,p-1}=1$, we obtain the standard HPs as the basis functions, i.e., $\phi_{p}^{\sf HP}(x[n])=|x[n]|^{2(p-1)}x[n]$. By tuning $c_{p,k}$, different classes of HPs can be generated as basis functions. Ignoring the higher-order approximation error in (3), we can express the PA function in (3) by the sum of the $\frac{P+1}{2}$ overparameterized nonlinear basis functions $\phi_{p}(x[n])$ as

Plugging (5) into (2), we can rewrite the PA output using the parallel HP with memory length $M$ as

where $h_{p}^{\sf PA}[m]$ is the impulse response of the PA for the $p$th order nonlinear input $\phi_{p}(x[n])$. The latter equality comes from our overparameterization approach, which will allow us to capture the input-dependent PA memory response effects $h_{p}^{\sf PA}[m]$ for $p\in\{1,2,\dots,\frac{P+1}{2}\}$.

Let $\{h_{\sf SI}[q]\}$ be the length-$Q$ impulse response of the SI propagation channel. Under an idealized LNA and ADC¹¹1PAs are often assumed the dominant sources of nonlinear SI [33, 22]; our model readily accommodates other sources of nonlinearity, nonetheless., the received baseband SI signal is approximately

Invoking (7) into (8), this approximated SI is equivalently

where $\{h_{p}[\ell]\}$ is the FIR filter containing the combined effects of $h_{p}^{\sf PA}[m]$ and $h_{\sf SI}[q]$, with total memory $L=M+Q-1$. Our approximation of the SI signal in (9) contains two sets of unknown parameters:

1. 1.

   The coefficients for constructing the nonlinear basis functions $\Theta_{1}=\left\{{\bf c}_{1},\ldots,{\bf c}_{\frac{P+1}{2}}\right\}$, with $|\Theta_{1}|=\frac{(P+1)(P+3)}{8}$.

2. 2.

   The coefficients of the effective FIR filter, $\Theta_{2}=\{\{h_{1}[\ell]\},\{h_{2}[\ell]\},\dots,\{h_{\frac{P+1}{2}}[\ell]\}\}$, with $|\Theta_{2}|=\frac{(P+1)L}{2}$.

We denote our approximation of the effective nonlinear SI channel in (1) by ${\hat{f}}(\cdot;\Theta_{1},\Theta_{2}):\mathbb{C}^{L}\rightarrow\mathbb{C}$ with parameters $\Theta_{1}$ and $\Theta_{2}$.

We emphasize that this model-based function approximation method differs from the model-free function approximation techniques with deep neural networks (DNNs) [38, 39], which approximate SI without exploiting any prior model knowledge. Our model-based method, on the other hand, leverages select parameters based on knowledge of established PA models and the FIR filters of linear time-invariant systems. Ultimately, this proves both numerically and experimentally to make our proposed approach more adaptive and effective in canceling SI, as we will see in Sections VI and VII.

Let $e[n]$ be the error between the true received SI $y[n]$ and the SI constructed by our approximation with input ${\bf x}_{L}[n]$:

With this function approximation in hand, our goal is to find the parameter sets $\Theta_{1}$ and $\Theta_{2}$ which minimize the empirical squared error

where $\mathcal{T}=\{n_{1},n_{1}+1,\ldots,n_{2}\}\subset\mathbb{Z}$ is an index set capturing the time interval of interest. In general, finding a jointly optimal solution of $\Theta_{1}$ and $\Theta_{2}$ is very challenging due to the non-convexity of $J(\Theta_{1},\Theta_{2})$ with respect to both parameter sets. In light of this difficulty, we instead optimize them in a disjoint manner using classical adaptive filter theory.

Let ${\bm{\phi}}_{p}({\bf x}_{L}[n];{\bf c}_{p})=\left[\phi_{p}(x[n];{\bf c}_{p}),\phi_{p}(x[n-1];{\bf c}_{p}),\ldots,\right.$ $\left.\phi_{p}(x[n-L+1];{\bf c}_{p})\right]$ be the filter input vector generated by the $p$th order basis function and ${\bf h}_{p}[n]=\left[h_{p}[n],h_{p}[n-1],\ldots,h_{p}[n-L+1]\right]$ be the filter response vector of the $p$th order basis function. By concatenating the input and filter response vectors, we can also define

and ${\bf h}[n]=\left[{\bf h}_{1}[n],\ldots,{\bf h}_{\frac{P+1}{2}}[n]\right]^{\top}\in\mathbb{C}^{\frac{L(P+1)}{2}\times 1}$. Then, for given $\Theta_{1}=\left\{{\bf c}_{1},\ldots,{\bf c}_{\frac{P+1}{2}}\right\}$, the approximation of SI in (9) can be rewritten in a linear form with respect to ${\bf h}[n]$ as

In the next section, we delineate the methodology for optimizing the basis parameter $\Theta_{1}$ of $\bm{\phi}$, leveraging the characteristics of the transmitted signals $\{\mathbf{x}[n]\}$.

Using the moments of the transmit signal $x[n]$, our goal is to construct a set of orthonormal basis functions $\phi_{1}(x[n];\mathbf{c}_{1}),\phi_{2}(x[n];\mathbf{c}_{2}),\ldots,\phi_{\frac{P+1}{2}}(x[n];\mathbf{c}_{\frac{P+1}{2}})$ that satisfy the following conditions:

Before proceeding, we first define a moment vector $\bm{\mu}_{2a}^{2b}$ and a Hankel matrix $\mathbf{M}_{p}$ as

and

where $a\leq b$, $p>1$, and $\mu_{2p}=\mathbb{E}[|x[n]|^{2p}]$. Now, we provide a theorem for finding polynomial coefficeints which satisfy the orthonormal condition (14), given known moments $\boldsymbol{\mu}_{2a}^{2b}$.

From Theorem 1, we can construct the coefficients of a polynomial satisfying (14) when matrix $\mathbf{M}_{p}$ is nonsingular. The process of sequentially finding the coefficients of orthogonal polynomials using the matrix equation is similar to the orthogonalization algorithm in traditional linear algebra. Extending the notion of a projection in a vector space, this method orthonormalizes in a function space, where the inner product is defined as in Fig. 2.

The advantage of this method is that it can find the orthonormal polynomials with only knowledge of the moments. However, one drawback arises from the need to perform several successive matrix inversions for higher-order orthonormal polynomials. In the following subsection, we propose an equivalent yet more efficient method which reduces the complexity of such successive matrix inversions.

In light of the computational costs of successive matrix inversions, we now propose an orthonormal polynomial construction method involving two steps: i) recursive computation of the orthogonal basis functions and ii) basis function normalization.

Recursive Computation: From Theorem 1, we construct the matrix equation satisfying the orthogonal condition in (14). To obtain the basis function coefficients, we need to solve the matrix equations in (18) for each $p\in\{1,2,\ldots,\frac{P+1}{2}\}$. The computational complexity for this operation scales considerably as the maximum nonlinear order $\frac{P+1}{2}$ increases. To reduce the computational complexity in constructing the orthonormal basis functions, we present a recursive algorithm using the structure of the Hankel matrix ${\bf M}_{p}$. In essence, we will use $\mathbf{M}_{p}$ to construct ${\bf M}_{p+1}$. Let $\mathbf{u}_{p}\triangleq\boldsymbol{\mu}_{2p}^{4p-4}$ and define a scalar $r_{p}\triangleq\mu_{4p-2}$. Then, $\mathbf{M}_{p+1}$ can be constructed as

Using the fact that ${\bf M}_{p}$ is the Schur complement [40] of $\mathbf{M}_{p+1}$, the inverse of $\mathbf{M}_{p+1}$ can be recursively computed using the inverse of $\mathbf{M}_{p}$ as

where $\tilde{\mathbf{u}}=\mathbf{M}_{p}^{-1}\mathbf{u}_{p}$ and $s_{p}=r_{p}-\mathbf{u}_{p}^{\top}\tilde{\mathbf{u}}_{p}$. By harnessing the structure of the Hankel matrix ${\bf M}_{p}$, we can compute a set of the orthogonal basis functions in this computationally efficient manner. This recursive computation of $\{\mathbf{M}_{p}^{-1}\}$ is summarized in Algorithm 1.

Normalization: Once the orthogonal basis functions are computed by the proposed recursive method, it is necessary to normalize the polynomial basis function to have unit norm. With this normalization, we obtain the orthonormal basis functions $\phi^{\sf AOP}_{p}(x;\hat{\mathbf{c}}_{p})$ as

where $\hat{\mathbf{c}}_{p}$ is an orthonormal coefficient vector of the $p$th order polynomial. The normalization coefficient $z$ in (19) is computed by

where $\tilde{c}_{p,k}=\sum_{i}[\mathbf{c}_{p}]_{i}[\mathbf{c}_{p}]_{k-i}$ is the self-convolution of $\mathbf{c}_{p}$.

Calculating coefficients for the $p$th orthonormal polynomial using conventional techniques such as Gaussian elimination involves a computational complexity of $\mathcal{O}(\frac{2}{3}p^{3})$. Employing our proposed method reduces this complexity to $\mathcal{O}(3p^{2})$.

Let $X$ be a complex Gaussian distributed random variable with mean zero and variance $\sigma_{x}^{2}$, denoted as $\mathcal{CN}(0,\sigma_{x}^{2})$, which is perhaps the most widely used distribution in communication theory. Recall, the only information necessitated by our algorithm to construct an orthonormal polynomial are the moments of the distribution. The $m$th-order moment, $\mathbb{E}[|X|^{m}]$, has a closed-form expression for even $m$ [41], given by

For the particular case of $\sigma_{x}^{2}=1$, the first three orthonormal polynomials obtained by our proposed algorithm are then

We point out that the orthonormal polynomials of complex Gaussian random variables have been studied in prior literature under the moniker of It$\hat{\mbox{o}}$-Hermite polynomials [42].

The random variable $X$ uniformly distributed over the interval $[-k,k]$, where $k>0$, has a probability density function (PDF) given by

It can be shown straighforwardly that $\mathbb{E}[X^{2m}]=\frac{1}{2m+1}k^{2m}$. Since the matrix $\mathbf{M}_{p}$ in (16) is non-singular when populated with moments of the uniform distribution, the corresponding polynomial coefficients can be obtained directly from Algorithm 1. For the particular case of $k=1$, the first three orthonormal polynomials obtained by Algorithm 1 are:

Note that it is well known that the orthonormal polynomials of the uniform distribution are the Legendre polynomials [43], with (IV-B) being a scaled version of such.

Suppose $X$ is an exponential random variable parameterized by $\lambda>0$, having PDF

The moment generating function (MGF) is well known to be $M_{X}(t)=\frac{\lambda}{\lambda-t}$. Moments of this distribution can be derived from the MGF as $\mathbb{E}[X^{m}]=\frac{\partial^{m}M_{X}(t)}{\partial t^{m}}|_{t=0}$ .

Note that distributions with non-zero mean, such as this, may yield multiple valid orthonormal polynomials, depending on if one considers only even degrees or both even and odd degrees. In the former case, the first polynomial is denoted as $\phi_{1}(x)=x$, whereas in the latter case, the degree of the first polynomial is $0$, leading it to be denoted as $\Phi_{0}=c_{0,0}$, where $\Phi_{p}$ is extended basis function defined as

The proposed algorithm can be directly extended to construct this expanded basis, $\Phi_{p}(\cdot)$, by including both the odd and even moments when generating the Hankel matrix (16). For the exponential distribution where $\lambda=1$, the orthonormal polynomials for the basis $\phi_{p}$ (having only even degrees) are

and $\Phi_{p}$ (having both even and odd degrees) are

This example illustrates that, even with the same distribution, more than one valid orthonormal polynomial can exist, depending on the form of the basis function. Note that the latter polynomial set in (IV-C) is well known as the Laguerre polynomial [44].

Other particularly relevant signal distributions to consider are those of 4QAM, 16QAM, 64QAM, and 256QAM, corresponding to digital modulation techniques used widely in modern communication systems. Let the constellation (set) of possible modulation symbols be defined as $\mathcal{S}=\{s_{1},s_{2},\dots,s_{i}\}$, where $i\in\{4,16,64,256\}$ is the modulation order. When symbols are drawn uniformly from a given constellation, the $m$th moment has the closed-form

For the particular case of 4QAM, the Hankel matrix in (16) is a rank-$1$ matrix. This is perhaps most obvious when the four symbols are on the unit circle, in which case they all have a second-order moment (or power) of $1$ and thus the higher-order moments are also $1$. Put simply, orthonormal polynomials of the third-order or higher do not exist for 4QAM. That is, the higher-order polynomials in (4) are

Thus, when transmitting 4QAM symbols, the nonlinearity is captured solely by the first-order polynomial $x$. Orthonormal polynomials for higher-order QAM constellations are summarized in Table I.

As a first step, our proposed SIC algorithm estimates the moments of the transmit symbols $x[n]$. The $p$th moment of the transmit symbol $x[n]$ is estimated by its sample average:

This estimation converges to the true moment as the sample size $N$ increases by the law of large numbers, i.e.,

Even for modest $N$, this sample average estimator can often closely approximate the true moment since estimation errors reduce linearly with sample size [45]. Then, the basis functions have orthogonality by constructing the coefficients according to Algorithm 1 as ${\bm{\phi}}^{\sf AOP}({\bf x}[n];\Theta_{1})=\left[\phi^{\sf AOP}_{1}(x[n],{\hat{\bf c}}_{1}),\ldots,\phi^{\sf AOP}_{\frac{P+1}{2}}(x[n],\hat{{\bf c}}_{\frac{P+1}{2}})\right]^{\top}$. Theorem 2 is provided below to show that SI can indeed be canceled using the proposed orthonormal polynomial basis functions.

To derive the weights $\{\mathbf{h}_{k}[n]\}$ which satisfy Theorem 2, we introduce an adaptive method with our proposed orthonormal basis function. The LMS algorithm is a stochastic approximation of the iterative steepest descent based implementation of the Wiener filter and is applicable when the SI and the input signal are jointly wide-sense stationary (JWSS) [24, 46, 47]. This stochastic approximation involves a simple update equation that can be implemented in practical systems with low computational complexity. Using the classical LMS algorithm [24], the linear filter parameter coefficients are updated as

where $\mu$ is the step size.

Remark 1 on Computational Complexity: The computational cost associated with implementing orthonormal polynomials in the AOP-LMS algorithm is $\mathcal{O}\left(\frac{P^{3}}{8}\right)$, where the majority of this complexity involves constructing the orthonormal basis functions.

Filter Adaptation: Based on adaptive filter theory [24], the performance of the LMS algorithm may vary depending on the conditioning of the covariance matrix of the input basis functions. That is, in an environment in which the distribution changes with time, using a fixed set of basis functions may cause degradation in SIC performance. This motivates the need for a method which generates orthonormal basis functions that adapt to changes in the distribution of the input signal. Fig. 3 depicts a block diagram of one such method. In this proposed method, the estimated moments of the transmit signal are updated by comparing the previous moment estimates to the current estimates, updating them based on some specified condition. Technically speaking, a change in the moments necessitates an update of the basis functions, but one could capture slight changes in the signal distribution by instead updating $\{\mathbf{h}_{k}[n]\}$. This avoids recomputing the polynomial basis functions and instead leverages the computational simplicity of adaptive algorithms, such as LMS.

As one example of such an approach, Algorithm 2 describes a method which uses nonlinear LMS to regularly update the estimated moments only in a predetermined interval. In practice, this interval could correspond to the duration of time slots defined by wireless standards, such as 5G and Wi-Fi. We define $N_{\mathrm{max}}$ as the maximum number of samples used to estimate the transmit signal’s statistics, $N_{\mathrm{int}}$ as the sample interval at which statistical information is fixed, and $n_{\mathrm{s}}$ as the start sample for statistical information collection.

LUT Adaptation: As mentioned before, in practical wireless systems, transmit signals are non-stationary due to link adaptation, whereby a transmitter adapts its MCS according to the channel strength/quality. Since the set of possible MCSs are known a priori, the transmit signal’s statistics and thus the orthonormal polynomials can be pre-computed for each MCS. Storing these pre-computed polynomials in a look-up table (LUT) and then referencing them at run-time can reduce computational complexity. Our proposed method employing such a technique is depicted in Fig. 4.

In Table II, we compare existing digital SIC algorithms against our proposed technique in terms of model complexity and performance. While the HP-W-LMS algorithm and our proposed AOP technique are on par with one another in terms of complexity, we will see in the next section that ours offers superior robustness/adaptation to changes in the transmit signal distribution and in the SI propagation channel.