Noise Parameter Identification in Ensemble of Atomic Clocks

Introduction

Precise timekeeping acts as the invisible backbone of modern infrastructure, synchronizing everything from global satellite navigation to high-frequency trading and distributed telecommunications [1], [2]. To guarantee extreme stability and robustness against individual hardware failures, high-quality time scales rely on ensembles of atomic clocks. While classical aggregation methods successfully reduce overall frequency instability, the next frontier in precision timekeeping demands a shift toward dynamic, model-based estimation approaches.

Parametric models such as state-space models can significantly enhance the stability and integrity of a time scale by incorporating prior knowledge of clock behavior. However, the efficacy of these model-based algorithms is fundamentally bounded by the accuracy of their underlying stochastic models. Real-world atomic clocks exhibit complex, often correlated noise profiles that are notoriously difficult to isolate and quantify. Therefore, the robust identification of noise parameters within an ensemble is a critical prerequisite. Overcoming this parameter identification challenge enables optimal weighting strategies and the timely detection of anomalies, bridging the gap between theoretical clock models and maximum real-world performance.

Although stochastic clock modeling and noise identification have been studied for decades, reliable estimation remains difficult in practical scenarios characterized by limited data records, correlated observations, and incomplete observability of individual clocks. Moreover, several theoretical aspects of the problem remain open, particularly regarding identifiability, sub-space unobservability, estimator consistency, or the statistical characterization of estimation uncertainty under realistic operating conditions. Furthermore, existing approaches predominantly address isolated aspects of the problem, seldom offering a unified framework that is simultaneously comprehensive, computationally efficient, and practically applicable.

Related Works

An ensemble of atomic clock is an underlying component of many applications that range from timekeeping to positioning. A successful operation of an ensemble relies on several algorithms that provide anomaly monitoring, clock combination, or clock steering. Pioneering clock combination algorithms are adaptive and only little information about clocks in the ensemble is required to initialize these algorithms [3]. On the other hand, there are algorithms which assume that a complete model of the clock ensemble is available [2], [4]. Therefore, there has been a concentrated effort devoted to clock ensemble identification. Since the the underlying dynamic is given by law of physics, the main effort is concentrated on estimating the parameters of the noises.

The parameter estimation method is tightly related to a selected model. Although there are some works focusing on ARIMA models [5], [6], most of the works related to the clock ensemble identification are based on the state space models [7] that are physically interpretable. A particular attention was devoted to the estimation of the frequency drift [8] from time deviations, the first order difference of time deviations (i.e. average normalized frequency deviations), and the second order difference of time deviations.

One approach to noise parameter identification is based on the maximum likelihood estimation method [9, 10, 11], where estimates of unknown model parameters are found to maximize the likelihood of measurements [12]. The likelihood of measurement is written using a chain rule as the product of measurement predictive probability density functions that are obtained using the Kalman filter. Missing measurements can be handled rather easily and apparent measurement outliers can be removed while constructing the likelihood. A potential drawback of the maximum likelihood estimation is the necessity of relying on an iterative optimization method as the likelihood is a nonlinear function of model parameters. To avoid negative variance estimates, either explicit constrained optimization is performed or variances are replaced by positive functions of auxiliary variables.

Another approach to noise parameter identification is based on the three-cornered-hat method that allows variances of three independent random variables to be estimated even if only samples of their pair-wise differences are available. This idea was applied to Allan variance (AVAR) plots under the assumption of noise-free measurements and independent internal noises of clocks in [13]. If more than three clocks are available, triads are formed and parameter estimated. Later on, the three-cornered-hat method has been extended to m-cornered-hat method [14] where a comparison with maximum likelihood method is given for a rather simplified ensemble model. The m-cornered-hat method also suffers from possibility of getting negative variance estimates. Several reasons for this issue are discussed in literature. The variance estimate uncertainty due to insufficient number of measurements is suspected in [13]. In [15], the measurement noise that influences the AVAR plot for short averaging intervals is suspected as a possible reason and the Groslambert covariance (Allan covariance (ACOV)) [16] is used instead of the AVAR under the assumption that the measurement noises are uncorrelated to make covariance estimates independent of measurement noise. A more in-depth discussion of possible reasons was given in [17]. It was argued that most commonly considered reasons in literature are not plausible and conclude that the clock noises should be admitted to be correlated. However, it makes the variance estimation problem ambiguous. These issue is discussed for the three-cornered-hat method in [17] and this discussion is extended to the m-cornered-hat method in [18] and [19].

Contribution

This paper addresses the problem of identifying noise parameters in a clock ensemble following the standard approach of utilizing differential measurements between a selected pivot clock and the remaining clocks. The main objective is to develop an effective identification method tailored to this measurement configuration, enabling reliable estimation of ensemble noise characteristics under realistic observational constraints.

Two noise parameter estimation methods are presented in this paper. The first one represents the extension of an underlying correlation approach for noise identification presented in [20]. Note that this approach was considered in the context of the clock ensemble identification in [21] but with unrealistic assumptions regarding the direct measurement of time deviations. The second presented methods builds upon the Hadamard covariance (HCOV) and m-cornered-hat method.

This is similar to the concept in [15], [22], where, however, the authors use Allan covariance and highlight its advantages over the conventional three-cornered-hat method. Both designed identification methods represent an original approach; they are broadly applicable to a wide class of ensemble atomic clocks. Moreover, their design is transparent and builds on the strong theoretical foundation of linear estimation.

A Clock and Ensemble Model Definition

The paper follows the standard modeling on atomic clock [7]. It is assumed that the $i$-th clock of an ensemble can be described at a time step $k = 0, 1, 2 \ldots, N$ by a second-order discrete-time state-space model²

where $\mathbf{x}_k^{(i)} = [x_{1,k}^{(i)}, x_{2,k}^{(i)}]^{\mathrm{T}} \in\mathbb{R}^{2}$ is a state consisting of a time deviation $x_{1,k}^{(i)} \in\mathbb{R}$ and a random walk process $x_{2,k}^{(i)} \in\mathbb{R}$, $\boldsymbol{\mu}^{(i)} \in\mathbb{R}^{2}$ is a drift term, $\mathbf{w}_k^{(i)} \in\mathbb{R}^{2}$ is a zero mean white Gaussian state noise with covariance matrix (CM) $\mathbf{Q}^{(i)} \in\mathbb{R}^{2 \times2}$, and $\Delta_k^{(i)} \in\mathbb{R}$ is a time deviation. The matrices $\mathbf{F}_1$, $\mathbf{Q}^{(i)}$ and vectors $\mathbf{h}_1$, $\boldsymbol{\mu}^{(i)}$ are given as

