# Bias-Free Estimation of the Auto- and Cross-Covariance and the Corresponding
Power Spectral Densities from Gappy Data

Nils Damaschke
Institute of General Electrical Engineering, University of Rostock
Faculty of Computer Science and Electrical Engineering
Albert-Einstein-Straße 2, 18059 Rostock, Germany
<EMAIL_ADDRESS>
Volker Kühn
Institute of Communications Engineering, University of Rostock
Faculty of Computer Science and Electrical Engineering
Richard-Wagner-Straße 31, 18119 Rostock-Warnemünde, Germany
<EMAIL_ADDRESS>
Holger Nobach
Max Planck Institute for Dynamics and Self-Organization
Am Faßberg 17, 37077 Göttingen, Germany
<EMAIL_ADDRESS>
###### Abstract

Signal processing of uniformly spaced data from stationary stochastic
processes with missing samples is investigated. Besides randomly and
independently occurring outliers also correlated data gaps are investigated.
Non-parametric estimators for the mean value, the signal variance, the
autocovariance and cross-covariance functions and the corresponding power
spectral densities are given, which are bias-free, independent of the spectral
composition of the data gaps. Bias-free estimation is obtained by averaging
over valid samples only from the data set. The procedures abstain from
interpolation of missing samples. An appropriate bias correction is used for
cases where the estimated mean value is subtracted out from the data. Spectral
estimates are obtained from covariance functions using Wiener-Khinchin’s
theorem.

##  1 Introduction

In normal operation, measurement instruments usually deliver a continuous
stream of equidistantly spaced samples of an observed, physical quantity.
Various reasons exist, which lead the normal operation to fail. There might be
boundary conditions, under which the measurement system cannot operate. This
includes cases, where the quantity under observation temporarily is not
accessible. For distributed measurement systems, the communication channels
might be temporarily disturbed. Also re-configuration or maintenance
interruption might be necessary for the instrument. In many cases, the
measurement principle includes a signal pre-processing, which may fail under
certain conditions and lead to outliers or gaps in the data stream.

Unfortunately, signal processing algorithms for statistical analysis usually
require a continuous stream of equidistantly sampled data, typically arranged
in blocks of a defined duration. A widely used way to circumvent the problem
caused by missing samples is to interpolate the signal and fill the gaps with
predicted values, see Vio et al. (2000). The interpolation scheme can mime the
statistical properties obtained from the valid parts of the signal and missing
values can be predicted. Then the statistical quantities like the mean value,
the variance, the covariance function or the power spectral density are
derived from the reconstructed signal consisting of a mixture of originally
valid samples and the interpolated ones. Examples for this principle are
Kalman interpolation in audio reconstruction by Niedźwiecki and Cisowski
(1993); Preihs et al. (2012), the adaptive filter-bank approach by Stoica et
al. (2000) or the Karhunen-Loève procedure resp. proper orthogonal
decomposition for gappy data by Everson and Sirovich (1995) or in turbulence
measurements by Venturi and Karniadakis (2004). However, even the best
interpolation in terms of the minimum prediction error, understood as the
minimum mean square error between the interpolated signal and the true signal,
will lead to a significant dynamic error depending on the probability of
invalid data points. The prediction will ideally replace missing samples by
something like the expectation of all possible continuations of the signal at
the respective time instance. Then the interpolation itself would be bias-
free, however, it will suppress parts of the fluctuations of the true signal.
Finally, statistical properties of the partially interpolated signal deviate
from those of the original signal. For rare and short gaps in the data
sequence this may work sufficiently, for more and longer gaps, the error may
easily become unacceptable. Note that this holds for any interpolation scheme,
even for those, which perfectly mimic the spectral composition of the signal
under observation.

Plantier et al. (2012) derive the expected dynamic error particularly for
sample-and-hold interpolation, see also Nobach (2019). 111Plantier et al.
(2012) state that Eq. (14) (in that publication) would be an approximation.
However, it is the explicit and exact formulation for the signal model used by
Nobach et al. (1998), namely the random sampling of a continuous time. The
respective formulation for the equidistant sampling used by Plantier et al.
(2012) is $c_{0}(\alpha)=1-\frac{2}{\alpha}+\frac{2}{\alpha^{2}}$. For an
assumed covariance function of the uninterrupted signal, the covariance
function after interpolation of missing data points is predicted, using the
statistical properties of the occurrence of these missing data points. The
primary covariance function estimated from the interpolated signal then can be
improved by a deconvolution derived from this correspondence. Under ideal
conditions the deconvolution entirely inverts the dynamic error caused by the
interpolation. Therefore, the idea to consider and to invert the influence of
the interpolation step, is a promising method to obtain bias-free estimates of
the covariance function and the spectrum. However, the derivation of the
mapping matrix depends on specific statistical properties of the data gaps,
namely a random occurrence of independent invalid samples. If the statistics
of the data gaps change, then the procedure needs substantial modification.
Therefore, it is no universal solution.

Sujbert and Orosz (2016) use no interpolation. Instead the expected spectrum
of the discrete Fourier transform has been derived for the signal with missing
samples or data gaps of various statistical characteristics. From these
expectations, a deconvolution could have been developed to improve the
estimation of the spectra, similar to the deconvolution after interpolation by
Plantier et al. (2012). Instead, a procedure has been used avoiding the data
gaps by rejecting all values past the first missing sample by means of zero
padding. Depending on the probabilities of the occurrence of the first missing
sample, the spectra get different resolution and the superposition of many
spectra from individual data blocks becomes smeared. However, this method is
very inefficient in using the information available, since significant amount
of valid samples get rejected. It is limited to small amounts of missing
samples anyway, since series of valid samples become too short for increasing
amounts of missing samples.

There exists a wide variety of direct spectral estimators optimized for
spectral estimation from a limited (typically small) number of unevenly
sampled observations of signals Lomb (1976); Scargle (1982); Masry (1983);
Ferraz-Mello (1986); Foster (1996); Mathias et al. (2004); Rivoira and Fleury
(2004); Stoica and Sandgren (2006); Babu and Stoica (2010); Damaschke et al.
(2018). They are widely counted as direct spectral estimators, since the
amplitudes of the spectrum are obtained directly at any frequency from the
sinusoidal fits. All these estimators potentially are able to process also
signals with missing samples. However, the spectral composition of
equidistantly sampled signals with independently missing samples deviates from
random sampling in continuous time, not to mention correlated data gaps, which
the above methods cannot handle accurately. Since Lomb-Scargle’s method is
widely used as a benchmark, it is included in the comparison below to prove it
as biased.

Missing data samples in equidistant data streams have also been investigated
broadly by Jones (1962a, b); Parzen (1963); Scheinok (1965); Bloomfield
(1970); Jones (1971, 1972); Ghazal and Elhassanein (2006); Munteanu et al.
(2016), including also specific cases of correlated data gaps. These
derivations strictly depend on the specific cases of missing data and are not
robust against changes in spectral content of the data gaps. This also holds
for Robinson (1977); Jones (1980); Dunsmuir and Robinson (1981a, b); Robinson
(1985), who performed parametric estimation. This way effective process
identification is possible in a limited search space. However, bias-free
estimation is not possible for unknown or changing spectral content of the
data gaps.

In the present article a non-parametric and bias-free method is introduced,
which is i) simple to realize, ii) efficient in using the available
information and iii) universal for various statistical properties of the data
gaps. It is a combination of three methods a) taken from Jones (1971), except
for any spectral or time windowing, b) using derivations as in Vogelsang and
Yang (2016) adapted to the two-sided autocovariance function including
weighting and transferred to the cross-covariance case (see details in Nobach
(2023)) and c) mentioned only briefly as an appropriate means for spectral
estimation in Bartlett (1948). This combination allows bias-free estimation of
the variance of the signal with data gaps, its covariance function and the
corresponding power spectrum, independent of the spectral content of the data
gaps, including those cases, where the mean value is estimated and subtracted
out from the data. No pre-knowledge about the characteristics of the signal
under investigation is needed. The procedure performs the averaging as it
occurs in the definition of the statistical properties, taking into account
only valid data samples and masking the invalid ones. Interpolation of missing
portions of the signal is not required. Neither interpolation of the data gaps
nor any kind of reconstruction of the signal are intended. Regularization
attempts or model-based estimation are not taken into account.

The invalid samples are assumed a priori known, given by an additional flag
(weight) for each data sample indicating whether the sample is valid or
invalid. This way, the samples with invalid values can be excluded from the
averaging process. It is assumed here, that the data gaps (randomly and
independently occurring outliers as well as correlated data gaps) are not
containing information from the data and that they occur randomly. While the
first assumption is essential for the sufficiency of the data ensemble, no
requirements are involved concerning the statistical characteristics or the
occurrence of the data gaps. In extreme cases, the gaps could also occur with
a static pattern in consecutive realizations of the signal. The only
requirement to the occurrence of data gaps is that pairs of data points exist
at all required time lags for estimating the covariance function. The programs
used here are available as supplementary material to the present article.

Note, that the introduced processing methods are suitable for signals from
stationary stochastic processes only and for their statistical analysis. They
are not suited for signals with time-dependent statistical properties, e.g.
single pulses. For those, missing information cannot be restored by the
methods presented here. The method is just using that information, which is
available. Note further, that the weighted averages lead to different
discrimination of the data ensemble at different lag times. This results in a
correlation matrix, which potentially may violate the non-negative
definiteness. As a consequence, negative values occur in the corresponding
power spectral density. Since the introduced procedures yield bias-free
estimates for both, the covariance function as well as the spectrum, averages
over multiple estimates of these functions will converge towards the true
functions of the underlying process. Furthermore, the procedures are shown to
be consistent. Therefore, the estimated functions will also converge towards
the true functions if applied to single but longer data sets, without losing
information between block boundaries. The ultimate solution of course is
regularization. Since this inevitably introduces a bias in both the covariance
function and the corresponding spectrum, regularization is not considered in
the present article, where bias-free estimation has priority.

##  2 Processing method

