A discrete time domain approach on time delay estimation


ACTA IMEKO 
ISSN: 2221-870X 
March 2019, Volume 8, Number 1, 77 - 86 

 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 77 

A discrete time domain approach on time delay estimation 

Jorma Kekalainen 

 

 

Section: RESEARCH PAPER  

Keywords: Time delay estimation; criterion function; index delay; cross-correlation; reliability checking 

Citation: Jorma Kekalainen, A discrete time domain approach on time delay estimation, Acta IMEKO, vol. 8, no. 1, article 11, March 2019, identifier: IMEKO-
ACTA-08 (2019)-01-11 

Section Editor: Konrad Jedrzejewski, Warsaw University of Technology, Poland 

Received April 26, 2018; In final form February 19, 2019; Published March 2019 

Copyright: © 2019 IMEKO. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Corresponding author: Jorma Kekalainen, e-mail: jormakekalainen@outlook.com 

1. INTRODUCTION 

Time delay estimation plays a significant role in numerous 
engineering applications, specifically signal processing [1], [2]. 
The problem of estimating time delays between noisy signals has 
attracted considerable attention over the years due to the variety 
of applications in areas such as acoustics, sonar, radar, flow, and 
velocity measurements. The accurate estimation of time delay is 
often complicated by conditions of poor signal-to-noise ratio, a 
changing time delay, multipath propagation, and dispersion. 
Various time delay estimation procedures have been proposed 
and implemented over the decades. The general goal has been to 
minimise the variance of the time delay estimates in the presence 
of uncorrelated noise at the receiving ends.  

For time delay estimation in linear systems, there are three 
basic operations involved: 1) filtration of input signals in the case 
of non-white signal and noise spectral characteristics in order to 
accentuate the frequency bands with good signal-to-noise ratio; 
2) comparison (or correlation) of the two filtered waveforms in 
order to apply the detection by means of threshold; and 3) 
computation of time delay estimate. 

Time delay estimation procedures based on cross correlation, 
smoothed coherence transforms, unit impulse response, and 
maximum likelihood estimates can generally be viewed as the 
inverse Fourier transforms of normalised (weighted) cross-

spectral density functions that yield delay estimates in time 
domain terms. Many of the methods proposed for time delay 
estimation have been shown to be related by means of the 
generalised cross-correlation (GCC) approach, which involves 
filtering the received signals and estimating the time delay as the 
time lag where the cross-correlation function of filtered signals is 
at its maximum [3]. 

The GCC function and the estimate of time delay can be 
expressed in the well-known form 