where $T_s \in\mathbb{R}^{++}$ is a sampling period in seconds, $d^{(i)} \in\mathbb{R}$ is a frequency drift, and $q_1^{(i)} \in\mathbb{R}^{++}$ and $q_2^{(i)} \in\mathbb{R}^{++}$ denote intensities of original continuous-time stochastic processes that represent the white noise frequency modulation and random walk frequency modulation, respectively.

The model of ensemble of $n$ clocks can be defined as

where $\mathbf{x}_k \in\mathbb{R}^{2n}$ is an ensemble state, $\boldsymbol{\mu} \in\mathbb{R}^{2n}$ is ensemble drift, $\mathbf{w}_k \in\mathbb{R}^{2n}$ is an ensemble state noise with zero mean and CM $\mathbf{Q} \in\mathbb{R}^{2n \times2n}$, $\mathbf{z}_k \in\mathbb{R}^{n_z}$ is a measurement containing $n_z = n - 1$ differences of time deviations of a pivot clock and all other clocks, and $\mathbf{v}_k \in\mathbb{R}^{n_z}$ is a white Gaussian measurement noise with zero mean and CM $\mathbf{R} \in\mathbb{R}^{n_z \times n_z}$. The vectors and matrices are given as

where $\mathbf{C}=[-\mathbf{1}_{n_z},\mathbf{I}_{n_z}]\in\mathbb{R}^{n_z\times n}$ is a matrix that specifies the first clock in the ensemble to be the pivot clock.

Problem Formulation

In principle, the problem does not permit an unambiguous estimation of the frequency drift of all clocks. Therefore, a common approach is to assume zero drift for the pivot clock and estimate the remaining drifts relative to it, i.e., the estimate the drift differences between the pivot and the other clocks in the ensemble. Thus, let’s define a vector of unknown model parameters as

where $\boldsymbol{\theta}\in\mathbb{R}^{n_\theta}$ with $n_\theta=(n^2+5n-2)/2$ denoting the number of the unknown parameters, $\mathbf{q}_1^{\mathsf{T}}=[q_1^{(1)},q_1^{(2)},\ldots,q_1^{(n)}]\in\mathbb{R}^n$, $\mathbf{q}_2^{\mathsf{T}}=[q_2^{(1)},q_2^{(2)},\ldots,q_2^{(n)}]\in\mathbb{R}^n$, $\mathbf{r}^{\mathsf{T}}=\operatorname{vecu}(\mathbf{R})\in\mathbb{R}^{n_r}$ collects elements of the upper triangular CM $\mathbf{R}$ with $n_r=(n^2-n)/2$ and $\mathbf{d}^{\mathsf{T}}=[d^{(2)},d^{(3)},\ldots,d^{(n)}]\in\mathbb{R}^{n-1}$. The vector $\boldsymbol{\theta}_a\in\mathbb{R}^{n_{\theta_a}}$ with $n_{\theta_a}=n(n+3)/2$ collects parameters $q_1^{(i)}$, $q_2^{(i)}$ and $R_{ij}$.

Goal: Given a model $(\mathbf{F},\mathbf{H},T_s)$ and an available sequence of measurements $\mathbf{z}_{0:N}$, the unknown model parameters $\boldsymbol{\theta}$ are estimated as

where $J:\mathbb{R}^{n_\theta\times(N+1)n_z}\mapsto\mathbb{R}^{++}$ is a designer-selected criterion.

Identification Methods

In this section two methods for estimating the model parameters are presented. Contrary to the standard maximum likelihood methods, both proposed methods aim to avoid solving large non-linear optimization problems by using a suitable transformation of available measurements.

Method based on Hadamard covariance (HCOV)

The basic idea is to use HCOV, whose unique properties enable (i) transformation into a linear optimization problem and (ii) elimination of drift effects, allowing the drift estimation problem to be treated separately.

HCOV Transformation

The HCOV can be defined for measurement sequences $\{z_{i,k}\}_{k=0}^{\infty}$ and $\{z_{j,k}\}_{k=0}^{\infty}$ as

where $\tau_m=mT_s\in\mathbb{R}^{++}$ is an averaging time with $m\in\mathbb{N}$. The ensemble model (4)–(5) and the theory of stochastic processes can be employed (see Appendix A) to get an analytical expression that relates the model parameters $\boldsymbol{\theta}$ to the HCOV $\sigma_{i,j}^{2}(\tau_m)$. If $i\ne j$, the HCOV is given as

If $j=i$ the HCOV reduces to Hadamard variance (HVAR) and it is given as

The expressions (15)–(16) show that for a fixed $\tau_m$, the HCOV/HVAR are linear functions of all model parameters except for the drifts on which they are independent.

Regression model

Assuming that HCOV estimates $\hat{\sigma}_{i,j}^{2}\in\mathbb{R}^{t}$ are computed using (14) where the expectation operator is evaluated as the sample mean with measurements $\mathbf{z}_{0:N}$ for $t$ different averaging times $\overline{\boldsymbol{\tau}}=[\bar{\tau}_1,\bar{\tau}_2,\ldots,\bar{\tau}_t]\in\mathbb{R}^{t}$, where $\bar{\tau}_p=m_pT_s$ with $m_p\in\{1,2,\ldots,N\}$. Then, it is possible to use (15)–(16) to write the following linear regression model

where $\mathbf{z}_a\in\mathbb{R}^{n_{z_a}}$ is a vector given as

with $n_{z_a}=n(n-1)t/2$. The regression matrix $\boldsymbol{\Phi}\in\mathbb{R}^{n_{z_a}\times n_{\theta_a}}$ is defined as

where

and $n_f=n_z(n_z-1)/2$. The HCOV estimation error $\boldsymbol{\epsilon}\in\mathbb{R}^{n_{z_a}}$ has zero mean. Since the true CM of $\boldsymbol{\epsilon}$ is difficult to obtain, it is approximated by a diagonal CM $\mathbf{R}_{\theta_a}\in\mathbb{R}^{n_{z_a}\times n_{z_a}}$ where each diagonal element represents the corresponding variance. Specifically, the variance of HCOV estimate $\hat{\sigma}_{i,j}^{2}(\bar{\tau}_p)$ can be obtained [23] from the properties of the Wishart distribution or as the variance of the sample covariance, where it is assumed that the signals are jointly Gaussian. Then, it can be written

