See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/228976901

Time-delay estimation using Farrow-based fractional-delay FIR filters: filter

approximation vs. estimation errors

Article · January 2006

CITATIONS READS

11 408

3 authors, including:

Håkan Johansson Per Bertil Löwenborg

Chalmers University of Technology Signal Procession Devices
187 PUBLICATIONS 2,572 CITATIONS 69 PUBLICATIONS 1,144 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Per Bertil Löwenborg on 30 May 2014.

The user has requested enhancement of the downloaded file.

 14th European Signal Processing Conference (EUSIPCO 2006), Florence, Italy, September 4-8, 2006, copyright by EURASIP

TIME-DELAY ESTIMATION USING

FARROW-BASED FRACTIONAL-DELAY FIR FILTERS:
FILTER APPROXIMATION VS. ESTIMATION ERRORS
Mattias Olsson, Håkan Johansson, and Per Löwenborg
Div. of Electronic Systems, Dept. of Electrical Engineering, Linköpings universitet
Campus Valla, SE-581 83 Linköping, Sweden
phone: +46-13-286698, email: {matol, hakanj, perl}@isy.liu.se
web: www.es.isy.liu.se

ABSTRACT The cost function F for d0=0.25 and w0=π/4

600
This paper provides error analysis regarding filter approxi-
500
mation errors versus estimation errors when utilizing Farrow-
400
based fractional-delay filters for time-delay estimation. Fur-
300

F
ther, a new technique is introduced which works on batches
200
of samples and utilizes the Newton-Raphson technique for
100
finding the minimum of the corresponding cost function.
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
d
1. INTRODUCTION
The need for estimating time-delays between two signals Figure 1: Example of the cost function F for a sinus input
arises in many different fields, including biomedicine, com- with d0 = 0.25.
munications, geophysics, radar, and ultrasonics. In [1],
a technique utilizing Farrow-based digital fractional-delay
(FD) filters [2] was introduced for this purpose. The use of
Farrow-based FD filters has two major advantages over other where d0 is the unknown difference in delay between the sig-
delay estimation techniques working in the digital domain. nals, T is the sampling period and e1 and e2 are uncorrelated
First, it is eminently suitable to handle delays that are frac-
additive noise1 . We assume that the delay d0 is a fraction of
tions of the sampling interval. This is in contrast to cross
T and that any integer sample delay has already been taken
correlation-based methods that require additional interpola-
care of in a proper manner.
tion [3]. Second, it can handle general bandlimited signals.
This is in contrast to techniques that assume a known input Now, let v(n) act as a reference signal and let the other
signal, like a sinusoidal signal [4]. signal x(n) pass through a FD filter generating the output
In [1], the idea of using Farrow-based FD filters for y(n, d) [see (7)]. The average squared difference function
delay-estimation was proposed. However, no analysis was (ASDF) F(d) for a certain delay d over a batch of N samples
provided as to filter approximation error versus estimation er- can then be written as
ror. Such an analysis will be provided in this paper. Further-
more, a new technique is introduced that works on batches of 1 n0 +N−1
samples and utilizes the Newton-Raphson technique for find- F(d) =
N n=n∑ (y(n, d) − v(n))2 . (3)
ing the minimum of the corresponding cost function. Since 0
the fractional delay of Farrow-based FD filters is governed
by only one parameter, analytical derivatives can be derived
An estimate db of the unknown fractional delay d0 in the ref-
for this purpose. Thereby, the problems associated with the
erence signal can then be computed by minimizing F as
use of numerical derivatives are avoided.
Following this introduction, Section 2 will provide a
short introduction to time-delay estimation using FD filters, db= arg min F(d). (4)
followed by a presentation of the delay-estimation technique d
in Section 3 and an error analysis in Section 4. In Section 5
we write FD filter design and in Section 6 we verify the error An example of this cost function F for a sinusoidal input
analysis. Finally some conclusions are drawn. can be seen in Fig. 1. Ideally, the function would be equal
to 2 sin 2 ( ω (d−d 0)
), which is approximately square for small
2. TIME-DELAY ESTIMATION USING FD FILTERS 2
d, however, if noise, magnitude errors, delay errors, etc, are
Two (or more) discrete-time signals, originally coming from present, it will deviate from the square shape.
one source, might experience different delays. We model this
as 1 Due to lack of space, the effect of noise is not included in this paper.

