text:
\begin{abstract}
Full-duplex transmission comprises the ability to transmit and receive at the same time on the same frequency band. It allows for more efficient utilization of spectral resources, but raises the challenge of strong self-interference (SI). Cancellation of SI is generally implemented as a multi-stage approach. This work proposes a novel adaptive SI cancellation algorithm in the digital domain and a comprehensive analysis of state-of-the-art adaptive cancellation techniques. Inspired by recent progress in acoustic echo control, we introduce a composite state-space model of the nonlinear SI channel in cascade structure.
We derive a SI cancellation algorithm that decouples the identification of linear and nonlinear elements of the composite state. They are estimated separately and consecutively in each adaptation cycle by a Kalman filter in DFT domain. We show that this adaptation can be supported by a-priori signal orthogonalization and decoding of the signal-of-interest (SoI).
In our simulation results, we analyze the performance by evaluating residual interference, system identification accuracy and communication rate. Based on the results, we provide recommendations for system design.
In case of input orthogonalization, our Kalman filter solution in cascade structure delivers best performance with low computational complexity. In this configuration, the performance lines up with that of the monolithic (parallel) Kalman filter or the recursive-least squares (RLS) algorithms. We show that the Kalman-based algorithm is superior over the RLS under time-variant conditions if the SoI is decoded and in this way the covariance information required by the Kalman filter can be provided to it.

\end{abstract}

\section{Introduction}
The sheer amount of data conveyed over the wireless medium requires novel concepts to be integrated into 5G (and beyond) transmission standards. A promising technique for more efficient spectrum utilization is RF (radio frequency) in-band full-duplex communication, where transmission and reception of signals are run simultaneously on the same frequency band~\cite{src:sexton20175g}. Theoretically such an approach supports an improved capacity of the wireless link~\cite{src:sabharwal2014band,src:aijaz2017simultaneous}. However, some considerable signal power from the transmitter inevitably leaks as self-interference (SI) into the receiver chain, and therefore any signal-of-interest (SoI) from a distant communication node cannot be reliably recovered at the receiver. Cancellation of the strong SI is a challenging task. Several landmark papers have successfully demonstrated the feasibility of such an approach~\cite{src:jain2011practical,src:riihonen2011hybrid,src:duarte2012experiment}. Common state-of-the-art SI cancellation is generally implemented as a multi-stage operation. Full-duplex prototypes have combined passive SI suppression, active analog cancellation, and digital cancellation~\cite{src:khandani2013two,src:duarte2014design,src:debaillie2014analog,src:ahmed2015all,src:askar2016agile}.

The availability of SI cancellation allows for a wide range of full-duplex applications. It can increase the communication rate in a wireless network by using full-duplex relays~\cite{src:avestimehr2008approximate,src:nunn2017antenna}, reduce end-to-end delay~\cite{src:kariminezhad2017fullduplex}, provide enhanced security at the physical layer by friendly jamming~\cite{src:zheng2013improving} or key agreement~\cite{src:vogt2016practical}, and it implies the potential of energy-harvesting~\cite{src:zeng2015full, src:bi2016accumulate}.

To remove the impact of SI in the digital domain, a variety of digital SI channel estimation metrics have been studied such as least-squares (LS)~\cite{src:bharadia2013full}, minimum-mean squared error (MMSE)~\cite{src:day2012fullduplex} or maximum-likelihood (ML) estimation~\cite{src:masmoudi2016maximum}. Non-uniform sampling can avoid undesired SI in Orthogonal Frequency-Division Multiplexing (OFDM) systems~\cite{src:bernhardt2018self}. Under time-variant conditions, adaptive algorithms generally can cover the task of adjusting an SI estimate to gradual changes in the SI channel~\cite{src:heino2015recent}. The linear SI channel model has been extended to nonlinear forms, where the Hammerstein polynomial model is the most popular~\cite{src:bharadia2014full}. It takes the form of a parallel, multi-channel representation with monolithic structure. Alternatively, the nonlinear SI channel can be characterized by an artificial neural network~\cite{src:stimming2018nonlinear}. Beyond the nonlinear SI model, previous work established enhancements like signal orthogonalization, least-mean square (LMS) adaptation~\cite{src:korpi2015adaptive,src:ferrand2017multi} or recursive-least square (RLS) algorithms~\cite{src:lemos2015fullduplex,src:emara2017nonlinear}. Adaptive approaches can be used in hybrid digital/analog cancellation designs~\cite{src:kiayani2018adaptive}. Many prototypes use an analog cancellation stage to remove the line-of-sight component, which is the dominant and, in most cases, the time-invariant part of the SI. However, moving objects close to the antenna can significantly change the multipath characteristics of the SI channel due to power backscattering~\cite{src:everett2014passive} and therefore introduce a significant time-variant contribution to the SI. The general problem of adaption to time-variant SI channels requires further research in the digital domain.

The SI is a challenge not only in wireless communication. Consider hands-free voice communication, where microphone and loudspeaker signals require full-duplex operations. Here, similar to the SI in wireless full-duplex applications, the unacceptable acoustic echo signal needs to be effectively removed. It was, however, found
that the eventual echo cancellation system involves a number of challenges, such as the fast adaptivity to time-variant loudspeaker-room-microphone systems~\cite{src:benesty2001advances,Haensler97,Haensler2006}.
The standard for acoustic echo control is the LMS and normalized least means-square (NLMS) type adaptive-filter algorithm in single- and multi-channel configurations~\cite{Haykin2002,Enzner2014}. Similarly, RLS-based approaches of the multi-channel type have been introduced~\cite{src:benesty2001advances}.  As a compromise between NLMS and RLS in terms of fast adaptation and computational complexity, approaches based on affine projection algorithms (APA)~\cite{src:gay1995the} have been proposed and applied to nonlinear echo cancellation~\cite{src:glicacho2012nonlinear}. State-space modeling of the time- and frequency-domain echo channel~\cite{Enzner06,Enzner2010} has been introduced.

Variants of the Kalman filter algorithm were provided in multichannel~\cite{Malik2011a} and nonlinear configurations \cite{MalikEnzner2012b,src:malik2013variational}.

Both the domains of acoustic echo control and wireless SI cancellation share many aspects of the system modeling. Thus, we can transfer concepts and insights from the acoustic echo cancellation. Most prominently, the nonlinear echo/SI models of both the acoustic loudspeaker and the RF power amplifiers are very similar. In the acoustic echo channel, however, impulse responses can have up to a few thousand of significant taps. In the wireless domain, the delay spread is much shorter, but the accurate depiction of fractional signal delays still requires a significant amount of taps~\cite{src:laakso96splitting}. There is a notable difference between acoustic echo control and wireless SI cancellation: in the handling of the SoI. In acoustics, the SoI is usually modeled as a random process and therefore contributes to the observation noise during adaptation. In wireless communications, however, the SoI is deliberately encoded and thus can be removed before adaptation.

In this work, we set our focus on the digital part of SI cancellation by an adaptive algorithm. Consider the system overview of~\figref{fig:overview}.
\def\antenna{
   -- +(0mm,4.0mm) -- +(2.625mm,7.5mm) -- +(-2.625mm,7.5mm) -- +(0mm,4.0mm)
}
\begin{figure}
	\centering
	\begin{tikzpicture}[
	sysBlock/.style={draw, fill=white, rectangle, minimum height=2em, minimum width=3.5em},
	nlBlock/.style={draw, fill=white, rectangle, minimum height=2em, minimum width=4em},
	shadowBlock/.style={draw, drop shadow, fill=white, rectangle},
	triangle/.style = {draw, regular polygon, regular polygon sides=3 },

	border rotated/.style = {shape border rotate=90}
	]

	\node (inTx) {$x_{\timeVar}$};
	\node (outRx) at ($(inTx)+(0,-3)$) {};

	\node[shadowBlock, minimum height=4em](dec) at ($(inTx)!0.5!(outRx)+(0.75,0.25)$) {Decoder};
	\myadd{add2}{$(dec)+(1.5,0)$}
	\node[shadowBlock, minimum height=4em, minimum width=4em, align=center, fill=red!10] at ($(add2)+(1.5,0)$) (adapt) {Adaptive\\algorithm};
	\myadd{add1}{$(adapt|-outRx)+(0,0)$}
	\node[dot] (dot1) at ($(adapt|-inTx)$) {};
	\node[dot] (dot2) at ($(add2|-outRx)+(0,0)$) {};

	\node  at ($(add1.east)+(0.1,0.2)$) {$-$};
	\node  at ($(add2.west)+(-0.1,0.2)$) {$-$};

	\node[shadowBlock, align=center, minimum height=10em] (analog) at ($(adapt)+(1.6,-0.25)$) {Tx/Rx\\path};
	\node (ant) at ($(analog)+(1.1,0)$) {};
	\draw[very thick] (ant.center) \antenna;
	\draw[very thick] (analog.east) -- (ant.center);

	\node[shadowBlock, align=center, minimum height=2em] (dist) at ($(analog)+(2.3,-1)$) {Distant\\node};
	\node (antDist) at ($(dist)+(0,0.5)$) {};
	\draw[very thick] (antDist.center) \antenna;
	\draw[very thick] (dist.north) -- (antDist.center);

	\draw[thick,->] (inTx) -- (dot1);
	\draw[thick,->] (dot1) -- (analog.west|-dot1);
	\draw[thick,->] (dot1) -- (adapt);

	\draw[thick,->] (analog.west|-add1) -- node[below] {$y_{\timeVar}$} (add1);
	\draw[thick,->] (adapt.south-|add1) -- node[right] {$\hat{x}_{\text{si},\timeVar}$} (add1);
	\draw[thick,->] (dot2) -- (add2);
	\draw[thick,->] (add1) -- node[above] {$e_{\timeVar}$} (dot2);
	\draw[thick,->] (dot2) -| (dec);
	\draw[thick,->] (dec) -- node[below] {$\hat{d}_{\timeVar}^h$} (add2);
	\draw[thick,->] (add2) -- node[above] {$\tilde{e}_{\timeVar}$} (adapt);

	\draw[dashed,->] ($(antDist)+(-0.2,0.2)$) -- node[below,pos=0.75] {$d_{\timeVar}^h$} ($(ant)+(0.2,0.2)$);

	\end{tikzpicture}
\caption{Proposed full-duplex SI estimation and cancellation.}
\label{fig:overview}
\end{figure}
A distant node sends a SoI~$d^h_k$ at time $k$ to the receiver, while the local transmitter conveys $x_k$. At the local receiver, an adaptive algorithm reconstructs an estimated SI~$\hat{x}_{\text{si},\timeVar}$ and subtracts its contribution from the received signal $y_k$, which leaves the error signal $e_k$. Next, this signal $e_k$ is fed into a decoder which retrieves the decoded SoI~$\hat{d}^h_k$. By removing $\hat{d}^h_k$ from the error signal $e_k$, the residual signal $\tilde{e}_k$ is given to the algorithm for adaptation. Inspired by recent progress in the area of acoustic echo control, we intend to tailor acoustic echo cancellation methods to wireless SI cancellation by employing a nonlinear, composite state-space system model in cascade structure.

The state-space representation allows more dedicated treatment of time-variant SI channels, since its a-priori knowledge is embedded directly into the system model.

This leads to a better control over the adaptation process and an improved performance of the solution.

Unlike in the parallel nonlinear SI channel model, the cascade structure is reducing the number of variables to be estimated and simultaneously the susceptibility to over-fitting~\cite{src:malik2013variational}.

We linearize the model by decoupling linear SI path and nonlinear coefficients, which are then estimated separately and iteratively by linear algorithms in each adaptation cycle.
The natural algorithmic solution to the estimations based on MMSE criterion turns out to be the Kalman filter. This Kalman filter algorithm is derived in DFT domain. Frequency domain approaches, while often explicitly tailored to OFDM systems~\cite{src:sohaib2017alow}, have the potential of reducing computational complexity while maintaining estimation accuracy~\cite{src:komatsu2017freq}.

In previous work, adaptive algorithms with input orthogonalization have been proposed~\cite{src:korpi2015adaptive,src:emara2017nonlinear} to meet the aforementioned assumptions. We maintain this feature to decouple the nonlinear basis function for improved speed of convergence. We pursue a systematic treatise of that design feature and explore its impact on the performance.
Furthermore, we integrate the decoding of the signal-of-interest (SoI), which has been conveyed by a distant node. This approach is beneficial, since in wireless full-duplex communication, the SoI contains structure which is a-priori known and therefore can be exploited. In the development of adaptive SI cancellation algorithms, the a-priori knowledge of SoI signal statistics have been rather neglected so far.

We study the performance by considering certain metrics like the signal-to-residual-interference-and-noise ratio, system identification accuracy and the communication rate. Both the time convergence behavior and the global performance with respect to input signal-to-interference-and-noise ratio are considered. We provide comparisons of the proposed algorithm in exact and approximated form to other approaches like Kalman-based, NLMS and RLS algorithms from the literature. This leads to a comprehensive analysis of certain design options  (such as input orthogonalization, complexity or model structure) on the performance of the novel and state-of-the-art adaptive algorithms.

We show that the temporal variations impose a fundamental performance limitation, regardless of whether input orthogonalization is applied or not. The simulation results indicate that the decoding of the SoI is generally beneficial to both speed of convergence and cancellation performance, especially under time-variant conditions.

Throughout the paper, we print column vectors and matrices as bold lower-case and upper-case letters, respectively. The operators $\expv\left[\cdot\right]$, $\tr\left[\cdot\right]$, $\left|\cdot\right|$, $\left\lVert\cdot\right\rVert_2$, $\left(\cdot\right)^T$, $\left(\cdot\right)^H$ denote expectation, trace of a matrix, absolute value, Euclidean norm, matrix transpose and Hermitian transpose, respectively. The operator $\diag\left[ \mybold{x}\right]$ represents a diagonal matrix with the elements of vector $\mybold{x}$ on its main diagonal. Signal vectors in DFT domain are marked by underline $\underline{\mybold{x}}$. The term $\mybold{x}\circ\mybold{y}$ denotes the Hadamard product, i.e., the element-wise multiplication of $\mybold{x}$ and $\mybold{y}$. The identity matrix of size $M\times M$ is given by $\eye{M}$.

The paper is organized as follows. Section~\ref{sec:systemmodel} introduces the system model with more details on the Tx/Rx path of~\figref{fig:overview}, its DFT domain representation and the state-space model. We derive the nonlinear adaptive algorithm in Section~\ref{sec:algorithm} and discuss useful approximations in Section~\ref{sec:approximations}. In Section~\ref{sec:results}, we define suitable metrics and evaluate the performance by considering computational complexity, time convergence behavior, global performance and various use cases. Finally, Section~\ref{sec:conclusion} concludes the paper.

