Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments

A sequential Bayesian estimation algorithm for multipath mitigation is presented, with an underlying movement model that is especially designed for dynamic channel scenarios. In order to facilitate efficient integration into receiver tracking loops, it builds upon complexity reduction concepts that previously have been applied within maximum likelihood (ML) estimators. To demonstrate its capabilities under different GNSS signal conditions, simulation results are presented for both BPSK-modulated and BOC-(1,1) modulated navigation signals.


INTRODUCTION
A major error source within global navigation satellite systems (GNSSs) comes from multipath, the reception of additional signal replica due to reflections, which introduce a bias into the estimate of the delay lock loop (DLL) of a conventional navigation receiver. For efficient removal of this bias, it is possible to formulate advanced maximum likelihood (ML) estimators that incorporate the echoes into the signal model and are capable of achieving the theoretical limits given by the Cramér Rao bound. The drawback of ML estimator techniques is that the parameters are assumed to be constant during the time of observation. Independent estimates are obtained for successive observation intervals, whose length has to be adapted to the dynamics of the channel.
In this paper, we consider the important practical case of dynamic channel scenarios and assess how the timedelay estimation can be improved if information is available about the temporal evolution of the channel parameters. Our approach is based on Bayesian filtering, the optimal and wellknown framework to address such dynamic state estimation problems [1]. Sequential Monte Carlo (SMC) methods [2,3] are used for computing the posterior probability density functions (PDFs) of the signal parameters. Figure 1 gives an overview of the relationships between different multipath mitigation and estimation approaches. In fact, we have chosen to discriminate approaches according to their primary objective. The left column represents the class of techniques that attempt to mitigate the effect of multipath in different ways. This can for example be achieved by modifications of the antenna response, either by means of hardware design or with signal processing techniques (e.g., beamforming). The majority of the remaining mitigation techniques are in some way aligning the more or less traditional receiver components (e.g., the early/late correlator) to the signal received in the multipath environment.

A COMPARISON OF VARIOUS MULTIPATH MITIGATION APPROACHES
The probably most simple technique is the adjustment of the correlator spacing applied in the Narrow Correlator [4]. Other well-known examples of this category are the Strobe Correlator [5], the Gated Correlator [6], or the Pulse Aperture Correlator [7]. To incorporate new signal forms (such as BOC), these methods need "tuning" in order to suffer as little as possible from multipath. On the other hand, multipath estimation techniques (right column) treat multipath (in particular the delay of the paths) as something to be estimated from the channel observations, so that its  effects can be trivially removed at a later processing stage. For the estimation techniques, we have differentiated between static and dynamic approaches, according to the underlying assumption of the channel dynamics. Examples for static multipath estimation are those belonging to the family of ML estimators, often using different efficient maximization strategies over the likelihood function [8][9][10][11][12][13]. For static channels without availability of prior information, the ML approach is optimal and performs significantly better than other techniques, especially if the echoes have short delay. An estimator based on sequential importance sampling (SIS) methods (particle filtering) for static multipath scenarios has been considered in [14], which has the advantage that prior channel knowledge can be incorporated.
As a first step towards addressing dynamic channels, one can incorporate ML estimators in receiver loops or formulate quasisequential estimators [15,16]. Finally, dynamic estimators that target the computation of the posterior PDF conditioned on the received channel output sequence at the receiver can be applied on a per single range basis or operate in the position domain. In this paper, we concentrate on dynamic estimators applied per each range. The sequential Monte Carlo approach has also been suggested in the communications field for estimation of time-varying channel responses in spread spectrum systems [17,18].

