Wrapped Particle Filtering for Angular Data

Particle filtering is probably the most widely accepted methodology for general nonlinear filtering applications. The performance of a particle filter critically depends on the choice of proposal distribution. In this paper, we propose using a wrapped normal distribution as a proposal distribution for angular data, i.e. data within finite range $(-\pi, \pi]$ . We then use the same method to derive the proposal density for a particle filter, in place of a standard assumed Gaussian density filter such as the unscented Kalman filter. The numerical integrals with respect to wrapped normal distribution are evaluated using Rogers-Szegő quadrature. Compared to using the unscented filter and similar approximate Gaussian filters to produce proposal densities, we show through examples that wrapped normal distribution gives a far better filtering performance when working with angular data. In addition, we demonstrate the trade-off involved in particle filters with local sampling and global sampling (i.e. by running a bank of approximate Gaussian filters vs running a single approximate Gaussian filter) with the former yielding a better filtering performance than the latter at the cost of increased computational load.


I. INTRODUCTION
The popular data-based analytical problem is a latent state estimation or filtering, which determines the internal or hidden states of a dynamic system by recursively combining noisy measurements and model-based prediction [1]. Some popular scientific areas involving filtering applications are target tracking [1], biomedical modeling [2], mathematical finance [3], and industrial diagnosis and prognosis [4], [5].
The filtering literature mainly begins with the linear Kalman filter (KF) [6], developed in the sixties. In the sixties itself, a preliminary extension of the linear KF, named extended Kalman filter (EKF) [7], [8], was introduced for nonlinear dynamical systems. Although the linear KF is optimal (in the sense of being conditional mean estimator) for linear Gaussian systems, the nonlinear EKF has no such optimality property due to the linear approximation involved for the nonlinearities. This approximation often results in poor The associate editor coordinating the review of this manuscript and approving it for publication was Shanying Zhu . accuracy and numerical instability. Later, advanced approximate Gaussian filtering techniques were introduced, which include the unscented KF [9], the Gauss-Hermite filter [10], and others [11], [12] for nonlinear systems to improve the accuracy and the numerical stability while still staying in the approximate Gaussian filtering framework. However, all approximate Gaussian filters work poorly in the presence of severe nonlinearities and/or non-Gaussian data [5], [8]. Particle filtering [13] was the next major advancement, which uses recursive Monte Carlo simulation, with the probability density being adjusted using measurements at each time step. This adapts well to nonlinear and non-Gaussian systems. Other recent development include heuristics for filtering with irregular measurements [14], [15] and assumed Gaussian density filtering with non-Gaussian noise [16].
A specific, but practically important filtering problem is the filtering with angular data [17], [18], where the state and measurement appear on circular path. An example is the angle estimation problem of a robot arm [19]. Applications in robotics often require closed-loop position control of the VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ robot arm, requiring a measurement of the current angle of the arm as feedback. This, in turn, involves estimate of an angular quantity from noisy sensor measurements. Modelling circular data in terms of angles generally leads to a nonlinear filtering problem. For nonlinear filtering, two methodologies are popular in literature, Gaussian filtering [11], [12] and particle filtering [13], [20]. The particle filtering can be significantly more accurate at the cost of increased computation and is the main focus of this paper.
The particle filtering represents the desired internal states of dynamical systems in terms of their posterior conditional probability density functions (PDF) and characterizes the PDFs as weighted sum of particles [20]. For generating the particles and the associated probability weights, particle filtering uses a representative PDF, often known as the proposal or importance density [13]. In the state-of-art particle filtering, the proposal density is approximated as an appropriate Gaussian density [21], [22], [23] and is accordingly characterized by its mean and its covariance. The proposal density may be defined locally or globally [21], [24], [25]. A globally defined proposal density generates all the particles at each time step [24]. However, a locally defined proposal density is used for generating a single particle and we require a bank of local proposal densities for generating all the particles. Some of the popular particle filtering techniques, classified in terms of the underlying proposal density used, include the particle filter with transition density as its proposal (PF), unscented particle filter (UPF) [21], cubature particle filter (CPF) [22], and Gauss-Hermite particle filter (GHPF) [23].
The particle filtering accuracy depends on the appropriateness of the proposal density and the re-sampling techniques. The choice of proposal density is further affected by: i) whether its shape is close to the unknown posterior PDF and ii) whether any multivariate integrals involved in computing its moments can be computed accurately in real time. The literature primarily uses Gaussian shape for the proposal density, while many contributions appeared in the literature by introducing different numerical approximation techniques [21], [22], [23], [25], [26], [27], [28]. Similarly, the literature also witnesses many developments by advancing the re-sampling techniques [5], [29]. As we will observe in the later parts of this paper, our objective is to improve the shape of the proposal density to better match the anticipated shape of the unknown PDFs. This will also require to modify the numerical approximation technique for computing the moments involved.
In all the above cases, Gaussian proposal density is over the entire real line as its support. However, Gaussian distribution fails to provide a close approximation of the proposal density, if the variables themselves are constrained to a smaller support, e.g. (−π , π] [30]. Thus we need a proposal density specifically designed for angular data. Academic literature witnesses some preliminary developments [19], [31], [32], [33]. An early development [31] utilized a truncated Fourier series with wrapped normal (WN) distribution. However, finite-length truncation of Fourier series affects estimation accuracy. Some of the later developments [32], [33] are designed only for linear system models in angular data. In a further development, Kurz et al. introduced a nonlinear circular filtering method through a series of publications [17], [19], [34]. However, this method is designed for univariate systems, whereas many of the real-life filtering problems are multivariate. Considering the several limitations of the existing methods, efficient angular filtering is still challenging.
In this paper, we propose a novel particle filtering method for handling nonlinear multivariate angular filtering problems. We represent the unknown angular proposal density appearing over (−π, π] with WN distribution, which is a counterpart of the ordinary normal distribution in the range (−π, π] [35]. The parameters of the wrapped normal distribution are the mean and the covariance. In the proposed method, the computation of the mean and the covariance involves intractable integrals of the form 'nonlinear function × wrapped normal distribution'. We use univariate Rogers-Szegő quadrature rule [36], [37] for approximating such intractable integrals and also extend it to the multivariate case using the product rule. Subsequently, we name the proposed filtering method as Rogers-Szegő particle filter (RSPF). We develop and test the proposed RSPF for both the locally and globally generated proposal densities. Further, we simulate the proposed filters for two angular filtering problems and validate its improved accuracy relative to the existing Gaussian proposal density filters from the simulation results.
Summarizing the above discussion, we highlight the following contributions of our paper in comparison to the existing literature on filtering.
• The proposed RSPF considers wrapped normally distributed proposal density, instead of the Gaussian proposal density which is most commonly used in the existing particle filtering algorithms.
• The proposed RSPF utilizes Rogers-Szegő quadrature rule for the numerical approximation of the integrals. To authors' knowledge, this is the first use of this quadrature rule in particle filtering context.
• We demonstrate through comprehensive numerical examples that the proposed RSPF can accurately handle angular data, where the existing filters underperform.
• We formulated the proposed RSPF for both the locally and globally generated proposal density. To authors' knowledge, this is the first explicit comparison of the two different particle filtering paradigms in the literature, which offer a compromise between estimation accuracy and computational cost for the same chosen proposal density. The rest of the paper are organized as follows. Section II discusses the problem formulation, followed by the discussion on the proposed RSPF in Section III. The simulation results are explained and illustrated in Section IV, while Section V provides conclusions of our work.