\section{System model}
\label{sec:systemmodel}
Consider the system model of~\figref{fig:systemmodel}, which depicts a nonlinear cascaded SI channel.
\begin{figure}
	\centering
	\begin{tikzpicture}[
	sysBlock/.style={draw, drop shadow, fill=white, rectangle, minimum height=2em, minimum width=3em},
	nlBlock/.style={draw, drop shadow, fill=white, rectangle, minimum height=2em, minimum width=4.5em},
	chanBlock/.style={draw, drop shadow, fill=white, rectangle, minimum height=3em, minimum width=2em}
	]

	\node (in1) {$x_{\timeVar}$};
	\node[dot] (dot2) at ($(in1)+(1,0)$) {};

	\node[nlBlock] at ($(dot2)+(1.35,0.8)$) (nlfunc0) {$\phi_0(\cdot)$};
	\node[nlBlock] at ($(dot2)+(1.35,0)$) (nlfunc1) {$\phi_1(\cdot)$};
	\node at ($(dot2)+(1.35,-0.6125)+(0,0)$) (colon1) {$\colon$};
	\node[nlBlock] at ($(dot2)+(1.35,-1.25)$) (nlfuncN) {$\phi_{N-1}(\cdot)$};

	\mymod{mult1}{$(nlfunc0)+(2.5,0)$}
	\mymod{mult2}{$(nlfunc1)+(2.5,0)$}
	\node at ($(colon1)+(2.5,0)+(0,0)$) (colon1) {$\colon$};
	\mymod{mult3}{$(nlfuncN)+(2.5,0)$}

	\node at ($(mult1)+(-0.75,-0.6)$) (a0) {$a_{0,\timeVar}$};
	\node at ($(mult2)+(-0.75,-0.6)$) (a1) {$a_{1,\timeVar}$};
	\node at ($(mult3)+(-0.75,-0.6)$) (aN) {$a_{N-1,\timeVar}$};

	\myadd{add1}{$(mult2)+(1.5,0)$}

	\node[chanBlock] at ($(add1)+(1,-1)$) (sichan) {$\mybold{w}_{\timeVar}$};

	\draw[thick,->] (in1) to (dot2);

	\draw[thick,->] (dot2) |- (nlfunc0);
	\draw[thick,->] (dot2) -- (nlfunc1);
	\draw[thick,->] (dot2) |- (nlfuncN);
	\draw[thick,->] (a0) -- (mult1);
	\draw[thick,->] (a1) -- (mult2);
	\draw[thick,->] (aN) -- (mult3);
	\draw[thick,->] (nlfunc0) -- (mult1);
	\draw[thick,->] (nlfunc1) -- (mult2);
	\draw[thick,->] (nlfuncN) -- (mult3);
	\draw[thick,->] (mult1) -| (add1);
	\draw[thick,->] (mult2) -- (add1);
	\draw[thick,->] (mult3) -| (add1);

	\myadd{add3}{$(sichan)+(0,-2)$}
	\myadd{add4}{$(add3)+(-1.5,0)$}
	\node (noise) at ($(add4)+(0,-0.75)$) {$n_{\timeVar}$};
	\node[sysBlock] at ($(add4)+(-3,0)$) (comchan) {$\mybold{h}_{\timeVar}$};
	\node (distant) at (in1 |- add4) {$d_{\timeVar}$};
	\node (out) at ($(add3)+(0,-1)$) {$y_{\timeVar}$};

	\node (labelSi) at ($(dot2) + (2,1.5)$) {Nonlinear cascade SI channel};
	\begin{pgfonlayer}{background}
		\node [fill=red!10,fit={($(labelSi.west) + (0,0.1)$) ($(sichan) + (0.55,-0.9)$)}] {};
	\end{pgfonlayer}

	\node (labelWc) at ($(comchan) + (-0.8,0.6)$) {Wireless channel};

	\draw[dashed] ($(labelWc.west) + (0,0.25)$) rectangle ($(comchan) + (1,-1)$);

	\draw[thick,->] (add1) -| (sichan);
	\draw[thick,->] (sichan) to node[right,pos=0.8] {$x_{\text{si},\timeVar}$} (add3);

	\draw[thick,->] (noise) to (add4);

	\draw[thick,->] (add4) to node[above] {$s_{\timeVar}$} (add3);
	\draw[thick,->] (distant) to (comchan);
	\draw[thick,->] (comchan) to node[above] {$d^{h}_{\timeVar}$} (add4);
	\draw[thick,->] (add3) to (out);
	\end{tikzpicture}
	\caption{Cascade model of the nonlinear SI channel with input $x_{\timeVar}$ and the signal-of-interest $d_{\timeVar}$.}
	\label{fig:systemmodel}
\end{figure}
The local transmitter creates an input signal $x_{\timeVar}$ to be conveyed to a distant receiver. Due to hardware impairments like I/Q imbalance in the quadrature mixer or high-order harmonics of the power amplifier, the signal is distorted by nonlinear components. This effect is especially severe at high transmitter output powers~\cite{src:korpi2014full}. We model the nonlinear contribution as $N^{\text{th}}$-order memoryless expansion, where the $i^{\text{th}}$ component consists of a nonlinear basis function $\phi_i\left(\cdot\right)$ and the coefficient $a_{i,\timeVar}$. Now, the linear SI channel is modeled by a linear, time-variant finite-impulse response (FIR) filter.  Let the aforementioned FIR filter have $\firLengthSI$ leading coefficients given by
\begin{align}
\label{eq:siChanDef}
\mybold{w}_{\timeVar}=\left[ w_{\timeVar,0}, w_{\timeVar,1}, \ldots, w_{\timeVar,\firLengthSI-1}\right]^T.
\end{align}
Alternatively, the nonlinear SI channel can be modeled in parallel structure, which is shown in~\figref{fig:parallelmodel}.
\begin{figure}
	\centering
	\begin{tikzpicture}[
	nlBlock/.style={draw, drop shadow, fill=white, rectangle, minimum height=1.75em, minimum width=4.5em},
	chanBlock/.style={draw, drop shadow, fill=white, rectangle, minimum height=1.75em, minimum width=4.5em}
	]

	\node (in1) {$x_{\timeVar}$};
	\node[dot] (dot2) at ($(in1)+(1,0)$) {};

	\node[nlBlock] at ($(dot2)+(1.35,0.8)$) (nlfunc0) {$\phi_0(\cdot)$};
	\node[nlBlock] at ($(dot2)+(1.35,0)$) (nlfunc1) {$\phi_1(\cdot)$};
	\node at ($(dot2)+(1.35,-0.6125)+(0,0)$) (colon1) {$\colon$};
	\node[nlBlock] at ($(dot2)+(1.35,-1.25)$) (nlfuncN) {$\phi_{N-1}(\cdot)$};

	\node[chanBlock] at ($(nlfunc0)+(2.5,0)$) (sichan0) {$\mybold{w}_{0,\timeVar}$};
	\node[chanBlock] at ($(nlfunc1)+(2.5,0)$) (sichan1) {$\mybold{w}_{1,\timeVar}$};
	\node at ($(colon1)+(2.5,0)+(0,0)$) (colon1) {$\colon$};
	\node[chanBlock] at ($(nlfuncN)+(2.5,0)$) (sichanN) {$\mybold{w}_{N-1,\timeVar}$};

	\myadd{add1}{$(sichan1)+(1.5,0)$}

	\draw[thick,->] (in1) to (dot2);

	\draw[thick,->] (dot2) |- (nlfunc0);
	\draw[thick,->] (dot2) -- (nlfunc1);
	\draw[thick,->] (dot2) |- (nlfuncN);
	\draw[thick,->] (nlfunc0) -- (sichan0);
	\draw[thick,->] (nlfunc1) -- (sichan1);
	\draw[thick,->] (nlfuncN) -- (sichanN);
	\draw[thick,->] (sichan0) -| (add1);
	\draw[thick,->] (sichan1) -- (add1);
	\draw[thick,->] (sichanN) -| (add1);

	\node (out) at ($(add1)+(1,0)$) {$x_{\text{si},\timeVar}$};
	\draw[thick,->] (add1) to (out);
	\end{tikzpicture}
	\caption{Parallel model of the nonlinear SI channel.}
	\label{fig:parallelmodel}
\end{figure}
In that case, each of the nonlinear basis functions is followed by a different FIR filter with coefficients $\mybold{w}_{i,\timeVar}$ for $0\leq i\leq N-1$. However, the parallel approach significantly increases the degrees of freedom to be estimated. Conventional SI channel models require a large number of filter taps, which has been identified as an obstacle for practical systems~\cite{src:nadh2017ataylor}. Specifically, such an architecture might lead to ambiguous solutions. Thus, we focus on the cascade modeling of~\figref{fig:systemmodel} in our work. Then, the SI signal of~\figref{fig:systemmodel} can be written in a generalized linear model~\cite{src:bishop2006pattern}
\begin{align}
x_{\text{si},\timeVar}
&= \sum^{\firLengthSI-1}_{l=0} w_{\timeVar,l} \sum_{i=0}^{N-1}a_{i,\timeVar}\phi_i\left(x_{\timeVar-l}\right), \nonumber\\
\label{eq:nldef}
&= \sum_{i=0}^{N-1}a_{i,\timeVar}\sum^{\firLengthSI-1}_{l=0} w_{\timeVar,l} \phi_i\left(x_{\timeVar-l}\right),
\end{align}
where $w_{\timeVar,l}$ is taken from~\eqref{eq:siChanDef}.

At a distant node, the transmitter broadcasts a SoI, denoted by $d_{\timeVar}$ with zero mean and power $\expv\left[|d_{\timeVar}|^2\right]=P_d$ over a wireless channel. Similar to the SI channel, the wireless channel is modeled as linear FIR filter with $\firLengthWC$ coefficients
\begin{align}
\label{eq:wChanDef}
\mybold{h}_{\timeVar}=\left[ h_{\timeVar,0}, h_{\timeVar,1}, \ldots, h_{\timeVar,\firLengthWC-1}\right]^T
\end{align}
with the same number of coefficients $L$ as in the case of the linear SI channel for convenience. Hence, at the receiver, the SoI is designated as
\begin{align}
\label{eq:rxSoi}
d^{h}_{\timeVar} = \sum_{l=0}^{\firLengthWC-1}h_{\timeVar,l}d_{\timeVar-l}.
\end{align}
The signal $s_{\timeVar}$ contains all elements at the receiver that are not related to the SI, such as the SoI~\eqref{eq:rxSoi} and the additive noise $n_{\timeVar}$, i.e.,
\begin{align}
s_{\timeVar} &= d^{h}_{\timeVar} + n_{\timeVar} \notag \\
\label{eq:notSIrelatedTime}
&= \sum_{l=0}^{\firLengthWC-1}h_{\timeVar,l}d_{\timeVar-l} + n_{\timeVar}
\end{align}
and hence we have the overall observed signal at the receiver
\begin{align}
y_{\timeVar} = x_{\text{si},\timeVar} + s_{\timeVar}.
\end{align}

Next, we introduce a vector notation for all signals in time domain. We group consecutive samples in time into frames. This is depicted in \figref{fig:frames}. Suppose that there are frames consisting of $M$ samples, where each frame is indexed by $\frameVar$. We form a new frame on every $R^{\text{th}}$ sample, such that we have

an overlap of neighboring frames by $\firLengthSI=M-R$ samples.
\begin{figure}
	\centering
	\newcommand*{\timeBlockLength}{1.75}
	\begin{tikzpicture}
	[
	timeBlock/.style={draw, drop shadow, fill=white, rectangle, minimum height=2em, minimum width=\timeBlockLength cm},
	dimen/.style={<->,>=latex,thin,every rectangle node/.style={fill=white,midway}},
	dimenright/.style={->,>=latex,thin,every rectangle node/.style={fill=white,right}}
	]
	\pgfmathsetmacro{\dotsGap}{1.25}
	\node[timeBlock] (t1) at (0,0) {\scriptsize $\frameVar R-M+1$};

	\node[] (c1) at ($(t1)+(\dotsGap,0)$) {$\dots$};
	\node[timeBlock] (t2) at ($(c1)+(\dotsGap,0)$) {\scriptsize $(\frameVar-1)R$};
	\node[timeBlock] (t3) at ($(t2)+(\timeBlockLength,0)$) {\scriptsize $(\frameVar-1)R+1$};
	\node[] (c2) at ($(t3)+(\dotsGap,0)$) {$\dots$};
	\node[timeBlock] (t4) at ($(c2)+(\dotsGap,0)$) {\scriptsize $\frameVar R$};

	\draw (t3.south west) -- ++(0,-0.5) coordinate (D1) -- +(0,-5pt);
	\draw (t4.south east) -- ++(0,-0.5) coordinate (D2) -- +(0,-5pt);
	\draw [dimen] (D1) -- (D2) node {$R$};
	\draw (t1.south west) -- ++(0,-1) coordinate (D3) -- +(0,-5pt);
	\draw (t4.south east) -- ++(0,-1) coordinate (D4) -- +(0,-5pt);
	\draw [dimen] (D3) -- (D4) node {$M$};

	\draw (t1.north west) -- ++(0,0.25) coordinate (time1) -- +(0,5pt);
	\draw [dimenright] (time1) -- ($(time1)+(1,0)$) node {$\timeVar$};

	\end{tikzpicture}
	\caption{Frame of $M$ samples with frame shift $R$ at index $\frameVar$.}
	\label{fig:frames}
\end{figure}
On this basis, we can express the convolution of the transmitted signal and the SI channel by the overlap-save method. Let

\begin{align}
\label{eq:nlBasisVec}
\mybold{\phi}_{i,\frameVar} = \left[ \phi_{i}\left(x_{\frameVar R-M+1}\right),\phi_{i}\left(x_{\frameVar R-M+2}\right),\ldots,\phi_{i}\left(x_{\frameVar R}\right) \right]^T
\end{align}
be the $\frameVar^{\text{th}}$ frame vector for the $i^{\text{th}}$ nonlinear basis function, applied to the input $x_{\timeVar}$.
Furthermore, let
\begin{align}
\label{eq:nlBasisMatrix}
\mybold{\Phi}_{\frameVar}=\left[ \mybold{\phi}_{0,\frameVar}, \mybold{\phi}_{1,\frameVar},\ldots, \mybold{\phi}_{N-1,\frameVar} \right]
\end{align}
be a $M\times N$ matrix with the signal of all nonlinear basis functions from~\eqref{eq:nlBasisVec}, stacked as column vectors.

We assume the SI channel to be slowly varying and therefore constant within a single frame, similar to block fading. Then, we form a FIR channel vector of size $M\times 1$ by appending $R$ zeros to~\eqref{eq:siChanDef}, and get
\begin{align}
\label{eq:siChanVec}
\mybold{w}_{\frameVar} = \left[ \mybold{w}^T_{\frameVar R}, \mybold{0}^T_{R\times 1} \right]^T.
\end{align}

Similarly, we define a $M\times 1$ frame vector for the desired $d_{\timeVar}$
\begin{align}
\mybold{d}_{\frameVar} &= \left[ d_{\frameVar R-M+1},d_{\frameVar R-M+2},\ldots,d_{\frameVar R} \right]^T

\end{align}
and a $R\times 1$ frame vector for the receiver noise
\begin{align}
\mybold{n}_{\frameVar}=\left[n_{(\frameVar-1)R+1},n_{(\frameVar-1)R+2},\ldots,n_{\frameVar R} \right]
\end{align}
and we define $\mybold{y}_{\frameVar}$, $\mybold{s}_{\frameVar}$ and $\mybold{x}_{\text{si},\frameVar}$ similarly as $\mybold{n}_{\frameVar}$.

\subsection{DFT-domain representation}
\label{sec:dft}
Next, we transform selected signals into DFT domain and express the convolutions of~\eqref{eq:nldef} and~\eqref{eq:rxSoi} as multiplications.
Let
\begin{align}
\underline{\mybold{\Phi}}_{\frameVar} &=\mybold{F}_M\mybold{\Phi}_{\frameVar}

 = \left[ \underline{\mybold{\phi}}_{0,\frameVar},
\label{eq:basisFrameMatrix}
\underline{\mybold{\phi}}_{1,\frameVar},\ldots, \underline{\mybold{\phi}}_{N-1,\frameVar} \right]
\end{align}
denote the $M\times N$ DFT-domain representations of the basis frame, where $\mybold{F}_M$ is the $M\times M$ DFT matrix. Similarly, let
\begin{align}
\label{eq:siChanVecDft}
\underline{\mybold{w}}_{\frameVar} &=\mybold{F}_M\mybold{w}_{\frameVar}
\end{align}
be the DFT-domain transform of the linear SI channel coefficients from~\eqref{eq:siChanVec}. It is likewise assumed that the nonlinear coefficients $a_{i,\timeVar}$ only vary slowly over time, thus they are approximately constant within a frame. We write
\begin{align}
	\label{eq:nlCoeffsVecEqual}
	\underline{a}_{i,\frameVar} = a_{i,(\frameVar-1)R} = a_{i,(\frameVar-1)R+1} = \ldots = a_{i,\frameVar R}