SIGNAL MODEL
Assume that the complex valued baseband-equivalent received signal is equal to (1) where c(t) is a delta-train code sequence that is modulated on a pulse g(t), N m is the maximum number of allowed paths reaching the receiver (to restrict the modeling complexity), e i (t) is a binary function that controls the activity of the i th path, and a i (t) and τ i (t) are their individual complex amplitudes and time delays, respectively. The signal is disturbed by additive white Gaussian noise n(t). Grouping blocks of L samples at times (m + kL)T s , m = 0, . . . , L − 1, together into vectors z k , k = 0, 1, . . . , n whilst assuming the parameter functions e i (t), a i (t), and τ i (t) being constant within the corresponding time interval and equal to e i,k , a i,k , and τ i,k , this can be rewritten as In the compact form on the right hand side, the samples of the delayed pulses g(τ i,k ) are stacked together as columns of the matrix G(τ k ) = [g(τ 1,k ), . . . , g(τ Nm,k )], C is a matrix representing the convolution with the code, and the delays and amplitudes are collected in the vectors τ k = [τ 1,k , . . . , τ Nm,k ] T and a k = [a 1,k , . . . , a Nm,k ] T , respectively. Furthermore, for concise notation we use E k = diag [e k ] whilst the elements of the vector e k = [e 1,k , . . . , e Nm,k ] T , e i,k ∈ [0, 1], determine whether the i th path is active or not by being either e i,k = 1 corresponding to an active path or e i,k = 0 for a path that is currently not active. The term s k denotes the signal hypothesis and is completely determined by the channel parameters τ k , a k, and e k . Using (2), we can write the associated likelihood function as The likelihood function will play a central role in the algorithms discussed in this paper; its purpose is to quantify the conditional probability of the received signal conditioned on the unknown signal (specifically the channel parameters).

Efficient likelihood computation
In [11], a general concept for the efficient representation of the likelihood (3) was presented, which is applicable to many of the existing ML multipath mitigation methods.
The key idea of this concept is to formulate (3) through a vector z c,k resulting from an orthonormal projection of the observed signal z k onto a smaller vector space, so that z c,k is a sufficient statistic according to the Neyman-Fisher factorization [19] and hence suitable for estimating s k . In other words, the reduced signal comprises the same information as the original signal itself. In practice, this concept becomes relevant as the projection can be achieved by processing the received signal (2) with a bank of correlators and a subsequent decorrelation of the correlator outputs. A variant of this very general concept, applied in [13], has also been referred to as the Signal Compression Theorem in [20] for a set of special projections that do not require the step of decorrelation due to the structure of the used correlators. For instance, unlike the correlation technique used in [8], the one suggested in [13] already projects onto an orthogonal and thus uncorrelated subspace, similar to the code matched correlator technique proposed in [11]. Due to complexity reasons, all practically relevant realizations of ML estimators [8,13] operate in a projected space, namely after correlation. The corresponding mathematical background will be discussed below, including also interpolation of the likelihood and elimination of complex amplitudes as further methods for complexity reduction.

Data compression
As explained above, the large vector containing the received signal samples z k is linearly transformed into a vector z c,k of much smaller size. Following this approach, the likelihood according to (2) can be rewritten as with the compressed received vector z c,k and the compressed signal hypothesis s c,k : and the orthonormal compression matrix Q c , which needs to fulfill to minimize the compression loss. According to [11], the compression can be two-fold so that we can factorize into a canonical component decomposition, given by an L×N cc matrix Q cc , and a principal component decomposition, given by an N cc × N pc matrix Q pc . In [11], two choices for Q cc are proposed Signal matched, where the elements of the vector τ b define the positions of the individual correlators. The noise-free outputs of the corresponding correlator banks are illustrated in Figure 2. To decorrelate the bank outputs (CG(τ b )) H y and C(τ b ) H y, as mentioned above, the whitening matrix R cc can be obtained from a QR decomposition of CG(τ b ) and C(τ b ), respectively. Apart from practical implementation issues, both correlation methods given by (8) are equivalent from a conceptual point of view. For details on the compression through Q pc, the reader is referred to [11].

Interpolation
In order to compute (4) independently of the sampling grid, advantage can be made of interpolation techniques. Using the discrete Fourier transformation (DFT), with Ψ being the DFT matrix and Ψ −1 being its inverse counterpart (IDFT), we get with Ω(τ k ) being a matrix of column-wise stacked vectors with Vandermonde structure [10,11], such that the element at row p and column q computes with N g is the length of the pulse g in samples. The advantage of the interpolation is that it can take place in the reduced space. The most costly computations in (9) can be carried out in precalculations as the matrix M sc , whose row dimension corresponds to the dimension of the reduced space and whose column dimension is N g , is constant.