This process involves two key approximations: substituting true (co)variances by their estimated values and reducing the degrees of freedom $\nu$ to account for correlations in the measurement sequence. As $\nu$ depends on the underlying noise process [24], adopting a random-walk frequency modulation (RWFM) model serves as a conservative choice with

The estimate of the unknown parameters $\boldsymbol{\theta}_a$ is obtained by solving the linear minimization problem (13) with the criterion

Note that the number of averaging times $t$ and their particular values represent design parameters.

3. Drift Estimation: Now let us return to the problem of drift estimation. The several approaches can be discussed:

• A regression analysis on the phase time-deviation sequence or its first and second differences using quadratic, linear, or constant models, respectively [24].

• Modification of HCOV method, where HCOV is replaced by Allan covariance [22]. This approach leads to a nonlinear optimization problem and, moreover, does not enable determination of the drift sign.

• Kalman filter based estimation, where the state-space model (1)–(5) can be re-casted such that the second order clock model (1)-(2) is replaced by the third order model, where the third state is the unknown constant drift.

• Measurement difference method (MDM) approach (see below § 4-B4).

In the considered scenario of long and correlated data records, estimation of the drift $d^{(i+1)}$ is performed by the first option. Concretely, a regression analysis is applied to the second difference of the $i$-th component of $\boldsymbol{z}_k$

Then, for given the state-space model (1)–(2), the appropriate regression model is

Using the analogical procedure from Appendix A, it can be easily shown that $\{\varepsilon_{i,k}\}$ is a moving average process of order 2 and satisfies

It follows that the original strongly correlated integration process $\boldsymbol{z}_k^{(i)}$ was transformed into a process with a very short correlation length, which significantly simplifies the following linear optimization problem

where $R_d^{(i)}\in\mathbb{R}^{(N-2)\times(N-2)}$ is an error CM. It has a band sparse structure with the $(I,J)$ term given as

where $\delta_{i,j}$ is a Kronecker delta operator. Note that (27)–(29) can be realized by replacing the true noise parameters with their estimates. As a result, the drift estimates depend on $\hat{\theta}_a$. Finally, the HCOV method is summarized in Algorithm 1.

Algorithm 1 The estimation algorithm based on HCOV

Input: $\boldsymbol{z}_{1:N}$, $\bar{\tau}$,

1: for $i=1,\ldots,n_z$ and $j=i,\ldots,n_z$ do

2: Estimate $\hat{\sigma}_{i,j}^2(\bar{\tau})$ as in (14)

3: Compute $\operatorname{var}(\hat{\sigma}_{i,i}^2(\bar{\tau}))$ and $\operatorname{var}(\hat{\sigma}_{i,j}^2(\bar{\tau}))$ as in (22)

4: Store the results from Steps 2–3

5: end for

6: Construct $\boldsymbol{z}_a$, $R_{\theta_a}$ and $\Phi_a$ as in (17)–(23).

7: Estimate $q_1^{(i)}$, $q_2^{(i)}$, $R_{ij}$ as in (13), (24).

8: for $i=1,\ldots,n_z$ do

9: Construct $\Delta^2\boldsymbol{z}_k^{(i)}$, $R_d$ as in (25), (27)–(29) and (31)

10: Estimate $d^{(i+1)}$ using the procedure (30)

11: end for

Output: Estimated parameters $\hat{\theta}$ and their CMs.

Measurement Difference Method

The unique properties of MDM stem from the definition of the residue, which can be calculated from available measurements and shown to be a linear function of the immeasurable state and measurement noises [25].

1. The Sampling Period and Resampled Model: Based on an initial analysis of the measurements and insights from the ACOV/HCOV techniques, it can be shown that the sampling period affects the extent to which individual parameters are reflected in the measurements and, consequently, the quality of the estimates of these parameters. This phenomenon is presented in Section VI-B10. Therefore, let’s define the resampled model of ensemble of $n$ clocks from Section IV as

for $k=0,1,\ldots,(\bar{N}-1)$, where $\bar{N}=([N/m]+1)$ is a total number of downsampled measurements, $\tau_m=mT_s\in\mathbb{R}^{+ +}$ is a resampling period² in seconds, $m\in\mathbb{N}$ determines the multiple of the original sampling period $T_s$, and the relations also hold:

2. Augmented Measurement Vector: Let’s define the augmented measurement vector as [25]

for $k=0,\ldots,(\bar{N}-L)$, where $L\ge1$ is a user-defined parameter, $\boldsymbol{Z}_k\in\mathbb{R}^{Ln_z}$, $\boldsymbol{W}_k\in\mathbb{R}^{(L-1)n_x}$, $\boldsymbol{V}_k\in\mathbb{R}^{Ln_z}$ $\mathcal{O} \in\mathbb{R}^{Ln_z \times n_x}$, $\Gamma\in\mathbb{R}^{Ln_z \times(L-1)n_x}$, $\mathbf{U} \in\mathbb{R}^{(L-1)n_x}$, and $\Upsilon_{\mathrm{d}} \in\mathbb{R}^{(L-1)n_x \times(n-1)}$ are defined as

Residue Definition

The augmented measurement vector $\mathbf{Z}_k$ is known and depends on the unknown state and measurement noises whose mean values and covariance matrices (CMs) are sought. To eliminate the state $\mathbf{x}_k$, define a non-zero annihilation matrix $\operatorname{am}(\mathcal{O}) \in\mathbb{R}^{n_{a\mathcal{O}} \times Ln_z}$ of the matrix $\mathcal{O}$ such that

where $n_{a\mathcal{O}}=Ln_z-\operatorname{rank}(\mathcal{O})$ and $\operatorname{rank}(\mathcal{O})$ denotes a rank of the matrix $\mathcal{O}$. Then, the residue $\widetilde{\mathbf{Z}}_k \in\mathbb{R}^{n_{a\mathcal{O}}}$ is defined as

where $\mathcal{A}\triangleq\operatorname{am}(\mathcal{O})[\Gamma,\mathbf{I}_{Ln_z}] \in\mathbb{R}^{n_{a\mathcal{O}}\times n_{\mathcal{E}}}$, $\mathcal{B}\triangleq\operatorname{am}(\mathcal{O})\Gamma\Upsilon_{\mathrm{d}} \in\mathbb{R}^{n_{a\mathcal{O}}\times n_z}$, $\boldsymbol{\mathcal{E}}_k\triangleq\begin{bmatrix}\mathbf{W}_k^{\mathsf{T}}&\mathbf{V}_k^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{n_{\mathcal{E}}}$ with $n_{\mathcal{E}}=(L-1)n_w+Ln_v$.

Frequency Drift Identification