Table 1: Nomenclature ${\mathbf{x}}$, ${\mathbf{y}}$ | signals (data sequences) without data gaps
---|---
$x_{i}$, $y_{i}$ |  $i$-th values of the sequences ${\mathbf{x}}$ and ${\mathbf{y}}$ resp.
| (similar use of indices for other sequences)
$t_{i}$ | time instance of the $i$-th sample
$\Delta\tau$ | fundamental sampling interval
$N$ | record length (number of values)
$N_{\mathrm{x}}$, $N_{\mathrm{y}}$ | number of values in the sequences ${\mathbf{x}}$ and ${\mathbf{y}}$ resp.
${\mathbf{z}}$ | signal with data gaps (valid values from signal ${\mathbf{x}}$
| and individual faulty values)
${\mathbf{w}}$ | sequence of weights ($w_{i}\in\\{0;1\\}$), indicating
| validity of values in sequence ${\mathbf{z}}$
${\mathbf{z}}_{\mathrm{x}}$,${\mathbf{z}}_{\mathrm{y}}$ | signals with data gaps from ${\mathbf{x}}$ and ${\mathbf{y}}$ resp.
${\mathbf{w}}_{\mathrm{x}}$,${\mathbf{w}}_{\mathrm{y}}$ | weights indicating valid values in ${\mathbf{z}}_{\mathrm{x}}$ and ${\mathbf{z}}_{\mathrm{x}}$
| resp.
$\mu$ | true mean value
$\mu_{\mathrm{x}}$, $\mu_{\mathrm{y}}$ | true mean values of signals ${\mathbf{x}}$ and ${\mathbf{y}}$ resp.
$\bar{z}$ | estimate of the mean value $\mu$ (obtained from ${\mathbf{z}}$)
$\bar{z}_{\mathrm{x}}$, $\bar{z}_{\mathrm{y}}$ | estimates of the mean values $\mu_{\mathrm{x}}$ and $\mu_{y}$ resp.
| (obtained from ${\mathbf{z}}_{\mathrm{x}}$ and ${\mathbf{z}}_{\mathrm{y}}$)
$\sigma^{2}$ | true signal’s variance
$s^{2}$ | estimate of signal’s variance
$\hat{s}^{2}$ | corrected estimate of signal’s variance
$\sigma_{\bar{z}}^{2}$ | true variance of the mean estimator $\bar{z}$
${\mathbf{\gamma}}$ | true autocovariance function of signal ${\mathbf{x}}$
${\mathbf{\gamma}}_{\mathrm{xy}}$ | true cross-covariance function of signals ${\mathbf{x}}$ and ${\mathbf{y}}$
${\mathbf{C}}$, ${\mathbf{C}}_{\mathrm{xy}}$ | estimates of ${\mathbf{\gamma}}$ and ${\mathbf{\gamma}}_{\mathrm{xy}}$ resp.
$\hat{\mathbf{C}}$, $\hat{\mathbf{C}}_{\mathrm{xy}}$ | corrected estimates of ${\mathbf{\gamma}}$ and ${\mathbf{\gamma}}_{\mathrm{xy}}$ resp.
$\tau_{k}$ |  $k$-th lag time
$K$ | number of lag times considered
$K_{1}$, $K_{2}$ | interval limits of lag times considered
${\mathbf{A}}$, ${\mathbf{A}}_{\mathrm{xy}}$ | mapping matrix (for the covariance)
$a_{kj}$, $a_{\mathrm{xy},kj}$ | elements of the mapping matrix
$\alpha^{\prime}$ | probability that a sample is valid
${\mathbf{\varepsilon}}$, ${\mathbf{\varepsilon}}_{\mathrm{xy}}$ | systematic error (bias) of the autocovariance and
| that of the cross-covariance resp.
${\mathbf{S}}$, ${\mathbf{S}}_{\mathrm{xy}}$ | empirical estimates of the power spectral density
| and the cross-power spectral density resp.
$f_{j}$ |  $j$-th frequency
$D$, ${\mathbf{Z}}$, ${\mathbf{W}}\\!$, ${\mathbf{G}}$, ${\mathbf{H}}$ | auxiliary values and sequences

Assuming a time-limited and equidistantly sampled signal ${\mathbf{x}}$
yielding a set of $N$ values $x_{i}$ with $i=0\ldots N-1$. Sampling instances
are the times $t_{i}=i\Delta\tau$ with the fundamental sampling interval
$\Delta\tau$. The observation ${\mathbf{z}}$ of this signal literally having
also $N$ values $z_{i}$ is assumed to include missing values. The
corresponding weights $w_{i}\in\\{0;1\\}$ provide information on the validity
of the observed values $z_{i}$. For $w_{i}=1$ the value $z_{i}$ is correct,
where $z_{i}=x_{i}$. For $w_{i}=0$ the value $z_{i}$ is faulty, where $z_{i}$
is arbitrary. Further errors within the measurement process are not
investigated.

###  2.1 Mean value and variance

The mean $\mu$ of the true signal ${\mathbf{x}}$ is defined as

$\mu=\left\langle{\mathbf{x}}\right\rangle$ (1)

with the expectation $\langle\cdot\rangle$. Considering appropriate weights
$w_{i}$, the mean value of the measured data set ${\mathbf{z}}$ with missing
data can be obtained based on the valid samples as

$\bar{z}=\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}\right)$ (2)

with

$D=\sum\limits_{i=0}^{N-1}w_{i}.$ (3)

If the appropriate weights $w_{i}$ correspond to the validity of the samples
only and do not include other information about the data, then this mean
estimator is bias-free.

$\langle\bar{z}\rangle=\mu$ (4)

The variance of the original signal ${\mathbf{x}}$ is defined as

$\sigma^{2}=\left\langle\left({\mathbf{x}}-\mu\right)^{2}\right\rangle$ (5)

with the true mean value $\mu$ of the signal. A suitable estimator for the
variance using the measured and gappy data set ${\mathbf{z}}$ considering the
weights $w_{i}$ is

$s^{2}=\frac{1}{D}\left[\sum\limits_{i=0}^{N-1}w_{i}\left(z_{i}-\bar{z}\right)^{2}\right].$
(6)

Here, the true mean value $\mu$ has been replaced by the estimate $\bar{z}$.
Even if the mean estimate $\bar{z}$ is bias-free, the variance estimate
$s^{2}$ is asymptotically bias-free only. The bias decreases with increasing
number of (valid) observations. An appropriate correction for short data sets
is given in Sec. 3.

###  2.2 Autocovariance function and spectrum

The definition of the autocovariance function ${\mathbf{\gamma}}$ of the true
signal ${\mathbf{x}}$ is

$\gamma_{k}=\left\langle\left(x_{i}-\mu\right)\left(x_{i+k}-\mu\right)\right\rangle$
(7)

with the expectation $\langle\cdot\rangle$ and the true mean value $\mu$ of
signal. The autocovariance function is also equidistantly sampled with the
sampling interval $\Delta\tau$. Considering the weights $w_{i}$, the empirical
autocovariance function ${\mathbf{C}}$ can be obtained directly from the
measured and gappy data set ${\mathbf{z}}$ through

$C_{k}=\frac{Z_{k}}{W_{k}}$ (8)

for all lag times $\tau_{k}=k\Delta\tau$ with

$\displaystyle Z_{k}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}(z_{i}-\bar{z})(z_{i+k}-\bar{z})$
(9) $\displaystyle W_{k}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}$ (10)

with $I_{1}=\max(0,-k)$, $I_{2}=\min(N,N-k)-1$ and with the estimated mean
value $\bar{z}$ as above. To accelerate the computing, the values $Z_{k}$ and
$W_{k}$ can also be obtained as the $k$-th values of the records

$\displaystyle{\mathbf{Z}}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\left\\{\left|\mathrm{FFT}\left\\{\underline{\mathbf{w}}\cdot\left(\underline{\mathbf{z}-\bar{z}}\right)\right\\}\right|^{2}\right\\}$
(11) $\displaystyle{\mathbf{W}}\\!$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\left\\{\left|\mathrm{FFT}\left\\{\underline{\mathbf{w}}\right\\}\right|^{2}\right\\}$
(12)

with the element-wise product $\cdot$, where the fast discrete Fourier
transform (FFT) and its inverse (IFFT) are utilized. The records
$\underline{\mathbf{z}-\bar{z}}$ and $\underline{\mathbf{w}}$ are the mean-
subtracted, measured and gappy signal ${\mathbf{z}}-\bar{z}$ and the sequence
of weights ${\mathbf{w}}$, both expanded by zero padding, concatenating $N$
zeros to each of the two sequences. This way the values $Z_{k}$ and $W_{k}$
are obtained for all $k=-N\ldots N-1$ in one step.

Following Wiener-Khinchin’s theorem, see Khintchine (1934) — the transcription
of the name from Cyrillic to Latin letters is ambiguous — the power spectral
density can be obtained from the autocovariance function through the discrete
Fourier transform. Unfortunately, the estimation of the power spectral density
from the autocovariance function obtained from a single data set has an
unacceptable high estimator variance. A common means to reduce the variance is
a subdivision of the data set into shorter blocks. The average of the power
spectral densities of the data blocks then has a significantly smaller
variance. This method is known as Bartlett’s method, see Bartlett (1948,
1950). A disadvantage of Bartlett’s method is that correlations between the
samples at the end of one block and the beginning of the next are not counted.
Furthermore, the wrap-around error may be increased if the assumption is made
that the signal respectively the block is periodic. For too short blocks this
may lead to significant deviations. In contrast, for longer blocks the
reduction of the estimator variance becomes less effective. With Welch’s
method, see Welch (1967), where the statistical functions from overlapping
blocks get averaged, correlations between block boundaries are counted.
However, this also partially generates redundancy. This fact is taken into
account by applying windowing functions to the data blocks prior to their
statistical analysis. Unfortunately, the additional modulation of the data
introduces an additional bias to the estimates of the statistical functions.
However, this is computationally costly.