II. PROBLEM FORMULATION
Our aim is to develop a particle filtering method suitable for angular data. The underlying state space model is represented as where θ k ∈ D n and y k ∈ D d , with D ∈ (−π, π], are state and measurement variables, respectively at k th instant with k ∈ {1, 2, · · · }. φ k : D n → D n and ψ k : D n → D d are general nonlinear functions, representing the state dynamics and the measurement equation, respectively. Finally, η k ∈ D n and ϑ k ∈ D d represent the process and measurement noises, respectively. Our objective is to recursively estimate θ k ∀k ∈ {1, 2, · · · } from a sequentially received set of measurements y k ∀k ∈ {1, 2, · · · }. The estimation of θ k from y k requires characterizing the PDF P(θ k |y 1:k ) analytically. In this regard, the particle filter [20] approximates P(θ k |y 1:k ) as a weighted summation of stochastically generated sample points, also known as particles. Particles are sampled from an appropriately chosen proposal density q(θ k |y 1:k ) for representing the unknown PDF P(θ k |y 1:k ). We denote θ i k and ω i k , ∀i ∈ {1, 2, · · · , N }, as the i th particle and the associated weight, respectively, at k th instant. The weights are normalized, i.e., N i=1 ω i k = 1. Subsequently, the desired PDF P(θ k |y 1:k ) is approximated as where P(y k |θ i k ) represents the likelihood function, δ(·) represents dirac delta function, and N represents the number of particles. If q(θ k |y 1:k ) is the same as the transition density P(θ k |θ k−1 ), then the weights can be recursively updated as ω i k ∝ ω i k−1 P(y k /θ i k ). However, using transition density as proposal ignores the information inferred from the latest measurement, and can result in poor approximation of P(θ k |y 1:k ) [22]. Using the proposal density resulting from approximate Gaussian filters such as the unscented Kalman filter, i.e., the UKF [21] or the cubature Kalman filter, i.e., the CKF [22], include the latest measurements information during generation of particles. There are two fundamentally different approaches for generating particles. In local sampling, [21] we approximate q(θ i k |θ i 0:k−1 , y 1:k ) = P(θ i k |θ i k−1 , y 1:k ) ≈ N (θ i k|k , P i k|k ) for each particle, where N (·) denotes the Gaussian distribution andθ i k|k , P i k|k , respectively, represent i th posterior mean and the posterior covariance estimates at time k. This update of density for each individual particle effectively involves running a bank of N individual Gaussian filters. Alternatively, [24], we can approximate q(θ i k |θ i 0:k−1 , y 1:k ) = P(θ i k |θ i k−1 , y 1:k ) ≈ N (θ k|k , P k|k ) for all the particles. This requires a single Gaussian proposal at each time step and is termed as global sampling. The authors have not come across explicit comparison of global and local approaches to proposal density generation. Most of the subsequent discussion in this paper follows local sampling strategy. However, we also compare our results using the global sampling strategy with those using the local sampling strategy in our simulation examples.
Note that Gaussian proposal density characterizes any data over entire real line. This is inappropriate characterization of the angular data bounded within (−π, π] and can result in poor accuracy for angular filtering. Here, we propose a novel particle filtering algorithm for angular filtering by characterizing the angular PDFs over (−π, π] instead of the entire real line. As discussed above, we also illustrate the developed particle filtering method for both the locally and globally defined proposal densities.