Consequently, the mean of the residue $\widetilde{\mathbf{Z}}_k$ is a linear function of the unknown frequency drift $\mathbf{d}$, i.e.,

Although the residue mean $\mathcal{R}_{\widetilde{\mathbf{Z}}}$ is unknown, its sample-based estimate $\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}}=\frac{1}{N-L+1}\sum_{k=0}^{N-L}\widetilde{\mathbf{Z}}_k$ can be computed from measurements and (38) can be modified as

The drift estimates can be obtained, assuming full-rank matrix $\mathcal{B}$, by the least squares (LS) method from (39) as follows

The CM of the estimate $\operatorname{cov}(\widehat{\mathbf{d}})$ depends on the resampling period $\tau_m$; a discussion of this is later in § 4-B10. To obtain the resampling period that leads to minimum CM of the estimate, the estimate of the residue $\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}}$ and, consequently, the drift $\widehat{\mathbf{d}}$ is split into $S>1$ mutually-uncorrelated⁴ sub-estimates as

for $p=1,2,\ldots,S$ where $M=\left\lfloor\frac{N-L+1}{S}\right\rfloor$. The sample-based estimate of the covariance of the drift estimate is computed as

The estimate of the covariance (42) can be computed for $t$ different resampling period $\bar{\tau}=[\bar{\tau}_1,\bar{\tau}_2,\ldots,\bar{\tau}_t]\in\mathbb{R}^{t}$, where $\bar{\tau}_p=m_pT_s$ with $m_p\in\left\{1,2,\ldots,\left\lfloor\frac{N}{L(S+1)-2}\right\rfloor\right\}$, and the one for which the covariance is minimal, in a certain sense, can be used to estimate the drift $\widehat{\mathbf{d}}$ (40). Furthermore, the CM of the estimated drift can be estimated from (42) using the relation

where $C(\bar{N},L)=\frac{1}{\bar{N}}+\sum_{i=1}^{\min(\bar{N}-1,L-1)}\frac{2(\bar{N}-i)}{\bar{N}^{2}}\frac{L-i}{L}$ is described in the Appendix B.

Relation of Residue and Noises Covariances

The residue definition (37) leads to the following key properties: i) $\widetilde{\mathbf{Z}}_k$ is a linear function of the unavailable state and measurement noises (and drift), and ii) it can be computed from the measurements.

Consequently, the residue correlation is a linear function of the CMs of the state and measurement noises, i.e.,

where $\mathcal{A}^{\otimes2}\triangleq\mathcal{A}\otimes\mathcal{A}$ denotes the Kronecker power, $\mathcal{R}_{\widetilde{\mathbf{Z}}^{2}}\in\mathbb{R}^{n_{a\mathcal{O}}^{2}}$ is the vectorised residual correlation containing the CM elements of the residue $\widetilde{\mathbf{Z}}_k$, $\mathcal{R}_{\mathcal{E}_2}\in\mathbb{R}^{n_{\mathcal{E}}^{2}}$ contains elements of the CM of the augmented state and measurement noise vector $\boldsymbol{\mathcal{E}}_k$ and will be specified below.

Noise CM Parametrisation via Structure-Defining Matrices

The noise covariance $\mathcal{R}_{\mathcal{E}_2}$ contains multiple copies of unique elements of the noise CMs $\bar{\mathbf{Q}}$ and $\mathbf{R}$, due to either the construction of the residue, application of the Kronecker algebra, or the CM symmetry. For efficient estimation of the noise CMs, it is necessary to express the unique elements of the noise CMs $\bar{\mathbf{Q}},\mathbf{R}$ from the vector $\mathcal{R}_{\mathcal{E}_2}$ in (44). An appealing parametrisation, introduced in [26], defines $\bar{\mathbf{Q}},\mathbf{R}$ using a set of known (or user-defined) structure-defining matrices $\mathbf{B}_{q_1}^{(i)}\in\mathbb{R}^{2n\times2n}$, $\mathbf{B}_{q_2}^{(i)}\in\mathbb{R}^{2n\times2n}$ and $\mathbf{B}_{r}^{(i)}\in\mathbb{R}^{n_z\times n_z}$ and unknown (and sought) parameters $\mathbf{q}_1,\mathbf{q}_2,\mathbf{r}$ as

where $\mathbf{B}_{q_1}^{(i)}$, $\mathbf{B}_{q_2}^{(i)}$ and $\mathbf{B}_{r}^{(i)}$ are defined as

where $\mathcal{I}(j)$ is a function that returns a zero matrix of dimension $n_z\times n_z$ containing ones only at the positions where the $i$-th unique element $r_i$ of the matrix $\mathbf{R}$ is located. For example, for the first unique element $r_1=R_{11}$, the matrix is $\mathcal{I}(1)=\operatorname{diag}([1\ \mathbf{0}_{1\times(n_z-1)}])$, and the element $r_{n_z+2}=R_{13}=R_{31}$, the matrix is $\mathcal{I}(n_z+1)=\operatorname{diag}\left(\left[\begin{array}{ccc}0&0&1\\0&0&0\\1&0&0\end{array}\right],\mathbf{0}_{n_z-3\times n_z-3}\right)$.

Residue Correlation

Then, the noise CMs parametrisation (46) allows to incorporate a priori knowledge of the noise structure using the minimum number of noise CM parameters gathered in the vector $\boldsymbol{\theta}_a$. The vectors $\mathcal{R}_{\mathcal{E}_2}$ and $\boldsymbol{\theta}_a$ are related as follows

where the notation $(\cdot)_{\mathrm{V}}$ stands for vectorisation of the matrix and the replication matrices $\mathbf{\Upsilon}_{q_1}\in\mathbb{R}^{n_a^2\mathbin{\times}n}$, $\mathbf{\Upsilon}_{q_2}\in\mathbb{R}^{n_a^2\mathbin{\times}n}$, and $\mathbf{\Upsilon}_{r}\in\mathbb{R}^{n_a^2\mathbin{\times}n_r}$ are defined as

Residue Correlation Unification

Similarly to the noise covariance $\mathcal{R}_{\mathcal{E}_2}$ in (49), the residue correlation $\mathcal{R}_{\widetilde{\mathbf{Z}}_2}$ also contains redundant copies of the same elements. Therefore, unique elements of the vector $\mathcal{R}_{\widetilde{\mathbf{Z}}_2}$ gathered in the vector denoted as $\mathcal{R}_{\widetilde{\mathbf{Z}}_2}^{\mathrm{U}}\in\mathbb{R}^{n_{a_O}(n_{a_O}+1)/2}$ are selected using a unification matrix denoted as $\widetilde{\mathbf{\Xi}}_{\widetilde{\mathbf{Z}}_2}\in\mathbb{R}^{(n_{a_O}(n_{a_O}+1)/2)\times n_{a_O}^2}$. Therefore, the following relationship holds

