Binary-encoding-based quantized Gaussian particle filtering for nonlinear systems with bit-flipping

Abstract

This paper investigates the state estimation problem for nonlinear systems with binary-encoding-based quantization, bit-flipping and heavy-tailed noise. Due to the limitation of the communication channel bandwidth in a network environment, the measurement data collected by sensors must be quantized based on a limited number of bits and converted into binary codes for transmission, thus generating quantization errors. Meanwhile, during data transmission, the binary code output from quantization is affected by channel noise, signal distortion, signal interference, and other factors. It is prone to bit-flipping (i.e. 0 is flipped to 1, and 1 is flipped to 0), which leads to communication errors. For the problems of quantization error, heavy-tailed noise and bit-flipping, this paper aims to achieve the state estimation of such nonlinear systems by improving the particle filter algorithm, with the aid of the Bayesian formula and the Monte Carlo simulation method. To address the issue of particle degradation in particle filtering, this study employs a Gaussian approximation of the posterior probability density as the proposal distribution. In the design of the proposal distribution, a normalized innovation sequence is introduced to mitigate quantization errors and bit-flip effects. Based on this formulation, the mean, covariance of the proposal distribution, and likelihood function are analytically derived, yielding a complete proposal distribution. Subsequently, particles are sampled from the constructed proposal distribution, and an updated particle weighting scheme is rigorously derived. The optimal state estimate is obtained through a weighted summation of the particles. Numerical simulations are conducted to validate the efficacy of the proposed algorithm.

Introduction

With the rapid development of modern science and technology, nonlinear systems play an increasingly important role in many fields such as aerospace, artificial intelligence, communication, and control systems.^1–14 Since the dynamic characteristics of nonlinear systems are complex and highly uncertain, it is often accompanied by the influence of unknown factors. Traditional linear state estimation methods are inadequate satisfy practical requirements. Because of this, scholars have carried out a lot of explorations to face the state estimation problem of nonlinear systems.

The Kalman Filter (KF) algorithm is one of the classical methods for state estimation problems.^15,16 Using the minimum mean square error method, the KF provides an effective recursive calculation method to estimate the system state. Buck et al. ¹⁷ Sunahara et al. ¹⁸ and Zheng et al. ¹⁹ proposed the well-known Extended Kalman Filter (EKF) algorithm based on the KF algorithm. However, the KF algorithm requires that the system be linear, which limits its application. As an early extension of the KF, the EKF adapts the linear framework to nonlinear systems via first-order Taylor series expansion and linearization. However, in the case of strong nonlinearity, neglected high-order terms often introduce significant errors, leading to filter divergence. Additionally, the computation of the Jacobian matrix is cumbersome. To mitigate parameter estimation bias arising from the linearization of nonlinear systems, the Unscented Kalman Filter (UKF) algorithm based on the Unscented Transform (UT) principle was proposed.^20,21 The core of the UKF algorithm is to approximate a probability distribution. However, EKF and UKF in the state estimation field require high accuracy of noise distribution and measurement information. To overcome the limitations of the above filtering methods regarding system model characteristics and noise distribution, subsequent researchers proposed particle filter algorithms. Particle filtering ²² is based on Monte Carlo sampling and approximates the posterior probability distribution of the system state using a large number of weighted particles. Theoretically, particle filtering can effectively solve the state estimation problem of arbitrary nonlinear and non-Gaussian systems. Moreover, since the Monte Carlo simulation approach relies on the law of large numbers, it can be readily extended to high-dimensional scenarios. After years of development, particle filtering algorithms have been gradually improved, which has led to the formation of a complete theoretical framework. However, particle filtering still faces several challenges, such as particle degradation, sample impoverishment, high computational complexity, and real-time performance issues. ²³ To address these challenges, researchers across various disciplines have proposed a range of optimization methods, thereby significantly advancing the development of particle filtering techniques.

As a fundamental technique in modern communication and control systems, signal quantization has recently been employed extensively in fields such as nonlinear control, distributed learning, and networked systems.^24–30 Owing to limited communication channel bandwidth in networked environments, analog sensor data are typically quantized into binary codes before transmission to a remote estimation center. Improper quantization can significantly degrade system performance and even cause divergence due to quantization errors. Therefore, signal quantization has become a focus area for many researchers. For example, dynamic quantizers designed for discrete-time nonlinear systems can adjust quantization parameters in real time using the adaptation parameter μ ³¹ to handle to signal variations. In parallel, they synchronously optimize controller and quantizer parameters to suppress the effect of quantization error through matrix inequalities. However, traditional deterministic truncated quantization methods suffer from intrinsic limitations, frequently resulting in considerable quantization errors. As demonstrated in,^32–34 recent studies have extensively investigated probabilistic quantizers as a more robust and adaptive alternative to address this challenge. The core operating mechanism of probabilistic quantizers ³⁵ involves partitioning the probability space into multiple quantization intervals. The continuously varying measurements are subsequently mapped to corresponding quantized values based on a predefined probabilistic mapping rule. Unlike traditional deterministic quantization methods, probabilistic quantizers not only focus on the magnitude of input values but also consider the probability of their occurrence. As a result, the statistical characteristics of the original data are better preserved during the quantization process. Notably, each quantized output is a random variable whose expected value equals the corresponding original signal value. This property fundamentally ensures the unbiasedness of the quantization process and provides a solid guarantee for subsequent signal processing and system operation. ³⁴ Shen et al. ³⁶ focused on the H₂ control of linear systems with multiple quantization channels and designed a composite controller to meet performance requirements. Another work was concerned with the mismatched quantized $H_{\infty }$ output-feedback control of fuzzy Markov jump systems via a dynamic guaranteed cost triggering scheme. ³⁷ A brief study investigated dynamic event-triggered iterative filter and fault compensation control for PMSG-based wind turbine systems under deception attack. ³⁸

The bit-flipping mechanism is a fundamental technique that alters the values of bits in binary data, that is, converting 0–1 and vice versa, to achieve specific functional objectives. The bit-flipping mechanism is a core error correction technique for information processing systems. In recent years, it has seen significant development in quantum computing, classical coding, and adversarial communication.^39–42 In the context of state estimation, most existing studies on quantization-based methods implicitly assume that the quantized output bit string can be transmitted without error. In practical scenarios, bit-flipping errors frequently occur during signal transmission due to channel noise, distortion, interference, and transceiver clock synchronization mismatches. These errors inevitably lead to communication failures, as illustrated in Figure 1. Given this, these challenges have attracted considerable research attention, and notable progress has been made. In linear discrete-time-varying systems, Liu et al. ³⁵ derived the Binary Coded Quantized Kalman Filter (BQKF) algorithm to model the effects of bit-flipping during quantized transmission.

Figure 1.

Probabilistic quantization and bit-flipping schematic diagram.

The main contributions are summarized as follows.

• (1) This paper investigates the state estimation problem of nonlinear systems with binary-encoding-based quantization, bit-flipping, and heavy-tailed noise.

• (2) This paper employs the normalized innovation sequence as an intermediate variable to jointly handle probabilistic quantization errors and bit-flipping errors. Furthermore, the particle filtering algorithm is applied to mitigate the interference caused by heavy-tailed noise.

• (3) An improved particle filtering method is proposed and applied to the state estimation problem of complex nonlinear systems. Numerical simulations are conducted to verify the effectiveness of the proposed algorithm.

The remainder of this article is organized as follows. Section II introduces the necessary definitions and assumptions. Section III focuses on improving the particle filtering algorithm and applying it to estimate the system state. Section IV conducts simulation experiments based on the algorithm developed in Section III to demonstrate its effectiveness. Finally, some conclusions are drawn in Section V.

Notation: $\mathbb{R}$ denotes the real number set. $\mathrm{E}[·]$ denotes the mathematical expectation. X ^T and X⁻¹ represent the transpose and inverse of the matrix X, respectively. $\left|X\right|$ refers to the spectral radius of the matrix X. $\mathbb{P}(·\mid ·)$ denotes the conditional probability. $\mathrm{Cov}(·)$ is the covariance matrix. $\mathrm{diag}\{a_1,a_2,\dots ,a_n\}$ stands for a diagonal matrix.

Problem formulation

Consider the following nonlinear dynamic system,

$$x_{k+1}=f\left(x_k\right)+w_k,$$

(1)

$$y_k=C_k{x}_k+v_k,$$

(2)

where (1) is the state transition equation, $x_k\in {\mathbb{R}}^n$ denotes the system state, $w_k\in {\mathbb{R}}^n$ is a process noise that follows a Gaussian distribution with zero mean and covariance $\sum _w\ge 0$ . $f(·)$ denotes a known nonlinear function where k denotes the discrete-time index. Equation (2) is the measurement equation. $y_k\in {\mathbb{R}}^{\xi }$ is the measurement vector from the sensors. $v_k\in {\mathbb{R}}^{\xi }$ is the measurement noise, which follows a heavy-tailed distribution with zero mean and covariance $\sum _v\ge 0$ (Figure 2). $C_k\in {\mathbb{R}}^{\xi \times r}$ is the known measurement matrix, k represents the moment.

Figure 2.

Heavy-tailed distribution image.

Probabilistic quantization model

During the transmission of sensor data, a probabilistic quantizer is employed to convert the original information into a finite-length binary code ³⁵ :

$$Q(·)\stackrel{\mathrm{\Delta }}{=}{\left[\begin{array}{cccc}{Q}_1(·)\hfill & Q_2(·)\hfill & \cdots \hfill & Q_{\xi }(·)\hfill \end{array}\right]}^T,$$

where $Q_k(·)$ is the quantizer processing the input scalar. For any vector $b_k\in {\mathbb{R}}^{\xi }$ , we have

$$Q\left(b_k\right)=[Q_1\left(b_{k,1}\right)Q_2\left(b_{k,2}\right)\cdots Q_{\xi }\left(b_{k,\xi }\right){]}^T,$$

in which b_k,n denotes the nth component of the vector b _k , with $n=1,2,\cdots ,\xi$ .