Amplitude elimination
In a further step, we reduce the number of parameters by optimizing (4) for a given set of τ k and e k with respect to the complex amplitudes a k , which can be achieved through a closed-form solution. Using and obtaining S + c,k by removing zero columns from S c,k, one yields the corresponding amplitude values of the active paths: 4 International Journal of Navigation and Observation  As we have introduced a potential source of performance loss by eliminating the amplitudes and thus practically are disregarding their temporal correlation, we propose to optimize (4) using with the adjustable averaging coefficient Q. When evaluating (4), we substitute s c,k by where the elements of the vector a k that are indicated to have an active path (a k,i : i→e k,i = 1) are set equal to the corresponding elements of a + k . All other elements (a k,i : i→e k,i = 0) can be set arbitrarily as their influence is masked by the zero elements of e k .

Review of the ML Concept
The concept of ML multipath estimation has drawn substantial research interest since the first approach was proposed in [8]. Despite being treated differently in various publications, the objective is the same for all ML approaches, namely to find the signal parameters that maximize (3) for a given observation z k : The signal parameters are thereby assumed as being constant throughout the observation period k. Different maximization strategies exist, which basically characterize the different approaches. Despite offering great advantages for theoretical analysis, the practical advantage of the generic ML concept is questionable due to a number of serious drawbacks.
(i) The ML estimator assumes that the channel is static for the observation period and is not able to exploit its temporal correlation throughout the sequence k = 1, . . . , n. Measured channel scenarios have shown significant temporal correlation [21].
(ii) Despite being of great interest in practice, the estimation of the number of received paths is often not addressed. The crucial problem here is to correctly estimate the current number of paths to avoid over determination, since an overdetermined estimator will tend to use the additional degrees of freedom to match the noise by introducing erroneous paths.
There exist various techniques based on model selection that can be employed to estimate the number of paths [22] but they suffer from the problem of having to dynamically adjust the decision thresholds. Typically, only a single hypothesis is tracked, which in practice causes error event propagation.
(iii) The ML estimator does only provide the most likely parameter set for the given observation. No reliability information about the estimates is provided. Consequently, ambiguities and multiple modes of the likelihood are not preserved by the estimator.
ML estimators require that the estimated parameters remain constant during the observation period. Due to data modulation and phase variations in practice, this period, which is often referred to as the coherent integration time, is limited to a range of 1 millisecond-20 milliseconds. To reach a sufficient noise performance with a ML estimator in practice, it is required to extend its observation interval approximately to the equivalent averaging time of a conventional tracking loop, which is commonly in the order of several hundred coherent integration periods. To overcome this problem, the observations are forced to be quasicoherent by aiding the ML estimator with a phase locked loop (PLL) and a data removal mechanism [8].

Optimal solution
In Section 3, we have established the models of the underlying time-variant processes. The problem of multipath mitigation now becomes one of sequential estimation of a hidden Markov process. We want to estimate the unknown channel parameters based on an evolving sequence of received noisy channel outputs z k . The channel process for each range of a satellite navigation system can be modeled as a first-order Markov process if future channel parameters given the present state of the channel and all its past states depend only on the present channel state (and not on any past states). We also assume that the noise affecting successive channel outputs is independent of the past noise values; so each channel observation depends only on the present channel state.
Intuitively, we are exploiting not only the channel observations to estimate the hidden channel parameters (via the likelihood function), but we are also exploiting our prior knowledge about the statistical dependencies between successive sets of channel parameters. We know from channel Measurements (observed) State transition PDF Hidden states measurements that channel parameters are time varying but not independent from one time instance to the next; for example, an echo usually experiences a "life-cycle" from its first occurrence, then a more or less gradual change in its delay and phase over time, until it disappears [21]. These measurements also show that the common channel models considered in communication systems [17,18] are not adequately reflecting the properties that are crucial for highresolution signal delay estimation as required in navigation systems. Now that our major assumptions have been established, we may apply the concept of sequential Bayesian estimation. The reader is referred to [23] which gives a derivation of the general framework for optimal estimation of temporally evolving (Markovian) parameters by means of inference, and we have chosen similar notation. The entire history of observations (over the temporal index k) can be written as Similarly, we denote the sequence of parameters of our hidden Markovian process by Here, x k represents the characterization of the hidden channel state, including the parameters that specify the signal hypothesis s k given in (2). Our goal is to determine the posterior PDF of every possible channel characterization given all channel observations: p(x k | Z k ) (see Figure 3). Once we have evaluated this posterior PDF, we can either determine that channel configuration that maximizes it-the so-called maximum a posteriori (MAP) estimate; or we can choose the expectation-equivalent to the minimum mean square error (MMSE) estimate. In addition, the posterior distribution itself contains all uncertainty about the current range and is thus the optimal measure to perform sensor data fusion in an overall positioning system.
It can be shown that the sequential estimation algorithm is recursive, as it uses the posterior PDF computed for time instance k − 1 to compute the posterior PDF for instance k (see Figure 4). For a given posterior PDF at time instance k −1, p(x k−1 | Z k−1 ), the prior PDF p(x k | Z k−1 ) is calculated in the so-called prediction step by applying the Chapman-Kolmogorov equation: with p(x k | x k−1 ) being the state transition PDF of the Markov process. In the update step, the new posterior PDF for step k is obtained by applying Bayes' rule to p(x k | z k , Z k−1 ) yielding the normalized product of the likelihood p(z k | x k ) and the prior PDF: The term p(z k | x k ) = p(z k | s c,k = s c,k ) follows from (4) and represents the probability of the measured channel output (often referred to as the likelihood value), conditioned on a certain configuration of channel parameters at the same time step k. The denominator of (19) does not depend on x k , and so it can be computed by integrating the numerator of (19) over the entire range of x k (normalization).
To summarize so far, the entire process of prediction and update can be carried out recursively to calculate the posterior PDF (19) sequentially, based on an initial value of p(x 0 | z 0 ) = p(x 0 ). The evaluation of the likelihood function p(z k | x k ) is the essence of the update step. Similarly, maximizing this likelihood function (i.e., ML estimation) would be equivalent to maximizing p(x k | Z k ) only in the case that the prior PDF p(x k | Z k−1 ) does not depend on Z k−1, and when all values of x k are a priori equally likely. Since these conditions are not met, evaluation of p(x k | Z k ) entails all the above steps.

