Robust Hammerstein Adaptive Filtering under Maximum Correntropy Criterion

Abstract

The maximum correntropy criterion (MCC) has recently been successfully applied to adaptive filtering. Adaptive algorithms under MCC show strong robustness against large outliers. In this work, we apply the MCC criterion to develop a robust Hammerstein adaptive filter. Compared with the traditional Hammerstein adaptive filters, which are usually derived based on the well-known mean square error (MSE) criterion, the proposed algorithm can achieve better convergence performance especially in the presence of impulsive non-Gaussian (e.g., α-stable) noises. Additionally, some theoretical results concerning the convergence behavior are also obtained. Simulation examples are presented to confirm the superior performance of the new algorithm.

1. Introduction

Nonlinear system identification is still an active research area [1]. Although linear systems have established a solid theory [2], most practical systems (e.g., hands-free telephone systems) may be more adequately represented as a nonlinear model. One of the main challenges for nonlinear system identification is the choice of an appropriate nonlinear filtering structure that accurately captures the characteristics of the underlying nonlinear system. A common structure used in nonlinear modeling is the block-oriented representation. The Wiener model and the Hammerstein model are two typical block-oriented nonlinear models [3]. Specifically, the Wiener model consists of a cascade of a linear time invariant (LTI) filter followed by a static nonlinear function, indicated as a linear-nonlinear (LN) model [4,5,6], and the Hammerstein model consists of a cascade of a static nonlinear function follow by a LTI filter, known as a nonlinear-linear (NL) model [7,8,9,10,11,12,13,14,15,16,17,18,19]. Other nonlinear models include neural networks (NNs) [20], Volterra adaptive filters (VAFs) [21], kernel adaptive filters (KAF) [22,23,24,25], among others.

Hammerstein filters can accurately model many real-world systems and, as a consequence, they have been successfully used in various applications of engineering [26,27,28,29]. Due to its simplicity and efficiency, the mean square error (MSE) criterion has been widely applied in Hammerstein adaptive filtering [30]. Adaptive algorithms under MSE usually perform very well when the desired signals are disturbed by Gaussian noises. However, when the desired signals are disturbed by non-Gaussian noises, especially in the presence of large outliers (observations that significantly deviate from the bulk of data), the performance of the MSE based algorithms may deteriorate rapidly. Actually, MSE is rather sensitive to outliers. In most practical situations, heavy-tailed impulsive noises may occur, which often cause large outliers. For instance, different types of artificial noises in electronic devices, atmospheric noises, and lighting spikes in natural phenomena, can be described as an impulsive noise [31,32].

In this work, instead of using the MSE criterion, we apply the maximum correntropy criterion (MCC) to develop a robust Hammerstein adaptive filtering algorithm. Correntropy is a nonlinear similarity measure between two signals [33,34]. The MCC aims at maximizing the similarity (measured by correntropy) between the model output and the desired response such that the adaptive model is as close as possible to the unknown system. It has been shown that, the MCC in terms of the stability and accuracy, is very robust with respect to impulsive noises [33,34,35,36,37,38,39]. Compared with the traditional Hammerstein adaptive filtering algorithms based on the MSE criterion, the new algorithm can achieve better performance especially in the presence of impulsive non-Gaussian noises.

The organization of the rest of the paper is as follows. In Section 2, after briefly introducing the correntropy, we derive a Hammerstein adaptive filtering algorithm under MCC criterion. In Section 3, we carry out the convergence analysis. In Section 4, we present simulation examples to demonstrate the superior performance of the proposed algorithm. Finally, we give the conclusion in Section 5.

2. Hammerstein Adaptive Filtering under the Maximum Correntropy Criterion

Figure 1 shows the structure of a Hammerstein adaptive filter under MCC criterion, where the filter consists of a polynomial memoryless nonlinearity followed by a linear FIR filter. This structure has been commonly used in Hammerstein adaptive filtering [8,9,27]. As shown in Figure 1, under the MCC criterion, the parameters of the linear and nonlinear parts are adjusted to maximize the correntropy between the model output and desired response.

Figure 1. Structure of a Hammerstein adaptive filter under maximum correntropy criterion (MCC) criterion.

Figure 1. Structure of a Hammerstein adaptive filter under maximum correntropy criterion (MCC) criterion.

2.1. Correntropy

Correntropy is a nonlinear similarity measure between two signals. Given two random variables X and Y, the correntropy is [33,34,35,36,37,38,39]

$$V(X,Y)=E[\mathrm{\kappa }(X,Y)]={\displaystyle \int \mathrm{\kappa }(x,y)f_{XY}(x,y)}dxdy$$

(1)

where E[·] denotes the expectation operator, κ(·,·) is a shift-invariant Mercer kernel, and f_XY(x, y) stands for the probability density function (PDF) of (X, Y). The most widely used kernel in correntropy is the Gaussian kernel, given by

$${\mathrm{\kappa }}_{\sigma }(x,y)=\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{e^2}{2{\sigma }^2})$$