To simplify the analysis, assume that b_k,n belongs to the interval $[-W,W]$ , where W is a positive scalar and L is the number of bits output by the encoder. The quantization range is evenly divided into $2^L-1$ subintervals, corresponding to 2 ^L points, that is $\mathcal{U}\stackrel{\mathrm{\Delta }}{=}\{{\tau }_1,{\tau }_2,\dots ,{\tau }_{2^L}\}$ . Denote the length of each subinterval by $\mathrm{\Delta }={\tau }_{i+1}-{\tau }_i$ , where $i=1,2,\cdots ,2^L-1$ . It follows that $\mathrm{\Delta }=2W/(2^L-1)$ .

When ${\tau }_i\le b_{k,n}\le {\tau }_{i+1}$ , the value b_k,n is quantized according to the following probabilistic rule

$$\begin{array}{c}\{\begin{array}{c}P\{Q_n\left(b_{k,n}\right)={\tau }_i\}=1-r_{k,n}\hfill \\ \multicolumn{1}{c}{P\{Q_n\left(b_{k,n}\right)={\tau }_{i+1}\}=r_{k,n},}\end{array}\hfill \end{array}$$

where $n=1,2,\cdots ,\xi$ , $r_{k,n}\stackrel{\mathrm{\Delta }}{=}(b_{k,n}-{\tau }_i)/\mathrm{\Delta }$ , and $r_{k,n}\in [0,1]$ . Based on the binary encoding scheme, the quantized output is expressed as

$$Q_n(b_{k,n}=-W+\sum _{l=1}^L{h}_{k,n}^{\left(l\right)}{2}^{l-1}\mathrm{\Delta },n=1,2,\dots ,\xi .$$

Therefore, $Q_n\left(b_{k,n}\right)$ can be encoded into the following bit string:

$$D_{k,n}\stackrel{\mathrm{\Delta }}{=}\{h_{k,n}^{\left(1\right)},h_{k,n}^{\left(2\right)},\dots ,h_{k,n}^{\left(L\right)}\}$$

where $h_{k,n}^{\left(l\right)}\in \{0,1\}$ , $n=1,2,\cdots ,\xi$ .

Bit-Flipping model

During the transmission of the bit string D_k,n over a memoryless Binary Symmetric Channel (BSC), bit-flipping may occur due to channel noise and other complex factors. Thus, the final received bit string can be expressed as ³⁵ :

$$\begin{array}{c}{D}_{k,n}^{\prime }\stackrel{\mathrm{\Delta }}{=}\{{\tilde{h}}_{k,n}^{\left(1\right)},{\tilde{h}}_{k,n}^{\left(2\right)},\dots ,{\tilde{h}}_{k,n}^{\left(L\right)}\},\hfill \\ \multicolumn{1}{c}{{\tilde{h}}_{k,n}^{\left(l\right)}={\theta }_{j\to i}^{\left(l\right)}(1-h_{k,n}^{\left(l\right)})+(1-{\theta }_{j\to i}^{\left(l\right)})h_{k,n}^{\left(l\right)}}\end{array}$$

where ${\tilde{h}}_{k,n}^{\left(l\right)}\in \{0,1\}$ , $n=1,2,\cdots ,\xi$ . When bit-flipping occurs, ${\theta }_{j\to i}^{\left(l\right)}=1$ , which means that ${\tilde{h}}_{k,n}^{\left(l\right)}=1-h_{k,n}^{\left(l\right)}$ ; when no flipping occurs, ${\theta }_{j\to i}^{\left(l\right)}=0$ , which implies ${\tilde{h}}_{k,n}^{\left(l\right)}=h_{k,n}^{\left(l\right)}$ . Additionally, define the total number of flipped bits as ${\theta }_{\mathrm{ij}}=\sum _{l=1}^L{\theta }_{j\to i}^{\left(l\right)}$ , and let p denote the probability of bit-flipping. The signal can be decoded using the bit string ${D\prime }_{k,n}$ as follows

$$m_{k,n}=-W+\sum _{l=1}^L{\tilde{h}}_{k,n}^{\left(l\right)}{2}^{l-1}\mathrm{\Delta },n=1,2,\cdots ,\xi .$$

The signals received over the time interval $[0,k]$ can be represented as

$${\mathcal{M}}_k\stackrel{\mathrm{\Delta }}{=}\{m_1,m_2,\dots ,m_k\},m_k=[m_{k,1},m_{k,2},\dots ,m_{k,\xi }{]}^T.$$

Definitions and assumptions

• Definition 1. Define the normalized innovation sequence ³⁵

$${\tilde{y}}_k\stackrel{\mathrm{\Delta }}{=}{S}_{k|k-1}^{-0.5}(y_k-{\hat{y}}_{k|k-1}).$$

The normalized innovation ${\tilde{y}}_k$ is used as the input to a probabilistic quantizer and is further transmitted over a memoryless binary symmetric channel (BES) to the remote estimator. ${\hat{y}}_{k|k-1}\stackrel{\mathrm{\Delta }}{=}\mathbb{E}[y_k\mid {\mathcal{M}}_{k-1}]$ represents the predicted value of the measurement y _k based on prior information, and $S_{k|k-1}\stackrel{\mathrm{\Delta }}{=}\mathrm{Cov}(y_k\mid {\mathcal{M}}_{k-1})$ is the corresponding predicted covariance. In the sequel, ${\tilde{y}}_k$ will be treated as the input to the quantization process.

• Assumption 1. The process noise and the measurement noise are mutually independent, which implies that ³⁵

$$E\left\{\left[\begin{array}{c}{w}_k\hfill \\ \multicolumn{1}{c}{v_k}\end{array}\right]\left[\begin{array}{cc}{w}_k^T\hfill & v_k^T\hfill \end{array}\right]\right\}=\left[\begin{array}{cc}{\Sigma }_w\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\Sigma }_v\hfill \end{array}\right].$$

• Assumption 2. The indicators ${\theta }_{j\to i}^{\left(l\right)}$ are white noise and mutually independent. ³⁵

• Remark: This paper makes substantial improvements to the particle filtering algorithm. Unlike, ⁴⁰ where the measurements are deterministic and not affected by communication processing. The observations in this paper are first probabilistically quantized, then transmitted over a memoryless binary symmetric channel, and further affected by random bit-flipping during transmission. Probabilistic quantization discretizes the measurements and significantly reduces the amount of data to be transmitted, while bit-flipping inevitably introduces transmission errors. Therefore, compared with, ⁴⁰ this paper not only accounts for more realistic observation and communication effects, but also provides systematic enhancements to the design of the proposal distribution and the weight-update mechanism in particle filtering.

Research objectives of this paper

The normalized innovation ${\tilde{y}}_k$ is subject to probabilistic quantization, from which the quantized value $Q\left({\tilde{y}}_k\right)$ is obtained by the estimator. During this process, unavoidable quantization errors are introduced. Moreover, when $Q\left({\tilde{y}}_k\right)$ is transmitted over a BSC, additional bit-flipping may occur due to channel noise, interference, and other complex disturbances, ultimately resulting in a received signal denoted by m _k , and causing signal distortion.

Quantization errors, bit-flip errors, and heavy-tailed noise adversely affect the accuracy and reliability of signal transmission. As a result, traditional state estimation algorithms often fail to achieve satisfactory performance in such nonlinear systems. To overcome these challenges, this paper develops improved particle filtering algorithms to enhance estimation accuracy and system stability.

State estimation based on particle filter

Design of the improved particle filter algorithm

Assume that the particle set at time k−1 is $\{x_{k-1}^i,w_{k-1}^i{\}}_{i=1}^{N_p}$ , where N _p is the total number of particles, and $w_{k-1}^i$ represents the associated weight of particle $x_{k-1}^i$ .

Based on the state transition model (1), the predicted particle ${\hat{x}}_{k|k-1}^i$ at time k can be obtained

$${\hat{x}}_{k|k-1}^i=f\left(x_{k-1}^i\right).$$

Accordingly, the predicted mean ${\hat{x}}_{k|k-1}$ and covariance $P_{k|k-1}$ can be computed

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{x}}_{k|k-1}=\sum _{i=1}^{N_p}{w}_{k-1}^i{\hat{x}}_{k|k-1}^i,}\\ \multicolumn{1}{c}{P_{k|k-1}=\sum _{i=1}^{N_p}{w}_{k-1}^i({\hat{x}}_{k|k-1}^i-{\hat{x}}_{k|k-1}){({\hat{x}}_{k|k-1}^i-{\hat{x}}_{k|k-1})}^T.}\end{array}$$

Consequently, the predicted observation ${\hat{y}}_{k|k-1}$ and its covariance $S_{k|k-1}$ can be derived as follows:

$$\begin{array}{c}{\hat{y}}_{k|k-1}=E\{y_k\mid {\mathcal{M}}_{k-1}\}=E\{C_k{x}_k+v_k\mid {\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{=C_k{\hat{x}}_{k|k-1},}\\ \multicolumn{1}{c}{S_{k|k-1}=\mathrm{Cov}\{y_k\mid {\mathcal{M}}_{k-1}\}=\mathrm{Cov}\{C_k{x}_k+v_k\mid {\mathcal{M}}_{k-1}\}}\\ \multicolumn{1}{c}{=C_k{P}_{k|k-1}{C}_k^T+{\Sigma }_v.}\end{array}$$

Furthermore, ${\tilde{\mathrm{y}}}_k$ can be calculated according to the definition of the normalized innovation sequence ${\tilde{y}}_k\stackrel{\mathrm{\Delta }}{=}{S}_{k|k-1}^{-0.5}(y_k-{\hat{y}}_{k|k-1})$ , where ${\tilde{y}}_k=[{\tilde{y}}_{k,1},{\tilde{y}}_{k,2},\dots ,{\tilde{y}}_{k,\xi }{]}^T$ is a ξ-dimensional vector that follows a standard Gaussian distribution.

Taking ${\tilde{y}}_k$ as the quantization input, according to the quantization rule, we can obtain

$$Q_n\left({\tilde{y}}_{k,n}\right)={\tau }_j=-W+\sum _{l=1}^L{h}_{k,n}^{\left(l\right)}{2}^{l-1}\mathrm{\Delta },$$

in which $n=1,2,\cdots ,\xi$

However, during the transmission over a memoryless BSC, τ _j may undergo bit-flipping. Therefore, the final received signal is

$$m_{k,n}={\tau }_i=-W+\sum _{l=1}^L{\tilde{h}}_{k,n}^{\left(l\right)}{2}^{l-1}\mathrm{\Delta }.$$

The vector $m_k\stackrel{\mathrm{\Delta }}{=}[m_{k,1},m_{k,2},\cdots ,m_{k,\xi }{]}^T$ is formed by m_k,n, where $n=1,2,\cdots ,\xi$ . Accordingly, the set of received vectors is given by ${\mathcal{M}}_k=\{m_1,m_2,\cdots ,m_k\}$ .

Based on Bayes’ rule, ⁴³ the Monte Carlo method is used to approximate the prior probability density function

$$\begin{array}{l}p\left(x_k|{\mathcal{M}}_{k-1}\right)=\int p\left(x_k|x_{k-1}\right)p\left(x_{k-1}|{\mathcal{M}}_{k-1}\right)d{x}_k\\ \approx \sum _{i=1}^{N_p}{w}_{k-1}^{i}p\left(x_k|x_{k-1}^i\right).\end{array}$$

Combined with Bayes’ formula, the posterior probability density function can be obtained

$$p(x_k|{\mathcal{M}}_k)=B_{k}p(m_k|x_k)p(x_k|{\mathcal{M}}_{k-1}),$$

(3)

where $p(x_k|{\mathcal{M}}_k)$ is the posterior probability density, B _k is the normalization constant, and $p(m_k|x_k)$ is the likelihood function.

Subsequently, according to the Monte Carlo approximation, the posterior probability density $p(x_k\mid {\mathcal{M}}_k)$ can be approximated as follows

$$\begin{array}{c}p(x_k|{\mathcal{M}}_k)=B_{k}p(m_k|x_k)p(x_k|{\mathcal{M}}_{k-1})\hfill \\ \multicolumn{1}{c}{\approx \sum _{i=1}^{N_p}{w}_{k-1}^i\frac{p(m_k|x_k)p(x_k|x_{k-1}^i)}{q(x_k|x_{k-1},{\mathcal{M}}_k)}q(x_k|x_{k-1},{\mathcal{M}}_k)}\\ \multicolumn{1}{c}{\approx \sum _{i=1}^{N_p}{w}_k^i\delta (x_k-x_k^i),}\end{array}$$

(4)

where $q(x_k\mid x_{k-1},{\mathcal{M}}_k)$ denotes the proposal distribution, $\delta (·)$ is the Dirac delta function. $w_k^i$ represents the weight associated with the particle $x_k^i$ at time k.

In particle filtering, particles $x_k^i$ can be drawn from the proposal distribution $q(x_k|x_{k-1},{\mathcal{M}}_k)$ , which is described in more detail in the next section. Similar to the idea of the Gaussian particle filter, ⁴³ this paper will adopt the Gaussian approximation of the posterior probability density function as the proposal distribution. After particles are drawn from the proposal distribution $q(x_k|x_{k-1},{\mathcal{M}}_k)$ , their weights are updated according to the method described in subsequent sections. In summary, the method proposed in this paper is shown in Table 1.

Table 1.

Algorithm flow table.

Algorithm: midrule
Step 1 Initialization
Step 2 Calculation proposal distribution $q(x_k\mid x_{k-1},{\mathcal{M}}_k)$ .
Step 3 Sampling from the proposal distribution $x_k^i~q(x_k\mid x_{k-1},{\mathcal{M}}_k)$ .
Step 4 Calculate the particle weight $w_k^i$
Step 5 ${\hat{x}}_k=\sum _{i=1}^{N_p}{w}_k^i{x}_k^i$ .
Step 6 $P_k=\sum _{i=1}^{N_p}{w}_k^i({\hat{x}}_k-x_k^i){({\hat{x}}_k-x_k^i)}^T$ .

Proposal Distribution Design

The proposal distribution is a probability distribution used to generate new particles during the iterative process of particle filtering.

In this work, we adopt a Gaussian approximation of the posterior probability density function as the proposal distribution

$$q(x_k|x_{k-1},{\mathcal{M}}_k)=p_G(x_k|{\mathcal{M}}_k).$$

The iteration of the probability density is simplified to that of the mean and covariance. In addition, the joint probability density of the state x _k and the measurement ${\mathcal{M}}_k$ is assumed to follow a Gaussian distribution. Under these assumptions, the Gaussian approximation of the posterior probability density can be expressed as

$$q(x_k|x_{k-1},{\mathcal{M}}_k)=p_G(x_k|{\mathcal{M}}_k)=\mathcal{N}(x_k;{\hat{x}}_{k|k},P_{k|k})$$

where

$$\begin{array}{c}{\hat{x}}_{k|k}={\hat{x}}_{k|k-1}+K_k{g}_k,\hfill \\ P_{k|k}=P_{k|k-1}-K_k{G}_k{K}_k^T,\hfill \end{array}$$

(5)

$$K_k=P_{k|k-1}{C}_k^T{S}_{k|k-1}^{-0.5},$$

(6)

specifically, the innovation vector and its covariance are defined as

$$\begin{array}{c}{\hat{y}}_{k|k-1}=E\{y_k\mid {\mathcal{M}}_{k-1}\}=E\{C_k{x}_k+v_k\mid {\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{=C_k{\hat{x}}_{k|k-1},}\\ \multicolumn{1}{c}{S_{k|k-1}=\mathrm{Cov}\{y_k\mid {\mathcal{M}}_{k-1}\}=\mathrm{Cov}\{C_k{x}_k+v_k\mid {\mathcal{M}}_{k-1}\}}\\ \multicolumn{1}{c}{=C_k{P}_{k|k-1}{C}_k^T+{\Sigma }_v.}\end{array}$$

Additionally, the coefficients α _i and β _i are computed as:

$$\begin{array}{c}\begin{array}{c}{g}_k={[{\alpha }_{m_{k,1}},{\alpha }_{m_{k,2}},\cdots ,{\alpha }_{m_{k,\xi }}]}^T,\hfill \\ \multicolumn{1}{c}{G_k=\mathrm{diag}\{{\beta }_{m_{k,1}},{\beta }_{m_{k,2}},\cdots ,{\beta }_{m_{k,\xi }}\},}\end{array}\hfill \end{array}$$

furthermore, we have:

$${\alpha }_i=\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}}\frac{F\left({\tau }_j\right)}{H\left({\tau }_j\right)},$$

$${\beta }_i={\alpha }_i^2-\sum _{i=1}^{2^L}{p}^{{\theta }_{\mathrm{ji}}}{(1-p)}^{L-{\theta }_{\mathrm{ji}}}\frac{G\left(\tau \right)}{H\left(\tau \right)}$$

with

$$\begin{array}{c}H\left(\tau \right)=1_{\{\tau \ne W\}}[h(\tau +\mathrm{\Delta })-h\left(\tau \right)+f\left(\tau \right)\mathrm{\Delta }]\hfill \\ \multicolumn{1}{c}{+1_{\{\tau \ne -W\}}[h(\tau -\mathrm{\Delta })-h\left(\tau \right)-f\left(\tau \right)\mathrm{\Delta }],}\\ \multicolumn{1}{c}{F\left(\tau \right)=1_{\{\tau \ne W\}}[f(\tau +\mathrm{\Delta })-f\left(\tau \right)+g\left(\tau \right)\mathrm{\Delta }]}\\ \multicolumn{1}{c}{+1_{\{\tau \ne -W\}}[f(\tau -\mathrm{\Delta })-f\left(\tau \right)-g\left(\tau \right)\mathrm{\Delta }],}\\ \multicolumn{1}{c}{G\left(\tau \right)=1_{\{\tau \ne W\}}[g(\tau +\mathrm{\Delta })-(1-\mathrm{\Delta })g\left(\tau \right)]}\\ \multicolumn{1}{c}{+1_{\{\tau \ne -W\}}[g(\tau -\mathrm{\Delta })-(1+\mathrm{\Delta }\tau )g\left(\tau \right)].}\end{array}$$

The design of the proposal distribution mainly consists of two parts: designing its mean and designing its covariance. Among these, the derivation process for calculating the mean of the proposal distribution is presented in the section “Solving for the Mean in the Proposal Distribution”, while the derivation of its covariance is presented in Section “Solving for the Covariance in the Proposal Distribution”.

• Lemma 1 ³⁵ If there exist d mutually independent events ${\mathrm{\Gamma }}_1,\dots ,{\mathrm{\Gamma }}_d$ , then the following conclusion holds

$$\mathbb{P}\{\mathrm{\Gamma }\mid \underset{i=1}{\stackrel{d}{\cup }}{\mathrm{\Gamma }}_i\}=\sum _{i=1}^d\frac{\mathbb{P}\{{\mathrm{\Gamma }}_i\}}{\sum _{j=1}^d\mathbb{P}\{{\mathrm{\Gamma }}_j\}}\mathbb{P}\{\mathrm{\Gamma }\mid {\mathrm{\Gamma }}_i\}.$$

Proof. Observe that the following identity holds

$$\mathbb{P}\{\mathrm{\Gamma }\mid \underset{i=1}{\stackrel{d}{\cup }}{\mathrm{\Gamma }}_i\}=\frac{\mathbb{P}\{\bigcup _{i=1}^d{\mathrm{\Gamma }}_i\mid \mathrm{\Gamma }\}\mathbb{P}\{\mathrm{\Gamma }\}}{\mathbb{P}\left\{\bigcup _{i=1}^d{\mathrm{\Gamma }}_i\right\}}.$$

Since the events $\{{\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2,\dots ,{\mathrm{\Gamma }}_d\}$ are mutually independent, we have $\mathbb{P}\left\{{\cup }_{i=1}^d{\mathrm{\Gamma }}_i\right\}=\sum _{i=1}^d\mathbb{P}\{{\mathrm{\Gamma }}_i\}$ and $\mathbb{P}\{{\cup }_{i=1}^d{\mathrm{\Gamma }}_i\mid \mathrm{\Gamma }\}=\sum _{i=1}^d\mathbb{P}\{{\mathrm{\Gamma }}_i\mid \mathrm{\Gamma }\}$ . Therefore,

$$\begin{array}{c}\begin{array}{c}\mathbb{P}\{\mathrm{\Gamma }\mid \underset{i=1}{\stackrel{d}{\cup }}{\mathrm{\Gamma }}_i\}=\sum _{i=1}^d\frac{\mathbb{P}\{{\mathrm{\Gamma }}_i\mid \mathrm{\Gamma }\}\mathbb{P}\{\mathrm{\Gamma }\}}{\sum _{j=1}^d\mathbb{P}\{{\mathrm{\Gamma }}_j\}}\hfill \\ \multicolumn{1}{c}{=\sum _{i=1}^d\frac{\mathbb{P}\{{\mathrm{\Gamma }}_i\}}{\sum _{j=1}^d\mathbb{P}\{{\mathrm{\Gamma }}_j\}}·\frac{\mathbb{P}\{{\mathrm{\Gamma }}_i\mid \mathrm{\Gamma }\}\mathbb{P}\{\mathrm{\Gamma }\}}{\mathbb{P}\{{\mathrm{\Gamma }}_i\}}}\\ \multicolumn{1}{c}{=\sum _{i=1}^d\frac{\mathbb{P}\{{\mathrm{\Gamma }}_i\}}{\sum _{j=1}^d\mathbb{P}\{{\mathrm{\Gamma }}_j\}}\mathbb{P}\{\mathrm{\Gamma }\mid {\mathrm{\Gamma }}_i\}.}\end{array}\hfill \end{array}$$

The proposal distribution is designed by determining its mean and covariance. The derivations for both are provided in the subsequent sections. Prior to the derivations, two lemmas are presented in this paper.

• Lemma 2 ³⁵ : If there exist σ- algebras ${\mathcal{C}}_1$ and ${\mathcal{C}}_2$ such that ${\mathcal{C}}_1\subset {\mathcal{C}}_2$ and $E\{\left|\chi \right|\}<\infty$ , then

$$E\{\mathrm{\Gamma }\mid {\mathcal{C}}_1\}=E\{E\{\mathrm{\Gamma }\mid {\mathcal{C}}_2\}\mid {\mathcal{C}}_1\}=E\{E\{\mathrm{\Gamma }\mid {\mathcal{C}}_1\}\mid {\mathcal{C}}_2\}.$$

Solving for the mean in the proposal distribution

We start to calculate the mean ${\hat{x}}_{k|k}=E\{x_k|{\mathcal{M}}_k\}$ in the proposal distribution. According to Lemma 2, it is not difficult to obtain

$$\begin{array}{c}{\hat{x}}_{k|k}=E\{x_k|{\mathcal{M}}_k\}=E\{x_k|{\mathcal{M}}_{k-1},m_k\}\hfill \\ \multicolumn{1}{c}{=E\{E\{x_k|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}|{\mathcal{M}}_{k-1},m_k\}.}\end{array}$$

(7)

Since the normalized innovation sequence ${\tilde{y}}_k$ follows a standard Gaussian distribution, the conditional distribution $p({\tilde{y}}_k\mid {\mathcal{M}}_{k-1})$ satisfies

$$E\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}=0,$$

$$\mathrm{Cov}\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}=I.$$

In addition, the joint distribution of the state x _k and the normalized innovation sequence ${\tilde{y}}_k$ is as follows

$$\begin{array}{c}\begin{array}{c}\mathrm{Cov}\{x_k,{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{=E\{(x_k-{\hat{x}}_{k|k-1}){[S_{k|k-1}^{-0.5}(y_k-{\hat{y}}_{k|k-1})]}^T|{\mathcal{M}}_{k-1}\}}\\ \multicolumn{1}{c}{=E\{(x_k-{\hat{x}}_{k|k-1}){(x_k-{\hat{x}}_{k|k-1})}^T|{\mathcal{M}}_{k-1}\}{C}_k^T{S}_{k|k-1}^{-0.5}}\\ \multicolumn{1}{c}{+E\{(x_k-{\hat{x}}_{k|k-1})v_k^T|{\mathcal{M}}_{k-1}\}{S}_{k|k-1}^{-0.5}}\\ \multicolumn{1}{c}{=P_{k|k-1}{C}_k^T{S}_{k|k-1}^{-0.5}.}\end{array}\hfill \end{array}$$

Since the state x _k and the measurement noise v _k are independent, it follows that $E\{(x_k-{\hat{x}}_{k|k-1})v_k^T|{\mathcal{M}}_{k-1}\}=0$ .

The covariance matrix $\mathrm{Cov}\{{\tilde{y}}_k,x_k|{\mathcal{M}}_{k-1}\}$ is symmetric. Therefore $\mathrm{Cov}\{{\tilde{y}}_k,x_k|{\mathcal{M}}_{k-1}\}=S_{k|k-1}^{-0.5}{C}_k{P}_{k|k-1}^T$ .

According to the formula $K_k=\mathrm{Cov}\{x_k,{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}/\mathrm{Cov}\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}$ , we can obtain

$$K_k=P_{k|k-1}{C}_k^T{S}_{k|k-1}^{-0.5}.$$

Therefore, the conditional expectation $E\{x_k|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}$ can be calculated as follows

$$\begin{array}{c}{\hat{x}}_k^o\stackrel{\mathrm{\Delta }}{=}E\{x_k|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}\hfill \\ \multicolumn{1}{c}{=E\{x_k|{\mathcal{M}}_{k-1}\}+\mathrm{Cov}\{x_k,{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}}\\ \multicolumn{1}{c}{\times \mathrm{Cov}{\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}}^{-1}({\tilde{y}}_k-E\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\})}\\ \multicolumn{1}{c}{={\hat{x}}_{k|k-1}+K_k{\tilde{y}}_k.}\end{array}$$

(8)

Combined with equation (7), by evaluating the expectation $E\{{\hat{x}}_k^o|{\mathcal{M}}_{k-1},m_k\}$ , we can obtain ${\hat{x}}_{k|k}$

$$\begin{array}{c}\begin{array}{c}{\hat{x}}_{k|k}=E\{{\hat{x}}_k^o\mid {\mathcal{M}}_{k-1},m_k\}\hfill \\ \multicolumn{1}{c}{={\hat{x}}_{k|k-1}+K_{k}E\{{\tilde{y}}_k\mid {\mathcal{M}}_{k-1},m_k\}.}\end{array}\hfill \end{array}$$

(9)

Then, we only need to calculate $E\{{\tilde{y}}_k|{\mathcal{M}}_{k-1},m_k\}$

Note that the conditional expectation $E\{{\tilde{y}}_k|{\mathcal{M}}_{k-1},m_k\}$ can be decomposed into components $E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},m_{k,n}\}$ , where $n=1,2,\cdots ,\xi$ . To compute each component, we define the event ${\mathcal{F}}_j^i$ , which denotes the event that the channel input $Q_n\left({\tilde{y}}_{k,n}\right)={\tau }_j$ (abbreviated as $Q_{k,n}={\tau }_j$ in the following text) is received as $m_{k,n}={\tau }_i$ due to bit-flipping during the transmission process. In practice, each received symbol $m_{k,n}={\tau }_i(n=1,2,\cdots ,\xi )$ can result from one of 2 ^L bit-flip patterns, corresponding to the following event set:

$$\{{\mathcal{F}}_1^i,{\mathcal{F}}_2^i,\cdots ,{\mathcal{F}}_{2^L}^i\}.$$

According to Lemma 1, we derive the expectation

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},m_{k,n}={\tau }_i\}\hfill \\ \multicolumn{1}{c}{=\sum _{j=1}^{2^L}\frac{\mathbb{P}\{{\mathcal{F}}_j^i|{\mathcal{M}}_{k-1}\}}{\sum _{s=1}^{2^L}\mathbb{P}\{{\mathcal{F}}_s^i|{\mathcal{M}}_{k-1}\}}}\\ \multicolumn{1}{c}{\times E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\},}\end{array}\hfill \end{array}$$

(10)

where $n=1,2,\dots ,\xi$

Since ${\mathcal{F}}_j^i$ is independent of ${\mathcal{M}}_{k-1}$ , we have $\mathbb{P}\{{\mathcal{F}}_j^i|{\mathcal{M}}_{k-1}\}=\mathbb{P}\{{\mathcal{F}}_j^i\}$ .

In addition, since ${\mathcal{F}}_j^i$ can be represented by a bit string of length L, and the bit-flipping are mutually independent, the flipping of each bit can be considered a Bernoulli trial. Therefore, ${\mathcal{F}}_j^i$ can be viewed as an L Bernoulli trial. According to the bit-flipping model introduced in Section 2, if the l-th bit of τ _j flips, then ${\theta }_{j\to i}^{\left(l\right)}=1$ ; otherwise, ${\theta }_{j\to i}^{\left(l\right)}=0$ . Therefore, the probability of the event ${\mathcal{F}}_j^i$ is given by

$$\mathbb{P}\{{\mathcal{F}}_j^i\}=p^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}},$$

(11)