Sequential estimation using particle filters
The optimal estimation algorithm relies on evaluating the integral (18), which is usually a very difficult task, except for certain additional restrictions imposed on the model and the noise process. So, very often a suboptimal realization of a Bayesian estimator has to be chosen for implementation. In this paper, we use a sequential Monte Carlo (SMC) filter, in particular a sampling importance resampling particle filter (SIR-PF) according to [23]. In this algorithm, the posterior density at step k is represented as a sum and is specified by a set of N p particles:  where each particle with index j has a state x j k and has a weight w j k . The sum over all particles' weights is one. In SIR-PF, the weights are computed according to the principle of importance sampling where the so-called proposal density is chosen to be p(x k | x k−1 = x j k−1 ), and with resampling at every time step. For N p →∞, the approximate posterior approaches the true PDF.
The key step, in which the measurement for instance k is incorporated, is in the calculation of the weight w j k which for the SIR-PF can be shown to be the likelihood function: p(z k | x j k ). The characterization of the channel process enters in the algorithm, when at each time instance k, the state of each particle x i k is drawn randomly from the proposal distribution, that is, from p(x k | x j k−1 ).

Exploiting linear substructures
If there exist linear substructures in the model, it is possible to reduce the computational complexity of the filter by means of marginalization over the linear state variables [24], also known as Rao-Blackwellization [25]. In a marginalized filter, particles are still used to estimate the nonlinear states, while for each of the particles the linear states can be estimated analytically. In our case, since the measurement z k is a linear function of the complex amplitudes a k , the likelihood function can be factorized as where the function g(z k , τ k , e k , a k ) is Gaussian with respect to a k , and f (z k , τ k , e k ) is proportional to p(z k | s c,k ). Marginalization of (21) over the linear state variables gives Assuming that the amplitudes are block-wise independent (Q = 1), it follows that the weights of the SIR particle filter are equal to p(z k | x j k ) = p(z k | s c,k = s j c,k ), which is given by (4).
If there exists a temporal correlation of the amplitudes, an optimal marginalized filter requires the implementation of a separate Kalman filter for each of the particles. As we do not want to increase the complexity per particle, we take the correlation into account by adjustment of Q in (13).