A powerful alternative without the necessity of block subdivision is the
reduction of the spectral resolution in a post-processing step. Reducing the
spectral resolution corresponds to a shorter support of the covariance
function. For a random process with arbitrarily long but finite memory the
autocovariance function is zero at longer lag times. In this case the
autocovariance function can be shortened to the extent of the longest lasting
correlation without losing information. Assuming that the empirical
autocovariance function ${\mathbf{C}}$, computed for $K$ lag times with $K\ll
2N$ (The factor 2 comes from previous zero padding of the signal.) is
negligible outside the interval $k=-\lfloor
K/2\rfloor\ldots\lfloor(K-1)/2\rfloor$ with $\lfloor i\rfloor$ being the
largest integer smaller or equal to $i$. Then the corresponding empirical
power spectral density ${\mathbf{S}}$ can be obtained from the shorter
autocovariance function ${\mathbf{C}}$ via the discrete (fast) Fourier
transform (FFT)

${\mathbf{S}}=\Delta\tau\cdot\mathrm{FFT}\\{{\mathbf{C}}\\}$ (13)

with

$S_{j}=\Delta\tau\cdot\sum\limits_{k=-\lfloor
K/2\rfloor}^{\lfloor(K-1)/2\rfloor}C_{k}\exp(-2\pi\mathbf{i}f_{j}\tau_{k})$
(14)

for the frequencies $f_{j}=j/(K\Delta\tau)$ with $j=-\lfloor
K/2\rfloor\ldots\lfloor(K-1)/2\rfloor$ and with the imaginary unit
$\mathbf{i}$. Due to the shorter support of the autocovariance function, the
spectral resolution is reduced accordingly, leading also to a significantly
lower estimator variance of the spectrum without introducing new errors from
too short a block subdivision or any modulation of the signal’s amplitude. The
advantage of this method compared to usual block subdivision is that the
correlations of all samples up to a certain maximum lag time are considered
from the entire data set or, at least, from sufficiently long blocks. Further
suppression of the wrap-around error e.g. via the application of a window
function is not needed. This way, one also avoids the spectrum becoming
smeared by additional modulation of the signal. This method is identical to
Blackman and Tukey (1958a, b) used with a rectangular window applied to the
covariance function estimated. Furthermore, this method has been previously
investigated by and briefly mentioned in Bartlett (1948), yielding results
comparable to block averaging. Since longer data records need to be Fourier
transformed, truncation of primarily long covariance estimates is
computationally more expensive. However, it is superior in efficiently using
the information available.

Since the discrete Fourier transform is linear, the above transform ensures
that a bias-free estimate of the covariance function also leads to a bias-free
estimate of the spectrum, except for possible aliasing errors. Note, that this
is valid only under the condition that all correlations of the investigated
process are within the support of the estimated covariance function. For
longer lasting correlations and too short support of the covariance estimate,
the individual values of the estimate still may be bias-free. However, the
spectrum gets smeared if correlations outside the investigated support are
missing. Remaining aliasing errors are not further investigated, since they
occur with too a low temporal resolution of the fundamental sampling process.
Instead, the estimation from time series with invalid samples is called bias-
free, if the estimates have an expectation identical to the equidistantly
sampled process without missing data. This includes appropriate reproduction
of possible aliasing errors.

###  2.3 Cross-covariance function and spectrum

For two signals ${\mathbf{x}}$ and ${\mathbf{y}}$, the cross-covariance
function ${\mathbf{\gamma}}_{\mathrm{xy}}$ is defined as

$\gamma_{\mathrm{xy},k}=\left\langle\left(x_{i}-\mu_{\mathrm{x}}\right)\left(y_{i+k}-\mu_{\mathrm{y}}\right)\right\rangle$
(15)

with the expectation $\langle\cdot\rangle$ and the true mean values
$\mu_{\mathrm{x}}$ of signal ${\mathbf{x}}$ and $\mu_{\mathrm{y}}$ of signal
${\mathbf{y}}$ respectively. From the two data sets $z_{\mathrm{x},i}$ with
$i=0\ldots N_{\mathrm{x}}-1$ and $z_{\mathrm{y},i}$ with $i=0\ldots
N_{\mathrm{y}}-1$ including their appropriate weights $w_{\mathrm{x},i}$ and
$w_{\mathrm{y},i}$, the empirical cross-covariance function
${\mathbf{C}}_{\mathrm{xy}}$ can be estimated directly as

$C_{\mathrm{xy},k}=\frac{Z_{\mathrm{xy},k}}{W_{\mathrm{xy},k}}$ (16)

for all lag times $\tau_{k}=k\Delta\tau$ with

$\displaystyle Z_{\mathrm{xy},k}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=I_{1}}^{I_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}(z_{\mathrm{x},i}-\bar{z}_{\mathrm{x}})(z_{\mathrm{y},i+k}-\bar{z}_{\mathrm{y}})$
(17) $\displaystyle W_{\mathrm{xy},k}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=I_{1}}^{I_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}$
(18)

with $I_{1}=\max(0,-k)$, $I_{2}=\min(N_{\mathrm{x}},N_{\mathrm{y}}-k)-1$ and
with the estimated mean values $\bar{z}_{\mathrm{x}}$ and
$\bar{z}_{\mathrm{y}}$ obtained from the data sets $z_{\mathrm{x}}$ and
$z_{\mathrm{y}}$ and the appropriate sequences of weights $w_{\mathrm{x}}$ and
$w_{\mathrm{y}}$ using the weighted mean estimator in Eq. (2). Using zero
padding and the (inverse) discrete (fast) Fourier transform (FFT resp. IFFT),
the values $Z_{\mathrm{xy},k}$ and $W_{\mathrm{xy},k}$ can also be obtained as
the $k$-th values of the records

$\displaystyle{\mathbf{Z}}_{\mathrm{xy}}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\left\\{\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{x}}}\cdot\left(\underline{\mathbf{z}_{\mathrm{x}}-\bar{z}_{\mathrm{x}}}\right)\right\\}^{\ast}\cdot\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{y}}}\cdot\left(\underline{\mathbf{z}_{\mathrm{y}}-\bar{z}_{\mathrm{y}}}\right)\right\\}\right\\}$
(19) $\displaystyle{\mathbf{W}\\!}_{\mathrm{xy}}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\left\\{\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{x}}}\right\\}^{\ast}\cdot\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{y}}}\right\\}\right\\}$
(20)

with the element-wise product $\cdot$ and the conjugate complex ∗. The
underlined terms are the sequences after zero padding. These are the sequences
${\mathbf{z}}_{\mathrm{x}}-\bar{z}_{\mathrm{x}}$ and
${\mathbf{w}}_{\mathrm{x}}$ extended by $N_{\mathrm{y}}$ concatenated zeros
and ${\mathbf{z}}_{\mathrm{y}}-\bar{z}_{\mathrm{y}}$ and
${\mathbf{w}}_{\mathrm{y}}$ extended by $N_{\mathrm{x}}$ concatenated zeros,
yielding sequences of identical length.

The empirical cross-covariance function ${\mathbf{C}}_{\mathrm{xy}}$, computed
for $K$ lag times $k=K_{1}\ldots K_{2}$, can finally be transformed into the
corresponding empirical power spectral density ${\mathbf{S}}_{\mathrm{xy}}$ by
the discrete (fast) Fourier transform (FFT)

${\mathbf{S}}_{\mathrm{xy}}=\Delta\tau\cdot\mathrm{FFT}\\{{\mathbf{C}}_{\mathrm{xy}}\\}$
(21)

with

$S_{\mathrm{xy},j}=\Delta\tau\cdot\sum\limits_{k=K_{1}}^{K_{2}}C_{\mathrm{xy},k}\exp(-2\pi\mathbf{i}f_{j}\tau_{k})$
(22)

for the frequencies $f_{j}=j/(K\Delta\tau)$ with $j=-\lfloor
K/2\rfloor\ldots\lfloor(K-1)/2\rfloor$ and with the imaginary unit
$\mathbf{i}$. Requirements according the choice of the investigated interval
of the covariance function are identical to those for the autocovariance case
above, namely no correlations outside the investigated interval.

##  3 Correction for short data sets (Bessel’s correction)

The estimators above are asymptotically bias-free only. For short data
records, the amount of valid data may become insufficient and the bias may
become significant. For the estimators introduced above an appropriate
correction is given below. However, other estimators could also benefit from
this correction in the case of short data records. Furthermore, the following
corrections are applied here to gappy data sets only with weights
$w_{i}\in\\{0;1\\}$. However, the procedures are suitable also for other
choices of weights including non-binary values.

###  3.1 Variance estimate

Let the mean estimator of Eq. (2) have the estimator variance
$\sigma_{\bar{z}}^{2}$. Since the variance of a sum of correlated variables is
the sum of all pair-wise covariances, the variance of the sum in the numerator
of Eq. (2) is $\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}w_{i}w_{j}\gamma_{j-i}$
involving the unknown true autocovariance function
$\gamma_{k}={\mathbf{\gamma}}(\tau_{k})$ with $\tau_{k}=k\Delta\tau$. The
denominator in Eq. (2) applies twice, to each factor of the pair-wise
covariances, yielding finally the variance of the mean estimator

$\sigma_{\bar{z}}^{2}=\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\gamma_{j-i}\right).$
(23)

The expectation of the variance estimation as in Eq. (6) with the estimated
mean subtracted from the data is (for the derivation see App. A)

$\langle
s^{2}\rangle=\sigma^{2}-\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\gamma_{j-i}\right)$
(24)

with the true variance $\sigma^{2}$ of the data and again with the true
autocovariance function ${\mathbf{\gamma}}$. The deviation from the correct
variance is exactly the variance of the mean estimator $\sigma_{\bar{z}}^{2}$.