{
𝑅(𝜏) = F−1{𝑊(𝑓)𝐺12(𝑓)}

𝜏𝐷𝑒𝑙𝑎𝑦 = Arg{Max
𝜏
[𝑅(𝜏)]}

 (1) 

where F-1 is the inverse Fourier transform, W(f) = H1(f)H2(f)* is 
a weighting function determined by the filter frequency 
responses (* is used to denote the complex conjugate), and G12 
is the cross power density spectrum of the received sensor signals 
1 and 2. The Arg symbol denotes an argument operator, that 
picks the time lag or delay argument τ of R. 

Cross-correlation is mainly suggested in applications in which 
some knowledge of time delays is required, such as the velocity 
measurement of solid or non-solid material flow (e.g. steel and 
paper strips; road and rail vehicles; and mixtures of solids, liquids, 
and gas) and range measurement by sonar and radar [4]-[10]. 

ABSTRACT 
A time delay estimation method based on the discrete time domain approach is introduced here. In this dual-channel time delay 
estimation model, the criterion function compares the time differences of time sequences between channels, not the magnitude values 
of time functions as in the conventional cross-correlation method. An estimation task is formulated as an extreme value problem in 
discrete index space. Using the index delay giving extreme value to the criterion function, it is possible to find the best estimate for time 
delay distribution in the meaning of that criterion. Using this method, the estimated delay distribution and criterion function are clearly 
separated. Thus, there are no theoretical problems in the determination of the average time delay or velocity in the non-constant or 
changing time delay case as long as a sufficient statistical similarity (correlation) exists between channel signals.   
The theoretical values of several criterion functions and the probability of occurrence of an anomalous estimate with the cross-
covariance criterion function are derived. A basic performance analysis of the estimation method is presented. Some potential real-time 
supervision methods based on the use of criterion functions in the detection of the possible unreliability of the time delay estimate are 
outlined. 

mailto:jormakekalainen@outlook.com


 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 78 

If the thresholding device in the basic model of Figure 1 is 
moved to the front of the cross-correlator, a simple correlation 
scheme is obtained [11]. Threshold crossing can be used as a 
criterion to save data, which is then processed by correlation and 
delay estimation. This algorithm, which allows for fast polarity 
correlation, was proposed by Henry [12]. 

Henry's idea is based on the fact that information used by a 
polarity cross-correlator comprises the actual zero-crossing times 
of the signals [13]. Henry noted that the zero crossings are 
randomly distributed in time; hence, the observation of these 
zero crossings gives data with a sufficient time resolution without 
the need for sampling the data at close intervals that were 
previously thought necessary to give an adequate resolution. 

If polarity information is replaced by crossing over times, then 
the time delay estimate is the time interval between the 
corresponding crossing over times on different channels. To find 
the crossing over times corresponding to each other in both 
channels, some criterion function is needed for the decision-
making process involved in the detection of the time delay 
estimate. The aim of this paper is to present the basic idea and 
analysis of a method that utilises pre-processed time signals and 
some simple criterion functions in time delay estimation. 

2. TIME DIFFERENCE DISTRIBUTION METHOD 

2.1 Time differences and criterion functions 

To determine an estimate for the time delay between two 
axially separated sensors, time differences are formed using the 
same sensor arrangement as in conventional cross-correlation 
methods, and the utilisation of information included in the zero-
crossing times of the signals or any time moment is actually based 
on some other detectable and randomly distributed 
characteristics (vortices in turbulent flow, gas bubbles, clumps of 
particles, macroscopic objects, surface patterns, injected 
chemical or radioisotope  markers, etc.) in flow. As a result, two 
primary time sequences {t} and {T} are obtained from one 
channel (input) and from the other channel (output). Of course, 
in any case, a sufficient statistical similarity (correlation) between 
channel signals with some time delay is an essential prerequisite 
for the successful utilisation of the method. 
The time difference distribution (TDD) method for solving time 
delay estimation problems is based on the use of any criterion 
function that is able to utilise that statistical similarity between 
the channel signals and to fix the time elements corresponding 
each other in sequences {t} and {T}. Using the index delay that 
gives the extreme value (maximum or minimum value depending 
on the criterion type) to the criterion function, it is possible to 
select one collection of (t, T)-pairs and to calculate the best 

estimate for time delay distribution (j , j=1, ... ) within the 
meaning of that criterion. With the help of estimated time delay 

distribution, the velocity distribution (L/j , j=1, ... ) can be 
calculated when the distance (L) between the axially separated 

sensors is known. Of course, sufficient statistical similarity 
(correlation) between channel signals sets an upper limit for the 
distance L. 
All suitable criterion functions with the previous time difference 
formulation of the time delay estimation problem have been 
inherently defined in the discrete (index) domain. Next, mainly 
as an example, we present the simplest possible criterion 
functions. The theoretical counterparts of these functions are 
presented in section 2.2. 

The first criterion function is called the sample variance delay 
function êvd. As a criterion of the best estimate, we use the 
minimum value of the sample variance delay function 

Min
𝑘
{�̂�𝑣𝑑(𝑘)} = Min

𝑘
{⟨∆2(𝑖 + 𝑘,𝑖)⟩ − ⟨∆(𝑖 + 𝑘,𝑖)⟩2} (2) 

where (i+k,i) = T(i+k)-t(i), T(i+k) is the (i+k)th element of 

{T}, t(i) is the ith element of {t}, k is an index delay, and  
denotes a sample mean. 

The second criterion function is called the sample average 
delay difference function êadd 

Min
𝑘
{�̂�𝑎𝑑𝑑(𝑘)}

= Min
𝑘
{⟨|∆(𝑖 + 𝑘,𝑖) − ∆(𝑖 + 𝑘 − 1,𝑖 − 1)|⟩} 

(3) 

where  is the absolute value. 
The third criterion function is the sample variance delay 

difference function êvdd 

Min
𝑘
{�̂�𝑣𝑑𝑑(𝑘)}

= Min
𝑘
{
⟨|∆(𝑖 + 𝑘,𝑖) − ∆(𝑖 + 𝑘 − 1,𝑖 − 1)|2⟩

−⟨|∆(𝑖 + 𝑘,𝑖) − ∆(𝑖 + 𝑘 − 1,𝑖 − 1)|⟩2
} . 

(4) 

The fourth criterion function is correlation coefficient êcor 

Max
𝑘
{�̂�𝑐𝑜𝑟(𝑘)}

= Max
𝑘
{
⟨∆𝑇(𝑖 + 𝑘)∆𝑡(𝑖)⟩ − ⟨∆𝑇(𝑖 + 𝑘)⟩⟨∆𝑡(𝑖)⟩

√𝑉�̂�𝑟{∆𝑇}𝑉�̂�𝑟{∆𝑡}
} 

(5) 

where 

𝑉â𝑟{∆𝛵} = ⟨∆𝛵2(𝑖 + 𝑘)⟩ − ⟨∆𝛵(𝑖 + 𝑘)⟩2 
𝑉â𝑟{∆𝑡} = ⟨∆𝑡2(𝑖)⟩ − ⟨∆𝑡(𝑖)⟩2 
∆𝑇(𝑖 + 𝑘) =  𝑇(𝑖 + 𝑘) − 𝑇(𝑖 + 𝑘 − 1) 
∆𝑡(𝑖) =  𝑡(𝑖) − 𝑡(𝑖 − 1). 
As the denominator and product of averages T(i+k) 

and t(i) are constants independent of index delay k in 
stationary sequences, the maximisation of the correlation 
coefficient is equivalent to the maximisation of the cross-
correlation function 

Max
𝑘
{�̂�𝑐𝑐(𝑘)} = Max

𝑘
{⟨∆𝑇(𝑖 + 𝑘)∆𝑡(𝑖)⟩} (6) 

where êcc(k) is sample cross-correlation function of time intervals. 
The minimisation of the quadratic delay difference function 

Min
𝑘
{�̂�𝑞𝑑𝑑(𝑘)} = Min

𝑘
{⟨|∆𝑇(𝑖 + 𝑘) − ∆𝑡(𝑖)|2⟩} (7) 

in stationary sequences is also equivalent to the maximisation of 
the correlation coefficient and the cross-correlation function.  

2.2 Theoretical values of some criterion functions  

In the following, stationary FIFO is assumed, meaning here that 
the detectable objects in the flow do not mix in the space 

between the sensors  undistorted delay distribution case. 
In fluid mechanics, a plug flow is a simple model of the velocity 

 

crosscorrelator threshold delay estimator 

H 
1 
(f) 

H (f) 
2 

 

Figure 1. Block diagram model for dual-channel delay estimation with the 
generalised correlation method. 



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 79 

profile of a fluid flowing in a pipe. In plug flow, the velocity of 
the fluid is assumed to be constant across any cross-section of 
the pipe perpendicular to the axis of the pipe. In an ideal plug 
flow, the fluid is not mixed in the axial direction in the space 
between the sensors, and any infinitesimal fluid plug element that 
enters the space between sensors at time t exits the space between 

the sensors at time T=t+, where  is the delay time of any plug 
between the sensors. The delay time distribution function is, 

therefore, the Dirac delta function at . However, any real plug 
flow has a delay time distribution that is a ‘narrow’ pulse around 
the mean delay time. Plug flow is a special type of previously 
assumed FIFO flow.  

The expectation of the sample variance delay function êvd (in 

the stationary FIFO case) for independent random variables t 

and  can be presented in the form 

𝑒𝑣𝑑(𝑘) = 𝐸{�̂�𝑣𝑑} = 𝐸{⟨∆
2(𝑖 + 𝑘,𝑖)⟩ − ⟨∆(𝑖 + 𝑘,𝑖)⟩2}

= |𝑘 − 𝑘0|𝑉𝑎𝑟{∆𝑡} + 𝑉𝑎𝑟{𝜏} 
(8) 

where E is the expectation operator, Var{t} is the variance of 

the time interval in the sequence {t}, Var{} is the variance of 

time delay , and k0 is a fixed index delay that minimises the 
criterion function êvd. 

There is a minimum of evd 

𝑀𝑖𝑛
𝑘=𝑘0

{𝑒𝑣𝑑(𝑘)} = 𝑉𝑎𝑟{𝜏} . (9) 

Therefore, moving (k-k0) index units from the optimum index 
delay k0, the value of the criterion function increases by an 
amount of time interval variance weighted with the absolute 
value of (k-k0) index units, as we see in Figure 2. 

The sample criterion function previously denoted as êadd can 
be presented in the form 

�̂�𝑎𝑑𝑑(𝑘)
= ⟨|∆𝜏(|𝑘 − 𝑘0| + 𝑖) + ∆𝑡(|𝑘 − 𝑘0| + 𝑖)
− ∆𝑡(𝑖)|⟩ ≥ ⟨|∆𝜏(|𝑘 − 𝑘0| + 𝑖)|⟩ 

(10) 

where the time delay difference is (s) = (s)-(s-1). 

Equality in equation (10) is true if time interval t() between the 
sequential elements in sequence {t} is constant or index delay 
k=k0, which minimises the value of criterion function êadd. 

 
The theoretical value of êadd(k) is 

𝑒𝑎𝑑𝑑(𝑘) = ∫ |𝑧|∫ 𝑓∆[∆𝑡](𝑧 − ∆𝜏)𝑓∆𝜏(∆𝜏)𝑑∆𝜏𝑑𝑧
∞

−∞

∞

−∞

 (11) 

where 

𝑓∆[∆𝑡](𝑧 − ∆𝜏) = ∫ 𝑓∆𝑡(∆𝑡)𝑓∆𝑡(𝑧 − ∆𝜏 − ∆𝑡)𝑑∆𝑡
∞

−∞
, 

𝑓∆𝜏(∆𝜏) = ∫ 𝑓𝜏(𝜏)𝑓𝜏(∆𝜏 − 𝜏)𝑑𝜏
∞

−∞
, 

𝑧 = ∆[∆𝑡] + ∆𝜏, 
f∆t is a probability density function of t, and f is a probability 

density function of  . 
The minimum value of eadd is 

Min
𝑘=𝑘0

{𝑒𝑎𝑑𝑑(𝑘)} = ∫ |𝑧|
∞

−∞

∫ 𝑓𝜏(𝜏)𝑓𝜏(𝑧 − 𝜏)𝑑𝜏𝑑𝑧
∞

−∞

 (12) 

where z = . 
Figure 3 shows the coincidence of the theoretical and sample 
values of the criterion function êadd. 

The sample variance delay difference function êvdd can be 
presented in the form 

�̂�𝑣𝑑𝑑(𝑘)
= ⟨|Δ𝜏(|𝑘 − 𝑘0| + 𝑖) + Δ𝑡(|𝑘 − 𝑘0| + 𝑖) − Δ𝑡(𝑖)|

2⟩
− ⟨|Δ𝜏(|𝑘 − 𝑘0| + 𝑖) + Δ𝑡(|𝑘 − 𝑘0| + 𝑖) − Δ𝑡(𝑖)|⟩

2 . 

(13) 

The theoretic counterpart of êvdd with the previous notes is 

𝑒𝑣𝑑𝑑(𝑘)

= [∫ 𝑧2
∞

−∞

∫ 𝑓Δ[Δ𝑡](𝑧 − Δ𝜏)𝑓Δ𝜏(Δ𝜏)𝑑Δ𝜏𝑑𝑧
∞

−∞

]

− [∫ |𝑧|∫ 𝑓Δ[Δ𝑡](𝑧 − Δ𝜏)𝑓Δ𝜏(Δ𝜏)𝑑Δ𝜏𝑑𝑧
∞

−∞

∞

−∞

]

2

. 

(14) 

The minimum value of evdd is 

Min
𝑘=𝑘0

{𝑒𝑣𝑑𝑑(𝑘)}

= [∫ 𝑧2
∞

−∞

∫ 𝑓𝑧(𝜏)𝑓𝜏(𝑧 − 𝜏)𝑑𝜏𝑑𝑧
∞

−∞

]

− [∫ |𝑧|∫ 𝑓𝑧(𝜏)𝑓𝜏(𝑧 − 𝜏)𝑑𝜏𝑑𝑧
∞

−∞

∞

−∞

]

2

. 

(15) 

 
 

4 

6 

8 

10 

12 

14 

6 7 8 9 10 11 12 13 14 15 16 

Index delay 

Values of criterion function êvd 

 

Figure 2. Theoretical and sample values of the square root of the variance 
delay function êvd. All 11 theoretical values presented in this index delay 
window are combined with the line. Sample distributions in the simulation 

experiments are approximately Gaussian tN(30,5²) and N(300,5²). The 
sample size is 1000. Theoretical distributions are exactly Gaussian with the 
previous parameters. The bar denotes the variation range of the sample 
values in 30 experiments. 

 

5 

6 

7 

8 

9 

6 7 8 9 10 11 12 13 14 15 16 

Index delay 

Values of criterion function êadd 

 

Figure 3. Theoretical and sample values of the average delay difference 
function êadd with the same statistics as in Figure 2. Theoretical values 
computed from the equation (11) and equation (12) are combined with the 
line. The bar denotes the variation range of the sample values in 30 
experiments. 



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 80 

 

Figure 4 shows the coincidence of the theoretical and sample 

values of the criterion function.  

The theoretical counterpart of the correlation coefficient can 

be presented in the form of normalised cross-covariance 

function 

𝑒𝑐𝑜𝑟(𝑘) =
𝐶𝑜𝑣{Δ𝑇(𝑖 + 𝑘),Δ𝑡(𝑖)}

√𝑉𝑎𝑟{Δ𝑇(𝑖 + 𝑘)}𝑉𝑎𝑟{Δ𝑡(𝑖)}
 . (16) 

The maximum value of ecor is 

Max
𝑘=𝑘0

{𝑒𝑐𝑜𝑟(𝑘)} = √
𝑉𝑎𝑟{Δ𝑡}

𝑉𝑎𝑟{Δ𝑡} + 2𝑉𝑎𝑟{𝜏}
 (17) 

if t and  are independent. 

Equivalently. cross-correlation function can be presented in 

the form 

𝑒𝑐𝑐(𝑘) = 𝑅Δ𝑡Δ𝑡(𝑘) + 𝑅Δ𝑡Δ𝜏(𝑘) (18) 

where Rtt is the autocorrelation function of t, and Rt is the 

cross-correlation function of t and . 

The maximum value of the cross-correlation criterion 

function for independent t and  is 

Max
𝑘=𝑘0

{𝑒𝑐𝑐(𝑘)} = 𝑉𝑎𝑟{Δ𝑡} + [𝐸{Δ𝑡}]
2 (19) 

because E{t}0 here, it is convenient to work with the cross-

covariance form of the criterion function 

𝑒𝑐𝑜𝑣(𝑘) = 𝑒𝑐𝑐(𝑘) − [𝐸{Δ𝑡}]
2 , (20) 

then 

Max
𝑘=𝑘0

{𝑒𝑐𝑜𝑣(𝑘)} = 𝑉𝑎𝑟{Δ𝑡} . (21) 

The quadratic delay difference function eqdd can be presented 

in the form 

𝑒𝑞𝑑𝑑(𝑘) = 𝐸{⟨|Δ𝑇(𝑖 + 𝑘) − Δ𝑡(𝑖)|
2⟩}

= 2[𝑉𝑎𝑟{Δ𝑡} + 𝑉𝑎𝑟{𝜏}] 
(22) 

and the minimum value of the criterion function is 

Min
𝑘=𝑘0

{𝑒𝑞𝑑𝑑(𝑘)} = 2𝑉𝑎𝑟{𝜏} (23) 

if t and  are independent. 

Because the maximisation of ecor, ecc, and ecov is equivalent to the 

minimisation of eqdd in stationary systems, Figure 5 only presents 

the theoretical and sample values of the criterion function ecor. 

Theoretical values of the previous criterion functions 

experimentally approached by the sample criterion functions 

were derived, assuming stationary stochastic processes 

generating independent time differences and FIFO-type order 

relation. 

3.  PERFORMANCE OF TIME DIFFERENCE DISTRIBUTION 
METHOD 

3.1 Probability of occurrence of an anomalous estimate 

Predicting the performance of time delay estimation methods, 
it is often possible to set a measure of performance by 
establishing bounds on performance. For the time delay 
estimation problem, the Cramer-Rao Lower Bound (CRLB) is 
commonly used as a performance standard. The CRLB yields a 
lower bound on the variance of any unbiased time delay estimate 
as a function of the signal and noise power spectra and the 
coherent processing time [3], [14]-[16]. 

Concerning cases of practical interest, there is a theorem that 
provides that the maximum likelihood estimate can be made 
arbitrary, close to the CRLB, for sufficiently long observation 
times [17]. Often, it is not possible to obtain observation times 
that are sufficiently long. Then, random variations in the values 
of the criterion function can mask correct extreme values, and 
large estimation errors or so-called anomalous estimates may 
occur. Therefore, the actual performance is much worse than the 
CRLB predicted by the linear analysis for a given observation 
time and signal-to-noise ratio [18], [19]. 

Let us estimate the reliability of the time delay estimate of this 
TDD method, with a finite observation time. The maximum type 
of the criterion function probability of the occurrence of an 
anomaly can be described as 

Values of criterion function êvdd 

 

4 

5 

6 

7 

6 7 8 9 10 11 12 13 14 15 16 

Index delay 

Values of criterion function êvdd 

 

 

Figure 4. Theoretical and sample values of the square root of evdd with the 
same statistics as in Figure 2. Theoretical values computed from the 
equations (14) and (15) are combined with line. Bar denotes a variation range 
of sample values in thirty experiments. 

 

-0,1 

0,0 

0,1 

0,2 

0,3 

0,4 

0,5 

0,6 

6 7 8 9 10 11 12 13 14 15 16 

Index delay 

Values of criterion function êcor 

 

 

Figure 5. Theoretical and sample values of the criterion function ecor with the 
same statistics as in Figure 2. Theoretical values computed from equation 
(16) and equation (17) are combined with the line. The bar denotes the 
variation range of the sample values in 30 experiments. 



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 81 

𝑃𝑎 = ∫ 𝑓�̂�(𝑚)(𝑥)
∞

−∞

[
 
 
 
 

1 − ∏ ∫ 𝑓�̂�(𝑗+𝑚)(𝑦)𝑑𝑦
𝑥

−∞

𝜔
2

𝑗=−
𝜔
2
 

(𝑗≠0) ]
 
 
 
 

𝑑𝑥 (24) 

where fê(m)(x) and fê(j+m)(y) denote the density functions of 
independent random variables ê(m) possessing correct index 

delay, and ê(j+m) possessing incorrect index delay if j0. The 
range of criterion functions (delay windows) over which the 

extreme value is searched is +1 ( is an even number). 
The probability of occurrence of an anomaly with the 

maximum criterion function type with identical Gaussian 
densities fê(k)(y) is in the normalised form 

𝑃𝑎 = 1 −
1

√2𝜋
∫ 𝑒

−
𝜃2

2

∞

−∞

{
1

√2𝜋
∫ 𝑒

−
𝑡2

2 𝑑𝑡
𝑧(𝜃)

−∞

}

𝜔

𝑑𝜃 (25) 

and 

𝑧(𝜃) =
𝜃√𝑉𝑎𝑟{�̂�(𝑚)} + 𝐸{�̂�(𝑚)} − 𝐸{�̂�(𝑘)}

√𝑉𝑎𝑟{�̂�(𝑘)}
 (26) 

where k=j+mm. 

Assumption of a Gaussian distribution of ê() is reasonable 
according to the central limit theorem in view of the large 
number of random variables summed to produce criterion 

function ê(). Thus, the quantities needed for densities of ê() are 

their means E{ê()} and variances Var{ê()}. Given the means 
and variances of the sample criterion function, the probability of 
occurrence of an anomaly can be numerically evaluated. 

As a special example of the maximum criterion function type 
is the cross-correlation function (or, in fact, the cross-covariance 
function) because it turns out that whenever the zero mean 

assumption is unrealistic (here E{t} = t0), it is more 
convenient to work with covariance functions than with 
correlation functions. Assuming a Gaussian distribution, and 

since Côvxy() is the unbiased estimate of Covxy() or 

E{Côvxy()}=Covxy(), the only other quantity needed for the 

density of Covxy() is its variance. The variance of the sample 
cross-correlation function is given by Van Trees [17] as 

𝑉𝑎𝑟{�̂�𝑥𝑦(𝑠)}

≈
1

𝑛
∑ [𝑅𝑥𝑥(𝑟)𝑅𝑦𝑦(𝑟)

∞

𝑟=−∞

+ 𝑅𝑥𝑦(𝑟 + 𝑠)𝑅𝑦𝑥(𝑟 − 𝑠)] 

(27) 

and the variance of sample cross-covariance function is 

𝑉𝑎𝑟{𝐶�̂�𝑣𝑥𝑦(𝑠)} ≈
1

𝑛
∑ [

𝐶𝑜𝑣𝑥𝑥(𝑟)𝐶𝑜𝑣𝑦𝑦(𝑟)

+𝐶𝑜𝑣𝑥𝑦(𝑟 + 𝑠)𝐶𝑜𝑣𝑦𝑥(𝑟 − 𝑠)
]

∞

𝑟=−∞

 (28) 

Now, the simple model considered here for pre-processed 
discrete-time signals received at two spatially separate sensors is 
given by 

{
𝑥(𝑖) = Δ𝑡(𝑖)

𝑦(𝑖) = Δ𝑇(𝑖) = Δ𝑡(𝑖 − 𝑚) + Δ𝜏(𝑖 − 𝑚)
 (29) 

and so 

{
 
 

 
 Cov𝑥𝑥

(𝑟) = Cov∆𝑡∆𝑡(𝑟)

Cov𝑦𝑦(𝑟) = Cov∆𝑇∆𝑇(𝑟) = Cov∆𝑡∆𝑡(𝑟) + Cov∆𝜏∆𝜏(𝑟)

Cov𝑥𝑦(𝑟 + 𝑠) = Cov∆𝑡∆𝑇(𝑟 + 𝑠) = Cov∆𝑡∆𝑡(𝑟 + 𝑠 − 𝑚)

Cov𝑦𝑥(𝑟 − 𝑠) = Cov∆𝑇∆𝑡(𝑟 − 𝑠) = Cov∆𝑡∆𝑡(𝑟 − 𝑠 + 𝑚)

 (30) 

Substituting equation (30) into equation (28) yields 

𝑉𝑎𝑟{𝐶�̂�𝑣𝑥𝑦(𝑠)}

≈
1

𝑛
∑ {

𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟)[𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟) + 𝐶𝑜𝑣Δ𝜏Δ𝜏(𝑟)]

+𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟 + 𝑠 − 𝑚)𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟 − 𝑠 + 𝑚)
}