\end{align}
and group the nonlinear coefficients~\eqref{eq:nlCoeffsVecEqual} in vector form
\begin{align}
	\label{eq:nlCoeffsVec}
\underline{\mybold{a}}_{\frameVar}= \left[\underline{a}_{0,\frameVar},\underline{a}_{1,\frameVar},\ldots,\underline{a}_{N-1,\frameVar}\right]^T.
\end{align}
Now we are equipped to express the SI signal of~\eqref{eq:nldef} in DFT domain. The nonlinear basis frames are given by $\underline{\mybold{\Phi}}_{\frameVar}\underline{\mybold{a}}_{\frameVar}$ based on~\Cref{eq:basisFrameMatrix,eq:nlCoeffsVec}. They are multiplied with the linear SI channel coefficients~\eqref{eq:siChanVecDft}, since the convolution is transformed into a multiplication in DFT domain. Thus, we construct
\begin{align}
\label{eq:obsCompact}
\underline{\mybold{x}}_{\text{si},\frameVar}
&=  \diag\left[\underline{\mybold{\Phi}}_{\frameVar}\underline{\mybold{a}}_{\frameVar}\right] \underline{\mybold{w}}_{\frameVar} \\
\label{eq:obsDetail}
&\hintedrel[eq:xsiDFT1]{=} \left(\sum_{i=0}^{N-1}\underline{a}_{i,\frameVar}\diag\left[ \underline{\mybold{\phi}}_{i,\frameVar}\right]\right) \underline{\mybold{w}}_{\frameVar},
\end{align}
where we introduce~$(\hintref{eq:xsiDFT1})$ to prepare a form that later on allows to efficiently infer the coefficients $\underline{a}_{i,\frameVar}$ by a linear algorithm.
The corresponding $R\times 1$ time-domain signal of~\eqref{eq:obsDetail}, in line with~\eqref{eq:nldef}, is
\begin{align}
\mybold{x}_{\text{si},\frameVar} &= \mybold{\Upsilon}^T\mybold{F}^{-1}_M\underline{\mybold{x}}_{\text{si},\frameVar} \\
\label{eq:siTime}
&=  \left(\sum_{i=0}^{N-1}\underline{a}_{i,\frameVar}\mybold{\Upsilon}^T\mybold{F}^{-1}_M\diag\left[ \underline{\mybold{\phi}}_{i,\frameVar}\right]\right) \underline{\mybold{w}}_{\frameVar},
\end{align}
where the matrix
\begin{align}
\mybold{\Upsilon}^T =
\begin{bmatrix}
    \zeros{R}{\firLengthSI} & \eye{R}
\end{bmatrix}
\end{align}
removes the first $\firLengthSI=M-R$ entries from a vector since these entries are contaminated by aliasing. Next, we address the SoI. In DFT domain, we have
\begin{align}
\label{eq:soiDFT}
\underline{\mybold{d}}_{\frameVar}&=\mybold{F}_M\mybold{d}_{\frameVar} \\
\label{eq:wChanDFT}
\underline{\mybold{h}}_{\frameVar}&=\mybold{F}_M\mybold{h}_{\frameVar}
\end{align}
for the SoI and the coefficients of the wireless channel.

Similar to~\eqref{eq:siTime}, we express the convolution of SoI and the wireless channel as multiplication of~\eqref{eq:soiDFT} and~\eqref{eq:wChanDFT} in DFT domain, and then transform it back into time domain
\begin{align}
\label{eq:soiPlusNoiseTime}
\mybold{s}_{\frameVar} &=  \mybold{\Upsilon}^T\mybold{F}^{-1}_M\diag\left[ \underline{\mybold{d}}_{\frameVar}\right] \underline{\mybold{h}}_{\frameVar}+\mybold{n}_{\frameVar}.
\end{align}

Next, we transform the overall received signal $\mybold{y}_{\frameVar} = \mybold{x}_{\text{si},\frameVar} + \mybold{s}_{\frameVar}$ back into DFT domain by first prepending $\firLengthSI$~zeros and then applying the DFT matrix. Thus we have
\begin{align}
\underline{\mybold{y}}_{\frameVar}&= \mybold{F}_M\mybold{\Upsilon}\mybold{y}_{\frameVar} \nonumber\\
&= \mybold{F}_M\mybold{\Upsilon}\mybold{x}_{\text{si},\frameVar} + \mybold{F}_M\mybold{\Upsilon}\mybold{s}_{\frameVar} \\
\label{eq:obsDFT}
&\hintedrel[eq:obsDFT1]{=} \left(\sum_{i=0}^{N-1} \underline{a}_{i,\frameVar}\underline{\mybold{C}}_{i,\frameVar}\right) \underline{\mybold{w}}_{\frameVar} + \underline{\mybold{s}}_{\frameVar},
\end{align}
where~$(\hintref{eq:obsDFT1})$ uses~\eqref{eq:siTime} and we employed the definition
\begin{align}
\label{eq:CDef}
 \underline{\mybold{C}}_{i,\frameVar} = \mybold{F}_M\mybold{\Upsilon}\mybold{\Upsilon}^T\mybold{F}^{-1}_M\diag\left[ \underline{\mybold{\phi}}_{i,\frameVar}\right]
\end{align}
and the DFT domain representation
\begin{align}
\label{eq:soiPlusNoiseDft}
\underline{\mybold{s}}_{\frameVar}&=\mybold{F}_M\mybold{\Upsilon}\mybold{s}_{\frameVar}.
\end{align}
The received signal without SI contribution~\eqref{eq:soiPlusNoiseDft} can be alternatively expressed by using~\eqref{eq:soiPlusNoiseTime} by
\begin{align}
\underline{\mybold{s}}_{\frameVar}&=

\underline{\mybold{d}}^{h}_{\frameVar} + \underline{\mybold{n}}_{\frameVar},
\end{align}
where the term $\underline{\mybold{n}}_{\frameVar}$ denotes the DFT domain representation of the additive receiver noise and

$\underline{\mybold{d}}^{h}_{\frameVar}$ is the DFT domain equivalent of the SoI at the receiver.

\subsection{State-space model}
We propose state-space models for both the linear SI channel~\eqref{eq:siChanVecDft} and the nonlinear coefficients~\eqref{eq:nlCoeffsVec}.

The time-variant nature of the nonlinear SI channel is specifically modeled by first-order Markov models of the linear SI channel and the nonlinear coefficients, i.e.,
\begin{align}
\label{eq:siChanMarkov}
\underline{\mybold{w}}_{\frameVar} &=
\underline{A}^{w}\underline{\mybold{w}}_{\frameVar-1}+\underline{\mybold{\Delta}}^{w}_{\frameVar},
 \\
\label{eq:nlCoeffsMarkov}
\underline{\mybold{a}}_{\frameVar} &=
\underline{\mybold{A}}^{a}\underline{\mybold{a}}_{\frameVar-1}+\underline{\mybold{\Delta}}^{a}_{\frameVar},
\end{align}
where the parameter $\underline{A}^{w}$ with $0\leq|\underline{A}^{w}|\leq 1$ and the diagonal matrix $\underline{\mybold{A}}^{a}$ denote the transition factor between consecutive channel realizations and nonlinear coefficients over time, respectively. We assume a scalar $\underline{A}^{w}$, since the linear SI channel coefficients are assumed to be similar in transient behavior. On the other hand, the matrix $\underline{\mybold{A}}^{a}$ reflects different degrees of temporal variations among the nonlinear coefficients. The Gaussian system noise variables $\underline{\mybold{\Delta}}^{w}_{\frameVar}\sim\cgauss{\mybold{0},\underline{\psi}^{\Delta w}_{\frameVar}\eye{M}}$ and $\underline{\mybold{\Delta}}^{a}_{\frameVar}\sim\cgauss{\mybold{0},\underline{\mybold{\Psi}}^{\Delta a}_{\frameVar}}$ are independent with variance $\underline{\psi}^{\Delta w}_{\frameVar}$ and covariance matrix $\underline{\mybold{\Psi}}^{\Delta a}_{\frameVar}$. By definition, we assume the top path from~\figref{fig:systemmodel} with index $i=0$ represents the linear component of the SI, thus, without loss of generality, we have $\phi_{0}(x_{\timeVar})=x_{\timeVar}$ and fix $\underline{a}_{0,\frameVar}=1$, $\underline{\Delta}^{a}_{0,\frameVar}=0$ in the following.

Let $\underline{\mybold{R}}^{w}_{\frameVar}=\expv\left[ \underline{\mybold{w}}_{\frameVar}\underline{\mybold{w}}^{H}_{\frameVar} \right]$ be the autocorrelation matrix of the random SI channel coefficients,

thus, from~\eqref{eq:siChanMarkov}, we have
\begin{align}
\label{eq:siChanCorr}
\tr\left[\underline{\mybold{R}}^{w}_{\frameVar}\right]=\left|\underline{A}^{w}\right|^{2}\tr\left[\underline{\mybold{R}}^{w}_{\frameVar-1}\right]+M\underline{\psi}^{\Delta w}_{\frameVar}.
\end{align}

Apparently, the system model is characterized by a significant number of parameters. Since these parameters are difficult to determine, we reduce their number by using certain relations between them. For instance, from a practical perspective, it is reasonable to assume that the statistical properties of the linear SI channel remain similar over time. Thus, we assume
\begin{align}
\label{eq:siChanCorrTrEqual}
\tr\left[\underline{\mybold{R}}^{w}_{\frameVar}\right] &= \tr\left[\underline{\mybold{R}}^{w}_{\frameVar+1}\right] = \ldots = \tr\left[\underline{\mybold{R}}^{w}_{\infty}\right].

\end{align}

and thus we find the covariance of the system noise as
 \begin{align}
 \label{eq:siChanSystemNoiseCovLongTerm}
\underline{\psi}^{\Delta w} =
\frac{1}{M}\tr\left[\underline{\mybold{R}}^{w}_{\infty}\right]
 \left(1-\left|\underline{A}^{w}\right|^{2}\right).
 \end{align}
The steady-state parameter $\tr\left[\underline{\mybold{R}}^{w}_{\infty}\right]$ determines the power of the linear SI channel components. In practice, this can be obtained from a-priori knowledge, since the relation of SI power level to the SoI is at least approximately known.

Furthermore, the \textit{coherence time} reflects the variability of the communication channel in time, and therefore is regarded as main characteristic of temporal variations. To best of our knowledge, the coherence time has not been systematically studied for SI channels so far. However, we believe it to be on the same order as the wireless channel, as indicated in some previous work~\cite{src:jain2011practical}. To integrate the notion of a channel coherence time into our framework, we introduce a coherence frame index $\frameVar^{w}_{\text{coh}}$, which denotes the time when the correlation of subsequent channel variables has dropped by half, thus
\begin{align}
\label{eq:channelCoherence}
\left| A^{w} \right|^{\frameVar^{w}_{\text{coh}}} = \frac{1}{2}.
\end{align}

For the zero-mean nonlinear coefficients $\underline{\mybold{a}}_{\frameVar}$, we define the diagonal covariance matrix $\underline{\mybold{R}}^{a}_{\frameVar}= \expv\left[ \underline{\mybold{a}}_{\frameVar}\underline{\mybold{a}}^{H}_{\frameVar} \right]$ with the $i^{\text{th}}$ diagonal element $\underline{p}^{a}_{i,\frameVar}=\expv\left[ \left|\underline{a}_{i,\frameVar}\right|^2 \right]$. From~\eqref{eq:nlCoeffsMarkov}, we can derive
\begin{align}
\label{eq:nlCorr}
\underline{\mybold{R}}^{a}_{\frameVar}=\underline{\mybold{A}}^{a}\underline{\mybold{R}}^{a}_{\frameVar-1}\underline{\mybold{A}}^{a^H}+\underline{\mybold{\Psi}}^{\Delta a}_{\frameVar}.
\end{align}

With the elements of $\underline{\mybold{R}}^{a}_{\frameVar}$ approximately constant over time,
\begin{align}
\underline{\mybold{R}}^{a}_{\frameVar}=\underline{\mybold{R}}^{a}_{\frameVar+1}=\ldots=\underline{\mybold{R}}^{a}_{\infty},
\end{align}
the system noise covariance matrix is expressed by
\begin{align}

\underline{\mybold{\Psi}}^{\Delta a}_{\frameVar} =
\underline{\mybold{R}}^{a}_{\infty}
 \left(\mybold{I}- \underline{\mybold{A}}^{a}\underline{\mybold{A}}^{a^H}\right),
\end{align}
since $\underline{\mybold{A}}^{a}$ is diagonal.
Similar to~\eqref{eq:channelCoherence}, we define a coherence time $\frameVar^{a}_{i,\text{coh}}$ for the $i^{\text{th}}$ nonlinear coefficient by
\begin{align}
\label{eq:nlCoherence}
\left| A_i^{a} \right|^{\frameVar^{a}_{i,\text{coh}}} = \frac{1}{2},
\end{align}
where $A_i^{a}$ denotes the $i^{\text{th}}$ element on the main diagonal of $\underline{\mybold{A}}^{a}$.

\section{State-space nonlinear adaptive algorithm}
\label{sec:algorithm}
Consider the proposed SI estimation and cancellation algorithm depicted in~\figref{fig:algorithm}.
\begin{figure}
	\begin{tikzpicture}[
	sysBlock/.style={draw, fill=white, rectangle, minimum height=2em, minimum width=3.5em},
	nlBlock/.style={draw, fill=white, rectangle, minimum height=2em, minimum width=4em},
	shadowBlock/.style={draw, drop shadow, fill=white, rectangle},
	triangle/.style = {draw, regular polygon, regular polygon sides=3 },

	border rotated/.style = {shape border rotate=90}
	]

	\node (inTx) {$\underline{\mybold{x}}_{\frameVar}$};

	\node[shadowBlock, align=center, minimum height=7em] (kalpred) at ($(inTx) +(1,-1.8)$) {Kalman\\prediction\\$\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}$,\\$\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar-1}$};
	\node (inRx) at ($(kalpred)+(-1,-2)$) {$\underline{\mybold{y}}_{\frameVar}$};

	\node[shadowBlock, align=center, text width=1.3cm, minimum height=7em] (sirec) at ($(kalpred) + (2,0)$) {SI\\recon\-struc\-tion};

	\node[shadowBlock, align=center, text width=2cm] (ortho) at ($(inTx)+(1.8,0)$) {Ortho\-go\-n\-ali\-za\-tion};
	\myadd{add1}{$(inRx-|sirec)$}
	\node  at ($(add1.east)+(0.1,0.2)$) {$-$};

	\node[shadowBlock, align=center, minimum height=7em] (dec) at ($(sirec)+(1.75,0)$) {Decoder};
	\myadd{add2}{$(dec)+(1.5,0)$}
	\node  at ($(add2.west)+(-0.1,0.2)$) {$-$};

	\node[shadowBlock, align=center, minimum height=7em] (kalupd) at ($(add2)+(1.5,0)$) {Kalman\\update\\$\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}$,\\$\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar}$};

	\node[dot] (dot1) at ($(sirec|-ortho)+(0.25,0)$) {};
	\node[dot] (dot2) at ($(dec|-inRx)$) {};

	\draw[thick,->] (inTx) to (ortho);
	\draw[thick,->] (ortho) -- (dot1);
	\draw[thick,->] (dot1) -- (sirec.north-|dot1);
	\draw[thick,->] (dot1) -| (kalupd);

	\draw[thick,->] (kalpred) -- (sirec);
	\draw[thick,->] (inRx) -- (add1);
	\draw[thick,->] (sirec) -- node[left] {$\underline{\hat{\mybold{x}}}_{\text{si},\frameVar}$} (add1);
	\draw[thick,->] (add1) -- node[above] {$\underline{\mybold{e}}_{\frameVar}$} (dot2);
	\draw[thick,->] (dot2) -- (dec);
	\draw[thick,->] (dot2) -| (add2);

	\draw[thick,->] (dec) -- node[below] {$\underline{\hat{\mybold{d}}}^h_{\frameVar}$} (add2);
	\draw[thick,->] (add2) -- node[above] {$\underline{\tilde{\mybold{e}}}_{\frameVar}$} (kalupd);