x(n) = xa (nT ) + e1 (1) It will not affect the filter approximation error vs. estimation error which is
in focus in this paper. The noise can be reduced to a level that is negliable
v(n) = xa (nT − d0 T ) + e2 (2) compared to the filter approximation error by increasing the batch length N.
 14th European Signal Processing Conference (EUSIPCO 2006), Florence, Italy, September 4-8, 2006, copyright by EURASIP

update equation is here

F ′ (dbn )
dbn+1 = dbn + (9)
F ′′ (dbn )

where F ′ (d) and F ′′ (d) are the derivatives of F(d). The

principle of the estimator is depicted in Fig. 3. Compared
to the common Least-Mean-Square (LMS) algorithm we do
Figure 2: The modified Farrow-based FD Filter where Gk (z) not need to specify a step length since it is computed explic-
are linear-phase FIR filters. In an implementation delays itly as 1/F ′′ (dbn ). For a perfectly quadratic function this step
must added to make sure the delay is equal for all the subfil- length is optimal and only one iteration is needed. However,
ters. in a real situation, a few more iterations are needed. Each
new iteration can use a new batch of samples N or the same
batch, depending on the amount of memory at hand and if
the estimation is run on-line or off-line.
3. PROPOSED TIME-DELAY ESTIMATOR To be able to calculate the next iterative estimate in (9),
3.1 Farrow-Based FD Filters the first and second derivatives of F(n, d) with respect to d
are needed, which can be calculated as
The desired frequency response of an FD filter is
1 n0 +N−1
Hdes (e jω T , d) = e− j(D+d)ω T , |ω T | < ωc T < π (5) F ′ (d) = 2 ∑ (y(n, d) − v(n)) y′ (n, d) (10)
N n=n 0
where D is an integer delay and d is a subsample (fractional)
delay. In practice an approximation of (5), designed for fre- and
quencies up to ωc T , is used.
The FD filters in this paper make use of the modified Far- 1 n0 +N−1
F ′′ (d) = 2 ∑ (y(n, d) − v(n))y′′ (n, d) + (y′ (n, d))2 .
row structure shown in Fig. 2 [5], [6] which is a modifica- N n=n 0
tion of the original Farrow structure [2] in that the subfilters (11)
Gk (z) are linear-phase FIR filters. The advantage of using When a Farrow-based FD-filter is used, the derivatives of (7)
linear-phase filters is that their impulse responses are sym- can be calculated analytically as
metrical and can be implemented with fewer multiplications.
The overall transfer function can thus be written as L

L
y′ (n, d) = ∑ kd k−1 yk (n) (12)
k k=1
H(z) = ∑d Gk (z) (6)
k=0
and
L
where Gk (z) are linear-phase FIR filters of either odd or even y′′ (n, d) = ∑ k(k − 1)d k−2 yk (n), (13)
order, say Mk , and with symmetric (anti-symmetric) impulse k=2
responses gk (n) for k even (k odd), i.e., gk (n) = gk (Mk − n)
[gk (n) = −gk (Mk − n)] for k even (k odd) [7]. When d is respectively. In Fig. 4, a straightforward implementation of
fixed, the overall filter approximates an allpass filter with the the derivatives (10) and (11) can be seen. Note that the sub-
fractional delay d, provided that the subfilters have been de- filters Gk (z) only have to be used once.
signed in a proper manner [5], [6]. Since a Farrow-based FD filer is an approximation of the
The output y(n, d) from the FD filter in Fig. 2 can be ideal response in (5), the approximation errors will affect the
written as estimator performance. In the next section we will investi-
L gate the performance of the algorithm when the FD filter is
y(n, d) = ∑ d k yk (n). (7) not ideal.
k=0

where 4. APPROXIMATION ERROR VS. ESTIMATION