Choice of appropriate channel process
To exploit the advantages of sequential estimation for our task of multipath mitigation/estimation, we must be able to describe the actual channel characteristics (channel parameters) so that these are captured by p(x k | x k−1 ). In other words, the model must be a first-order Markov model, and all transition probabilities must be known. In our approach, we approximate the channel as follows.
(i) The channel is totally characterized by a direct path (index i = 1) and at most N m − 1 echoes, (ii) each path has complex amplitude a i,k and relative delay and Δτ i,k = τ i,k − τ 1,k ; where echoes are constrained to have delay τ i,k ≥ τ 1,k , that is, Δτ i,k ≥ 0, (iii) the different path delays follow the process: (iv) each parameter α i,k that specifies the speed of the change of the path delay follows its own process: (v) the magnitudes and phases of the individual paths, represented by the complex amplitudes a i,k , are eliminated according to (12) and (14) for the computation of the likelihood (4), (vi) each path is either "on" or "off," as defined by channel parameter e i,k ∈ {1 ≡ "on, 0 ≡ "off }, (vii) where e i,k follows a simple two-state Markov process with a-symmetric crossover and same-state probabilities: The model implicitly incorporates three i.i.d noise sources: Gaussian n τ and n α as well as the noise process driving the state changes for e i,k . These sources provide the randomness of the model. The parameter K defines how quickly the speed of path delays can change (for a given variance of n α ). Finally, Δt is the time between instances k − 1 and k. We assume all model parameters (i.e., K, Δt, noise variances, and the "on"/"off " Markov model) to be independent of k (see Figure 5). The hidden channel state vector x k can therefore be represented as τ 1,k , Δτ 2,k , . . . , Δτ Nm,k , α 1,k , . . . , α Nm,k , e i,k , . . . , e Nm,k T .
(27) Note that the model implicitly represents the number of paths as a time-variant parameter. When applied to our particle filtering algorithm, drawing from the proposal density is simple. Each particle stores the above-channel parameters of the model, and then the new state of each particle is drawn randomly from p(x k | x j k−1 ) which corresponds to drawing values for n α and n τ as well as propagating the "on"/"off " Markov model, and then updating the channel parameters for k according to (23)-(26).

Model matching
It is important to point out that a sequential estimator is only as good as its state transition model matches the real world situation. The state model needs to capture all relevant hidden states with memory and needs to correctly model their dependencies, while adhering to the first order Markov condition. Furthermore, any memory of the measurement noise affecting the likelihood function p(z k | x k ) must be explicitly contained as additional states of the model x, so that the measurement noise is i.i.d.
The channel state model is motivated by channel modeling work for multipath prone environments such as the urban satellite navigation channel [21,26]. In fact, the process of constructing a channel model in order to characterize the channel for signal level simulations and receiver evaluation comes close to our task of building a first-order Markov process for sequential estimation. For particle filtering, the model needs to satisfy the condition that one can draw states with relatively low computational complexity. Adapting the model structure and the model parameters to the real channel environment is a task for current and future work. It may even be possible to envisage hierarchical models in which the selection of the current model itself follows a process. In this case, a sequential estimator will automatically choose the correct weighting of these models according to their ability to fit the received signal.

Integration into a receiver
For receiver integration, the computational complexity of the filtering algorithm is crucial. From a theoretical point of view, it is desirable to run the sequential filter clocked corresponding to the coherent integration period of the receiver and with a very large number of particles. From the practical point of view, however, it is desirable to reduce the sequential filter rate to the navigation rate and to minimize the number of particles. Existing ML approaches can help here to achieve a flexible complexity/performance tradeoff, as strategies already developed to extend the observation   periods of ML estimators can be used directly to reduce the rate of the sequential filtering algorithm.

PERFORMANCE EVALUATION
To demonstrate the capabilities of the SMC-based Bayesian time delay estimator proposed in Section 4, we have carried out computer simulations for both static and dynamic channel environments. The signal-to-multipath ratio (SMR) was chosen constant and equal to 6 dB, while the relative phases are changing according to Δϕ i,k = 2πΔτ i,k f c with f c = 1575.42 MHz being the frequency of the L1 carrier. The Bayesian estimator uses a time interval of Δt = 1 millisecond corresponding to the duration of a codeword. The amplitude averaging coefficient is set to Q = 10, and signal compression is applied with N cc = 41 (code-matched correlators) N pc = 25. The channel parameters σ 2 i,τ , σ 2 i,α , K, and p(e i,k | e i,k−1 ) are selected to fit the statistics of a real channel according to [21]. The SIR-PF uses the minimum mean square error (MMSE) criterion to estimate the parameters x k from the posterior.