If the variance of the mean estimation is known beforehand, then a bias-free
estimate of the data variance is222If all samples are independent and all
weights are one, this leads to the bias-free estimator
$\hat{s}^{2}=\frac{1}{N-1}\cdot\sum_{i=0}^{N-1}\left(z_{i}-\bar{z}\right)^{2}$,
where the division by $N-1$ instead of $N$ is widely known as Bessel’s
correction for the variance estimate for independent samples. Similar
corrections can be made to estimates of the autocovariance function or the
cross-covariance function derived from two different data sets. Unfortunately,
this requires considering that the data samples are correlated — why one would
otherwise calculate the covariance function?

$\hat{s}^{2}=s^{2}+\sigma_{\bar{z}}^{2}.$ (25)

For unknown a variance of the mean estimation, first the autocovariance
estimate must be corrected. Therefore, the recipe to obtain a bias-free
estimate of the variance is given at the end of Sec. 3.2, first introducing
the correction for the autocovariance estimate.

###  3.2 Autocovariance function and spectrum

Vogelsang and Yang (2016) derived a procedure of an exact correction for short
data sets for the autocovariance estimation. This procedure has been further
developed in Nobach (2023), for the two-sided autocovariance function,
correcting the autocovariance estimate normalized with the number of pairs of
values instead of a constant, including sample weights and also extended to
the cross-covariance case.

The autocovariance estimator in Eq. (8) has the expectation

$\langle C_{k}\rangle=\gamma_{k}+\varepsilon_{k}$ (26)

with the true autocovariance function ${\mathbf{\gamma}}$ and the bias (for
the derivation see App. B)

$\displaystyle\varepsilon_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\gamma_{j-i}\right)$
(27)
$\displaystyle-\frac{1}{DW_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N-1}w_{i}w_{i+k}w_{j}(\gamma_{j-i}+\gamma_{i+k-j})\right]\quad$

with $I_{1}=\max(0,-k)$, $I_{2}=\min(N,N-k)-1$. The value is constant for
uncorrelated data, otherwise it varies with $k$. The first term again is the
variance $\sigma_{\bar{z}}^{2}$ of the mean estimator. Since the true
covariance function ${\mathbf{\gamma}}$ is unknown in real measurements, the
prediction cannot be made directly. However, the relation between the true
covariance function and its estimate is linear. Therefore, one can built a
matrix ${\mathbf{A}}$, mapping a hypothetical covariance function
${\mathbf{\gamma}}$ onto the expectation of the estimated one
$\langle{\mathbf{C}}\rangle$.

$\langle{\mathbf{C}}\rangle={\mathbf{A}}{\mathbf{\gamma}},$ (28)

If the matrix ${\mathbf{A}}$ has the elements $a_{kj}$ then the prediction of
the estimated covariance at lag time $\tau_{k}$ is

$\langle C_{k}\rangle=\sum_{j=K_{1}}^{K_{2}}a_{kj}\gamma_{j}.$ (29)

The range $K_{1}\ldots K_{2}$ of covariances considered should again include
the full range of occurring correlations, such that all true covariances
outside this interval can be neglected. If one sorts the summands in Eq. (26)
and Eq. (27) for an increasing index of the true autocovariance
${\mathbf{\gamma}}$, one obtains the elements of the matrix

$a_{kj}=\delta_{k-j}+\frac{W_{j}}{D^{2}}-\frac{G_{kj}+H_{kj}}{DW_{k}}$ (30)

with

$\delta_{k-j}=\left\\{\begin{array}[]{ll}1&\mbox{for $k-j=0$}\\\
0&\mbox{otherwise}\end{array}\right.$ (31)

and with $W_{j}$ and $W_{k}$ being the $j$-th and the $k$-th values of the
record ${\mathbf{W}}\\!$ above and

$\displaystyle G_{kj}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=\max(0,-j,-k)}^{\min(N,N-j,N-k)-1}w_{i}w_{i+j}w_{i+k}$
(32) $\displaystyle H_{kj}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=\max(0,-j,k-j)}^{\min(N,N-j,N+k-j)-1}w_{i}w_{i+j}w_{i+j-k}.$
(33)

Alternatively, the values $G_{kj}$ and $H_{kj}$ can be obtained as the $j$-th
values of the records

$\displaystyle{\mathbf{G}}_{k}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\Big{\\{}\mathrm{FFT}\left\\{\underline{\mathbf{w}}\cdot\underline{\mathbf{w}}_{\,+k}\right\\}^{\ast}\cdot\mathrm{FFT}\left\\{\underline{\mathbf{w}}\right\\}\Big{\\}}$
(34) $\displaystyle{\mathbf{H}}_{k}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\Big{\\{}\mathrm{FFT}\left\\{\underline{\mathbf{w}}\right\\}^{\ast}\cdot\mathrm{FFT}\left\\{\underline{\mathbf{w}}\cdot\underline{\mathbf{w}}_{\,-k}\right\\}\Big{\\}}$
(35)

with the element-wise product $\cdot$ and the conjugate complex ∗, involving
again the (fast) discrete Fourier transform (FFT) and its inverse (IFFT). The
record $\underline{\mathbf{w}}$ again is the sequence of weights
${\mathbf{w}}$, expanded by zero padding, concatenating $N$ zeros. The records
$\underline{\mathbf{w}}_{\,+k}$ and $\underline{\mathbf{w}}_{\,-k}$ are the
sequence of weights, expanded by zero padding and then shifted by $+k$ or $-k$
time steps respectively. While $\underline{\mathbf{w}}$ at index $i$ has the
weight $w_{i}$, $\underline{\mathbf{w}}_{\,+k}$ has the weight $w_{i+k}$ at
index $i$ and $\underline{\mathbf{w}}_{\,-k}$ has the weight $w_{i-k}$.

The inverse ${\mathbf{A}}^{-1}$ of the matrix ${\mathbf{A}}$ applied to the
primary estimate ${\mathbf{C}}$ yields an improved, bias-free estimate
$\hat{\mathbf{C}}$ of the covariance

$\hat{\mathbf{C}}={\mathbf{A}}^{-1}{\mathbf{C}}.$ (36)

For given $N$ samples $z_{i}$, the covariance function after zero padding has
$2N-1$ non-zero values $C_{k}$ in the range $-(N-1)\ldots N-1$. Unfortunately,
the appropriate matrix ${\mathbf{A}}$ then is singular and the inverse does
not exist for this case. The inverse can be calculated only, if the covariance
function is limited to the range $K_{1}\ldots K_{2}$ with $-(N-1)<K_{1}\leq
K_{2}<N-1$. The improved covariance estimate then is bias-free, as long as the
true covariance of the original signal is zero outside the reduced interval of
lag times $\tau_{K_{1}}\ldots\tau_{K_{2}}$. This coincides with the
requirement that the interval of investigated lag times is larger than the
longest correlation lasts and the observation interval of the signal is longer
than the largest lag time investigated.

The improved estimate $\hat{\mathbf{C}}$ of the covariance function then can
be used to derive the variance of the mean estimator $\sigma_{\bar{z}}^{2}$
following Eq. (23), where the true covariance ${\mathbf{\gamma}}$ is replaced
by the improved, bias-free estimate $\hat{\mathbf{C}}$, and finally to improve
the estimation $\hat{s}^{2}$ of the variance following Eq. (25). Note, that
for the weights being only zero or one, the relations $C_{0}=s^{2}$ and
$\hat{C}_{0}=\hat{s}^{2}$ hold. For other weights, $C_{0}$ may differ from
$s^{2}$ as well as $\hat{C}_{0}$ from $\hat{s}^{2}$.

Finally, the corresponding spectrum is obtained again via the discrete Fourier
transform using Eq. (13).

###  3.3 Cross-covariance function and spectrum

The cross-covariance estimator in Eq. (16) has the expectation

$\langle
C_{\mathrm{xy},k}\rangle=\gamma_{\mathrm{xy},k}+\varepsilon_{\mathrm{xy},k}$
(37)

with the true cross-covariance function ${\mathbf{\gamma}}_{\mathrm{xy}}$ and
the bias (for the derivation see App. C)

$\displaystyle\varepsilon_{\mathrm{xy},k}$ $\displaystyle=$
$\displaystyle\frac{1}{D_{\mathrm{x}}D_{\mathrm{y}}}\left(\sum\limits_{i=0}^{N_{\mathrm{x}}-1}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},j}\gamma_{\mathrm{xy},j-i}\right)$
(38)
$\displaystyle-\frac{1}{D_{\mathrm{y}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{y},j}\gamma_{\mathrm{xy},j-i}\right]$
$\displaystyle-\frac{1}{D_{\mathrm{x}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{x},j}\gamma_{\mathrm{xy},i+k-j}\right]\quad$

with $I_{1}=\max(0,-k)$, $I_{2}=\min(N_{\mathrm{x}},N_{\mathrm{y}}-k)-1$. The
value is constant for uncorrelated data and only if the weights are all
identical for the two data sets, otherwise it varies with $k$. Here

$\displaystyle D_{\mathrm{x}}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},i}$ (39)
$\displaystyle D_{\mathrm{y}}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=0}^{N_{\mathrm{x}}-1}w_{\mathrm{y},i}$ (40)

and $W_{\mathrm{xy},k}$ as above are used for short notation. The matrix
${\mathbf{A}}_{\mathrm{xy}}$, mapping a hypothetical covariance function
${\mathbf{\gamma}}_{\mathrm{xy}}$ onto the expectation of the estimated one
${\mathbf{C}}_{\mathrm{xy}}$ via

$\langle{\mathbf{C}}_{\mathrm{xy}}\rangle={\mathbf{A}}_{\mathrm{xy}}{\mathbf{\gamma}}_{\mathrm{xy}}$
(41)

can be used to predict the estimated covariance at lag time $\tau_{k}$ as

$\langle
C_{\mathrm{xy},k}\rangle=\sum_{j=K_{1}}^{K_{2}}a_{\mathrm{xy},kj}\gamma_{\mathrm{xy},j}$
(42)

with the elements $a_{\mathrm{xy},kj}$ of the matrix
${\mathbf{A}}_{\mathrm{xy}}$. The range $K_{1}\ldots K_{2}$ of covariances
considered should again include the full range of occurring correlations, such
that all true covariance outside this interval can be neglected.