∞

𝑟=−∞

 . 
(31) 

Suppose that s=m 

𝑉𝑎𝑟{𝐶�̂�𝑣𝑥𝑦(𝑠 = 𝑚)}

≈
1

𝑛
{
𝐶𝑜𝑣Δ𝑡Δ𝑡(0)[𝐶𝑜𝑣Δ𝑡Δ𝑡(0) + 𝐶𝑜𝑣Δ𝜏Δ𝜏(0)]

+𝐶𝑜𝑣Δ𝑡Δ𝑡(0)𝐶𝑜𝑣Δ𝑡Δ𝑡(0)
}

≈
1

𝑛
[2𝑉𝑎𝑟{Δ𝑡} + 𝑉𝑎𝑟{Δ𝜏}]𝑉𝑎𝑟{Δ𝑡}  

(32) 

because 

{
𝐶𝑜𝑣Δ𝑡Δ𝑡(0) = 𝑉𝑎𝑟{Δ𝑡}

𝐶𝑜𝑣Δ𝜏Δ𝜏(0) = 𝑉𝑎𝑟{Δ𝜏}
 . (33) 

Correspondingly, if sm, 

𝑉𝑎𝑟{𝐶�̂�𝑣𝑥𝑦(𝑠 ≠ 𝑚)}

≈
1

𝑛
∑ {

𝐶𝑜𝑣Δ𝑡Δ𝑡
2 (𝑟) + 𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟)𝐶𝑜𝑣Δ𝜏Δ𝜏(𝑟)

+𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟 + 𝑠 − 𝑚)𝐶𝑜𝑣Δ𝑡Δ𝑡(𝑟 − 𝑠 + 𝑚)
}