Noise Parameter Estimation

Consequently, the correlation of the residue is a linear function of the noise parameters and equations (47) and the kron-square of the state noise mean vector (drift) (49), i.e.,

The residue correlation $\mathcal{R}_{\widetilde{\mathbf{Z}}_2}$ (59) is unknown. However, its sample-based estimate $\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}_2}^{\mathrm{U}}=\widetilde{\mathbf{\Xi}}_{\widetilde{\mathbf{Z}}_2}\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}_2}=\frac{\widetilde{\mathbf{\Xi}}_{\widetilde{\mathbf{Z}}_2}}{\overline{N}-L+1}\sum_{k=0}^{\overline{N}-L}\widetilde{\mathbf{Z}}_k^{\otimes2}$ can be computed and (50) can be modified as

Analogously to the drift $\mathbf{d}$ and its estimate (45), the noise parameters and squared drift can directly be estimated by the LS method from (51). However, in this straightforward solution the MSE of a subset of estimated parameter will be significantly higher than the MSE of the remaining parameters due to the “uneven” influence of the sampling period, that is analyzed Section VI-B10. To equalize the MSE for all estimated parameters, we, therefore, propose a two-step procedure for sequential noise parameters estimation based on two purposely selected (down)-sampled periods.

The first step, it can be shown that for a small resampling period, e.g., $\tau_m\leq100\mathrm{s}$, the influence of the noise parameters $\mathbf{q}_2$ and $\mathbf{d}$ is negligible, and the equation (61) can be considered in the following form

The noise parameters $\mathbf{q}_1$ and $\mathbf{r}$ estimates can be obtained by LS from (52) as follows

The second step, the influence of the noise parameters $\mathbf{q}_1$, $\mathbf{r}$, and $\mathbf{d}$ can be eliminated from equation (51) by using the available estimates $\widehat{\mathbf{q}}_1$, $\widehat{\mathbf{r}}$ (63) and $\widehat{\mathbf{d}}$ (45), and subsequently estimate the last remaining parameter $\mathbf{q}_2$ as follows

The CM of the estimates $\widehat{\mathbf{q}}_1$, $\widehat{\mathbf{r}}$, and $\widehat{\mathbf{q}}_2$ depends on the resampling period $\tau_m$. The technique for obtaining the resampling period $\tau_m$ that minimizes the CM of the estimates (similar to the estimation of the drift (40)) is to use $S$ mutually uncorrelated sub-estimates as

for $p=1,2,\ldots,S$, where $\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}_2}^{p,\mathrm{U}}=\widetilde{\mathbf{\Xi}}_{\widetilde{\mathbf{Z}}_2}\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}_2}^{p}$, $\widehat{\mathcal{R}}_{\widetilde{\mathbf{Z}}_2}^{p}=\frac{1}{\overline{M}-L+1}\sum_{k=0}^{\overline{M}-L}\widetilde{\mathbf{Z}}_{k+(p-1)\overline{M}}^{\otimes2}$ and $\overline{M}=\left\lfloor\frac{\overline{N}-L+1}{S}\right\rfloor$.

Using the sub-estimates (55), the appropriate resampling period can be found for the estimate of the pair $\left[\widehat{\mathbf{q}}_1\ \widehat{\mathbf{r}}\right]$ (53) as well as for the estimate of $\widehat{\mathbf{q}}_2$ (54). As a consequence, each of the chosen resampling period minimizes the CM of the respective estimate error, which can be calculated as

where

Influence of Resampling Period on Noise Parameter Estimation

Notes, and Algorithm: Let’s focus on observations that determined why individual estimates are calculated in these particular ways. In this section consider a set of AHM clocks as described in § 6-A. The observations presented below are illustrative; however, they have been valid for various model parameters as well.

Observation I: To estimate $\mathrm{d}$ (45), the choice of the resampling period $\tau_m$ is crucial, as it can significantly affect the variance of the estimate, as shown in Figure 1. However, it can be seen the wide range between $\tau_m=60\,\mathrm{s}$ and $\tau_m=41000\,\mathrm{s}$ where the STD of the estimate differs from the minimum by 1%.

Observation II: The appropriate resampling periods $\tau_m$ determined based on sub-estimates by (47) and (57) are very close to the resampling periods that lead to the minimum variance of the estimate calculated both theoretically and using Monte Carlo (MC) simulations i.e., $\tau_m=4800\,\mathrm{s}$. This can be seen in Figure 1 for the frequency drift identification.

Observation III: The error of the estimated CM of the estimate based on (48) and (56) (compared to the CM estimated based on MC simulations) for different value of measurements $N$ (with $T_s=5\,\mathrm{s}$) is shown in Figure 2. The average error of CMs for $\hat{q}_1$ and $\hat{r}$ is 6%, independently of the total number of measurements $N$ with $T_s=5\,\mathrm{s}$. In contrast, the error of CMs for $\hat{q}_2$ and $\hat{\mathrm{d}}$ decreases as the value of $N$ increases.

Observation IV: The resampling time $\tau_m$ determines the ratio of the contributions of the unknown parameters $q_1$, $q_2$, $r$, and $\mathrm{d}$ (60) to the correlation value $\widetilde{R}_{\widetilde{Z}_2}$, as shown in Figure 3.

Therefore, for a small resampling period $\tau_m$ (Area 1), the correlation $\widetilde{R}_{\widetilde{Z}_2}$ is formed exclusively by the parameters $q_1$ and $r$. As a result, the parameters $q_1$ and $r$ can be estimated by (63) individually, without knowledge of the remaining parameters $q_2$ and $\mathrm{d}$. For longer resampling period $\tau_m$ (Area 2), the correlation $\widetilde{R}_{\widetilde{Z}_2}$ is composed of all parameters. However, since the estimates $\mathrm{d}$ (45) and $q_1$, $r$ (63) are already available, their effects can be removed from the equation, and the remaining parameter $q_2$ can be estimated by (64).

Observation V: The appropriate resampling period for both $q_1$, $r$ parameters and $q_2$ parameter yields the minimum variance of these estimates (see Observation II), corresponds to the resampling period at which the parameters have a significant contributions (see Observation IV) i.e., $\tau_m=8\,\mathrm{s}$ for $q_1$, $r$ and $\tau_m=36000\,\mathrm{s}$ for $q_2$.