Sorting the summands in Eq. (37) and Eq. (38) for an increasing index of the
true cross-covariance ${\mathbf{\gamma}}_{\mathrm{xy}}$ yields

$a_{\mathrm{xy},kj}=\delta_{k-j}+\frac{W_{\mathrm{xy},j}}{D_{\mathrm{x}}D_{\mathrm{y}}}-\frac{G_{\mathrm{xy},kj}}{D_{\mathrm{y}}W_{\mathrm{xy},k}}-\frac{H_{\mathrm{xy},kj}}{D_{\mathrm{x}}W_{\mathrm{xy},k}}$
(43)

with $\delta_{k-j}$ as above and with $W_{\mathrm{xy},j}$ and
$W_{\mathrm{xy},k}$ being the $j$-th and the $k$-th value of the record
${\mathbf{W}}\\!_{\mathrm{xy}}$ above and

$\displaystyle G_{\mathrm{xy},kj}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=\max(0,-j,-k)}^{\min(N_{\mathrm{x}},N_{\mathrm{y}}-j,N_{\mathrm{y}}-k)-1}w_{\mathrm{x},i}w_{\mathrm{y},i+j}w_{\mathrm{y},i+k}$
(44) $\displaystyle H_{\mathrm{xy},kj}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=\max(0,-j,k-j)}^{\min(N_{\mathrm{x}},N_{\mathrm{y}}-j,N_{\mathrm{x}}+k-j)-1}w_{\mathrm{x},i}w_{\mathrm{y},i+j}w_{\mathrm{x},i+j-k}$
(45)

Alternatively, $G_{\mathrm{xy},kj}$ and $H_{\mathrm{xy},kj}$ can be obtained
as the $j$-th values of the records

$\displaystyle{\mathbf{G}}_{\mathrm{xy},k}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\left\\{\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{x}}}\cdot\underline{\mathbf{w}_{\mathrm{y}}}{}_{\,,+k}\right\\}^{\ast}\cdot\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{y}}}\right\\}\right\\}\qquad$
(46) $\displaystyle{\mathbf{H}}_{\mathrm{xy},k}$ $\displaystyle=$
$\displaystyle\mathrm{IFFT}\left\\{\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{x}}}\right\\}^{\ast}\cdot\mathrm{FFT}\left\\{\underline{\mathbf{w}_{\mathrm{y}}}\cdot\underline{\mathbf{w}_{\mathrm{x}}}{}_{\,,-k}\right\\}\right\\},\qquad$
(47)

involving again the (fast) discrete Fourier transform (FFT) and its inverse
(IFFT) and where again, $\underline{\mathbf{w}_{\mathrm{x}}}$ and
$\underline{\mathbf{w}_{\mathrm{y}}}$ are the weights
${\mathbf{w}}_{\mathrm{x}}$ and ${\mathbf{w}}_{\mathrm{y}}$ after zero padding
and $\underline{\mathbf{w}_{\mathrm{x}}}{}_{\,,-k}$ and
$\underline{\mathbf{w}_{\mathrm{y}}}{}_{\,,+k}$ are additionally shifted by
$-k$ or $+k$ time steps respectively. While
$\underline{\mathbf{w}_{\mathrm{x}}}$ and
$\underline{\mathbf{w}_{\mathrm{y}}}$ at index $i$ have the weights
$w_{\mathrm{x},i}$ and $w_{\mathrm{y},i}$,
$\underline{\mathbf{w}_{\mathrm{x}}}{}_{\,,-k}$ has the weight
$w_{\mathrm{x},i-k}$ at index $i$ and
$\underline{\mathbf{w}_{\mathrm{y}}}{}_{\,,+k}$ has the weight
$w_{\mathrm{y},i+k}$.

The inverse ${\mathbf{A}}_{\mathrm{xy}}^{-1}$ of the matrix
${\mathbf{A}}_{\mathrm{xy}}$ applied to the estimate
${\mathbf{C}}_{\mathrm{xy}}$ yields an improved, bias-free estimate
$\hat{\mathbf{C}}_{\mathrm{xy}}$ of the cross-covariance

$\hat{\mathbf{C}}_{\mathrm{xy}}={\mathbf{A}}_{\mathrm{xy}}^{-1}{\mathbf{C}}_{\mathrm{xy}}.$
(48)

For given $N_{\mathrm{x}}$ samples $z_{\mathrm{x},i}$ and $N_{\mathrm{y}}$
samples $z_{\mathrm{y},i}$, the covariance function after zero padding has
$N_{\mathrm{x}}+N_{\mathrm{y}}-1$ non-zero values $c_{\mathrm{xy},k}$ in the
range $-(N_{\mathrm{x}}-1)\ldots N_{\mathrm{y}}-1$. Unfortunately, the
appropriate matrix ${\mathbf{A}}_{\mathrm{xy}}$ then is singular and the
inverse does not exist for this case. The inverse can be calculated only, if
the covariance function is limited to the range $K_{1}\ldots K_{2}$ with
$-(N_{\mathrm{x}}-1)<K_{1}\leq K_{2}<N_{\mathrm{y}}-1$. The improved
covariance estimate then is bias-free, as long as the true covariance of the
original signal is zero outside the reduced interval of lag times
$\tau_{K_{1}}\ldots\tau_{K_{2}}$. This coincides with the requirement that the
interval of investigated lag times is larger than the longest correlation
lasts and the observation interval of the signal is longer than the largest
lag time investigated.

Finally, the corresponding spectrum is obtained again via the discrete Fourier
transform using Eq. (21).

##  4 Simulation

To demonstrate the ability of the estimation routines to derive consistent
estimates of the covariance functions and the power spectral densities, a
moving-average stochastic process is generated from noise and analyzed in
Monte-Carlo runs. This includes both, bias-free estimation and estimator
variance decreasing with increasing amount of information within the signal.
Both aspects can be measured by means of the root mean square error (RMS),
shown in Sec. 4.2. However, since systematic and random errors get mixed, the
plots of the RMS error alone are not specific enough to identify individual
features of the various estimators. For systematic errors, which are small in
absolute magnitude but significant in comparison to a small true value, or for
small random errors around a biased value, the characteristic dynamics of the
investigated process may get affected significantly, e.g. low-pass filtered.
Characteristic dynamic features of the process under observation may get
concealed this way. Therefore, the next section investigates systematic errors
first, before RMS errors are investigated separately. The empirical mean
estimates are shown in direct comparison to the expected values of the
simulation. This way, the main features of the estimators get more obvious,
especially dynamic errors.

###  4.1 Systematic errors

The simulation uses an autoregressive linear stochastic process of order 100,
yielding an artificial spectrum with an exponentially increasing slope and
with an artificial dip in the observed frequency range. The simulated
autoregressive process has an infinite impulse response and, therefore,
clearly mismatches the assumption of having no correlations outside the
observation interval. This mismatch has been chosen on purpose, to demonstrate
that the introduced estimation methods, especially the correction for short
data sets, are robust against cases, where processes with remaining long
lasting correlations are observed in smaller time windows. The parameters of
the process and the values of the correlation functions and the spectrum are
provided as supplementary data with the present article. Each run of the
simulation process generates two signals with such spectral characteristic
with a total length of $100\,\mathrm{tu}$ (time units). The signals have a
mean value of $8\,\mathrm{au}$ (amplitude units), a variance of
$4\,\mathrm{au^{2}}$, a cross-covariance of $3\,\mathrm{au^{2}}$ and a time
delay between the corresponding signals of $10\,\mathrm{tu}$. Then for each
sample of the two primary signals a weight (1 for a valid sample and 0 for an
invalid one) is chosen from a random process to mimic the data gaps. The first
simulation generates individual samples marked as invalid, independently from
each other, with a probability of $50\,\mathrm{\%}$. The second simulation
generates series of invalid samples, where the state of validity changes with
a probability of $10\,\mathrm{\%}$ at each time step. This procedure also
yields $50\,\mathrm{\%}$ invalid samples on average, where the length of valid
data or that of invalid data has an exponential distribution with a mean of
ten samples and the sequence of weights gets correlated. For invalid data
points, the samples in the data series are set to the constant value of
$-1\,\mathrm{au}$ for visualization purposes and to cross check that the
results get not affected by these faulty values. In Fig. 1 individual
realizations of the signals and the weights are shown. While in Fig. 1a only
individual samples are marked independently as invalid, in Fig. 1b longer
sequences of invalid samples can be seen. However, in both cases
$50\,\mathrm{\%}$ of the samples get marked as invalid on average.

Figure 1: Single realizations of the signals with (a) independent outliers and
(b) with long gaps ($\mathrm{au}$ \- amplitude unit, $\mathrm{tu}$ \- time
unit, $+$ data $z_{\mathrm{x}}$, $\times$ weights $w_{\mathrm{x}}$,
$\bigtriangledown$ data $z_{\mathrm{y}}$, $\bigtriangleup$ weights
$w_{\mathrm{y}}$)

The data sets have been processed to auto- and cross-covariance functions for
the range of lag times between $-25\,\mathrm{tu}$ and $24\,\mathrm{tu}$
($K=50$) respectively $-20\,\mathrm{tu}$ and $29\,\mathrm{tu}$ ($K_{1}=-20$,
$K_{2}=29$). The corresponding power spectral densities have been obtained in
the range between $-0.5\,\mathrm{tu^{-1}}$ and $0.48\,\mathrm{tu^{-1}}$.

Figure 2: Empirical mean from averages over 1 000 realizations with
independent outliers and 100 samples each: (a) autocovariance function, (b)
cross-covariance function, (c) power spectral density (autospectrum), (d)
magnitude of the cross power spectral density (cross-spectrum, $\mathrm{au}$
\- amplitude unit, $\mathrm{tu}$ \- time unit, $+$ sample-and-hold
interpolation w/o deconvolution, $\times$ sample-and-hold interpolation w
deconvolution, $\square$ Lomb-Scargle’s method, $\bigtriangledown$ only valid
points w/o correction for short data sets, $\bigtriangleup$ only valid points
w correction for short data sets, — simulation).

