Modeling and Identification of Nonlinear Effects in Massive MIMO Systems Using a Fifth-Order Cumulants-Based Blind Approach

Abstract

Pre-processing associated with massive multiple input-multiple output (MIMO) systems can lead to signals with high envelope fluctuations, which are very prone to nonlinear effects, especially when massive MIMO schemes are combined with orthogonal transform multiplexing (OFDM) modulations. If the nonlinear characteristics that affect a given system are known, we can design appropriate receivers that take into account the nonlinear effects introduced by the transmitter. Cubic systems are particularly important, not only because they can approximate many nonlinear effects (e.g., due to the power amplifier or clipping effects), but also because many more complex nonlinear characteristics in communication schemes can be replaced by equivalent lower-order nonlinear characteristics in general, and cubic characteristics in particular. To compensate the effects at the receiver side (e.g., by using the so-called Bussgang receivers), we need to estimate the nonlinear operation that was introduced at the transmitter, and this should be done blindly, without the need of training symbols. The paper contains a description of a mathematical approach for modeling and identification of nonlinear kernels in cubic systems. Based on theoretical tools of HOC in cubic systems, we build a new formula which relates the second- and fifth-order cumulants. Our performance results indicate that the proposed approach allows an accurate identification, yielding the desired kernels via fifth-order cumulants, and ensures a very good convergence, outperforming existing adaptive methods. This is achieved blindly, by exploiting the maximum information of the output system, making it suitable for many practical nonlinear effects.

1. Introduction

In digital communication, strong nonlinear effects throughout the transmission chain can have a negative impact on the overall system’s performance [1]. Pre-processing associated with massive MIMO systems can lead to signals with envelope fluctuations, which are very prone to nonlinear effects, especially when massive MIMO schemes are combined with OFDM modulations [1,2,3,4,5,6,7,8,9,10]. Several works focused on MIMO systems and their important parameters such as mutual coupling [11,12,13,14]. Pre-processing techniques employed at the transmitter side of MIMO systems, such as the ones used by MIMO singular value decomposition (MIMO-SVD) techniques, can lead to signals with large envelope fluctuations and high peak to-average power ratio (PAPR) [3]. The same occurs with the Zero Forcing Transmitter (ZFT), Maximum Ratio Transmitter (MRT), and Equal Gain Transmitter (EGT), when used in m-MIMO implementation with pre- or post-processing [4]. On the other hand, several authors have qualified nonorthogonal multiple access (NOMA) as the main candidate to support fifth-generation (5G) wireless communication. This has been investigated in several works (see, e.g., [15,16,17]). By combining signals with different power levels, NOMA schemes can also be very prone to nonlinear effects.

If the nonlinear characteristics that affect a given system are known, we can design appropriate receivers that take into account the nonlinear effects introduced by the transmitter. In fact, efficient receivers that jointly perform detection and estimation and cancellation of nonlinear distortion effects, can have excellent performance/complexity trade-offs [2,8]. On the other hand, taking into account nonlinear distortion effects in the receiver design can have significant performance improvements [18], sometimes even outperforming the linear case [10,18,19,20]. As an alternative, we can employ receivers such as the Bussgang receiver, which estimate and cancel nonlinear distortion [21,22,23]. All these techniques require the knowledge of the nonlinear characteristic used by the transmitter, which can be obtained blindly using higher-order cumulants (HOC) theory. We propose to study the issue of nonlinear system identification using the fifth-order cumulants domain.

One of the first uses of HOC in signal processing has been devoted to the problem of blind system identification. These systems, which can be linear or nonlinear, are excited by an unobservable, zero-mean, independent and identically distributed (i.i.d), stationary, ergodic and non-Gaussian signal under the presence of Gaussian noise. Linear systems identification has been very deeply investigated in the literature (see, e.g., [24,25,26,27,28]). The primary existing blind approaches are based on second-order cumulants (usually called autocorrelation function). The second- and third-order cumulants were used by Alshebeili et al. to identify the nonminimum phase (NMP) of finite impulse response (FIR) systems [24]. Likewise, Giannakis et al. [25] presented an identification method by using both second- and higher-order statistics of an NMP output system and many approaches devoted to the same issue have been reported in the literature (see, e.g., [26]). However, these second-order-based methods are not able to identify correctly the linear system’s model and their performance was degraded when the output is affected by additive Gaussian noise. To overcome these inconveniences, other linear algebraic methods involving only HOC of output system have been reported in the literature (see, e.g., [29,30,31,32,33,34]), which present several advantages such as effectiveness to reject additive Gaussian noise and ability to reconstruct the phase. Other authors showed, both analytically and by simulations, that better results can be obtained using these approaches. The case when the noise is non-Gaussian was studied deeply by Sadler et al. by using HOC tools [35]. A binary random signal was also used in the issue of blind identification of linear systems using symbolic dynamics [36,37].

Most of these existing blind algorithms are focused on linear system models due of their simplicity. However, many practical systems are nonlinear, and these nonlinearities cannot be ignored without serious performance degradation. This applies to many practical fields, including signal processing, prediction, communication channels, biomedicine, control of industrial processes, detection, equalization and many others [38,39,40,41].

