Estimating and detecting random processes on the unit circle

Consider a detection system with two hypothesis:

where $Z_{t}$ is the noisy measurements and $X_{t}=e^{i\Theta_{t}}\in S^{1}$ is the circle-valued signal with $\Theta_{t}$ obeying an arbitrary random process. Measurement function $h(X_{t})$ in our case, is given by

with $\Re[\cdot]$ denoting the real part. $\{W_{t},0\leq t\leq T\}$ is the standard Wiener process, which represents the measurement noise and is independent of $X_{t}$. $\sigma_{0}$ denotes the variance of the measurement noise and is assumed known. In this paper, we consider three specific random processes for $\Theta_{t}$:

where $q_{\Theta}$ and $q_{w}$ are known functions of phase, $\Theta$, and frequency, $w$. They represent the noise variance for phase and frequency, respectively. $w_{0}$ is the known initial frequency and $\phi_{0}$ is the random initial phase, uniformly distributed across $[0,2\pi)$. $\{B_{t},0\leq t\leq T\}$ and $\{\tilde{B}_{t},0\leq t\leq T\}$ are two independent standard Wiener processes, which are also independent of $X_{t}$ and $\{W_{t},0\leq t\leq T\}$.

The Neyman-Pearson (EC) detector is the optimal detector in the sense that it maximizes the detection probability (Pd) with given false alarm probability (Pf). In order to form this type of detector, the log likelihood ratio $\log\Lambda_{t}$ has to be evaluated, which is given by

with $\mathcal{F}_{t}$ the filtration of $\sigma$-algebras generated by $\{Z_{\tau},0\leq\tau\leq t\}$, given in (2) with $\mathcal{F}_{s}\subseteq\mathcal{F}_{t}$, $\forall s\leq t$.

From (7) and (8), the EC detector is formed in two steps: 1.) calculate the conditional mean estimate $\hat{X}_{t}$; 2.) form the log likelihood ratio $\log\Lambda_{T}$ at the termination time $T$ and compare it with a pre-specified threshold to claim $H_{0}$ or $H_{1}$ correspondingly.

Given the filtered probability space $(\Omega,\mathcal{F},\mathbb{P};\{\mathcal{F}_{t}\})$ with $\Omega$ the underlying sample space and $\mathbb{P}$ the probability measure defined on $\mathcal{F}$, suppose the signal and measurement processes are given by

with $x_{t}\in\mathbb{R}^{d}$ and $b:\mathbb{R}^{d}\to\mathbb{R}^{d}$, $\sigma:\mathbb{R}^{d}\to\mathbb{R}^{d\times p}$, $h:\mathbb{R}^{d}\to\mathbb{R}$, $a\overset{\triangle}{=}\frac{1}{2}\sigma\sigma^{\rm{T}}$ with superscript $\rm{T}$ being the transpose of a matrix. The filtering problem can be formulated as: for a general measurable function $g(x_{t})$, calculate $\hat{g}(x_{t})$ by

where $p(x,t|\mathcal{F}_{t})$ denotes the conditional density of $x$ at time $t$. The recursive formula for evolving $p(x,t|\mathcal{F}_{t})$ is described by the Kushner equation (Kushner (1967))

with $\mathcal{L}^{*}$ the adjoint of the infinitesimal generator $\mathcal{L}$, given respectively by

To uniquely characterize the conditional density, we pair (11) with an arbitrary test function $g(x)$ with compact support, where the duality pairing is denoted by $\langle\cdot,\cdot\rangle$. We then have

which further implies

here $\mathcal{L}$ is defined in (13). Both $\hat{}$ and $\langle\cdot\rangle$ denote the conditional mean with respect to the conditional density $p(x,t|\mathcal{F}_{t})$. $\hat{g}(x_{t})\overset{\triangle}{=}\Big{\langle}p(x,t|\mathcal{F}_{t}),g(x)\Big{\rangle}$ is defined in (10).

From the previous section, we observe the duality relationship between the conditional density and the estimated test function. In order to fully characterize $p(x,t|\mathcal{F}_{t})$, an infinite number of statistics (or test functions $g(\cdot)$) are required. In this paper, we set the test functions to be the moment functions of the signal $X_{t}$, given as $X_{n}(t)\overset{\triangle}{=}e^{in\Theta_{t}}$ for $n=0\cdots,N-1$, considering firstly, $X_{n}(t)$ can be regarded as the characteristic function of $\Theta$, which completely determines the properties of the probability distribution of $\Theta$. Secondly, $X_{n}(t)$ spans the eigenspace of the generator $\mathcal{L}$, as discussed in the following section. In doing so, solving an infinite dimensional Kushner equation boils down to solving $N$’s ODEs.

As we can see from (15), the filtering process is composed of two terms: the one-step prediction (the first term) and the update (the second term). We will derive the solutions for $\hat{X}_{n}$ in terms of these two properties.

In order to obtain the infinitesimal generator of the signal process, we need to write the signal dynamics in Ito differential form, for three scenarios. In particular, we have for scenario (i) and (ii)

Comparing (17) with (16), a constant rotation term $inw_{0}$ is added to the drift coefficient. Considering the martingale process in both circumstances: $inq_{\Theta}^{1/2}X_{n}dB_{t}$, with the property $\mathbf{E}[\int_{s}^{t}inq_{\Theta}^{1/2}X_{n}(\tau)dB_{\tau}|\mathcal{F}_{s}]=0$, we have a separate infinitesimal generator for each $X_{n}$, where $\mathcal{L}_{n}=inw_{0}-\frac{q_{\Theta}n^{2}}{2}$, with $w_{0}=0$ in scenario (i).