(2)

where e = x − y, and σ stands for the kernel bandwidth. In this work, without being mentioned otherwise, the kernel function is a Gaussian kernel. In practical situations, the join distribution of X and Y is usually unknown and only a finite number of data {(d(i), y(i))}${}_{i=1}^K$ are available. In these cases, one can use a sample mean estimator of the correntropy:

$${\widehat{V}}_{N,\sigma }(X,Y)=\frac{1}{K}{\displaystyle \sum _{i=1}^K{\mathrm{\kappa }}_{\sigma }(d(i)-y(i))}$$

(3)

The optimization cost under MCC is thus

$$\max J_{MCC}={\displaystyle \sum _{i=1}^K{\mathrm{\kappa }}_{\sigma }(e(i))}$$

(4)

where e(i) = d(i) − y(i). We can evaluate the sensitivity (derivative) of the MCC cost J_MCC with respect to the error e(i),

$$\frac{\partial J_{MCC}}{\partial e(i)}=\frac{-1}{\sqrt{2\pi }{\sigma }^3}\cdot \exp \left(\frac{-e^2(i)}{2{\sigma }^2}\right)e(i)$$

(5)

The derivative curves of −J_MCC for different kernel widths are illustrated in Figure 2. As one can see, when the magnitude of error is very large, the derivative will become rather small especially for a smaller kernel width. Therefore, the MCC training is insensitive (hence robust) to a large error.

Figure 2. Derivative curves of −J_MCC with respect to e(i) for different kernel widths.

Figure 2. Derivative curves of −J_MCC with respect to e(i) for different kernel widths.

2.2. Hammerstein Adaptive Filtering

Assuming that the input-output mapping of the memoryless polynomial nonlinearity is

$$s(n)=p_{1}x(n)+p_2{x}^2(n)+\cdots +p_M{x}^M(n)$$

(6)

where M and p_M denote the polynomial order and the m-th order coefficient, Expression (6) can be rewritten as

$$s(n)={\text{p}}^T(n){\text{x}}_p(n)$$

(7)

where x_p(n) = [x(n) x²(n)···x^M(n)]^T is the polynomial regressor, and p(n) = [p₁ p₂···p_M]^T is the polynomial coefficient vector. The output of the FIR filter can be expressed as

$$y(n)={\text{w}}^T(n)\text{s}(n)$$

(8)

where w(n) = [w₀ w₁···w_N−1]^T is the FIR weight vector, and s(n) = [s(n) s(n − 1)···s(n − N + 1)]^T is the FIR input vector, with N being the FIR memory size. Let X(n) = [x_p(n) x_p(n − 1)···x_p(n − N + 1)]^T. Then we have

$$\text{s}(n)={\mathrm{{\rm X}}}^T(n)\text{p}(n)$$

(9)

Combining Equations (8) and (9) yields

$$y(n)={\text{w}}^T(n)\text{s}(n)={\text{w}}^T(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)$$

(10)

Assume that the unknown system that needs to be identified is also a Hammerstein system with parameter vectors ${\mathbf{p}}^{*}={[p_1^{*} p_2^{*}···p_M^{*}]}^T$ and ${\mathbf{w}}^{*}={[w_0^{*} w_1^{*}···w_{N-1}^{*}]}^T$. Then, the desired signal can be expressed as

$$d(n)={\text{w}}^{*T}{\mathrm{{\rm X}}}^T(n){\text{p}}^{\ast }+v(n)$$

(11)

where v(n) stands for an additive disturbance noise. The error signal can then be calculated as e(n) = d(n) − w^T(n)X^T(n)p(n). In the following, we derive an adaptive algorithm to estimate the Hammerstein parameter vectors using MCC instead of MSE as an optimization criterion. Let us consider the following cost function

$$\begin{array}{l}{J}_{MCC}(\mathbf{\text{p,w}})={\displaystyle \sum _{j=n-L+1}^n{\mathrm{\kappa }}_{\sigma }(d(j),y(j))}\\ =\frac{1}{\sqrt{2\pi }\sigma }{\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}\end{array}$$

(12)

where e(j) = d(j) − y(j), and L denotes the sliding data length. Then, a steepest ascent algorithm for estimating the polynomial coefficient vector can be derived as follows:

$$\begin{array}{l}\frac{\partial J_{MCC}(\mathbf{\text{p,w}})}{\partial \mathbf{\text{p}}(n)}=\frac{1}{\sqrt{2\pi }{\sigma }^3}\cdot {\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}e(j)\cdot \frac{\partial e(j)}{\partial \mathbf{\text{p}}(j)}\\ =\frac{1}{\sqrt{2\pi }{\sigma }^3}\cdot {\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}e(j)\cdot \mathrm{{\rm X}}(j)\mathbf{\text{w}}(j)\end{array}$$

(13)