Fig. 2 shows the empirical mean values of the estimated autocovariance
function, that of the cross-covariance function, the auto-power spectral
density and the cross-power spectral density obtained from 1 000 generated
data sets for the case of independent outliers. To simplify the identification
of dynamic errors only one half of the spectra are shown in log-log plots. The
introduced estimation routines have been used for the analysis, including the
correction for short data sets. For comparison, the results of Lomb-Scargle’s
method are also shown as well as those of the sample-and-hold interpolation.
For the latter one, additionally, the results of the deconvolution with the
mapping matrix following Plantier et al. (2012) are given. Interpolation
smears the autocovariance function, yielding deviations, if estimated without
the deconvolution (Fig. 2a). More obvious this dynamic error becomes with the
autospectrum in Fig. 2c, where a frequency dependent deviation from the
reference spectrum occurs, which is characteristic for the interpolation
scheme used. Lomb-Scargle’s method has a significant systematic error of the
autocovariance estimate at lag time zero corresponding to an offset of the
autospectrum. Principally, this offset can be estimated and corrected using

$\hat{\mathbf{S}}=\mathbf{S}-\frac{\Delta\tau}{D}\left(\frac{1}{\alpha^{\prime}}-1\right)\sum\limits_{i=0}^{N-1}w_{i}\left(z_{i}-\bar{z}\right)^{2}$

where $\mathbf{S}$ is the primary, biased estimate of Lomb-Scargle’s method,
$\alpha^{\prime}$ is the probability of a sample to be valid and
$\hat{\mathbf{S}}$ is the improved estimate of the power spectral density. For
details, see Damaschke et al. (2021). However, this offset correction is
suited only for randomly occurring and independent data points missing.
Interpolation with deconvolution and the estimators using only valid samples
show no such significant deviations from the expected spectrum (Fig. 2c).
However, for all estimators without the correction for short data sets, there
remains a deviation of the autocovariance estimates (Fig. 2a) if the estimated
mean value is removed from the data. The additional correction for short data
sets is able to correct this error. However, it is available only for the
estimators using only valid samples. With an equivalent correction, also
sample-and-hold interpolation with deconvolution as well as Lomb-Scargle’s
method could obtain improved estimates of the autocovariance and the spectrum.
However, this holds only for strictly independent outliers, since the
deconvolution procedure implies this assumption and also Lomb-Scargle’s method
has a significant dynamic error for correlated data gaps.

For the cross-covariance and the cross-spectral density the deconvolution
after interpolation has not been derived by Plantier et al. (2012). Without
this deconvolution, the dynamic error will always affect the results for
interpolation leading to similar conclusions as for the autocovariance and
spectrum. No such deviation from the expected values can be identified for the
new estimation methods except for the remaining deviation due to the short
record length if no appropriate correction is used. The additional correction
for short data sets is able to correct this error also for the cross-
covariance case.

The bias due to the short data sets mainly consists of an offset. The linear
term in the estimate of the cross-covariance becomes not obvious in the
present simulation. However, in other spectral compositions it becomes more
obvious and may even lead to an over-estimation of the cross-covariance at
certain lag times. The bias of the autocovariance function is symmetric and,
therefore, has no linear term. Higher-order contributions exist. However, they
are significantly smaller.

Due to the asymmetry of the cross-covariance function the cross-spectral
density becomes complex. In Fig. 2d only the magnitude of the mean complex
power spectral density is shown. In the log-log plot of the spectra (both,
auto and cross) only one half of the spectrum is shown and also the value at
frequency zero is not shown. Since the bias in the covariance functions mainly
consists of an offset, this bias does not contribute significantly to the
spectrum at frequencies different than zero. Therefore, changes due to the
bias-correction for short data sets contribute only little to the values of
the spectra shown in the diagrams.

The correction will be complete only under the assumption, that no correlation
exists outside the observed interval. Since in the present simulation, this is
not entirely fulfilled, the correction is not exact. However, the deviation
from this assumption is small enough and therefore the correction is
sufficient with no obvious bias remaining. Furthermore, this evidences that
the correction is robust against (small) remaining correlations outside the
observed correlation interval.

Figure 3: Empirical mean from averages over 1 000 realizations with long data
gaps and 100 samples each: (a) autocovariance function, (b) cross-covariance
function, (c) power spectral density (autospectrum), (d) magnitude of the
cross power spectral density (cross-spectrum, $\mathrm{au}$ \- amplitude unit,
$\mathrm{tu}$ \- time unit, $+$ sample-and-hold interpolation w/o
deconvolution, $\times$ sample-and-hold interpolation w deconvolution,
$\square$ Lomb-Scargle’s method, $\bigtriangledown$ only valid points w/o
correction for short data sets, $\bigtriangleup$ only valid points w
correction for short data sets, — simulation).

Fig. 3 shows the mean values of the estimated functions for the simulation
with long data gaps, where the occurrence of invalid samples is correlated. In
this case, the model for the deconvolution after sample-and-hold interpolation
is not suited. Therefore, both, the autocovariance function as well as the
autospectral density (Fig. 3a and c), show significant deviations for the
interpolation method even with deconvolution. The deconvolution is not able to
compensate the dynamic error due to the violation of the assumptions for the
derivation of the deconvolution matrix. Without the deconvolution,
interpolation always has such a dynamic error. Lomb-Scargle’s method also has
such a dynamic error, as can be noticed in the autocorrelation function as an
additional smearing of the peak around zero lag time. Correspondingly, the
error in the spectrum is not simply an offset any more and the above
correction of the spectral estimate will fail here. On the contrary, the
methods based on the average over valid samples only, show no such dynamic
errors and the correction for short data sets finally yields bias-free
estimates. This holds for both, the autocovariance and spectrum estimation as
well as for the cross-covariance and spectrum.

###  4.2 Root mean square error

With increasing record length, more and more pairs of valid samples of a
certain lag time get available. The averages converge accordingly and the
estimator variance of the covariance function decreases for each and for all
lag times. As long as the support of the covariance function considered for
the transform into a spectrum remains unchanged, also the estimator variance
of the spectrum decreases accordingly. This also holds for bias correction for
short data sets, since the correction consists of a linear transformation of
the primary estimates. Finally, one can conclude that the introduced method
for bias-free estimation is consistent. This holds independent of the spectral
characteristics of the data gaps. Therefore, it is not necessary to
investigate random outliers and longer data gaps in comparison. Since both
cases are representative, from here on only the results for the more rich case
of data sets with long data gaps are presented. Furthermore, this is the case,
where other methods potentially fail.

Figure 4: Empirical root mean square (RMS) error from averages over 1 000
realizations with long data gaps and 100 samples each (valid and invalid
samples together): (a) autocovariance function, (b) cross-covariance function,
(c) power spectral density (autospectrum), (d) magnitude of the cross power
spectral density (cross-spectrum, $\mathrm{au}$ \- amplitude unit,
$\mathrm{tu}$ \- time unit, $+$ sample-and-hold interpolation w/o
deconvolution, $\times$ sample-and-hold interpolation w deconvolution,
$\square$ Lomb-Scargle’s method, $\bigtriangledown$ only valid points w/o
correction for short data sets, $\bigtriangleup$ only valid points w
correction for short data sets).

Fig. 4 shows the empirical RMS errors obtained from the previous simulation
with long data gaps. Since the new estimator is bias-free, its RMS error can
be smaller than that of the biased estimators where ever the latter ones have
a strong enough bias. In comparison to Lomb-Scarge’s estimator this is obvious
for the autocovariance estimate in Fig. 4a for small lag times, where strong
correlations exist and a strong bias has been found for Lomb-Scarge’s
estimator (cf. Figs. 2a and 3a). For all other lag times, the estimators all
have about the same RMS errors with minor differences. However, sample-and-
hold interpolation without deconvolution owes its small RMS error only to the
strong damping of this method, corresponding to its significant dynamic error
(cf. Figs. 2c and 3c). Coincidentally, this yields values tending towards the
correct ones since correlations are required to vanish at larger lag times for
the correction procedure for short data sets. The estimators with valid data
points only have an almost constant RMS error for all lag times. Including
correction for short data sets increases the RMS error slightly. This also
holds for the cross-correlation and cross-spectrum (Figs. 4b and d). Due to
damping, the biased estimates of the interpolation finally yield lower RMS
errors than those of the bias-free estimators with valid samples only.

Figure 5: Empirical root mean square (RMS) error from averages over 1 000
realizations with long data gaps and 10 000 samples each (valid and invalid
samples together): (a) autocovariance function, (b) cross-covariance function,
(c) power spectral density (autospectrum), (d) magnitude of the cross power
spectral density (cross-spectrum, $\mathrm{au}$ \- amplitude unit,
$\mathrm{tu}$ \- time unit, $+$ sample-and-hold interpolation w/o
deconvolution, $\times$ sample-and-hold interpolation w deconvolution,
$\square$ Lomb-Scargle’s method, $\bigtriangledown$ only valid points w/o
correction for short data sets, $\bigtriangleup$ only valid points w
correction for short data sets).

The advantage of bias-free estimation becomes obvious only for longer data
records, especially for correlated data gaps, as shown in Fig. 5 for data sets
with 10 000 samples each. All other simulation parameters have held unchanged.
While bias-free estimators with increasing amount of data converge towards the
correct values yielding smaller and finally vanishing RMS errors, biased
estimators cannot fall below a certain value. Therefore, the RMS errors of the
new estimators easily fall below the RMS errors of the other (biased)
estimators. This holds at least at those lag times, where these estimators
yield biased estimates, while Lomb-Scargle’s method still has low RMS errors
at large lag times. However, the RMS error of the spectrum (Fig. 5c) uniquely
demonstrates the advantage of the new methods using the valid samples only.
The absence of any bias for these methods finally also leads to smaller RMS
errors compared to biased estimators, at least for an increasing amount of
data.

##  5 Conclusion

