Nonlinear Filtering of Weak Chaotic Signals

Abstract

In recent years, the application of nonlinear filtering for processing chaotic signals has become relevant. A common factor in all nonlinear filtering algorithms is that they operate in an instantaneous fashion, that is, at each cycle, a one moment of time magnitude of the signal of interest is processed. This operation regime yields good performance metrics, in terms of mean squared error (MSE) when the signal-to-noise ratio (SNR) is greater than one and shows moderate degradation for SNR values no smaller than −3 dB. Many practical applications require detection for smaller SNR values (weak signals). This chapter presents the theoretical tools and developments that allow nonlinear filtering of weak chaotic signals, avoiding the degradation of the MSE when the SNR is rather small. The innovation introduced through this approach is that the nonlinear filtering becomes multimoment, that is, the influence of more than one moment of time magnitudes is involved in the processing. Some other approaches are also presented.

1. Introduction

The detection of chaotic (stochastic) weak signals is relevant (among others) for applications such as biomedical telemetry [1, 2], seismological signal processing [3], underwater signal processing [4], interference modeling [5], etc. Effective detection of weak and rather weak chaotic signals (−3 dB or less) is a challenge whose solution can improve, for example, the link budget (communication distance). Among different approaches to this problem, one can mention techniques such as stochastic resonance [4], instantaneous spectral cloning [6], etc. The problem in this chapter is addressed from the standpoint of nonlinear filtering techniques which earlier was designed to operate with signal-to-noise ratio (SNR) values bigger than one or at least rather close to one (with an acceptable slight degradation as the SNR approaches −3 dB [7]. Far down −3 dB, the performance of the available filtering methods drops down sharply and becomes ineffective. One of the possible explanations for this issue is that current nonlinear filtering algorithms can be considered as one moment in the sense that they operate in an instantaneous fashion, that is, during each operation cycle, they process an instantaneous one moment of time magnitude of the received aggregate signal; in the next cycle, a new instantaneous one moment of time magnitude is processed and so on. This is precisely the operation rule for all known optimum algorithms and their quasi-optimum versions as well, for instance, the extended Kalman filter (EKF) [7], but it can also be found in strategies such as unscented Kalman filter (UKF), Gauss-Hermite filter (GHF), and quadrature Kalman filter (QKF), among others. One of the goals of this chapter is to describe the detection of weak chaotic signals applying the principles of noninstantaneous filtering in a block way, that is, multimoment filtering theory [8], through a real-time implementation in a digital signal processing (DSP) block. Moreover, some space of this chapter will be dedicated to the conditionally optimum approach for the nonlinear filtering methods as well, together with some asymptotic methods.

Theoretically, for many cases, the chaos might be represented as an output signal of dissipative continuous dynamic systems (strange attractors) [9]:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\left(t\right)\right),\mathit{x}\in R^n,\mathit{x}\left(t_0\right)={\mathit{x}}_0,$$E1

where f(•) = [f₁(x), …f_n(x)]^T is a differentiable vector function.

According to the idea of Kolmogorov, the equations for strange attractors (1) can be successfully transformed in the equivalent stochastic form as a stochastic differential equation (SDE) [9, 10]:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\left(t\right)\right)+\mathbf{\epsilon \xi }\left(t\right).$$E2

The influence of a weak external source of white noise is denoted by ξ(t), and the noise intensities are given in a matrix form ε = [ε_ij]^nxn.

Note that a stationary distribution W_st(x) exists even when the weak white noise component is tending to zero [11, 12, 13].

Nonlinear filtering of chaotic desired signals comes up naturally when SDE (2) is used as model of chaos. This follows straight from the classical theory of nonlinear filtering for Markov processes, proposed more than 50 years ago [14, 15] and extensively developed in subsequent studies [16, 17, 18, 19, 20, 21], although those methods are still under development.

From the practical implementation point of view, the nonlinear filtering strategies are approximate (see the references above). This follows from the fact that, in general, there is no analytical solution for the a posteriori probability density functions when one attempts solving the Stratonovich-Kushner equations (SKE).

In the following, some of the numerous nonlinear filtering approximate approaches that have been developed will be presented.

2. Nonlinear filtering for Markovian processes

Let assume that filtering of the following received signal is required:

$$\mathit{y}\left(t\right)=\mathit{s}\left(t,\mathit{x}\left(t\right)\right)+{\mathit{n}}_0\left(t\right),$$E3

where s (⋅) is a vector function of the “message dependent” desired signal (which is subject of filtering) of dimension “m,” the received signal is denoted by the vector y(t) (also of dimension “m”), and n₀ is a vector of the white additive noises characterized by the intensity matrix N₀(m × m). The following SDE is used to model the signal s (⋅) as an n-dimensional Markov diffusion process [22]:

$$\dot{\mathit{x}}=\mathit{g}\left(t,\mathit{x}\right)+\mathbf{\xi }\left(t\right).$$E4

Strictly speaking, Eqs. (4) and (2) are the same SDE, and the vector function g (⋅) substitutes f (⋅) in (2); for (4), D denotes the correspondent matrix of intensities for ξ(⋅).

Under this assumption ([14, 22] and so on), one can use the so-called Fokker-Planck-Kolmogorov (FPK) equation in order to solve the a priori probability density function (a priori PDF), for x(t):

$$\frac{\partial W_{\mathit{PR}}\left(\mathit{x},t\right)}{\partial t}=-\sum _{i=1}^n\frac{\partial }{\partial x_i}\left[g_i\left(t,\mathit{x}\right)W_{\mathit{PR}}(\mathit{x},t)\right]+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n\frac{{\partial }^2}{\partial x_i\partial x_j}\left[D_{\mathit{ij}}{W}_{\mathit{PR}}\left(\mathit{x},t\right)\right],$$E5

where W_PR(x,t₀) = W₀(x).

The Eq. (5) can be rewritten in another form [21, 23] as well:

$$\frac{\partial W_{\mathit{PR}}\left(\mathit{x},t\right)}{\partial t}=-\mathit{div}\mathbf{\pi }\left(\mathit{x},t\right),$$E6

or

$$\frac{\partial W_{\mathit{PR}}\left(\mathit{x},t\right)}{\partial t}=L_{\mathit{PR}}\left\{W_{\mathit{PR}}\left(\mathit{x},t\right)\right\},$$E7