\end{tikzpicture}
\caption{Proposed nonlinear SI estimation and cancellation.}
\label{fig:algorithm}
\end{figure}
All signals are represented here in DFT domain. The adaptive algorithm has access to the transmitted signal $\underline{\mybold{x}}_{\frameVar}$. Initially, the input is subject to an orthogonalization process, where the nonlinear contributions are transformed to uncorrelated equivalents. Then, the adaptation is performed in an iterative process over time. At each step $\frameVar$, we first acquire a predicted Kalman estimate of the linear SI channel and the nonlinear coefficients.

We intend to cancel as much SI as possible, thus the algorithm produces a reconstruction $\underline{\hat{\mybold{x}}}_{\text{si},\frameVar}$ of the SI signal. After subtracting~$\underline{\hat{\mybold{x}}}_{\text{si},\frameVar}$  from the received signal $\underline{\mybold{y}}_{\frameVar}$, the algorithm outputs an error signal $\underline{\mybold{e}}_{\frameVar}$. This error signal is then provided to the decoder which produces a reconstructed SoI $\underline{\hat{\mybold{d}}}^h_{\frameVar}$. Finally, after $\underline{\hat{\mybold{d}}}^h_{\frameVar}$ is subtracted from the error signal $\underline{\mybold{e}}_{\frameVar}$, the residual error $\underline{\tilde{\mybold{e}}}_{\frameVar}$ supports the Kalman update estimations for both the linear SI channel and the nonlinear coefficients.

\subsection{Orthogonalization}
\label{sec:ortho}
The nonlinear basis functions $\phi_i(\cdot)$ do not necessarily represent an orthogonal basis for the SI signal. Therefore, we intend to change the basis by a matrix transformation. There are many approaches available that achieve this orthogonalization, such as eigenvalue decomposition~\cite{src:korpi2015adaptive}, Cholesky factorization~\cite{src:emara2017nonlinear} or Gram-Schmidt method~\cite{src:proakis2007digital}.

We describe the statistical connection between different basis functions by a time-invariant $N\times N$ autocorrelation matrix. Using~\eqref{eq:basisFrameMatrix}, we have

\begin{align*}
	\mybold{R}^{\Phi} &= \expv\Big[ \underline{\mybold{\Phi}}^T_{\frameVar} \underline{\mybold{\Phi}}^*_{\frameVar}\Big].
\end{align*}

Next, we choose a $N\times N$ transform matrix $\mybold{G}$, which is not necessarily unitary, but full-rank, and such that $\mybold{G}\mybold{R}^{\Phi}\mybold{G}^H$ is diagonal.

We apply the transformation to each basis frame by
\begin{align}
	\underline{\tilde{\mybold{\Phi}}}_{\frameVar} &= \underline{\mybold{\Phi}}_{\frameVar} \mybold{G}^T \nonumber\\
	\label{eq:orthoPhi}
	&= \left[ \underline{\tilde{\mybold{\phi}}}_{0,\frameVar}, \underline{\tilde{\mybold{\phi}}}_{1,\frameVar},\ldots, \underline{\tilde{\mybold{\phi}}}_{N-1,\frameVar} \right].
\end{align}
Next, starting from~\eqref{eq:obsCompact}, we have
\begin{align}
\underline{\mybold{x}}_{\text{si},\frameVar}
&=  \diag\left[\underline{\mybold{\Phi}}_{\frameVar}\mybold{G}^{T}\left(\mybold{G}^T\right)^{-1}\underline{\mybold{a}}_{\frameVar}\right] \underline{\mybold{w}}_{\frameVar} \nonumber\\
&\hintedrel[eq:obsOrtho1]{=} \diag\left[\underline{\tilde{\mybold{\Phi}}}_{\frameVar}\underline{\tilde{\mybold{a}}}_{\frameVar}\right]\underline{\mybold{w}}_{\frameVar} \nonumber\\
&= \diag\left[\underline{\tilde{\mybold{\Phi}}}_{\frameVar}\underline{\tilde{\mybold{a}}}_{\frameVar}\frac{1}{\underline{\tilde{a}}_{0,\frameVar}}\right]\underline{\tilde{a}}_{0,\frameVar}\underline{\mybold{w}}_{\frameVar} \nonumber\\
&\hintedrel[eq:obsOrtho2]{=}  \left(\sum_{i=0}^{N-1}\underline{\check{a}}_{i,\frameVar}\diag\left[ \underline{\tilde{\mybold{\phi}}}_{i,\frameVar}\right]\right) \underline{\tilde{\mybold{w}}}_{\frameVar},
\end{align}
where~$(\hintref{eq:obsOrtho1})$ uses~\eqref{eq:orthoPhi} and $\underline{\tilde{\mybold{a}}}_{\frameVar}=\left(\mybold{G}^T\right)^{-1}\underline{\mybold{a}}_{\frameVar}$, and~$(\hintref{eq:obsOrtho2})$ applies normalization with $\underline{\check{a}}_{i,\frameVar}=\underline{\tilde{a}}_{i,\frameVar}/\underline{\tilde{a}}_{0,\frameVar}$ and $\underline{\tilde{\mybold{w}}}_{\frameVar}=\underline{\tilde{a}}_{0,\frameVar}\underline{\mybold{w}}_{\frameVar}$. Apparently, the change of basis has an effect on the outcome of both the estimated linear SI channel and the nonlinear coefficients as they represent the inner state of the system. However, the overall SI reconstruction is left unchanged. In the following, for the sake of simplicity, we drop the ``tilde'', but orthogonalization is always applied unless noted otherwise.

\subsection{Predictions}
\label{sec:predictions}
Consider the state equations for the linear SI channel~\eqref{eq:siChanMarkov} and the nonlinear coefficients~\eqref{eq:nlCoeffsMarkov}. We intend to derive estimators for $\underline{\mybold{w}}_{\frameVar}$ and $\underline{\mybold{a}}_{\frameVar}$. Let
\begin{align}
\label{eq:siChanEstimatorPred}
\underline{\mybold{w}}_{\frameVar} &= \hat{\underline{\mybold{w}}}_{\frameVar|\frameVar-1} + \underline{\mybold{n}}^w_{\frameVar|\frameVar-1}\\
\label{eq:nlCoeffsEstimatorPred}
\underline{\mybold{a}}_{\frameVar} &= \hat{\underline{\mybold{a}}}_{\frameVar|\frameVar-1} + \underline{\mybold{n}}^a_{\frameVar|\frameVar-1},
\end{align}
where $\hat{\underline{\mybold{w}}}_{\frameVar|\frameVar-1}$ and $\hat{\underline{\mybold{a}}}_{\frameVar|\frameVar-1}$ are estimators for the linear SI channel and the nonlinear coefficients given all past observations $\left(\underline{\mybold{y}}_{0},\ldots,  \underline{\mybold{y}}_{\frameVar-2}, \underline{\mybold{y}}_{\frameVar-1}\right)$, respectively. The zero-mean a-priori errors $\underline{\mybold{n}}^w_{\frameVar|\frameVar-1}$ and $\underline{\mybold{n}}^a_{\frameVar|\frameVar-1}$ have the covariances $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}$ and $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$, respectively.

Similarly, we define $\hat{\underline{\mybold{w}}}_{\frameVar|\frameVar}$ and $\hat{\underline{\mybold{a}}}_{\frameVar|\frameVar}$ to be estimators for the linear SI channel and the nonlinear coefficients given the present and all past observations $\left(\underline{\mybold{y}}_{0},\ldots,  \underline{\mybold{y}}_{\frameVar-1}, \underline{\mybold{y}}_{\frameVar}\right)$, respectively. The zero-mean a-posteriori errors
have the covariances $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}$ and $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar}$, respectively. Furthermore, consider the following assumptions:
\begin{enumerate}
\item The linear SI channel $\underline{\mybold{w}}_{\frameVar}$ and the nonlinear coefficients $\underline{\mybold{a}}_{\frameVar}$ are independent.
\item The joint probability density functions of both the estimators
$\hat{\underline{\mybold{w}}}_{\frameVar|\frameVar-1}$,
$\hat{\underline{\mybold{a}}}_{\frameVar|\frameVar-1}$,
$\hat{\underline{\mybold{w}}}_{\frameVar|\frameVar}$,
$\hat{\underline{\mybold{a}}}_{\frameVar|\frameVar}$ and the a-priori, a-posterori errors

are circular-symmetric, complex Gaussian distributions.
\item The a-priori, a-posteriori errors

are independent.
\end{enumerate}
The adaptive algorithms for the estimation of the linear SI channel and of the nonlinear coefficients are derived in the following. We show the overall adaptive algorithm in~\algref{alg:overall}. Following the lines of~\cite{src:scharf1991statistical}, the predictions of the estimates $\hat{\underline{\mybold{w}}}_{\frameVar|\frameVar-1}$, $\hat{\underline{\mybold{a}}}_{\frameVar|\frameVar-1}$ and the estimation error covariances $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}$, $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$ are computed first. Next, the information encoded in the SoI is retrieved and its contribution is removed from the received signal. This is explained in more detail in Subsection~\ref{sec:decoding}. Finally, the update steps of the linear SI channel and the nonlinear coefficients are computed. This is presented in Subsections~\ref{sec:estLin} and~\ref{sec:estNl}, respectively.
\begin{algorithm}
\caption{Iterative estimation of linear SI channel~$\underline{{\mybold{w}}}_{\frameVar}$ and nonlinear coefficients~$\underline{\mybold{a}}_{\frameVar}$}
\label{alg:overall}
\begin{algorithmic}
\REQUIRE Initialization of $\underline{\hat{\mybold{w}}}_{0|0}$, $\underline{\mybold{P}}^{w}_{0|0}$, $\underline{\hat{\mybold{a}}}_{0|0}$ and $\underline{\mybold{P}}^{a}_{0|0}$
\REPEAT

\STATE $\frameVar \leftarrow \frameVar + 1$
\PREDICT{linear SI channel estimate}
\STATE $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}=\underline{A}^{w} \underline{\hat{\mybold{w}}}_{\frameVar-1|\frameVar-1}$
\ENDBLOCK
\PREDICT{linear SI channel estimation error covariance}
\STATE
$\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1} = \left|\underline{A}^{w}\right|^{2}\underline{\mybold{P}}^{w}_{\frameVar-1|\frameVar-1}

+\underline{\mybold{\Psi}}^{\Delta w}_{\frameVar}$
\ENDBLOCK
\PREDICT{coefficients estimate}
\STATE $\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar-1}=\underline{\mybold{A}}^{a} \underline{\hat{\mybold{a}}}_{\frameVar-1|\frameVar-1}$
\ENDBLOCK
\PREDICT{coefficients estimation error covariance}
\STATE
$\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1} = \underline{\mybold{A}}^{a}\underline{\mybold{P}}^{a}_{\frameVar-1|\frameVar-1}\underline{\mybold{A}}^{a^H}
+\underline{\mybold{\Psi}}^{\Delta a}_{\frameVar}$
\ENDBLOCK
\STATE{\textbf{Decode}~the information from the SoI and subtract its contribution from the received signal, see Section~\ref{sec:decoding}}
\UPDATE{linear SI channel estimate from prediction $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}$, see~\algref{alg:linearChannel} in Subsection~\ref{sec:estLin}}
\STATE Requires all predictions and provides $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}$ and $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}$
\ENDBLOCK
\UPDATE{nonlinear coefficients estimate from prediction $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$, see~\algref{alg:nonlinearCoeffs} in Subsection~\ref{sec:estNl}}
\STATE Requires all predictions, the updated $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}$ and $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}$, and provides $\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar}$ and $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar}$
\ENDBLOCK
\UNTIL{forever}
\end{algorithmic}
\end{algorithm}

\subsection{Decoding}
\label{sec:decoding}
After predicting the estimates of both the linear SI channel and the nonlinear coefficients, the SI contribution can be reconstructed based on~\eqref{eq:obsDFT} for the current time frame $\frameVar$ by
\begin{align}
\label{eq:reconstructedSI}
\underline{\hat{\mybold{x}}}_{\text{si},\frameVar}= \sum_{i=0}^{N-1}\underline{\hat{a}}_{i,\frameVar|\frameVar-1}\underline{\mybold{C}}_{i,\frameVar}\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}.
\end{align}
We define the error signal
\begin{align}
\label{eq:errorSignal}
\underline{\mybold{e}}_{\frameVar} &= \underline{\mybold{y}}_{\frameVar} - \underline{\hat{\mybold{x}}}_{\text{si},\frameVar}.
\end{align}
At this point the error signal~\eqref{eq:errorSignal} comprises many sources of observation noise, such as the residual SI, the SoI and independent additive noise. Large observation noise generally reduces the convergence speed or increases the misalignment of adaptive algorithms~\cite{Haykin2002}, but here we can reduce its contribution. Unlike in acoustic echo cancellation, where the acoustic signal is modeled as a random process, the SoI in wireless communication contains structure from encoded information. By decoding the information from the SoI and subtracting the corresponding SoI signal from~\eqref{eq:errorSignal}, we have better conditioning for subsequent updates of the SI path.

Let $\underline{\hat{\mybold{d}}}^{h}_{\frameVar}$ be the SoI generated from decoded information in DFT domain. Then,

we have the residual error
\begin{align}
\label{eq:residualError}
\underline{\tilde{\mybold{e}}}_{\frameVar} = \underline{\mybold{y}}_{\frameVar} - \underline{\hat{\mybold{x}}}_{\text{si},\frameVar} - \underline{\hat{\mybold{d}}}^{h}_{\frameVar}.
\end{align}

It comprises a potential decoding error, the residual SI and independent noise.

\subsection{Update estimation of linear SI channel}
\label{sec:estLin}
During the $\frameVar^{\text{th}}$~step of adaption, we need to know the estimate of the nonlinear coefficients from the previous step. Therefore, in the following, we assume that the estimates $\underline{\hat{a}}_{i,\frameVar|\frameVar-1}$ are given. Recall the observation equation~\eqref{eq:obsDFT}. After the decoding, we subtract the reconstructed SoI and have
\begin{align}
\underline{\tilde{\mybold{y}}}_{\frameVar} &= \underline{\mybold{y}}_{\frameVar} - \underline{\hat{\mybold{d}}}^{h}_{\frameVar} \notag \\
\label{eq:obsDec}
&\hintedrel[eq:y1]{=} \left(\sum_{i=0}^{N-1} \underline{a}_{i,\frameVar}\underline{\mybold{C}}_{i,\frameVar}\right) \underline{\mybold{w}}_{\frameVar} + \underline{\tilde{\mybold{s}}}_{\frameVar} \\
&\hintedrel[eq:y2]{=}
\left(\sum_{i=0}^{N-1} \underline{\hat{a}}_{i,\frameVar|\frameVar-1}\underline{\mybold{C}}_{i,\frameVar}\right) \underline{\mybold{w}}_{\frameVar} + \underline{\tilde{\mybold{s}}}^{w1}_{\frameVar},
\label{eq:obsLinear}
\end{align}
where~$(\hintref{eq:y1})$ introduces the noise term $\underline{\tilde{\mybold{s}}}_{\frameVar} = \underline{\mybold{s}}_{\frameVar} - \underline{\hat{\mybold{d}}}^{h}_{\frameVar}$