ERROR
yk (n) = gk (n) ∗ x(n). (8)
We model the non-ideal FD filter as
In this paper, it is assumed that the order M = max{Mk }
is even. The reason is that the delay of H(z) is then an in- ˜
H(e jω , d) = (1 + δ (ω , d))e− jω (D+d+d(ω ,d)) (14)
teger, D = M/2, when d = 0. This is suitable for the time-
delay estimation problem. The filter is normally optimized where δ (d, ω ) denotes the magnitude error and d(d, ˜ ω ) de-
for |d| ≤ 0.5 to cover one sampling interval. notes the delay error of the FD filter for a certain fractional
delay d and a certain angular frequency ω . D is the integer
3.2 The Estimator delay of the FD-filter. This is a formulation which is general
To find the minimum of F with respect to d, the well-known for all types of FD filters, not just the Farrow-based FD filter
Newton Raphson (NR) algorithm is used. The algorithm used in this paper, and the resulting expressions in this sec-
is iterative and tends towards the closest zero of a one- tion is hence independent of the FD filter type. In Section 5,
dimensional function, in this case the derivative of F. The expressions specific for Farrow-based FD filters are derived.
 14th European Signal Processing Conference (EUSIPCO 2006), Florence, Italy, September 4-8, 2006, copyright by EURASIP

Figure 3: The principle of the proposed estimator.

2
Accumulator F’(d)

2
(.)2 Accumulator F’’(d)
–

y0(n) v(n)
G0(z)
y(n,d)

y1(n)
G1(z)
Adder Adder Adder
x(n) y2(n)
G2(z)

L(L–1)dL–2
3x2xd
LdL–1

y3(n)
2d
d2

d3

dL
d

G3(z)

yL(n)
GL(z)

Figure 4: A straightforward implementation of the derivatives.

Now, if we assume that x(n) is sinusoidal (for the choice later, d + d˜ will be close to d0 . Using this fact we can rewrite
of signal see the conclusions), the output from the FD filter FN′ as
at a frequency ω0 is
δ sin(N ω0 ) 
˜
y(n, d) = (1 + δ ) sin(ω0 (n − d − d)) (15) FN′ (ω , d)|d+d=d
˜ 0 =− (1 + δ )ω0 (1 + d˜′ )
N sin(ω0 )