∞

𝑟=−∞

≈
1

𝑛
[𝐶𝑜𝑣Δ𝑡Δ𝑡

2 (0) + 𝐶𝑜𝑣Δ𝑡Δ𝑡(0)𝐶𝑜𝑣Δ𝜏Δ𝜏(0)]

≈
1

𝑛
𝑉𝑎𝑟{Δ𝑡}[𝑉𝑎𝑟{Δ𝑡} + 𝑉𝑎𝑟{Δ𝜏}] . 

(34) 

According to equation (26) and the previous discussion, we 
obtain for the cross-covariance criterion function 

𝑧(𝜃) ≈
𝜃√2(

𝜎Δ𝑡
𝜎Δ𝜏

)
2

+ 1 + (
𝜎Δ𝑡
𝜎Δ𝜏

)√𝑛

√(
𝜎Δ𝑡
𝜎Δ𝜏

)
2

+ 1

=
𝜃√2[(

𝜎Δ𝑡
𝜎𝜏
)
2

+ 1] + (
𝜎Δ𝑡
𝜎𝜏
)√𝑛

√(
𝜎Δ𝑡
𝜎𝜏
)
2

+ 2

 

(35) 

where =2, ²=Var{} and t²=Var{t}. 
If we compare equation (25) and equation (35) to the results 