Blind nonlinear system identification has received great attention in the literature (see, e.g., [42,43,44,45]). In the nonlinear domain, discrete Volterra representations were used to describe the input–output relation, thus that characterization, analysis and synthesis are easily amenable [42]. Stathaki et al. [43] proposed to identify the Volterra kernels by using the third cumulant with a constraint of the second-order cumulant. It has been demonstrated that this constrained optimization approach can achieve excellent results for nonlinear system identification. An inconvenience of this method is its sensitivity to additive Gaussian noise since it is based on second-order cumulant, as pointed out in [44]. To overcome this issue, Chen et al. [44] used a fourth-order-based method link as alternative solution of the second-order cumulant to smooth out the additive Gaussian noise and improve the convergence by proposing a hybrid optimization method of simplex and genetic algorithms. In the same frame, Glentis et al. [38] contributed to Volterra system identification by developing an efficient adaptive transversal algorithm, with the advantages of low complexity and modular structure. There are also many existing works in the context of nonlinear Volterra system identification (see, e.g., [46,47,48]). In [47], the authors used the third-order moment (TOM) domain for blind identification of nonlinear Volterra systems with studies on the computational complexity and demonstrated via simulation that the proposed approach can be applied in a practical scenario to identify blindly the primary path of active-noise-control. A class of truncated sparse Volterra systems was identified blindly by using the second-order moment (SOM) or TOM sequences of the system outputs in [48]. The authors of [42] investigated the third-order Volterra nonlinear systems with a novel approach based only second-order statistics and proved that was able to remove white noise with any distribution offering significant reduction in the computational burden. A complex random sequence input coupled with output HOC was employed in [40] to identify blindly the linear quadratic Volterra systems by examining different cases and the simulations were also used to illustrate the performance of the proposed approaches. The authors in [45] proposed an extension of the main concept behind the error whitening criterion (EWC) in the linear case to the unbiased identification by analytically examining the true parameters of order-2 Volterra series models of nonlinear dynamical systems. It was emphasized by the authors that this extension does not apply to higher-order Volterra systems. A hydroturbine shaft system modeled as a quadratic Volterra system has been identified blindly using third-order cumulants and a reversely recursive method by Bai and Zhang in [46]. The proposed method was tested using engineering applications and three numerical experiments, and their applicability was demonstrated by the blind identification of the hydroturbine shaft system. Furthermore, identification of the quadratic nonlinear systems has been reported in the past (see, e.g., [49,50,51,52]). In [49], the properties of higher-order moment sequences and the calculation of the higher-order spectral moments of an i.i.d for various frequencies were used to develop algorithms to identify time-invariant quadratic nonlinear systems. A concrete application of higher-order statistics theory with more experimental discussion was developed by Chow et al. [50], to identify machine faults through vibration measurements using higher-order statistics-based methods, such as a nonparametric phase-bearing and a parametric linear or nonlinear modeling approach. Another application was developed in [51] to model the delay of video-packets transmission by identification of quadratic nonlinear systems using fourth-order cumulants-based blind approaches, which were compared to the Levenberg–Marquardt algorithm, least mean square (LMS) and recursive least square (RLS). A HOC was also used by Zidane et al. [52] to develop both extending approaches, developed for the linear and nonlinear case, which uses fourth-order only and combined third- and fourth-order cumulants, to identify the diagonal parameters of nonlinear quadratic systems. Numerical simulations were presented to verify the theoretical results with comparison to the exiting method in terms of the normalized mean-squared error (NMSE) and the fluctuations around the true kernels. Motivated by reductions in computational requirements and the mathematical tractability of nonlinear system identification issue, Ralston et al. [41] used a Hammerstein series to build a new method to identify a specific time invariant nonlinear system when the input is a non-Gaussian stationary signal.

As previously indicated, several authors considered the issue of nonlinear system identification using HOC theory, but most works focused on second-order or quadratic system identification, and the order of cumulants was limited to fourth-order cumulants. The fifth-order received little attention up to now. Indeed, crosscumulants up to the fifth-order were used by Koukoulas et al. [39] to identify the input–output quadratic system. Cubic system identification has also received so far little attention. Indeed, the Volterra kernels of the cubic system identification method using higher-order moments were developed by Tseng and Powers in [53]. The authors showed that the proposed method reduces the complexity of Volterra kernel identification compared to the non-Gaussian and non-i.i.d. input case. The authors of [54] used the TOM from the extending SOM to identify the nonlinear cubic system. It was shown that more statistical knowledge can be obtained in the third-order statistics domain for blind system identification. In [55], mixed third- and fourth-order cumulants were used to extend linear HOC-based methods to nonlinear systems.

In general, when we increase the nonlinearity degree in the model, the performances of methods relying on cumulants with small orders can degrade substantially in the presence of Gaussian noise. Therefore, we increase the order of the proposed approach by considering the fifth-order. In this framework, we consider cubic nonlinear systems, which can be used to describe a wide range of nonlinear systems, requiring only a small number of kernels in the Volterra series. The main contributions of this work are the development of new fifth-order cumulants which exploit the maximum information of the output cubic system without the need of training signals, i.e., blindly. The methodology of the proposed blind approach is based on the main theoretical tools of HOC using the link between Fourier transform of second-order and fifth-order cumulants in the case of cubic systems which are established to build a new formula which relates the second- and fifth-order cumulants by using the inverse Fourier transform. Based on this formula, which relates the second- and fifth-order cumuants, we can propose an approach based only on fifth-order under some assumptions for the kernels of systems and the properties of the autocorrelation function. To analyze the convergence and test the effectiveness of the proposed blind approach, a set of simulations is carried out for various signal to noise ratio (SNR) levels and the proposed approach is compared with adaptive existing algorithms [56] that use the input–output relations to identify the kernels of cubic systems.

The outline of this paper is organized as follows. We start in the Section 2 with the problem formulation and some assumptions regarding the system model representation. Section 3 describes the main theoretical tools of HOC in the case of cubic systems. Theoretical development based using polyspectra are provided in the Section 4 to propose a blind approach based only on fifth-order cumulants. A set of performance results is supplied to support the theoretical development of the proposed blind approach, which is compared with an adaptive algorithm [56] in Section 5.

2. Mathematical Definitions

The characteristic function of the vector $X={(x_1,x_2,\dots ,x_k)}^T$ composed of k real random variables $x_i$ is defined by

$${\Phi }_X\left(V\right)=E\left[exp\left(j{V}^{T}X\right)\right],$$

(1)

where $V^T=(v_1,v_2,\dots ,v_k)$ is vector composed of k random variables. $V^{T}X=v_1{x}_1+v_2{x}_2+\dots +v_k{x}_k$ represents the scalar product of V and X and $E[.]$ is the mathematical expectation operator.

The coefficients of the kth-order moments are defined in the Taylor series development of the characteristic function ${\Phi }_X\left(V\right)$ as

$$\begin{array}{ccc}\hfill {\Phi }_X\left(V\right)& =& \sum _{0\le m\le s}\frac{j^m}{n_1!n_2!\dots n_k!}Mom(x_1^{n_1},x_2^{n_2},\dots ,x_k^{n_k})\hfill \\ & \times & v_1^{n_1}{v}_2^{n_2}\dots v_k^{n_k}+{o\left(\right|V|}^s),\hfill \end{array}$$

(2)

where