III. ROGERS-SZEGŐ PARTICLE FILTER
In this section, we introduce the proposed RSPF, which has potential to be far more accurate for angular data appeared over (−π, π]. We approximate the unknown proposal density with wrapped normal distribution 37, defined over (−π, π]. Before proceeding further, we briefly discuss this distribution in the next sub-section.

A. WRAPPED NORMAL DISTRIBUTION
The wrapped normal distribution may be obtained by wrapping the horizontal axis of an ordinary normal distribution curve around a unit circle [19]. Similar to a normal distribution, wrapped normal distribution is also completely characterized by the mean and the variance. If θ is wrapped normally distributed with mean µ and variance σ 2 , then the distribution of θ is given as [17], The zero-mean and unit-variance wrapped normal distribution is plotted in Fig. 1. We refer to [17], [38], [39] for more detailed discussion on the wrapped normal distribution. VOLUME 10, 2022

B. PROPOSAL DENSITY APPROXIMATION
As stated in [38] and [17], the wrapped normal distribution follows circular central limit theorem for angular data. Thus, we consider the following approximations.
• The prior and posterior distributions of θ i k are approximated as wrapped normally distributed, i.e., where f WN (·) represents multivariate wrapped normal distribution. Here,θ i k|k , P i k|k , respectively, represent the i th posterior mean and the posterior covariance of the state at time k.
k|k−1 and P yyi k|k−1 denote the mean and the covariance of y i k , respectively. • The desired proposal density is approximated as and f WN (0, R k ), respectively, where Q k and R k are the covariances of the respective noises.
• Additionally, we assume that the noises, η k and ϑ k are independent of each other as well as serially independent.
From the above assumptions, the proposal density can be characterized by meanθ i k|k and covariance P i k|k for each particle, which corresponds to local sampling in our earlier discussion. As in any standard Gaussian filtering algorithm, we obtainθ i k|k and P i k|k in two steps, prediction and update. The prediction characterizes the PDF P(θ i k |y k−1 ) by computing the meanθ i k|k−1 and covariance P i k|k−1 . The update characterizes the desired PDF P(θ i k |y k ) by determining the meanθ i k|k and the covariance P i k|k approximately, using the formulae for conditional mean and the conditional variance for Gaussian distributions, which is a step similar to other assumed Gaussian density filters [11], [39]. We provide the computational aspects of these two steps in the following discussion.