Routines for non-parametric estimation of the mean value, the variance, the
auto- and cross-covariance function as well as the corresponding power
spectral densities from signals with invalid samples have been introduced. The
information about the validity of the samples is assumed a priori known. The
introduced estimation methods are based on ensemble averages over valid
samples only, allowing bias-free estimates. Both sources of systematic errors
are tackled, data gaps as well as too short data sets. The routines allow
bias-free estimation, independent of the dynamic characteristics of the gaps.
Furthermore, the estimators are consistent, yielding vanishing RMS errors for
increasing amounts of data. This is the advantage compared to other, biased
estimates, which on the contrary have remaining RMS errors even for infinite
an amount of data.

For independent outliers, Lomb-Scargle’s method yields a significant offset of
the spectral estimate. For correlated data gaps, a dynamic error adds.
Interpolation without deconvolution shows such a dynamic error independent of
the spectral composition of the sampling scheme. Dynamic errors can miss-lead
the interpretation of the process under investigation. Appropriate
deconvolution can invert and correct this dynamic error. However, detailed
information is required about the sampling scheme in correspondence with the
interpolation. So far, appropriate deconvolution matrices are available only
for simple case like purely random outliers with sample-and-hold
interpolation. For any other combination (different interpolation or
correlated sampling schemes), the available deconvolution procedures are not
suited and need significant adaptation.

Bias-free estimation, especially a posteriori bias corrections, typically come
with an increase of estimation uncertainties. Increasing the amount of data
can efficiently reduce estimation uncertainties. Appropriate RMS errors vanish
due to the absence of systematic errors. Therefore, bias-free estimation
outplays its advantages tending for longer data sets. For those, the bias-
correction from Sec. 3 can be omitted, because the simple averages from Sec. 2
are asymptotically bias-free and for long enough data sets they are
sufficient. However, bias correction for short data sets can be applied by
default also to longer data sets to ensure bias-free estimates independent
from the amount of data. There is no risk of an overcorrection or making the
bias worse. In the present study, the procedures for bias-free estimation have
been applied to gappy data sets with weights of one for valid samples and zero
for invalid samples. However, the procedures including the correction for
short data sets are suitable also for other weights, including non-binary
values. The programs used here are available with the present article as
supplementary material.

##  Appendix A Derivation of Eqs. (24) and (25)

From Eq. (6) follows

$s^{2}=\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}^{2}\right)-\frac{2}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}\bar{z}\right)+\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}\bar{z}^{2}\right)$
(49) $\displaystyle=$
$\displaystyle\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}^{2}\right)-\frac{2}{D}\left\\{\sum\limits_{i=0}^{N-1}w_{i}z_{i}\left[\frac{1}{D}\left(\sum\limits_{j=0}^{N-1}w_{j}z_{j}\right)\right]\right\\}$
(50)
$\displaystyle+\frac{1}{D}\left\\{\sum\limits_{i=0}^{N-1}w_{i}\left[\frac{1}{D}\left(\sum\limits_{j=0}^{N-1}w_{j}z_{j}\right)\right]^{2}\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}^{2}\right)-\frac{2}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}z_{i}z_{j}\right)+\left[\frac{1}{D}\left(\sum\limits_{j=0}^{N-1}w_{j}z_{j}\right)\right]^{2}$
(51) $\displaystyle=$
$\displaystyle\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}^{2}\right)-\frac{2}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}z_{i}z_{j}\right)+\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}z_{i}z_{j}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{D}\left(\sum\limits_{i=0}^{N-1}w_{i}z_{i}^{2}\right)-\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}z_{i}z_{j}\right)$
(53) $\displaystyle=$
$\displaystyle\frac{1}{D^{2}}\left[\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\left(z_{i}^{2}-z_{i}z_{j}\right)\right].$
(54)

The expectation of $s^{2}$ then is

$\langle
s^{2}\rangle=\frac{1}{D^{2}}\left\\{\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\left[\left(\sigma^{2}+\mu^{2}\right)-\left(\gamma_{j-i}+\mu^{2}\right)\right]\right\\}$
(55)
$=\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\sigma^{2}\right)-\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\gamma_{j-i}\right)$
(56) $=\sigma^{2}-\sigma_{\bar{z}}^{2}$ (57)

Finally, $s^{2}+\sigma_{\bar{z}}^{2}$ is a bias-free estimate of the variance
$\sigma^{2}$.

##  Appendix B Derivation of Eqs. (26) and (27)

From Eq. (8) follows

$\displaystyle C_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{W_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}z_{i}z_{i+k}\right]-\frac{\bar{z}}{W_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}\left(z_{i}+z_{i+k}\right)\right]+\bar{z}^{2}\qquad$
(58) $\displaystyle=$
$\displaystyle\frac{1}{W_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}z_{i}z_{i+k}\right]-\frac{1}{DW_{k}}\left(\sum\limits_{j=0}^{N-1}w_{j}z_{j}\right)\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{i}w_{i+k}\left(z_{i}+z_{i+k}\right)\right]$
(59)
$\displaystyle+\left[\frac{1}{D}\left(\sum\limits_{j=0}^{N-1}w_{j}z_{j}\right)\right]^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{W_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}z_{i}z_{i+k}\right]-\frac{1}{DW_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N-1}w_{i}w_{i+k}w_{j}\left(z_{i}+z_{i+k}\right)z_{j}\right]$
(60)
$\displaystyle+\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}z_{i}z_{j}\right)$

with $I_{1}=\max(0,-k)$ and $I_{2}=\min(N,N-k)-1$.

The expectation of $C_{k}$ then is

$\displaystyle\langle C_{k}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{W_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}w_{i}w_{i+k}\left(\gamma_{k}+\mu^{2}\right)\right]$
(61)
$\displaystyle-\frac{1}{DW_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N-1}w_{i}w_{i+k}w_{j}\left(\gamma_{j-i}+\gamma_{j-(i+k)}+2\mu^{2}\right)\right]$
$\displaystyle+\frac{1}{D^{2}}\left[\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\left(\gamma_{j-i}+\mu^{2}\right)\right]$
$\displaystyle=$
$\displaystyle\gamma_{k}+\frac{1}{D^{2}}\left(\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}w_{i}w_{j}\gamma_{j-i}\right)$
(62)
$\displaystyle-\frac{1}{DW_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N-1}w_{i}w_{i+k}w_{j}\left(\gamma_{j-i}+\gamma_{i+k-j}\right)\right]$
$\displaystyle=$
$\displaystyle\gamma_{k}+\sigma_{\bar{z}}^{2}-\frac{1}{DW_{k}}\left[\sum\limits_{i=I_{1}}^{I_{2}}\sum\limits_{j=0}^{N-1}w_{i}w_{i+k}w_{j}\left(\gamma_{j-i}+\gamma_{i+k-j}\right)\right].\quad$
(63)

##  Appendix C Derivation of Eqs. (37) and (38)

From Eq. (16) follows

$\displaystyle C_{\mathrm{xy},k}$ $\displaystyle=$
$\displaystyle\frac{1}{W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{x},i}z_{\mathrm{y},i+k}\right]-\frac{\bar{z}_{\mathrm{y}}}{W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{x},i}\right]$
(64)
$\displaystyle-\frac{\bar{z}_{\mathrm{x}}}{W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{y},i+k}\right]+\bar{z}_{\mathrm{x}}\bar{z}_{\mathrm{y}}$
$\displaystyle=$
$\displaystyle\frac{1}{W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{x},i}z_{\mathrm{y},i+k}\right]$
(65)
$\displaystyle-\frac{1}{D_{\mathrm{y}}W_{\mathrm{xy},k}}\left(\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{y},j}z_{\mathrm{y},j}\right)\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{x},i}\right]$
$\displaystyle-\frac{1}{D_{\mathrm{x}}W_{\mathrm{xy},k}}\left(\sum\limits_{j=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},j}z_{\mathrm{x},j}\right)\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{y},i+k}\right]$
$\displaystyle+\frac{1}{D_{\mathrm{x}}D_{\mathrm{y}}}\left(\sum\limits_{j=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},j}z_{\mathrm{x},j}\right)\left(\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{y},j}z_{\mathrm{y},j}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}z_{\mathrm{x},i}z_{\mathrm{y},i+k}\right]$
(66)
$\displaystyle-\frac{1}{D_{\mathrm{y}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{y},j}z_{\mathrm{x},i}z_{\mathrm{y},j}\right]$
$\displaystyle-\frac{1}{D_{\mathrm{x}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}\sum\limits_{j=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{x},j}z_{\mathrm{y},i+k}z_{\mathrm{x},j}\right]$
$\displaystyle+\frac{1}{D_{\mathrm{x}}D_{\mathrm{y}}}\left(\sum\limits_{i=0}^{N_{\mathrm{x}}-1}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},j}z_{\mathrm{x},i}z_{\mathrm{y},j}\right)$

with $I_{1}=\max(0,-k)$ and $I_{2}=\min(N_{\mathrm{x}},N_{\mathrm{y}}-k)-1$.

The expectation of $C_{\mathrm{xy},k}$ then is