of Chan, Yansouni, and Crozier in [20], we obtain formal 
equivalence between the signal-to-noise ratio (SNR) in the 

conventional cross-correlation method and the (t/)-ratio or 

the (t/)-ratio in the TDD method. 
According to equation (25) and equation (35), Pa is constant 

when (t/)n is constant, and (t/)²<<1. This means 

that the distance between two Pa curves in the (t/)-scale with 

large sample sizes of n1(≫ 1) and n2(≫ 1) is approximately 
10 log(n2/n1) in decibels, as we clearly see in Figure 6, in which 
theoretical Pa curves computed with equation (25) and equation 
(35), and the corresponding experimental relative frequencies of 
anomalous estimates are presented. For example, the probability 
curves with sample sizes of 1000 and 3000 correspond to each 

other with a 5 dB shift (approximately) in respect to the (t/)-



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 82 

ratio, as we see in Figure 6. Each relative frequency point in 
Figure 6 was obtained from a series of 60 simulation experiments 
with sample sizes of 300, 1000, and 3000 per experiment. This is 
in good agreement with the theory presented previously. 

3.2 Threshold effect 

It is well known that below a certain signal-to-noise ratio, the 
estimation variance increases rapidly and strongly compared to 
the theoretical variance predicted by linear and high signal-to -
noise ratio analysis. This is due to the threshold effect caused by 
anomalous estimates. Of course, the CRLB is no more accurate 
an estimate of performance for time delay estimators when the 
probability of occurrence of an anomaly is sufficiently increased 
from zero. Several studies have been conducted to find a bound 
that is tighter than the CRLB and that would predict 
performance more accurately in the case of small signal-to-noise 
ratios [18]-[23]. For example, it has been evidenced that the Ziv-
Zakai Lower Bound (ZZLB) developed in [21] is a tighter bound 
than CRLB for time delay estimators. In fact, ZZLB is very 
nearly the greatest lower bound for the dual-sensor delay 
estimation problem [22]. 

Generally, the variance in the time delay estimate is 
characterised by three distinct regions [24]. If the probability of 
occurrence of an anomaly in the TDD method is increased and 

approaches /(+1), prior information (e.g. the known 
maximum delay) limits the variance. This is equivalent to a pure 
guess of the index delay in the selected delay window. Obviously, 
then, we can simply assume a uniform distribution of anomalies 
between the limits of the delay window as Ianniello does in [18]. 
If the probability of occurrence of an anomaly approaches zero, 
the variance of the time delay estimate approaches CRLB. 
Between these limits, there is a transition region from the prior 
information limit to the linear mode or CRLB. The theoretical 
variance of the time delay estimate in the presence of anomalous 
estimates can be determined from Ianniello’s correlator 
performance estimate (CPE) model [18], which is modified here 
to  

𝑉𝑎𝑟{⟨𝜏⟩} = (1 − 𝑃𝑎)𝑉𝑎𝑟{⟨𝜏⟩}𝑚𝑖𝑛 + 𝑃𝑎
[(𝜔 + 1)𝐸{Δ𝑡}]2

12
 (36) 

where Var{ }min is a variance of  sample mean   below the 
threshold. 

In normalised form, the variance of   is 

𝑉𝑎𝑟{< 𝜏 >}𝑛𝑜𝑟𝑚 =
(1 − 𝑃𝑎)𝑉𝑎𝑟{< 𝜏 >}𝑚𝑖𝑛

[(𝜔 + 1)𝐸{Δ𝑡}]2

12

+ 𝑃𝑎 . (37) 

An interval between two Var{< >}norm curves in (t/)-
scale with large sample sizes of n1 and n2 is (approximately) 
10 log(n2/n1) in decibels. According to the previous discussion of 
constant Pa values, Pa-curves intersect the corresponding (same 

sample size) Var{< >}min curves (linear parts of Figure 7) with 

an approximately constant value of normalised Var{< >} as 

long as (t/)²<<1 or (t/)²<<2, which necessarily means 
that the sample size is large because the constant threshold value 
of Pa is much smaller than one. 

Thus, the threshold variance, or threshold time delay error, is 
approximately constant independent of the sample size. The 

threshold (t/)-ratio, at which the variance of the time delay 
estimate begins to deviate from the minimum variance 

corresponding to CRLB, is approximately proportional to 1/n 
or to –10 log(n) in decibels. These results are equivalent to the 
results of Scarborough, Tremblay, and Carter [24], [25]. They 
observed that the value of variance at which ZZLB begins to 
deviate from CRLB remains essentially constant, independent of 
the coherent processing time, with a large bandwidth-time 
product. The threshold signal-to-noise ratio is approximately 
inversely proportional to the square root of the coherent 
processing time for the large bandwidth-time product. 

Figure 7 demonstrates the sample sizes of 300, 1000, 3000, 
and 10000 for which the theoretical threshold variances stay 
essentially constant with large sample sizes. Figure 7 shows that 
a sample size of 100 samples is not sufficiently large. The 

threshold value condition (t/)
2<<2 is strongly violated. 

Sample variances and theoretical variances of < > correspond 

quite well. Only the experimental threshold (t/)-ratio seems 
to be decreased in the order of 5 dB compared to the theoretical 
threshold. Clearly, this effect is caused by the small number of 
experiments. For example, the probability that at least one 
anomalous estimate occurs in the nearest neighbourhood of the 

threshold (10-6 in Figure 7) with 60 experiments is 0.00006 
assuming the experiment outcome ‘anomalous estimate’ follows 
the binomial distribution. 

The threshold effect has particular significance concerning 
coherent vs. incoherent signal processing for time delay 
estimation. Incoherent processing is an attractive alternative to 
coherent processing in many time delay applications of 
conventional correlation techniques because it becomes 
necessary to compensate for the effect of variations in time delay 
(e.g. due time-varying time delay, source or receiver motion) for 
coherent processing times other than short ones [26], [27], [28]. 

In incoherent processing, the whole sample size is divided 
into sections. The sections are processed individually to obtain 
one estimate for each section. These time delay estimates are 
then averaged to obtain the final estimate at the end of the whole 
observation interval. The coherent processing sample size with 
incoherent processing is reduced to n/N samples for each of the 

Pa values as a function of (t/) 

 
t/ (dB) 

Figure 6. Probability curves and corresponding relative frequencies of 

occurrence of an anomalous estimate with =10 and sample sizes 300 (• • 

•), 1000 (* * *) and 3000 () samples. Each relative frequency point was 
obtained from series of sixty simulation experiments with sample sizes of 
300, 1000 and 3000 per experiment. 

0,0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

-40 -35 -30 -25 -20 -15 -10 -5 0



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 83 

N sections, although the total observation sample size is still n 
samples. 

If the actual (t/)-ratio is larger than the (t/)n/N-ratio, 

which is the threshold (t/)-ratio for the coherent processing 
sample size of n/N samples, the performance of the incoherent 
processor (according to previous threshold discussions) 
coincides with the coherent processor and with small error 
performance. However, if 

(
𝜎Δ𝑡
𝜎𝜏
)
𝑛

< (
𝜎Δ𝑡
𝜎𝜏
) < (

𝜎Δ𝑡
𝜎𝜏
)
𝑛
𝑁

 (38) 

the incoherent processor exhibits significantly poorer 
performance than the coherent processor. If n/N is large, the 

threshold (t/)-ratio for the coherent processing sample size 
of n/N samples is 

(
𝜎Δ𝑡
𝜎𝜏
)
𝑛
𝑁