where ${\theta }_{\mathrm{ij}}=\sum _{l=1}^L{\theta }_{j\to i}^{\left(l\right)}$ . Obviously, the total probability over all L-fold Bernoulli trials equals 1, which means $\sum _{l=1}^{2^L}\mathbb{P}\{{\mathcal{F}}_l^i\}=1$ . By substituting (11) into (10), we can obtain

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}\mid {\mathcal{M}}_{k-1},m_{k,n}={\tau }_i\}\hfill \\ \multicolumn{1}{c}{=\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}}}\\ \multicolumn{1}{c}{\times E\{{\tilde{y}}_{k,n}\mid {\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\},}\end{array}\hfill \end{array}$$

(12)

where $n=1,2,\cdots ,\xi$ . Moreover, the information set $\{{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\}$ can be equivalently characterized as the union of the two sets $\{{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\}$ and $\{{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}$ , for $j=2,\cdots ,2^L-1$ .

${\mathcal{U}}_j$ denotes the interval $[{\tau }_j,{\tau }_{j+1}]$ ; ${\mathcal{L}}_k^j$ indicates that ${\tilde{y}}_{k,n}\in {\mathcal{U}}_j$ is quantized to τ _j (i.e. the left endpoint of ${\mathcal{U}}_j$ ); ${\mathcal{R}}_k^j$ indicates that ${\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1}$ is quantized to τ _j (i.e. the right endpoint of ${\mathcal{U}}_{j-1}$ ). This holds for all $n=1,2,\cdots ,\xi$ .

According to Lemma 1, we have

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\}\hfill \\ \multicolumn{1}{c}{=\frac{p_{{\mathcal{R}}_k^j}}{p_{{\mathcal{L}}_k^j}+p_{{\mathcal{R}}_k^j}}E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}}\\ \multicolumn{1}{c}{+\frac{p_{{\mathcal{L}}_k^j}}{p_{{\mathcal{L}}_k^j}+p_{{\mathcal{R}}_k^j}}E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\},}\end{array}\hfill \end{array}$$

(13)

where $n=1,2\cdots ,\xi$ , we use $p_{{\mathcal{R}}_k^j}$ as a shorthand for the probability $\mathbb{P}\{{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}$ , and $p_{{\mathcal{L}}_k^j}$ for $\mathbb{P}\{{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\}$ . These expressions are valid for $j=2,\cdots ,2^L-1$ ; $n=1,2,\cdots ,\xi$ .

In the subsequent calculations, we first evaluate the probabilities $p_{{\mathcal{L}}_k^j}$ and $p_{{\mathcal{R}}_k^j}$ . According to Bayes’ theorem, the probability $p_{{\mathcal{L}}_k^i}$ can be factorized into the product of the following three components

$$\begin{array}{c}{p}_{{\mathcal{L}}_k^j}=\mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{\times \mathbb{P}\{{\tilde{y}}_{k,n}\in {\mathcal{U}}_j|{\mathcal{M}}_{k-1}\}\times \mathbb{P}\{{\mathcal{M}}_{k-1}\}.}\end{array}$$

Given that $p({\tilde{y}}_{k,n}\mid {\mathcal{M}}_{k-1})~\mathcal{N}(y_{k,n};0,1)$ , it follows that

$$\mathbb{P}\{{\tilde{\mathrm{y}}}_{k,n}\in {\mathcal{U}}_j|{\mathcal{M}}_{k-1}\}=\phi \left({\tau }_{j+1}\right)-\phi \left({\tau }_j\right).$$

Moreover, the conditional probability of ${\mathcal{L}}_k^j$ can be calculated by the following formula

$$\begin{array}{c}\mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{=\underset{\mathbb{R}}{\int }p({\tilde{y}}_{k,n}|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1})}\\ \multicolumn{1}{c}{\times \mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1}\}d{\tilde{y}}_{k,n}.}\end{array}$$

Given the specific situation of ${\tilde{y}}_{k,n}$ from the perspective of probability calculation ${\tilde{y}}_{k,n}\in {\mathcal{U}}_j$ and ${\mathcal{M}}_{k-1}$ are redundant, so we have

$$\begin{array}{c}\multicolumn{1}{c}{\mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1}\}=\mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n}\}}\\ \multicolumn{1}{c}{=\frac{{\tau }_{j+1}-{\tilde{y}}_{k,n}}{\mathrm{\Delta }}.}\end{array}$$

In addition, we also have

$$\begin{array}{c}\multicolumn{1}{c}{p({\tilde{y}}_{k,n}|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1})=}\\ \multicolumn{1}{c}{\{\begin{array}{c}\frac{\mathcal{N}({\tilde{y}}_{k,n};0,1)}{\phi \left({\tau }_{j+1}\right)-\phi \left({\tau }_j\right)},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j.\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise}.}\end{array}}\end{array}$$

Therefore, we can derive the following result

$$\begin{array}{c}\begin{array}{c}\mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1}\}=\underset{{\mathcal{U}}_j}{\int }\frac{\mathcal{N}({\tilde{y}}_{k,n};0,1)}{\phi \left({\tau }_{j+1}\right)-\phi \left({\tau }_j\right)}\hfill \\ \multicolumn{1}{c}{\times \frac{{\tau }_{j+1}-{\tilde{y}}_{k,n}}{\mathrm{\Delta }}d{\tilde{y}}_{k,n},}\end{array}\hfill \end{array}$$

in which $n=1,2\cdots ,\xi$ .

To simplify the notation, we define ${\overline{p}}_{{\mathcal{L}}_k^j}$ as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{p}}_{{\mathcal{L}}_k^j}\stackrel{\mathrm{\Delta }}{=}\mathbb{P}\{{\tilde{y}}_{k,n}\in {\mathcal{U}}_j|{\mathcal{M}}_{k-1}\}\mathbb{P}\{{\mathcal{L}}_k^j|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{M}}_{k-1}\}}\\ \multicolumn{1}{c}{={\int }_{{\mathcal{U}}_j}\frac{{\tau }_{j+1}-{\tilde{y}}_{k,n}}{\mathrm{\Delta }}\mathcal{N}({\tilde{y}}_{k,n};0,1)d{\tilde{y}}_{k,n}}\\ \multicolumn{1}{c}{=\frac{\exp (-\frac{{\tau }_{j+1}^2}{2})-\exp (-\frac{{\tau }_j^2}{2})}{\sqrt{2\pi }\mathrm{\Delta }}}\\ \multicolumn{1}{c}{+\frac{{\tau }_{j+1}[\phi \left({\tau }_{j+1}\right)-\phi \left({\tau }_j\right)]}{\mathrm{\Delta }}.}\end{array}$$

Accordingly, the probability $p_{{\mathcal{L}}_k^j}$ can be expressed as

$$p_{{\mathcal{L}}_k^j}={\overline{p}}_{{\mathcal{L}}_k^j}\mathbb{P}\{{\mathcal{M}}_{k-1}\}.$$

In a similar manner, we define

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{p}}_{{\mathcal{R}}_k^j}\stackrel{\mathrm{\Delta }}{=}\frac{\exp (-\frac{{\tau }_{j-1}^2}{2})-\exp (-\frac{{\tau }_j^2}{2})}{\sqrt{2\pi }\mathrm{\Delta }}}\\ \multicolumn{1}{c}{+\frac{{\tau }_{j-1}[\phi \left({\tau }_{j-1}\right)-\phi \left({\tau }_j\right)]}{\mathrm{\Delta }},}\end{array}$$

with the corresponding probability $p_{{\mathcal{R}}_k^j}={\overline{p}}_{{\mathcal{R}}_k^j}\mathbb{P}\{{\mathcal{M}}_{k-1}\}$ .

Furthermore, we denote ${\overline{p}}_{{\mathcal{R}}_k^j}+{\overline{p}}_{{\mathcal{L}}_k^j}$ as $H\left({\tau }_j\right)/\mathrm{\Delta }$ , where the function $H\left(\tau \right)$ is given by

$$H\left(\tau \right)=h(\tau -\mathrm{\Delta })+h(\tau +\mathrm{\Delta })-2h\left(\tau \right),$$

in which the auxiliary function $h\left(\tau \right)$ is defined as

$$h\left(\tau \right)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{{\tau }^2}{2})+\tau \phi \left(\tau \right).$$

Next, we calculate the expectations $E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}$ and $E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\}$ . According to Bayes’ formula, we can derive

$$\begin{array}{c}p({\tilde{y}}_{k,n}|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j,{\mathcal{M}}_{k-1})\hfill \\ \multicolumn{1}{c}{=\{\begin{array}{c}\frac{p({\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1})\mathbb{P}\{{\mathcal{L}}_n^i|{\tilde{y}}_{n,k}\}}{\mathbb{P}\{{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j|{\mathcal{M}}_{k-1}\}},\mathrm{if}{\tilde{y}}_{k,n}\in {\mathcal{U}}_j.\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise}.}\end{array}}\\ \multicolumn{1}{c}{=\{\begin{array}{c}\frac{1}{{\overline{p}}_{{\mathcal{L}}_k^j}}\frac{{\tau }_{j+1}-{\tilde{y}}_{k,n}}{\mathrm{\Delta }}\mathcal{N}({\tilde{y}}_{k,n};0,1),\mathrm{if}{\tilde{y}}_{k,n}\in {\mathcal{U}}_j.\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise}.}\end{array}}\end{array}$$

According to the formula for calculating the mean in probability theory, we can obtain

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}|{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j,{\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{=\frac{1}{{\overline{p}}_{{\mathcal{L}}_k^j}\mathrm{\Delta }}\underset{{\mathcal{U}}_j}{\int }{\tilde{y}}_{k,n}({\tau }_{j+1}-{\tilde{y}}_{k,n})\mathcal{N}({\tilde{y}}_{k,n};0,1)d{\tilde{y}}_{k,n}}\\ \multicolumn{1}{c}{=\frac{1}{{\overline{p}}_{{\mathcal{L}}_k^j}\mathrm{\Delta }}[\phi \left({\tau }_j\right)-\phi \left({\tau }_{j+1}\right)+\frac{\mathrm{\Delta }}{\sqrt{2\pi }}\exp (-\frac{{\tau }_j^2}{2})].}\end{array}\hfill \end{array}$$

(14)

Applying the same reasoning, we derive

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}|{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j,{\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{=\frac{1}{{\overline{p}}_{{\mathcal{R}}_k^j}\mathrm{\Delta }}[\phi \left({\tau }_j\right)-\phi \left({\tau }_{j-1}\right)-\frac{\mathrm{\Delta }}{\sqrt{2\pi }}\exp (-\frac{{\tau }_j^2}{2})].}\end{array}\hfill \end{array}$$

(15)

By substituting (14) and (15) into (13), we can simplify the expression as follows

$$E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\}=\frac{F\left({\tau }_j\right)}{H\left({\tau }_j\right)},$$

(16)

where $n=1,2,\cdots ,\xi$ , $F\left(\tau \right)=f(\tau -\mathrm{\Delta })+f(\tau +\mathrm{\Delta })-2f\left(\tau \right)$ , and $f\left(\tau \right)=-\phi \left(\tau \right)$ .

In particular, when j= 1, the expectation $E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\}$ can be simplified to

$$E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\},$$

and when j = 2 ^L , it reduces to

$$E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}.$$

Hence, the function $F\left(\tau \right)$ can be expressed as

