# On the Accuracy of Phase Extraction from a Known-Frequency Noisy Sinusoidal Signal

Emmanuel Dervieux, Florian Tilquin, Alexis Bisiaux, Wilfried Uhring

► To cite this version:

Emmanuel Dervieux, Florian Tilquin, Alexis Bisiaux, Wilfried Uhring. On the Accuracy of Phase Extraction from a Known-Frequency Noisy Sinusoidal Signal. APSIPA Transactions on Signal and Information Processing, 2024, 13 (1), pp.e35. 10.1561/116.20240042. hal-04984002

HAL Id: hal-04984002

https://hal.science/hal-04984002v1

Submitted on 10 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Original Paper

# On the Accuracy of Phase Extraction
from a Known-Frequency Noisy
Sinusoidal Signal

Emmanuel Dervieux¹,²*, Florian Tilquin¹, Alexis Bisiaux¹ and
Wilfried Uhring²

¹Biosency, Cesson-Sévigné, France

²ICube, University of Strasbourg and CNRS, Strasbourg, France

ABSTRACT

Accurate phase extraction from sinusoidal signals is a crucial task
in various signal processing applications. While prior research
predominantly addresses the case of asynchronous sampling with
unknown signal frequency, this study focuses on the more specific
situation where synchronous sampling is possible, and the signal's
frequency is known. In this framework, a comprehensive analysis
of phase estimation accuracy in the presence of both additive and
phase noises is presented. A closed-form expression for the asymp-
totic Probability Density Function (PDF) of the resulting phase
estimator is validated by simulations depicting Root Mean Square
Error (RMSE) trends in different noise scenarios. This estimator
is asymptotically efficient, converging rapidly to its Cramèr-Rao
Lower Bound (CRLB). Three distinct RMSE behaviours were
identified based on SNR, sample count (N), and noise level: (i)
saturation towards a random guess at low Signal to Noise Ratio
(SNR), (ii) linear decrease with the square roots of N and SNR at
moderate noise levels, and (iii) saturation at high SNR towards
a noise floor dependent on the phase noise level. By quantifying
the impact of sample count, additive noise, and phase noise on

<footnote>*Corresponding author: emmanuel.dervieux@biosency.com.</footnote>
phase estimation accuracy, this work provides valuable insights for designing systems requiring precise phase extraction, such as phase-based fluorescence assays or system identification.

**Keywords:** Spectral estimation, phase estimation, synchronous sampling, Fourier analysis

# 1 Introduction

Spectral estimation plays a critical role in signal processing by characterising a signal's spectral attributes, including amplitudes and phase shifts. This topic has gathered considerable research interest for decades due to its numerous applications in various fields, such as telecommunications, radar, seismology, and power grid analysis [23, 29, 33, 24]. In the general case, the frequencies of interest $f_i$ of the signal under study are a priori unknown. Thus, it is exceedingly unlikely that given a sampling frequency $f_s$ and a sampling length $N$, the numbers $f_i \cdot N/f_s$ are integers. This condition—known as “asynchronous sampling”—leads to the infamous picket fence and spectral leaking effects [12], which may be mitigated by an appropriate windowing function choice [30], the use of all-phase Discrete Fourier Transform (DFT) [38, 14], or both [34], for example.

There are certain cases, however, for which the signal under study is purely sinusoidal with a known frequency. This scenario arises when characterising linear systems, which may be fed a sinusoidal excitation signal of known frequency, amplitude, and phase, while recording their output. The analysis of the attenuation and phase shift induced by the system at hand can then yield useful information. For instance, in the context of frequency-based Dual Lifetime Referencing (f-DLR) [17], the phase shift between a fluorescence excitation signal of known frequency and the re-emitted one can be used to accurately measure the concentration of a variety of analytes [37, 1, 32, 35]. In this situation—known as “synchronous sampling”—the number of samples taken, as well as the sampling and excitation frequencies $f$ and $f_s$, can be chosen so that $f \cdot N/f_s$ is an integer, which suppresses the above-mentioned deleterious effects [12, 10].

Yet, as far as we are aware, no comprehensive study has been conducted to characterise the achievable accuracy of phase estimation in such a synchronous sampling scenario. In this paper, we present theoretical developments leading to a closed-form expression of the asymptotic Probability Density Function (PDF) of the phase estimate of a noisy sinusoidal signal in the presence of both phase and additive noises. The presented derivations are supported by simulations results, showing the resulting phase Root Mean Square Error
(RMSE) at different noise levels. We then show that the derived phase estimator is asymptotically efficient with a fast convergence. Finally, we discuss its asymptotic behaviour in the case of very high or very low noise levels and sample numbers.

## 2 Problem Formulation

In the remainder of this document, the objective is always to retrieve the phase $\varphi$ of a real discrete signal of length $N$ defined as

$$s_n = A_s \cdot \cos \left( \frac{2 \cdot \pi \cdot f_0 \cdot n}{f_s} + \varphi + p_n \right) + x_n, \text{ with } \begin{cases} x_n \stackrel{iid}{\sim} \mathcal{N}(0, \sigma_x^2) \\ p_n \stackrel{iid}{\sim} \mathcal{N}(0, \sigma_p^2) \\ n \in [0; N-1] \end{cases} \quad (1)$$

with $f_0$ the frequency of the signal itself, $f_s$ its sampling frequency—always chosen such that $f_s > 2 \cdot f_0$, the Nyquist frequency—and $A_s$ its amplitude. The $x_n$ and $p_n$ random variables—of variances $\sigma_x^2$ and $\sigma_p^2$—represent additive measurement noise and sampling-induced phase noise, respectively. Of note, it is also considered that $A_s$ and $\varphi$, though unknown, remain constant throughout the acquisition duration $N/f_s$. A representative illustration of the issue at hand, involving most of the parameters introduced above, may be seen in Figure 1.

Typically, in an f-DLR sensing scheme, $A_s$ and $\varphi$ would correspond to:

1. the intensity of the collected light: a function of the quantum yield of the involved fluorophores, of their concentrations, and of the illumination and light collection parameters, and

2. the phase shift: function of the ratio of the different fluorophores species, conveying the concentration of the analyte of interest.

Hence, it is of particular importance to accurately estimate $\varphi$, and to characterise the influence of $\sigma_p$, $\sigma_x$, and $N$ on its RMSE, since it will directly influence the reachable accuracy on the measurement of a given analyte's concentration.

In the remainder of this article, we adopt the following notations: $x_n$ refers to the n-th element of a given vector $X$, $\square^T$ is the transpose operator, $\mathbb{0}_N$ and $\mathbb{1}_N$ stand for the zero and unit vector in $\mathbb{R}^N$, respectively, $\mathbb{R}(z)$ and $\mathbb{F}(z)$ stand for the real and imaginary parts of a given complex number $z$, while $|z|$ and $\arg(z)$ stand for its modulus and argument. $\triangleq$ means “per definition”, $\square$ is the complex conjugate operator, $\mathcal{N}$ and $\mathcal{CN}$ stand for the normal and complex normal distributions, respectively, and $x \perp y$ denotes
<figure><img src="image_2.png" /><figcaption>Figure 1: Temporal representation of the problem: the objective is to estimate the phase $\varphi$ of an idealised sinusoidal signal of amplitude $A_s$ (●) from a noisy measurement of the latter (↔). The noise originates from two distinct sources: (i) an additive noise $x_n$ which makes the sampled signal depart from its ideal counterpart, and (ii) a phase noise $p_n$, which randomly shifts the sampling times (···) in comparison to an ideal sampling at frequency $f_s$ (■). The symbol $\propto$ denotes proportionality, not to be confused with the Greek letter $\alpha$ (α), used later on in this paper.</figcaption></figure>

the independence between two random variable $x$ and $y$. Finally, the Signal to Noise Ratio (SNR) of the measurement is defined as

$$ \text{SNR} = \frac{A_s^2}{2 \cdot \sigma_x^2} \quad \text{and} \quad \text{SNR}_{\text{dB}} = 10 \cdot \log_{10}(\text{SNR}) \qquad (2) $$

### 3 Characterisation of the DFT Distribution

This paper focuses on $\varphi$ estimation through the study of the DFT of the above-presented noisy signal. Indeed, we demonstrate in Section 4 that an unbiased and efficient estimator of $\varphi$—denoted as $\hat{\varphi}$—can be derived by taking the argument of the signal's DFT at frequency $f_0$. In order to derive the PDF of $\hat{\varphi}$, the PDF of this DFT must thus be known first.

To this end, let us first consider the k-th index of the N-points DFT of the above-mentioned signal [7], i.e. its DFT at frequency $f_0$:

$$ D_N = \sum_{n=0}^{N-1} s_n \cdot e^{-i \cdot \frac{2 \cdot \pi \cdot k \cdot n}{N}} \qquad (3) $$
with $k$, $N$, $f_0$ and $f_s$ chosen such that $\frac{k}{N} = \frac{f_0}{f_s}$ (synchronous sampling hypothesis). Let $\alpha = \frac{2\cdot\pi\cdot f_0}{f_s}$ and

$$S = \begin{pmatrix} s_0 \\ s_1 \\ \vdots \\ s_{N-1} \end{pmatrix}, X = \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{pmatrix}, Z = \begin{pmatrix} e^{i \cdot p_0} \\ e^{i \cdot p_1} \\ \vdots \\ e^{i \cdot p_{N-1}} \end{pmatrix}, \mathcal{U}_{\alpha} = \begin{pmatrix} 1 \\ e^{-i \cdot \alpha} \\ \vdots \\ e^{-i \cdot \alpha \cdot (N-1)} \end{pmatrix} \quad (4)$$

Then $D_N$ may be rewritten as

$$D_N = \sum_{n=0}^{N-1} s_n \cdot e^{-i \cdot \alpha \cdot n} = \mathcal{U}_{\alpha}^{\intercal} \cdot S \quad (5)$$

Let also define $\tilde{S}$ as $S = \tilde{S} + X$. Using Euler's formula, each element $\tilde{s}_n$ of $\tilde{S}$ may then be expressed as

$$\begin{aligned} \tilde{s}_n &= A_s \cdot \frac{e^{i \cdot \left( \frac{2 \cdot \pi \cdot f_0 \cdot n}{f_s} + \varphi + p_n \right)} + e^{-i \cdot \left( \frac{2 \cdot \pi \cdot f_0 \cdot n}{f_s} + \varphi + p_n \right)}}{2} \\ &= \frac{A_s}{2} \cdot \left( e^{i \cdot (\alpha \cdot n + \varphi + p_n)} + e^{-i \cdot (\alpha \cdot n + \varphi + p_n)} \right) \end{aligned} \quad (6)$$