comprising the residual SI and additive noise with covariance matrix $\underline{\mybold{\Psi}}^{\tilde{s}}_{\frameVar}$. We assume that $\underline{\mybold{\Psi}}^{\tilde{s}}_{\frameVar}$ can be characterized from a-priori information on the noise level. The term~$(\hintref{eq:y2})$ is due to both~\eqref{eq:nlCoeffsEstimatorPred} and the zero-mean augmented noise term
\begin{align}
\label{eq:augNoiseLinear}
\underline{\tilde{\mybold{s}}}^{w1}_{\frameVar} = \left(\sum_{i=0}^{N-1} \underline{n}^{a}_{i,\frameVar|\frameVar-1}\underline{\mybold{C}}_{i,\frameVar}\right) \underline{\mybold{w}}_{\frameVar} + \underline{\tilde{\mybold{s}}}_{\frameVar}.
\end{align}
The covariance matrix of the augmented noise $\underline{\tilde{\mybold{s}}}^{w1}_{\frameVar}$ is
\begin{align}
\underline{\mybold{\Psi}}^{\tilde{s}^{w1}}_{\frameVar}
&= \expv\left[ \underline{\tilde{\mybold{s}}}^{w1}_{\frameVar}\underline{\tilde{\mybold{s}}}^{w1^H}_{\frameVar} \right] \nonumber \\

&=
\sum_{i=0}^{N-1}\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\underline{\mybold{C}}_{i,\frameVar}

\underline{\mybold{R}}^{w}_{\frameVar}
\underline{\mybold{C}}^{H}_{i,\frameVar}
+\underline{\mybold{\Psi}}^{\tilde{s}}_{\frameVar},

\label{eq:covAugNoiseLinear}
\end{align}

where we have $\underline{p}^{a}_{i,\frameVar|\frameVar-1}=\expv\left[ |\underline{n}^{a}_{i,\frameVar|\frameVar-1}|^2 \right]$ and $\underline{\mybold{R}}^{w}_{\frameVar}=\expv\left[\underline{\mybold{w}}_{\frameVar}\underline{\mybold{w}}^H_{\frameVar}\right]$. Due to the product term of the nonlinear estimation error and the channel vector, the overall term~\eqref{eq:augNoiseLinear} is non-Gaussian. As a consequence, the derivation of exact Kalman filter equations is not possible, since it requires jointly Gaussian distributions within the state-space model. Therefore, in the following, we approximate~\eqref{eq:augNoiseLinear} by an independent Gaussian random vector $\underline{\tilde{\mybold{s}}}^{w2}_{\frameVar}$, which has the same second-order moment as $\underline{\tilde{\mybold{s}}}^{w1}_{\frameVar}$.

Thus, the modified observation vector is
\begin{align}
\underline{\tilde{\mybold{y}}}_{\frameVar}
&\approx
\left(\sum_{i=0}^{N-1} \underline{\hat{a}}_{i,\frameVar|\frameVar-1}\underline{\mybold{C}}_{i,\frameVar}\right) \underline{\mybold{w}}_{\frameVar} + \underline{\tilde{\mybold{s}}}^{w2}_{\frameVar} \nonumber\\
&=
\underline{\mybold{C}}^{\hat{a}}_{\frameVar|\frameVar-1} \underline{\mybold{w}}_{\frameVar} + \underline{\tilde{\mybold{s}}}^{w2}_{\frameVar}
\label{eq:obsLinearGauss}
\end{align}
with
\begin{align}
\underline{\mybold{C}}^{\hat{a}}_{\frameVar|\frameVar-1} = \sum_{i=0}^{N-1} \underline{\hat{a}}_{i,\frameVar|\frameVar-1}\underline{\mybold{C}}_{i,\frameVar}.
\label{eq:bwTau}
\end{align}
We are now prepared to provide the Kalman gain and update equations~\cite{src:scharf1991statistical}, based on the observation equation~\eqref{eq:obsLinearGauss}, and by using~\eqref{eq:covAugNoiseLinear} and~\eqref{eq:bwTau}, and the definitions of error covariance matrices from Section~\ref{sec:predictions}. The computations for the $\frameVar^{\text{th}}$ step are depicted in~\algref{alg:linearChannel}.

\begin{algorithm}
\caption{Obtain linear SI channel estimate at step $\frameVar$}
\label{alg:linearChannel}
\begin{algorithmic}
\REQUIRE known $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}$, $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}$, $\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar-1}$ and $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$

\GAIN{}
\STATE

$\underline{\mybold{K}}^{w}_{\frameVar} =
\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{a}^H}_{\frameVar|\frameVar-1}\cdot$\newline
$\phantom{\underline{\mybold{K}}^{w}_{\frameVar} =
\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}
}
\left(
\underline{\mybold{C}}^{\hat{a}}_{\frameVar|\frameVar-1}
\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{a}^H}_{\frameVar|\frameVar-1}+
\underline{\mybold{\Psi}}^{\tilde{s}^{w2}}_{\frameVar}
\right)^{-1}
$
\ENDBLOCK
\UPDATE{linear SI channel estimation}
\STATE $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}=
\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}+
\underline{\mybold{K}}^{w}_{\frameVar}
\left(
\underline{\tilde{\mybold{y}}}_{\frameVar} - \underline{\mybold{C}}^{\hat{a}}_{\frameVar|\frameVar-1}\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}
\right)$
\ENDBLOCK
\UPDATE{linear SI channel estimation error covariance}
\STATE
$\underline{\mybold{P}}^{w}_{\frameVar|\frameVar} =
\left(
\eye{M}-\underline{\mybold{K}}^{w}_{\frameVar}\underline{\mybold{C}}^{\hat{a}}_{\frameVar|\frameVar-1}
\right)
\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}$
\ENDBLOCK
\end{algorithmic}
\end{algorithm}

\subsection{Update estimation of nonlinear coefficients}
\label{sec:estNl}
We first provide an alternative expression for the observation~\eqref{eq:obsDec} at the receiver. The adaptation with respect to the nonlinear coefficients~$\underline{\mybold{a}}_{\frameVar}$ requires another linear form of the observation. Thus, we express~\eqref{eq:obsDec} alternatively by using~\eqref{eq:nlCoeffsVec}
\begin{align}
\label{eq:obsAlt}

\underline{\tilde{\mybold{y}}}_{\frameVar}&= \underline{\mybold{C}}^{w}_{\frameVar} \underline{\mybold{a}}_{\frameVar} + \underline{\tilde{\mybold{s}}}_{\frameVar},
\end{align}
where we define

\begin{align}
\underline{\mybold{C}}^{w}_{\frameVar} =
\left[
\underline{\mybold{C}}_{0,\frameVar}\underline{\mybold{w}}_{\frameVar},
\ldots,
\underline{\mybold{C}}_{N-1,\frameVar}\underline{\mybold{w}}_{\frameVar}
\right].
\end{align}
Compared to the original form~\eqref{eq:obsDec}, we have essentially swapped the roles of~$\underline{\mybold{a}}_{\frameVar}$ and the linear SI channel~$\underline{\mybold{w}}_{\frameVar}$. Let
\begin{align}
\label{eq:y2dft}
\underline{\tilde{\mybold{y}}}_{\frameVar}=

\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar} \underline{\mybold{a}}_{\frameVar} + \underline{\tilde{\mybold{s}}}^{a1}_{\frameVar}
\end{align}
where we use~\eqref{eq:obsAlt} and the definition
\begin{align}
\label{eq:cwtau}
\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar} &=
\left[
\underline{\mybold{C}}_{0,\frameVar}\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar},
\underline{\mybold{C}}_{1,\frameVar}\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar},
\ldots,
\underline{\mybold{C}}_{N-1,\frameVar}\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}
\right] \\
&=
\begin{bmatrix}
\underline{\mybold{c}}^{\hat{w}}_{0,\frameVar|\frameVar} &
\underline{\mybold{c}}^{\hat{w}}_{1,\frameVar|\frameVar} &
\ldots &
\underline{\mybold{c}}^{\hat{w}}_{N-1,\frameVar|\frameVar}
\end{bmatrix}.
\label{eq:frameColsDef}
\end{align}

Furthermore, the zero-mean augmented noise here is
\begin{align}
\label{eq:SignalNoiseA}
\underline{\tilde{\mybold{s}}}^{a1}_{\frameVar} &=
\left[
\begin{matrix}
\underline{\mybold{C}}_{0,\frameVar}\underline{\mybold{n}}^{w}_{\frameVar|\frameVar}
& \underline{\mybold{C}}_{1,\frameVar}\underline{\mybold{n}}^{w}_{\frameVar|\frameVar}
& \ldots
\end{matrix} \right.\nonumber\\
&\qquad\left.
\begin{matrix}
& \underline{\mybold{C}}_{N-1,\frameVar}\underline{\mybold{n}}^{w}_{\frameVar|\frameVar}
\end{matrix}
\right]
\underline{\mybold{a}}_{\frameVar}

 + \underline{\tilde{\mybold{s}}}_{\frameVar}.
\end{align}
The covariance matrix of the augmented noise $\underline{\tilde{\mybold{s}}}^{a1}_{\frameVar}$ is
\begin{align}
\underline{\mybold{\Psi}}^{\tilde{s}^{a1}}_{\frameVar}
&=
\sum_{i=0}^{N-1}

\underline{p}^{a}_{i,\frameVar}
\underline{\mybold{C}}_{i,\frameVar}
\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}
\underline{\mybold{C}}^{H}_{i,\frameVar}
+\underline{\mybold{\Psi}}^{\tilde{s}}_{\frameVar}.

\label{eq:SignalNoiseACov}
\end{align}

with $\underline{p}^{a}_{i,\frameVar}=\expv\left[ |\underline{a}_{i,\frameVar}|^2 \right]$ and $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}=\expv\left[ \underline{\mybold{n}}^{w}_{\frameVar|\frameVar}\underline{\mybold{n}}^{w^{H}}_{\frameVar|\frameVar} \right]$.

As a starting point,  it is noteworthy to say that the estimation of the linear SI channel provided in Section~\ref{sec:estLin} depends on the predictions of $\underline{\mybold{w}}_{\frameVar}$ and $\underline{\mybold{a}}_{\frameVar}$ only, while here the updated $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}$ is available.

Similar to~\eqref{eq:augNoiseLinear}, we approximate~\eqref{eq:SignalNoiseA} by an independent Gaussian random vector $\underline{\tilde{\mybold{s}}}^{a2}_{\frameVar}$ with the same second-order moment as $\underline{\tilde{\mybold{s}}}^{a1}_{\frameVar}$. Thus, we have
\begin{align}
\label{eq:yaDef}
\underline{\tilde{\mybold{y}}}_{\frameVar}
&\approx
\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar} \underline{\mybold{a}}_{\frameVar}
 + \underline{\tilde{\mybold{s}}}^{a2}_{\frameVar}.
\end{align}
On this basis, we provide the Kalman filter equations, based on the system equation~\eqref{eq:nlCoeffsMarkov} and the observation equation~\eqref{eq:yaDef}, and by considering~\eqref{eq:SignalNoiseACov}. The computations for the $\frameVar^{\text{th}}$ step are depicted in~\algref{alg:nonlinearCoeffs}.
\begin{algorithm}
\caption{Obtain nonlinear coefficients estimate at step $\frameVar$}
\label{alg:nonlinearCoeffs}
\begin{algorithmic}
\REQUIRE known $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}$, $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}$, $\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar-1}$ and $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$

\GAIN{}
\STATE
$\underline{\mybold{K}}^{a}_{\frameVar} =
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}
\left(
\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar}
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}+
\underline{\mybold{\Psi}}^{\tilde{s}^{a2}}_{\frameVar}
\right)^{-1}$
\ENDBLOCK
\UPDATE{coefficients estimate}
\STATE $\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar}=
\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar-1}+
\underline{\mybold{K}}^{a}_{\frameVar}
\left(
\underline{\tilde{\mybold{y}}}_{\frameVar} - \underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar}\underline{\hat{\mybold{a}}}_{\frameVar|\frameVar-1}
\right)$
\ENDBLOCK
\UPDATE{coefficients estimation error covariance}
\STATE
$\underline{\mybold{P}}^{a}_{\frameVar|\frameVar} =
\left(
\eye{M}-\underline{\mybold{K}}^{a}_{\frameVar}\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar}
\right)
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$
\ENDBLOCK
\end{algorithmic}
\end{algorithm}

\section{Approximations}
\label{sec:approximations}
In this section, we discuss three approximations of the algorithm derived in Section~\ref{sec:algorithm} such as an intra-channel approximation (Section~\ref{sec:intraChannelApprox}), an overlap-save diagonalization (Section~\ref{sec:overlapSaveDiag}) and an nonlinear diagonalization (Section~\ref{sec:nonlinearDiag}).

\subsection{Intra-channel approximation}
\label{sec:intraChannelApprox}
The computation of noise covariances~\eqref{eq:covAugNoiseLinear} and~\eqref{eq:SignalNoiseACov} during the Kalman update step requires a-priori knowledge of the second moments $\underline{p}^{a}_{i,\frameVar}$, $\underline{\mybold{R}}^{w}_{\frameVar}$ of the linear channel and nonlinear coefficients, respectively. However, this information might not be available, and therefore we approximate the second moments by their estimated counterparts:
\begin{align}
\underline{\mybold{R}}^{w}_{\frameVar} &=
\expv\left[ \underline{\mybold{w}}_{\frameVar}\underline{\mybold{w}}^H_{\frameVar} \right] \nonumber \\
&\approx
\diag\left[
\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar-1}
\circ
\underline{\hat{\mybold{w}}}^{*}_{\frameVar|\frameVar-1}
\right]
+
\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}. \\
\underline{p}^{a}_{i,\frameVar} &=
\expv\left[
|\underline{a}_{i,\frameVar}|^2
\right] \nonumber \\
&\approx
|\underline{\hat{a}}_{i,\frameVar|\frameVar-1}|^2+
\underline{p}^{a}_{i,\frameVar|\frameVar-1},
\end{align}
where we approximate the linear channel and nonlinear coefficient covariance by considering the instantaneous representations of the variables~\eqref{eq:siChanEstimatorPred} and~\eqref{eq:nlCoeffsEstimatorPred} instead of the average.

Furthermore, we assume the error covariance matrices $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}$, $\underline{\mybold{P}}^{w}_{\frameVar|\frameVar}$ and the augmented noise covariance matrix of~\eqref{eq:covAugNoiseLinear} to be diagonal.