$$F\left(\tau \right)=\{\begin{array}{c}\phi \left(\tau \right)-\phi (\tau +\mathrm{\Delta })+\frac{\mathrm{\Delta }}{\sqrt{2\pi }}\exp (-\frac{{\tau }^2}{2}),\tau =-W.\hfill \\ \multicolumn{1}{c}{\phi \left(\tau \right)-\phi (\tau -\mathrm{\Delta })-\frac{\mathrm{\Delta }}{\sqrt{2\pi }}\exp (-\frac{{\tau }^2}{2}),\tau =W.}\end{array}$$

In addition, the function $H\left(\tau \right)$ is given by

$$H\left(\tau \right)=\{\begin{array}{c}\frac{1}{\sqrt{2\pi }}[\exp (-\frac{{(\tau +\mathrm{\Delta })}^2}{2})-\exp (-\frac{{\tau }^2}{2})]\hfill \\ \multicolumn{1}{c}{+(\tau +\mathrm{\Delta })[\phi (\tau +\mathrm{\Delta })-\phi \left(\tau \right)],\tau =-W.}\\ \multicolumn{1}{c}{\frac{1}{\sqrt{2\pi }}[\exp (-\frac{{(\tau -\mathrm{\Delta })}^2}{2})-\exp (-\frac{{\tau }^2}{2})]}\\ \multicolumn{1}{c}{+(\tau -\mathrm{\Delta })[\phi (\tau -\mathrm{\Delta })-\phi \left(\tau \right)],\tau =W.}\end{array}$$

Let ${\alpha }_i\stackrel{\mathrm{\Delta }}{=}E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_{k-1},m_{k,n}={\tau }_i\}$ denote the expected value associated with the received symbol τ _i at time k, where $n=1,2,\cdots ,\xi$ . Based on (12) and (16), we obtain ${\alpha }_i=\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}}\frac{F\left({\tau }_j\right)}{H\left({\tau }_j\right)}$ .

Therefore, (9) can be rewritten as

$${\hat{x}}_{k|k}={\hat{x}}_{k|k-1}+K_k{g}_k,$$

where $g_k=[{\alpha }_{m_{k,1}},{\alpha }_{m_{k,2}},\cdots ,{\alpha }_{m_{k,L}}{]}^T$ .

Solving for the covariance in the proposal distribution

We next compute the covariance $P_{k|k}\stackrel{\mathrm{\Delta }}{=}E\{(x_k-{\hat{x}}_{k|k}){(x_k-{\hat{x}}_{k|k})}^T|{\mathcal{M}}_k\}$ .

According to (8) and (9) the estimation error $e_k=x_k-{\hat{x}}_{k|k}$ can be expressed as

$$\begin{array}{c}\begin{array}{c}{e}_k=x_k-{\hat{x}}_k^o+{\hat{x}}_k^o-{\hat{x}}_{k|k}\hfill \\ \multicolumn{1}{c}{=x_k-{\hat{x}}_k^o+K_k({\tilde{y}}_k-E\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}).}\end{array}\hfill \end{array}$$

(17)

By Lemma 2, it follows that

$$\begin{array}{c}\begin{array}{c}{P}_{k|k}=E\{e_k{e}_k^T|{\mathcal{M}}_k\}\hfill \\ \multicolumn{1}{c}{=E\{E\{e_k{e}_k^T|{\mathcal{M}}_{k-1},{\tilde{\mathrm{y}}}_k\}|{\mathcal{M}}_{k-1},m_k\}.}\end{array}\hfill \end{array}$$

(18)

Since under the condition $E\{(x_k-{\hat{x}}_k^o)|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}=E\{x_k|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}-{\hat{x}}_k^o={\hat{x}}_k^o-{\hat{x}}_k^o=0$ and $\{{\mathcal{M}}_{k-1},{\tilde{y}}_k\}$ , ${\tilde{y}}_k$ and $E\{{\tilde{y}}_k|{\mathcal{M}}_k\}$ are deterministic functions, we can substitute (17) into (18) and expand to obtain

$$\begin{array}{c}\begin{array}{c}E\{e_k{e}_k^T|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}\hfill \\ \multicolumn{1}{c}{=E\{(x_k-{\hat{x}}_k^o){(x_k-{\hat{x}}_k^o)}^T|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}}\\ \multicolumn{1}{c}{+K_k({\tilde{y}}_k-E\{{\tilde{\mathrm{y}}}_k|{\mathcal{M}}_k\}){({\tilde{y}}_k-E\{{\tilde{y}}_k|{\mathcal{M}}_k\})}^T{K}_k^T.}\end{array}\hfill \end{array}$$

(19)

Meanwhile, since the first term after the equal sign in (19) is the covariance of $p(x_k|{\mathcal{M}}_{k-1},{\hat{y}}_k)$ , it can be obtained by the following formula

$$\begin{array}{c}\begin{array}{c}{P}_k^o\stackrel{\mathrm{\Delta }}{=}\mathrm{Cov}\{x_k|{\mathcal{M}}_{k-1},{\tilde{y}}_k\}=\mathrm{Cov}\{x_k|{\mathcal{M}}_{k-1}\}\hfill \\ \multicolumn{1}{c}{-\mathrm{Cov}\{x_k,{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}\mathrm{Cov}{\{{\tilde{y}}_k|{\mathcal{M}}_{k-1}\}}^{-1}}\\ \multicolumn{1}{c}{\times \mathrm{Cov}\{{\tilde{y}}_k,x_k|{\mathcal{M}}_{k-1}\}=P_{k|k-1}-K_k{K}_k^T.}\end{array}\hfill \end{array}$$

(20)

Combining (18), (19), and (20), we obtain

$$\begin{array}{c}\begin{array}{c}{P}_{k|k}=P_k^o+K_{k}E\{({\tilde{y}}_k-E\{{\tilde{y}}_k|{\mathcal{M}}_k\})\hfill \\ \multicolumn{1}{c}{\times {({\tilde{y}}_k-E\{{\tilde{y}}_k|{\mathcal{M}}_k\})}^T|{\mathcal{M}}_k\}{K}_k^T}\\ \multicolumn{1}{c}{=P_{k|k-1}-K_k{K}_k^T+K_k\mathrm{Cov}\{{\tilde{y}}_k|{\mathcal{M}}_k\}{K}_k^T}\\ \multicolumn{1}{c}{=P_{k|k-1}-K_k(I-\mathrm{Cov}\{{\tilde{y}}_k|{\mathcal{M}}_k\})K_k^T,}\end{array}\hfill \end{array}$$

where

$$\mathrm{Cov}\{{\tilde{y}}_k|{\mathcal{M}}_k\}=E\{{\tilde{y}}_k{\tilde{y}}_k^T|{\mathcal{M}}_k\}-E\{{\tilde{y}}_k|{\mathcal{M}}_k\}E\{{\tilde{y}}_k|{\mathcal{M}}_k{\}}^T.$$

Since the normalized innovation sequence ${\tilde{y}}_{k,n}$ is a sequence of independent and identically distributed random variables following $p\left({\tilde{y}}_{k,n}\right)=\mathcal{N}({\tilde{y}}_{k,n};0,1)$ , the off-diagonal variables of the covariance $\mathrm{Cov}\{{\tilde{y}}_k|{\mathcal{M}}_k\}$ are all 0. For the diagonal components, we have $\mathrm{Var}\{{\tilde{y}}_{k,n}|{\mathcal{M}}_k\}=E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_k\}-{(E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_k\})}^2$ , where the first moment $E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_k\}$ is obtained from (12), with $n=1,2,\cdots ,\xi$ . Note that the second moment $E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_k\}$ can also be calculated similarly.

Next, we calculate the moment $E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_k\}$

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_{k-1},m_{k,n}={\tau }_i\}\hfill \\ \multicolumn{1}{c}{=\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}}}\\ \multicolumn{1}{c}{\times E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\},}\end{array}\hfill \end{array}$$

and for $j=2,\cdots ,2^L-1$ ,we have

$$\begin{array}{c}\begin{array}{c}E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_{k-1},Q_{k,n}={\tau }_j\}\hfill \\ \multicolumn{1}{c}{=\frac{p_{{\mathcal{R}}_k^j}}{p_{{\mathcal{L}}_k^j}+p_{{\mathcal{R}}_k^j}}E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}}\\ \multicolumn{1}{c}{+\frac{p_{{\mathcal{L}}_k^j}}{p_{{\mathcal{L}}_k^j}+p_{{\mathcal{R}}_k^j}}E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_{k-1},{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\}.}\end{array}\hfill \end{array}$$

Finally, similar to the calculation method of the first moment $E\{{\tilde{y}}_{k,n}|{\mathcal{M}}_k\}$ , we can obtain

$$E\{{\tilde{y}}_{k,n}^2|{\mathcal{M}}_{k-1},m_k={\tau }_i\}=\frac{\Lambda \left({\tau }_i\right)}{H\left({\tau }_i\right)},$$

where $n=1,2,\cdots ,\xi$ , $\Lambda \left(\tau \right)=\lambda (\tau -\mathrm{\Delta })+\lambda (\tau +\mathrm{\Delta })-2\lambda \left(\tau \right)$ , and $\lambda \left(\tau \right)=\frac{2}{\sqrt{2\pi }}\exp (-\frac{{\tau }^2}{2})+\tau \phi \left(\tau \right)$ .

Next, we define ${\beta }_i\stackrel{\mathrm{\Delta }}{=}1-\mathrm{Var}\{{\tilde{y}}_{k,n}|{\mathcal{M}}_k\}$ , and we derive

$$\begin{array}{c}{\beta }_i=1-\{\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ji}}}{(1-p)}^{L-{\theta }_{\mathrm{ji}}}\frac{\Lambda \left(\tau \right)}{H\left(\tau \right)}-{\alpha }_i^2\}\hfill \\ \multicolumn{1}{c}{={\alpha }_i^2-\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ji}}}{(1-p)}^{L-{\theta }_{\mathrm{ji}}}\frac{\Lambda \left(\tau \right)-H\left(\tau \right)}{H\left(\tau \right)}}\\ \multicolumn{1}{c}{={\alpha }_i^2-\sum _{i=1}^{2^L}{p}^{{\theta }_{\mathrm{ji}}}{(1-p)}^{L-{\theta }_{\mathrm{ji}}}\frac{G\left(\tau \right)}{H\left(\tau \right)},}\end{array}$$

where $G\left(\tau \right)=g(\tau -\mathrm{\Delta })+g(\tau +\mathrm{\Delta })-2g\left(\tau \right)$ , $g\left(\tau \right)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{{\tau }^2}{2})$ . Particularly, when i= 1 or $i=2^L$ , it follows that