The differential equation for $X_{n}$ in scenario (iii) can be obtained in the same manner

However, as shown here, the randomness appears in the drift coefficient. Now it is difficult to interpret the generator by simply writing $\mathcal{L}_{n}(w)=inw(t)-\frac{q_{\Theta}n^{2}}{2}$ with $w(t)$ a random process. Instead, $\mathcal{L}_{n}(w)$ should be interpreted in the mean sense by marginalizing out the random variable $w(t)$. In other words, we have the infinitesimal generator computed by $\mathcal{L}_{n}=\mathbf{E}_{w}[\mathcal{L}_{n}(w)|\mathcal{F}_{t}]=\int_{w\in\mathbb{R}}\mathcal{L}_{n}(w)p(w,t|\mathcal{F}_{t})dw$, with $d\hat{X}_{n}-\mathcal{L}_{n}\hat{X}_{n}dt$ still a martingale process. Alternatively, we introduce the new state variables as $X_{mn}\overset{\triangle}{=}w^{m}e^{in\Theta}$, with $m=0,\dots,M-1$ and $n=0,\dots N-1$, where superscript $m$ and $n$ denote the $m$th and $n$th power, respectively. This is an “augmented” state space by noting that $X_{n}=X_{mn}|_{m=0}$. For predicting $X_{mn}$, likewise, we write down the differential equation for $X_{mn}$ as

with $\mathcal{L}_{mn}(X_{mn})=\frac{q_{w}m(m-1)}{2}X_{m-2,n}+inX_{m+1,n}-\frac{q_{\Theta}n^{2}}{2}X_{mn}$ for each $X_{mn}$. Here by augmenting the state space from $X_{n}$ to $X_{mn}$, the time-variant generator is replaced by a time-invariant generator at the cost of higher computational complexity.

In all three scenarios, we have measurement function $h(X_{t})=\Re[X_{t}]$ from (3). This is a considerable advantage since $h$ can be regarded as a linear operator in both the subspace $\{X_{n}\}_{n=0}^{N-1}$ and $\{X_{mn}\}_{m=0,n=1}^{M-1,N-1}$ by noting

with $()^{*}$ denoting the complex conjugation. Combining (16), (17), (19) with (20), (21) and substituting into (15), we can compute $\hat{X}_{n}(t)$ for (i)&(ii) and $\hat{X}_{mn}(t)$ for (iii) recursively by solving a sequence of stochastic ODEs.

The log-likelihood ratio, defined in (7) can then be updated recursively as well by substituting in $\hat{X}_{t}$ from (15) without too much effort, considering $\hat{X}_{t}=\hat{X}_{n}(t)|_{n=1}$ for (i)&(ii) and $\hat{X}_{t}=\hat{X}_{mn}(t)|_{m=0,\;n=1}$ for (iii).

Since scenario (i) is analogous to (ii) with $w_{0}=0$, we only do experiments for (ii) and (iii). The Neyman-Pearson (EC) detector is composed of a conditional mean estimator and a log-likelihood ratio detector, hence we need to check the performance of both estimation and detection. For estimation, we generate the synthetic signal $X_{\rm syn}=A_{\rm syn}e^{i\Theta_{\rm syn}}$ with $\Theta_{\rm syn}$ given in (5) and (6). Noisy measurements are generated from (2). For detection, we randomly generate measurements from (1) or (2), where signal $X_{\rm syn}=A_{\rm syn}e^{i\Theta_{\rm syn}}$ is generated as above. We fix termination time $T$ and compare $\log\Lambda_{T}$ with a pre-specified threshold and claim detection if $\log\Lambda_{T}$ is greater than the threshold, and vice versa. We simulate both the signal process and the measurement process numerically using the Euler-Maruyama method with sampling interval $\Delta t=0.1$s. The parameters are listed in Table 1.

Two realizations of EC estimated $\Re[\hat{X}_{t}]$, with phase dynamics described in (ii) and (iii) are shown in Fig. 1 and Fig. 2, compared with EKF estimated $\Re[\hat{X}_{t}]$. From both plots, with SNR = 0.1, the signal is submerged in large noise. However, the EC estimated signal (blue) is closer to the synthetic signal (yellow) compared with EKF estimated signal (red). When frequency experiences small wandering in scenario (iii), as displayed in Fig. 2, by estimating the product $X_{mn}$, we can still extract $\hat{X}_{t}$ with higher accuracy than EKF estimated signal.

To quantify the performance of the approximate EC detector, receiver operating characteristic (ROC) curves with different SNRs under scenario (ii) with $q_{\Theta}=10^{-1}$ and $w_{0}=0.012$Hz are displayed in Fig. 3, compared with the EKF detector as well. In Fig. 4, ROCs of the EC and EKF detector under scenario (iii) with $q_{\Theta}=10^{-1}$, $q_{w}=10^{-8}$ and $w_{0}=0.012$Hz are plotted. From Fig. 3, as SNR drops from 0.1 (top) to 0.0816 (middle) and then to 0.0577 (bottom), Pd of the EC detector (blue) descends from 0.9 to 0.6, then to 0.4 at Pf $=10^{-2}$. Fig. 4 exhibits a similar trend, with higher Pd as SNR gets larger. Comparing two top panels between Fig. 3 and Fig. 4, when the frequency is slowly wandering, with $q_{w}=10^{-8}$ in Fig. 4, EC experiences a slight degradation, with Pd drops from 0.9 to 0.8 at SNR = 0.1 and Pf = $10^{-2}$. Throughout all six panels, the EC detector advantages over the EKF detector with higher Pd across the Pf range, especially when Pf $<10^{-2}$.