\subsection{Overlap-save diagonalization}
\label{sec:overlapSaveDiag}
The computation of the Kalman gain in \algref{alg:linearChannel} requires the inversion of an $M\times M$ matrix, which is both computationally intensive and numerically unstable. To simplify~\eqref{eq:CDef}, we apply the intra-channel diagonalization from Section~\ref{sec:intraChannelApprox} and use the approximations~\cite{src:benesty2001advances,Enzner06} for each nonlinear contribution
\begin{align}
\label{eq:approx1}
\underline{\mybold{C}}_{i,\frameVar}
&\approx
\frac{R}{M}\diag\left[ \underline{\mybold{\phi}}_{i,\frameVar}\right] \\
\label{eq:approx2}
\underline{\mybold{C}}_{i,\frameVar}\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}\underline{\mybold{C}}^{H}_{j,\frameVar}
&\approx
\frac{R}{M}\diag\left[ \underline{\mybold{\phi}}_{i,\frameVar}\right]\underline{\mybold{P}}^{w}_{\frameVar|\frameVar-1}\diag\left[ \underline{\mybold{\phi}}^{*}_{j,\frameVar}\right].
\end{align}
The approximation errors of~\eqref{eq:approx1},~\eqref{eq:approx2} vanish for $R\to M$~\cite{src:malik2013variational}. The diagonalization step essentially ignores the overlap-save constraint during fast convolution, however, this effect is small if the difference $M-R$ is rather low. Considering the aforementioned changes, the prediction, Kalman gain and update steps from~\algref{alg:linearChannel} are diagonalized.

\subsection{Nonlinear diagonalization}
\label{sec:nonlinearDiag}
The Kalman gain of \algref{alg:nonlinearCoeffs} uses the inversion of an $M\times M$ matrix. We intend to reduce computational complexity of that operation by diagonalizing that matrix. Since $\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}$ and $\underline{\mybold{\Psi}}^{s^{a2}}_{\frameVar}$ are already diagonal by the approximations of Sections~\ref{sec:intraChannelApprox} and~\ref{sec:overlapSaveDiag}, only the $M\times N$ frame matrix $\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar}$ has to be considered. Recall the definition~\eqref{eq:frameColsDef}. Each of the columns represents a nonlinear basis frame after convolution with the linear SI channel. Thus, we can write
\begin{align}
\label{eq:nldiagInv}
&\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar}
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}+
\underline{\mybold{\Psi}}^{\tilde{s}^{a2}}_{\frameVar} = \\
&\qquad=
\sum_{i=0}^{N-1}
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}
\underline{\mybold{c}}^{\hat{w}^H}_{i,\frameVar|\frameVar}+
\underline{\mybold{\Psi}}^{\tilde{s}^{a2}}_{\frameVar}, \notag \\
\label{eq:nlDiagInvApprox}
&\qquad
\hintedrel[eq:nldiagInv1]{\approx}
\sum_{i=0}^{N-1}
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}
\underline{\mybold{c}}^{\hat{w}^H}_{i,\frameVar|\frameVar}+
\underline{\sigma}^{2}_{\frameVar}\eye{M},
\\
&\qquad
\hintedrel[eq:nldiagInv2]{=}
\mybold{U} \diag\left[ \lambda_0, \lambda_1, \ldots, \lambda_{M-1} \right] \mybold{U}^H
\end{align}
where in~$(\hintref{eq:nldiagInv1})$, $\underline{\mybold{\Psi}}^{\tilde{s}^{a2}}_{\frameVar}$ is approximated by $\underline{\sigma}^{2}_{\frameVar}\eye{M}$ and $\underline{\sigma}^{2}_{\frameVar}$ denotes the largest element of $\underline{\mybold{\Psi}}^{\tilde{s}^{a2}}_{\frameVar}$. Furthermore, in~$(\hintref{eq:nldiagInv2})$, $\mybold{U}$ denotes the matrix of eigenvectors and $\lambda_i$ are the eigenvalues. The basis frames $\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}$ are not necessarily orthogonal. However, if we assume they are approximately orthogonal, then the basis frames are eigenvectors of
$
\underline{\mybold{C}}^{\hat{w}}_{\frameVar|\frameVar}
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}
$, and, since $\underline{\sigma}^{2}_{\frameVar}\eye{M}$ is just a scaled identity matrix, they are also eigenvectors of~\eqref{eq:nlDiagInvApprox} with eigenvalues
\begin{align*}
\lambda_i =
\begin{cases}
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\underline{\mybold{c}}^{\hat{w}^H}_{i,\frameVar|\frameVar}
\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}+
\underline{\sigma}^{2}_{\frameVar} &
0\leq i\leq N-1 \\
\underline{\sigma}^{2}_{\frameVar} &
N\leq i\leq M-1.
\end{cases}
\end{align*}
Furthermore, the first $N$ columns of $\mybold{U}$ are normalized versions of the $\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}$. Next, for the Kalman gain, we have
\begin{align*}
\underline{\mybold{K}}^{a}_{\frameVar} &=
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}
\left(
\sum_{i=0}^{N-1}
\underline{p}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}
\underline{\mybold{c}}^{\hat{w}^H}_{i,\frameVar|\frameVar}+
\underline{\sigma}^{2}_{\frameVar}\eye{M}
\right)^{-1} \\
&=
\underline{\mybold{P}}^{a}_{\frameVar|\frameVar-1}
\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}
\left(
\mybold{U} \diag\left[ \frac{1}{\lambda_0}, \frac{1}{\lambda_1}, \ldots, \frac{1}{\lambda_{M-1}} \right] \mybold{U}^H
\right).
\end{align*}
Let $\underline{\mybold{k}}^{a^{T}}_{i,\frameVar}$ be the $i$th row of  $\underline{\mybold{K}}^{a}_{\frameVar}$, then we have
\begin{align}
\label{eq:kalmanApprox}
\underline{\mybold{k}}^{a^{T}}_{i,\frameVar}=
\frac{
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
}{
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\underline{\mybold{c}}^{\hat{w}^H}_{i,\frameVar|\frameVar}
\underline{\mybold{c}}^{\hat{w}}_{i,\frameVar|\frameVar}+
\underline{\sigma}^{2}_{\frameVar}.
}
\underline{\mybold{c}}^{\hat{w}^H}_{i,\frameVar|\frameVar},
\end{align}
since the $i$th row of $\underline{\mybold{C}}^{\hat{w}^H}_{\frameVar|\frameVar}$ is orthogonal to all columns of $\mybold{U}$ except the $i$th one. The computation of~\eqref{eq:kalmanApprox} is much less expensive in terms of computational complexity. Finally, we apply the overlap-save diagonalization of~\eqref{eq:approx1} to~\eqref{eq:kalmanApprox}, and get
\begin{align*}

\underline{\mybold{k}}^{a^{T}}_{i,\frameVar}\approx
\frac{
\frac{R}{M}
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\left(
\underline{\mybold{\phi}}_{i,\frameVar}
\circ
\underline{\mybold{\hat{w}}}_{i,\frameVar|\frameVar}
\right)^H
}{
\frac{R}{M}
\underline{p}^{a}_{i,\frameVar|\frameVar-1}
\left(
\underline{\mybold{\phi}}_{i,\frameVar}
\circ
\underline{\mybold{\hat{w}}}_{i,\frameVar|\frameVar}
\right)^H
\left(
\underline{\mybold{\phi}}_{i,\frameVar}
\circ
\underline{\mybold{\hat{w}}}_{i,\frameVar|\frameVar}
\right)
+
\underline{\sigma}^{2}_{\frameVar}.
}.
\end{align*}

\section{Results}
\label{sec:results}

In this section, we evaluate the performance of the algorithms derived in this work by comparing them to other approaches in the literature. We refer to~\algref{alg:linearChannel},~\algref{alg:nonlinearCoeffs} from Section~\ref{sec:algorithm} as the \emph{exact} Kalman algorithm in cascade structure. In addition, the \emph{approximated} Kalman algorithm is taken from Section~\ref{sec:approximations}. From the literature, we employ the Kalman algorithm in parallel structure~\cite{MalikEnzner2012b} and both the time-domain nonlinear NLMS and RLS in parallel structure.

The section is organized as follows. First, we define the performance metrics in Section~\ref{sec:metrics}. Then, the computational complexity is analyzed in Section~\ref{sec:complexity}. By performing simulations, we look into time convergence behavior in Section~\ref{sec:convergence}, global performance of the proposed algorithms in Sections~\ref{sec:global} and use cases and~\ref{sec:usecases}.

\subsection{Performance metrics}
\label{sec:metrics}
We define the key metrics that serve as performance indicators for the adaptive SI cancellation algorithms. First, we introduce definitions according to the mean-squared error (MSE) principle. The average signal to residual-interference-and-noise ratio (SRINR) is given as follows
\begin{align}
\label{eq:srinrDef}
\text{SRINR}_{\frameVar}=
\frac{
\expv\left[\left| \mybold{d}^{h}_{\frameVar} \right|^2\right]
}{
\expv\left[\left| \mybold{e}_\frameVar-\mybold{d}^{h}_{\frameVar} \right|^2\right]
}
,
\end{align}
where $\mybold{d}^{h}_{\frameVar}$ denotes the desired signal at the receiver, and $\mybold{e}_\frameVar$ is the error signal before decoding~\eqref{eq:errorSignal} in time domain.

The wireless channel~\eqref{eq:wChanDef} is assumed to be random, thus the process~\eqref{eq:srinrDef} is non-ergodic, and therefore we define the ``global'' SRINR to be the average in time for $\frameVar\to\infty$ over many realizations of the wireless channel.

The capability of system identification can be measured by the system distance, i.e., the power of the estimation error compared to the power of the system variables. In the domain of Kalman filtering, the internal system state is the unknown quantity to be observed. To evaluate the quality of state estimation, we define the system distance for the linear SI channel coefficients as follows:
\begin{align}
\label{eq:sysDistWDef}
\text{SysDist}^{w}_{\frameVar} &=
\frac{
\expv\left[
\left( \mybold{w}_{\frameVar}-\hat{\mybold{w}}_{\frameVar|\frameVar} \right)^H  \left( \mybold{w}_{\frameVar}-\hat{\mybold{w}}_{\frameVar|\frameVar} \right)
\right]
}{
\expv\left[ \mybold{w}^H_{\frameVar} \mybold{w}_{\frameVar}\right]
}
,
\end{align}
where $\mybold{w}_{\frameVar}$ is the true SI channel vector per frame, as given by~\eqref{eq:siChanVec}, and $\hat{\mybold{w}}_{\frameVar|\frameVar}$ is the corresponding estimate of $\underline{\hat{\mybold{w}}}_{\frameVar|\frameVar}$ from~\algref{alg:linearChannel} in time domain. Similarly, the system distance for the $i$th nonlinear coefficient is defined as follows:
\begin{align}
\label{eq:sysDistADef}
\text{SysDist}^{a}_{i,\frameVar} &=
\frac{
\expv\left[
\left| \underline{a}_{i,\frameVar}-\underline{\hat{a}}_{i,\frameVar|\frameVar} \right|^2
\right]
}{
\expv\left[ \left| \underline{a}_{i,\frameVar} \right|^2 \right]
}
,
\end{align}
where $\underline{a}_{i,\frameVar}$ and $\underline{\hat{a}}_{i,\frameVar|\frameVar}$ are the nonlinear coefficient from~\eqref{eq:nlCoeffsVecEqual} and its updated estimate at frame $\frameVar$, respectively. We need a fair comparison of the algorithms based on either cascaded (see~\figref{fig:systemmodel}) or parallel (see~\figref{fig:parallelmodel}) channel structure. The coefficients $\hat{\mybold{w}}_{0,\frameVar|\frameVar}$ of the parallel structure correspond to the linear SI channel from the cascade model. The nonlinear coefficients $\underline{\hat{a}}_{i,\frameVar}$ in case of a parallel structure are hence obtained by a Least-Squares solution
\begin{align}
\label{eq:nlLSEst}
\underline{\hat{a}}_{i,\frameVar|\frameVar}=\frac{\hat{\mybold{w}}^H_{0,\frameVar|\frameVar}\hat{\mybold{w}}_{i,\frameVar|\frameVar}}{\hat{\mybold{w}}^H_{0,\frameVar|\frameVar}\hat{\mybold{w}}_{0,\frameVar|\frameVar}}.
\end{align}

Next, we turn to a metric of digital communications. To simplify the analysis, we assume the wireless channel to be constant over time. Furthermore, since after SI cancellation the residual SI can be approximated by white noise~\cite{src:day2012fullduplex},  we assume that the residual SI process is independent and ergodic. Thus, the information capacity with power constraint $P_d$ at the distant node is given by~\cite{src:cover2006elements}
\begin{align}
C &=\max_{f(d):\expv\left[|d|^2\right]<P_d}\:I\left(d;y\right) \nonumber\\
&\hintedrel[eq:ratelower1]{\geq}\max_{f(d):\expv\left[|d|^2\right]<P_d}\:I\left(d;y_G\right) \nonumber\\
\label{eq:ratelower}
&\hintedrel[eq:ratelower2]{=} \log\left( 1+\frac{P_d\left\lVert \mybold{h} \right\rVert^2_2}{\sigma^2_{\tilde{e}}} \right) =: R
\end{align}
where in~$(\hintref{eq:ratelower1})$ we approximate the received signal $y$ by a circular-symmetric complex Gaussian counterpart $y_G$, since Gaussian noise serves as the worst-case noise~\cite{src:shomorony2012Is}. We arrive at~$(\hintref{eq:ratelower2})$, since in a Gaussian point-to-point channel, it is known that a Gaussian input is optimal. Here, $\sigma^2_{\tilde{e}}$ denotes the variance of the residual SI plus independent noise.

\subsection{Complexity}
\label{sec:complexity}
The cascade model of the nonlinear SI channel (recall~\figref{fig:systemmodel}) is essentially motivated by reducing computational complexity of the adaptive cancellation algorithm, which is primarily determined by the number of multiplications and divisions. In our analysis, we focus on the principal complexity only and do not evaluate the cost of each calculation step in detail. Thus, we identify those parts of the algorithms that exhibit the most significant impact on the complexity, depending on the (potentially large) frame shift $R$, the FIR filter length $\firLengthSI$ and on the order of the nonlinear expansion $N$. \tabref{tab:complexity} shows the results. In order to compare time-domain and DFT domain SI estimation and cancellation algorithms, we refer to the complexity as per sample. In all cases, DFT domain approaches are normalized to the frame shift $R=M-\firLengthSI$. This is reducing complexity by an order of magnitude if $R$ is growing while $\firLengthSI$ is kept fixed.

Now, consider the frame shift $R$. The exact Kalman algorithm in cascade structure, derived in Section~\ref{sec:algorithm}, does have at least quadratic complexity, since it requires the inversion and multiplication of $M\times M$ matrices. On the other hand, considering the approximations of Section~\ref{sec:approximations}, the DFT operation becomes dominant, and therefore we have logarithmic complexity for the Kalman algorithm with nonlinear diagonalization. Similar reasoning holds for the Kalman algorithm in parallel structure after submatrix or full diagonalization~\cite{MalikEnzner2012b}. The complexity of the time-domain algorithms is determined by the the FIR filter size $\firLengthSI$. The classic NLMS algorithm in time domain exhibits linear complexity in $\firLengthSI$, while the parallel RLS algorithm in time domain has quadratic complexity due to the update step. However, a direct comparison of time-domain and DFT domain complexities is more involved and depends on the priorities of Kalman filter design. Certain designs might increase $\firLengthSI$ for finer channel estimation (and thus keep $R$ essentially constant), while others might use larger frame shifts $R$ to diminish the impact of the overlap-save approximation.