where π(x, t) is a probabilistic “flow” with components:

$$\mathbf{\pi }\left(\mathit{x},t\right)=g_i\left(\mathit{x},t\right)W_{\mathit{PR}}\left(\mathit{x},t\right)-\frac{1}{2}\sum _{j=1}^n\frac{\partial }{\partial x_j}\left[D_{\mathit{ij}}{W}_{\mathit{PR}}\left(\mathit{x},t\right)\right].$$E8

In Eqs. (5)–(8), {D_ij} denote diffusion coefficients of the Markov process and ${\left\{g_i\left(\mathit{x},t\right)\right\}}_1^n$ are the correspondent drift coefficients, and both of them will be used in the Stratonovich sense [14, 22]; L_PR{⋅} denotes a FPK linear operator.

The integrodifferential equation for the a posteriori probability density function W_PS(x, t) is given by any of the two equivalent expressions (see [14]:

$$\frac{\partial W_{\mathit{PS}}\left(\mathit{x},t\right)}{\partial t}=L_{\mathit{PR}}\left\{W_{\mathit{PS}}\left(\mathit{x},t\right)\right\}+\frac{1}{2}\left[\mathit{F}\left(\mathit{x},t\right)-\underset{-\infty }{\overset{\infty }{\int }}\mathit{F}(\mathit{x},t\left)W_{\mathit{PS}}\right(\mathit{x},t)d\mathit{x}\right]W_{\mathit{PS}}\left(\mathit{x},t\right)$$E9

or

$$\frac{\partial W_{\mathit{PS}}\left(\mathit{x},t\right)}{\partial t}=-\mathit{div}\widehat{\mathbf{\pi }}\left(\mathit{x},t\right)+\frac{1}{2}\left[\mathit{F}\left(\mathit{x},t\right)-<F(\mathit{x},t)>\right]W_{\mathit{PS}}\left(\mathit{x},t\right),$$E10

where $<F\left(\mathit{x},t\right)>$ denotes the averaging of F(x,t) given by $<F\left(\mathit{x},t\right)>=\underset{-\infty }{\overset{\infty }{\int }}\mathit{F}\left(\mathit{x},t\right)W_{\mathit{PS}}\left(\mathit{x},t\right)d\mathit{x}$, $\widehat{\pi }\left(\mathit{x},t\right)$ is (5), W_PR(x, t) is substituted by W_PS(x, t), and

$$\mathit{F}\left(\mathit{x},t\right)={\left[\mathit{y}\left(t\right)-\frac{1}{2}s(\mathit{x},t)\right]}^T{\mathit{N}}_0^{-1}\left[\mathit{y}\left(t\right)-\frac{1}{2}s(\mathit{x},t)\right].$$E11

The combination of Eqs. (9)–(11) is known as the Stratonovich-Kushner nonlinear equations (SKE), and they have an appealing physical sense: the first term in (9) represents the dynamics of the a priori data of x(t). For the second term, the analysis of observations is used to drive the innovation of the a priori data.

Using any optimization criteria, one can get $\widehat{\mathit{x}}\left(t\right)$ (the optimum estimation of x(t)) which comes as a solution of (9), when y(t) is the input signal, that is, filtering of x(t).

Here, one has to note that Eq. (9) turns into FPK (6) if the intensity of the additive noises N₀ grows (the first term in (9) is dominant), and as a consequence, the filtering accuracy diminishes drastically. In the opposed scenario (large signal-to-noise ratio), the W_PS(x, t) tends to the unimodal Gaussian PDF [14, 20].

Note that the time evolution of W_PS(x, t) is completely described by the SKE but, as it was mentioned earlier, does not provide exact analytical solutions. There are very few exceptions: linear SDE (4) which yields the well-known Kalman filtering algorithm [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], the Zakai approach [25], and so on. Due to this, the nonlinear filtering algorithms are practically always approximate. As it was mentioned before, during almost 50 years of intensive research, the bibliography for nonlinear filtering algorithms has become enormous; in the next section, we will consider only few of those works taking into account the following considerations:

• the models applied for filtering of chaos correspond to the equations for Rössler, Chua, and Lorenz strange attractors with n = 3, that is, low dimensional;

• the algorithms for nonlinear filtering have to be of reduced computational complexity in order to satisfy real-time application requirements;

• the algorithms for nonlinear filtering, according to the aim of the material of the chapter, have to be able to perform satisfactorily in scenarios with low or very low signal-to-noise ratios (SNR), although the Gaussian assumption for W_PS(x) is not always valid;

• $$s\left(\mathit{x}\left(t\right),t\right)\cong \mathit{x}\left(t\right);$$E12

• All D_ij are equal to zero, except D₁₁ ≅ D₁ [11].

2.1. Approximate approaches for nonlinear filtering

For the sake of simplicity, it is “easier” to approximate the a posteriori PDF W_PS(x, t) than the nonlinearity at (4) and (9) [16, 17, 19]. In this sense, let us just list some of the approximate approaches for W_PS(x, t):

• Integral or global approximations for W_PS(x, t) [20];

• Functional approximations for W_PS(x, t) [16, 21];

• Higher Order Statistics (HOS) approximations for W_PS(x, t), and so on;

• Gaussian approximations: extended Kalman filter (EKF) [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]; unscented Kalman filter (UKF) [19]; quadrature Kalman filter (QKF) [17]; iterated Kalman filter (IKF), etc.

It is hardly feasible to give a complete overview of all those methods; moreover, not all of them are adequate, taking into account the observations introduced at the end of the previous section.

Let us start with the extended Kalman filter (EKF): considering W_PS(x, t) as a three-dimensional Gaussian PDF- ${\widehat{W}}_G\left(\mathit{x},t\right)$, from (9), it is possible to obtain the following equations for per-component of the mean estimates ${\left\{{\widehat{x}}_i\right\}}_1^3$ and for estimates of the elements of the a posteriori covariance matrix ${\left\{{\widehat{R}}_{\mathit{ij}}\right\}}_{i,j=1}^3$:

$$\begin{array}{c}\hfill {\dot{\widehat{x}}}_i=\underset{-\infty }{\overset{\infty }{\int }}\left({\widehat{\mathbf{\pi }}}^{\mathbf{T}}\left(\mathit{x},t\right)grad{x}_i\right)d\mathit{x}+\frac{1}{2}\left[\underset{-\infty }{\overset{\infty }{\int }}{x}_i\mathit{F}\left(\mathit{x},t\right){\widehat{W}}_G(\mathit{x},t)d\mathit{x}-{\widehat{x}}_i\underset{-\infty }{\overset{\infty }{\int }}\mathit{F}(\mathit{x},t\left){\widehat{W}}_G\right(\mathit{x},t)d\mathit{x}\right]\hfill \\ \hfill {\dot{\widehat{R}}}_{\mathit{ij}}=\underset{-\infty }{\overset{\infty }{\int }}\left({\widehat{\mathbf{\pi }}}^{\mathbf{T}}\left(\mathit{x},t\right)grad{{\displaystyle \stackrel{\circ }{x}}}_i{{\displaystyle \stackrel{\circ }{x}}}_j\right)d\mathit{x}+\frac{1}{2}\left[\underset{-\infty }{\overset{\infty }{\int }}{{\displaystyle \stackrel{\circ }{x}}}_i{{\displaystyle \stackrel{\circ }{x}}}_j\mathit{F}\left(\mathit{x},t\right){\widehat{W}}_G(\mathit{x},t)d\mathit{x}-{\widehat{R}}_{\mathit{ij}}\underset{-\infty }{\overset{\infty }{\int }}\mathit{F}(\mathit{x},t\left){\widehat{W}}_G\right(\mathit{x},t)d\mathit{x}\right],\hfill \end{array}$$E13

where ${{\displaystyle \stackrel{\circ }{x}}}_i=x_i-{\widehat{x}}_i$ and ${{\displaystyle \stackrel{\circ }{x}}}_j=x_j-{\widehat{x}}_j$.

The matrix form [14, 15, 16, 20] can be used to represent Eq. (13); however, for some specific applications, per-component representation (13) could be more adequate (see the following).

It is reasonable to assume convergence to the stationary values $\overline{R_{\mathit{ij}}}$ for $\forall {\widehat{R}}_{\mathit{ij}}\left(t\right)$ when t→∞, and as a result, the second equation in (13) can be expressed as a system of nonlinear algebraic equations, with standard numerical solutions. This consideration is relevant for real-time scenarios, as it significantly simplifies the implementation of the related EKF algorithms.

Functional approximation for W_PS(x, t) is, as it was described in [16, 21],

$$W_{\mathit{PS}}\left(\mathit{x},t\right)=\prod _{i=1}^3{W}_{\mathit{PS}}\left(x_i\right)\left[1+\sum _{q=2}^3\sum _{j=1}^{q-1}\frac{R_{\mathit{qj}}}{R_{\mathit{qq}}{R}_{\mathit{ji}}}\left(x_q-{\widehat{x}}_q\right)\left(x_j-{\widehat{x}}_j\right)\right].$$E14

From (14), we see that the functional approximation for the PDF is sufficiently non-Gaussian (marginal W_PS(x_i) is arbitrary), but for “joint” characterization of the vector $\widehat{\mathit{x}}$, only elements of the a posteriori covariance matrix ${\widehat{R}}_{\mathit{ij}}$ are considered.

It can be shown that the equations for ${\left\{{\widehat{x}}_i\right\}}_1^n$ and $\left\{{\widehat{R}}_{\mathit{ij}}\right\}$ coincide with those in (13), and the unique difference would be that one has to apply in (13) the approximation for W_PS(x, t) instead of ${\widehat{W}}_G\left(\mathit{x},t\right)$. The resulting integrals can be solved either through the Gauss-Hermit quadrature formula [17, 18] or analytically.

The integral or Global approximation for W_PS(x, t) is another approach for approximate solution. Maybe the experienced reader already noticed that the last two approximations for W_PS(x, t) can be considered as “local” as they offer maximum of W_PS(x, t), estimation of $\left\{{\widehat{x}}_i\right\}$, and $\left\{{\widehat{R}}_{\mathit{ij}}\right\}$.

For conditions of significantly large SNR, this is sufficient, but for low SNR, one has to find a different approach, known as integral approximation. This strategy was suggested as an adequate approximation of W_PS(x, t) together with the PDF’s “tails,” that is, for the whole span of x.

Let us suppose that W_PS(x, t) can be characterized as:

$$W_{\mathit{PS}}\left(\mathit{x},t\right)=W_{\mathit{PS}}\left(\mathit{x},\mathbf{\alpha }\left(t\right)\right).$$E15

Here α is an unknown vector of approximation parameters. As an approximation criterion for PDF, it is possible to use the Kullback measure; thus, one might obtain the following equation for the unknown vector α:

$${\displaystyle \dot{\mathbf{\alpha }}}=⟨L_{\mathit{PR}}^{+}\left\{\mathit{h}\left(\mathit{x},t\right)\right\}⟩+{\mathit{V}}^{-1}\left(t\right)⟨\mathit{h}\left(\mathit{x},t\right)\mathit{F}(\mathit{x},t)⟩,$$E16

where $\mathit{h}\left(\mathit{x},t\right)=\frac{\text{∂}ln{W}_{\mathit{PS}}\left(\mathit{x},\mathbf{\alpha }\left(t\right)\right)}{\text{∂}\mathbf{\alpha }}$, $\mathit{V}\left(t\right)=\underset{-\infty }{\overset{\infty }{\int }}{\left[\frac{\text{∂}ln{W}_{\mathit{PS}}\left(\mathit{x},\mathbf{\alpha }\left(t\right)\right)}{\text{∂}\mathbf{\alpha }}\right]}^T{W}_{\mathit{PS}}\left(\mathit{x},\mathbf{\alpha }\left(t\right)\right)\mathit{dx}=〈\frac{{\text{∂}}^2{W}_{\mathit{PS}}\left(\mathit{x},\mathbf{\alpha }\left(t\right)\right)}{\text{∂}\mathbf{\alpha }\text{∂}{\mathbf{\alpha }}^{\mathbf{{\rm T}}}}〉$, $L_{\mathit{PR}}^{+}\left\{•\right\}$ is a self-adjoint operator to the FPK operator [22].

Now, as an integral approximation of W_PS(x, α(t), let us choose the so-called “Dynkin PDF” with α(t) is the vector of sufficient statistics for W_PS(⋅):

$$W_{\mathit{PS}}\left(\mathit{x},\mathbf{\alpha }\left(t\right)\right)=Cexp\left\{\sum _{p=1}^K{\alpha }_p\left(t\right){\mathbf{\varphi }}_p\left(\mathit{x}\right)+{\mathbf{\varphi }}_0\left(\mathit{x}\right)\right\},$$E17

where {φ_p(x)} are orthogonal multidimensional operators: Laguerre, Hermite, and so on.

One can notice that there is a significant coincidence between (17) and the orthogonal series characterization of W_PS(x, α(t) [22]: even though both apply series of orthogonal functions, in (17), it is not used for W_PS(x, α(t)) but for its monotonical transform ln{W_PS(x, α(t)}. So, the coefficients {α_p(t)} can be expressed by means of the cumulants of W_PS(x) [22]. Thanks to this, instead of searching for a solution of (17), hardly possible in an analytically way, one can search directly equations for the cumulants (HOS) of W_PS(x, t) [16, 26].

Here, the HOS approach will be presented because the last problem was addressed in the cited references. It is worth noticing the following: for real-time scenarios when n > 1, equations for HOS and Eq. (16) are significantly complex; for n = 1, both strategies are equivalent [26].

3. Multimoment filtering of chaos

As it follows from the material of Section 2, all the algorithms are “one-moment” in the sense that they are operating only with the data at each time instant, that is, they are tracking instantaneously one moment magnitude of the received aggregate signal. As it was shown at [27], the adequate filtering algorithm (for the one-moment case) is an Extended Kalman Filter (EKF).

This choice is more or less expected, due to the experience which is already known from the available references (see above). EKF shows rather good performance for the filtering of chaotic signals: the mean squared error (MSE) is less than 1% when SNR is about −3 dB, and for SNR bigger than −3 dB, the results are much better.

In this regard, a question arises: is it possible to improve this approach in the sense of getting still rather good MSE’s for successively lower thresholds of the SNR with an algorithm of reasonable complexity? The following material attempts to prove that the answer is “yes,” if one can apply some additional information from the received aggregate signal taken on several sequential time instants.

It means that the information has to be considered in the block manner by aggregating data, in our case, for several time instants ([8, 16, 27], and so on.). The difference between the following approach and that from the cited references is precisely the aggregated data obtained for many time instants: multimoment algorithms are carried out through the generalization of the Stratonovich-Kushner equations (SKE) for the corresponding multimoment data, and therefore, in the following, all heuristics for the simplification, considered as Generalized SKE (GSKE), are not arbitrary but can be taken as generalized heuristics from the “standard” one-moment SKE (see below). This gives a “hope” to achieve the abovementioned improvement for the SNR threshold with less complex tools.

It follows from the fact that, as it was shown in [8] (see also the references therein), the GSKE comes from the same structure as its one-moment prototype. So the way of its simplification (except for the limiting of the number of time instants) in order to get a quasi-optimum algorithm, could be done in a similar way as for the one-moment case: approximation of the a posteriori PDF (characteristic function) in SKE with a minimum set of significant parameters. Moreover, there is an additional way to improve the accuracy of the quasi-optimum solution for the GSKE: assume this quasi-optimum algorithm as a “given structure,” as it was proposed in [16] and also considered in the following.

3.1. Generalization of SKE for the multimoment case

In the same way, as it was underlined earlier, the chaos is “generated” by the equation:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right),\mathit{x}\in R^n,\mathit{x}\left(t_0\right)={\mathit{x}}_0,$$E18

where f(•) = [f₁(x), …f_n(x)]^T is a differentiable vector function and it can be considered as a degenerated Markov process from the following stochastic equation:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right)+\mathbf{\epsilon \xi }\left(t\right),$$E19

where ξ(t) is a vector of “weak” external white noise with the related positively defined matrix of “intensities” ε = [ε_ij]^n×n.

In the following, one can consider both the ordinary differential equation (ODE) (18) and the stochastic differential equation (SDE) (19) when the noise intensities tend equally to zero. Adding the ε term in (19) guarantees the existence of a stationary PDF for x(t) as well, no matter how small the elements of ε might be [28]. So, one can suppose that this stationary PDF, W_ST(x), is known a priori. For our case in practical sense, one can deal actually only with the stationary PDF, which we assumed is modeled by means of a chaotic process (concretely let us say the first component, x₁(t), of certain attractor model). Certainly W_ST(x₁) can easily be obtained from W_ST(x). If the two PDFs coincide in terms of certain fitness criteria, then only for simplicity in the subsequent developments, the SDE (19) can be substituted by its statistically equivalent one-dimensional SDE with the same W_ST(x₁):

$${\dot{x}}_1=f\left(x_1\right)+\sqrt{\epsilon }\xi \left(t\right),$$E20

where $f\left(x_1\right)=\frac{\epsilon }{2}\frac{d}{{\mathit{dx}}_1}ln{W}_{\mathit{ST}}\left(x_1\right)$ and ε in (20) can be considered here as a “scale factor” and can be chosen by equalizing the average powers of real x₁(t) and solution of (20). Formally, there is no need for all those operations, but then the reader has to be extremely concentrated with “multiindex” definitions: one index for the number of applied components of the attractor and another index for the time instant, that is, x^m(t_i), which might cause confusion in further developments, as x₁(t) is an observable component whose dynamics depends on other “nonobservable” components. For those reasons, in the following, (20) will be considered as a model of the desired signal for filtering.

Let us introduce the following notation for the time instants (time moments): t₁ < t₂ < t₃ … <t_n and x_i = x(t_i), $i=\overline{1,n}$. Then, ${\left\{x\left(t_i\right)\right\}}_1^n$ forms a vector x(t) = [x(t₁),…, x(t_n)]^T and W_n(x, t) ≅ W_n(x₁, …x_n; t₁, …t_n); W_n(x,t) is an a priori PDF for x(t). As it follows from ([16], ch. 5):

$$\frac{\partial W_n\left(\mathit{x},t\right)}{\partial t_i}=L_i\left\{W_n\left(\mathit{x},t\right)\right\}$$E21

where $L_i\left\{•\right\}=-\frac{\text{∂}}{\text{∂}{x}_i}{K}_1\left(x_i\right)+\frac{1}{2}\frac{{\text{∂}}^2}{\text{∂}{x}_i^2}{K}_2\left(x_i\right)$ is the FPK operator [16] with K₁(x_i) = f₁(x_i), K₂(x_i) = ε². It is easy to show that by consecutive differentiation one can obtain:

$$\frac{{\partial }^n{W}_n\left(\mathit{x},t\right)}{\partial t_1\dots \partial t_n}=\prod _{i=1}^n{L}_i\left\{W_n\left(\mathit{x},t\right)\right\},$$E22

$$L_{\mathit{PR}}\left\{•\right\}=\prod _{i=1}^n{L}_i\left\{•\right\}.$$E23

Certainly, the adjoint operator [16, 22] for the multimoment case is:

$${\widehat{L}}_{\mathit{PR}}\left\{•\right\}=\prod _{i=1}^n{\widehat{L}}_i\left\{•\right\},$$E24

where ${\widehat{L}}_i=K_1\left(x_i\right)\frac{\text{∂}}{\text{∂}{x}_i}+\frac{K_2\left(x_i\right)}{2}\frac{{\text{∂}}^2}{\text{∂}{x}_i^2}$ is a Kolmogorov operator [22].

Let us then introduce the a posteriori n-dimensional PDF W_ps(y|x,t) ≅ W_ps(x,t) for the multimoment case. Then, repeating formally the development for the SKE, but in this case generalized for the “n” time case (in the same way as it was done at [27]), one can get:

$$\frac{{\partial }^n{W}_{\mathit{ps}}\left(\mathit{x},t\right)}{\partial t_1\dots \partial t_n}=L_{\mathit{PR}}\left\{W_{\mathit{ps}}\left(\mathit{x},t\right)\right\}+\frac{1}{2}\left[\mathit{F}\left(\mathit{x},t\right)-\underset{R^n}{\int }\mathit{F}(\mathit{x},t\left)W_{\mathit{ps}}\right(\mathit{x},t)d\mathit{x}\right]W_{\mathit{ps}}\left(\mathit{x},t\right)$$E25

with t = [t₁, …, t_n]^T,

$$\mathit{F}\left(\mathit{x},t\right)=\frac{{\left[\mathit{y}\left(t\right)-\frac{1}{2}\mathit{x}\left(t\right)\right]}^T}{N_0}\left[\mathit{y}\left(t\right)-\frac{1}{2}\mathit{x}\left(t\right)\right],$$E26

where y(t) = [y(t₁), …, y(t_n)]^T is the vector of ${\left\{x\left(t_i\right)\right\}}_1^n$ taken from y(t) = x(t) + n(t) and n(t) is the AWGN with intensity N₀.

Analyzing (25) by comparing it with the standard form of the SKE (see Eqs. (9) and (10) in part II), one can see that there is a total “structural” identity! The same matter takes place for the a posteriori cumulants [16, 27], that is:

$$\frac{\partial {\kappa }_{r_1,\dots r_n}^{\mathit{ps}}\left(t\right)}{\partial t_1\dots \partial t_n}={\left(-j\right)}^k\frac{{\partial }^k}{\partial {\lambda }_1^{r_1},\dots \partial {\lambda }_n^{r_n}}\left\{{\left[\frac{M\left\{\widehat{L}exp\left\{j{\mathbf{\lambda }}^T\mathit{x}\right\}\mathit{F}(\mathit{x},t)\right\}}{M\left\{exp\left\{j{\mathbf{\lambda }}^T\mathit{x}\right\}\right\}}\right]}_{\lambda =0}+{\left[\frac{M\left\{exp\left\{j{\mathbf{\lambda }}^T\mathit{x}\right\}\mathit{F}(\mathit{x},t)\right\}}{M\left\{exp\left\{j{\mathbf{\lambda }}^T\mathit{x}\right\}\right\}}\right]}_{\lambda =0}\right\}$$E27

where r₁+ r₂ + … + r_n = k, k = 1, 2, ….

One can see from (25) and (26) that those algorithms are rather complex for implementation in real-time regime. So, in addition to the one-moment SKE, they have to be modified in order to get the quasi-optimum solution.

3.2. Quasi-optimum solutions. Generalized EKF

One has to know that “quasi-optimum” solutions (for any problem) are based on some heuristics and those heuristics have to be reasonable and based on previous experience in solving similar problems. In the case of multimoment filtering, the analogies can be the following (of course implicit considerations for complexity have always to be taken into account):

• The priority will be given to the quasi-linear approximation for nonlinear functions in the same way as it was assumed for the “standard” one-moment filtering.

• All algorithms for block processing show that there is in some sense a reasonable block length for the processed data. Taking into account the complexity limits and that the covariance function of the chaos initially drops rather fast [29], let us take first n = 2.

• The approximation of the a posteriori PDF (characteristic function) has to apply the minimum set of first cumulants; one has to remind that, as the order of cumulants grows, their significance for PDF approximation vanishes [22];

Taking these observations into account, let us take n = 2, that is, two-moment filtering case, then [16, 22]:

$${\theta }_{\mathit{ps}}\left(\mathbf{\lambda }\right)\cong exp\left\{\sum _{s=1}^2\frac{\mathit{js}}{s!}\sum _{r_1,r_2}^S{\kappa }_s\left(t_1,t_2\right){\mathbf{\lambda }}^{r_1}{\mathbf{\lambda }}^{r_2}\right\},$$E28

and cumulants are:

$${\kappa }_{r_1,\dots r_n}^{\left(\mathit{ps}\right)}\left(t\right)={\left(-j\right)}^k{\left[\frac{{\partial }^k}{\partial {\lambda }_1^{r_1},\dots \partial {\lambda }_n^{r_n}}ln{\theta }_{\mathit{ps}}\left(\mathbf{\lambda }\right)\right]}_{\lambda =0}.$$

Another assumption is that the a posteriori process is supposed to be stationary; then, the one-moment cumulants for t₁ and t₂ have to be the same and the only mutual cumulant taken into account might be κ₁₁(t₁, t₂). Next, for each moment “t₁” and “t₂” one-moment cumulants can be calculated applying Gaussian approximation for the a posteriori PDF, and for the two-moment case, the “functional approximation” could be applied. In a rigorous sense, the a posteriori variance ${\kappa }_2^{\mathit{ps}}$ has to be evaluated as ${\kappa }_2^{\mathit{ps}}\left(t_1,t_2\right)$, considering the covariance among time instants “t₁” and “t₂”; in the following, the heuristic strategy will be introduced, which avoids the cumbersome calculations.

One can obtain the first two-moment cumulants:

$$\begin{array}{l}{\dot{\kappa }}_1=<K_1\left(x\right)>+\frac{1}{2}<\mathit{xF}\left(x,t\right)>-\frac{{\kappa }_1}{2}<\mathit{xF}\left(x,t\right)>\\ {\dot{\kappa }}_2=<2{\mathit{xK}}_1\left(x\right)>-{\kappa }_1<K_1\left(x\right)>+\epsilon +\frac{1}{2}<{\left(x-{\kappa }_1\right)}^{2}F\left(x,t\right)>-\frac{{\kappa }_2}{2}<\mathit{xF}\left(x,t\right)>,\end{array}$$E29

where < > is a symbol for the averaging procedure, $F\left(x,t\right)=\frac{1}{N_0}\left[y\left(t\right)-\frac{1}{2}x\left(t\right)\right]$, and K₁(x) is the drift coefficient for (19).

One has to notice that at (29) κ₁(t) is an estimation of the filtered signal (in our case, it is a chaotic signal); κ₂(t) is a measure of the filtering accuracy. As it can be is seen from (29), those equations were written without any intention for linearization, that is, they are presented in a generalized form. For the quasi-linear algorithms, it is well known [27] that κ₂(t)/N₀ is the main part for the “averaging coefficient” of the second element in the first quasi-linear equation of (29), that is, it is an averaging value for the instantaneous information actualization from the entering desired signal plus noise. Thus, if one can reduce κ₂ through the two-moment processing, the accuracy of the quasi-linear method will grow and the challenge stated before will be almost solved. To achieve the latter, one can take into account that κ₂(t) in the stationary regime is oscillating around its stationary value $\overline{{\kappa }_2\left(t\right)}=\underset{t\to \infty }{lim}{\kappa }_2\left(t\right)$ which is commonly assumed as an accuracy measure in the one-moment case.

The value of $\overline{{\kappa }_2\left(t\right)}$ can be diminished applying the information from κ₁₁(t) also in the stationary case, that is, ${\overline{\kappa }}_{11}=\underset{t_1,t_2\to \infty }{lim}{\kappa }_{11}\left(t_1,t_2\right)$; then, it is known that ${\widehat{\overline{\kappa }}}_2={\overline{\kappa }}_2\left(1-{\overline{\kappa }}_{11}^2\right)$ and it is always less than ${\overline{\kappa }}_2$, if and only if the ${\overline{\kappa }}_{11}\ge 0$; In this way, ${\widehat{\overline{\kappa }}}_2$ can be used as a new weighting coefficient in (29). To find κ₁₁(t₁,t₂) from (27), some cumbersome developments are required which finally yield to:

$$\frac{\partial {\kappa }_{11}\left(t_1,t_2\right)}{\partial t_1\partial t_2}=<K_1\left(x_1\right)K_1\left(x_2\right)>+<\frac{x_1{x}_2}{2}F\left(\mathit{x},t\right)>-\frac{{\kappa }_{11}\left(t_1,t_2\right)}{2}<F\left(\mathit{x},t\right)>$$E30

and

$${\overline{\kappa }}_{11}=\frac{2\left[<\frac{x_1{x}_{2}F\left(\mathit{x},t\right)}{2}+2K_1\left(x_1\right)K_1\left(x_2\right)>\right]}{<F\left(\mathit{x},t\right)>}.$$E31

First we would like to stress here that, as we are interested in covariance calculation, it is necessary to preserve the notations x(t₁) = x₁ and x(t₂) = x₂. Second, we want to “improve” the stationary value ${\overline{\kappa }}_2$ evaluated for the one-moment case through its indirect dependence on ${\overline{\kappa }}_{11}$ as if it was “evaluated” for the two-moment case.

Thus in doing so, the direct calculation of the quasi-linear algorithm for the two-moment case is bypassed (see (29) and (30)). For applications in real time, the formal calculus is almost impossible. Instead, we simplified it with a formal “ignorance” of the two-moment features. There might be for sure a compromise between the complexity and the improvement attempt for the “classic” EKF.

In order to avoid some additional complexities for the calculation of (31), let us make the following assumption: introduce the SNR of the filtering in the way: $h^2=\frac{{{\displaystyle \overline{\kappa }}}_1^2}{N_0}<<1$, that is, weak signal case. In this regard [16, 27], the a priori data are the main influence, that is, approximately only <K₁(x₁) K₁(x₂) > can be applied. Or one can simply apply a Gaussian approximation for the second equation in (29) for the stationary regime (${\dot{\kappa }}_2\cong 0$)

$$2\left[\overline{K_1\text{'}\left({\kappa }_1\right)}+{\overline{\kappa }}_{11}\right]\frac{{\overline{\kappa }}_2}{2}+\epsilon +\frac{1}{4}\overline{F\text{'}\text{'}\left({\kappa }_1\right)}{\overline{\kappa }}_2^2=0.$$E32

In the case h² < 1, it is possible to achieve:

$${\overline{\kappa }}_2\sim \frac{1}{2\left[\overline{K\text{'}\left({\kappa }_1\right)}+{\overline{\kappa }}_{11}\right]},$$E33

and if ${\overline{\kappa }}_{11}$ > 0, and K′(κ₁) ≥ 0, κ₂ is always reduced compared with the one-moment approach. Formula (33) can be seen as another illustration about the usefulness of the heuristic approximation proposed above. Then, to evaluate the order of the ${\overline{\kappa }}_{11}$, let us apply for averaging of <K₁(x₁) K₁(x₂)> the functional approximation of W_ps(x₁, x₂) in the way:

$$W_{\mathit{ps}}\left(x_1,x_2\right)=\frac{1}{2\pi {\overline{\kappa }}_2}exp\left[-\frac{{\left(x_1-{\kappa }_1\right)}^2}{2{\overline{\kappa }}_2}\right]exp\left[-\frac{{\left(x_2-{\kappa }_1\right)}^2}{2{\overline{\kappa }}_2}\right]\left[1+{\overline{\kappa }}_{11}\left(x_1-{\kappa }_1\right)\left(x_2-{\kappa }_1\right)\right].$$E34

As an approximate result, one can substitute (34) in (33), assume h² < 1 and see that the normalized value ${\overline{\kappa }}_{11}$ has the same order as h², that is, ${\overline{\kappa }}_{11}$ ∼ O(h²). This is an important consideration because usually the pure chaos has a low covariance interval [29] and one can obtain a very small MSE for two time instants t₁ and t₂ arbitrarily close. In this sense and fixing SNR ∼ 0.5 and MSE ∼ 0.1%, an equivalent MSE can be reached using the two-moment approach but with an SNR threshold 30% lower than for the one-moment case. Let us be emphatic and say that the approximation ${\overline{\kappa }}_{11}$ ∼ O(h²) is valid just for h² < 1, and calculation of ${\widehat{\overline{\kappa }}}_2\sim {\overline{\kappa }}_2\left(1-{\overline{\kappa }}_{11}^2\right)$ has to be updated instantaneously because h² is varying in the interval 0 ≤ h² < 1.

Of course this calculation is quite approximated and true superiority for the two-moment case of the modified quasi-linear strategy has to be verified by computer experiments. Anyway it is a strong sign indicating that the use of the two-moment strategy can be very opportunistic if and only if one can find strategies to reduce the computational complexity, for example, the generalized extended Kalman filter (GEKF) algorithm.

Finally, let us reiterate that the GEKF is yet a one-moment strategy for quasi-optimum filtering, but internally makes processing of the statistical features of the chaotic data (input) through the multimoment (two-moment) apparatus. That is why this modified GEKF improved accuracy in comparison with the standard EKF. In the following in order to additionally improve the accuracy of this one-moment modified EKF, it is convenient to apply the principles of the theory of so-called “conditionally optimum filtering” proposed in ([16], ch. 9), taking this generalized EKF as the “tolerance” or “admitted” filter.

4. Conditionally optimum filtering approach

The ideas and methods for conditionally optimum filtering are rather simple and are thoroughly described at ([16], ch. 9). So, let us first present the basic idea of this method. In the general case, the conditional optimum filter for the optimum estimation of the desired signal x(t) in presence of AWGN n(t) can be presented in the form [16]:

$${\dot{\kappa }}_1=\mathit{\alpha \xi }\left(y,{\kappa }_1,t\right)+\mathit{\beta \eta }\left(y,{\kappa }_1,t\right)y\left(t\right),$$E35

where κ₁(t) is a filtered signal; y(t) = x(t) + n(t); n(t) is the AWGN with intensity N₀; α, β are some time-dependent coefficients which have to be found.

The representation (35) is a generalized representation of the filtering algorithms where ${\dot{\kappa }}_1$ is the expectation of the filtered signal. It is clear as well [16] that this form is valid also for the quasi-optimum nonlinear filtering algorithms. In the previous part, a modified EKF algorithm was proposed for the two-time-moment case, which shows rather opportunistic improvement of the filtering accuracy, applying some heuristics related to the simplified implementation of the two-moment principle of filtering. Sure those simplifications do not allow taking full advantage of the application of the two-moment principle. Once again, this simplification is reasonable for diminishing the dimension of the filtering algorithm in order to make it practical for real-time applications. Therefore, the hope for further improvement of the characteristics of this modified EKF might be based on further optimization in the framework of conditional optimality [16].

In the theory of conditional optimality, the structure of the filter is already chosen (in our case, it is the GEKF) and the only chance for further accuracy improvement is to optimize the coefficients α(t) and β(t) in order to minimize the MSE. The structure which was chosen initially is a so-called admitted structure which actually belongs to a class of the admitted filters. The next step is to minimize the MSE. The minimization of the MSE is a strategy in which the admitted filter makes an optimal transition at the moment “s” (s > t, s → t) from an initial stage, at moment “t,” to a new stage at the moment “s” with the minimum MSE. The algorithm of such kind of filter is “conditionally optimum” according to Ref. [16].

Hereafter we are not going to present all the material related to this approach as it was comprehensively described at ([16], ch. 9), we will only apply the necessary final formulas from there. Unfortunately, full use of the abovementioned approach is not possible (as we will see in the following), and so, we will present some developments that allow to obtain the coefficients α(t) and β(t) successfully.

4.1. Approach to find unknown coefficients α(t) and β(t)

It is possible to present an admitted structure of the conditionally optimum filter from (29) in two equivalent forms:

$${\dot{\kappa }}_1=\alpha \left[K_1\left({\kappa }_1\right)+\frac{{\widehat{{\displaystyle \overline{\kappa }}}}_2}{2}{K_1}^{\text{'}\text{'}}\left({\kappa }_1\right)\right]+\beta \frac{{\widehat{\overline{\kappa }}}_2}{N_0}\left[y\left(t\right)-{\kappa }_1\left(t\right)\right]$$E36

$${\dot{\kappa }}_1=\alpha \left[K_1\left({\kappa }_1\right)+\frac{{\widehat{\overline{\kappa }}}_2}{2}{K_1}^{\text{'}\text{'}}\left({\kappa }_1\right)-\frac{{\widehat{\overline{\kappa }}}_2{\kappa }_1}{N_0}\right]+\beta \frac{{\widehat{\overline{\kappa }}}_2}{N_0}y\left(t\right),$$E37

where, as it was proposed earlier,

$${\widehat{\overline{\kappa }}}_2={\overline{\kappa }}_2\left(1-{\overline{\kappa }}_{11}^2\right).$$E38

Then, from (36) and (37), one has

$$\xi \left(t\right)=K_1\left({\kappa }_1\right)+\frac{{\widehat{\overline{\kappa }}}_2}{2}{K_1}^{\text{'}\text{'}}\left({\kappa }_1\right),\eta \left(t\right)=\frac{{\widehat{\overline{\kappa }}}_2}{N_0}\left[y\left(t\right)-{\kappa }_1\left(t\right)\right]$$E39

$$\xi \left(t\right)=K_1\left({\kappa }_1\right)+\frac{{\widehat{\overline{\kappa }}}_2}{2}{K_1}^{\text{'}\text{'}}\left({\kappa }_1\right)-\frac{{\widehat{\overline{\kappa }}}_2{\kappa }_1}{N_0},\eta \left(t\right)=\frac{{\widehat{\overline{\kappa }}}_2}{N_0}y\left(t\right).$$E40

One can see that in this regard, α and β are weighting coefficients of a priori information related to the desired chaotic signal and a posteriori data. This issue was thoroughly commented in [27]. For SNR < 1, the weight of ξ(t) obviously prevails, because a posteriori data are strongly corrupted by the additive noise. Nevertheless, taking into account that ${\widehat{\overline{\kappa }}}_2$ is rather small for the modified EKF, in the following, ${\widehat{\overline{\kappa }}}_2$ (which is actually the MSE) will be considered as a “small parameter” in all the approximations.

In order to follow all definitions and notations from ([16], ch. 9), one has to use the Ito form in all the equations:

$$\begin{array}{l}\mathit{dy}=\mathit{Xdt}+{\mathit{dW}}_1={\varphi }_1\left(y,x,t\right)\mathit{dt}+{\psi }_1\left(y,x,t\right){\mathit{dW}}_1\\ \mathit{dx}=f\left(x\right)\mathit{dt}+{\mathit{dW}}_2={\varphi }_1\left(x,t\right)\mathit{dt}+\psi \left(x,t\right){\mathit{dW}}_2,\end{array}$$E41

where {W_i(t)} are independent Wiener processes, i = 1, 2. It is obvious that:

$$\begin{array}{l}{\phi }_1\left(x,t\right)=f\left(x\right)={\kappa }_1\left(x\right)\\ {\phi }_1\left(y,x,t\right)=x\\ {\psi }_1\left(y,x\right)=1\\ \psi \left(x,t\right)=1\end{array}$$E42

Then, from ([16], ch. 9)

$${\widehat{x}}_s-{\widehat{x}}_t\cong {\kappa }_s-{\kappa }_t={\mathit{\alpha \xi }}_t\mathrm{\Delta t}+{\mathit{\beta \eta }}_t\left({\varphi }_{1_t}\mathrm{\Delta }t+{\varphi }_{1_t}\mathrm{\Delta }W\right).$$E43

Unbiased conditions for the optimum estimation from (43) are [16]:

$$\alpha <{\xi }_t>+<{\eta }_t{\varphi }_{1_t}>-<{\varphi }_1>=0.$$E44

Taking ξ_t and η_t according to its definitions from (40), it is easy to get from (44):

$$\alpha m_1+\beta m_2=m_0,$$E45

where m₀ = <φ_t>, m₁ = <ξ_t>, m₂ =  < η_tφ_1t>.

Taking into account (42) with conditions ${\widehat{\overline{\kappa }}}_2$ < 1 and assuming that K₁(κ₁) ≈ K₁´´(κ₁) ≈ 0¹, finally one gets:

$$\frac{\beta }{\alpha }=\frac{{\kappa }_2}{<x^2>}.$$E46

The next step, as it was proposed in ([16], ch. 9), is focused on checking the correlation conditions for the error (κ_s−x_s) with the vector [ξΔt, ηΔy] which yields to [16]:

$$\beta ={\kappa }_{0_2}{\kappa }_{22}^{-1},$$E47

where

$${\kappa }_{0_2}=<\left(x_t-{\kappa }_{1_t}\right)x_t\frac{{\widehat{\overline{\kappa }}}_2}{N_0}{y}_t>+<{\eta }_t\cdot y_t>\frac{{\widehat{\overline{\kappa }}}_2}{N_0},{\kappa }_{22}={\left(\frac{{\widehat{\overline{\kappa }}}_2}{N_0}\right)}^2<y_t^2\cdot {\eta }_t>.$$E48

From the second equation in (48), it follows that β → ∞ which is a clear absurd. So, why this happened and what is wrong? Is the approach in ([16], ch. 9) wrong? Definitively, no. It is possible to show that the estimate κ₁ is unbiased and decorrelated with both components ξ(t) and η(t), but for our special case, the condition that κ₂₂ (a matrix in the general case) has to be invertible is violated. Opposed as it was stated in ([16], ch. 9), the approach is not working.

The solution might be found from direct calculation of (x−κ₁) from the SDE of chaos and (29) and by minimization of <(x-κ₁)² > by α or β.

4.2. Direct evaluation of the MSE and its minimization

As a first step, let us calculate the difference between the solution of (20) and (39) by applying (46):

$$\left(x-{\kappa }_1\right)=\underset{0}{\overset{T}{\int }}\left\{\left[K_1\left(x\right)-\alpha K_1\left({\kappa }_1\right)\right]-\frac{{\mathit{\alpha \kappa }}_1{\widehat{\overline{\kappa }}}_{2}n\left(t\right)}{<x^2>N_0}\right\}\mathit{dt}.$$E49

Let us take the second power of (49) and make a statistical average. One has to notice that the second power of (49) is a double integral and <n(t₁) n(t₂)> = N₀δ(t₂−t₁). Then, applying finally the assumption ${\widehat{\overline{\kappa }}}_2$ < 1, one can get for the MSE:

$$\mathit{MSE}\approx <K_1^2\left(x\right)>+{\alpha }^2<K_1^2\left({\kappa }_1\right)>-2\alpha <K_1\left(x\right)K_1\left({\kappa }_1\right)>+\frac{{\alpha }^2{\widehat{\overline{\kappa }}}_2}{<x^2>N_0}.$$E50

Looking for the minimum of (50) in terms of “α”, one easily finds:

$$\alpha =\frac{<K_1\left({\kappa }_1\right)K_1\left(x\right)>}{<K_1^2\left({\kappa }_1\right)>+\frac{{\widehat{\overline{\kappa }}}_2}{<x^2>}}.$$E51

Assuming that still ${\widehat{\overline{\kappa }}}_2$ is a “small parameter,” it follows that α ≈ 1 and $\beta \cong \frac{{\kappa }_1}{<x^2>}\cong \mathrm{O}\left(\frac{1}{{\kappa }_1}\right)$. In this regard,

$$\mathit{MSE}\sim \frac{{\left({\widehat{\overline{\kappa }}}_2\right)}^2}{<x^2>}.$$E52

Comparing Eq. (52) with the MSE of the one-moment filtering which is κ₂, one can see that the conditional optimum filtering might significantly improve the MSE with the same SNR or significantly diminish the SNR threshold for a fixed MSE.

The authors consider that the two-moment filtering of chaos together with the conditionally optimum principle is a very opportunistic approach to significantly improve the MSE for chaos filtering.