Note I: The noise can be arbitrarily distributed, and the assumption of Gaussian distributed noises are not necessary.

Note II: The technique for finding the appropriate resampling period from set $\bar{\tau}$ can also be applied to $\ell$ different values of $L$ i.e., $\bar{L}=[\bar{L}_1,\bar{L}_2,\ldots,\bar{L}_{\ell}]\in\mathbb{R}^{\ell}$, where $\bar{L}_p\in\mathbb{N}$, to find the appropriate one for which the noise parameters (as a whole) have minimum variance. The appropriate value is $L=8$ as is shown in Figure 4.

Algorithm 2 The MDM estimation algorithm

Input: $z_{1:N}$, $\bar{\tau}$, $\bar{L}$, $S$

1. for $L\in\bar{L}$ do

2. for $\tau_m\in\bar{\tau}$ do

3. Compute residue $\widetilde{Z}_k$ (42).

4. Estimate $\hat{\mathbf{d}}^p$ as in (46) and $\begin{bmatrix}\hat{q}_1^p\\\hat{r}^p\end{bmatrix}$ (65).

5. Estimate $\operatorname{cov}(\hat{\mathrm{d}})$ as in (48) and $\operatorname{cov}\left(\begin{bmatrix}\hat{q}_1\\\hat{r}\end{bmatrix}\right)$ (67).

6. end for

7. Estimate $\hat{\mathrm{d}}$ (45) for such a $\tau_m$ that $\operatorname{cov}(\hat{\mathrm{d}})$ is minimal.

8. Estimate $\begin{bmatrix}\hat{q}_1\\\hat{r}\end{bmatrix}$ (63) for such a $\tau_m$ that $\operatorname{cov}\left(\begin{bmatrix}\hat{q}_1\\\hat{r}\end{bmatrix}\right)$ is minimal.

9. for $\tau_m\in\bar{\tau}$ do

10. Estimate $\hat{q}_2$ (66) with $\hat{q}_1$, $\hat{r}$, and $\hat{\mathrm{d}}$ from steps 7–8.

11. Estimate $\operatorname{cov}(\hat{q}_2)$ (68).

12. end for

13. Estimate $\hat{q}_2$ (64) for such a $\tau_m$ that $\operatorname{cov}(\hat{q}_2)$ is minimal.

14. end for

15. Through all $L$, pick estimates $\begin{bmatrix}\hat{q}_1\\\hat{r}\end{bmatrix}$, $\hat{\mathrm{d}}$, and $\hat{q}_2$ from steps 7, 8, and 13, for which their CMs are minimal.

Output: Estimated parameters $\hat{\boldsymbol{\theta}}$ and their CMs.

Results

This section presents illustration and comparison of both identification methods for simulated and real data.

Simulated data

The ensemble consists of four active hydrogen maser (AHM) clocks with true parameters listed in tbl:I and tbl:II. Other parameters considered are: $T_s=5\,\mathrm{s}$, $N=$ $6.312\mathrm{E}6$ (i.e., this corresponds to a one-year data period). The parameter for HCOV is $t=20$, where averaging times $\bar{\tau}_p$ are evenly spaced in a log-space ranging 5 s to $3.15\mathrm{E}7$ s and the parameters for MDM are $S=100$, $t=400$, where resampling period $\bar{\tau}_p$ are evenly spaced in a log-space ranging from 5 s to $5\cdot\frac{6\,312\,000}{101\cdot\bar{L}^{-2}}$ s, and $\bar{L}=[4,5,6,7,8,9,10,11]$ i.e., $\ell=8$.

The resulting average and standard deviations of the estimates for both methods (HCOV and MDM) based on 100 Monte Carlo simulations are summarized in Table I (PDF p. 9) and Table II. Both methods seems to provide noise parameter estimates of almost identical quality. A more significant difference can be observed in the uncertainty estimation for parameter $q_2$, where the HCOV method provides optimistic CM, whereas MDM provided CM corresponding to those based on MC simulations. The graphical representation of the results is provided in Figure 5 and Figure 6. The figures present the Allan variance (AVAR) of the individual clocks computed analytically using the true parameters (blue line) and using the average values of the estimated parameters from MC simulations (dashed red line). The confidence intervals (CIs) of the AVAR are shown as shaded areas. The AVARs computed using MC average of parameter estimates are close to the AVARs computed using true parameters. In both cases, a wide confidence interval is obtained for the pivot clock for averaging times greater then $1\mathrm{E}4$ s. This behavior can be attributed to (i) the assumption of zero drift for the pivot clock, (ii) the low value of the parameter $q_2$ for the pivot clock, and (iii) the limited length of the measurement record. Although the width of a confidence intervals could be further decreased by increasing the duration of measurement, a substantial reduction

		$q_1^{(i)}$ ($\times1\mathrm{E}{-27}$)				$q_2^{(i)}$ ($\times1\mathrm{E}{-35}$)				$d^{(i)}$ ($\times1\mathrm{E}{-21}$)
	$i$	1	2	3	4	1	2	3	4	2	3	4
	true	1	1.5	5	7	0.1	2	1.5	2.5	8	7.5	3
HCOV	MC-Based $\mathbb{E}(\hat{\theta})$	1	1.5	5.001	7	0.133	1.869	1.328	2.197	7.982	7.414	3.051
	MC-Based $\operatorname{std}(\hat{\theta})$	0.003	0.004	0.006	0.008	0.048	0.251	0.206	0.448	0.761	0.810	0.916
	Method provided $\operatorname{std}(\hat{\theta})$	0.002	0.004	0.007	0.009	0.043	0.093	0.092	0.131	0.815	0.696	0.876
MDM	MC-Based $\mathbb{E}(\hat{\theta})$	1	1.5	4.999	7	0.089	2.012	1.495	2.488	8.137	7.431	3.187
	MC-Based $\operatorname{std}(\hat{\theta})$	0.004	0.004	0.008	0.011	0.070	0.115	0.126	0.199	0.902	0.688	0.991
	Method provided $\operatorname{std}(\hat{\theta})$	0.004	0.004	0.008	0.011	0.077	0.153	0.151	0.199	0.842	0.755	0.944

		$(\times1\mathrm{E}{-27})$	$R_{11}$	$R_{12}$	$R_{13}$	$R_{22}$	$R_{23}$	$R_{33}$
	true	9	6	5	8.7	4	9.5