≈ (
𝜎Δ𝑡
𝜎𝜏
)
𝑛

√𝑁 . (39) 

Thus, significant performance gains could be obtained for 

small (t/)-ratios by increasing the coherent processing time. 

Note that the analogy between the (t/)-ratio in TDD and the 
signal-to-noise ratio in GCC continues, as we can see that 
comparing equation (39) to the following equation holds for 
large bandwidth-time products [24] 

𝑆𝑁𝑅𝑇
𝑁
≈ 5𝑙𝑜𝑔(𝑁) + 𝑆𝑁𝑅𝑇 (dB) (40) 

where SNRT is the threshold SNR for a coherent processing time 

of T, and SNRT/N is the threshold SNR for a coherent processing 

time of T /N.  
In the TDD method, it is not necessary to compensate for the 
effect of moderate time variations in the time delay because the 
constant index delay giving the extreme value to the criterion 
function also contains information on the variable time delay and 
determines the best estimate for time delay distribution within 
the meaning of that criterion. In the conventional correlation 
method with coherent processing, there is only a constant time 
delay argument, which maximises the cross-correlation function 
without noticing a variable time delay. Of course, if the time delay 
is strongly and constantly changing, it may destroy a sufficient 
statistical similarity (correlation) between the time sequences in 
TDD with coherent processing so that only random noise and 

sensor error effects remain. Thus, the cessation of the error 
propagation coherent processing time in TDD is also bound. 
Then, a prerequisite for the successful utilisation of incoherent 

processing is the condition (t/)>(t/)n/N. 

3.2 Distortion of time delay distribution in flow systems 
disobeying the FIFO order relation 

The theoretical values of the previous criterion functions were 

derived assuming stationary stochastic processes generating 

independent time differences and FIFO-type order relations 

between sensors. 
If the flow system is not a FIFO system, time delay 

distribution, coherently computed with the help of information 
on two time sequences obtained from two sensors a fixed 
distance apart, is a distorted version of true time delay 
distribution. Such a situation may, for example, be in traffic 
applications (a multilane motorway with unlimited speed as a 
pathological example) where the distance between the sensors is 

sufficiently long that this relation t <<  is true and the delay 
time distribution of the vehicles is sufficiently wide that this 

relation  << t is not true, the consequence being that the 
time delay distribution detected by the measurement system is 
distorted (without identification of the vehicles and their arrival 
and departure times).  

The results of applying the theory presented in the previous 
sections are fully applicable if the distortion is small. Obviously, 
the distortion is small if the following expression 

𝑃0 = ∫ 𝑓Δ𝑡(Δ𝑡)∫ 𝑓Δ𝜏(Δ𝜏)𝑑Δ𝜏𝑑Δ𝑡
𝛥𝑡

−𝛥𝑡

∞

−∞

 (41) 

where 

𝑓Δ𝜏(Δ𝜏) = ∫ 𝑓𝜏(𝜏)𝑓𝜏(Δ𝜏 − 𝜏)
∞

−∞

𝑑𝜏 

is near to one. P0 describes the probability per sample that the 
FIFO order relation is not violated. Quite arbitrary, we could 
require that the condition of the small distortion is valid if P0 is 
in the order of 0.9 or more. 

To obtain a simple approximation, we replace ft with Dirac's 

delta function D() 

𝑓Δ𝑡(Δ𝑡) = 𝛿𝐷(Δ𝑡 − 𝐸{Δ𝑡}) . (42) 

Then 

𝑃0 = 𝐹Δ𝜏(𝐸{Δ𝑡}) − 𝐹Δ𝜏(−𝐸{Δ𝑡}) (43) 

where F() is a cumulative distribution function of random 

variable . Assuming Gaussian variable  equation (43) 
becomes the simple form 

𝑃0 = 2𝑒𝑟𝑓 (
𝜇Δ𝑡

𝜎𝜏√2
) (44) 

where erf() is the error function  

𝑒𝑟𝑓(𝑥) =
1

√2𝜋
∫ 𝑒

−
𝑦2

2 𝑑𝑦
𝑥

0
. 

Now, the condition of the previously determined small distortion 

holds if µt/2.33. 

For Gaussian variables t and  

𝑃0 = 2 𝑒𝑟𝑓 (
𝜇Δ𝑡

√𝜎Δ𝑡
2 + 2𝜎𝜏

2
) (45) 

Var{< >}norm    (dB ) 

 
t/ (dB) 

Figure 7. Normalised theoretical and sample variances of <> as a function of 

the (t/)-ratio. The sample size is the parameter. The theoretical variance 
curves are numerically computed using equation (25), equation (35), and 

equation (37) with =10, t=30, and t=1. Each sample variance is 
calculated from 60 simulation experiments using sample sizes n=100 (* * *), 

n=300 (• • •), and n=1000 (). 

-70

-60

-50

-40

-30

-20

-10

0

-30 -25 -20 -15 -10 -5 0 5 10
-70

-60

-50

-40

-30

-20

-10

0

-30 -25 -20 -15 -10 -5 0 5 10

n=10000 n=3000 n=1000 n=300 n=100



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 84 

and the condition of small distortion holds if 

𝜇𝛥𝑡

√𝜎𝛥𝑡
2 + 2𝜎𝜏

2
≥ 1.65 . 

As a measure of distortion, we could use the following 
expression 

𝐷 = ∑𝑖∫ 𝑓𝛥𝑡(𝛥𝑡)∫ 𝑓𝛥𝜏(𝛥𝜏)𝑑𝛥𝜏
−𝑖𝛥𝑡

−(𝑖+1)𝛥𝑡

∞

−∞

∞

𝑖=1

𝑑𝛥𝑡 . (46) 

D describes the average amount of FIFO order violations per 

sample. For Gaussian variables t and  equation (46) can be 
obtained after manipulation to form 

𝐷 = ∑𝑖{𝑒𝑟𝑓[𝐴(𝑖 + 1)] − 𝑒𝑟𝑓[𝐴(𝑖)]}

∞

𝑖=1

 (47) 

where 

𝐴(𝑖) =
𝑖𝜇Δ𝑡

√𝑖𝜎Δ𝑡
2 + 2𝜎𝜏

2
 .  

In Figure 8, we see that the condition for the small distortion 

holds if the distortion percentage (100 % ∙  D) is about 5 % or 
smaller. 

The expectation of the sample variance of the estimated time 
delay is reduced due to distortion 

𝐸{𝑉�̂�𝑟{𝜏}} = 𝑃0𝑉𝑎𝑟{𝜏} (48) 

where P0 is obtained from equation (41). Therefore, the previous 
condition for the small distortion means that the standard 
deviation of the time delay is underestimated at approximately 5 
% or less if the idealised FIFO model is used. In Figure 8, the 

reduction of E{Vâr{}} from its true value Var{} with 

Gaussian variables t and  as a function of /µt presents 

holding t as a constant. Figure 8 shows rather good harmony 
between the experimental results and the theoretical curves 
computed from equation (45), equation (47), and equation (48). 
Each sample variance and distortion percentage were computed 
from the experiment of 1000 time samples.  

Of course, the reduced sample variance of the time delay due 
to distortion will bias the estimate of mean transport velocity. 
However, the estimate of mean transport delay is still unbiased 
because of the compensation mechanism due to the inherent 
symmetry of distortion effects. However, violation of the FIFO 
order or mixing in the space between the sensors effectively 
reduces the sufficient statistical similarity (correlation) between 
the channel signals, which is an essential prerequisite for the 
successful utilisation of the method. Therefore, a sufficient 
statistical similarity between channel signals sets such an upper 
limit for the distance between sensors that the previous condition 
for the small distortion cannot disobey much. 