The remainder of this section is organised as follows: the expected value and variance of $D_N$ are computed in Sections 3.1 and 3.2, respectively, while its asymptotic PDF is derived in Section 3.3.

## 3.1 Expected Value of $D_N$

The expected value of $D_N$ is given by

$$\mathbf{E}[D_N] = \mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}] + \mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot X] \quad (7)$$

Regarding $\mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}]$,

$$\begin{aligned} \mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S} &= \frac{A_s}{2} \left( e^{i \cdot \varphi} \sum_{n=0}^{N-1} e^{i \cdot p_n} + e^{-i \cdot \varphi} \sum_{n=0}^{N-1} e^{-2 \cdot i \cdot \alpha \cdot n} e^{-i \cdot p_n} \right) \\ &= \frac{A_s}{2} \cdot \left( e^{i \cdot \varphi} \cdot \mathbb{1}_N^{\intercal} \cdot Z + e^{-i \cdot \varphi} \cdot (\mathcal{U}_{\alpha}^2)^{\intercal} \cdot \bar{Z} \right) \end{aligned} \quad (8)$$

hence

$$\mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}] = \frac{A_s}{2} \cdot \left( e^{i \cdot \varphi} \cdot \mathbb{1}_N^{\intercal} \cdot \mathbf{E}[Z] + e^{-i \cdot \varphi} \cdot (\mathcal{U}_{\alpha}^2)^{\intercal} \cdot \mathbf{E}[\bar{Z}] \right) \quad (9)$$
We thus need to compute $\mathbf{E}[Z]$. For a given $p$ amongst $p_n$, let us consider $\mathbf{E}[e^{i \cdot p}]$ first, and let $[g : v \mapsto e^{i \cdot v}]$. Then

$$ \mathbf{E}[e^{i \cdot p}] = \mathbf{E}[g(p)] = \int_{-\infty}^{+\infty} g(v) \cdot f_p(v) dv \quad (10) $$

thanks to the Law Of The Unconscious Statistician (LOTUS) [31], wherein $f_p$ is the probability density function of $p$—namely a centred normal distribution of variance $\sigma_p^2$, see Equation 1—given by

$$ \left[ f_p : v \mapsto \frac{1}{\sigma_p \cdot \sqrt{2 \cdot \pi}} \cdot e^{-\frac{v^2}{2 \cdot \sigma_p^2}} \right] \quad (11) $$

We thus have

$$ \begin{aligned} \mathbf{E}[e^{i \cdot p}] &= \int_{-\infty}^{+\infty} \frac{1}{\sigma_p \cdot \sqrt{2 \cdot \pi}} \cdot e^{-\frac{v^2}{2 \cdot \sigma_p^2} + i \cdot v} dv \\ &= \frac{e^{-\frac{\sigma_p^2}{2}}}{\sigma_p \cdot \sqrt{2 \cdot \pi}} \int_{-\infty}^{+\infty} e^{\left( \frac{i \cdot v}{\sigma_p \cdot \sqrt{2}} + \frac{\sigma_p}{\sqrt{2}} \right)^2} dv \end{aligned} \quad (12) $$

Using the change of variable $u = h(v)$ with

$$ \left[ h : v \mapsto \frac{i \cdot v}{\sigma_p \cdot \sqrt{2}} + \frac{\sigma_p}{\sqrt{2}} \right] \quad (13) $$

we thus have $du = \frac{i \cdot dv}{\sigma_p \cdot \sqrt{2}}$, with erfi being the imaginary error function,

$$ \begin{aligned} \mathbf{E}[e^{i \cdot p}] &= \frac{e^{-\frac{\sigma_p^2}{2}}}{\sigma_p \cdot \sqrt{2 \cdot \pi}} \cdot \lim_{v \to +\infty} \int_{-h(v)}^{+h(v)} e^{u^2 \cdot \frac{\sigma_p \cdot \sqrt{2}}{i}} du \\ &= \frac{e^{-\frac{\sigma_p^2}{2}}}{2 \cdot i} \cdot \lim_{v \to +\infty} \underbrace{\frac{2}{\sqrt{\pi}} \int_{-h(v)}^{+h(v)} e^{u^2} du}_{= [\text{erfi}(u)]_{-h(v)}^{+h(v)}} \\ &= \underbrace{2 \cdot i}_{=2 \cdot i} \end{aligned} \quad (14) $$

Back to $\mathbf{E}[\mathcal{U}_{\alpha}^{I} \cdot \tilde{S}]$, since

$$ \mathbf{E}[\overline{Z}] = \overline{\mathbf{E}[Z]} = \mathbf{E}[Z] = 1_N \cdot e^{-\frac{\sigma_p^2}{2}} \quad (15) $$

we thus have

$$ \mathbf{E}[\mathcal{U}_{\alpha}^{I} \cdot \tilde{S}] = \frac{A_s}{2} \cdot \left( e^{i \cdot \varphi} \cdot N \cdot e^{-\frac{\sigma_p^2}{2}} + e^{-i \cdot \varphi} \cdot e^{-\frac{\sigma_p^2}{2}} \cdot \sum_{n=0}^{N-1} u_{\alpha,n}^2 \right) \quad (16) $$
Noting that $\sum_{n=0}^{N-1} u_{\alpha,n}^2 = 0$ because $e^{-2 \cdot i \cdot \alpha}$ is an N-th root of unity, comes

$$ \mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}] = \frac{A_s \cdot N}{2} \cdot e^{i \cdot \varphi} \cdot e^{-\frac{\sigma_p^2}{2}} \quad (17) $$

Regarding $\mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot X]$,

$$ \mathbf{E}[\mathcal{U}_{\alpha}^{\intercal} \cdot X] = \mathcal{U}_{\alpha}^{\intercal} \cdot \underbrace{\mathbf{E}[X]}_{=0} = 0 \quad (18) $$

Finally, the expected value of $D_N$ comes to be

$$ \mathbf{E}[D_N] = \frac{A_s \cdot N}{2} \cdot e^{i \cdot \varphi} \cdot e^{-\frac{\sigma_p^2}{2}} \quad (19) $$

## 3.2 Variance of $D_N$

The variance of $D_N$ may also be calculated in a similar manner, starting with

$$ \mathrm{Var}(D_N) = \mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}) + \mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot X) \quad (20) $$

Regarding $\mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S})$,

$$ \begin{aligned} \mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}) &= \mathrm{Var} \left[ \frac{A_s}{2} \cdot \left( e^{i \cdot \varphi} \cdot \mathbb{1}_N^{\intercal} \cdot Z + e^{-i \cdot \varphi} \cdot (\mathcal{U}_{\alpha}^2)^{\intercal} \cdot \overline{Z} \right) \right] \\ &= \frac{A_s^2}{4} \cdot \left( \sum_{n=0}^{N-1} \mathrm{Var}(z_n) + \mathrm{Var}(\overline{z_n}) \right) \end{aligned} \quad (21) $$

We thus need to compute $\mathrm{Var}(z_n)$, i.e. $\mathrm{Var}(e^{i \cdot p_n})$.

$$ \begin{aligned} \mathrm{Var}(e^{i \cdot p}) &= \mathbf{E} \left[ |e^{i \cdot p}|^2 \right] - \left| \mathbf{E} [e^{i \cdot p}] \right|^2 \\ &= 1 - e^{-\sigma_p^2} = \mathrm{Var}(e^{-i \cdot p}) \end{aligned} \quad (22) $$

as $\mathbf{E}[e^{-i \cdot p}] = \overline{\mathbf{E}[e^{i \cdot p}]}$. Thus

$$ \mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot \tilde{S}) = \frac{N \cdot A_s^2}{2} \cdot (1 - e^{-\sigma_p^2}) \quad (23) $$

Regarding $\mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot X)$,

$$ \begin{aligned} \mathrm{Var}(\mathcal{U}_{\alpha}^{\intercal} \cdot X) &= \sum_{n=0}^{N-1} \mathrm{Var}(u_{\alpha,n} x_n) = \sum_{n=0}^{N-1} \left( \mathbf{E} \left[ |u_{\alpha,n} x_n|^2 \right] - \left| \mathbf{E} [u_{\alpha,n} x_n] \right|^2 \right) \quad (24) \\ &= \sum_{n=0}^{N-1} \mathbf{E}[x_n^2] = N \cdot \sigma_x^2 \end{aligned} $$
Finally, the variance of $D_N$ comes to be

$$ \mathrm{Var}(D_N) = N \cdot \left( \frac{A_s^2}{2} \cdot (1 - e^{-\sigma_p^2}) + \sigma_x^2 \right) \quad (25) $$

### 3.3 Distribution of $D_N$

Now that the expected value and variance of $D_N$ are known, the next step is to study its PDF. To do so, we focus on a reduced version of $D_N—\tilde{D}_N—$defined as $\tilde{D}_N = \frac{2 \cdot D_N}{A_s \cdot N}$, with

$$ \beta_p = e^{-\frac{\sigma_p^2}{2}} \qquad \begin{aligned} \mathbf{E}[\tilde{D}_N] &= \beta_p \cdot e^{i \cdot \varphi} \\ \mathbf{Var}(\tilde{D}_N) &= \frac{2}{N} \cdot \left(1 - \beta_p^2 + \frac{1}{\text{SNR}}\right) \end{aligned} \quad (26) $$

We then proceed in two steps: at first, the convergence in law of $\Re(\tilde{D}_N)$ and $\Im(\tilde{D}_N)$ towards normal distributions is demonstrated. Then, the asymptotic independence of the latter two quantities is shown. These two demonstrations establish that $\tilde{D}_N$ converges in law<fnref n="1" /> towards a complex normal distribution [18, pp. 540–559], which is a crucial requirement for the forthcoming developments (see Section 4). However, before delving any deeper into this two-step demonstration, we can further simplify the issue at hand, observing that

$$ \tilde{D}_N = \underbrace{\frac{2 \cdot \mathcal{U}_\alpha^\intercal \cdot \tilde{S}}{A_s \cdot N}}_{\tilde{D}_s} + \underbrace{\frac{2 \cdot \mathcal{U}_\alpha^\intercal \cdot X}{A_s \cdot N}}_{\tilde{D}_x} \quad (27) $$

Since $X$ is stationary and ergodic, it readily follows that $\tilde{D}_x$ converges in distribution toward a complex normal distribution [27, 4]. Since $\tilde{D}_s$ and $\tilde{D}_x$ are independent, the two above-mentioned steps thus only have to be performed for $\tilde{D}_s$.