1) PREDICTION
This step computesθ i k|k−1 and P i k|k−1 , respectively, as [23] 2) UPDATE This step computesθ i k|k and P i k|k , respectively, as [23] θ k|k−1 , and K i k are determined, respectively, as [23] Eqs. (5) to (11) give the steps needed to generate the i th mean-covariance pair θ i k|k , P i k|k , given θ i k−1|k−1 , P i k−1|k−1 and the measurement y k . Note that these equations use the same approximate equations for prior and posterior distributions as in the case of other assumed Gaussian density filters [12], apart from one crucial difference that wrapped normal density is used in place of normal density for calculating moments.
The characterization of the proposal density using the above equations requires computing integrals of the form where F : D n → D n is a general nonlinear function. As closed-form solution is not available for such integrals in general, the proposed RSPF introduces Rogers-Szegő quadrature rule [37], [40] for approximating such integrals. The integrals occurring here are different from those which occur in traditional proposal densities or in approximate Gaussian filters [11], [12] or in traditional proposal densities for particle filters [21], [22], [23], since these involve a Gaussian density (rather than a wrapped normal density).

C. APPROXIMATION OF INTEGRALS
In this section, we introduce the Rogers-Szegő quadrature rule for approximating the desired intractable integral I n (F). While the Rogers-Szegő quadrature rule is applicable only for univariate systems, the desired integral I n (F) is multivariate. We additionally use the product rule [10] for extending the univariate Rogers-Szegő quadrature rule to multivariate case.
We will first approximate the standard form (zero-mean and unity-covariance) of I n (F), given as 90290 VOLUME 10, 2022 where I n represents the identity matrix. Further, the multivariate integral I n 0 (F) can be expressed as where the univariate integrals I 1 0 (F s ) ∀s ∈ {1, 2, · · · , n} can be expressed as with The univariate Rogers-Szegő quadrature rule is designed for approximating I 1 0 (F s ) given in Eq. (15).
The solution of the recurrence relation Eq. (18) corresponding to f WN (θ s ; 0, 1) gives the desired monic sequences in terms of τ -binomial coefficient as 37 where τ -binomial coefficient is , with 0 < τ < 1 and After γ m (Z) is obtained from Eq. (19), the m th Rogers-Szegő polynomial is obtained as where From [41], τ = e −1 corresponding to the desired weight function f WN (θ s ; 0, 1). Utilizing the above interpretations, 37 derives and states that the N s number of desire Rogers-Szegő quadrature points can be obtained as the phase of complex roots of the polynomial where |δ| = 1 and with j ∈ {1, 2, · · · , N s }. Let us denote the roots of A N s (Z) as λ j , ∀j ∈ {1, 2, · · · , N s }. Then, the j th Rogers-Szegő quadrature point is obtained as χ j = λ j ∀j ∈ {1, 2, · · · , N s }. Subsequently, the weight W j associated with χ j can be obtained in terms of monic sequence γ N s as 37 where Re[·] represents the real part of complex number and γ denotes the first derivative of γ . After χ j and W j are obtained the univariate integral I 1 0 (F s ) given in Eq. (15) can be approximated as Table-1 illustrates the complex values of λ, the angular points χ, and the associated weights W for eight-points univariate Rogers-Szegő quadrature rule with τ = e −1 . Please note that τ = e −1 is an essential requirement for the proposed RSPF as mentioned above.