Static multipath scenario
The capability of multipath mitigation techniques is commonly assessed by showing the systematic error due to a single multipath replica plotted as a function of the relative multipath delay in a static channel scenario. In Figure 6, the root mean square error (RMSE) is shown for the proposed sequential estimator, implemented as a SIR-PF with N s = 2000 particles. For comparison, the performance of conventional DLLs with Narrow Correlator  activity tracking. The performance of a single path estimator is comparable to that of a DLL with infinitesimal correlator spacing and shows a considerable bias over a large delay range. The estimator with fixed two-path model successfully mitigates the multipath bias for delays greater than 30 m. However, for smaller delays it shows an increasing variance and is outperformed by the single path estimator. The estimator with path activity tracking is capable of combining the advantages of both models. From the posterior, it is possible to calculate the estimated average probability P(N m,k = 2 | Z k ) of a two-path model, which is shown in Figure 7, and indicates the transition between the models: for small delays the two paths essentially merge to a single one. Note that in these simulations, the model parameters of the sequential estimator are still the ones designed for the dynamic channel and not optimal for this static scenario. The RMSE as function of the number of particles N s is shown in Figure 8 for a relative multipath delay of 14 m, which corresponds to the worst case in Figure 6.

Dynamic scenario
The dynamic multipath channel scenario with up to N m = 3 paths, used in the following simulations and depicted in Figure 9, has been generated according to the movement model, whereby the parameters K = 25000, σ nα = 10 −10 , σ nτ = 10 −8 , and p onoff = p offon = 0.0001 were chosen to resemble a typical urban satellite navigation channel environment [26]. Figure 10 shows for this multipath scenario at a carrierto-noise ratio of C/N 0 = 45 dB-Hz the BPSK and the BOC(1,1) performance of a DLL with a narrow correlator spacing of 10 −7 seconds, corresponding to a tenth of a code chip. A DLL loop bandwidth of 1 Hz was selected, which appeared to deliver the best DLL performance for this channel. For a fair comparison of the different modulation techniques, both signals were generated with the same C/A code sequence of length 1023, the same receiver bandwidth of 5 MHz (one-sided), a sampling rate of 1/T s = 10.23 MHz, and a coherent integration time of LT s = 1 millisecond. Although the results show an improvement for the BOC(1,1) modulation, there still remains a considerable error due to the multipath. This behavior is confirmed by a comparison with the multipath error envelope, depicted in Figure 11. While in some multipath delay regions, the error is reduced by the BOC(1,1) modulation (see, e.g, the segment marked by the dashed box in Figure 10), more sophisticated multipath mitigating techniques are required to reduce the errors due to short range multipath.
The simulation results for the particle filter-based estimator with N s = 20000 particles are given in Figures 12(a) and  Figure 12: Dynamic scenario: Performance of the sequential estimator with particle filtering (black) compared to the DLL (grey). The sequential estimator is clearly superior for multipath mitigation.
12(b) for the same channel as in Figure 10. They show that the performance could be improved substantially, resulting in a reduction of the root mean square (RMS) error from 3.77 m to 0.769 m for BPSK modulation and from 2.61 m to 0.694 m for BOC(1,1) modulation, respectively.
The RMSE as function of the number of particles N s is shown in Figure 13.

CONCLUSIONS
We have demonstrated how sequential Bayesian estimation techniques can be applied to the multipath mitigation problem in a navigation receiver. The proposed approach is characterized by code-matched, correlator-based signal compression together with interpolation techniques for efficient likelihood computation in combination with a particle filter realization of the prediction and update recursion. The considered movement model has been adapted to dynamic multipath scenarios and incorporates the number of echoes as a time-variant hidden channel state variable that is tracked together with the other parameters in a probabilistic fashion. A further advantage compared to ML estimation is that the posterior PDF at the output of the estimator represents reliability information about the desired parameters and preserves the ambiguities and multiple modes that may occur within the likelihood function. Simulation results for BPSKand BOC-(1,1) modulated signals show that in both cases significant improvements can be achieved compared to a DLL with narrow correlator. In this work, we have employed two methods to reduce complexity: the signal compression to facilitate the computation of the likelihood function as well as a simple form of Rao-Blackwellization to eliminate the complex amplitudes from the state space. Further work will concentrate on additional complexity reduction techniques such as more suitable proposal functions or particle filtering algorithms such as the auxiliary particle filter that are possibly more efficient with respect to the number of particles when applied to our problem domain.