#### 3.3.1 Convergence in Law Towards a Normal Distribution

Let us consider the real part of $\tilde{D}_s$

$$ \begin{aligned} \Re(\tilde{D}_s) &= \frac{2}{N} \cdot \sum_{n=0}^{N-1} \cos(\alpha \cdot n) \cdot \cos(\alpha \cdot n + \varphi + p_n) \\ &= \frac{2}{N} \cdot \sum_{n=0}^{N-1} V_n, \quad \text{with } V_n = \cos(\alpha \cdot n) \cdot \cos(\alpha \cdot n + \varphi + p_n) \end{aligned} \quad (28) $$

<footnote n="1">¹Of note, convergence in law is sometimes also referred to as “convergence in distribution” or “weak convergence” [36, p. 18].</footnote>
We will use Lyapunov's Central Limit Theorem (L-CLT) to demonstrate the convergence of $\sum V_n$ towards a normal distribution. To do so, we will first show that there exists a positive $\delta$ such that

$$ \lim_{N \to +\infty} \underbrace{\frac{1}{s_N^{2+\delta}} \sum_{n=0}^{N-1} \mathbf{E} \left[ |V_n - \mathbf{E}[V_n]|^{2+\delta} \right]}_{\gamma_N} = 0 \quad (29) $$

wherein $s_N^2 = \sum_{n=0}^{N-1} \mathbf{Var}(V_n)$, and $\mathbf{Var}(V_n) = \cos^2(\alpha \cdot n) \underbrace{\mathbf{Var}(\cos(t_n + p_n))}_{=\mathbf{E}(\dots)^2 - \mathbf{E}(\dots)^2}$, with $t_n = \alpha \cdot n + \varphi$. Where

$$ \mathbf{E}(\dots)^2 = \frac{1}{2} + \frac{e^{-2 \cdot \sigma_p^2}}{2} \cdot \cos(2 \cdot t_n), \quad \text{and} \quad \mathbf{E}(\dots)^2 = e^{-\sigma_p^2} \cdot \frac{1 + \cos(2 \cdot t_n)}{2} \quad (30) $$

and thus

$$ \begin{aligned} \mathbf{Var}(V_n) &= \cos^2(\alpha \cdot n) \cdot \left[ \frac{1}{2} + e^{-\sigma_p^2} \cdot \left( \frac{\cos(2 \cdot t_n)}{2} \cdot e^{-\sigma_p^2} - \frac{1 + \cos(2 \cdot t_n)}{2} \right) \right] \\ &= \frac{\cos^2(\alpha \cdot n)}{2} \cdot (1 - e^{-\sigma_p^2}) \cdot (1 - e^{-\sigma_p^2} \cdot \cos(2 \cdot t_n)) \\ &\geq \frac{\cos^2(\alpha \cdot n)}{2} \cdot (1 - e^{-\sigma_p^2})^2 \end{aligned} \quad (31) $$

Back to $s_n$,

$$ s_n^2 \geq \frac{(1 - e^{-\sigma_p^2})^2}{2} \underbrace{\sum_{n=0}^{N-1} \cos^2(\alpha \cdot n)}_{=N/2^{(\frac{\delta}{2})}} \geq \frac{N \cdot (1 - e^{-\sigma_p^2})^2}{4} \quad (32) $$

wherein $(\frac{\delta}{2})$ comes from the facts that $\cos^2 x = \frac{1+\cos(2 \cdot x)}{2}$, and that $e^{-2 \cdot i \cdot \alpha}$ is an N-th root of unity (see Equations 16-17). Then, $\forall \delta > 0$

$$ \gamma_N = \frac{1}{s_N^{2+\delta}} \sum_{n=0}^{N-1} \underbrace{\mathbf{E} \left[ |V_n - \mathbf{E}[V_n]|^{2+\delta} \right]}_{\le 2^{2+\delta}} \le \frac{N \cdot 2^{2+\delta}}{\left( \frac{N}{4} \cdot (1 - e^{-\sigma_p^2})^2 \right)^{1+\frac{\delta}{2}}} \xrightarrow{N \to +\infty} 0 \quad (33) $$

Since $V_n$ are independent and of finite variance, according to L-CLT, we thus have

$$ \frac{1}{s_N} \sum_{n=0}^{N-1} (V_n - \mathbf{E}[V_n]) \xrightarrow{N \to +\infty} \mathcal{N}(0, 1) \quad (34) $$

wherein $d \to \infty$ denotes convergence in distribution.
Hence, since $\Re(\tilde{D}_s) = \frac{2}{N} \cdot \sum_{n=0}^{N-1} V_n$,

$$ \Re(\tilde{D}_s) \xrightarrow[N \to +\infty]{d} \mathcal{N} \left( \frac{2}{N} \cdot \sum_{n=0}^{N-1} \mathbf{E}[V_n], \left( \frac{2 \cdot s_N}{N} \right)^2 \right) \quad (35) $$

A similar train of thought can be followed to also demonstrate the asymptotic normality of $\Im(\tilde{D}_s)$.

### 3.3.2 Asymptotic Independence

Demonstrating the complex normality of $\tilde{D}_s$ then only requires to demonstrate that $\Re(\tilde{D}_s) \perp \Im(\tilde{D}_s)$. To do so, it suffices to show that [28, Th. 4.5.1]:

(i) $\Re(\tilde{D}_s)$ and $\Im(\tilde{D}_s)$ follow a bivariate normal distribution, and that

$$ (ii) \operatorname{Cov}(\Re(\tilde{D}_s), \Im(\tilde{D}_s)) = 0. $$

#### Bivariate normality

$$ \forall (a, b) \in \mathbb{R}^2 \text{ let} $$

$$ \begin{aligned} T &= a \cdot \Re(\tilde{D}_s) + b \cdot \Im(\tilde{D}_s) \\ &= \sum_{n=0}^{N-1} \frac{a \cdot \cos(\alpha \cdot n) + b \cdot \sin(\alpha \cdot n)}{N} \cdot \cos(\alpha \cdot n + \varphi + p_n) \end{aligned} \quad (36) $$

It can be shown—as was done in the previous section with $\Re(\tilde{D}_s)$—that $T$ also converges in law towards a normal distribution. Thus, by definition, $\Re(\tilde{D}_s)$ and $\Im(\tilde{D}_s)$ follow a bivariate normal distribution.

#### Covariance

Let

$$ C = \operatorname{Cov}(\Re(\tilde{D}_s), \Im(\tilde{D}_s)) = \underbrace{\mathbf{E}[\Re(\tilde{D}_s) \cdot \Im(\tilde{D}_s)]}_{C_L} - \underbrace{\mathbf{E}[\Re(\tilde{D}_s)] \cdot \mathbf{E}[\Im(\tilde{D}_s)]}_{C_R} \quad (37) $$

and $\forall n, \Phi_n = \varphi + p_n$. Then

$$ \begin{aligned} C_L &= \mathbf{E} \left[ \left( \sum_{l=0}^{N-1} \frac{\cos(\alpha \cdot l) \cdot \cos(\alpha \cdot l + \Phi_l)}{N} \right) \cdot \left( \sum_{l=0}^{N-1} \frac{\sin(\alpha \cdot m) \cdot \cos(\alpha \cdot m + \Phi_m)}{N} \right) \right] \\ &= \sum_{l,m=0}^{N-1} \left( \frac{\cos(\alpha \cdot l) \cdot \sin(\alpha \cdot m)}{N^2} \cdot \mathbf{E} \left[ \cos(\alpha \cdot l + \Phi_l) \cdot \cos(\alpha \cdot m + \Phi_m) \right] \right) \end{aligned} \quad (38) $$
while

$$C_R = \sum_{l,m=0}^{N-1} \left( \frac{\cos(\alpha \cdot l) \cdot \sin(\alpha \cdot m)}{N^2} \cdot \mathbf{E}[\cos(\alpha \cdot l + \Phi_l)] \cdot \mathbf{E}[\cos(\alpha \cdot m + \Phi_m)] \right) \quad (39)$$

Since $(p_n)$ are independent, $(\Phi_n)$ are also independent and $\forall l \neq m$

$$\mathbf{E}[\cos(\alpha \cdot l + \Phi_l) \cdot \cos(\alpha \cdot m + \Phi_m)] = \mathbf{E}[\cos(\alpha \cdot l + \Phi_l)] \cdot \mathbf{E}[\cos(\alpha \cdot m + \Phi_m)] \quad (40)$$

The corresponding terms in $C_L$ and $C_R$ thus cancel out and

$$C = \frac{1}{N^2} \sum_{n=0}^{N-1} \cos(\alpha \cdot n) \cdot \sin(\alpha \cdot n) \cdot \underbrace{[\mathbf{E}[\cos^2(\alpha \cdot n + \Phi_n)] - \mathbf{E}[\cos(\alpha \cdot n + \Phi_n)]^2]}_{= \mathrm{Var}(\cos(\alpha \cdot n + \Phi_n))} \quad (41)$$

Using Popoviciu's inequality on variances to bound $\mathrm{Var}(\cos(...))$ yields

$$C \le \frac{1}{N^2} \sum_{n=0}^{N-1} \underbrace{|\cos(\alpha \cdot n)|}_{\le 1} \cdot \underbrace{|\sin(\alpha \cdot n)|}_{\le 1} \cdot \underbrace{\mathrm{Var}(\cos(...))}_{\le 1} \le \frac{1}{N} \xrightarrow[N \to +\infty]{} 0 \quad (42)$$

Thus $\Re(\tilde{D}_s) \perp \Im(\tilde{D}_s)$ asymptotically, and we finally demonstrated that $\tilde{D}_s$—and thus $\tilde{D}_N$—converges in law towards a complex normal distribution with increasing values of $N$. This can be rephrased as