2) MULTIVARIATE EXTENSION OF ROGERS-SZEGŐ QUADRATURE RULE
In the above discussion, we introduced the univariate Rogers-Szegő rule for approximating the univariate integral I 1 0 (F s ). We now extend the univariate Rogers-Szegő rule for approximating the multivariate integral I n 0 (F). The product rule states that given the univariate Rogers-Szegő quadrature points {χ 1 , χ 2 , · · · , χ N s } and the associated weights {W 1 , W 2 , · · · , W N s }, we can approximate I n 0 (F) as where s 1 , s 2 , · · · , s n ∈ {1, 2, · · · , N s }. Finally, we can approximate the desired intractable integral I n corresponding to f WN (θ ;θ, P), as We can rewrite the above expression as where N p = (N s ) n denotes the number of the multivariate Rogers-Szegő quadrature points, while j = χ s 1 χ s 2 , · · · , χ s n and W j = n j=1 W s j , with s 1 , s 2 , · · · , s n ∈ {1, 2, · · · , N s }, represent the multivariate quadrature points and weights, respectively. Moreover, S represents the Cholesky decomposition of P, i.e., P = SS T .
We can use the product rule-based multivariate Rogers-Szegő quadrature rule for approximating the intractable integrals appeared through Eqs. (5) to (11). Subsequently, we can construct the proposal density q(θ i k |θ i 0:k−1 , y 1:k ) ≈ f WN (θ i k ;θ i k|k , P i k|k ) by determining the meanθ i k|k and the covariance P i k|k from Eqs. (7) and (8), respectively. Considering that the multivariate quadrature points and weights, i.e., j and W j , ∀j ∈ {1, 2, · · · , N p }, are available, we refer to [23] for analytical steps used for implementing the prediction and update steps discussed in Section III-B.

D. PARTICLE FILTERING WITH WRAPPED PROPOSAL DISTRIBUTION
We now introduce the proposed RSPF using the wrapped normal distribution for handling angular data. The RSPF uses the proposal density determined through the subsections III-A-III-B in particle filtering. As discussed previously, we develop two filters, one for local sampling and another for global sampling. We abbreviate the two filters as LRSPF and GRSPF, respectively. The previous discussions of this section provides the formulations for the local sampling, i.e. for the LRSPF. The formulation for GRSPF is similar, although simpler due to a single filter at each time step.