$$\mathbf{\text{p}}(n+1)=\mathbf{\text{p}}(n)+{\mathrm{\mu }}_p\frac{\partial J_{MCC}(\mathbf{\text{p,w}})}{\partial \mathbf{\text{p}}(n)}=\mathbf{\text{p}}(n)+\frac{{\mathrm{\mu }}_p}{\sqrt{2\pi }{\sigma }^3}\cdot {\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}e(j)\cdot \mathrm{{\rm X}}(j)\mathbf{\text{w}}(j)$$

(14)

In a similar way, we propose the following weight update equation for the coefficients of the FIR filter:

$$\begin{array}{l}\frac{\partial J_{MCC}(\mathbf{\text{p,w}})}{\partial \mathbf{\text{w}}(n)}=\frac{1}{\sqrt{2\pi }{\sigma }^3}\cdot {\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}e(j)\cdot \frac{\partial e(j)}{\partial \mathbf{\text{w}}(j)}\\ =\frac{1}{\sqrt{2\pi }{\sigma }^3}\cdot {\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}e(j)\cdot {\mathrm{{\rm X}}}^T(j)\mathbf{\text{p}}(j)\end{array}$$

(15)

$$\mathbf{\text{w}}(n+1)=\mathbf{\text{w}}(n)+{\mathrm{\mu }}_w\frac{\partial J_{MCC}(\mathbf{\text{p,w}})}{\partial \mathbf{\text{w}}(n)}=\mathbf{\text{w}}(n)+\frac{{\mathrm{\mu }}_w}{\sqrt{2\pi }{\sigma }^3}\cdot {\displaystyle \sum _{j=n-L+1}^n\exp \left(\frac{-e^2(j)}{2{\sigma }^2}\right)}e(j)\cdot {\mathrm{{\rm X}}}^T(j)\mathbf{\text{p}}(j)$$

(16)

In Equations (14) and (16), μ_p and μ_w are, respectively, step-sizes for polynomial nonlinearity subsystem and FIR subsystem. In this work, for simplicity we consider only the stochastic gradient based algorithm (i.e., L = 1). In this case, we have

$$\mathbf{\text{p}}(n+1)=\mathbf{\text{p}}(n)+{\mathrm{\eta }}_p\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)e(n)\cdot \mathrm{{\rm X}}(n)\mathbf{\text{w}}(n)$$

(17)

$$\mathbf{\text{w}}(n+1)=\mathbf{\text{w}}(n)+{\mathrm{\eta }}_w\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)e(n)\cdot {\mathrm{{\rm X}}}^T(n)\mathbf{\text{p}}(n)$$

(18)

where ${\mathrm{\eta }}_p=\frac{{\mathrm{\mu }}_p}{\sqrt{2\pi }{\sigma }^3}$, and ${\mathrm{\eta }}_w=\frac{{\mathrm{\mu }}_w}{\sqrt{2\pi }{\sigma }^3}$. The above update equations are referred to as the Hammerstein adaptive filtering algorithm under MCC criterion, whose pseudocodes are presented in Algorithm 1. The proposed algorithm is in form similar to the traditional Hammerstein adaptive filters under MSE criterion [7], but the step-sizes are different.

Algorithm 1: Hammerstein adaptive filtering Algorithm under MCC.

Parameters setting: μp, μw, σ Initialization: p(0), w(0)

For n = 1, 2, … do (1) $\mathbf{\text{s}}(n)={\mathrm{{\rm X}}}^T(n)\mathbf{\text{p}}(n)$ (2) $y(n)={\mathbf{\text{w}}}^T(n)\mathbf{\text{s}}(n)$ (3) $e(n)=d(n)-y(n)$ (4) $\text{p}(n+1)=\text{p}(n)+{\eta }_p\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)e(n)\cdot \mathrm{{\rm X}}(n)\text{w}(n)$ (5) $\text{w}(n+1)=\text{w}(n)+{\eta }_w\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)e(n)\cdot {\mathrm{{\rm X}}}^T(n)\text{p}(n)$ End for

3. Convergence Analysis

3.1. Stability Analysis

Using the Taylor series expansion of the error e(n + 1) around the instant n and keeping only the linear term, we have [4,7,40]

$$e(n+1)=e(n)+\frac{\partial e(n)}{\partial \text{p}(n)}{|}_{\text{w}(n)=const}\mathrm{\Delta }\text{p}(n)+\frac{\partial e(n)}{\partial \text{w}(n)}{|}_{\text{p}(n)=const}\mathrm{\Delta }\text{w}(n)+h.o.t$$

(19)

where h.o.t denotes higher-order terms. Combining Equations (11), (17) and (18), we can obtain

$$\frac{\partial e(n)}{\partial \text{p}(n)}=-\mathrm{{\rm X}}(n)\text{w}(n)$$

(20)

$$\frac{\partial e(n)}{\partial \text{w}(n)}=-{\mathrm{{\rm X}}}^T(n)\text{p}(n)$$

(21)