The number of nonlinear basis functions $N$ affects the computational complexity of the various algorithms in different ways. In cascade structure, the complexity generally scales linearly with respect to $N$. On the other hand, in parallel structures, Kalman or RLS algorithms take the correlation between all sub-channels into account, and therefore require cubic and quadratic complexity in $N$, respectively. Especially the cubic term might be greater than $\log_2M$ even for small $N$ and thus can have a significant impact. Kalman algorithm with full diagonalization~\cite{MalikEnzner2012b} or the NLMS allow to reduce the complexity to linear scale, since the correlation between parallel channels is neglected.
\begin{table}
\begin{center}
  \begin{tabular}{ | l | l |}
  \hline

    Kalman, Cascade, Exact~\ref{sec:algorithm} & $\mathcal{O}\left(N\frac{(R+L)^3}{R}\right)$ \\ \hline
    Kalman, Cascade, Approx.~\ref{sec:nonlinearDiag} & $\mathcal{O}\left(N\frac{(R+L)\log_2(R+L)}{R}\right)$  \\ \hline
    Kalman, Parallel, Sub. diag.~\cite{MalikEnzner2012b}& $\mathcal{O}\left(N^3\frac{R+L}{R}\right) + \mathcal{O}\left(N\frac{(R+L)\log_2(R+L)}{R}\right)$ \\ \hline
    Kalman, Parallel, Full diag.~\cite{MalikEnzner2012b} & $\mathcal{O}\left(N\frac{(R+L)\log_2(R+L)}{R}\right)$ \\ \hline
    NLMS, Parallel, Time domain & $\mathcal{O}(NL)$  \\ \hline
    RLS, Parallel, Time domain & $\mathcal{O}(N^2L^2)$  \\ \hline
  \end{tabular}
\end{center}
\caption{Comparison of computational complexity with respect to frame shift $R$, FIR filter length $\firLengthSI$ and expansion order $N$ for different SI estimation and cancellation algorithms.}
\label{tab:complexity}
\end{table}

\subsection{Time convergence behavior}
\label{sec:convergence}
In this Section, we analyze the performance of the adaptive estimations over the course of time. We evaluate the performance of the algorithm by simulations.  We choose the parameters $M=64$, $R=56$ and select $N=3$ basis functions. The basis function are chosen according to the following emphases: (i) the linear part $\phi_0(x_{\timeVar})=x_{\timeVar}$ of the SI, then (ii) the widely-linear component $\phi_1(x_{\timeVar})=x^*_{\timeVar}$~\cite{src:korpi2014widely} and finally, (iii), the $3^{\text{rd}}$ order nonlinear component $\phi_2(x_{\timeVar})=x^2_{\timeVar}x^*_{\timeVar}$ with initial nonlinear coefficients $\underline{a}_{1,0}=\underline{a}_{2,0}=-10$~dB. The additive noise level is set to $-50$~dB. For the SI signal and the SoI, we generate zero-mean, standard normal i.i.d. samples. The Kalman filter parameters are perfectly matched to the simulated time-variant SI channel conditions. The NLMS step size is $10^{-2}$, while the RLS forgetting factor is matched to the Markov model of~\eqref{eq:siChanMarkov}.

 Let the \emph{input SINR} be the ratio of the SoI power to the SI-and-noise power at the receiver before any signal processing occurs. The input SINR is fixed to $-15$~dB, thus, the SI is much stronger than the SoI. Consider~\figref{fig:convergence}, which shows the time convergence behavior of SI cancellation algorithms over the frame index $\frameVar$.
\begin{figure*}
\ifshowTikzPictures
\begin{tikzpicture}
   \begin{groupplot}[group style={
                      group name=framePlots,
                      horizontal sep=3cm,
                      vertical sep=2cm,
                      group size= 2 by 2},
                      height=5cm,
                      width=8cm,
                      xmin=0,
                      xmax=400,
                      grid=major]
        \pgfmathtruncatemacro{\markRepeatFaster}{\markRepeat/2}
        \nextgroupplot[ylabel={SRINR [dB]}, xlabel={$\frameVar$}]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=\markRepeatFaster,mark phase=\markPOne]
			  table[]{\resSrinrVsSinr NonOrtho-Static-1.tsv};\label{plots:plot1a}

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=\markRepeat,mark phase=\markPFour]
			  table[]{\resSrinrVsSinr NonOrtho-Static-4.tsv};\label{plots:plot4d}
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=\markRepeat,mark phase=\markPFive]
			  table[]{\resSrinrVsSinr NonOrtho-Static-5.tsv};\label{plots:plot5e}

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=\markRepeat,mark phase=\markPSeven]
			  table[]{\resSrinrVsSinr NonOrtho-Static-7.tsv};\label{plots:plot7g}
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=\markRepeat,mark phase=\markPEight]
			  table[]{\resSrinrVsSinr NonOrtho-Static-8.tsv};\label{plots:plot8h}

        \nextgroupplot[ylabel={System distance [dB]}, xlabel={$\frameVar$}, ymin=-45, ymax=20]

			\addplot [color=blue, mark=asterisk, mark options={solid, blue}, mark repeat=\markRepeat,mark phase=\markPThree]
			  table[]{\resSysDistWVsSinr NonOrtho-Static-3.tsv};
			\addplot [color=cyan, mark=star, dashdotted, mark repeat=\markRepeat,mark phase=\markPFour]
			  table[]{\resSysDistWVsSinr NonOrtho-Static-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=\markRepeat,mark phase=\markPFive]
			  table[]{\resSysDistWVsSinr NonOrtho-Static-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=\markRepeat,mark phase=\markPSeven]
			  table[]{\resSysDistWVsSinr NonOrtho-Static-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=\markRepeat,mark phase=\markPEight]
			  table[]{\resSysDistWVsSinr NonOrtho-Static-8.tsv};

        \nextgroupplot[ylabel={System distance [dB]}, xlabel={$\frameVar$}, ymin=-45, ymax=20]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=\markRepeatFaster,mark phase=\markPOne]
			  table[]{\resSysDistAOneVsSinr NonOrtho-Static-1.tsv};

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=\markRepeat,mark phase=\markPFour]
			  table[]{\resSysDistAOneVsSinr NonOrtho-Static-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=\markRepeat,mark phase=\markPFive]
			  table[]{\resSysDistAOneVsSinr NonOrtho-Static-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=\markRepeat,mark phase=\markPSeven]
			  table[]{\resSysDistAOneVsSinr NonOrtho-Static-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=\markRepeat,mark phase=\markPEight]
				  table[]{\resSysDistAOneVsSinr NonOrtho-Static-8.tsv};

        \nextgroupplot[xlabel={Rounds}, ylabel={System distance [dB]}, xlabel={$\frameVar$}, ymin=-45, ymax=20]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=\markRepeatFaster,mark phase=\markPOne]
			  table[]{\resSysDistATwoVsSinr NonOrtho-Static-1.tsv};

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=\markRepeat,mark phase=\markPFour]
			  table[]{\resSysDistATwoVsSinr NonOrtho-Static-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=\markRepeat,mark phase=\markPFive]
			  table[]{\resSysDistATwoVsSinr NonOrtho-Static-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=\markRepeat,mark phase=\markPSeven]
			  table[]{\resSysDistATwoVsSinr NonOrtho-Static-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=\markRepeat,mark phase=\markPEight]
			  table[]{\resSysDistATwoVsSinr NonOrtho-Static-8.tsv};
    \end{groupplot}

\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]framePlots c1r1.south) {\subcaption{Signal-to-residual-interference-and-noise ratios~\eqref{eq:srinrDef}.\label{fig:srinrNonOrtho}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]framePlots c2r1.south) {\subcaption{System distance~\eqref{eq:sysDistWDef} to the linear SI channel $\underline{\mybold{w}}_{\frameVar}$\label{fig:sysdistWNonOrtho}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]framePlots c1r2.south) {\subcaption{System distance~\eqref{eq:sysDistADef} to the nonlinear coefficient $\underline{a}_{1,\frameVar}$\label{fig:sysdistA1NonOrtho}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]framePlots c2r2.south) {\subcaption{System distance~\eqref{eq:sysDistADef} to the nonlinear coefficient $\underline{a}_{2,\frameVar}$\label{fig:sysdistA2NonOrtho}}};

\path (framePlots c1r1.north west|-current bounding box.north)--
      coordinate(legendposFramePlots)
      (framePlots c2r1.north east|-current bounding box.north);
\matrix[
    matrix of nodes,
    anchor=south,
    draw,
    inner sep=0.2em,
    nodes={anchor=west}
  ]at([yshift=1ex]legendposFramePlots)
  {
    \ref*{plots:plot1a}& Kalman, Cascade, Exact~\ref{sec:algorithm} &[5pt]

    \ref*{plots:plot4d}& Kalman, Cascade, Approximated~\ref{sec:nonlinearDiag} & [5pt]
    \ref*{plots:plot5e}& Kalman, Parallel, Sub. diag.~\cite{MalikEnzner2012b}& \\

    \ref*{plots:plot7g}& NLMS, Parallel, time domain & [5pt]
    \ref*{plots:plot8h}& RLS, Parallel, time domain  \\
    };

\end{tikzpicture}
\fi
\vspace{-0.25cm}
\caption{Time convergence performance of SI cancellation algorithms over the frame index $\frameVar$ with non-orthogonalized inputs with static SI path}
\label{fig:convergence}
\end{figure*}
The inputs are not orthogonalized. Both the linear SI channel and the nonlinear coefficients remain constant over time. The SRINR results of~\figref{fig:srinrNonOrtho} demonstrate that without prior orthogonalization, the algorithms with strong inherent decorrelation feature, such as the Kalman in parallel structure with submatrix diagonalization or the RLS, are advantageous. Conversely, the NLMS is much slower in convergence. The Kalman approaches with cascade structure are bounded from above by the RLS and from below by the NLMS. Interestingly, the performance difference between the exact Kalman approach in cascade structure and its approximated counterpart

are almost negligible. Sometimes the exact calculation performs slightly worse. This has been identified as a particularity of the iterative design~\algref{alg:linearChannel} and~\algref{alg:nonlinearCoeffs}, which does not jointly estimate the linear SI channel and the nonlinear coefficients. The results of the system distance~(\Cref{fig:sysdistWNonOrtho,fig:sysdistA1NonOrtho,fig:sysdistA2NonOrtho}) essentially support the aforementioned findings. It can be observed that the system distances of the linear SI channel coefficients $\underline{\mybold{w}}_{\frameVar}$ of~\figref{fig:sysdistWNonOrtho} and the first nonlinear coefficient $\underline{a}_{1,\frameVar}$ in~\figref{fig:sysdistA1NonOrtho} are ranked consistently. In contrast to that, the distance of second nonlinear coefficient $\underline{a}_{2,\frameVar}$ in~\figref{fig:sysdistA2NonOrtho} is much larger, especially for the NLMS. This is due to the fact that the input signal and its $3$rd order component are correlated, and thus the performance deteriorates if the algorithms lack a decorrelation feature. The system distances in~\Cref{fig:sysdistWNonOrtho,fig:sysdistA1NonOrtho,fig:sysdistA2NonOrtho} of the Kalman approach in parallel structure with submatrix diagonalization and the RLS are almost identical. But the SRINR result (\figref{fig:srinrNonOrtho}) indicates that the Kalman algorithm does not precisely achieve the same level as the RLS. However, this observation can be explained. Both algorithms are devised in parallel structure. Thus, the nonlinear coefficients are computed by~\eqref{eq:nlLSEst}, which essentially is averaging over the channel coefficients and thus removing potential estimation errors. The parallel Kalman algorithm is slightly more susceptible to errors than the RLS since it is formulated in DFT domain with overlap-save approximation similar to the method described in sub-section~\ref{sec:overlapSaveDiag}. These error phenomena are compensated by reducing the number of estimated variables, such as by using~\eqref{eq:nlLSEst}.

\subsection{Global performance}
\label{sec:global}
For the global performance, we evaluate SRINR (after cancellation) and system distances with respect to certain input SINR (before cancellation), which is illustrated in~\figref{fig:global}.
\begin{figure*}
\ifshowTikzPictures
\begin{tikzpicture}
   \begin{groupplot}[group style={
                      group name=globalPlots,
                      horizontal sep=3cm,
                      vertical sep=2.5cm,
                      group size= 2 by 1},
                      height=5cm,
                      width=8cm,
                      xmin=-20,
                      xmax=20,
                      grid=major]

        \nextgroupplot[ylabel={SRINR [dB]}, xlabel={Input SINR [dB]}, ymin=-5, ymax=25]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=2,mark phase=1]
			  table[]{\resSrinrConvVsSinr NonOrtho-Static-1.tsv};\label{plots:global1}

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=2]
			  table[]{\resSrinrConvVsSinr NonOrtho-Static-4.tsv};\label{plots:global4}
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=2,mark phase=1]
			  table[]{\resSrinrConvVsSinr NonOrtho-Static-5.tsv};\label{plots:global5}

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=2,mark phase=1]
			  table[]{\resSrinrConvVsSinr NonOrtho-Static-7.tsv};\label{plots:global7}
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=2]
			  table[]{\resSrinrConvVsSinr NonOrtho-Static-8.tsv};\label{plots:global8}
			\addplot [color=black, mark options={none}, thick] coordinates {(0,0) (20,20)};
			\node[text width=2cm,align=center] at ($(axis cs:14,2.5)+(0,0)$) {Treating-interference-as-noise threshold};
			\addplot [color=black, mark options={none}, thick] coordinates {(-20,20) (20,20)};
			\node[text width=2cm,align=center] at ($(axis cs:-10,23)+(0,0)$) {SNR};
			\addplot [color=black, mark options={none}, thick] coordinates {(-20,0) (0,0)};
			\node[text width=2cm,align=center] at ($(axis cs:-10,3)+(0,0)$) {min. SRINR};

        \nextgroupplot[ylabel={System distance [dB]}, xlabel={Input SINR [dB]}, ymin=-50,ymax=20]
			\addplot [color=red, mark=+, mark options={solid, red},mark repeat=2, mark phase=1]
			  table[]{\resSysDistWConvVsSinr NonOrtho-Static-1.tsv};

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2, mark phase=2]
			  table[]{\resSysDistWConvVsSinr NonOrtho-Static-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=2,mark phase=1]
			  table[]{\resSysDistWConvVsSinr NonOrtho-Static-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=2,mark phase=1]
			  table[]{\resSysDistWConvVsSinr NonOrtho-Static-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=2]
			  table[]{\resSysDistWConvVsSinr NonOrtho-Static-8.tsv};

    \end{groupplot}

\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]globalPlots c1r1.south) {\subcaption{Signal-to-residual -interference-and-noise ratios~\eqref{eq:srinrDef}\label{fig:srinrConvNonOrtho}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]globalPlots c2r1.south) {\subcaption{System distance~\eqref{eq:sysDistWDef} to the linear SI channel $\underline{\mybold{w}}_{\frameVar}$\label{fig:sysdistWConvNonOrtho}}};

\path (globalPlots c1r1.north west|-current bounding box.north)--
      coordinate(legendposglobalPlots)
      (globalPlots c2r1.north east|-current bounding box.north);
\matrix[
    matrix of nodes,
    anchor=south,
    draw,
    inner sep=0.2em,
    nodes={anchor=west}
  ]at([yshift=1ex]legendposglobalPlots)
  {
    \ref*{plots:global1}& Kalman, Cascade, Exact~\ref{sec:algorithm} &[5pt]

    \ref*{plots:global4}& Kalman, Cascade, Approximated~\ref{sec:nonlinearDiag} & [5pt]
    \ref*{plots:global5}& Kalman, Parallel, Sub. diag.~\cite{MalikEnzner2012b}&\\

    \ref*{plots:global7}& NLMS, Parallel, time domain &
    \ref*{plots:global8}& RLS, Parallel, time domain  \\
    };