$$\begin{array}{ccc}\hfill Mom(x_1^{n_1},x_2^{n_2},\dots ,x_k^{n_k})& =& j^{-m}[{\left(\frac{\partial }{\partial v_1}\right)}^{n_1}{\left(\frac{\partial }{\partial v_2}\right)}^{n_2}\hfill \\ & \dots & {\left(\frac{\partial }{\partial v_k}\right)}^{n_k}{\Phi }_X\left(V\right){]}_{V=0}\hfill \end{array}$$

(3)

The kth-order moments of these random variables are given by

$$Mom(x_1^{n_1},x_2^{n_2},\dots ,x_k^{n_k})=E[x_1^{n_1}{x}_2^{n_2}\dots x_k^{n_k}],$$

(4)

The cumulant-generating function is defined as logarithm of the first characteristic function

$${\Psi }_X\left(V\right)=ln\left({\Phi }_X\left(V\right)\right).$$

(5)

The mth-order cumulant of these random variables are defined as the coefficient of $(v_1,v_2,\dots ,v_m)$ in the Taylor series development of the cumulant-generating function

$$\begin{array}{ccc}\hfill {\Psi }_X\left(V\right)& =& \sum _{0\le m\le s}\frac{j^m}{n_1!n_2!\dots n_k!}Cum(x_1^{n_1},x_2^{n_2},\dots ,x_k^{n_k})\hfill \\ & \times & v_1^{n_1}{v}_2^{n_2}\dots v_k^{n_k}+{o\left(\right|V|}^s),\hfill \end{array}$$

(6)

where

$$\begin{array}{ccc}\hfill Cum(x_1^{n_1},x_2^{n_2},\dots ,x_k^{n_k})& =& j^{-m}[{\left(\frac{\partial }{\partial v_1}\right)}^{n_1}{\left(\frac{\partial }{\partial v_2}\right)}^{n_2}\hfill \\ & \dots & {\left(\frac{\partial }{\partial v_k}\right)}^{n_k}{\Psi }_X\left(V\right){]}_{V=0}\hfill \end{array}$$

(7)

The operator $Cum[·]$ stands for the nth-order joint cumulant of the random variables $x_1^{n_1},x_2^{n_2},\dots ,x_k^{n_k}$. Thus, considering a zero-mean random process $\left\{y\right(t\left)\right\}$, we have

$$\begin{array}{c}\hfill C_{n,\Psi }=Cum\left[y\left(t\right),y(t+{\tau }_1),\dots ,y(t+{\tau }_{n-1})\right]\end{array}$$

(8)

where, $\Psi =\left\{y\left(t\right),y(t+{\tau }_1),\dots ,y(t+{\tau }_{n-1})\right\}$, and, hence, the involved real random variables are time-shifted samples of the process $\left\{y\right(t\left)\right\}$. Notice that, for stationary processes, the nth-order statistics depend only on the $n-1$ time lags ${\tau }_1=t_1-t_0,\dots ,{\tau }_{n-1}=t_{n-1}-t_0$. This allows us to introduce the following notations:

$$\begin{array}{c}\hfill C_{n,y}({\tau }_1,{\tau }_2,\dots ,{\tau }_{n-1})=Cum\left[y\left(t\right),y(t+{\tau }_1),\dots ,y(t+{\tau }_{n-1})\right]\end{array}$$

(9)

where $\left\{y\right(t\left)\right\}$ is a zero-mean stationary random process.

3. Problem

Cubic systems can approximate many nonlinear effects (e.g., due to the power amplifier or quantization effects), and we find that more complex nonlinear characteristics can be replaced by equivalent lower-order nonlinear characteristics in general, and cubic characteristics in particular [21]. Indeed, Figure 1 describes the block diagram of the nonlinear system identification problem which is represented by a Volterra series with only a small number of kernels. The system, which we propose to identify in this work, is a cubic nonlinearity that has the form