3.3 Some discussions on reliability-checking possibilities 

Because the TDD method is based on the utilisation of 
criterion functions, there are several simple ways, in smart 
measurement systems, of implementing automatic reliability 
checking or supervisory operations even in real time. The trivial 
checking operation is based on magnitude testing of the extreme 
value of the criterion function. If the extreme value exceeds or 
passes under some predetermined range, it is reasonable to doubt 
the reliability of results. Of course, a great disadvantage of this 
simple checking method is the determination of the acceptable 
range of extreme values, which may be problematic. 

Obviously, it is possible to utilise the theoretical relationship 
between the extreme value of the criterion function and the 

(t/)-ratio to fix the predetermined range. For example, 

choosing the lowest acceptable limit of (t/), we get the 

lowest acceptable maximum value of the criterion function and 
vice versa. Fixing the sample size and index delay window, we 
get a clear probability measure Pa for a criterion of reliability. For 
example, the probability of an anomalous estimate as a function 
of the maximum value of the normalised cross-covariance 
criterion function ecor is 

𝑃𝑎 = 1 −
1

√2𝜋
∫ 𝑒

−
𝜃2

2

∞

−∞

{
1

√2𝜋
∫ 𝑒

−
𝑡2

2 𝑑𝑡
𝑧(𝜃)

−∞

}

𝜔

𝑑𝜃 , 

𝑧(𝜃) = 𝜃√𝑀𝑎𝑥
𝑘
{𝑒𝑐𝑜𝑟(𝑘)}

2 + 1 + √𝑛𝑀𝑎𝑥
𝑘
{𝑒𝑐𝑜𝑟(𝑘)} . 

(49) 

Figure 9 shows several Pa-curves with some sample sizes.  
The second possibility concerning the evaluation unreliability 

of measurement results is to study the sensitivity of the criterion 
function in surroundings of selected index delay. Obviously, we 
can think that a result is reliable if the extreme value differs 
sufficiently from nearby values. Determination of a sufficient 
difference between the extreme value and other values of the 
criterion function is also problematic. A quantitative measure of 
the quality of the time delay estimator in the delay window of 

(+1) index values ( is even) could be 

�̂� =
�̂�(𝑘0) −

1
𝜔
∑ �̂�(𝑘0 + 𝑖)
𝜔
2
|𝑖|=1

√1
𝜔
∑ [�̂�(𝑘0 + 𝑖) −

1
𝜔
∑ �̂�(𝑘0 + 𝑖)
𝜔
2
|𝑖|=1

]

2𝜔
2
|𝑖|=1

 
(50) 

where k0 is the estimated index delay. A similar quality measure 

(‘peakedness’) was presented by Mars and Van Arragon using the 
concept of the average mutual amount of information in the 
estimation of time delay in nonlinear systems [29]. 

We can further expand this idea. For example, in the case of 
the cross-covariance criterion function, a comparison of 

 

P0 D 

0,0 

0,1 

0,2 

0,3 

0,4 

0,5 

0,6 

0,7 

0,8 

0,9 

1,0 

0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0 2,2 

 D 
 

i 

s 

t 

o 

r 

t 

i 

o 

n 

  

 

E{Vâr{}}/Var{} 

/t 

Figure 8. Reduction of the detected variance of transport delay and the 

distortion of delay distribution with Gaussian variables (t=30, t=5). Each 
sample variance is computed from 1000 samples. Theoretical curves are 
computed from equation (45), equation (47), and equation (48). Black spots 
denote the variation range of the sample values in 20 simulation 
experiments.  



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 85 

equation (26), equation (35), and equation (50) shows that 
measure û is an estimate of 

𝑢 =
𝐸{�̂�(𝑘0)} − 𝐸{�̂�(𝑘)}

√𝑉𝑎𝑟{�̂�(𝑘)}
≈

√𝑛(
𝜎𝛥𝑡
𝜎𝜏
)

√(
𝜎𝛥𝑡
𝜎𝜏
)
2

+ 2

 (51) 

where kk0. 

In this way, we can transform this measure of ‘peakedness’ to the 
probability of occurrence of an anomaly, substituting 

𝑧(𝜃) = 𝜃√
𝑢2

𝑛
+ 1 + 𝑢 (52) 

into equation (49). 
The third possibility concerning supervision of reliability is 

more practical and is ready to use. This method utilises several 
different criterion functions at the same time. An agreement of 
criterion functions concerning index delay is the measure of 
reliability. The advantage of this method is that it is not necessary 
to determine any threshold values or ranges. 

Simultaneous use of several different criterion functions in 
reliability checking is based on the presumption that it is very 

unlikely that all different criterion functions have the extreme 
value in the same incorrect index delay. This presumption seems 
to be correct, as we see in Table 1, in which each figure describes 
the relative number of undetected errors calculated in 60 
experiments. Due to the error detection capability of combined 
criterion functions, the relative amount of the undetected error 
strongly decreases compared to the single reference criterion 
function (the main criterion function), which in Table 1, is a 
cross-covariance function êcov. 

The fourth possible method for the reliability evaluation of 
results is also practical and ready to use in smart measurement 
systems. This method utilises coherent and incoherent 
processing at the same time, continually comparing the 
magnitudes of extreme values and detecting optimum index 
delay changes, which may reveal errors caused by strongly and 
constantly changing time delays or all types of sensor errors, 
which may destroy a sufficient statistical similarity (correlation) 
between time sequences. An agreement on the optimum index 
delay is the measure of reliability. The advantage of this method 
is that it is not necessary to determine any threshold values or 
ranges. 

4. CONCLUSION 

In this time delay estimation model, the criterion function 
compares the time differences of time sequences between 
channels, not the amplitudes of time functions as in the 
conventional cross-correlation method. Using the index delay, 
which gives the extreme value to the criterion function, it is 
possible to calculate the best estimate for time delay distribution 
within the meaning of that criterion. With the help of estimated 
time delay distribution, the velocity distribution can be calculated 
when the distance between sensors is known. Using this method, 
the estimated delay distribution and criterion function are clearly 
separated. There is no need to make symmetry assumptions for 
the criterion function or to assume the criterion function as an 
estimate of the time delay distribution. Thus, there are no 
theoretical problems in the determination of average time delay 
or velocity in the non-constant or changing time delay case as 
long as a sufficient statistical similarity (correlation) exists 
between the channel signals.  

The theoretical values of several criterion functions and the 
probability of occurrence of an anomalous estimate with the 
cross-covariance criterion function are derived. The threshold 
effect caused by anomalous estimates, coherent vs. incoherent 
processing, and the distortion of estimated delay distribution are 
discussed. Some possibilities concerning reliability supervision 
based on the use of criterion functions are outlined. Potential 
applications of the time difference distribution method are the 
delay and velocity measurements in stochastic flow, or the 
transport processes and time delay estimation problems in radar, 
sonar, and wireless communication systems.  

The basis and an analysis of the estimation method using 
simple data were presented here. Topics for additional studies 
could be, for example, favourable subjects of practical 
application and real-time reliability-checking methods based on 
simultaneous use of different criterion functions and parallel 
incoherent and coherent processing in smart measurement 
systems.  

REFERENCES 

[1] G. C. Carter, Special issue on time delay estimation, IEEE Trans. 
Acoust., Speech, Signal Processing ASSP-29, 3 (1981). 

Pa (dB below 1) 

 
Max{ecor} 

n=100 n=300 n=1000

n=3000 n=10000  

Figure 9. The probability of an anomalous estimate as a function of the 
maximum value of normalised cross-covariance criterion function computed 

by equation (49), and =10. 