$$G\left(\tau \right)=\{\begin{array}{c}g(\tau +\mathrm{\Delta })-(1-\mathrm{\Delta }\tau )g\left(\tau \right),\tau =-W.\hfill \\ \multicolumn{1}{c}{g(\tau -\mathrm{\Delta })-(1+\mathrm{\Delta }\tau )g\left(\tau \right),\tau =W.}\end{array}$$

Define $G_k\stackrel{\mathrm{\Delta }}{=}\mathrm{diag}\{{\beta }_{k,1},{\beta }_{k,2},\cdots ,{\beta }_{k,\xi }\}$ . Therefore, it is easy to obtain

$$P_{k|k}=P_{k|k-1}-K_k{G}_k{K}_k^T.$$

Update of particle weights

The particle weight can be determined by calculating the ratio of the true posterior probability density to the proposal distribution. From (4), the particle weight update is given by the following formula

$$w_k^i=w_{k-1}^i\frac{p(m_k|x_k)p(x_k|x_{k-1}^i)}{q(x_k|x_{k-1},{\mathcal{M}}_k)}.$$

(21)

Since the proposal distribution $q(x_k|x_{k-1},{\mathcal{M}}_k)$ has been determined in section Solving for the Covariance in the Proposal Distribution, we now focus on computing the likelihood function $p(m_k|x_k)$ .

It follows from probability theory that

$$\begin{array}{c}\begin{array}{c}p(m_k|x_k)=p(\underset{n=1}{\overset{\xi }{\cup }}\{m_{k,n}={\tau }_i\}|x_k)\hfill \\ \multicolumn{1}{c}{=\underset{n=1}{\overset{\xi }{\Pi }}p(m_{k,n}={\tau }_i|x_k),}\end{array}\hfill \end{array}$$

where $n=1,2,\cdots ,\xi$ . According to the relationship between the set $\{m_{k,n}={\tau }_i\}$ and the set ${\cup }_{j=1}^{2^L}\{Q\left({\tilde{y}}_{k,n}\right)={\tau }_j,{\mathcal{F}}_j^i\}$ we can get

$$\begin{array}{c}\begin{array}{c}p(m_{k,n}={\tau }_i|x_k)=p(\underset{j=1}{\overset{2^L}{\cup }}\{Q\left({\tilde{\mathrm{y}}}_{k,n}\right)={\tau }_j,{\mathcal{F}}_j^i\}|x_k)\hfill \\ \multicolumn{1}{c}{=\sum _{j=1}^{2^L}\mathbb{P}\{{\mathcal{F}}_j^i\}\mathbb{P}\{Q\left({\tilde{\mathrm{y}}}_{k,n}\right)={\tau }_j|x_k\},}\end{array}\hfill \end{array}$$

where the probability $\mathbb{P}\{{\mathcal{F}}_j^i\}=p^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}}$ , $n=1,2,\cdots ,\xi$ . Based on Bayes’ theorem, we can obtain

$$\begin{array}{l}\mathbb{P}\left\{Q\left({\tilde{y}}_{k,n}\right)\right|x_k\}={\int }_{\mathbb{R}}p\left({\tilde{y}}_{k,n},Q\left({\tilde{y}}_{k,n}\right)|x_k\right)\mathrm{d}{\tilde{y}}_{k,n}\\ ={\int }_{\mathbb{R}}\mathbb{P}\left\{Q\left({\tilde{y}}_{k,n}\right)\right|{\tilde{y}}_{k,n}\}p\left({\tilde{y}}_{k,n}|x_k\right)\mathrm{d}{\tilde{y}}_{k,n}.\end{array}$$

(22)

It is easy to obtain

$$p({\tilde{y}}_{k,n}|x_k)=\mathcal{N}({\tilde{y}}_{k,n};{\mu }_{k,n},{\sigma }_{k,n}^2),n=1,2,\cdots ,\xi$$

where

$${\mu }_{k,n}=\frac{C_k(x_k-{\hat{x}}_{k|k-1})}{\sqrt{C_k{P}_{k|k-1}{C}_k^T+{\Sigma }_v^2}},$$

$${\sigma }_{k,n}=\frac{{\Sigma }_v}{\sqrt{C_k{P}_{k|k-1}{C}_k^T+{\Sigma }_v^2}}.$$

In addition, since the set $\{Q\left({\tilde{y}}_{k,n}\right)={\tau }_j\}$ is the union of the sets $\{{\tilde{y}}_{k,n}\in {\mathcal{U}}_j,{\mathcal{L}}_k^j\}$ and $\{{\tilde{y}}_{k,n}\in {\mathcal{U}}_{j-1},{\mathcal{R}}_k^j\}$ , (22) can be written as

$$\begin{array}{c}\multicolumn{1}{c}{\begin{array}{c}\multicolumn{1}{c}{\mathbb{P}\{Q\left({\tilde{y}}_{k,n}\right)={\tau }_j|x_k\}={\int }_{{\mathcal{U}}_{j-1}}\frac{{\tilde{y}}_{k,n}-{\tau }_{j-1}}{\mathrm{\Delta }}}\\ \multicolumn{1}{c}{\times \mathcal{N}({\tilde{y}}_{k,n};{\mu }_{k,n},{\sigma }_{k,n}^2)d{\tilde{y}}_{k,n}}\\ \multicolumn{1}{c}{+{\int }_{{\mathcal{U}}_j}\frac{{\tau }_{j+1}-{\tilde{y}}_{k,n}}{\mathrm{\Delta }}}\\ \multicolumn{1}{c}{\times \mathcal{N}({\tilde{y}}_{k,n};{\mu }_{k,n},{\sigma }_{k,n}^2)d{\tilde{y}}_{k,n},}\end{array}}\end{array}$$

(23)

in which $j=2,\cdots ,2^L-1$ , $n=1,2,\cdots ,\xi$ . We define the function ${\psi }_{k,n}(\tau ;x_k)\stackrel{\mathrm{\Delta }}{=}\psi \left(\frac{\tau -{\mu }_{k,n}}{{\sigma }_{k,n}}\right)$ , where the function $\psi \left(\tau \right)={\left(2\pi \right)}^{-0.5}\times \exp (-\frac{{\tau }^2}{2})+\tau \phi \left(\tau \right)$ and the function ${\mathrm{\Psi }}_{k,n}(\tau ;x_k)\stackrel{\mathrm{\Delta }}{=}{\psi }_{k,n}(\tau -\mathrm{\Delta };x_k)+{\psi }_{k,n}(\tau +\mathrm{\Delta };x_k)-2\times {\psi }_{k,n}(\tau ;x_k)$ are also defined. From equation (23), we can obtain

$$\mathbb{P}\{Q\left({\tilde{y}}_{k,n}\right)={\tau }_j|x_k\}=\frac{{\sigma }_{k,n}^2}{\mathrm{\Delta }}{\mathrm{\Psi }}_{k,n}({\tau }_j;x_k),n=1,2,\cdots ,\xi .$$

Particularly, when j= 1, the expression can be simplified to the first integral term in (23). When $j=2^L$ , $\mathbb{P}\{Q\left({\tilde{y}}_{k,n}\right)={\tau }_j|x_k\}$ can be simplified to the second integral term in (23). Therefore, we have

