Estimation of subsample time delay differences in narrowbanded ultrasonic echoes using the Hilbert
transform correlation
Anders Grennberg and Magnus Sandell
December 1993
Abstract
In many areas the time delay of arrival (TDOA) is desired. Conventional estimators usually have the property of breaking down for low SNR. In the case of narrowbanded signals we propose a fast and simple method to estimate the small time delays. The method is shown to have the same accuracy as conventional methods for small delays and high SNR while it performs better for low SNR. It is based on using the Hilbert transform in correlation between two signals. Apart from using it for better resolution in TDOA, it may also be used in applications with narrowbanded signals where the measurements are repeatable, such as ultrasonic imaging and ECG. In ultrasonic applications, due to fluctuations in the insonified media, a small random time shift may be present causing the signals to be misaligned in time. Averaging signals under these conditions will result in a distortion of the signal shape. We propose an averaging method to avoid this and to accomplish a higher SNR without the distortion. Simulations and experiments from ultrasonic applications are presented.
1 Introduction
The time delay of arrival (TDOA) is required in many areas, e.g. ultrasonic imaging [16]. Estimators of the TDOA are designed in different ways [3, 5, 7, 10, 17], but they have a few things in common. Usually they have an SNR threshold under which they perform substantially less well [18]. For correlator based estimators this is due to peak ambiguity, i.e. the possibility of selecting the wrong peak of the correlation. These estimators also need some kind of interpolation to obtain subsample resolution, e.g. a parabolic fit [4, 14]. In [11], a comparison of different methods is made. Under the assumption that the signals are narrowbanded and that the time delays are small, we propose a novel method to estimate the TDOA. This method is shown to be fast and simple. Apart from the increased resolution in TDOA, it can be used for other purposes.
One such application is the problem of averaging. In areas where the measurements are repeatable, such as ultrasonic imaging, nondestructive evaluation (NDE) and ECG, the process of averaging can be used to improve the SNR. If the useful part of the signal is the same for all measurements, it will remain the same after the averaging while the disturbances (e.g. noise and artifacts) will be attenuated since they occur randomly in each separate measurement. However, if the useful part of the signal has a small random time shift due to physical reasons or inaccuracies in the measurement equipment, the averaging will distort this signal. This phenomenon can appear e.g. in ultrasonic imaging when the insonified medium is not homogeneous but varies slightly with time. Using the above mentioned method we can estimate these small random time shifts and design an averaging algorithm that time aligns (i.e. causes the useful part of the signals to have the same position in time) so that very little signal shape distortion occurs when the averaging is carried out. Another area of application is in ultrasonic pulse-echo imaging.
Airborne ultrasound suffers from inhomogeneities in the medium such as temperature gradients and air turbulence [8]. These effects can be reduced if a reference reflector, such as a line reflector across the ultrasonic field, is placed at a fixed distance from the transducer [9]. The main echo from the surface will not be disturbed in any major way and changes in the medium can be estimated since the distance to the reference reflector is known. Since these changes can be expected to be small and the reference echo, due to its limited size of the reflector, have a poor SNR, the algorithm described in this paper is very useful.
2 Derivation of the estimator
We model all received echoes as having the same signal shape and only differing in the arrival time. Thus, using r(t) as a reference echo we find that the received echo is
s(t) = r(t− θ) (1)
where the time delayθ is considered to be an unknown deterministic parameter. In many conventional time delay estimators [3, 5] the cross correlation between the reference echo r(t) and the received echo s(t) is considered. However, using the fact that the echo is narrowbanded and the time delay θ is assumed to be small, an alternative approach is to use the correlation with the Hilbert transform of r(t). The Hilbert transform is defined as [2]
ˇ
r(t) =H{r(t)} = 1 π
Z ∞
−∞
r(x)
x− tdx = h(t)∗ r(t) (2) where the integral is a Cauchy Principal Value (CPV) and ∗ denotes convolution. The Hilbert kernel is denoted by h(t) = −πt1. We observe that if r(t) ∈ L2(−∞, ∞) then it follows that ˇr(t) ∈ L2(−∞, ∞), see [6]. Now, consider the cross correlation between s(t) and ˇr(t) for time lag zero, i.e. the standard scalar product between s(t) and ˇr(t) and use Parseval’s identity:
Z ∞
−∞s(t)ˇr(t)dt = < s, ˇr > =
Z ∞
−∞r(t− θ) [h(t) ∗ r(t)] dt
=
Z ∞
−∞e−j2πfθR(f )H(f )R(f )df =
Z ∞
−∞ej2πf θH(f )|R(f)|2df (3) If the ultrasonic echoes are narrowbanded, i.e. their energy is concentrated in frequency intervals of size B around±f0, we may formulate (3) as
< s, ˇr >≈Z −f0+B/2
−f0−B/2 ej2πf θH(f )|R(f)|2df +
Z f0+B/2
f0−B/2 ej2πf θH(f )|R(f)|2df
≈ −je−j2πf0θZ −f0+B/2
−f0−B/2 |R(f)|2df + jej2πf0θ
Z f0+B/2
f0−B/2 |R(f)|2df
= Er
2 j³−e−j2πf0θ+ej2πf0θ´=−Ersin 2πf0θ =−Ersinω0θ (4) where Er is the energy of r(t) and ω0 = 2πf0. We have used the fact that H(f ), the Fourier transform of the Hilbert kernel, is j sgn(f ), see [2], and assumed the exponen- tional function to be constant over the integration interval. This gives us
sinω0θ ≈ − < s, ˇr >
E = −R−∞∞ s(t)ˇr(t)dt
R∞
−∞r2(t)dt (5)
This estimator can be expected to work well only for narrowbanded signals and small time delaysθ. The assumption that ej2πf θ is constant over the integration interval will not
be valid otherwise. However, to solve this problem we may design a ”hybrid” algorithm consisting of two parts. The first is a rough estimate of the time delay which is robust but not very accurate. Its resolution is typically a few samples which is sufficient since the ”fine tuning” can be done with the second part of the algorithm, i.e. our proposed method. As we work with sampled signals, we may give the discrete version of (5).
Assuming we haveN samples and denoting these samples s[n] = s(nTs), where Ts is the sample interval, we may approximate the integrals of (5) with sums and yield
θ =ˆ 1 ω0
sin−1
Ã−R−∞∞ s(t)ˇr(t)dt
R∞
−∞r2(t)dt
!
≈ 1 ω0
arcsin
Ã−Ts
PN−1
k=0 s[k]ˇr[k]
Ts
PN−1
k=0 r2[k]
!
(6) Using vectors and denoting them s = [s[0] s[1] . . . s[N − 1]]T, ˇr = [ˇr[0] ˇr[1] . . . ˇr[N − 1]]T and r = [r[0] r[1] . . . r[N− 1]]T, we have
θ =ˆ 1 ω0
arcsin
Ã
−sTˇr rTr
!
(7) Since r(t) is a reference echo, we may calculate −ˇr/rTr = b one time for all and arrive at the fast and simple time delay estimator
θ =ˆ 1 ω0
arcsin³sTb´ (8)
If the angular centre frequency ω0 is not known, it may be estimated by ”inducing our own time delay”. Introducing a constructed signal x(t) = r(t− Ts) we have θ = Ts. In vector notation this is equal to x = [0 r]. For this knownθ we can rewrite (7) to estimate ω0
ˆ ω0 = 1
Ts
arcsin³[0 r]T[b 0]´ (9)
3 Performance of the estimator
3.1 The influence of bandwidth
In the derivation of the estimator we assumed an infinitesimally small bandwidth. This is of course not the case for an arbitrary narrowbanded signal. To analyze the influence of bandwidth we assume that the signalr(t) has the following Fourier transform:
R(f ) =
( qEr
2B |f ± f0| < B
0 otherwise (10)
This signal has energyErand bandwidthB. In Appendix B we show that the bandwidth introduces a bias in the form of an multiplicative term:
sinω0θ =ˆ −R−∞∞ s(t)ˇr(t)dt
R∞
−∞r2(t)dt = sinc(Bθ) sin ω0θ (11) We see that in order to keep this bias low, the product of bandwidth and time-delay must be small. This gives us an idea of how large θ this estimator can be assumed to work well for.
3.2 Noise error
Assume that r(t) and s(t) are corrupted by independent additive white Gaussian noise with power spectral density ofN0/2. Since both r(t) and s(t) are bandlimited we can also assume that noise has a limited bandwidth. If this is not the case, the received signals can always be lowpass filtered without any loss of information. Note that the Hilbert transform of a white Gaussian process n(t) is also white and Gaussian:
Snˇ(f ) = Sn(f )|H(f)|2 =Sn(f ) (12) Here we let Sn(f ) denote the power spectral density of the stochastic process n(t), i.e.
the Fourier transform of the autocorrelation function Rn(τ ). If we include the noise in (5) we get
θ =ˆ 1 ω0
arcsin
Ã−R0T (s(t) + n(t)) (ˇr(t) + ˇn0(t)) dt
RT
0 (r(t) + n0(t))2dt
!
(13) where the reference echo is r(t) + n0(t) and the received echo s(t) + n(t). Since we are interested in θ and we work with discrete data, i.e. N samples of respective signal, we approximate the integrals in (13) with sums:
θ =ˆ 1 ω0
arcsin
Ã−Ts
PN−1
k=0 (s[k] + n[k])(ˇr[k] + ˇn0[k]) Ts
PN−1
k=0 (r[k] + n0[k])2
!
= 1 ω0
arcsin
µB A
¶
(14) The sampled white Gaussian processesn[k] and n0[k] have the variance σ2 =N0Fs, where Fs is the sample frequency, see Appendix A. Since this estimator is highly nonlinear, it is
difficult to analyze. However, since the mean is much larger than the standard deviation for the denominator, we may linearize the estimator in (14), see Appendix C, to obtain:
E{ˆθ} ≈ 1 ω0
arcsin
õB
µA
!
E{³θˆ− θ´2} ≈ µ2AVar{B} + µ2BVar{A}
µ2A(µ2A− µ2B)ω20 +
à 1 ω0
arcsin
õB
µA
!
− θ
!2
(15) where
µA = E{A} = Er+N Tsσ2 µB = E{B} = Ersinω0θ Var{A} = 4ErTsσ2+ 2N Ts2σ4
Var{B} = 2ErTsσ2+N Ts2σ4 (16) The mean of the estimate is
E{ˆθ} = 1 ω0
arcsin
µ Ersinω0
Er+N Tsσ2
¶
= 1 ω0
arcsin
µ SN R
SN R + 1sinω0θ
¶
(17) where we define
SN R = Er
N Tsσ2 (18)
This is not the standard definition of the SNR but rather the relationship between the energies of the signal and the noise. As we see from (17), the estimator can be expected to be asymptotically unbiased for large SNR but contain a bias for small SNR. As for the error variance, it can be shown from (15) that
lim
SN R→−∞dBE{(ˆθ − θ)2} = 1
N ω02 +θ2 lim
SN R→∞dBE{(ˆθ − θ)2} = 0 (19)
From this we see that the standard deviation of the error for very low SNR will be approximately equal to the delay θ. Generally, for an arbitrary SNR, we will have an increase in error variance for larger θ.
3.3 Numerical simulations
In order to verify the performance given by the previous section we simulated (7) numer- ically by using a synthetic echo in the form of
r(t) = Ate−at2cosω0t (20)
with the parameters A = 7 · 105, a = 1010, ω0 = 2π106 and the sample frequency fs = 20 MHz which corresponds to our equipment. The synthetic echo is shown in Fig.
-3 -2 -1 0 1 2 3
0 5 10 15 20 25 30 35 40 45 50
Time [ms]
Amplitude [V]
Figure 1: Synthetic echo used in simulations.
1. Since we usedfs = 20 MHz and the centre frequencyf0 = 1 MHz, we have 20 samples per period in the echo. For each value of the SNR, we simulated 1000 measurements and calculated the mean of the estimator (7). We added noise to both r(t) and s(t).
The simplified expression of the mean and variance (15) was compared to the mean and variance of the simulation. In Fig. 2 and 3, we show these curves as a function of SNR with θ = Ts. As we can see, the simplified expression (15) agrees very well with the simulations. The curves are in fact hard to tell apart. To see how the error variance varies with θ, we chose to calculate the performance of the estimator for three different values of the SNR: -10 dB, 10 dB and 30 dB. The error variance of the estimator is shown as a function of θ in Fig. 4. As we can see, the larger θ is, the larger the error variance becomes. This phenomenon stems from the fact that for large θ, the bandwidth of r(t) introduces a bias which can not be neglected. As a reference, we have also simulated a conventional time delay estimator. This type of estimator correlates the received echo with the reference echo and looks for the maximum of the cross correlation.
Its performance has been discussed in [7]. In order to obtain subsample resolution, a simple interpolation technique was used. A second order polynomial was fitted to the maximum correlation and its preceding and following point. The time lag that maximizes this polynomial was chosen as the subsample time lag that gives maximum correlation.
Why this works is easily seen using (4). If we correlate the received echo s(t) with r(t) instead of ˇr(t), the sin ω0θ-term will be replaced by a cos ω0θ-term. This implies that the cross correlation will have the shape of a cosine in the neighbourhood of its maximum. For smallx we can make the approximation cos x≈ 1 + x2/2, i.e. a second order polynomial, which explains why the interpolation works so well. In Fig. 5 we have compared this
SNR [dB]
E {q} [samples]^
0 0.2 0.4 0.6 0.8 1
-30 -20 -10 0 10 20 30 40 50
Figure 2: Simulated (solid) and theoretical (dashed) mean of the estimator for θ = Ts. Note that the two curves are difficult to distinguish between.
estimator with ours for the time delay corresponding to a sample interval (θ = Ts). The reason for the increase in the error variance for the conventional estimator is that it is based on peak detection [18]. For small SNR, the possibility of detecting the wrong peak increases. This is especially true for narrowbanded signals, which have a small difference in amplitude for the peaks in the autocorrelation. If we let the conventional estimator know that the time delay is small, i.e. it is not necessary to perform a global search of the maximum, this phenomenon disappears and the performance will be approximately the same as for our method. However, the complexity of the estimator is greater since the cross correlation for several time lags has to be calculated and an interpolation needs to be done to achieve subsample resolution.
SNR [dB]
E {(q-q) } [samples ]^ 2 2
10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101
-30 -20 -10 0 10 20 30 40 50
Figure 3: Simulated (solid) and theoretical (dashed) error variance of the estimator for θ = Ts
10-5 10-4 10-3 10-2 10-1 100 101 102
0 0.5 1 1.5 2 2.5 3 3.5 4
q [samples]
E {(q-q) } [samples ]^ 2 2
SNR = -10 dB
SNR = 10 dB
SNR = 30 dB
Figure 4: Simulated error variance of the estimator for different values of the SNR.
SNR [dB]
E {(q-q) } [samples ]^ 2 2
10-7 10-4 10-1 102 105 108
-30 -20 -10 0 10 20 30 40 50
Figure 5: Simulated error variance of the estimator (solid) compared with a conventional method (dashed) for θ = Ts.
4 Application example: Averaging algorithm
Using the estimator of small time delays (7), we can proceed to formulate an averaging algorithm. This method becomes useful when the echoes have small random time shifts.
These could originate from fluctuations in the insonified medium. In our line of research, which contains surface profiling by using airbourne ultrasound [16], these fluctuations take the form of temperature variations and air turbulence. A change in the air temper- ature of 0.1o C results in a 0.06 m/s change of sound speed [12]. Using our transducer with a radius of curvature of 45 mm, i.e. a total path length of 90 mm for a surface in the focal plane, this means a time delay of 45 ns = 0.95Ts in our case, where Ts is the sample interval. These small random time shifts have to be taken into account when an averaging process is applied. Otherwise we will introduce a distortion of the signal due to the misalignment of the echoes. Assume that we want to average over M echoes,
si(t) = r(t− θi) +ni(t) i = 0 . . . M − 1 (21) where we define s0(t) =∆ r(t) + n0(t). In (21), θi are small independent random time shifts and ni(t) white Gaussian additive noise which are independent of θi. The first thing is to estimate the individual time delays between the echoes. This can be done as described in section 2. We use the first echo s0(t) as the reference echo. To be able to average over the M echoes, we need to align the echoes in time. This has been done by using an interpolation formula, which is approximately an all-pass filter with linear phase (corresponding to the time delay). Since we have samples of the echoes the problem is expressed as having si(nTs) and wanting si(nTs+θi) = s0(nTs). We use the same notation as before, i.e. the samples are si[n] = si(nTs). A simple interpolation method is Lagrange’s formula [1] of the order 2N + 1
f (x + t∆) =
XN k=−N
hk(t)f (x− k∆) (22)
which uses the fact that f (x) is known in the equally spaced points xk = x− k∆. The weightshk depend on t which is the fraction of the distance ∆ for which we want to find f . The interpolation (22) can be seen as a filter with taps hk. We can now align si[n]
with s0[n] by forming s0i[n] = si[n]∗ h[n] where ∗ denotes convolution. By doing this for all echoes si[n] we can average over the M echoes:
s[n] = 1 M
MX−1 i=0
s0i[n] (23)
The averaged echos[n] will have the same position in time as s0[n] since this is used as a reference echo. In order to find an absolute position in time we average over all estimated time delays
θ = 1 M
MX−1 i=0
θˆi (24)
This gives us an estimation of the correct position in time of s[n]. To obtain the final averaged signal we interpolate with the Lagrange interpolation filter with taps hk which are determined by θ.
5 Experiments
Apart from the difficulties of defining signal shape distortion, there are also practical problems. A noise-less reference echo is usually not available and thus it is difficult to measure the effect of the signal shape distortion. Instead we chose to get another quantitative measure of the performance of the algorithm. One area where it can be used is airborne ultrasound. We used a focused transducer with a radius of curvature of 45 mm, a diameter of 20 mm and a centre frequency of approximately 1.07 MHz which makes it a strongly focused transducer [13]. It was designed by Hans W. Persson, Lund Institute of Technology, for the investigation presented in [16], and is acoustically adapted for air. We applied it to a phantom in the form of a massive brass block with smooth surfaces. The block was placed perpendicular to the transducer, the position of which was fixed throughout the whole experiment. We registered 500 echoes which took about 30 seconds with our equipment. These echoes, which should be identical since they were measured at the same spot, were then investigated. A measured echo is shown in Fig. 6 and the SNR was estimated to be approximately 30 dB.
-3 -2 -1 0 1 2 3
0 5 10 15 20 25 30 35 40 45 50
Time [ms]
Amplitude [V]
Figure 6: Measured ultrasonic echo from the phantom.
6 Results
If we plot a cross-section of the registered echoes, i.e. for a fixed value of the time t we plot how the amplitude varies for the different measurements, we can see how the misalignment effect can cause problems when we average the echoes. This cross-section plot is shown in Fig. 7. If the echoes did not have a random time shift, these curves
0 0.5 1 1.5 2 2.5
0 50 100 150 200 250 300 350 400 450 500
Measurement nr Amplitude [V]
t=t1
t=t2
Figure 7: Cross-section of the measured echoes for two values oft. The uppermost curve corresponds to t = t1 which is a peak in the echo and the lowermost tot = t2 which is a zero-crossing. Total time for the 500 measurements was 30 seconds.
would be straight horizontal lines apart from the measurement noise. The uppermost
0 0.2 0.4 0.6 0.8 1
0 50 100 150 200 250 300 350 400 450 500
Measurement nr q [samples]^
Figure 8: Estimated time shifts relative to the first echo for the measurements in Fig. 7.
curve corresponds to a point in time where the echo has a peak, i.e. small derivative, and the lowermost to a point in time when the echo has a zero crossing, i.e. large derivative.
Looking at the estimations of the time shifts relative to the first echo, Fig. 8, we can see that it has the same drift as the cross-section. The magnitude of the derivative of the echo in Fig. 6 at the zero crossings is estimated to be almost 20 V/µs. From Fig. 8 and Fig. 7 we see that the delay of one sample, i.e. 50 ns, at a zero crossing results in an amplitude difference of approximately 1 V. If we linearize the echo at the zero crossing we get an amplitude of 20 V/µs × 50 ns = 1 V. From Fig. 8 we can see that it has a low-frequency shape which corresponds to the fluctuations in the medium that appear during the total measurement time of 30 seconds. To see the effects of using our averaging algorithm, we look at the cross-section variance of s(t) for each point of time, i.e. v(t) = Var{s(t)}.
We compare the variance for the original echoes with the one where we estimated the individual time shifts and interpolated the echoes in order for them to be time aligned with the first echo. As we can see in Fig. 9, the alignment is very good. The peaks in
0 0.02 0.04 0.06 0.08 0.1
10 12 14 16 18 20 22 24
Time [ms]
Var{s(t)} [V ]2
Figure 9: Variance of the cross-section before (solid) and after (dashed) time alignment.
variance for the original echoes correspond to points of time where the echoes’ derivative is large, thus yielding large variations in the cross-section. We can also see in the plot of the variance after the time alignment that the variance has a constant amplitude and does not vary with time. The small variations that are still present are mainly due to additive noise.
7 Summary
We have proposed a new method for estimating small time delays in narrowbanded signals. It is based on correlation with the Hilbert transform of the reference echo instead of the original echo. The estimator has been shown to be simple, fast and with the same or better accuracy than conventional methods. This method could be used for increased resolution in time delay of arrival (TDOA) estimation. We have also proposed an averaging algorithm which can be used when the individual signals have a small random time shift between themselves. In this case we can have the advantage of using the above estimation method to align all signals in time before averaging them. The mean of the estimated time shifts can also be used to find a better absolute position in time for the averaged signal. Typical areas of application for these methods are ultrasonic applications and ECG.
Appendix A Sampled stochastic processes
A continuous-time white noise processn(t) with bandwidth Fs/2 has the autocorrelation function (acf)rn(τ ) = N20Fssinc(Fsτ ), where N0/2 is the amplitude of the power density spectrum and sinc(x) = sin(πx)/(πx). If this process is sampled with a frequency of Fs, i.e. n[k] = n(k/Fs) =n(kTs), the samples will constitute a discrete white noise process:
E{n[k]} = E{n(kTs)} = 0
E{n[k]n[k − m]} = E{n(kTs)n((k− m)Ts)} = rn(mTs) = N0
2 Fssinc(FsmTs) = N0
2 Fs
sinπm
πm (25)
The last expression will equal zero for all integersm6= 0. For m = 0, we have E{n[k]2} = E{n(kTs)2} = rn(0) =N0Fs/2. Thus, n[k] is a white noise process with zero mean and variance σ2 =N0Fs/2.
Appendix B The effect of bandwidth
Assuming that the narrowbanded signalr(t) has the Fourier transform R(f ) =
( qEr
2B |f ± f0| < B
0 otherwise (26)
where Er denotes the energy and B the bandwidth, we have, using Parseval’s theorem and the model s(t) = r(t− θ)
Z ∞
−∞s(t)ˇr(t)dt =
Z ∞
−∞e−j2πfθR(f )jsgn(f )R(f )df =
Z −f0+B/2
−f0−B/2 −jej2πf θ|R(f)|2df +
Z f0+B/2
f0−B/2 jej2πf θ|R(f)|2df =
−jEr
2B
"
ej2πf θ j2πθ
#−f0+B/2
−f0−B/2
+jEr
2B
"
ej2πf θ j2πθ
#f0+B/2 f0−B/2
=
−jEr
2Be−j2πf0θBsinc(Bθ) + jEr
2Bej2πf0θBsinc(Bθ) = jEr
2BBsinc(Bθ)2j sin(2πf0θ) =−Ersinc(Bθ) sin(2πf0θ) (27) Thus we have
−R−∞∞ s(t)ˇr(t)dt
R∞
−∞r2(t)dt = Ersinc(Bθ) sin ω0θ Er
= sinc(Bθ) sin ω0θ (28)
Appendix C Linearizing the estimator
In order to linearize the estimator (14) we use the Taylor series of the function f (x, y) = 1
ω0
arcsin
Ãx y
!
= f (x0, y0) + (x− x0)∂f
∂x |(x,y)=(x0,y0) +(y− y0)∂f
∂y |(x,y)=(x0,y0)+. . . (29) The partial derivatives are
∂f
∂x = 1
ω0
r 1
1−³xy´2 1 y = 1
ω0
√ 1
y2− x2
∂f
∂y = 1
ω0
r 1
1−³xy´2
−x
y2 =− 1 ω0
x y
√ 1
y2− x2 (30)
Introducing εB, which is a zero-mean stochastic variable with variance Var{B}, and µB = E{B} we can approximate the estimator (14) by putting B = µB + εB and correspondingly A = µA+εA and use the above linearization at the point (µB, µA)
θ =ˆ 1 ω0
arcsin
µB A
¶
≈ 1
ω0
arcsin
õB
µA
!
+εB
1 ω0
q
µ2A− µ2B
− εA
µB
ω0µA
q
µ2A− µ2B
(31)
The mean of ˆθ is now found to be
E{ˆθ} = 1 ω0
arcsin
õB
µA
!
(32) The error variance is
E{³θˆ− θ´2} = E
1 ω0
arcsin
õB
µA
!
− θ + εB
ω0
q
µ2A− µ2B
− µBεA
ω0µA
q
µ2A− µ2B
2
=
à 1 ω0
arcsin
õB
µA
!
− θ
!2
+ Var{B}
ω02(µ2A− µ2B) + µ2BVar{A}
ω02µ2A(µ2A− µ2B) =
à 1 ω0
arcsin
õB
µA
!
− θ
!2
+µ2AVar{B} + µ2BVar{A}
ω02µ2A(µ2A− µ2B) (33) if we assume that εA and εB are uncorrelated. The mean of B is found to be
E{B} = E
(
−Ts
XN k=0
s[k]ˇr[k] + s[k]ˇn0[k] + n[k]ˇr[k] + n[k]ˇn0[k]
)
≈ −Z ∞
−∞s(t)ˇr(t)dt≈ Ersinω0θ (34)
since n[k] and ˇn0[k] are independent and have zero mean. The mean of A follows from
E{A} = E
(
Ts NX−1
k=0
r2[k] + 2r[k]n0[k] + n20[k]
)
≈Z ∞
−∞r2(t)dt + TsN σ2 =Er+N Tsσ2 (35) For the variances ofA and B we have
Var{B} = E{(B − µB)2} = Ts2E
ÃN−1
X
k=0
s[k]ˇn0[k] + n[k]ˇr[k] + n[k]ˇn0[k]
!2
= Ts2E
(N−1 X
k=0 NX−1
l=0
s[k]ˇn0[k]s[l]ˇn0[l] + n[k]ˇr[k]n[l]ˇr[l] + n[k]ˇn0[k]n[l]ˇn0[l]+
2s[k]ˇn0[k]n[l]ˇr[l] + 2s[k]ˇn0[k]n[l]ˇn0[l] + 2n[k]ˇr[k]n[l]ˇn0[l]} = Ts2
µEr
Ts
σ2+Er
Ts
σ2+N σ2σ2
¶
= 2ErTsσ2+N Ts2σ4 (36) and
Var{A} = E{(A − E{A})2} = Ts2E
ÃN−1
X
k=0
³
2r[k]n[k] + n2[k]− σ2´!
2
= Ts2E
(N−1 X
k=0 NX−1
l=0
4r[k]n[k]r[l]n[l] + n2[k]n2[l] + σ4+ 4r[k]n[k]n2[l]− 4r[k]n[k]σ2− 2n2[k]σ2o= Ts2
µ
4Er
Ts
σ2+ (N2+ 2N )σ4+N2σ4 − 2N2σ2σ2
¶
= 4ErTsσ2 + 2N Ts2σ4 (37) Here we have used the fourth moment of a zero-mean Gaussian stochastic variable E{x4} = 3σx4, [15], giving us
E
(N−1 X
k=0 NX−1
l=0
n2[k]n2[l]
)
=E
X
k6=l
n2[k]n2[l] +X
k=l
n4[k]
=
(N2− N)σ2σ2+N 3σ4 = (N2+ 2N )σ4 (38) and with the third moment, E{x3} = 0, [15],
E
(N−1 X
k=0 NX−1
l=0
n[k]n2[l]
)
=E
X
k6=l
n[k]n2[l] +X
k=l
n3[k]
= 0 (39)
References
[1] M. Abramowitz and I.A. Stegun, Eds, Handbook of mathematical functions with formulas, graphs and mathematical tables, National Bureau of Standards Applied Mathematics Series 55, U.S. Government Printing Office, Washington, DC.
[2] R.N. Bracewell, The Fourier transform and its applications, 2nd ed., New York:McGraw-Hill, 1986.
[3] J.N. Bradley and R.L. Kirlin, ”Delay estimation by expected value”, IEEE Trans.
Acoust., Speech, Signal Processing, vol. 32, no. 1, pp. 19-27, 1984.
[4] R.E. Boucher and J.C. Hassab, ”Analysis of discrete implementation of generalized cross correlator”, IEEE Trans. Acoust., Speech, Signal Processing, vol 29, no. 3, pp.
609-611, 1981.
[5] G.C. Carter, Ed., Coherence and time delay estimation, IEEE Press, New York, 1993.
[6] H. Dym and H.P. McKean, Fourier series and integrals, New York:Academic Press, 1972.
[7] H. Eriksson, Modelling of waveform deformation and time-of-flight estimation, Li- centiate thesis, 1993:15L, Lule˚a University of Technology, May 1993.
[8] J.D. Fox, B.T. Khuri-Yakub and G.S. Kino, ”High-frequency acoustic wave measure- ments in air”, Proceedings of the IEEE 1983 Ultrasonics symposium, pp. 581-584, Atlanta, USA.
[9] A. Grennberg and M. Sandell, ”Improved range resolution in imaging with narrow- band airbourne ultrasound using a small reference reflector”, Submitted to EUSIPCO -94.
[10] I.A. Hein and W.D. O’Brien, ”Current time-domain methods for assessing tissue motion by analysis from reflected ultrasound echoes - A review”, IEEE Trans. Ul- trason., Ferroelec., Freq. Contr., vol. 40, no. 2, pp. 84-102, 1993.
[11] G. Jacovitti and G. Scarano, ”Discrete time techniques for time delay estimation”, IEEE Trans. Signal Processing, vol.41, no. 2, pp. 525-533, 1993.
[12] L.E. Kinsler, A.R. Frey, A.B. Coppens and J.V. Sanders, Fundamentals of acoustics, Wiley, New York, 1980.
[13] G. Kossoff, ”Analysis of focusing action of spherically curved transducers”, Ultra- sound in Med. and Biol., vol 5, pp. 359-65, 1979.
[14] R. Moddemeijer, ”On the determination of the position of extrema of sampled cor- relators”, IEEE Trans. Signal Processing, vol 39, no. 1, pp. 216-219, 1991.
[15] A. Papoulis, ”Probability, random variables and stochastic processes”, 2nd ed., Singapore:McGraw-Hill, 1984.
[16] N. Sundstr¨om, P.O. B¨orjesson, N.-G. Holmer, L. Olsson and H.W. Persson, ”Regis- tration of surface structures using airborne focused ultrasound”, Ultrasound in Med.
and Biol., vol. 17, no. 5, 1991.
[17] H.L. van Trees, Detection, estimation and modulation theory, Part 1, USA:Wiley, 1968.
[18] A. Zeira and P.M. Schultheiss, ”Thresholds and related problems in time delay esti- mation”, In Proceedings of international conference on acoustics, speech and signal processing, pp. 1261-1264, Toronto, Canada, May 1991.