HCOV	MC-Based $\mathbb{E}(\hat{\theta})$	8.999	5.999	4.999	8.698	3.999	9.497
	MC-Based $\operatorname{std}(\hat{\theta})$	0.013	0.011	0.011	0.018	0.013	0.025
	Method provided $\operatorname{std}(\hat{\theta})$	0.01	0.008	0.008	0.016	0.009	0.020
MDM	MC-Based $\mathbb{E}(\hat{\theta})$	9.001	6.001	4.999	8.703	4	9.500
	MC-Based $\operatorname{std}(\hat{\theta})$	0.014	0.013	0.012	0.023	0.010	0.028
	Method provided $\operatorname{std}(\hat{\theta})$	0.013	0.012	0.013	0.021	0.013	0.027

would require measurement duration on the order of years. Figure 7 and Figure 8 illustrate the relationship between the quality of individual noise parameter estimates and the number of measurement $N\in[1\mathrm{E}5,1\mathrm{E}8]$ (equivalent to the total observation period ranging from 6 days to 3 years). The estimation quality is quantified by a t-statistic, defined as $t_{stat}=\mathbb{E}(a)/\operatorname{std}(a)$. Both methods yield similar trends for most of the estimated parameters. Statistically significant estimates are obtainable for all parameters $q_1$ and $R_{ij}$ (assuming a confidence level of 0.95 where $t\geq2$) even for the minimum duration of $N=1\mathrm{E}5$. In contrast, for certain $q_2$ (clk1) parameters and drift (clk4), the estimates achieve statistical significance only when $N\geq2\mathrm{E}6$ (equivalent to a measurement period of approximately three months). It is consistent with previous studies [15], in which the estimation process may produce negative values for $q_2$, a result that is physically non-realizable.

The simulation results validate HCOV and MDM as efficient identification methods offering comparable estimation quality.

The relevant parameters in MDM, as the result of the calculation, are as follows: the resampling periods $\tau_m$ that minimize the covariance of the estimates (§ 4-B10 Observation II) are $\tau_m=5\mathrm{s}$ for $q_1$ and $r$ (see Figure 3); $\tau_m=12180\mathrm{s}$ for $q_2$ (see Figure 3); $\tau_m=9070\mathrm{s}$ for the estimate of $d$ (see Figure 1) and $\bar{L}=8$ (see Figure 4).

Real data

The historical records of measurements were provided by the ESTEC timing laboratory during year 2022 as part of the ESA project NAVISP-EL1-056. The ensemble consists of triplet AHM atomic clocks iMaser3000 [27]. The nominal

The sampling period of the phase difference measurements is $5\ \mathrm{s}$ and the dataset contains $N = 1.59\mathrm{E}6$ sampling instants (i.e., this corresponds to a three-months data period). The parameter for HCOV is $t = 20$, where averaging times $\bar{\tau}_p$ are evenly spaced in a log-space ranging $5\ \mathrm{s}$ to $7.95\mathrm{E}6\ \mathrm{s}$ and the parameters for MDM are $S = 100$, $t = 400$, where resampling period $\bar{\tau}_p$ are evenly spaced in a log-space ranging from $5\ \mathrm{s}$ to $5\cdot\left[\frac{1\,589\,695}{101\cdot\bar{L}-2}\right]\ \mathrm{s}$, and $\bar{L} = [4, 5, 6, 7, 8, 9, 10, 11]$ i.e., $\ell= 8$.

The obvious measurement outliers were removed from the real measurements before the parameter estimation methods were applied. The AVARs of individual clocks computed using estimated parameters (dashed lines) and AVAR obtained from data-sheet (solid line) are depicted in Figure 9 and 10. Contrary to the simulation scenario, the parameter estimates obtained by the proposed methods show significant discrepancies. For clocks 1 and 2, the parameter estimates differ by up to 20 %, whereas for clock 3 the differences span several orders of magnitude. These discrepancies may be attributed to (i) an insufficient amount of data (the available three-month record may be too short, leading to high estimation uncertainty) (ii) the presence of additional noise types in real measurements that are included in the model, and (iii) anomalies in measurements that were not fully removed during preprocessing. However, the results provide an indication of the achievable estimation accuracy under realistic conditions. This insight is particularly relevant for subsequent use of the ensemble model in combination and detection algorithms, which should be designed to remain robust in the presence of such uncertainties.

The appropriate parameters in MDM are as follows: the resampling periods $\tau_m$ that minimize the covariance of the estimates (§ 4-B10 Observation II) are: $\tau_m = 5$ s for $q_1$ and $r$ (see Figure 3); $\tau_m = 13\,680$ s for $q_2$ (see Figure 3); $\tau_m = 5600$ s

for the estimate of $d$ (see Figure 1) and $L = 5$.

Discussion

Although the motivations and underlying principles of the proposed methods differ, they share the following features:

	$i$ or $ij$	$q_1^{(i)}$ ($\times1\mathrm{E}{-27}$)			$q_2^{(i)}$ ($\times1\mathrm{E}{-36}$)			$d^{(i)}$ ($\times1\mathrm{E}{-21}$)		$R_{ij}$ ($\times1\mathrm{E}{-27}$)
		1	2	3	1	2	3	2	3	11	12	22
MDM	$\hat{\theta}$	2.398	0.990	4.435	0.254	4.420	0.140	7.756	0.727	8.227	6.151	7.729
HCOV	Method provided std$(\hat{\theta})$	0.009	0.020	0.0170	0.233	0.820	0.255	0.769	0.227	0.021	0.022	0.035
MDM	$\hat{\theta}$	1.554	1.085	3.817	5.701	26.654	127.564	7.792	0.844	9.374	7.492	10.133
HCOV	Method provided std$(\hat{\theta})$	0.019	0.013	0.017	5.230	8.110	17.867	2.310	3.626	0.055	0.044	0.058

• Independence on state, linear filtering: Both methods prevent the state dynamics from influencing noise estimates. HVAR can be interpreted alternatively as the output of a specific band-pass filter (rectangular averaging followed by the second-order differentiation) that eliminates the effect of the state (see Appendix A). In MDM, this is achieved through the transformation of the augmented measurement vector (37).

• Noise components decomposition: Figure 3, in the case of the MDM, and the HVAR plot, in the case of the HCOV method, can be seen as visualization of different noise types as a function of the resampling period $\tau_m$. However, while the HCOV method uses regression to fit the HCOV curve at selected points $\bar{\tau}$ using a model with unknown noise parameters, the MDM first systematically separates the individual noise components and then estimates the noise parameters for the appropriate determined resampling periods $\tau_m$.