Table 1. The relative amount of undetected errors (anomalous estimates) as 

a function of (t/)-ratio and corresponding reference probability of an 
anomalous estimate with cross-covariance criterion function Pacov. All 
distributions are approximately Gaussian. Each figure describing the relative 
amount of undetected error of criterion functions (ecov, eadd, evd) is calculated 
from 60 experiments. 

t/ 
in dB 

Pacov 
Relative amount of undetected errors 

ecov only 
Combination of 

 eadd, ecov, and evd 

n = 100 n = 300 n = 300 ω = 10 n = 100 n = 300 n = 100 n = 300 

0 -5 ≪ 0.001 0 0 0 0 

-5 -10 0.030 0.02 0.05 0 0 

-10 -15 0.296 0.27 0.38 0.02 0.07 

-15 -20 0.606 0.57 0.62 0.13 0.05 

-20 -25 0.769 0.73 0.72 0.08 0.10 

-25 -30 0.841 0.87 0.85 0.10 0.08 
 

-180

-160

-140

-120

-100

-80

-60

-40

-20

0

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0



 

ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 86 

[2] J. S. Bendat, A. G. Piersol, Engineering Applications of 
Correlation and Spectral Analysis, Wiley-Interscience, New York, 
1980. 

[3] C. H. Knapp, G. C. Carter, The generalized correlation method 
for estimation of time delay, IEEE Trans. Acoust., Speech, Signal 
Processing ASSP-24, 4 (1976) pp. 320-327. 

[4] M. H. Butterfield, G. F. Bryant, J. Dowsing, A new method of 
strip speed measurement using random waveform correlation, 
Trans. Soc. Instrum. Tech. 13 (1961) pp. 111-123. 

[5] T. Idowaga, T. Ono, A new method of vehicle-speed 
measurement using cross correlation, Proc. of the VII IMEKO 
World Congress, London, UK, 10 - 14 May 1976, Institute of 
Measurement and Control, London pp. BLV 218/1-9. 

[6] C. Zimmer, H. Ryser, H. Meyr, An instrument for non contact 
speed measurement by image correlation, Proc. of the VII 
IMEKO World Congress, London, UK, 10 - 14 May 1976, 
Institute of Measurement and Control, London, pp. BLV 217/1-
10. 

[7] H. Kashiwagi, Paper strip speed measurement using correlation 
techniques, Trans. Soc. Instrument and Contr. Engrs., Japan, 
1968, pp. 4304-312. 

[8] M. S. Beck, Correlation in instruments: cross correlation 
flowmeters, J. Phys. E: Sci. Instrum. 14 (1981) pp. 7-19. 

[9] J. W. Arthur, J. F. Dix, E. Harland, J. W. Widdowson, P. B. 
Denyer, J. Mavor, A. D. Milne, Large time-bandwidth product 
CCD correlators for sonar, IEE Conference Publication 
1979/180, pp. 62-68. 

[10] D. E. N. Davies, Signal processing for radar applications, IEE 
Conference Publication no. 1979/180, pp. 145-150. 

[11] W. H. Haas, C. S. Lindquist, A synthesis of frequency domain 
filters for time delay estimation, IEEE Trans. Acoust., Speech, 
Signal Processing ASSP-29 (1981) pp. 540-548. 

[12] R. M. Henry, An improved algorithm allowing fast on-line polarity 
correlation by microprocessor or minicomputer, IEEE 
Colloquium on Correlation Processing, May 1974, London, UK 
pp. 3/1-4. 

[13] L. R. Morris, The role of zero crossings in speech recognition and 
processing, Proc. of the IEEE Conf. Speech Comm. Processing, 
Boston, USA, 24 – 26 April 1972, pp. 446-450. 

[14] B. Friedlander, On the Cramer-Rao bound for time delay and 
Doppler estimation, IEEE Trans. Inform. Theory, vol. IT-30, 
May 1984, pp. 575-580. 

[15] W. R. Hahn, Optimum signal processing for passive sonar range 
and bearing estimation, J. Acoust. Soc. Amer. 58 (1975) pp. 201-
207. 

[16] R. L. Kirlin, J. N. Bradley, Delay estimation simulations and a 
normalized comparison of published results, IEEE Trans. 
Acoust., Speech, Signal Processing ASSP-30 (1982) pp. 508-511. 

[17] H. L. Van Trees, Detection, Estimation, and Modulation Theory, 
Part I, John Wiley and Sons, New York, 1968. 

[18] J. P. Ianniello, Time delay estimation via cross-correlation in the 
presence of large estimation errors, IEEE Trans. Acoust., Speech, 
Signal Processing ASSP-30 (1982) pp. 998-1003. 

[19] J. P. Ianniello, Large and small error performance limits for 
multipath time delay estimation, IEEE Trans. Acoust., Speech, 
Signal Processing ASSP-34 (1986) pp. 245-251. 

[20] Y. T. Chan, P. Yansouni, S. Crozier, Time delay estimation with 
low signal to noise ratios, Proc. of the IEEE International 
Conference on Acoustics, Speech, and Signal Processing ICASSP 
‘84, 19 – 21 March 1984, San Diego, USA, pp. 616-619. 

[21] A. Weiss, E. Weinstein, Fundamental limitations in passive time 
delay estimation - Part I: Narrow-band systems, IEEE Trans. 
Acoust., Speech, Signal Processing ASSP-31 (1983) pp. 472-486. 

[22] J. P. Ianniello, E. Weinstein, A. Weiss, Comparison of the Ziv-
Zakai lower bound on time delay estimation with correlator 
performance, Proc. of the IEEE International Conference on 
Acoustics, Speech, and Signal Processing ICASSP ‘83, Boston, 
USA, 14-16 April 1983, pp. 875–878. 

[23] J. P. Ianniello, Comparison of the Ziv-Zakai lower bound on 
multipath time delay estimation with autocorrelator performance, 
ICASSP 84, Proc. of the IEEE International Conference on 
Acoustics, Speech, and Signal Processing ICASSP ‘84, 19 – 21 
March 1984, San Diego, USA, pp. 620-623. 

[24] K. Scarborough, R. J. Tremblay, G. C. Carter, Performance 
predictions for coherent and incoherent processing techniques of 
time delay estimation, IEEE Trans. Acoust., Speech, Signal 
Processing ASSP-31 (1983) pp. 1191-1196. 

[25] K. Scarborough, G. C. Carter, R. J. Tremblay, Implications of 
threshold effects for coherent and incoherent processing 
techniques of time delay estimation, Proc. of the IEEE 
International Conference on Acoustics, Speech, and Signal 
Processing ICASSP ‘84, 19 – 21 March 1984, San Diego, USA, pp. 
612-615. 

[26] C. H. Knapp, G. C. Carter, Estimation of time delay in the 
presence of source or receiver motion, J. Acoust. Soc. Amer. 61 
(1977), pp. 1545-1549. 

[27] W. B. Adams, J. P. Kuhn, W. P. Whyland, Correlator 
compensation requirements for passive time delay estimation with 
moving source or receivers, IEEE Trans. Acoust., Speech, Signal 
Processing ASSP-28, 2 (1980), pp. 158-168. 

[28] J. P. Kuhn, S. E. Robison, D. W. Winfield, Time Delay Estimation 
Techniques for Moving Sources, Proc. of the 24th Midwest Symp. 
on Circuits and Systems, Albuquerque, USA, 29 June 1981, pp. 69-
75. 

[29] N. J. I. Mars, G. W. Van Arragon, Time delay estimation in 
nonlinear systems, IEEE Trans. Acoust., Speech, Signal 
Processing ASSP-29 (1981) pp. 619-621.