$$\mathrm{\Delta }\text{p}(n)={\eta }_p\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)e(n)\cdot \mathrm{{\rm X}}(n)\text{w}(n)$$

(22)

$$\mathrm{\Delta }\text{w}(n)={\eta }_w\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)e(n)\cdot {\mathrm{{\rm X}}}^T(n)\text{p}(n)$$

(23)

Substituting Equations (20)–(23) in Equation (19), and after simple manipulation, we have

$$e(n+1)=\left[1-{\eta }_p\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)\cdot {\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2-{\eta }_w\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)\cdot {\Vert {\mathrm{{\rm X}}}^T(n)\text{p}(n)\Vert }^2\right]e(n)$$

(24)

To ensure the stability of the proposed algorithm, we must assure that |e(n + 1)| ≤ |e(n)|, and hence

$$\left|1-{\eta }_p\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)\cdot {\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2-{\eta }_w\cdot \exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)\cdot {\Vert {\mathrm{{\rm X}}}^T(n)\text{p}(n)\Vert }^2\right|\le 1$$

(25)

which yields

$$0<{\eta }_p\cdot {\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2+{\eta }_w\cdot {\Vert {\mathrm{{\rm X}}}^T(n)\text{p}(n)\Vert }^2\le \frac{2}{\exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)}$$

(26)

Since $\exp \left(\frac{-e^2(n)}{2{\sigma }^2}\right)\le 1$, the following condition guarantees convergence:

$$0<{\eta }_p\cdot {\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2+{\eta }_w\cdot {\Vert {\mathrm{{\rm X}}}^T(n)\text{p}(n)\Vert }^2\le 2$$

(27)

Remark 1. The derived bound on step-sizes is only of theoretical importance as in general, Equation (27) cannot be verified in a practical situation. Similar theoretical results can be found in [7].

3.2. Steady-State Mean Square Performance

We denote e_pw(n) the a priori error of the whole system, e_p(n) the a priori error when only the nonlinear part is adapted while the linear filter is fixed, and e_w(n) the a priori error when only the linear filter is adapted while the nonlinear part is fixed. Let ${\mathrm{{\rm H}}}_p=\underset{n\to \infty }{\lim }E\left[e_p^2(n)\right]$, ${\mathrm{{\rm H}}}_w=\underset{n\to \infty }{\lim }E\left[e_w^2(n)\right]$, and ${\mathrm{{\rm H}}}_{pw}=\underset{n\to \infty }{\lim }E\left[e_{pw}^2(n)\right]$ be the steady-state excess mean square errors (EMSEs). In addition, we denote

$$f(e(i))=\exp \left(\frac{-e^2(i)}{2{\sigma }^2}\right)e(i)$$

(28)

Before evaluating the theoretical values of the steady-state EMSEs, we make the following assumptions:

(A)

The noise v(n) is zero-mean, independent, identically distributed, and is independent of the input X(n), $\widehat{\mathbf{\text{s}}}(n)$ and e(n).

(B)

The a priori errors e_p(n) and e_w(n) are zero-mean Gaussian, and independent of the noise v(n).

(C)

||X(n)w(n)||² and ||X^T(n)p(n)||² are asymptotically uncorrelated with f²(e(n)), that is

$$\underset{n\to \infty }{\lim }E\left[{\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2{f}^2(e(n))\right]=Tr(R_{WX})\underset{n\to \infty }{\lim }E\left[f^2(e(n))\right]$$

(29)

$$\underset{n\to \infty }{\lim }E\left[{\Vert {\mathrm{{\rm X}}}^T(n)\text{p}(n)\Vert }^2{f}^2(e(n))\right]=Tr(R_{XP})\underset{n\to \infty }{\lim }E\left[f^2(e(n))\right]$$

(30)

where $R_{WX}=E\left[(\mathrm{{\rm X}}(n)\text{w}(n)){(\mathrm{{\rm X}}(n)\text{w}(n))}^T\right]$ and $R_{XP}=E\left[({\mathrm{{\rm X}}}^T(n)\text{p}(n)){({\mathrm{{\rm X}}}^T(n)\text{p}(n))}^T\right]$ are the covariance matrices, and $Tr(\cdot )$ denotes the trace operator.

Remark 2. For the assumption (A), it is very common to assume that the noise is independent of the regression vector [41,42,43]. In addition, the noise is often restricted to be zero-mean, identically distributed [33,34,35]. As discussed in [44,45], the assumption (B) is reasonable for long adaptive filters. Since e_p(n) is the a priori error when only the nonlinear part is adapted while the linear filter is fixed, we have the approximation w^* ≈ w(n) such that w(n) is asymptotically uncorrelated with f²(e(n)). Due to the independent assumption (A), X(n) is also asymptotically uncorrelated with f²(e(n)). So ||X^T(n)w(n)||² is asymptotically uncorrelated with f²(e(n)). Similarly, ||X^T(n)p(n)||² is asymptotically uncorrelated with f²(e(n)). Therefore, the assumption (C) is rational.

When only the polynomial part with parameter vector p is adapted, the error e_p(n) is

$$e_p(n)={\text{w}}^{*T}{\mathrm{{\rm X}}}^T(n){\text{p}}^{\ast }-{\text{w}}^T(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)\approx {\text{w}}^T(n){\mathrm{{\rm X}}}^T(n)\tilde{\text{p}}(n)$$

(31)

where $\widetilde{\mathbf{\text{p}}}(n)={\mathbf{\text{p}}}^{*}-{\mathbf{\text{p}}}^T(n)$. In Equation (31), we use the approximation w^* ≈ w(n) at steady-state. From Equation (17), it follows easily that

$$\tilde{\text{p}}(n+1)=\tilde{\text{p}}(n)-{\eta }_p\cdot f(e(n))\cdot \mathrm{{\rm X}}(n)\text{w}(n)$$

(32)

Squaring both sides of Equation (32), we have

$${\Vert \tilde{\text{p}}(n+1)\Vert }^2={\Vert \tilde{\text{p}}(n)\Vert }^2-2{\eta }_p\cdot f(e(n))\cdot e_p(n)+{\eta }_p^2\cdot f^2(e(n))\cdot {\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2$$

(33)

Taking the expectations of the both sides of Equation (33) yields

$$E\left[{\Vert \tilde{\text{p}}(n+1)\Vert }^2\right]=E\left[{\Vert \tilde{\text{p}}(n)\Vert }^2\right]-2{\eta }_p\cdot E\left[f(e(n))\cdot e_p(n)\right]+{\eta }_p^2\cdot E\left[f^2(e(n))\cdot {\Vert \mathrm{{\rm X}}(n)\text{w}(n)\Vert }^2\right]$$

(34)

Assuming the filter is stable and attains the steady state, it holds

$$\underset{n\to \infty }{\lim }E\left[{\Vert \tilde{\text{p}}(n+1)\Vert }^2\right]=\underset{n\to \infty }{\lim }E\left[{\Vert \tilde{\text{p}}(n)\Vert }^2\right]$$

(35)

Combining Equations (34) and (35) and the above assumptions, we obtain

$$2\cdot \underset{n\to \infty }{\lim }E\left[f(e(n))\cdot e_p(n)\right]={\eta }_{p}Tr(R_{WX})\underset{n\to \infty }{\lim }E\left[f^2(e(n))\right]$$

(36)

In order to derive a theoretical value of the steady-state EMSE, we consider two cases below.

Case A. Gaussian Noise

Recalling that e(n) = e_p(n) + v(n), and assuming that the noise v(n) is zero-mean Gaussian, with variance ${\mathrm{ϛ}}_v^2$, we get [34]

$$\underset{n\to \infty }{\lim }E\left[f(e(n))\cdot e_p(n)\right]\text{=}\frac{{\sigma }^3{\mathrm{{\rm H}}}_p}{{({\sigma }^2+{\varsigma }_v^2+{\mathrm{{\rm H}}}_p)}^{3/2}}$$

(37)

where ${\mathrm{\sigma }}_e^2$ denotes the variance of the error, and ${\mathrm{\sigma }}_e^2$ = E[$e_p^2$(n)] + ${\mathrm{ϛ}}_v^2$. Similarly, we obtain [29]

$$\underset{n\to \infty }{\lim }E\left[f^2(e(n))\right]=\frac{{\sigma }^3({\mathrm{{\rm H}}}_p+{\varsigma }_v^2)}{{({\sigma }^2+2{\varsigma }_v^2+2{\mathrm{{\rm H}}}_p)}^{3/2}}$$

(38)

Substituting Equations (37) and (38) into Equation (36), we have

$$2\cdot \frac{{\sigma }^3{\mathrm{{\rm H}}}_p}{{({\sigma }^2+{\varsigma }_v^2+{\mathrm{{\rm H}}}_p)}^{3/2}}={\eta }_{p}Tr(R_{WX})\frac{{\sigma }^3({\mathrm{{\rm H}}}_p+{\varsigma }_v^2)}{{({\sigma }^2+2{\varsigma }_v^2+2{\mathrm{{\rm H}}}_p)}^{3/2}}$$

(39)

Therefore, the steady-state EMSE H_p satisfies

$${\mathrm{{\rm H}}}_p=\frac{{\eta }_p}{2}Tr(R_{WX})\frac{({\mathrm{{\rm H}}}_p+{\varsigma }_v^2){({\sigma }^2+{\varsigma }_v^2+{\mathrm{{\rm H}}}_p)}^{3/2}}{{({\sigma }^2+2{\varsigma }_v^2+2{\mathrm{{\rm H}}}_p)}^{3/2}}$$

(40)

Theorem 1. In a Gaussian noise environment and with the same step-size, the proposed nonlinear Hammerstein adaptive filter under MCC criterion has a smaller steady-state EMSE than under MSE criterion. As the kernel width increases, their values of the steady-state EMSE will become almost identical.

Proof. It can be shown that [34]

$${\mathrm{{\rm H}}}_{p-MSE}=\frac{{\eta }_{p}Tr(R_{WX}){\varsigma }_v^2}{2-{\eta }_{p}Tr(R_{WX})}$$

(41)

where H_p−MSE denotes the steady-state EMSE under MSE criterion. From Equation (40), we have

$${\mathrm{{\rm H}}}_p=\frac{\epsilon {\eta }_{p}Tr(R_{WX}){\varsigma }_v^2}{2-\epsilon {\eta }_{p}Tr(R_{WX})}$$

(42)

where $\epsilon =\frac{{({\sigma }^2+{\varsigma }_v^2+{\mathrm{{\rm H}}}_p)}^{3/2}}{{({\sigma }^2+2{\varsigma }_v^2+2{\mathrm{{\rm H}}}_p)}^{3/2}}$. Since $\epsilon <1$, it holds

$${\mathrm{{\rm H}}}_p<{\mathrm{{\rm H}}}_{p-MSE}$$

(43)

Further, as σ → ∞, we have H_p → H_p−MSE.

Case B. Non-Gaussian Noise

Taking the Taylor series expansion of f(e(n)) around v(n) yields

$$\begin{array}{l}f(e(n))=f(e_p(n)+v(n))\\ =f(v(n))+f^{\prime }(v(n))e_p(n)+\frac{1}{2}{f}^{″}(v(n))e_p^2(n)+h.o.t\end{array}$$

(44)

with

$$\begin{array}{l}{f}^{\prime }(v(n))=\exp \left(\frac{-v^2(n)}{2{\sigma }^2}\right)\left(1-\frac{v^2(n)}{{\sigma }^2}\right)\\ f^{″}(v(n))=\exp \left(\frac{-v^2(n)}{2{\sigma }^2}\right)\left(\frac{v^3(n)}{{\sigma }^4}-\frac{3v(n)}{{\sigma }^2}\right)\end{array}$$

(45)

Under the assumptions (A) and (B), we get [34]

$$\underset{n\to \infty }{\lim }E\left[f(e(n))\cdot e_p(n)\right]\approx E\left[f^{\prime }(v(n))\right]{\mathrm{{\rm H}}}_p$$

(46)

$$\underset{n\to \infty }{\lim }E\left[f^2(e(n))\right]\approx E\left[f^2(v(n))\right]+E\left[f(v(n))f^{″}(v(n))+{\left|f^{\prime }(v(n))\right|}^2\right]{\mathrm{{\rm H}}}_p$$

(47)

Substituting Equations (46) and (47) into Equation (36), we have

$${\mathrm{{\rm H}}}_p=\frac{{\eta }_{p}Tr(R_{WX})E\left[f^2(v(n))\right]}{2E\left[f^{\prime }(v(n))\right]-{\eta }_{p}Tr(R_{WX})E\left[f(v(n))f^{″}(v(n))+{\left|f^{\prime }(v(n))\right|}^2\right]}$$

(48)

Further, substituting Equation (45) into Equation (48), we obtain

$${\mathrm{{\rm H}}}_p=\frac{{\eta }_{p}Tr(R_{WX})E\left[\exp \left(\frac{-v^2(n)}{{\sigma }^2}\right)v^2(n)\right]}{2E\left[\exp \left(\frac{-v^2(n)}{2{\sigma }^2}\right)\left(1-\frac{v^2(n)}{{\sigma }^2}\right)\right]-{\eta }_{p}Tr(R_{WX})E\left[\exp \left(\frac{-v^2(n)}{{\sigma }^2}\right)\left(1+\frac{2v^4(n)}{{\sigma }^4}-\frac{5v^2(n)}{{\sigma }^2}\right)\right]}$$

(49)

When only the linear filter with parameter vector w(n) is adapted, we get

$$2\cdot \underset{n\to \infty }{\lim }E\left[f(e(n))\cdot e_w(n)\right]={\eta }_{p}Tr(R_{XP})\underset{n\to \infty }{\lim }E\left[f^2(e(n))\right]$$

(50)

where $e_w(n)\approx {\tilde{\text{w}}}^T(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)$, $\tilde{\text{w}}(n)={\text{w}}^{*}-\text{w}(n)$. For Gaussian noise case, we obtain

$${\mathrm{{\rm H}}}_w=\frac{{\eta }_w}{2}Tr(R_{XP})\frac{({\mathrm{{\rm H}}}_w+{\varsigma }_v^2){({\sigma }^2+{\varsigma }_v^2+{\mathrm{{\rm H}}}_w)}^{3/2}}{{({\sigma }^2+2{\varsigma }_v^2+2{\mathrm{{\rm H}}}_w)}^{3/2}}$$

(51)

In non-Gaussian environments, we have

$${\mathrm{{\rm H}}}_w=\frac{{\eta }_{w}Tr(R_{XP})E\left[\exp \left(\frac{-v^2(n)}{{\sigma }^2}\right)v^2(n)\right]}{2E\left[\exp \left(\frac{-v^2(n)}{2{\sigma }^2}\right)\left(1-\frac{v^2(n)}{{\sigma }^2}\right)\right]-{\eta }_{w}Tr(R_{XP})E\left[\exp \left(\frac{-v^2(n)}{{\sigma }^2}\right)\left(1+\frac{2v^4(n)}{{\sigma }^4}-\frac{5v^2(n)}{{\sigma }^2}\right)\right]}$$

(52)

Theorem 2. The H_pw satisfies the following condition

$${\mathrm{{\rm H}}}_{pw}\ge {\mathrm{{\rm H}}}_p+{\mathrm{{\rm H}}}_w$$

(53)

Proof. Using Equation (31), we derive

$$\begin{array}{l}{e}_{pw}(n)={\text{w}}^{*T}{\mathrm{{\rm X}}}^T(n){\text{p}}^{\ast }-{\text{w}}^T(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)\\ =\left({\text{w}}^{*T}{\mathrm{{\rm X}}}^T(n){\text{p}}^{\ast }-{\text{w}}^{*T}(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)\right)\text{+}\left({\text{w}}^{*T}{\mathrm{{\rm X}}}^T(n)\text{p}(n)-{\text{w}}^T(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)\right)\\ \approx {\text{w}}^T(n){\mathrm{{\rm X}}}^T(n)\tilde{\text{p}}(n){\tilde{\text{w}}}^T(n){\mathrm{{\rm X}}}^T(n)\text{p}(n)\\ \text{ =}{e}_p(n)\text{+}{e}_w(n)\end{array}$$

(54)

It follows that

$$\begin{array}{l}{\mathrm{{\rm H}}}_{pw}=\underset{n\to \infty }{\lim }E\left[{e^2}_{pw}(n)\right]\approx \underset{n\to \infty }{\lim }E\left[{\left(e_p(n)\text{+}{e}_w(n)\right)}^2\right]\\ \text{ =}\underset{n\to \infty }{\lim }E\left[e_p^2(n)\right]+2\underset{n\to \infty }{\lim }E\left[e_p(n)e_w(n)\right]+\underset{n\to \infty }{\lim }E\left[e_w^2(n)\right]\\ \text{ =}{\mathrm{{\rm H}}}_p+{\mathrm{{\rm H}}}_w+2{\mathrm{{\rm H}}}_{cross}\end{array}$$

(55)

where ${\mathrm{{\rm H}}}_{cross}=\underset{n\to \infty }{\lim }E\left[e_p(n)e_w(n)\right]$ stands for the cross-EMSE and H_cross ≥ 0 (H_cross = 0 when e_p(n) and e_w(n) are statistically independent and zero mean) [7]. Therefore, H_pw ≥ H_p + H_w, which completes the proof.

4. Simulation Results

Now, we present simulation results to demonstrate the performance of the Hammerstein adaptive filtering under MCC. In order to show the performance of the proposed algorithm in non-Gaussian noises, we adopt the alpha-stable distribution to generate the disturbance noise, whose characteristic function is [32,46]

$$f(\mathrm{t})=\exp \{\mathrm{j}\delta t-\gamma |\mathrm{t}{|}^{\alpha }[1+\mathrm{j}\beta \mathrm{sgn}(\mathrm{t})\mathrm{S}(\mathrm{t},\alpha )]\}$$

(56)

in which

$$\mathrm{S}(\mathrm{t},\alpha )=\{\begin{array}{cc}\tan \frac{\alpha \pi }{2}& if \alpha \ne 1\\ \frac{2}{\pi }\log \left|t\right|& if \alpha =1\end{array}$$

(57)

where α ϵ (0, 2] denotes the characteristic factor, −∞ < δ < +∞ is the location parameter, β ϵ [−1, 1] stands for the symmetry parameter, and γ > 0 is the dispersion parameter. The characteristic factor α measures the tail heaviness of the distribution. The smaller α is, the heavier the tail is. In addition, γ measures the dispersion of the distribution. The distribution is symmetric about its location δ when β = 0. Such a distribution is called a symmetric alpha-stable (SαS) distribution. The parameters vector of the noise model is defined as V = (α, β, γ, δ).

In the simulations below, the input signal considered is a colored signal obtained from the following equation:

$$x(n)=ax(n-1)+\sqrt{1-a^2}\xi (n)$$

(58)

with a = 0.95, and ξ(n) being a white Gaussian signal of unit variance. In addition, the coefficient vectors are initialized with the first coefficient equal to 1 and the others equal to zero [7].

4.1. Experiment 1

First, we consider an unknown Hammerstein system with parameter vectors p^* = [1, 0.6], w^* = [1, 0.6, 0.1, −0.2, −0.06, 0.04, 0.02, −0.03, −0.02, 0.01]. Thus, M = 2, N = 10. The kernel width σ is 1.0. The noise vector V is set at (1.2, 0, 0.6, 0), and the noise signal is shown in Figure 3. Simulation results are averaged over 100 independent Monte Carlo runs, and in each simulation, 15,000 iterations are run to ensure the algorithm will reach the steady state, and the steady-state MSE is obtained as an average over the last 2000 iterations. The step-sizes are set at μ_p = μ_w = 0.005 and μ_p = 0.01, μ_w = 0.01 for MSE and MCC, respectively. Figure 4 shows the average convergence curves under MCC and MSE. As we can see, the Hammerstein adaptive filtering under MCC criterion achieves faster convergence speed and lower steady-state testing MSE than under MSE criterion. Here the testing MSE is evaluated on a test set with 100 samples.

Figure 3. A typical sequence of the alpha-stable noise with V = (1.2, 0, 0.6, 0).

Figure 3. A typical sequence of the alpha-stable noise with V = (1.2, 0, 0.6, 0).

Figure 4. Convergence curves under maximum correntropy criterion (MCC) and mean square error (MSE) (for unknown system with polynomial nonlinearity).

Figure 4. Convergence curves under maximum correntropy criterion (MCC) and mean square error (MSE) (for unknown system with polynomial nonlinearity).

Second, we investigate the performance of the algorithms with different noise parameters. The steady-state MSEs with different γ (0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6) and different α (0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2.0) are shown in Figure 5 and Figure 6, respectively. We observe: (1) In most cases, the new algorithm performs better and achieves a lower steady-state MSE compared with the Hammerstein adaptive filtering under MSE criterion; (2) When α is close to 2.0, the Hammerstein adaptive filtering under MSE criterion can achieve better performance than under MCC criterion. The main reason for this is that, when α ≈ 2.0, the noise will be approximately Gaussian. Simulation results suggest that the proposed algorithm is particularly useful for identifying a Hammerstein system in non-Gaussian noises.

Figure 5. Steady-state mean square error (MSE) with different γ (α = 1.2).

Figure 5. Steady-state mean square error (MSE) with different γ (α = 1.2).

Figure 6. Steady-state mean square error (MSE) with different α (γ = 0.6).

Figure 6. Steady-state mean square error (MSE) with different α (γ = 0.6).

4.2. Experiment 2

The second experiment is drawn from [47]. The nonlinear dynamic system is composed of two blocks. The first block is a non-polynomial nonlinearity

$$s(n)=\sqrt[3]{x(n)}$$

(59)

while the second block is an FIR filter with weight vector

$$h={\left[10.750.50.250-0.25\right]}^T$$

(60)

The noise vector V is set at (1.0, 0, 0.8, 0) (see Figure 7 for a typical sequence of the noise), and the polynomial order M and the FIR memory size N are set at 3 and 6, respectively. Simulation results are averaged over 50 independent Monte Carlo runs, and in each simulation, 30,000 iterations are run to ensure the algorithm will reach the steady state, and the steady-state MSE is obtained as an average over the last 2000 iterations. The testing MSE is evaluated on a test set with 100 samples. Figure 8 demonstrates the convergence curves under MCC and MSE. For both adaptive filtering algorithms, the step-sizes are set at μ_p = 0.005, μ_w = 0.015. It can be seen that, the Hammerstein adaptive filter under MCC criterion performs better (say, with faster convergence speed and smaller mismatch error) than under MSE criterion.

Figure 7. A typical sequence of the alpha-stable noise with V = (1.0, 0, 0.8, 0).

Figure 7. A typical sequence of the alpha-stable noise with V = (1.0, 0, 0.8, 0).

Figure 8. Convergence curves under maximum correntropy criterion (MCC) and mean square error (MSE) (for unknown system with non-polynomial nonlinearity).

Figure 8. Convergence curves under maximum correntropy criterion (MCC) and mean square error (MSE) (for unknown system with non-polynomial nonlinearity).

Finally, we show the steady-state performance of the algorithms with different kernel widths σ (0.01, 1.0, 2.0, 3.0, 4.0, 5.0). Simulation results are shown in Figure 9. As we can see, the kernel width has significant influence on the performance of the proposed algorithm. In this example, the lowest steady-state MSE is obtained when σ = 1.0.

Figure 9. Steady-state mean square error (MSE) with different kernel widths.

Figure 9. Steady-state mean square error (MSE) with different kernel widths.

5. Conclusions

The MCC has been successfully applied in domains of machine learning and signal processing due to its strong robustness in impulsive non-Gaussian situations. In this work, we develop a robust Hammerstein adaptive filter under MCC criterion. Different from the traditional Hammerstein adaptive filtering algorithms, the new algorithm use the MCC instead of the well-known MSE as the adaptation criterion, which can achieve desirable performance especially in impulsive noises. Based on [7,31], we carry out the convergence analysis, and obtain some important theoretical results. Simulation examples confirm the excellent performance of the proposed algorithm. How to verify the derived theoretical results is an interesting topic for future study.