$\displaystyle\langle C_{\mathrm{xy},k}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}w_{\mathrm{x},i}w_{\mathrm{y},i+k}\left(\gamma_{\mathrm{xy},k}+\mu_{\mathrm{x}}\mu_{\mathrm{y}}\right)\right]$
(67)
$\displaystyle-\frac{1}{D_{\mathrm{y}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{y},j}\left(\gamma_{\mathrm{xy},j-i}+\mu_{\mathrm{x}}\mu_{\mathrm{y}}\right)\right]$
$\displaystyle-\frac{1}{D_{\mathrm{x}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}\sum\limits_{j=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{x},j}\left(\gamma_{\mathrm{xy},i+k-j}+\mu_{\mathrm{x}}\mu_{\mathrm{y}}\right)\right]$
$\displaystyle+\frac{1}{D_{\mathrm{x}}D_{\mathrm{y}}}\left[\sum\limits_{i=0}^{N_{\mathrm{x}}-1}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},j}\left(\gamma_{\mathrm{xy},j-i}+\mu_{\mathrm{x}}\mu_{\mathrm{y}}\right)\right]$
$\displaystyle=$
$\displaystyle\gamma_{\mathrm{xy},k}+\frac{1}{D_{\mathrm{x}}D_{\mathrm{y}}}\left(\sum\limits_{i=0}^{N_{\mathrm{x}}-1}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},j}\gamma_{\mathrm{xy},j-i}\right)$
(68)
$\displaystyle-\frac{1}{D_{\mathrm{y}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}\sum\limits_{j=0}^{N_{\mathrm{y}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{y},j}\gamma_{\mathrm{xy},j-i}\right]$
$\displaystyle-\frac{1}{D_{\mathrm{x}}W_{\mathrm{xy},k}}\left[\sum\limits_{i=N_{1}}^{N_{2}}\sum\limits_{j=0}^{N_{\mathrm{x}}-1}w_{\mathrm{x},i}w_{\mathrm{y},i+k}w_{\mathrm{x},j}\gamma_{\mathrm{xy},i+k-j}\right].\qquad$

## References

  * (1)
  * Babu and Stoica (2010) Babu, P. and Stoica, P. (2010), ‘Spectral analysis of nonuniformly sampled data — a review’, Digital Signal Processing 20, 359–378.
  * Bartlett (1948) Bartlett, M. S. (1948), ‘Smoothing periodograms from time-series with continuous spectra’, Nature 161, 686–687.
  * Bartlett (1950) Bartlett, M. S. (1950), ‘Periodogram analysis and continuous spectra’, Biometrika 37(1–2), 1–16.
  * Blackman and Tukey (1958a) Blackman, R. D. and Tukey, J. W. (1958a), ‘The measurement of power spectra from the point of view of communications engineering — part I’, The Bell System Technical Journal 37(1), 185–282.
  * Blackman and Tukey (1958b) Blackman, R. D. and Tukey, J. W. (1958b), ‘The measurement of power spectra from the point of view of communications engineering — part II’, The Bell System Technical Journal 37(2), 485–569.
  * Bloomfield (1970) Bloomfield, P. (1970), ‘Spectral analysis with randomly missing observations’, Journal of the Royal Statistical Society. Series B: Statistical Methodology 32(3), 369–380.
  * Damaschke et al. (2018) Damaschke, N., Kühn, V. and Nobach, H. (2018), ‘A fair review of non-parametric bias-free autocorrelation and spectral methods for randomly sampled data in laser Doppler velocimetry’, Digital Signal Processing 76, 22–33.
  * Damaschke et al. (2021) Damaschke, N., Kühn, V. and Nobach, H. (2021), ‘Bias correction for direct spectral estimation from irregularly sampled data including sampling schemes with correlation’, EURASIP J. Adv. Signal Process 2021, 7.
  * Dunsmuir and Robinson (1981a) Dunsmuir, W. and Robinson, P. M. (1981a), ‘Estimation of time series models in the presence of missing data’, Journal of the American Statistica Association 76(375), 560–568.
  * Dunsmuir and Robinson (1981b) Dunsmuir, W. and Robinson, P. M. (1981b), ‘Parametric estimators for stationary time series with missing observations’, Advances in Applied Probability 13(1), 129–146.
  * Everson and Sirovich (1995) Everson, E. and Sirovich, L. (1995), ‘Karhunen-Loève procedure for gappy data’, Journal of the Optical Society of America A 12, 1657–1664.
  * Ferraz-Mello (1986) Ferraz-Mello, S. (1986), ‘Estimation of periods from unequally spaced observations’, Astronomical Journal 86, 619–624.
  * Foster (1996) Foster, G. (1996), ‘Time series analysis by projection. I. statistical properties of Fourier analysis’, Astronomical Journal 111, 541–554.
  * Ghazal and Elhassanein (2006) Ghazal, M. A. and Elhassanein, A. (2006), ‘Periodogram analysis with missing observations’, Journal of Applied Mathematics and Computing 22(1–2), 209–222.
  * Jones (1962a) Jones, R. H. (1962a), ‘Spectral analysis with regularly missed observations’, The Annals of Mathematical Statistics 33(2), 455–461.
  * Jones (1962b) Jones, R. H. (1962b), ‘Spectral estimates and their distributions, part ii’, Scandinavian Actuarial Journal 1962(3–4), 135–153.
  * Jones (1971) Jones, R. H. (1971), ‘Spectrum estimation with missing observations’, Annals of the Institute of Statistical Mathematics 23(1), 387–398.
  * Jones (1972) Jones, R. H. (1972), ‘Aliasing with unequally spaced observations’, Journal of Applied Meteorology 11(2), 245–254.
  * Jones (1980) Jones, R. H. (1980), ‘Maximum likelihood fitting of arma models to time series with missing observations’, Technometrics 22(3), 389–395.
  * Khintchine (1934) Khintchine, A. (1934), ‘Korrelationstheorie der stationären stochastischen Prozesse’, Mathematische Annalen 109, 604–615.
  * Lomb (1976) Lomb, N. R. (1976), ‘Least-squares frequency analysis of unequally spaced data’, Astrophysics and Space Science 39, 447–462.
  * Masry (1983) Masry, E. (1983), Spectral and probability density estimation from irregularly observed data, in E. Parzen, ed., ‘Statistical Analysis of Irregularly Observed Time Series in Statistics, Lecture Notes in Statistics’, Springer, New York, pp. 224–250.
  * Mathias et al. (2004) Mathias, A., Grond, F., Guardans, R., Seese, D., Canela, M. and Diebner, H. H. (2004), ‘Algorithms for spectral analysis of irregularly sampled time series’, Journal of Statistical Software 11, 1–27.
  * Munteanu et al. (2016) Munteanu, C., Negrea, C., Echim, M. and Mursula, K. (2016), ‘Effect of data gaps: comparison of different spectral analysis methods’, Ann. Geophys. 34(4), 437–449.
  * Niedźwiecki and Cisowski (1993) Niedźwiecki, M. and Cisowski, K. (1993), Adaptive scheme for elimination of background noise and impulsive disturbances from audio signals, pp. 519–522.  Proceedings of the Quatrozieme Colloque GRETSI, Juan-les-Pins, France.
  * Nobach (2019) Nobach, H. (2019), ‘Note on nonparametric spectral analysis of wideband spectrum with missing data via sample-and-hold interpolation and deconvolution’, Digital Signal Processing 87, 19–20.
  * Nobach (2023) Nobach, H. (2023), ‘Practical realization of Bessel’s correction for a bias-free estimation of the auto-covariance and the cross-covariance functions’.  arXiv:2303.11047 [stat.ME].
  * Nobach et al. (1998) Nobach, H., Müller, E. and Tropea, C. (1998), ‘Efficient estimation of power spectral density from laser Doppler anemometer data’, Experiments in Fluids 24, 499–509.
  * Parzen (1963) Parzen, E. (1963), ‘On spectral analysis with missing observations and amplitude modulation’, Sankhya: The Indian Journal of Statistics, Series A 245, 383–392.
  * Plantier et al. (2012) Plantier, G., Moreau, S., Simon, L., Valière, J.-C., Duff, A. L. and Baillet, J. (2012), ‘Nonparametric spectral analysis of wideband spectrum with missing data via sample-and-hold interpolation and deconvolution’, Digital Signal Processing 22, 994–1004.
  * Preihs et al. (2012) Preihs, S., Stöter, F.-R. and Ostermann, J. (2012), Low delay error concealment for audio signals, in ‘46th AES Conference on Audio Forensics’.
  * Rivoira and Fleury (2004) Rivoira, A. and Fleury, G. A. (2004), ‘A consistent nonparametric spectral estimator for randomly sampled signals’, IEEE Transactions on Signal Processing 52, 2383–2395.
  * Robinson (1977) Robinson, P. M. (1977), ‘Estimation of a time series model from unequally spaced data’, Stochastic Processes and their Applications 6(1), 9–24.
  * Robinson (1985) Robinson, P. M. (1985), ‘Testing for serial correlation in regression with missing observations’, Journal of the Royal Statistical Society. Series B: Statistical Methodology 47(3), 429–437.
  * Scargle (1982) Scargle, J. D. (1982), ‘Studies in astronomical time series analysis. II. statistical aspects of spectral analysis of unevenly spaced data’, The Astrophysical Journal 263, 835–853.
  * Scheinok (1965) Scheinok, P. A. (1965), ‘Spectral analysis with randomly missed observations: The binomial case’, The Annals of Mathematical Statistics 36(3), 972–977.
  * Stoica et al. (2000) Stoica, P., Larsson, E. G. and Li, J. (2000), ‘Adaptive filter-bank approach to restoration and spectral analysis of gapped data’, The Astronomical Journal 120(4), 2163–2173.
  * Stoica and Sandgren (2006) Stoica, P. and Sandgren, N. (2006), ‘Spectral analysis of irregularly-sampled data: Paralleling the regularly-sampled data approach’, Digital Signal Processing 16(6), 712–734.
  * Sujbert and Orosz (2016) Sujbert, L. and Orosz, G. (2016), ‘FFT-based spectrum analysis in the case of data loss’, IEEE Trans. Instrumentation and Measurement 65, 968–976.
  * Venturi and Karniadakis (2004) Venturi, D. and Karniadakis, G. E. (2004), ‘Gappy data and reconstruction procedures for flow past a cylinder’, Journal of Fluid Mechanics 519, 315–336.
  * Vio et al. (2000) Vio, R., Strohmer, T. and Wamsteker, W. (2000), ‘On the reconstruction of irregularly sampled time series’, Publications of the Astronomical Society of the Pacific 112(767), 74–90.
  * Vogelsang and Yang (2016) Vogelsang, T. J. and Yang, J. (2016), ‘Exactly/nearly unbiased estimation of autocovariances of a univariate time series with unknown mean’, Journal of Time Series Analysis 37, 723–740.
  * Welch (1967) Welch, P. D. (1967), ‘The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms’, IEEE Transactions on Audio and Electroacoustics 15(2), 70–73.