1) LRSPF
The implementation of the proposed LRSPF comprises the following steps.
• Initialization: We initialize the filter at k = 0 with initial state θ 0 and initial covariance P 0|0 . Subsequently, we generate the initial set of particles as {θ i 0 } N i=1 ∼ P(θ 0 ) and the associated weights as {ω i 0 } N i=1 = 1/N . • Particles and weights calculation at k th instant (k ≥ 1): Construct the initialization for k th instant from the latest particles and covariance (assignθ .

(27)
• Weight normalization: The weights are normalized as and determine the effective sample size N eff as • Re-sampling: If N eff is below a preassigned threshold value N th , then we perform re-sampling [22] and a new set θ i k ∀i ∈ {1, 2, · · · , N } is generated from the current set of particlesθ i k ∀i ∈ {1, 2, · · · , N }. In general, we consider N th = 2/3N . Furthermore, a new updated weights are obtained as, ω i k = 1/N . • Estimation: Finally, the desired estimate and covariance are obtained as 90292 VOLUME 10, 2022 The LRSPF requires implementing the algorithm discussed through Eqs. (5) to (11) implemented once for each particle, which increases the computational demand.

2) GRSPF
In an alternate method to the LRSPF, a single proposal density is globally defined in order to reduce the computational demand. This single proposal density is used for generating all particles. The steps for implementing the GRSPF is provided below.
• Initialization: This step is same as the LRSPF, which can be followed from the previous discussion.
• Compute the mean and covariance of particles at k th instant (k ≥ 1): -Propagateθ k−1|k−1 and P k−1|k−1 through the algorithm discussed in Section III-B to obtainθ k|k and P k|k for single mean and covariances. • Compute particles: Generate the desired particles as {θ .
• The remaining steps are same as those discussed for LRSPF. For angular data, the proposed RSPF with wrapped normal proposal and Rogers-Szegő quadrature rule for integration is more accurate than the numerical approximation methods used in the state-of-art particle filters such as the UPF, CPF, and GHPF under assumed Gaussian density, as amply illustrated in the next section. The LRSPF is slightly more accurate than the GRSPF due to the local sampling and may be preferred if computational budget is available.
For angular data, the wrapped normal distribution is probably the most appropriately shaped distribution for which numerical approximation methods to the desired multivariate integrals are available. We have used this intuition to develop a new particle filtering algorithm. In the next section, we demonstrate its superior performance on simulation examples on multivariate nonlinear systems where data is inherently angular.

IV. SIMULATION EXAMPLES
In this section, we validate the performance of the proposed RSPF for two angular filtering problems. The simulation is performed in MATLAB on a PC with Intel Core i5, 7th gen processor running at 3.40 GHz, 32 GB RAM, and a 64-bit operating system.
We compare the performance of the proposed LRSPF and GRSPF with various existing filters, including the particle filter with transition density as its proposal (referred to as PF) [20], unscented particle filter or the UPF [21], cubature particle filter or the CPF [22], and Gauss-Hermite particle filter or the GHPF [23], and the existing nonlinear circular filter (CF) [19]. The existing filters, UPF, CPF, and GHPF are implemented with local sampling to characterize them with their best accuracy. Note that CF is applicable only for univariate systems. Thus, it will be included in comparison for the first simulation problem only, as the second problem is multivariate.
The proposed LRSPF, GRSPF, and GHPF are implemented with two univariate quadrature points. The κ parameter for implementing the UPF is assigned as κ = 4. Finally, CF is applied by considering a three-point wrapped Dirac mixture distribution.
The performance analysis and comparison are based on angular root mean square error (RMSE) between the true and estimated states. Note that the true and the estimated states may often fall beyond the angle range (−π, π] due to noise. In such cases, we perform aliasing to obtain the equivalent angle within (−π, π]. The angular root mean square error (RMSE) is expressed as, where θ i k andθ i k|k are the true and the estimated states at k th time-step and in i th simulation run, and MC is the number of Monte-carlo simulations.

A. PROBLEM 1
This is an angle estimation problem of a robot's arm [19], which is briefly discussed in the introduction section. The dynamic state space model of the system is given as [19] where d 1 and d 2 are constants. The initial true and estimated states are taken as θ 0 = 0 and θ 0|0 = π, respectively, while the initial variance is taken as P 0|0 = 2. We assign d 1 = 0.1, d 2 = 0.15, Q = 0.1, and R = diag(0.2, 0.2). The simulation is performed for 200 time-steps and angular RMSEs are computed by implementing 1000 Monte-Carlo simulations. Fig. 2 shows the angular RMSE plots and Table-2 shows the average angular RMSEs of the LRSPF, GRSPF and the existing filters using different numbers of particles. The figure and table together indicate that RMSE is significantly reduced for both the LRSPF and GRSPF, as compared to the existing choices of proposal densities, viz PF, UPF, CPF and GHPF. More specifically, Table 2 concludes that the average angular RMSEs of LRSPF and GRSPF (100 particles)  are reduced by 26.45% and 18.05%, respectively compared with the PF. We can observe a similar relative percentage of improvement for the proposed filters compared with the remaining existing filters. However, it is also interesting to note that the existing circular filter CF outperforms the other existing filters.
The RMSE of the GRSPF is higher than the LRSPF, as expected. However, even global sampling using wrapped normal distribution as a proposal (i.e. the GRSPF) leads to a better estimation accuracy than Gaussian density proposal filters using local sampling, even as local sampling leads to a significantly higher computational cost.
Relative computational times of different filters are listed in Table 3. Thus, the computational time for the proposed GRSPF is lower than all the algorithms which require local sampling, while it still delivers a superior estimation performance to PF, UPF, CPF and GHPF. CF is not a particle filtering algorithm and has a constant computational cost as well as constant RMSE, relative to the number of particles. Note that global sampling in UPF, CPF and GHPF leads to poorer accuracy than that using local sampling, and hence, the results of the global sampling for these algorithms are omitted. LRSPF further increases the accuracy of our  algorithm, at the expense of somewhat higher computational cost.

B. PROBLEM 2
The second problem considered is a general multivariate nonlinear angular estimation problem [14]. In this problem, the state dynamics are of oscillatory nature, while the measurement equation is a monotone increasing function of arguments (e.g. a positive quadratic form or its positive square root). Similar, system models often appear in sonar-based bearing measurements and GPS (Global positioning system)based information on the angle of arrival. This problem has been widely used in literature [11], [14], [44], [45] for validating the filtering performance. The state-space model  of this problem can be written as [14], We consider a three-dimensional system (θ k ∈ D 3 , y k ∈ D) and assign the initial true and estimated states as θ 0 = [0, −π, π] T andθ 0|0 = [−π, π, π/2] T , respectively, while initial error covariance is taken as P 0|0 = 2I n . The noise covariances are assigned as Q = diag([0.05, 0.05, 0.05]) and R = 0.1. The states are estimated for 200 time-steps and the results are evaluated by performing 1000 Monte-Carlo runs. Fig. 3 shows the angular RMSE plots of the three states obtained using the LRSPF, GRSPF and the existing filters. Moreover, Table-4 shows the average angular RMSEs obtained using each filter for varying number of particles. Note that the figures and table do not include the results for CF, as this problem is multivariate. As in the case of problem 1, LRSPF and GRSPF yield lower average angular RMSE than all the other filters examined. Furthermore, Table 4 concludes that the average angular RMSEs of LRSPF and GRSPF (100 particles) are reduced by 22% and 18.25%, respectively compared with the PF. We can conclude a similar relative percentage of accuracy improvement for the proposed filters compared to other filters as well. The order of computational times of all filters was observed similar to the order reported for the first problem. In this case as well, GRSPF gives a better accuracy than existing Gaussian proposal filters at a lower computational cost.

C. NOISE PARAMETER VARIATIONS
We further study the effect of varying process and measurement noise parameters on the performance of the proposed RSPF. In this regard, we define three different scenarios by varying the process and measurement noise covariances Q and R.  Fig. 4 for Problem 1 and in Fig. 5 for Problem 2. In all the cases, the average angular RMSE is lower for the LRSPF and the GRSPF as compared to the other filters.

D. DISCUSSION
The simulation results infer that the proposed RSPF outperforms the traditional PF, the existing extensions of the PF (such as the UPF, CPF, and GHPF), and the circular filter CF for angular data. Interestingly, the proposed RSPF could outperform all the existing forms of PF at lower computational demand. We also introduce two forms of the proposed RSPF, abbreviated as LRSPF and GRSPF. They can be used to achieve a different trade-off between the accuracy and computational demand, particularly if the number of particles is high. We also observe that the proposed RSPF outperforms all the existing filters for various noisy environments, which validates the improved accuracy of the proposed method under different practical scenarios.

V. CONCLUSION
The particle filtering is a popular and widely accepted nonlinear filtering method available in literature. A crucial deter-VOLUME 10, 2022 minant of its performance is the choice of proposal density. Most existing particle filters use Gaussian proposal and tend to perform poorly for estimation of angular quantities on a restricted domain. The proposed RSPF uses wrapped normal distribution instead of the ordinary normal distribution as the proposal density. Subsequently, it closely represents the proposal density for angular data. Further, we propose to use i) conditional density approximation similar to the Gaussian filtering case to derive posterior densities and ii) use Rogers-Szegő quadrature rule along with the product rule for approximating the integrals involved. We compare the performance of two variants of the new filter (with local sampling and with global sampling) with existing Gaussian filters. On two different simulation examples, the proposed filter is shown to outperform Gaussian proposal filters comprehensively. We show that even a global sampling variant of our filter is more accurate than the local sampling versions of existing Gaussian filters, even though global sampling leads to a significantly reduced computational cost. Local sampling version of the RSPF can lead to some what increased accuracy at the cost of a higher computational load. Additionally, the product rule suffers from the curse of dimensionality problem [46]. However, the practitioners may replace the product rule with the Smolyak rule [46] and adaptive sparse-grid method [47] to partially address this problem.

DATA ACCESS STATEMENT
This research did not use any experimentally generated data or data from any publicly available dataset. Model definitions (including the specified probability distributions) and parameter values (including the initialization parameters) provided in the paper are adequate for reproducing the exact qualitative behavior of the algorithms illustrated in the paper.