and the reference signal v(n) is sin(ω0 (2d0 + N − 1)) + δ ′ cos(2ω0 (2d0 + N − 1)) . (18)
v(n) = sin(ω0 (n − d0 )). (16) Since δ and δ ′ usually are very small, (18) will be small even
for a small N. However, when ω0 tends towards 0, (18) will
Henceforth, the dependence of δ and δ ′ on d and ω0 will be tend towards −δ ′ δ , independently of N.
omitted in the notation for the sake of clarity. The iterative Newton-Raphson algorithm in (9) will, if
After inserting (15) and (16) into (3) and some simplifi- the function is behaving well, converge towards the zero of
cations we arrive at F0′ . If we assume that N is so large that FN′ can be approxi-
F ′ (ω , d) = δ ′ [(1 + δ ) − cos(ω (d − d0 + d))]+
˜ mated as zero for all d and ω0 , it can be seen from (17) that
′
if δ ′ and d˜ are zero, F ′ will be zero when d = d0 , irrespective
+ (1 + δ )ω (1 + d˜ ) sin(ω (d − d0 + d))+
˜
of δ . On the other hand, if δ ′ is small, but not zero, the effect
1  of a constant magnitude error δ will affect the estimator. The
(1 + δ )ω0 (1 + d˜′ ) sin(ω0 (−(d0 + d + d)˜ + 2ϕ )−
effect of a delay error derivative d˜′ is normally neglible since
N
− (1 + δ )2 ω0 (1 + d˜′ ) sin(ω0 (−(d + d)
˜ + 2ϕ )+ it tends to be small compared to 1.
Unfortunately, it is impossible to directly from (17) ana-
′ ˜ + 2ϕ )−
+ δ cos(ω0 (−(d0 + d + d) lytically find the d that makes F0′ zero. It can be done numer-
 ically and to find an approximation of the estimation error we
˜ + 2ϕ ) sin(N ω0 )
−(1 + δ )δ ′ cos(ω0 (−(d0 + d) do a first-order Taylor expansion of the sine and cosine in F0′
sin(ω0 )
and write it as
′ ′
= F0 (ω , d) + FN (ω , d) (17)
ω02
F0′
where does not depend on the batch length N and con- FN′ F0′ ≈ δ ′ (δ + ˜ 2 )−(1+ δ )(1+ d˜′ )ω02 (d −d0 + d)
(d −d0 + d) ˜
2
tains the terms that do. At the minimum of F, as we will see (19)
 14th European Signal Processing Conference (EUSIPCO 2006), Florence, Italy, September 4-8, 2006, copyright by EURASIP

If 1 + d˜ and 1 + δ are approximated by 1, we can solve for Simulated estimator error, N=10000 Simulated estimator error, N=100

derr ≈ d − d0 and get

q x 10
−6
x 10
−5

1 1 ˜
derr ≈ − ′ + ω02 − 2(δ ′ )2 δ − d. (20) 3 3
δ ω0 δ ′ 2.5
2.5

This expression can be used to predict the final estimation er- 2 2

Simulated error
Simulated error
ror. The first part becomes small when ω02 is large compared 1.5 1.5

to 2(δ ′ )2 δ and the error will then be dominated by d.˜ For 1 1

small ω0 , if δ or δ ′ are not zero, the error deviates more and 0.5
more from d.˜ If d ′ = 0 the error becomes derr = d.˜
0.5

0
0
0.4
0.4
5. FILTER DESIGN 0.3
0 0.3 0
0.1 0.1
As was seen in (17) the main error source is the delay er- 0.2 0.2
0.2
0.2

ror d˜ which directly affects the estimate. However, when the

0.3 0.1 0.3
0.1
0.4 0.4
derivative of the magnitude error, δ ′ , is nonzero the magni- d 0 0.5 ω0/π d 0 0.5 ω /π
0

tude error δ and the frequency ω0 will come into effect. To

find an FD filter which is optimized for small estimation er- Figure 5: The simulated estimation error for N = 10000 and
rors, the expression in (20) can be used as cost function in N = 100 samples. 0.1 < ω0 < π /2.
an optimization problem. The optimization problem can for
example be formulated as a minimax problem

minimize ε subject to |derr (ω , d)| < ε (21)

which is then optimized over a range of angular frequen-

cies ωk and delays dk . The fminimax-routine in MATLAB
efficiently implements a sequential quadratic programming
method that is capable of solving the minimax constraint
problem in (21). However, the solution is not guaranteed to
be the global minimum as the problem is nonlinear.
The derivative of δ , which is needed to calculate the ex-
pected estimation error, can be found analytically by noting
that the magnitude of the Farrow-based FD filter can be writ-
ten as
⌊L/2⌋
|H(d, e jω T )| = ∑ d 2k G(2k)R (ω T )+
k=0
⌊(L+1)/2⌋
+j ∑ d 2k−1 G(2k−1)R (ω T ) Figure 6: The expected error. 0.1π < ω0 < π /2.
k=1
= |A(d) + jB(d)| = 1 + δ (22)

where G(2k)R and G(2k−1)R are the zero phase frequency re- 6. PERFORMANCE OF THE ESTIMATION ERROR
sponses for the even-order linear-phase FIR filters with im- PREDICTION
pulse responses gk (n). To verify the error analysis, a number of simulations were
Since the derivative of the magnitude is equal to the performed. An FD filter with L = 7 subfilters was optimized
derivative of the magnitude error we can calculate δ ′ as for minimal error up to the frequency ω0 = π /2. The result-
ing even-order subfilters had orders Mk ranging from 8 to 18.
A(d)A′ (d) + B(d)B′ (d)
δ′ = 2 p (23) The maximum magnitude error was δmax = 1.83 · 10−4 , the
A(d)2 + B(d)2 maximum delay error was d˜max = 1.168 · 10−5 and the max-
imum derivative of the delay error was d˜max
′ = 1.01 · 10−4 .
where A′ (d) and B′ (d) can be calculated as In Fig. 5 the simulated estimation error for N = 10000
⌊L/2⌋
and N = 100 samples can be seen. As predicted by (20) the
estimation error increases for decreasing frequencies when
A′ (d) = ∑ 2kd 2k−1 G(2k)R (ω T ) (24)
there is a magnitude error. The reason for the lower perfor-
k=1
mance for N = 100 samples is that N must be larger for the
and variance of FN′ to become small, especially when δ and δ ′ are
relatively large, which they tend to be for small ω0 and large
⌊(L+1)/2⌋ d. In Fig. 6 the expected error derr can be seen. Compared to
B′ (d) = ∑ (2k − 1)d 2k−2 G(2k−1)R (ω T ). (25) the simulated error in Fig. 5 a smaller frequency range has
k=1 been used to enhance the details.
 14th European Signal Processing Conference (EUSIPCO 2006), Florence, Italy, September 4-8, 2006, copyright by EURASIP

Difference between expected and simulated error, N=10000 because for small N noise and the variance of FN′ will dom-
inate, while for large N the filter approximation errors will
dominate.
−8
x 10
8 7. CONCLUSIONS
6 We have presented a novel method to estimate the delay er-
ror between two sets of samples using a FD filter. The idea
Difference

4 is to use an iterative, Newton-Raphson-based, estimator to

minimize the mean-squared-difference between a reference
2
signal and a signal with an unknown delay.
0 The effects of a nonideal FD filter with magnitude and
0.4 0.1 delay errors have been studied theoretically and in simula-
0.3 0.2 tions. An expression of the expected estimation error derr
0.2 0.3
was derived, which can be used to optimize the FD filter in
0.1 0.4
the estimator for a minumum estimation error. The expected
0 0.5
d
ω0/π estimation error derr is a bias, or in other words, the best one
we can achieve without noise, i.e when N → ∞.
Figure 7: The difference between the expected estimation er- In the analysis we have assumed a single sinusoid, but the
ror and the simulated error, N = 10000. 0.1π < ω0 < π /2. time delay estimator can be used for more general bandlim-
ited signals. However, in this case, the estimation error will
differ from the case where a sinusoidal is used. The analysis
done in this paper is still useful since it gives an insight into
errors caused by the nonideal interpolation performed by the
x 10
−6 fractional-delay filters. Additionally, in some applications
7
where the training sequence may be chosen freely we can
actually choose a sinusoid and achieve the limits computed
6
in this paper.
5 The time-delay estimation method can easily be extended
Difference

4 to more sets of samples, using one set as the reference, which

3
could be useful e.g. in the calibration of time-interleaved
analog-to-digital-converters (TIADCs).
2

1
REFERENCES
0
0.4
0.3
[1] S. R. Dooley and A. K. Nandi, “On explicit time delay
0.1
0.2 0.2 estimation using the Farrow structure,” Signal Process-
0.3
0.1
0.4 ing, vol. 72, pp. 53–57, 1999.
0 0.5
d
ω0/π [2] C. W. Farrow, “A continously variable digital delay el-
ement,” in Proc. IEEE Int. Symp. Circuits & Syst., June
Figure 8: The difference between the expected estimation er- 7–9 1988, vol. 3, pp. 2641–2645.
ror and the simulated error, N = 100. 0.1π < ω0 < π /2. [3] J. O. Smith and B. Friedlander, “Adaptive interpolated
time-delay estimation,” IEEE Trans. Aerospace, Elec-
tronic Syst., vol. AES-21, no. 2, pp. 180–199, Mar 1985.
[4] D. L. Maskell and G. S. Woods, “The discrete-time
When derr was derived we omitted FN′ , assuming that N quadrature subsample estimation of delay,” IEEE Trans.
is infinite or at least large enough. If the batch size N is Instrum. Meas., vol. 51, no. 1, pp. 133–137, Feb 2002.
decreased the expected estimation error will increase, which
can be seen in Fig. 7 and 8. As expected, the difference [5] J. Vesma and T. Saramäki, “Optimization and efficient
between the estimated error and the simulated error is in- implementation of FIR filters with adjustable fractional
creased when N becomes smaller so that FN′ no longer can delay,” in Proc. IEEE Int. Symp. Circuits Syst., June 9–
be approximated to be zero. For N = 100 samples the esti- 12 1997, vol. IV, pp. 2256–2259.
mator still performs quite well, see Fig. 5, but the expected [6] H. Johansson and P. Löwenborg, “On the design of ad-
error derr no longer works as a good estimation of the actual justable fractional delay FIR filters,” IEEE Trans. on Cir-
error since the expected error is almost as large as the actual cuits and Syst. II: Analog and Digital Signal Processing,
difference between the simulated and expected error. This is vol. 50, no. 4, pp. 164–169, Apr 2003.

View publication stats