$${\mathrm{\Psi }}_{k,n}(\tau ;x_k)=\{\begin{array}{c}{\psi }_{k,n}(\tau +\mathrm{\Delta };x_k)-{\psi }_{k,n}(\tau ;x_k)\hfill \\ \multicolumn{1}{c}{-\frac{\mathrm{\Delta }}{{\sigma }_{k,n}}\phi \left(\frac{\tau -{\mu }_{k,n}}{{\sigma }_{k,n}}\right),\tau =-W.}\\ \multicolumn{1}{c}{{\psi }_{k,n}(\tau -\mathrm{\Delta };x_k)-{\psi }_{k,n}(\tau ;x_k)}\\ \multicolumn{1}{c}{+\frac{\mathrm{\Delta }}{{\sigma }_{k,n}}\phi \left(\frac{\tau -{\mu }_{k,n}}{{\sigma }_{k,n}}\right),\tau =W.}\\ \multicolumn{1}{c}{{\psi }_{k,n}(\tau -\mathrm{\Delta };x_k)+{\psi }_{k,n}(\tau +\mathrm{\Delta };x_k)}\\ \multicolumn{1}{c}{-2{\psi }_{k,n}(\tau ;x_k),\mathrm{otherwise}.}\end{array}$$

If we define ${\Upsilon }_{k,n}\stackrel{\mathrm{\Delta }}{=}\sum _{j=1}^{2^L}{p}^{{\theta }_{\mathrm{ij}}}{(1-p)}^{L-{\theta }_{\mathrm{ij}}}{\mathrm{\Psi }}_{k,n}({\tau }_j;x_k)$ , then we have

$$p(m_{k,n}|x_k)=\frac{{\sigma }_{k,n}^2}{\mathrm{\Delta }}{\Upsilon }_{k,n},$$

in which $n=1,2,\cdots ,\xi$

Finally, the corresponding particle weights at time k can be obtained from (21).

According to the particle weight update in (21), a set of unnormalized weights can be calculated, and subsequently normalized such that their sum equals 1.

Therefore, at time k, particles are obtained through the proposal distribution $q(x_k|x_{k-1},{\mathcal{M}}_k)=p_G(x_k|{\mathcal{M}}_k)$ , and the weights of the corresponding particles are calculated according to (21). In this way, the particle set $\{x_k^i,w_k^i{\}}_{i=1}^{N_p}$ at time k is obtained. The mean ${\hat{x}}_k$ and covariance P _k of the posterior probability density can be calculated according to the following two formulas

$$\begin{array}{c}{\hat{x}}_k=\sum _{i=1}^{N_p}{w}_k^i{x}_k^i,\hfill \\ \multicolumn{1}{c}{P_k=\sum _{i=1}^{N_p}{w}_k^i({\hat{x}}_k-x_k^i){({\hat{x}}_k-x_k^i)}^T.}\end{array}$$

Here, ${\hat{x}}_k$ is the estimate of the state x _k in this paper, and it is also the estimate closest to the true value with the minimum mean-square error.

Convergence analysis

According to the given algorithm, the expectation of an arbitrary function $t(·)$ can be calculated by the following formula

$$E[t\left({\hat{x}}_k\right)]=\frac{\sum _{i=1}^{N_p}t\left(x_k^i\right)w_k^i}{\sum _{i=1}^{N_p}{w}_k^i}.$$

(24)

By substituting the particle weight update (21) into (24), we have

$$E[t\left({\hat{x}}_k\right)]=\frac{\sum _{i=1}^{N_p}t\left(x_k^i\right)w_{k-1}^i\frac{p(m_k|x_k)p(x_k|x_{k-1}^i)}{q(x_k|x_{k-1},{\mathcal{M}}_k)}}{\sum _{i=1}^{N_p}{w}_{k-1}^i\frac{p(m_k|x_k)p(x_k|x_{k-1}^i)}{q(x_k|x_{k-1},{\mathcal{M}}_k)}}.$$

According to the Monte Carlo approximation, when the number of particles N _p is sufficiently large, the discrete sum of particles can be transformed into an integral form. Therefore, it follows that

$$E[t\left({\hat{x}}_k\right)]\to \frac{\int t\left(x_k\right)\frac{p(m_k|x_k)p(x_k|{\mathcal{M}}_{k-1})}{q(x_k|x_{k-1},{\mathcal{M}}_k)}q(x_k|x_{k-1},{\mathcal{M}}_k)d{x}_k}{\int \frac{p(m_k|x_k)p(x_k|{\mathcal{M}}_{k-1})}{q(x_k|x_{k-1},{\mathcal{M}}_k)}q(x_k|x_{k-1},{\mathcal{M}}_k)d{x}_k}.$$

After further simplification, the expectation becomes

$$E[t\left({\hat{x}}_k\right)]=\frac{\int t\left(x_k\right)p(m_k|x_k)p(x_k|{\mathcal{M}}_{k-1})d{x}_k}{\int p(m_k|x_k)p(x_k|{\mathcal{M}}_{k-1})d{x}_k}.$$

(25)

Finally, substituting (3) into (25) and simplifying, we obtain

$$E[t\left({\hat{x}}_k\right)]=\int t\left(x_k\right)p(x_k|{\mathcal{M}}_k)d{x}_k.$$

According to the Central Limit Theorem, the Monte Carlo approximation converges as the number of particles increases (i.e. $N_p\to \infty$ ). Therefore, the estimated value ${\hat{x}}_k$ of the system state and the estimated value P _k of the covariance both converge to their respective minimum mean-square error estimates $E[p(x_k|{\mathcal{M}}_k)]$ and $E[(x_k-E[p(x_k|{\mathcal{M}}_k)]){(x_k-E[p(x_k|{\mathcal{M}}_k)])}^T|{\mathcal{M}}_k]$ .

• Remark ³⁵ : investigates the state estimation problem of linear systems with probabilistic quantization errors and bit-flipping scenarios. It has been shown that, according to the measurement ${\mathcal{M}}_k$ with quantization errors and bit-flipping, the estimation of the system state x _k can achieve good results. Meanwhile, an in-depth analysis of the estimation error range has also been carried out.

Specifically, based on Theorem 2 and Theorem 3 in the stability analysis (Section C, Section 4 of ³⁵ ), we can obtain that the upper bound of the covariance $P_{k|k-1}$ is $\overline{\mathrm{M}}$ , that is, $P_{k|k-1}\le \overline{\mathrm{M}}$ , where $\overline{\mathrm{M}}$ satisfies the algebraic equation $\overline{\mathrm{M}}=r_{\beta }\left(\overline{\mathrm{M}}\right)$ . Furthermore, from the covariance update equation in Section 3 of, ³⁵ the covariance P _k at time k has an upper bound. The Cramér–Rao (C-R) lower bound of P _k is given by the inverse of the Fisher information matrix $\left({\mathbf{J}}_k\right)$ , that is, $P_k\ge {\mathbf{J}}_k^{-1}$ . Therefore, the covariance P _k at time k is bounded. For more details, see reference. ³⁵ Similarly, it can be obtained that the estimation error of the improved particle filter algorithm in this paper is also bounded.

Thus, it can be shown that the algorithm designed in this paper has a favorable and stable estimation effect, and the obtained estimation result can converge to the estimate of the minimum mean square error.

Numerical simulation

In this section, the performance of the proposed method is evaluated through a target tracking example. Assume that there are four range-measuring sensors, whose positions are as follows: Sensor 1 is located at (0 cm, 0 cm), Sensor 2 at (0 cm, 120 cm), Sensor 3 at (120 cm, 0 cm), and Sensor 4 at (120 cm, 120 cm). The motion of the target is modeled using the Constant Turn (CT) model described below:

$$\begin{array}{c}{x}_k=\left[\begin{array}{cccc}1\hfill & \frac{\sin \left(\alpha T\right)}{\alpha }\hfill & 0\hfill & \frac{\cos \left(\alpha T\right)-1}{\alpha }\hfill \\ \multicolumn{1}{c}{0}& \cos \left(\alpha T\right)\hfill & 0\hfill & \sin \left(\alpha T\right)\hfill \\ \multicolumn{1}{c}{0}& \frac{1-\cos \left(\alpha T\right)}{\alpha }\hfill & 1\hfill & -\frac{\sin \left(\alpha T\right)}{\alpha }\hfill \\ \multicolumn{1}{c}{0}& -\sin \left(\alpha T\right)\hfill & 0\hfill & \cos \left(\alpha T\right)\hfill \end{array}\right]x_{k-1}+G{w}_k\hfill \\ \multicolumn{1}{c}{z_k^r=\sqrt{{(x_{p,k}-x_r)}^2+{(y_{p,k}-y_r)}^2}+v_k^r,}\end{array}$$

where $G={\left[\begin{array}{cccc}0\hfill & 0\hfill & T\hfill & 1\hfill \\ \multicolumn{1}{c}{T}& 1\hfill & 0\hfill & 0\hfill \end{array}\right]}^T$ , the angular velocity is $\alpha =\pi /20\mathrm{rad}/\mathrm{s}$ and the sampling period is $T=0.5\mathrm{s}$ . The system state is $x_k=[x_{p,k},x_{v,k},y_{p,k},y_{v,k}{]}^T$ , where x_p,k and y_p,k are the positions of the target, x_v,k and y_v,k are the velocities of the target, x _r and y _r represent the positions of the r-th sensor in the x-axis and y-axis. Meanwhile, the initial state is assumed to be $x_0=[60\mathrm{cm},0\mathrm{cm}/\mathrm{s},30\mathrm{cm},0\mathrm{cm}/\mathrm{s}{]}^T$ . Let the number of quantization bits be denoted as L= 3, the number of intervals is $M=2^L$ , and the probability of bit-flipping is set to p=0.01. For the particle filter, we set the number of particles to N _p = 1000 and the resampling threshold to $0.5\times N_p$ . The number of Monte Carlo simulation runs is 200. The variance of the process noise is set to ${\Sigma }_w=0.1\times \mathrm{diag}[1{\mathrm{cm}}^2,0.{5\mathrm{cm}}^2/{\mathrm{s}}^2]$ , and the measurement noise variance is specified as R=0.1. The measurement noise is as follows

$$v_k~0.8\mathcal{N}(0,R)+0.2\mathcal{N}(0,\lambda R),$$

where λ is the scale factor and $\lambda >>1$ .

Next, the Root Mean Square Error for position (RMSEp) estimation is defined ⁴³

$$\mathrm{RMSEp}=\frac{1}{\rho }\times \sum _{i=1}^{\rho }\sqrt{{(x_{p,k}-{\hat{x}}_{p,k})}^2+{(y_{p,k}-{\hat{y}}_{p,k})}^2},$$

(26)

where ρ= 200 is the number of Monte Carlo simulation runs.

With the aid of Matlab, the simulation results are presented in Figures 3 to 6.

Figure 3.

Sensor 1: normalized measurement before/after quantization and bit-flipping.

Figure 4.

x & y Axes: estimated versus true states.

Figure 5.

xoy Plane: estimated versus true states.

Figure 6.

RMSEp curve.

As can be observed from Figure 3, the amount of information after probabilistic quantization is significantly smaller than in the case without probabilistic quantization. Meanwhile, some transmission errors are introduced due to bit-flipping. Figures 4 and 5 intuitively demonstrate the dynamic comparison between the true system states and the estimated states obtained by the improved particle filtering algorithm. These figures present a multidimensional comparison of the true trajectories of the system state variables and the estimated outputs provided by the improved particle filter algorithm. Figures 4 and 5 clearly show that, under the combined effects of quantization errors, bit-flipping, and heavy-tailed noise, the estimated trajectories from the improved particle filter closely follow the actual state trajectories. Whether during the acceleration or deceleration phases of the target’s motion, or at abrupt turning points, the algorithm consistently tracks the underlying state with high accuracy.

Figure 6 provides a quantitative evaluation of the Root Mean Square Error (RMSE) associated with the improved particle filtering algorithm. By plotting the RMSE over time, the figure visually demonstrates the stability of the algorithm’s performance. As shown in Figure 6, the proposed algorithm exhibits strong real-time responsiveness and maintains a high level of estimation accuracy, even in complex environments and over extended operational periods.

In summary, the proposed algorithm demonstrates excellent estimation performance for nonlinear systems subject to quantization errors, bit-flipping, and heavy-tailed noise.

Conclusion

This paper investigates the problem of state estimation for nonlinear systems in the presence of binary-encoding-based quantization, bit-flipping, and heavy-tailed noise. First, considering the limitations of traditional estimation methods under such complex conditions, the particle filter is selected as the basis for algorithmic improvement. The state estimation is performed using Bayes’ formula combined with Monte Carlo simulation. To address the commonly encountered particle degeneracy issue in particle filtering, a Gaussian approximation of the posterior probability density function is adopted as the proposal distribution, thereby enhancing the algorithm’s stability. In the process of designing the proposal distribution, a normalized innovation sequence is defined to effectively handle the effects of quantization errors and bit-flipping. This sequence serves as a bridge for incorporating the effects of quantization and bit-flipping into the estimation process. Based on this construction, the mean, covariance, and likelihood function of the proposal distribution are derived, enabling the design of a complete proposal distribution. Subsequently, particles are sampled from this distribution, and an updated formula for particle weights is derived. Finally, the optimal state estimation is obtained by computing the weighted sum of the particles. Future work may also explore the integration of adaptive quantization bit allocation into the state estimation framework, aiming to balance communication bandwidth efficiency and estimation accuracy more flexibly under dynamic network conditions.