$$\begin{array}{c}\left\{\begin{array}{c}{y}_0\left(k\right)={\displaystyle \sum _{i=0}^q}h(i,i,i)x^3(k-i)\hfill \\ \\ y\left(k\right)=y_0\left(k\right)+n_G\left(k\right)\hfill \end{array}\right.\end{array}$$

(10)

where $\left\{y_0\left(k\right)\right\}$ is the cubic system output, and $\left\{y\right(k\left)\right\}$ is the observation output, which is contaminated by a zero-mean additive Gaussian noise $\left\{n_G\left(k\right)\right\}$, assumed independent of $\left\{y_0\left(k\right)\right\}$, excluding the effect of the noise where nth-order cumulants are superior of 2 (i.e., $C_{n,y_0}({\tau }_1,{\tau }_2,\dots ,{\tau }_n)=C_{n,y}({\tau }_1,{\tau }_2,\dots ,{\tau }_n)$) as indicated above. The input random signal, $\left\{x\right(k\left)\right\}$, of the model is the unobservable zero mean, independent identically distributed (i.i.d) stationary, ergodic and non-Gaussian with known distribution statistics ${\gamma }_{3,x}=E\left[x^3\left(k\right)\right]\ne 0$, ${\gamma }_{6,x}=E\left[x^6\left(k\right)\right]\ne 0$, ${\gamma }_{9,x}=E\left[x^9\left(k\right)\right]\ne 0$, ${\gamma }_{12,x}=E\left[x^{12}\left(k\right)\right]\ne 0$ and ${\gamma }_{15,x}=E\left[x^{15}\left(k\right)\right]\ne 0$. Unknown cubic kernels include of $\left\{h\right(i,i,i)\forall i=1,\dots ,q\}$, and ${\displaystyle \sum _i}\left|h(i,i,i)\right|<\infty$, which suggest that the system is stable causal implying bounded input–output.

The cubic system identification is more complex and we need even more information to calculate the cumulants of the generated output signal. The main objective of this investigation is to identify the kernels of the cubic nonlinear system in (10) by using the blind approach based on fifth-order cumulants of the system’s measured output, $C_{5,y}({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)$, and some knowledge of the properties of the input random signal $\left\{x\right(k\left)\right\}$.

4. Theoretical Tools of HOC

In this section, we focus on the theoretical development of HOC, which present the main relationships linking diagonal kernels of cubic systems and the cumulants of the output that received random signals up to the fifth-order for the purpose to use the HOC to identify the cubic system described by (10).

The starting point for all linear-quadratic blind methods is the Leonov–Shiryaev formula [57,58], which links the different order of cumulants to the moments. The latter formula allows expressing any cumulant as a function of moments of lower or equal orders and is well-known in the case of non-delayed cumulants [59]. It is assumed above that the output signal is stationary and ergodic and these cumulants do not depend on time but on the difference time between the instant of observation.

The rth-order cumulants of the output random signal are linked to the moments, where the order p is inferior or equal to r, by the following formula of Leonov and Shiryayev [57,58]:

$$\begin{array}{ccc}\hfill Cum[y_1,\dots ,y_r]& =& \sum {(-1)}^{k-1}(k-1)!E\left[\prod _{i\in v_1}{y}_i\right]E\left[\prod _{j\in v_2}{y}_j\right]\dots E\left[\prod _{k\in v_p}{y}_k\right],\hfill \end{array}$$

(11)

The summation extends over all ensembles $\{v_1,v_2,\dots ,v_p:1\le p\le r\}$ forming a partition $\{1,2,\dots ,r\}$. In this formula, k is the number of elements that compose the partition.

For the second-order terms, both partitions are possible (1, 2) and (1) (2). Thus, we have the following expression:

$$\begin{array}{ccc}\hfill Cum[y_1,y_2]& =& {(-1)}^0\left(0\right)!E\left[y_1{y}_2\right]+{(-1)}^{2-1}(2-1)!E\left[y_1\right]E\left[y_2\right]\hfill \\ & =& E\left[y_1{y}_2\right]-E\left[y_1\right]E\left[y_2\right].\hfill \end{array}$$

(12)

For the stationary random output signal, $\left\{y\right(t\left)\right\}$, its second-order cumulants, in the cubic system case, become

$$\begin{array}{ccc}\hfill C_{2,y}\left({\tau }_1\right)& =& Cum[y_1,y_2]=Cum[y\left(k\right),y(k+{\tau }_1)]\hfill \\ & =& E[\sum _{i=0}^{q}h(i,i,i)x^3(k-i)\sum _{j=0}^{q}h(j,j,j)\hfill \\ & \times & x^3(k+{\tau }_1-j)]-E\left[\sum _{i=0}^{q}h(i,i,i)x^3(k-i)\right]\hfill \\ & \times & E\left[\sum _{j=0}^{q}h(j,j,j)x^3(k+{\tau }_1-j)\right]\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill C_{2,y}\left({\tau }_1\right)& =& E[\sum _{i=0}^{q}h(i,i,i)x^3(k-i)h(i+{\tau }_1,i+{\tau }_1,\hfill \\ & & i+{\tau }_1)x^3(k-i)]-E\left[\sum _{i=0}^{q}h(i,i,i)x^3(k-i)\right]\hfill \\ & \times & E\left[\sum _{i=0}^{q}h(i+{\tau }_1,i+{\tau }_1,i+{\tau }_1)x^3(k-i)\right]\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill C_{2,y}\left({\tau }_1\right)& =& \sum _{i=0}^{q}h(i,i,i)h(i+{\tau }_1,i+{\tau }_1,i+{\tau }_1)\hfill \\ & \times & (E\left[x^3(k-i)\times x^3(k-i)\right]\hfill \\ & -& E\left[x^3(k-i)\right]E\left[x^3(k-i)\right]).\hfill \end{array}$$

(15)

Under the assumption that the input sequence $x\left(k\right)$ is i.i.d zero mean, stationary, non-Gaussian with ${\gamma }_{n,x}=E\left[x^n\left(k\right)\right]\ne 0$, $\forall n=3,6,9,12,15$ and (15), the second-order cumulants and the diagonal kernels of cubic systems are linked by the following expression

$$C_{2,y}\left({\tau }_1\right)=\left({\gamma }_{6,x}-{\gamma }_{3,x}^2\right)\sum _{i=0}^{q}h(i,i,i)h(i+{\tau }_1,i+{\tau }_1,i+{\tau }_1).$$

(16)

In the case of fifth-order cumulants, we have 52 possible partitions of 7 different types:

• 1 partition of type: (1, 2, 3, 4, 5) with k = 1;
• 5 partitions of type: (1) (2, 3, 4, 5) with k = 2;
• 10 partitions of type: (1, 2) (3, 4, 5) with k = 2;
• 10 partitions of type: (1) (2) (3, 4, 5) with k = 3;
• 15 partitions of type: (1) (2, 3) (4, 5) with k = 3;
• 10 partitions of type: (1) (2) (3) (4, 5) with k = 4;
• 1 partition of type: (1) (2) (3) (4) (5) with k = 5.

$$\begin{array}{ccc}\hfill Cum\left[y_1,y_2,y_3,y_4,y_5\right]& =& E\left[y_1{y}_2{y}_3{y}_4{y}_5\right]\hfill \\ & -& \left[5\right]E\left[y_1\right]E\left[y_2{y}_3{y}_4{y}_5\right]\hfill \\ & -& \left[10\right]E\left[y_1{y}_2\right]E\left[y_3{y}_4{y}_5\right]\hfill \\ & +& \left[10\right]2E\left[y_1\right]E\left[y_2\right]E\left[y_3{y}_4{y}_5\right]\hfill \\ & +& \left[15\right]2E\left[y_1\right]E\left[y_2{y}_3\right]E\left[y_4{y}_5\right]\hfill \\ & -& \left[10\right]6E\left[y_1\right]E\left[y_2\right]E\left[y_3\right]E\left[y_4{y}_5\right]\hfill \\ & +& 24E\left[y_1\right]E\left[y_2\right]E\left[y_3\right]\hfill \\ & \times & E\left[y_4\right]E\left[y_5\right].\hfill \end{array}$$

(17)

The fifth-order cumulants which relate the diagonal kernels of cubic systems can be expressed as (see Appendix A):

$$\begin{array}{ccc}\hfill C_{5,y}({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)& =& ({\gamma }_{15,x}-5{\gamma }_{3,x}{\gamma }_{12,x}-10{\gamma }_{6,x}{\gamma }_{9,x}\hfill \\ & +& 20{\gamma }_{3,x}^2{\gamma }_{9,x}+30{\gamma }_{3,x}{\gamma }_{6,x}^2-60{\gamma }_{3,x}^3{\gamma }_{6,x}\hfill \\ & +& 24{\gamma }_{3,x}^5)\sum _{i=0}^{q}h(i,i,i)h(i+{\tau }_1,i+{\tau }_1,\hfill \\ & & i+{\tau }_1)h(i+{\tau }_2,i+{\tau }_2,i+{\tau }_2)\hfill \\ & \times & h(i+{\tau }_3,i+{\tau }_3,i+{\tau }_3)\hfill \\ & \times & h(i+{\tau }_4,i+{\tau }_4,i+{\tau }_4)\hfill \end{array}$$

(18)

5. Proposed Blind Approach

The nth-order polyspectrum is defined as the $(n-1)$-dimensional discrete-time Fourier transform of the nth-order cumulants $C_{n,y}({\tau }_1,{\tau }_2,\dots ,{\tau }_n)$, i.e.,

$$\begin{array}{ccc}\hfill S_{n,y}({\omega }_1,{\omega }_2,\dots ,{\omega }_{n-1})& =& \sum _{{\tau }_1=-\infty }^{+\infty }\cdots \sum _{{\tau }_{n-1}=-\infty }^{+\infty }{C}_{n,y}({\tau }_1,{\tau }_2,\hfill \\ & \dots & ,{\tau }_{n-1})exp\left(-j\sum _{i=1}^{n-1}{\omega }_i{\tau }_i\right).\hfill \end{array}$$

(19)

The second-order polyspectrum, Fourier transform of the (16), (called power spectrum) is linked to the cubic system transfer function $H(\omega ,\omega ,\omega )$ by

$$\begin{array}{ccc}\hfill S_{2,y}\left(\omega \right)& =& T.F\left\{C_{2,y}\left(\tau \right)\right\}\hfill \\ & =& {\varphi }_{2}H(\omega ,\omega ,\omega )H(-\omega ,-\omega ,-\omega ),\hfill \end{array}$$

(20)

where, ${\varphi }_2={\gamma }_{6,x}-{\gamma }_{3,x}^2$.

The fifth-order polyspectrum, Fourier transform of Equation (18), of the cubic system output is linked to the system transfer function by

$$\begin{array}{ccc}\hfill S_{5,y}({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4)& =& T.F\left\{C_{5,y}({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)\right\}\hfill \\ & =& {\varphi }_{5}H\left({\omega }_1,{\omega }_1,{\omega }_1\right)H\left({\omega }_2,{\omega }_2,{\omega }_2\right)\hfill \\ & \times & H\left({\omega }_3,{\omega }_3,{\omega }_3\right)H\left({\omega }_4,{\omega }_4,{\omega }_4\right)\hfill \\ & \times & H\left(-\sum _{i=1}^4{\omega }_i,-\sum _{i=1}^4{\omega }_i,-\sum _{i=1}^4{\omega }_i\right),\hfill \end{array}$$

(21)

where ${\varphi }_5={\gamma }_{15,x}-5{\gamma }_{3,x}{\gamma }_{12,x}-10{\gamma }_{6,x}{\gamma }_{9,x}+20{\gamma }_{3,x}^2{\gamma }_{9,x}+30{\gamma }_{3,x}{\gamma }_{6,x}^2-60{\gamma }_{3,x}^3{\gamma }_{6,x}+24{\gamma }_{3,x}^5$.

Defining $\omega ={\displaystyle \sum _{i=1}^4}{\omega }_i$, (20) becomes

$$\begin{array}{ccc}\hfill S_{2,y}\left(\sum _{i=1}^4{\omega }_i\right)& =& {\varphi }_{2}H\left(\sum _{i=1}^4{\omega }_i,\sum _{i=1}^4{\omega }_i,\sum _{i=1}^4{\omega }_i\right)\hfill \\ & \times & H\left(-\sum _{i=1}^4{\omega }_i,-\sum _{i=1}^4{\omega }_i,-\sum _{i=1}^4{\omega }_i\right).\hfill \end{array}$$

(22)

From (21) and (22), we obtain

$$\begin{array}{c}\hfill S_{2,y}\left(\sum _{i=1}^4{\omega }_i\right)=\frac{{\varphi }_2}{{\varphi }_5}\frac{S_{5,y}({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4)H\left({\displaystyle \sum _{i=1}^4}{\omega }_i,{\displaystyle \sum _{i=1}^4}{\omega }_i,{\displaystyle \sum _{i=1}^4}{\omega }_i\right)}{H\left({\omega }_1,{\omega }_1,{\omega }_1\right)H\left({\omega }_2,{\omega }_2,{\omega }_2\right)H\left({\omega }_3,{\omega }_3,{\omega }_3\right)H\left({\omega }_4,{\omega }_4,{\omega }_4\right)}\end{array}$$

(23)

$$\begin{array}{ccc}& & S_{5,y}({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4)H\left(\sum _{i=1}^4{\omega }_i,\sum _{i=1}^4{\omega }_i,\sum _{i=1}^4{\omega }_i\right)=\hfill \\ & & {\eta }_{5,2}H\left({\omega }_1,{\omega }_1,{\omega }_1\right)H\left({\omega }_2,{\omega }_2,{\omega }_2\right)H\left({\omega }_3,{\omega }_3,{\omega }_3\right)H\left({\omega }_4,{\omega }_4,{\omega }_4\right)S_{2,y}\left(\sum _{i=1}^4{\omega }_i\right),\hfill \end{array}$$

(24)

where, ${\eta }_{5,2}=\frac{{\varphi }_5}{{\varphi }_2}=\frac{{\gamma }_{15,x}-5{\gamma }_{3,x}{\gamma }_{12,x}-10{\gamma }_{6,x}{\gamma }_{9,x}+20{\gamma }_{3,x}^2{\gamma }_{9,x}+30{\gamma }_{3,x}{\gamma }_{6,x}^2-60{\gamma }_{3,x}^3{\gamma }_{6,x}+24{\gamma }_{3,x}^5}{{\gamma }_{6,x}-{\gamma }_{3,x}^2}$.

Then, we use the inverse Fourier transform of (24) to demonstrate that the second-order, the fifth-order cumulants and the kernels of cubic systems are linked by

$$\begin{array}{ccc}& & \sum _{i=0}^{q}h(i,i,i)C_{5,y}({\tau }_1-i,{\tau }_2-i,{\tau }_3-i,{\tau }_4-i)\hfill \\ & & ={\eta }_{5,2}\sum _{i=1}^{q}h(i,i,i)h({\tau }_2-{\tau }_1+i,{\tau }_2-{\tau }_1+i,{\tau }_2-{\tau }_1+i)\hfill \\ & & \times h({\tau }_3-{\tau }_1+i,{\tau }_3-{\tau }_1+i,{\tau }_3-{\tau }_1+i)\times \hfill \\ & & \times h({\tau }_4-{\tau }_1+i,{\tau }_4-{\tau }_1+i,{\tau }_4-{\tau }_1+i)C_{2,y}({\tau }_1-i).\hfill \end{array}$$

(25)

Based on this, we can develop the following algorithm to determine the unknown kernels in the cubic system case. We start by taking ${\tau }_1=2q$, and we use the properties of the second-order cumulants of the stationary, ergodic systems which are $C_{2,y}\left(\tau \right)\ne 0$ for $-q\le \tau \le +q$, and zero elsewhere. Therefore, (25) yields

$$\begin{array}{ccc}& & \sum _{i=0}^{q}h(i,i,i)C_{5,y}(2q-i,{\tau }_2-i,{\tau }_3-i,{\tau }_4-i)=\hfill \\ & & {\eta }_{5,2}h(q,q,q)h({\tau }_2-q,{\tau }_2-q,{\tau }_2-q)\hfill \\ & & \times h({\tau }_3-q,{\tau }_3-q,{\tau }_3-q)\hfill \\ & & \times h({\tau }_4-q,{\tau }_4-q,{\tau }_4-q)C_{2,y}\left(q\right),\hfill \end{array}$$

(26)

The assumed system is causal and truncated, i.e., $h(i,i,i)=0$ if $i<0$ and $i>q$ with $h(0,0,0)=1$. Therefore, ${\tau }_2,{\tau }_3,{\tau }_4$ take values q and $2q$. If we take ${\tau }_2={\tau }_3=q$ in (26), we have

$$\begin{array}{ccc}& & \sum _{i=0}^{q}h(i,i,i)C_{5,y}(2q-i,q-i,q-i,{\tau }_4-i)\hfill \\ & & ={\eta }_{5,2}h(q,q,q)h^2(0,0,0)\hfill \\ & & \times h({\tau }_4-q,{\tau }_4-q,{\tau }_4-q)C_{2,y}\left(q\right),\hfill \end{array}$$

(27)

where $h(0,0,0)=1$.

Else if, ${\tau }_4=2q$, (27) becomes:

$$\begin{array}{c}\hfill h(q,q,q)C_{5,y}(q,0,0,q)={\eta }_{5,2}{h}^2(q,q,q)C_{2,y}\left(q\right),\end{array}$$

(28)

Using (27) and (28), we obtain the following expression which links the kernels of cubic system and fifth-order cumulants only:

$$\begin{array}{ccc}& & \sum _{i=0}^{q}h(i,i,i)C_{5,y}(2q-i,q-i,q-i,{\tau }_4-i)\hfill \\ & =& C_{5,y}(q,0,0,q)h({\tau }_4-q,{\tau }_4-q,{\tau }_4-q),\hfill \end{array}$$

(29)

where

$$q\le {\tau }_4\le 2q.$$

(30)

Thus, (29) and (30) are combined in the following matrix form

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}{C}_{5,y}(2q-1,q-1,q-1,q-1)& \dots & C_{5,y}(q,0,0,0)\\ C_{5,y}(2q-1,q-1,q-1,q)-\psi & \dots & C_{5,y}(q,0,0,1)\\ ⋮& \ddots & ⋮\\ C_{5,y}(2q-1,q-1,q-1,2q-1)& \dots & C_{5,y}(q,0,0,q)-\psi \end{array}\right]\\ \hfill \times \left[\begin{array}{c}h(1,1,1)\\ .\\ .\\ .\\ h(i,i,i)\\ .\\ .\\ .\\ h(q,q,q)\end{array}\right]=\left[\begin{array}{c}\psi -C_{5,y}(2q,q,q,q)\\ -C_{5,y}(2q,q,q,q+1)\\ .\\ .\\ .\\ .\\ .\\ .\\ -C_{5,y}(2q,q,q,2q)\end{array}\right],\end{array}$$

(31)

where, $\psi =C_{5,y}(q,0,0,q)$.

Equation (31) can be rewritten in other compact form as

$$MH=D,$$

(32)

where, $H={[h(1,1,1),\dots ,h(q,q,q)]}^T$ is a column vector of kernels of the cubic system with size q, $D={[\psi -C_{5,y}(2q,q,q,q),-C_{5,y}(2q,q,q,q+1),\dots ,-C_{5,y}(2q,q,q,2q)]}^T$ is a column vector of size $(q+1)$ and M the matrix of size $(q+1)\times q$ elements. The least square solution of (32) is expressed by

$$\widehat{H}={\left(M^{T}M\right)}^{-1}{M}^{T}D$$

(33)

6. Performance Results

To study the effectiveness and analyze the convergence of the proposed blind approach under various signal-to-noise levels (SNR), defined as $\mathrm{SNR}=E\left[y_0^2\left(k\right)\right]/E\left[n_G^2\left(k\right)\right]$. We consider 200 Monte Carlo runs and use two cubic systems.

To measure the precision of cubic kernel identification with respect to the true values, we define the NMSE for each run as

$$\mathrm{NMSE}=\frac{{\displaystyle \sum _{i=0}^q}{\left(h(i,i,i)-\widehat{h}(i,i,i)\right)}^2}{{\displaystyle \sum _{i=0}^q}{h}^2(i,i,i)},$$

(34)

where $\widehat{h}(i,i,i)$ are the estimated kernels in each run, and $h(i,i,i)$ are the true kernels in the cubic system.

System I: The third-order cubic system is

$$\begin{array}{c}\left\{\begin{array}{c}{y}_0\left(k\right)=x^3\left(k\right)-0.15x^3(k-1)-0.35x^3(k-2)\hfill \\ +0.90x^3(k-3)\hfill \\ \\ y\left(k\right)=y_0\left(k\right)+n_G\left(k\right)\hfill \end{array}\right.\end{array}$$

(35)

With reference to (10), in this case, we have $h(0,0,0)=1$, $h(1,1,1)=-0.15$, $h(2,2,2)=-0.35$ and $h(3,3,3)=0.90$.

Then, we try to increase the system order to test robustness of our approach. We consider the fifth-order cubic system given by

System II:

$$\begin{array}{c}\left\{\begin{array}{c}{y}_0\left(k\right)=x^3\left(k\right)-0.25x^3(k-1)+0.36x^3(k-2)\hfill \\ -0.688x^3(k-3)+0.56x^3(k-4)\hfill \\ +0.227x^3(k-5)\hfill \\ \\ y\left(k\right)=y_0\left(k\right)+n_G\left(k\right)\hfill \end{array}\right.\end{array}$$

(36)

With reference to (10), in this case, we have $h(0,0,0)=1$, $h(1,1,1)=-0.25$, $h(2,2,2)=0.36$, $h(3,3,3)=-0.688$, $h(4,4,4)=0.56$ and $h(5,5,5)=0.227$.

The simulation results of the estimated cubic kernels, mean and the standard deviations (Std. Dev.) of system I, over 200 Monte Carlo runs with non-Gaussian sample length input $N=2400$ for various SNR are shown in

Table 1. True and estimated (means ± Std. Dev.) cubic kernels of the system I under 200 Monte Carlo runs.

Approach	$\widehat{\mathit{h}}(\mathit{i},\mathit{i},\mathit{i})$ ± $\mathbf{Std}$	0 dB	8 dB	16 dB	24 dB
AlgCum5	$\widehat{h}(1,1,1)$ ± $Std$	$-0.1237$ ± $0.0554$	$-0.1246$ ± $0.0365$	$-0.1278$ ± $0.0347$	$-0.1297$ ± $0.0318$
	$\widehat{h}(2,2,2)$ ± $Std$	$-0.3442$ ± $0.0573$	$-0.3446$ ± $0.0293$	$-0.3464$ ± $0.0205$	$-0.3456$ ± $0.0198$
	$\widehat{h}(3,3,3)$ ± $Std$	$0.8905$ ± $0.0372$	$0.8964$ ± $0.0290$	$0.8982$ ± $0.0275$	$0.8980$ ± $0.0265$
	$\mathrm{NMSE}$	$4.1833\times {10}^{-4}$	$3.5208\times {10}^{-4}$	$2.6065\times {10}^{-4}$	$2.2156\times {10}^{-4}$
CIMLLAD [56]	$\widehat{h}(1,1,1)$ ± $Std$	$-0.1843$ ± $0.2376$	$-0.1353$ ± $0.1882$	$-0.1704$ ± $0.2087$	$-0.1300$ ± $0.1907$
	$\widehat{h}(2,2,2)$ ± $Std$	$-0.3880$ ± $0.2272$	$-0.3457$ ± $0.1874$	$-0.3788$ ± $0.2206$	$-0.3611$ ± $0.1974$
	$\widehat{h}(3,3,3)$ ± $Std$	$0.8836$ ± $0.2223$	$0.8629$ ± $0.1762$	$0.8950$ ± $0.2161$	$0.8773$ ± $0.1903$
	$\mathrm{NMSE}$	$15\times {10}^{-4}$	$8.2349\times {10}^{-4}$	$6.5195\times {10}^{-4}$	$5.2984\times {10}^{-4}$

. To analyze the convergence of this proposed blind approach to identify the nonlinear system I with high levels of noise (i.e., SNR = 0 dB), we represent the fluctuations around the means of the estimated nonlinear cubic kernels in Figure 2, Figure 3 and Figure 4, respectively. To illustrate this convergence and to test and validate it, these fluctuations are compared with those obtained using the adaptive algorithm such as CIM-LLAD ($\mu =0.01$, $\rho =0.0001$, $\tau =1.2$, $\sigma =0.02$, $M=1024$) developed in [56]. The proposed blind method has lower fluctuations around the true values of the identifying cubic kernels when compared with the CIM-LLAD algorithm of [56].

Table 1 compares the different techniques. Clearly, the standard deviations obtained with the proposed approach are smaller than for the CIM-LLAD algorithm proposed in [56]. Moreover, the proposed approach uses only the output and can identify the cubic system blindly without any information about the input model with a good precision. The main reason is that the proposed method uses the fifth-order cumulants which exploit the maximum information by calculating the cumulants of the output making it more robust when we use a small sample size of non-Gaussian signal input. Additionally, Gaussian noise had a minor effect on the estimated cubic kernels when using the proposed approach as indicated in Table 1.

The advantages of the proposed method are even more clear when we look to the NMSE, especially for high noise levels (i.e., SNR = 0 dB). This is well supported by the existing literature as indicated above. Indeed, the additive Gaussian noise will almost vanish in the higher cumulants domain (superior to second-order) and the proposed approach, which uses fifth-order cumulants can achieve excellent identification kernels in very noisy environments compared to the CIM-LLAD algorithm. To illustrate the differences in the performance between the proposed and existing methods, we calculate the corresponding NMSE for various SNR and input signal sample sizes. The corresponding identification results are depicted in Figure 5.

Clearly, the NMSE and standard deviation values obtained using the developed approach are much lower than those for the CIM-LLAD algorithm, for all considered values of SNR and input data sizes, being able to improve significantly the identification of cubic kernels of the system blindly, even in very noisy environments and/or small sample sizes (e.g., SNR = 0 dB, $N=400$).

To study the convergence of the proposed fifth-order cumulants-based approach more clearly, an additional set of simulation results is represented in Figure 6 and Figure 7. It describes through system I the behavior distribution of NMSE and indicates the effect of the Gaussian noise on the estimated cubic kernels in noisy environments for various SNR levels over 200 Monte-Carlo runs.

Figure 8 concerns the differences between the true and estimated output signal using the proposed blind approach (AlgCum5) and the CIM-LLAD method. We take MSE on the true and estimated outputs as criterion of comparison, i.e.,

$$\mathrm{MSE}=\frac{1}{N}\sum _{n=0}^N{\left[y\left(n\right)-\widehat{y}\left(n\right)\right]}^2,$$

(37)

where $y\left(n\right)$ is the true output signal of the cubic system considered and $\widehat{y}\left(n\right)$ is the estimated output signal of the cubic system to be identified.

We conclude from Figure 8 and the results summarized in

Table 2. MSE of different approaches.

Approach	System I	System II
AlgCum5	$0.0452$	$0.0762$
CIMLLAD [56]	$0.5655$	$0.7218$

using MSE criterion that no noticeable differences between true output signal of the system I and their corresponding estimates using the proposed blind approach. On the other hand, for the CIM-LLAD method, we have an MSE about 12 times higher.

The obtained results of nonlinear system II with non-Gaussian input sample of length $N=2400$ and different SNR levels are given in

Table 3. True and estimated (means ± Std. Dev.) cubic kernels of system II.

SNR	$\widehat{\mathit{h}}(\mathit{i},\mathit{i},\mathit{i})$ ± $\mathbf{Std}$	0 dB	8 dB	16 dB	24 dB
	$\widehat{h}(1,1,1)$ ± $Std$	$-0.2449$ ± $0.0591$	$-0.2562$ ± $0.0396$	$-0.2618$ ± $0.0451$	$-0.2587$ ± $0.0432$
	$\widehat{h}(2,2,2)$ ± $Std$	$0.3782$ ± $0.0446$	$0.3705$ ± $0.0332$	$0.3683$ ± $0.0373$	$0.3692$ ± $0.0372$
AlgCum5	$\widehat{h}(3,3,3)$ ± $Std$	$-0.6842$ ± $0.0556$	$-0.6862$ ± $0.0314$	$-0.6841$ ± $0.0328$	$-0.6890$ ± $0.0249$
	$\widehat{h}(4,4,4)$ ± $Std$	$0.5579$ ± $0.0322$	$0.5621$ ± $0.0200$	$0.5613$ ± $0.0283$	$0.5647$ ± $0.0248$
	$\widehat{h}(5,5,5)$ ± $Std$	$0.2407$ ± $0.0853$	$0.2373$ ± $0.0550$	$0.2308$ ± $0.0708$	$0.2334$ ± $0.0649$
	$\mathrm{NMSE}$	$2.7712\times {10}^{-4}$	$1.2968\times {10}^{-4}$	$1.1822\times {10}^{-4}$	$1.1068\times {10}^{-4}$
	$\widehat{h}(1,1,1)$ ± $Std$	$-0.2870$ ± $0.2373$	$-0.2790$ ± $0.2128$	$-0.2349$ ± $0.1949$	$-0.2611$ ± $0.2103$
	$\widehat{h}(2,2,2)$ ± $Std$	$0.3066$ ± $0.2354$	$0.3429$ ± $0.2117$	$0.3482$ ± $0.2072$	$0.3618$ ± $0.2099$
CIMLLAD [56]	$\widehat{h}(3,3,3)$ ± $Std$	$-0.6906$ ± $0.2161$	$-0.6923$ ± $0.2234$	$-0.6834$ ± $0.2013$	$-0.7022$ ± $0.2266$
	$\widehat{h}(4,4,4)$ ± $Std$	$0.5775$ ± $0.2157$	$0.5849$ ± $0.2041$	$0.5324$ ± $0.1990$	$0.5784$ ± $0.2219$
	$\widehat{h}(5,5,5)$ ± $Std$	$0.2104$ ± $0.2398$	$0.2176$ ± $0.2221$	$0.2273$ ± $0.1969$	$0.2444$ ± $0.2125$
	$\mathrm{NMSE}$	$24\times {10}^{-4}$	$9.1512\times {10}^{-4}$	$5.6523\times {10}^{-4}$	$4.7641\times {10}^{-4}$

. Clearly, we can make remarks similar to those made in the case of system I. In fact, there is still a significant improvement when compared to the CIM-LLAD method, regardless of the SNR.

Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 illustrate the convergence of the proposed blind approach in a very noise environment (SNR = 0 dB). From these figures and the Std. Dev. depicted in Table 3, it is clear that our approach converges to the true values of cubic kernels with little fluctuation, even if we use a small sample length input such as N = 1000. The convergence corresponding to CIM-LLAD algorithm has more fluctuations of the estimated kernels when the SNR varies from 0 to 24 dB values as shown in Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13, respectively, and Std. Dev. calculated in Table 3 which is still in the same level for various SNR. For example, the Std. Dev. of the first estimated kernels, $h(1,1,1)$, obtained using the proposed method is about 4 times lower than for CIM-LLAD, once again, with the additional advantage that we can estimate the kernels of the cubic system blindly (i.e., without any information about the excitation non-Gaussian signal), which is not the case of the adaptive algorithm that requires the input and output to identify the kernels of the model.

Figure 14, where we present the NMSE as function the SNR, gives us a good idea about the precision of the estimated kernels. For high noise level (SNR = 0 dB), the NMSE using the proposed approach is almost 9 times lower than for the CIM-LLAD. The reason is the same as in system I: the proposed approach uses only the fifth-order cumulants, which are able to remove most of the Gaussian noise, allowing a more accurate identification of the cubic kernels.

Figure 15 and Figure 16 concern the convergence of the proposed blind approach and the CIM-LLAD method, respectively. Figure 17 shows the differences between the true and estimated output signal using the proposed blind approach (AlgCum5) and the CIM-LLAD method. First, we identify the cubic kernels, then, we convolve it with the non-Gaussian input signal. We also computed the MSE, which is depicted in Table 2, once again, showing the clear advantages of the proposed blind approach.

7. Conclusions

The identification kernels of cubic systems using a new blind fifth-order-based approach was carried out in this paper. The proposed method allows an accurate identification approach and yields the desired kernels via fifth-order cumulants which exploit the maximum information of the output system, allowing estimating it blindly and being able to to remove the Gaussian noise almost completely. The proposed technique was compared with existing ones, showing excellent convergence capabilities and the advantage and accurate estimation of system parameters, even for low SNR and/or small input data size. In fact, it was shown that the proposed schemes achieve significantly better convergence with little fluctuations around the true values of the identifying cubic kernels, even in a very noisy environment, i.e., SNR = 0 dB, and input data size N = 800 compared to the adaptive method which presents more fluctuations.