$$\left\{ \begin{array}{l} \tilde{D}_N = x + i \cdot y \\ x \sim \mathcal{N}(\mu_x = \beta_p \cdot \cos(\varphi), \sigma^2) \\ y \sim \mathcal{N}(\mu_y = \beta_p \cdot \sin(\varphi), \sigma^2) \\ x \perp y \end{array} \right. \quad \text{wherein } \sigma^2 = \frac{1}{N} \cdot \left( 1 - \beta_p^2 + \frac{1}{\text{SNR}} \right) \quad (43)$$

### 3.3.3 Convergence in Practice

While this convergence is theoretically proven for $N \to +\infty$, its practical validity was studied for relatively small values of $N$. To do so, simulations were performed for SNR<sub>dB</sub>= 0..30 dB, $\sigma_p = 0.1..10^\circ$, $N = 10..100$. The multivariate normality of $\Re(\tilde{D}_N)$ and $\Im(\tilde{D}_N)$ was tested using Henze-Zirkler test [13] on $2.10^3 \tilde{D}_N$ outcomes. Their independence, was asserted using Hoeffding's D measures [11] between $\Re(\tilde{D}_N)$ and $\Im(\tilde{D}_N)$ using $10^5 \tilde{D}_N$ outcomes. In the case of the Henze-Zirkler test, simulations were repeated ten times, and their p-values-adjusted using the Benjamini-Hochberg correction [2] and combined with Fisher's method. The outcomes of these simulations were as follows. Henze-Zirkler multivariate test revealed no deviation from normality for $N$.
values above 20 using a significance level of 0.05. Hoeffding's D measures between $\Re(\tilde{D}_N)$ and $\Im(\tilde{D}_N)$, on their part, were found to be below $10^{-4}$ for as little as $N = 20$ measurement points, clearly demonstrating the practical independence of the two random variables. These two tests confirm that the convergence of $\tilde{D}_N$ towards a complex normal distribution is observed in practice across a wide range of conditions, even for relatively small $N$ values. Of note, this fast convergence is again demonstrated in Figure 3, further below, for as little as 20 samples.

# 4 Phase Estimation from the DFT

The phase $\varphi$ of the signal $S$ may be estimated by $\hat{\varphi} = \arg(\tilde{D}_N)$. Of paramount importance are thus the mean bias and RMSE of this estimator, i.e. the two quantities

$$ \mathbf{E}[\varphi - \hat{\varphi}] \quad \text{and} \quad \sqrt{\mathbf{E}[(\varphi - \hat{\varphi})^2]} \qquad (44) $$

In order to derive them, we can notice that $\hat{\varphi} = \arg(x + i \cdot y)$, and focus on the joint probability function of $(x, y)$, given by

$$ f(x, y) = \frac{1}{2 \cdot \pi \cdot \sigma^2} \cdot e^{-\frac{1}{2 \cdot \sigma^2} \cdot ((x - \mu_x)^2 + (y - \mu_y)^2)} \qquad (45) $$

Then—thanks to the LOTUS—$\hat{\varphi}$'s mean bias ($a = 1$) and RMSE ($a = 2$) may be computed using

$$ \mathbf{E}[(\varphi - \hat{\varphi})^a] = \iint_{\mathbb{R}^2} \arg(x + i \cdot y)^a \cdot f(x, y) \, dxdy \qquad (46) $$

While direct numerical calculations are presented in Section 4.3, another approach involving a switch to polar coordinates is first presented in the next section, allowing an in-depth comprehension of the influence of noises by means of meaningful illustrations.

## 4.1 Marginalisation in Polar Coordinates

A representation of the estimation of $\varphi$ in the complex plane can be seen in Figure 2. Naturally, different realisations of $S$ would lead to different $\tilde{D}_N$ values. Since $\tilde{D}_N$ follows a complex normal distribution, this translates into the small black dots (·) whose repartition is characteristic of a bivariate normal law, in the complex plane. This distribution is centred around its mean—$\mathbf{E}[\tilde{D}_N]$—of Cartesian coordinates $(\mu_x, \mu_y)_c$, represented as a large black dot (●). Interestingly, this centre is distinct from its position in the noiseless case, represented as a large blue dot (●), due to the $e^{-\sigma_p^2}$ term in
<figure><img src="image_3.png" /><figcaption>Figure 2: A representation of the estimation issue at hand in the complex plane. The (□)$_c$ and (□)$_p$ subscripts denote Cartesian and polar coordinates, respectively. See the text for further explanations.</figcaption></figure>

$\mathbf{E}[\tilde{D}_N]$, which is dragging the $\tilde{D}_N$ distribution towards the origin. A given $\tilde{D}_N$ realisation—depicted as a large red dot (●)—yields an estimation of $\varphi$, namely $\hat{\varphi}$.

Formally, the following change of variable from Cartesian to polar coordinates may be performed:

$$ \mathbf{F}: \begin{pmatrix} \mathbb{R}^+ \\ (r) \\ (\theta) \end{pmatrix} \mapsto \begin{pmatrix} \mathbb{R}^2 \\ (r \cdot \cos \theta) \\ (r \cdot \sin \theta) \end{pmatrix} \quad \text{and} \quad \mathbf{J}_{\mathbf{F}}(r, \theta) = r \quad (47) $$

wherein $\mathbf{J}_F$ is the Jacobian of $\mathbf{F}$. Expressing $f$ as a function of $r$ and $\theta$ yields

$$ \begin{align} f(r, \theta) &= \frac{1}{2 \cdot \pi \cdot \sigma^2} \cdot e^{-\frac{1}{2 \cdot \sigma^2} \cdot (r^2 - 2 \cdot r \cdot \beta_p \cdot \cos(\theta - \varphi) + \beta_p^2)} \nonumber \\ &= \frac{1}{2 \cdot \pi \cdot \sigma^2} \cdot e^{-\frac{(r - \beta_p \cdot \cos(\theta - \varphi))^2 + \beta_p^2 \cdot \sin^2(\theta - \varphi)}{2 \cdot \sigma^2}} \tag{48} \\ &= \frac{e^{-\beta_p^2 \cdot \frac{\sin^2(\theta - \varphi)}{2 \cdot \sigma^2}}}{2 \cdot \pi \cdot \sigma^2} \cdot e^{-\frac{1}{2 \cdot \sigma^2} \cdot (r - \beta_p \cdot \cos(\theta - \varphi))^2} \nonumber \end{align} $$

Applying the above-mentioned change of variables, taking care to replace $dxdy$ by $\mathbf{J}_F(r, \theta)drd\theta$ in the integrand leads to the following expression for the marginal density probability of $\theta$ under $\varphi$:

$$ \iint_{(x,y) \in \mathbb{R}^2} f(x,y)dxdy = \iint_{\theta \in [0; 2\pi[} f(r, \theta)r \, dr \, d\theta \quad (49) $$
and $f$ can be marginalised with respect to $r$ alone, leading to

$$ g_{\varphi}(\theta) \triangleq \int_{0}^{+\infty} f(\theta, r) \cdot r \, dr = \frac{e^{-\beta_p^2 \frac{\sin^2(\theta-\varphi)}{2 \cdot \sigma^2}}}{2 \cdot \pi \cdot \sigma^2} \cdot \underbrace{\int_{0}^{+\infty} r \cdot e^{-\frac{(r-\beta_p \cdot \cos(\theta-\varphi))^2}{2 \cdot \sigma^2}} dr}_{\mathcal{G}_{\theta}} \quad (50) $$

$\mathcal{G}_{\theta}$ can be further decomposed into $\mathcal{H}_{\theta}$ and $\mathcal{K}_{\theta}$, following

$$ \mathcal{G}_{\theta} = \underbrace{\int_{0}^{+\infty} (r - \beta_p \cdot \cos(\theta - \varphi)) \cdot e^{-\frac{(r - \beta_p \cdot \cos(\theta - \varphi))^2}{2 \cdot \sigma^2}} dr}_{\mathcal{H}_{\theta}} + \beta_p \cdot \cos(\theta - \varphi) \cdot \underbrace{\int_{0}^{+\infty} e^{-\frac{1}{2 \cdot \sigma^2} \cdot (r - \beta_p \cdot \cos(\theta - \varphi))^2} dr}_{\mathcal{K}_{\theta}} \quad (51) $$

And $\mathcal{H}_{\theta}$ and $\mathcal{K}_{\theta}$ can be explicitly calculated using appropriate change of variables as

$$ \mathcal{H}_{\theta} = \int_{-\beta_p \cdot \cos(\theta - \varphi)}^{+\infty} t \cdot e^{-\frac{t^2}{2 \cdot \sigma^2}} dt = -\sigma^2 \cdot \int_{-\beta_p \cdot \cos(\theta - \varphi)}^{+\infty} \frac{d}{dt} \left( e^{-\frac{t^2}{2 \cdot \sigma^2}} \right) dt = \sigma^2 \cdot e^{-\beta_p^2 \frac{\cos^2(\theta - \varphi)}{2 \cdot \sigma^2}} \quad (52) $$

and

$$ \begin{aligned} \mathcal{K}_{\theta} &= \int_{-\beta_p \cdot \cos(\theta - \varphi)}^{+\infty} e^{-\frac{t^2}{2 \cdot \sigma^2}} dt = \int_{-\beta_p \cdot \cos(\theta - \varphi)}^{+\infty} e^{-\left(\frac{t}{\sqrt{2} \cdot \sigma}\right)^2} dt = \sigma \cdot \sqrt{2} \cdot \int_{-\beta_p \cdot \frac{\cos(\theta - \varphi)}{\sigma \cdot \sqrt{2}}}^{+\infty} e^{-u^2} du \\ &= \sigma \cdot \sqrt{\frac{\pi}{2}} \cdot \frac{2}{\sqrt{\pi}} \cdot \underbrace{\int_{-\beta_p \cdot \frac{\cos(\theta - \varphi)}{\sigma \cdot \sqrt{2}}}^{+\infty} e^{-u^2} du}_{\text{erfc}(\dots)} = \sigma \cdot \sqrt{\frac{\pi}{2}} \cdot \text{erfc} \left( -\beta_p \cdot \frac{\cos(\theta - \varphi)}{\sigma \cdot \sqrt{2}} \right) \quad (53) \end{aligned} $$

and thus:

$$ \mathcal{G}_{\theta} = \sigma^2 \cdot e^{-\beta_p^2 \frac{\cos^2(\theta-\varphi)}{2 \cdot \sigma^2}} + \sigma \cdot \beta_p \cdot \sqrt{\frac{\pi}{2}} \cdot \cos(\theta - \varphi) \cdot \text{erfc} \left( -\beta_p \cdot \frac{\cos(\theta - \varphi)}{\sigma \cdot \sqrt{2}} \right) \quad (54) $$

Finally, putting it all together leads to

$$ g_{\varphi}(\theta) = \frac{e^{-\frac{\beta_p^2}{2 \cdot \sigma^2}}}{2 \cdot \pi} + \frac{\beta_p \cdot \cos(\theta - \varphi) \cdot e^{-\frac{\beta_p^2}{2 \cdot \sigma^2} \frac{\sin^2(\theta - \varphi)}{2 \cdot \sigma^2}}}{2 \cdot \sqrt{2 \cdot \pi} \cdot \sigma} \cdot \text{erfc} \left( -\beta_p \cdot \frac{\cos(\theta - \varphi)}{\sigma \cdot \sqrt{2}} \right) \quad (55) $$
$g_\varphi(\theta)$ is represented in Figure 3, along with the histogram of simulated $\hat{\varphi}$ values—i.e. $\arg(\tilde{D}_N)$ values, $\tilde{D}_N$ being computed from simulated $S$ vectors. $g_\varphi(\theta)$ resembles a normal distribution which would have been wrapped around the $[0, 2\pi]$ interval—although it is not a wrapped normal distribution nor a von Mises distribution strictly speaking (for further information on these distributions, see Collett et al. [6, 22, 21]). As could have been expected intuitively, the shape of $g_\varphi(\theta)$ narrows as $N$ or SNR$_{dB}$ increase, or as $\sigma_p$ decreases, as emphasised in the inset of Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Normalised representations of $g_\varphi(\theta)$ and of the histogram of simulated $\hat{\varphi}$ values ($10^6$ draws, $\varphi=60^\circ$). Main picture: $N = 20$, SNR$_{dB}=-10$ dB and $\sigma_p = 5^\circ$. The magnifying glass clearly shows how the histogram of simulated $\hat{\varphi}$ is close to the analytic expression of $g_\varphi(\theta)$. Inset: $N = 20$, SNR$_{dB}=0$ dB and $\sigma_p = 1^\circ$.</figcaption></figure>

## 4.2 Expressing the Error in the Polar Case

Since $g_\varphi(\theta)$ is the probability density function to measure a phase shift $\hat{\varphi} = \theta$ given a true phase shift $\varphi$, the mean bias and RMSE of the $\hat{\varphi}$ estimator, may be expressed as

$$ 
\begin{align} 
\mathbf{E}[\varphi - \hat{\varphi}] &= \int_{0}^{2\pi} \mathrm{err}_{\varphi}(\theta) \cdot g_{\varphi}(\theta) d\theta, \quad \text{and} \tag{56} \\ 
\mathrm{RMSE} &= \sqrt{\int_{0}^{2\pi} \mathrm{err}_{\hat{\varphi}}^{2}(\theta) \cdot g_{\varphi}(\theta) d\theta} \nonumber 
\end{align} 
$$
wherein err$_{ϕ}$ (θ) is the estimation error at a given polar angle θ, as represented in Figure 4, and is defined as

$$err_{\varphi}(\theta) = \min(d, 2\cdot\pi - d) \qquad (57)$$

with d ≡ θ − ϕ (mod 2 · π), that is

$$err_{\varphi}(\theta) = \pi - |\pi - \theta| \qquad (58)$$

<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of err on the trigonometric circle: err always represent the minimum angular distance between the true ϕ value and an angle θ.</figcaption></figure>

It can also be demonstrated that the RMSE and bias are independent of ϕ, as could have been expected intuitively. Only the RMSE case is detailed below, but a similar train of thought can be followed for the bias. Since the square root function is bijective on $R^+$, we only have to demonstrate that $\forall(\varphi, \varphi') \in [0, 2\cdot\pi]^2$

$$\underbrace{\int_{0}^{2\cdot\pi} \mathrm{err}_{\varphi}^{2}(\theta) \cdot g_{\varphi}(\theta) d\theta}_{\mathcal{J}_{\varphi}} = \underbrace{\int_{0}^{2\cdot\pi} \mathrm{err}_{\varphi'}^{2}(\theta) \cdot g_{\varphi'}(\theta) d\theta}_{\mathcal{J}_{\varphi'}} \qquad (59)$$

Then, with the change of variable $t = \theta - \varphi + \varphi'$,

$$\mathcal{J}_{\varphi} = \int_{-\varphi + \varphi'}^{2\cdot\pi - \varphi + \varphi'} \mathrm{err}_{\varphi}^{2}(t + \varphi - \varphi') \cdot g_{\varphi}(t + \varphi - \varphi') d\theta \qquad (60)$$

Using the definitions of err and g comes

$$\begin{aligned} \mathcal{J}_{\varphi} &= \underbrace{\int_{-\varphi+\varphi'}^{2\cdot\pi-\varphi+\varphi'} \mathrm{err}_{\varphi'}^{2}(t) \cdot g_{\varphi'}(t) dt}_{\mathcal{M}_{\varphi}} \\ &= \underbrace{\int_{-\varphi+\varphi'}^{0} \mathrm{err}_{\varphi'}^{2}(t) \cdot g_{\varphi'}(t) dt}_{\mathcal{M}_{\varphi}} + \int_{0}^{2\cdot\pi-\varphi+\varphi'} \mathrm{err}_{\varphi'}^{2}(t) \cdot g_{\varphi'}(t) dt \end{aligned} \qquad (61)$$
Using the change of variable $u = t + 2\cdot\pi$, and since err and g are $2\cdot\pi$-periodic,

$$ \mathcal{M}_{\varphi} = \int_{2\cdot\pi-\varphi+\varphi'}^{2\cdot\pi} \mathrm{err}_{\varphi'}^2(u-2\cdot\pi) \cdot g_{\varphi'}(u-2\cdot\pi) du = \int_{2\cdot\pi-\varphi+\varphi'}^{2\cdot\pi} \mathrm{err}_{\varphi'}^2(u) \cdot g_{\varphi'}(u) du \quad (62) $$

Leading to

$$ \begin{aligned} \mathcal{J}_{\varphi} &= \int_{0}^{2\cdot\pi-\varphi+\varphi'} \mathrm{err}_{\varphi'}^2(t) \cdot g_{\varphi'}(t) dt + \int_{2\cdot\pi-\varphi+\varphi'}^{2\cdot\pi} \mathrm{err}_{\varphi'}^2(u) \cdot g_{\varphi'}(u) du \\ &= \int_{0}^{2\cdot\pi} \mathrm{err}_{\varphi'}^2(t) \cdot g_{\varphi'}(t) dt = \mathcal{J}_{\varphi'} \end{aligned} \quad (63) $$

The RMSE is thus independent of $\varphi$, and taking $\varphi = 0$ leads to $err(\theta) = |\theta|$ $\forall \theta \in [-\pi, \pi]$, yielding finally

$$ \mathrm{RMSE} = \sqrt{\int_{-\pi}^{+\pi} \theta^2 \cdot g_0(\theta) d\theta} \quad (64) $$

However, due to the complexity of $g_0(\theta)$, we did not manage to derive a closed-form expression of the RMSE in the general case, and used numerical simulations to compute the latter as a function of $\sigma_x$—or, equivalently, the SNR—$\sigma_p$, and $N$.

The bias, on the other hand, can be readily computed since it can be shown similarly that:

$$ \mathbf{E}(\hat{\varphi} - \varphi) = \int_{-\pi}^{+\pi} \theta \cdot g_0(\theta) d\theta \quad (65) $$

and since $[\theta \mapsto \theta \cdot g_0(\theta)]$ is an odd function, it follows that the bias is null and that $\hat{\varphi}$ is an unbiased estimator of $\varphi$, which further justifies its choice as $\varphi$ estimator in the first place.

## 4.3 Numerical Calculation of the RMSE

Despite giving a better understanding of the issue at hand, and allowing one to clearly visualize the influence of $N$, $\sigma_p$ or $\sigma_x$ on the probability density function of $\hat{\varphi}$—as demonstrated by Figures 2 and 3—the polar coordinates considerations did not give a closed form expression for the RMSE of $\hat{\varphi}$. The latter may thus be calculated numerically using either Equation 46 (Cartesian case) or Equation 64 (polar case):

$$ \mathrm{RMSE} = \sqrt{\iint_{\mathbb{R}^2} \arg (x + i \cdot y)^2 \cdot f(x, y) \, dxdy} = \sqrt{\int_{-\pi}^{+\pi} \theta^2 \cdot g_0(\theta) d\theta} \quad (66) $$
While these two expression of the RMSE obviously yield the same quantity, the Cartesian approach is much more computationally expensive than the polar one, due to the double integration over the discretised, truncated $\mathbb{R}^2$ plane. Thus, all the results below were obtained using the polar approach. The influence of both the $\text{SNR}_{\text{dB}}$—and thus $\sigma_x$—and the phase noise $\sigma_p$ on the RMSE is depicted in Figure 5. Simulated RMSE values were also added to the latter figure, computed in a similar fashion to Figure 3. The closeness between these simulated values and the theoretical predictions consolidates the above-presented approach.

<figure><img src="image_6.png" /><figcaption>Figure 5: The RMSE on $\varphi$ estimation as a function of $\text{SNR}_{\text{dB}}$ for various $\sigma_p$ values. In the absence of phase noise—i.e. the black, continuous line (—)—a linear regression can be performed for $\text{SNR}_{\text{dB}} \ge -20$ dB, yielding the displayed equation and associated determination coefficient. As indicated, all calculations were performed taking $N=1000$ points. Simulated values were computed using only $10^2$ realisations so that some dispersion is perceivable. Indeed, using $10^6$ realisations as was done in Figure 3 led to marks indistinguishable from the plain red line.</figcaption></figure>

The evolution of the RMSE as a function of the $\text{SNR}_{\text{dB}}$ can be split into three behaviours:

- At very low $\text{SNR}_{\text{dB}}$ values—i.e. below -30 dB—the RMSE saturates, to reach approximately 100°. This corresponds to an exceedingly noisy case, for which $\hat{\varphi}$ is basically no more than a random guess on $[0; 2\pi[$. In this case, $g_{\varphi}(\theta)$ converges towards a uniform distribution, and the RMSE tends towards $\pi/\sqrt{3}$ rad ($\approx 104°$), as shown in Section 5.1.

- At higher $\text{SNR}_{\text{dB}}$ value and in the absence of phase noise—i.e. when $\text{SNR}_{\text{dB}} > -20$ dB and $\sigma_p=0$—the RMSE follows a linear relationship
with the SNR (or a log-linear relationship with $\text{SNR}_{\text{dB}}$, as indicates the equation on the graphic). In this case, $g_{\varphi}(\theta)$ converges towards a normal distribution centred around $\varphi$, and the RMSE tends towards $1/\sqrt{N \cdot \text{SNR}}$, as demonstrated in Section 5.2.

- Finally, when $\text{SNR}_{\text{dB}} \ge -20$ dB but the phase noise is significant—typically above $0.5^\circ$ in Figure 5—the RMSE becomes independent of $\text{SNR}_{\text{dB}}$ at some point as $\text{SNR}_{\text{dB}}$ increases, and becomes a function of $\sigma_p$ alone. This latter case is developed in Section 5.3.

The influence of $N$ on $\hat{\varphi}$ RMSE is presented in Figure 6. Starting with the left part of the figure, while the three above-mentioned behaviours can be observed irrespective of $N$ value, increasing $N$ by a factor $m$ has two main effects. First, it shifts the saturation threshold in case of extreme noise—i.e. for $\text{SNR}_{\text{dB}}$ values in the [-50,-20] dB range—by a factor $10 \cdot \log_{10}(m)$ (in dB) to the left. Then, it divides the RMSE at higher SNRs—i.e. for $\text{SNR}_{\text{dB}}$ values in the [-20,50] dB range—by a factor $\sqrt{m}$. Looking at the right part of Figure 6, $\hat{\varphi}$ RMSE appears to decrease linearly with increasing $N$ values, while increasing the SNR shifts the RMSE lines downwards.

<figure><img src="image_7.png" /><figcaption>Figure 6: The RMSE on $\varphi$ estimation as a function of $\text{SNR}_{\text{dB}}$ and $N$, with $\sigma_p=1^\circ$.</figcaption></figure>

Of note, the reader should bear in mind that the $\text{SNR}_{\text{dB}}$ thresholds given above, as well as the numerical values given below, are of course dependent on $N$ and $A_s$. As a general rule, $A_s$ is considered unitary in all our simulations, and $N=1000$ unless otherwise stated. Still, the same above-described behaviours would be observed with different $N$ and $A_s$ values, only the numerical values
stated in our developments would be altered. Of note, $\sigma_p$ values were also deliberately chosen relatively high for the sake of illustration. In practice, modern analogue-to-digital converters can feature phase noises in the $[10^{-2}-10^{-3}]^\circ$ range [16, 3].

Finally, Figure 7 depicts the evolution of $\hat{\varphi}$ RMSE with $\sigma_p$. Here, the behaviour observed in Figure 5 can be examined from another perspective: at low enough phase noise levels—i.e. when $\sigma_p$ is below approximately 0.01°—the RMSE is solely function of $N$ and the SNR. However, as soon as $\sigma_p$ increases, it starts to act as a noise floor, and increases the RMSE progressively. This effect is particularly pronounced at relatively high $N$ and SNR values, potentially ruining an otherwise excellent accuracy, as can be seen in the left part of the figure for $N = 10^5$ and SNR$_{dB} = 50$ dB.

<figure><img src="image_8.png" /><figcaption>Figure 7: The RMSE on $\varphi$ estimation as a function of $\sigma_p$, for various $N$ and SNR values.</figcaption></figure>

## 4.4 Estimation Efficiency

The efficiency of the performed estimation can be found by computing the Cramér-Rao Lower Bound (CRLB) of $\hat{\varphi}$. The likelihood $\mathcal{L}(\hat{\varphi}|\tilde{D}_N)$ of observing $\tilde{D}_N$ under a given $\hat{\varphi}$—hereafter noted $\mathcal{L}$ for the sake of conciseness—is given by

$$ \mathcal{L} = \prod_{z=\{x,y\}} \frac{1}{\sqrt{2 \cdot \pi \cdot \sigma^2}} \cdot e^{-\frac{(z-\mu_z)^2}{2 \cdot \sigma^2}} \quad (67) $$
The log-likelihood is then

$$ \log(\mathcal{L}) = -\log(2\pi\cdot\sigma^2) - \sum_{z=\{x,y\}} \frac{(z-\mu_z)^2}{2\cdot\sigma^2} \qquad (68) $$

and its first and second derivatives are given by

$$ \frac{\partial \log(\mathcal{L})}{\partial \hat{\varphi}} = \frac{1}{\sigma^2} \sum_{z=\{x,y\}} (z - \mu_z) \cdot \frac{\partial \mu_z}{\partial \hat{\varphi}} \qquad (69) $$

and

$$ \frac{\partial^2 \log(\mathcal{L})}{\partial \hat{\varphi}^2} = \frac{1}{\sigma^2} \sum_{z=\{x,y\}} \left[ (z - \mu_z) \cdot \frac{\partial^2 \mu_z}{\partial \hat{\varphi}^2} - \left( \frac{\partial \mu_z}{\partial \hat{\varphi}} \right)^2 \right] \qquad (70) $$

The Fisher information $\mathcal{I}(\hat{\varphi})$ may then be derived as

$$ \begin{aligned} \mathcal{I}(\hat{\varphi}) &= -\mathbf{E} \left[ \frac{\partial^2 \log(\mathcal{L})}{\partial \hat{\varphi}^2} \right] = \frac{1}{\sigma^2} \sum_{z=\{x,y\}} \left[ \left( \frac{\partial \mu_z}{\partial \hat{\varphi}} \right)^2 - \mathbf{E} \underbrace{[z - \mu_z] \cdot \frac{\partial^2 \mu_z}{\partial \hat{\varphi}^2}}_{=0} \right] \quad (71) \\ &= \frac{\beta_p^2}{\sigma^2} \end{aligned} $$

Finally leading to

$$ \text{CRLB} = \frac{1}{\mathcal{I}(\hat{\varphi})} = \frac{\sigma^2}{\beta_p^2} \qquad (72) $$

The convergence of $\hat{\varphi}$ towards its CRLB is illustrated in Figure 8. $\hat{\varphi}$ appears to be an asymptotically efficient estimator of $\varphi$ with a fast convergence rate, exhibiting 1-Efficiency values below $10^{-3}$ for as little as 1000 samples even in the presence of strong noise (SNR=0 dB). As a remainder, the efficiency of an unbiased estimator—which is the case for $\hat{\varphi}$—is defined as [28, p. 279]

$$ \text{Efficiency} = \frac{\text{CRLB}}{\text{Var}(\hat{\varphi})} = \frac{\text{CRLB}}{\text{RMSE}^2} \qquad (73) $$

# 5 Asymptotical Behaviours of the RMSE

## 5.1 Saturation in Case of Excessive Additive Noise

For SNR<sub>dB</sub> values below approximately -30 dB, the RMSE appears to be converging towards approximately 100°—see Figure 5. This phenomenon corresponds to an extremely noisy case, wherein the estimated $\varphi$ value is no
<figure><img src="image_9.png" /><figcaption>Figure 8: Asymptotical efficiency of $\hat{\varphi}$. $\sigma_p$ was set to 1°.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: $g_0(\theta)$ for different SNR<sub>dB</sub> values. For increasing SNR<sub>dB</sub> values, $g_0(\theta)$ converges towards a uniform distribution. $N=1000$, $\sigma_p=0$.</figcaption></figure>

better than a random guess on the $[0; 2\pi[$ interval. In this case, $g_0(\theta)$ converges toward a uniform distribution as SNR<sub>dB</sub> decreases, as illustrated in Figure 9.
Indeed, as SNR<sub>dB</sub> decreases, σ increases and we can make the following approximation

$$ g_0(\theta) = \frac{1}{2 \cdot \pi} \cdot \underbrace{\frac{e^{-\frac{\beta_p^2}{2\sigma^2}}}{\sigma \to +\infty}}_{1} + \underbrace{\frac{\beta_p \cdot \cos(\theta)}{2 \cdot \sqrt{2\cdot\pi} \cdot \sigma}}_{\substack{\xrightarrow{\sigma \to +\infty} 0 \\ \underbrace{|\dots| \le 1}_{\approx 0}}} \cdot \underbrace{\frac{e^{-\frac{\beta_p^2 \cdot \sin^2(\theta)}{2 \cdot \sigma^2}}}{<1}}_{\xrightarrow{\sigma \to +\infty} 1} \cdot \operatorname{erfc} \left( \underbrace{-\beta_p \cdot \frac{\cos(\theta)}{\sigma \cdot \sqrt{2}}}_{\xrightarrow{\sigma \to +\infty} 1} \right) \quad (74) $$

$$ g_0(\theta) \xrightarrow[\sigma \to +\infty]{\sigma \to +\infty} \frac{1}{2 \cdot \pi} $$

In other words, $g_0(\theta)$ can be approximated by a uniform distribution on the $[-\pi, \pi]$ interval. The Kullback-Leibler divergence between $g_0(\theta)$ and a uniform distribution can also be computed, and is presented in Figure 10. As expected, the divergence decreases steeply with a decreasing SNR<sub>dB</sub>, further confirming the above-mentioned convergence phenomenon.

<figure><img src="image_11.png" /><figcaption>Figure 10: The Kullback-Leibler divergence between $g_0(\theta)$ and a uniform distribution as a function of SNR<sub>dB</sub>. The inset shows the same data on a logarithmic scale. N=1000, σ<sub>p</sub>=0.</figcaption></figure>

Subsequently the RMSE becomes

$$ \text{RMSE} = \sqrt{\int_{-\pi}^{\pi} \theta^2 \cdot \frac{1}{2 \cdot \pi} d\theta} = \frac{\pi}{\sqrt{3}} \quad (75) $$
hence the saturation behaviour observed for the RMSE at high $\sigma_c$ values in Figure 5. The 100° plateau value noted above simply comes from the radian to degree conversion ($\pi/\sqrt{3}$ rad≈104°). That being said, contrary to the linear case given in the next section, this saturation phenomenon is of little practical interest since it corresponds to an extremely noisy case, which only yields random guesses as phase estimation. It was thus only presented here for the sake of completeness.

## 5.2 Linear Relationship With the SNR

The observed linear relationship between the RMSE and the SNR can be explained by the fact that, for increasing SNR values, $g_0(\theta)$ converges towards a normal distribution, as illustrated in Figure 11.

<figure><img src="image_12.png" /><figcaption>Figure 11: Normalised probability density functions: $g_0(\theta)$ is plotted with its associated simplified Gaussian model. For $\text{SNR}_{\text{dB}}$ above -10 dB, $g_0(\theta)$ becomes indistinguishable from its Gaussian approximation. $N=1000$, $\sigma_p=0$.</figcaption></figure>

Indeed, taking $\varphi = 0$ and high $\text{SNR}_{\text{dB}}$ values leads to an extremely narrow $g_0(\theta)$ function—as was already noted above in Figure 3—which is non-negligible only for very small $\theta$ deviations from zero. More formally:

$$ \begin{align} & \mathrm{Var}(\tilde{D}_N) \xrightarrow[\text{SNR}_{\text{dB}} \to +\infty]{\phantom{=}} 0 \quad (*) \nonumber \\ (*) & \implies \tilde{D}_N \xrightarrow[\text{SNR}_{\text{dB}} \to +\infty]{L_2} \mathbf{E}[\tilde{D}_N] \tag{76} \\ (*) & \implies \theta = \arg(\tilde{D}_N) - \varphi \xrightarrow[\text{SNR}_{\text{dB}} \to +\infty]{P} 0 \nonumber \end{align} $$
wherein $\underline{L_2}$ and $\underline{P}$ denote convergence in mean square and convergence in probability, respectively. By definition of the latter convergence, $g_0(\theta)$ is thus non-negligible only for very small deviations from zero with a non-null probability. Under these conditions, and in the absence of phase noise, $\beta_p = 1$ and $\sigma$ tends towards zero as $\text{SNR}_{\text{dB}}$ tends towards infinity. We can then make the following approximation:

$$ g_0(\theta) = \frac{e^{-\frac{\beta_p^2}{2\sigma^2}}}{2\pi} + \frac{\beta_p}{2\sqrt{2\cdot\pi}\sigma} \cdot \underbrace{\cos(\theta)}_{\approx 1} \cdot e^{-\beta_p^2 \frac{\sin^2(\theta)}{2\cdot\sigma^2}} \cdot \underbrace{\text{erfc}\left(-\beta_p \frac{\cos(\theta)}{\sigma\sqrt{2}}\right)}_{\approx 2} \quad (77) \\ \approx \frac{1}{\sqrt{2\cdot\pi}\sigma} \cdot e^{-\frac{\theta^2}{2\cdot\sigma^2}} $$

Thus, for high enough SNR values, $g_0(\theta)$ can be reasonably well approximated by a simple normal distribution of null mean, and variance $\sigma^2$. The fast convergence towards this approximation is further illustrated in Figure 12, wherein the Kullback-Leibler divergence and Bhattacharyya distance between $g_0(\theta)$ and its Gaussian approximation are represented as a function of $\text{SNR}_{\text{dB}}$. It appears that for $\text{SNR}_{\text{dB}}$ values above -10 dB, this distance becomes virtually null ($\approx 10^{-5}$).

<figure><img src="image_13.png" /><figcaption>Figure 12: The Kullback-Leibler divergence (KL) and Bhattacharyya distance (Bhat.) between $g_0(\theta)$ and its Gaussian approximation as a function of $\text{SNR}_{\text{dB}}$. The inset shows the same data on a logarithmic scale. $N=1000$, $\sigma_p=0$.</figcaption></figure>
As a reminder, the Bhattacharya distance quantifies the closeness between two distributions. Given two random variables P and Q with probability density function p(x) and q(x) the Bhattacharya distance and Kullback-Leibler divergence are defined as follows [15]:

(i) Bhattacharya distance:

$$d_B(P \| Q) = - \log \left( \int p(x) \cdot q(x) \, dx \right) \qquad (78)$$

(ii) Kullback-Leibler divergence:

$$d_{KL}(P \| Q) = \int p(x) \cdot \log \left( \frac{p(x)}{q(x)} \right) dx \qquad (79)$$

One may notice that the division by q(x) in $d_{KL}$ may be problematic in case of numerical applications if q(x) is near zero. Such an issue can be observed in Figure 12 for SNR<sub>dB</sub> values above roughly -8 dB: $d_{KL}$ cannot be computed even using 64-bits double-precision floats, because the Gaussian probability density function tends towards zero extremely fast as soon as $\theta$ deviates from zero. This is the reason why the Bhattacharya distance was introduced in the first place, so as to better cover the case of high SNR<sub>dB</sub> values.

Subsequently the RMSE becomes

$$\text{RMSE} = \sqrt{\int_{-\pi}^{\pi} \theta^2 \cdot \frac{1}{\sqrt{2 \cdot \pi \cdot \sigma}} \cdot e^{-\frac{\theta^2}{2 \cdot \sigma^2}} d\theta} = \sigma = \frac{1}{\sqrt{N \cdot \text{SNR}}} \qquad (80)$$

hence the linear relationship observed in Figure 6 between $\log_{10}(\text{RMSE})$ on the one hand, and SNR<sub>dB</sub> and N on the other hand. This result is especially interesting for practical applications. Indeed, a SNR<sub>dB</sub> above -10 dB can easily be reached in practice, ensuring the validity of the above-mentioned Gaussian approximation for $g_0(\theta)$, and the ensuing conclusions on the RMSE. Of particular interest, if this condition is fulfilled, **the expected RMSE on the phase estimation can be directly inferred from N and the SNR**, using Equation 80. Still, one should bear in mind that this conclusions holds only if $\beta_p \approx 1$, otherwise the developments presented in the next section should apply.
## 5.3 Saturation in the Presence of Phase Noise

When some amount of phase noise is present—i.e. $\beta_p \in ]0; 1[$—and for high enough SNR, the following approximation can be made, as was done in Equation 77:

$$g_0(\theta) \approx \frac{\beta_p}{\sqrt{2 \cdot \pi} \cdot \sigma} \cdot e^{-\frac{\beta_p^2 \cdot \theta^2}{2 \cdot \sigma^2}} \quad (81)$$

Again, $g_0(\theta)$ can be approximated by a simple normal distribution, but of variance $\sigma^2/\beta_p^2$, instead of $\sigma^2$ alone in the previous section. The convergence towards a Gaussian is nearly identical with that presented in Section 5.2 and thus Figures 11 and 12 were not reproduced for the sake of conciseness. Similarly to Equation 80 comes

$$\text{RMSE} = \frac{\sigma}{\beta_p} = \sqrt{\frac{1}{\beta_p^2 \cdot N} \cdot \left(1 - \beta_p^2 + \frac{1}{\text{SNR}}\right)} \quad (82)$$

which, in case of high enough SNR—i.e. $1 - \beta_p^2 \gg 1/\text{SNR}$—becomes

$$\text{RMSE} \approx \sqrt{\frac{1}{N} \cdot \left(\frac{1}{\beta_p^2} - 1\right)} \quad (83)$$

This result is also interesting because it provides a lower limit for the RMSE, even at exceedingly large SNR values: **the RMSE is ultimately limited by the phase noise, which acts as a noise floor.** This explains the saturations observed on the right part of Figure 5 and 6: the lower limits reached by the different curves with non-zero phase noises directly depend on their respective $\sigma_p$ values, following Equation 83. Most interestingly, Equation 82 gives a generic expression for the RMSE at reasonably high SNR values taking into account the joint influences of: (i) the number of points $N$, (ii) the amplitude of the phase noise—through $\beta_p$—and (iii) that of the additive noise—through the SNR.

# 6 Conclusion

This article presents a thorough analysis of the influence of additive and phase noises on the accuracy of the phase measurement of a known-frequency sinusoidal signal. More specifically, we focused on synchronous detection, a measurement scheme for which the number of collected samples $N$ on the one hand, and the sampling and probing frequency $f_s$ and $f_0$ on the other hand can be chosen so that $f_0 \cdot N/f_s$ is an integer. In this particular case, a closed-form expression of the PDF of the phase estimate $\hat{\varphi}$ could be derived, depending on $N$ and on the levels of additive and phase noises—$\sigma_x$ and $\sigma_p$,
respectively. $\hat{\phi}$ was also shown to be asymptotically efficient, with a fast convergence towards its CRLB, even using a limited number of samples in the presence of substantial noise levels.

When using the above-mentioned PDF to compute $\hat{\phi}$ RMSE, three main behaviours could be identified: (i) in case of excessive noise, the RMSE saturates towards the random guess situation, (ii) as the SNR increases, the RMSE decreases linearly with the square roots of the SNR and N, and (iii) as the SNR further decreases, the RMSE saturates again, reaching a noise floor caused by phase noise. While (i) is of little practical interest, (ii) and (iii) are of major importance in practical scenarios, since they can tell the experimenter whether lowering the RMSE should be achieved by taking more samples or by increasing the SNR—e.g. by increasing the emission power. Indeed, in case (iii), increasing the SNR is no use once the phase noise floor is reached, and only taking more samples can yield lower RMSE values.

In addition to its theoretical significance, this paper is thus also of practical value, allowing for informed decision-making when designing a phase-measuring apparatus. In particular, in the context of f-DLR mentioned in introduction, a compromise often has to be made between reducing the number of samples and the illumination power—which preserves the involved dyes from photobleaching, and saves power in case of battery-powered devices—and reducing the RMSE by increasing the two latter parameters—at the expense of power consumption, computing costs, and dye photobleaching. Of particular importance, we recently published a short communication demonstrating a real-world implementation of f-DLR applied to Carbon dioxide ($CO_2$) sensing [9]. In the latter work, we validated the theoretical results derived in the present paper in the case of a low phase noise and a relatively high SNR—i.e. Equation 80.

Finally, one possible limitation of this work resides in its simplified modelling of phase noise. Indeed, phase noise—studied through its influence on the power spectral density of a clock signal in the frequency domain [19], or by the distribution of the jitter it induces in the time domain [26]—can result from several phenomena that can be either purely random or partially deterministic [8, Chap. 2–4]. While a thorough review on the origins and statistical properties of phase noise is clearly out of the scope of this paper, and although some oscillators do exhibit Gaussian phase jitters in practice [20, 25, 5], this may not be the case in general [26]. As a result, future research could focus on extending our work to address non-Gaussian and/or non-i.i.d. phase noise.

## Acknowledgments

We are grateful to Yoshitate Takakura for his valuable insights in spectral estimation, and to Morgan Madec for his meticulous early review of this work.
References

[1] D. Atamanchuk, A. Tengberg, P. J. Thomas, J. Hovdenes, A. Apostolidis, C. Huber, and P. O. Hall, "Performance of a Lifetime-Based Optode for Measuring Partial Pressure of Carbon Dioxide in Natural Waters", *Limnology and Oceanography: Methods*, 12(2), 2014, 63–73, DOI: 10.4319/lom.2014.12.63.

[2] Y. Benjamini and Y. Hochberg, "Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing", *Journal of the Royal Statistical Society. Series B (Methodological)*, 57(1), 1995, 289–300, ISSN: 00359246, http://www.jstor.org/stable/2346101.

[3] C. E. Calosso, A. C. C. Olaya, and E. Rubiola, "Phase-Noise and Amplitude-Noise Measurement of DACs and DDSs", in *2019 Joint Conference of the IEEE International Frequency Control Symposium and European Frequency and Time Forum (EFTF/IFC)*, 2019 Joint Conference of the IEEE International Frequency Control Symposium and European Frequency and Time Forum (EFTF/IFC), 2019, 431–9, DOI: 10.1109/FCS.2019.8856100.

[4] C. Cerovecki and S. Hörmann, "On the CLT for Discrete Fourier Transforms of Functional Time Series", *Journal of Multivariate Analysis*, 154 (February), February 2017, 282–95, ISSN: 0047-259X, https://www.sciencedirect.com/science/article/pii/S0047259X16301592.

[5] G. Ciarpi, D. Monda, M. Mestice, D. Rossi, and S. Saponara, "Asymmetric 5.5 GHz Three-Stage Voltage-Controlled Ring-Oscillator in 65 nm CMOS Technology", *Electronics*, 12(3), 2023, ISSN: 2079-9292, DOI: 10.3390/electronics12030778.

[6] D. Collett and T. Lewis, "Discriminating Between the Von Mises and Wrapped Normal Distributions", *Australian Journal of Statistics*, 23(1), March 1981, 73–9, ISSN: 0004-9581, DOI: 10.1111/j.1467-842X.1981.tb00763.x, https://doi.org/10.1111/j.1467-842X.1981.tb00763.x.

[7] J. Cooley, P. Lewis, and P. Welch, "The Finite Fourier Transform", *IEEE Transactions on Audio and Electroacoustics*, 17(2), 1969, 77–85, DOI: 10.1109/TAU.1969.1162036.

[8] N. Da Dalt and A. Sheikholeslami, *Understanding Jitter and Phase Noise: A Circuits and Systems Perspective*, Cambridge: Cambridge University Press, 2018, ISBN: 9781107188570, DOI: 10.1017/9781316981238.

[9] E. Dervieux and W. Uhring, "Dual Lifetime Referencing for Accurate CO2 Sensing: An Experimental in vitro Validation", in *2024 22nd IEEE Interregional NEWCAS Conference (NEWCAS)*, 2024 22nd IEEE Interregional NEWCAS Conference (NEWCAS), 2024, 253–7, DOI: 10.1109/NewCAS58973.2024.10666315.
[10] A. Ferrero and R. Ottoboni, "High-Accuracy Fourier Analysis Based on Synchronous Sampling Techniques", *IEEE Transactions on Instrumentation and Measurement*, 41(6), 1992, 780-5, ISSN: 1557-9662, DOI: 10.1109/19.199406, https://doi.org/10.1109/19.199406.

[11] A. Fujita, J. R. Sato, M. A. A. Demasi, M. C. Sogayar, C. E. Ferreira, and S. Miyano, "Comparing Pearson, Spearman and Hoeffding's D Measure for Gene Expression Association Analysis", *Journal of Bioinformatics and Computational Biology*, 07(04), August 2009, 663-84, ISSN: 0219-7200, DOI: 10.1142/S0219720009004230, https://doi.org/10.1142/S0219720009004230.

[12] A. A. Girgis and F. M. Ham, "A Quantitative Study of Pitfalls in the FFT", *IEEE Transactions on Aerospace and Electronic Systems*, AES-16(4), 1980, 434-9, ISSN: 1557-9603, DOI: 10.1109/TAES.1980.308971, https://doi.org/10.1109/TAES.1980.308971.

[13] N. Henze and B. Zirkler, "A Class of Invariant Consistent Tests for Multivariate Normality", *Communications in Statistics - Theory and Methods*, 19(10), January 1990, 3595-617, ISSN: 0361-0926, DOI: 10.1080/03610929008830400, https://doi.org/10.1080/03610929008830400.

[14] X. Huang, Z. Wang, L. Ren, Y. Zeng, and X. Ruan, "A Novel High-Accuracy Digitalized Measuring Phase Method", in *2008 9th International Conference on Signal Processing, 2008 9th International Conference on Signal Processing*, 2008, 120-3, DOI: 10.1109/ICOSP.2008.4697084.

[15] T. Kailath, "The Divergence and Bhattacharyya Distance Measures in Signal Selection", *IEEE Transactions on Communication Technology*, 15(1), 1967, 52-60, DOI: 10.1109/TCOM.1967.1089532.

[16] W. Kester, "Converting Oscillator Phase Noise to Time Jitter", tech. rep. No. MT-008, Analog Devices, 2010.

[17] I. Klimant, C. Huber, G. Liebsch, G. Neurauter, A. Stangelmayer, and O. S. Wolfbeis, "Dual Lifetime Referencing (DLR) — A New Scheme for Converting Fluorescence Intensity into a Frequency-Domain or Time-Domain Information", in, *New Trends in Fluorescence Spectroscopy: Applications to Chemical and Life Sciences*, Berlin, Heidelberg: Springer Berlin Heidelberg, 2001, 257-74, ISBN: 978-3-642-56853-4, DOI: 10.1007/978-3-642-56853-4\_13.

[18] A. Lapidoth, *A Foundation in Digital Communication*, 2nd ed., Cambridge: Cambridge University Press, 2017, ISBN: 9781107177321, DOI: 10.1017/9781316822708.

[19] D. B. Leeson, "Oscillator Phase Noise: A 50-Year Review", *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, 63(8), 2016, 1208-25, ISSN: 1525-8955, DOI: 10.1109/TUFFC.2016.2562663.

[20] P. Lesage and C. Audoin, "Characterization of Frequency Stability: Uncertainty due to the Finite Number of Measurements", *IEEE Trans-
actions on Instrumentation and Measurement, 22(2), 1973, 157–61, ISSN: 1557-9662, DOI: 10.1109/TIM.1973.4314128.

[21] C. Ley and T. Verdebout, *Modern Directional Statistics*, Chapman and Hall / CRC, 2017, ISBN: 9781315119472, DOI: 10.1201/9781315119472.

[22] K. V. Mardia and P. E. Jupp, *Directional Statistics*, Wiley Online Library, 1999, ISBN: 978-0-471-95333-3.

[23] S. L. Marple, “A Tutorial Overview of Modern Spectral Estimation”, in *International Conference on Acoustics, Speech, and Signal Processing, International Conference on Acoustics, Speech, and Signal Processing*, 1989, 2152–2157 vol.4, DOI: 10.1109/ICASSP.1989.266889.

[24] S. Marple, *Digital Spectral Analysis: Second Edition*, Dover Books on Electrical Engineering, Dover Publications, 2019, ISBN: 9780486780528.

[25] C. Nazareth, P. Kaufmann, A. Alves Ferreira Júnior, M. Paiva, and J. W. Swart, “New Jitter Measurement Technique Using TDC Principle in a FPGA Component”, *Seminário Internacional de Metrologia Elétrica – VIII SEMETRO*, 2009.

[26] C.-K. Ong, D. Hong, K.-T. T. Cheng, and L.-C. Wang, “Random Jitter Extraction Technique in a Multi-Gigahertz Signal”, in *Proceedings Design, Automation and Test in Europe Conference and Exhibition, Vol. 1, Proceedings Design, Automation and Test in Europe Conference and Exhibition*, 2004, 286–91, DOI: 10.1109/DATE.2004.1268862.

[27] V. M. Panaretos and S. Tavakoli, “Fourier Analysis of Stationary Time Series in Function Space”, *The Annals of Statistics*, 41(2), April 2013, 568–603, DOI: 10.1214/13-AOS1086, https://doi.org/10.1214/13-AOS1086.

[28] D. Z. Robert V. Hogg Elliot Tanis, *Probability and Statistical Inference*, 9th, Pearson, 2014, ISBN: 9781292062358.

[29] L. Scharf and C. Demeure, *Statistical Signal Processing: Detection, Estimation, and Time Series Analysis, Addison-Wesley series in electrical and computer engineering*, Addison-Wesley Publishing Company, 1991, ISBN: 9780201190380.

[30] S. Schuster, S. Scheiblhofer, and A. Stelzer, “The Influence of Windowing on Bias and Variance of DFT-Based Frequency and Phase Estimation”, *IEEE Transactions on Instrumentation and Measurement*, 58(6), 2009, ISSN: 1557-9662, DOI: 10.1109/TIM.2008.2006131, https://doi.org/10.1109/TIM.2008.2006131.

[31] J. Soch, T. J. Faulkenberry, K. Petrykowski, and C. Allefeld, *The Book of Statistical Proofs*, version 2020, December 2020, DOI: 10.5281/zenodo.4305950.

[32] C. Staudinger, M. Strobl, J. P. Fischer, R. Thar, T. Mayr, D. Aigner, B. J. Müller, B. Müller, P. Lehner, G. Mistlberger, E. Fritzsche, J. Ehgartner, P. W. Zach, J. S. Clarke, F. Geißler, A. Mutzberg, J. D. Müller, E. P. Achterberg, S. M. Borisov, and I. Klimant, “A Versatile Optode System for Oxygen, Carbon Dioxide, and pH Measurements in Seawater With
Integrated Battery and Logger", *Limnology and Oceanography: Methods*, 16(7), 2018, 459–73, DOI: 10.1002/lom3.10260.

[33] P. Stoica and R. Moses, *Spectral Analysis of Signals*, Pearson Prentice Hall, 2005, ISBN: 9780131139565.

[34] T. Su, M. Yang, T. Jin, and R. C. C. Flesch, "Power Harmonic and Interharmonic Detection Method in Renewable Power Based on Nuttall Double-Windows All-Phase FFT Algorithm", *IET Renewable Power Generation*, 12(8), June 2018, 953–61, ISSN: 1752-1416, DOI: 10.1049/iet-rpg.2017.0115, https://doi.org/10.1049/iet-rpg.2017.0115.

[35] T. B. Tufan and U. Guler, "A Miniaturized Transcutaneous Carbon Dioxide Monitor Based on Dual Lifetime Referencing", in *2022 IEEE Biomedical Circuits and Systems Conference (BioCAS)*, 2022, 144–8, DOI: 10.1109/BioCAS54905.2022.9948600.

[36] A. van der Vaart and J. Wellner, *Weak Convergence and Empirical Processes: With Applications to Statistics*, Springer Series in Statistics, Springer, 1996, ISBN: 9780387946405.

[37] C. von Bültzingslöwen, A. McEvoy, C. Mcdonagh, B. Maccraith, I. Klimant, C. Krause, and O. Wolfbeis, "Sol-gel Based Optical Carbon Dioxide Sensor Employing Dual Luminophore Referencing for Application in Food Packaging Technology", *The Analyst*, 127 (December), December 2002, 1478–83, DOI: 10.1039/B207438A.

[38] H. Xiaohong, W. Zhaohua, and C. Guoqiang, "New Method of Estimation of Phase, Amplitude, and Frequency Based on all Phase FFT Spectrum Analysis", in *2007 International Symposium on Intelligent Signal Processing and Communication Systems, 2007 International Symposium on Intelligent Signal Processing and Communication Systems*, 2007, 284–7, DOI: 10.1109/ISPACS.2007.4445879.