• Transformation into a linear estimation problem(s): HCOV uses the analytical form of HCOV (15), (16), where unknown noise parameters appear linearly. In MDM, this is achieved by using a combination of a specific residue (37), its correlation (44), and the replication matrices $\Upsilon_d$ (35) and $\Upsilon_{q_1}$, $\Upsilon_{q_2}$, $\Upsilon_r$ (48).

• Resampling: Whereas HCOV and other Allan like statistics can be interpreted as a multiscale resampling variance, MDM uses explicit resampling for both measurements and the model (32).

• Bootstrap technique: Both methods use bootstrap techniques

as statistical resampling methods that estimate the statistical properties of $\hat{\theta}$ by repeatedly generating new datasets through observed data sampling. Note that the samples are selected deterministically in both cases.

Conclusions

The article presents two identification methods for determining ensemble clocks noise parameters. The performance was evaluated on both synthetic and real data. While the accuracy of the parameter estimation is primarily constrained by sample size, both methods prove to be computationally efficient, ensuring reliable characterization of noise characteristics under realistic conditions. The proposed methods yield results of similar quality when applied to synthetic data. In a real data case, the results differ, which is due to their different underlying concepts. The choice of method ultimately depends on whether the user prefers a design model fitted to the HCOV plot or a model based on the systematic decomposition of noise into its individual components.

HCOV (HVAR) derivation

HCOV in (14) can be rewritten using the state-space representations (1)–(5) and with the variables rearrangement as

This can be rewritten using the third difference of variable $g_k$ denoted as $\Delta^3 g_k=g_{k+3m}-3g_{k+2m}+3g_{k+m}-g_k$ as

Considering that the individual ensemble clocks are uncorrelated (as follows from the block diagonal form of the CM $\mathbf{Q}$) and using $\mathbb{E}(x\cdot y)=\operatorname{cov}(x,y)+\mathbb{E}(x)\mathbb{E}(y)$, the equation (59) can be simplified to the form

Before focusing on individual members in equation (60), it is convenient to introduce the relations

where $\mathbf{F}^{m}=\begin{bmatrix}1&mT_s\\0&1\end{bmatrix}$, white Gaussian noises $\mathbf{w}_{m}^{(i)}$, $\widetilde{\mathbf{w}}_{m}^{(i)}$ and $\widetilde{\mathbf{w}}_{m}^{(i)}$ are independent and have identical statistical properties, i.e., $\mathbf{w}_{m}^{(i)},\widetilde{\mathbf{w}}_{m}^{(i)},\widetilde{\mathbf{w}}_{m}^{(i)}\sim\mathcal{N}\left(\boldsymbol{\mu}_{m}^{(i)},\mathbf{Q}_{m}^{(i)}\right)$ with

Note that (61)–(63) follows from the fact that $\mathbf{x}_{k+m}$ can be understood as a change in the sampling period of the dynamic system (4)–(5) by a multiple of $m$.

Now, the individual members of expression (60) can be covered. The first one is concerned with $\Delta^{3}x_{1,k}^{(i)}$ that can be rewritten using equations (61)–(63) as follows

It is clear that $\left(\mathbf{F}^{3m}-3\mathbf{F}^{2m}+3\mathbf{F}^{m}-\mathbf{I}\right)$ is a zero matrix. Thus, $\Delta^{3}x_{1,k}^{(i)}$ depends solely on the noise characteristics as

One can easily obtain mean and variance of $\Delta^{3}x_{1,k}^{(i)}$ as

Next, let’s focus on $\Delta^{3}v_{i,k}$. Since

and $v_k$ is an independent sequence, for mean and variance of $\Delta^{2}v_{i,k}$ can be stated

Finally, the covariance $\operatorname{cov}(\Delta^{3}v_{i,k},\Delta^{3}v_{j,k})$ can be deduced as

Using (68), (69), (71) and (72), HCOV $\sigma_{i,j}^{2}(\tau_m)$ from (60) results in (15). Finally, HVAR $\sigma_{i,i}^{2}(\tau_m)$ given by (16) can be obtained analogously by the procedure (14) and (58)–(72).

Sample-Based Covariance Estimation

Consider the following linear equation

where $\alpha_k$ is a known measurement, $\mathbf{D}$ is a known matrix, $\boldsymbol{\omega}$ is a vector of unknown parameters, and $\boldsymbol{\pi}_k$ is zero mean noise with CM $\operatorname{cov}(\boldsymbol{\pi}_{k+r},\boldsymbol{\pi}_k)=\mathbf{\Pi}_r\ne0$ for $\operatorname{abs}(r)<L$ and $\mathbf{\Pi}_r=0$ for $\operatorname{abs}(r)\ge L$.

The estimate of $\boldsymbol{\omega}$ using LS for a total of $N$ available measurements is the following

The CM of the estimate (74) is determined by the following

It can be shown that the matrix inequality $2\mathbf{\Pi}_0\succeq(\mathbf{\Pi}_1+\mathbf{\Pi}_1^{\mathsf{T}})\succeq\cdots\succeq(\mathbf{\Pi}_L+\mathbf{\Pi}_L^{\mathsf{T}})=0$ holds i.e., the dependence $\operatorname{cov}(\boldsymbol{\pi}_{k+r},\boldsymbol{\pi}_k)$ decreases as $\operatorname{abs}(r)$ increases from 0 to $L$. Based on this, the user chose a linearly decreasing dependency as the most appropriate compromise between the absence of dependency (i.e., $\mathbf{\Pi}_0\ne0$ and $\mathbf{\Pi}_1=\mathbf{\Pi}_2=\cdots=\mathbf{\Pi}_L=0$) and full dependency (i.e., $2\mathbf{\Pi}_0=(\mathbf{\Pi}_1+\mathbf{\Pi}_1^{\mathsf{T}})=\cdots=(\mathbf{\Pi}_{L-1}+\mathbf{\Pi}_{L-1}^{\mathsf{T}})$ and $\mathbf{\Pi}_L=0$). Therefore the estimate $(\mathbf{\Pi}_r+\mathbf{\Pi}_r^{\mathsf{T}})=2\mathbf{\Pi}_0\frac{L-r}{L}$ for $0\le\operatorname{abs}(r)\le L$ is used. Using this estimate in (75b), the estimate of the CM of the estimate can be obtained by

The relation (76b) is useful when the estimate of the CM of the estimate $\operatorname{cov}(\widehat{\boldsymbol{\omega}}_N)$ for a total of $N$ measurements is available, and the CM of the estimate $\operatorname{cov}(\widehat{\boldsymbol{\omega}}_M)$ for a total of $M$ measurements is needed. The estimate can be obtained by