\end{tikzpicture}
\fi
\vspace{-0.25cm}
\caption{Global performance of SI cancellation algorithms over input SINR with non-orthogonalized inputs and static SI path.}
\label{fig:global}
\end{figure*}
Again, the input signals are not orthogonalized and both the linear SI channel and the nonlinear coefficients remain static. Consider~\figref{fig:srinrConvNonOrtho} first, which depicts the SRINR. The solid black line labeled ``SNR'' designates the maximum possible SRINR, since the noise floor is fixed at $-20$~dB with respect to the SoI. If the SRINR vs. input SINR graph crosses the ``treating-interference-as-noise threshold'', then it is no longer reasonable to use the algorithm and better to rather accept the SI as noise. The ``min. SRINR'' line is simply attained when the algorithm selects $\underline{\hat{\mybold{x}}}_{\text{si},\frameVar}=\underline{\mybold{y}}_{\frameVar}$ as output such that the received signal is completely suppressed. The results show that all algorithms can be used over the whole range of SINR scenarios. However, as already previously indicated, algorithms with already existing decorrelation property (Kalman in parallel structure with submatrix diagonalization or the RLS) have an advantage. This is especially prominent for the low SINR regime, where a great amount of cancellation is required. On the other hand, in the high SINR regime, the performance difference between the algorithms significantly diminishes. The system distance results of the linear SI channel coefficients $\underline{\mybold{w}}_{\frameVar}$ are shown in~\Cref{fig:sysdistWConvNonOrtho}. If we compare the performance of the various algorithms, the graphs are essentially reflecting the expected ``upside-down'' pattern compared to the SRINR figure, since a negative system distance naturally pays off in terms of SRINR. Here, however, the overall system distance is much better in the low SINR regime. This is intuitive, since at low SINR, the SI power is much larger than the SoI power, and thus the SI system components can be identified more accurately. Conversely, in the high SINR regime, the SI is buried in the SoI and the noise floor, and therefore the coefficient estimations are less reliable. While the detailed system identification is more difficult in the high SINR regime, the SRINR performance (denoting the overall SI estimation) is, however, not degraded.

 \subsection{Use cases}
\label{sec:usecases}
In this Section, we investigate the impact of time-variant channels on the performance of the SI estimation and cancellation algorithms. Static and time-variant channels represent fundamentally different use cases that significantly influence the proper choice of algorithms. As a performance metric, we apply the communication rate~\eqref{eq:ratelower}. The time-variant case is characterized by the (normalized) coherence time $\frameVar^{w}_{\text{coh}}=10^3$ for the linear SI channel and $\frameVar^{a}_{\text{coh}}=10^4$ for the nonlinear coefficients. We assume that the algorithm parameters are perfectly matched to the simulated channel conditions. Now consider~\figref{fig:rates}, where we illustrate the rate over the input SINR for different cases of temporal variations with or without orthogonalization.
 \begin{figure*}
\ifshowTikzPictures
\begin{tikzpicture}[spy using outlines=
                      	{rectangle, magnification=4, connect spies}

]
   \begin{groupplot}[group style={
                      group name=ratePlots,
                      horizontal sep=3cm,
                      vertical sep=2cm,
                      group size= 2 by 2},
                      height=5cm,
                      width=8cm,
                      xmin=-20,
                      xmax=20,
                      grid=major
                     ]

        \nextgroupplot[ylabel={Rate [bits/sample]}, xlabel={Input SINR [dB]}, ymin=0]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Static-1.tsv};\label{plots:rate1}

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr NonOrtho-Static-4.tsv};\label{plots:rate4}
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Static-5.tsv};\label{plots:rate5}

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Static-7.tsv};\label{plots:rate7}
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr NonOrtho-Static-8.tsv};\label{plots:rate8}
			\addplot [color=blue, mark=none, thick, mark options={solid, blue}]
			  table[]{\resRateVsSinr NonOrtho-Static-9.tsv};\label{plots:rate9}

        \nextgroupplot[ylabel={Rate [bits/sample]}, xlabel={Input SINR [dB]}, ymin=0]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr Orthogonal-Static-1.tsv};

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr Orthogonal-Static-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr Orthogonal-Static-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr Orthogonal-Static-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr Orthogonal-Static-8.tsv};
			\addplot [color=blue, mark=none, thick, mark options={solid, blue}]
			  table[]{\resRateVsSinr Orthogonal-Static-9.tsv};

        \nextgroupplot[ylabel={Rate [bits/sample]}, xlabel={Input SINR [dB]}, ymin=0]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Varying-1.tsv};

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr NonOrtho-Varying-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Varying-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Varying-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr NonOrtho-Varying-8.tsv};
			\addplot [color=blue, mark=none, thick, mark options={solid, blue}]
			  table[]{\resRateVsSinr NonOrtho-Varying-9.tsv};
			\coordinate (spypoint1) at (axis cs:-15,2);
			\coordinate (magnifyglass1) at (axis cs:12,2.25);

        \nextgroupplot[ylabel={Rate [bits/sample]}, xlabel={Input SINR [dB]}, ymin=0]
			\addplot [color=red, mark=+, mark options={solid, red}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr Orthogonal-Varying-1.tsv};

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr Orthogonal-Varying-4.tsv};
			\addplot [color=magenta, mark=x, mark options={solid, magenta}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr Orthogonal-Varying-5.tsv};

			\addplot [color=black, mark=diamond, mark options={solid, black}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr Orthogonal-Varying-7.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr Orthogonal-Varying-8.tsv};
			\addplot [color=blue, mark=none, thick, mark options={solid, blue}]
			  table[]{\resRateVsSinr Orthogonal-Varying-9.tsv};
			\coordinate (spypoint2) at (axis cs:-15,2);
			\coordinate (magnifyglass2) at (axis cs:12,2.25);

    \end{groupplot}

\spy [black, width=2.25cm, height=2.25cm] on (spypoint2)
	in node at (magnifyglass2);

\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]ratePlots c1r1.south) {\subcaption{Non-orthogonalized inputs, static SI path.\label{fig:rateNonOrthoStatic}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]ratePlots c2r1.south) {\subcaption{Orthogonalized inputs, static SI path.\label{fig:rateOrthoStatic}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]ratePlots c1r2.south) {\subcaption{Non-orthogonalized inputs, time-variant SI path.\label{fig:rateNonOrthoVarying}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]ratePlots c2r2.south) {\subcaption{Orthogonalized inputs, time-variant SI path.\label{fig:rateOrthoVarying}}};

\path (ratePlots c1r1.north west|-current bounding box.north)--
      coordinate(legendposratePlots)
      (ratePlots c2r1.north east|-current bounding box.north);
\matrix[
    matrix of nodes,
    anchor=south,
    draw,
    inner sep=0.2em,
    nodes={anchor=west}
  ]at([yshift=1ex]legendposratePlots)
  {
    \ref*{plots:rate1}& Kalman, Cascade, Exact~\ref{sec:algorithm} &[5pt]

    \ref*{plots:rate4}& Kalman, Cascade, Approximated~\ref{sec:nonlinearDiag} & [5pt]
    \ref*{plots:rate5}& Kalman, Parallel, Sub. diag.~\cite{MalikEnzner2012b}& \\

    \ref*{plots:rate7}& NLMS, Parallel, time domain & [5pt]
    \ref*{plots:rate8}& RLS, Parallel, time domain & [5pt]
    \ref*{plots:rate9}& Capacity &[5pt] \\

    };
\end{tikzpicture}
\fi
\vspace{-0.25cm}
\caption{Communication rate~\eqref{eq:ratelower} over the input SINR for different cases of input processing and temporal variations.}
\label{fig:rates}
\end{figure*}
In each figure, we have provided the capacity as the upper bound for all rates. In~\figref{fig:rateNonOrthoStatic}, we have depicted the rate for static conditions of both the linear SI channel and the nonlinear coefficients, while the input signals are not orthogonalized. This is essentially a similar situation shown as in~\figref{fig:srinrConvNonOrtho}, but now the performance gaps in the low SINR regime are even more pronounced. Algorithms with inherent decorrelation like the Kalman in parallel structure with submatrix diagonalization perform the best regardless of the input SINR, while the RLS actually almost achieves the capacity. The cascaded and computationally efficient Kalman algorithms proposed in this work have a middle position, while the simple NMLS exhibits worse performance.

However, if orthogonalization of the input signal is applied, the situation fundamentally changes, as shown in~\figref{fig:rateOrthoStatic}. Now all the algorithm performances are completely independent of the input SINR. Furthermore, the performance gaps are almost insignificant, and thus a system designer can switch to a different metric in order to choose the best algorithm, such as the computational complexity.

In the time-convergence behavior, where we can obverse the same evolution for all the algorithms except the NLMS, but we omit the figure due to space constraints. In addition, the system distance results (not shown here) are very similar for all studied algorithms.

\figref{fig:rateNonOrthoVarying} illustrates the communication rate in time-variant scenario, but without using input orthogonalization. Compared to the result of~\figref{fig:srinrConvNonOrtho}, it can clearly be seen that now the input SINR has a significant impact on the performance of all algorithms.

\figref{fig:rateOrthoVarying} shows the results under time-variant conditions for orthogonalized input signals. In this case, the orthogonalization does not seem to help either, since the temporal variations represent the overall performance limitation.

\subsection{Decoding}
\label{sec:decodingRes}
The decoding step has been introduced in Section~\ref{sec:decoding}. It takes place after the Kalman prediction, right before the update step is performed. Now, we evaluate the benefits of the information decoding to the SI estimation and cancellation. We essentially regard two relevant cases. First, we consider the situation that the algorithm does not decode the desired information at all. As a consequence, the SoI is considered as unknown noise to the adaptive algorithm, and thus the estimation performance is equal to the results of the previous Sections. Second, in the case of perfect decoding, the desired information is obtained without any errors, and thus the remaining residual error signal contains the residual SI only and some independent noise with $\text{SNR}=20$~dB.

Consider~\figref{fig:DecRateNonOrthoStatic}.
\begin{figure*}
\ifshowTikzPictures
\begin{tikzpicture}
   \begin{groupplot}[group style={
                      group name=decPlots,
                      horizontal sep=3cm,
                      vertical sep=2.5cm,
                      group size= 2 by 1},
                      height=5cm,
                      width=8cm,
                      xmin=-20,
                      xmax=20,
                      grid=major]

        \nextgroupplot[ylabel={Rate [bits/sample]}, xlabel={Input SINR [dB]}, ymin=0]

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Static-4.tsv};\label{plots:dec4}

			\addplot [color=orange, mark=square, mark options={solid, orange}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr PerfectDecNonOrtho-Static-4.tsv};\label{plots:dec6}
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr PerfectDecNonOrtho-Static-8.tsv};\label{plots:dec8}
			\addplot [color=green, mark=diamond, mark options={solid, green}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr PerfectDecNonOrtho-Static-10.tsv};\label{plots:dec10}
  			\addplot [color=blue, mark=none, thick, mark options={solid, blue}]
			 table[]{\resRateVsSinr Orthogonal-Static-11.tsv};\label{plots:dec11}

        \nextgroupplot[ylabel={Rate [bits/sample]}, xlabel={Input SINR [dB]}, ymin=0]

			\addplot [color=cyan, mark=star, dashdotted, mark repeat=2,mark phase=1]
			  table[]{\resRateVsSinr NonOrtho-Varying-4.tsv};

			\addplot [color=orange, mark=square, mark options={solid, orange}, mark repeat=2,mark phase=2]
			  table[]{\resRateVsSinr PerfectDecNonOrtho-Varying-4.tsv};
			\addplot [color=blue, mark=triangle, mark options={solid, blue}, mark repeat=2,mark phase=1]
						  table[]{\resRateVsSinr PerfectDecNonOrtho-Varying-8.tsv};
			\addplot [color=green, mark=diamond, mark options={solid, green}, mark repeat=2,mark phase=2]
						  table[]{\resRateVsSinr PerfectDecNonOrtho-Varying-10.tsv};
  			\addplot [color=blue, mark=none, thick, mark options={solid, blue}]
				 table[]{\resRateVsSinr Orthogonal-Varying-11.tsv};

    \end{groupplot}

\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]decPlots c1r1.south) {\subcaption{Communication rate~\eqref{eq:ratelower} for static SI path.\label{fig:DecRateNonOrthoStatic}}};
\node[text width=8cm,align=center,anchor=north] at ([yshift=-\subcaptionSpacing cm]decPlots c2r1.south) {\subcaption{Communication rate~\eqref{eq:ratelower} for time-variant SI path.\label{fig:DecRateNonOrthoVarying}}};

\path (decPlots c1r1.north west|-current bounding box.north)--
      coordinate(legendposdecPlots)
      (decPlots c2r1.north east|-current bounding box.north);
\matrix[
    matrix of nodes,
    anchor=south,
    draw,
    inner sep=0.2em,
    nodes={anchor=west}
  ]at([yshift=1ex]legendposdecPlots)
  {

    \ref*{plots:dec4}& Kalman, Cascade, Approximated~\ref{sec:nonlinearDiag}, No decoding &

    \ref*{plots:dec6}& Kalman, Cascade, Approximated~\ref{sec:nonlinearDiag}, Perfect decoding & \\
    \ref*{plots:dec8}& RLS, Parallel, time domain, No decoding & [5pt]
    \ref*{plots:dec10}& RLS, Parallel, time domain, Perfect decoding & [5pt] \\
    \ref*{plots:dec11}& Capacity &[5pt] \\
    };
\end{tikzpicture}
\fi
\vspace{-0.25cm}
\caption{Impact of decoding strategies on the SI estimation and cancellation performance.}
\label{fig:decoding}
\end{figure*}
It shows the communication rate~\eqref{eq:ratelower} for different input SINR in static~(\figref{fig:DecRateNonOrthoStatic}) and time-variant~(\figref{fig:DecRateNonOrthoVarying}) environments without input orthogonalization.

We keep our focus on Kalman algorithms in cascade structure with approximations~(Section~\ref{sec:nonlinearDiag}) and the time-domain RLS in parallel configuration. At first glance, it can be presumed that with perfect decoding, the performance is no longer dependent on the input SINR. However, this is only the case for the RLS algorithm in static environments. In the low SINR regime, the SI power is strong and thus the approximated algorithm is rather limited by estimation errors or favors certain levels of SI power.  Moreover, the probability of decoding errors is much higher in that regime and thus the potential rate gain is in fact even lower. This leads to the conclusion that in low SINR regime for static environments, the decoding step is of less relevance. However, the picture is different if the environment is time-variant, as illustrated in~\figref{fig:DecRateNonOrthoVarying}. The primary limiting factor is the temporal variations. The Kalman algorithm benefits twofold from the perfect decoding: First, the SoI is removed from the received signal before the SI channel update, and second, the signal covariance is adjusted to the noise floor. This improves adaptation to the time-variant context and makes it possible for the Kalman algorithm to be uniformly superior over the RLS in both low and high input SINR regimes. The actual rate performance of decoding with partial errors is expected to lie somewhere in between the two extreme cases of either perfect or no decoding at all. Thus,~\figref{fig:DecRateNonOrthoVarying} provide insights on the lower and upper performance bounds. Decoding increases the rate performance especially in moderate-to-high input SINR regimes.

\section{Conclusion}
\label{sec:conclusion}
In this work, we have proposed an adaptive SI cancellation algorithm for full-duplex communication in the digital domain. It is based on a state-space model of the nonlinear SI channel in cascade structure. The signal model is derived in DFT domain. Linear and nonlinear components, although statistically coupled by the system model, can be linearized and then estimated separately in a two-step iterative approach, each by a Kalman filter. The algorithm comprises input orthogonalization, prediction, decoding of the desired signal, and updating the estimations.

The state-space representation has introduced a more fine-grained handling of the time-variant nature of the SI channel. It underlines the significance of a-priori knowledge on the scale of temporal variations, which is represented by the state-space system model. Our simulation results indicate the following conclusion. While the choice of the underlying algorithmic principle like NLMS, RLS or Kalman do influence the performance of SI cancellation, the structural elements such as orthogonalization, parallel or cascade model, and the decoding have a much more significant impact. Especially the orthogonalization equalizes the performance gain for all algorithms. This work also demonstrates that the temporal variations are determining the fundamental performance limitation.