Signal Processing - November 2017 - 166

Appendix
The derivation of (4) proceeds in three steps: 1) derive a
polynomial expression approximating Figure S1's continuous p (t ) sinusoid in terms of known p [n] samples; 2) set
that expression's time derivative equal to zero; and 3)
replace t with f and solve for f.
Our derivation begins by assuming the largest of our known
p [n] samples is located at time t = 0. Approximating
Figure S1's p (t ) sinusoid with a second-order Taylor series
expression in the vicinity of t = 0, we begin by writing
	

p (t ) . p (0) + pl (0) t + 1 p m (0) t 2, (S1)
2

where pl(0) and p m(0) represent the first and second derivatives of p (t ) at time p (t ) = 0. Given the (S1) polynomial,
we next approximate the unknown pl(0) and p m(0) coefficients by using the central difference formula. Doing so we
write the first-order derivative pl(0) as
pl(0) .

	

p (h) - p (- h)
.(S2)
2h

Having an approximation of pl(0), we next approximate
the second-order derivative p m(0) using the first-order derivatives centered at the hypothetical p (- h/2) and p (h/2)
samples in Figure S1. Those first-order derivatives are
p (0) - p (-h)
and
h
p (h) - p (0)
pl(h/2) .
.
h

pl(-h/2) .

	

[0, 1) range. The computation of parameters offset and fr are shown as the firststage processing in Figure 3.
The neat trick of our proposed algorithm is that neither the p [n] sequence
nor the continuous p (t) signal need to
be computed. Our first-stage processing
produces a rough estimate of the angle
166

p (h)

−h
-h
2

p (−h)

0

f h
2

h

t

FIGURE S1. The sinusoidal p (t ) signal and the desired time value f.
Assuming the time between our known p [n] samples is
unity sets h = 1 and recalling that, for our p [n] samples,
p [-1] =-p [1], we can rewrite (S2) and (S3) as
pl(0) .

	

p (1) - p (- 1)
= p (1)(S4)
2

p (1) - 2p (0) + p (-1)
12
=-2p (0).
(S5)

p m(0) .

	

Substituting (S4) and (S5) as coefficients in (S1), our
desired approximation of p (t ) is
p (t ) . p (0) + p (1) t - p (0) t 2 .(S6)

	

That completes the first step of our derivation. As the second step
of our derivation we take the derivative of p (t ) to produce
pl (t ) . p (1) - 2p (0) t.(S7)

Setting (S7)'s pl(t ) = 0 gives us an approximation of the
time location of the maximum value of Figure S1's p (t ) signal. Doing so and defining that estimated time value as fr
we write
	

p m(0) .

p (h/2)

p (−h/2)

	

Given pl(- h/2) and pl(- h/2), we write our desired second-order derivative p m(0) as
pl(h/2) - pl(- h/2)
h
p (h) - p (0) p (0) - p (- h)
h
h

.
h
p (h) - 2p (0) + p (- h)
=
.
h2

p (t )
p (0)

0 = p (1) - 2p (0) fr .(S8)

Finally, solving (S8) for our desired expression for f in
terms of known p [n] sample values we arrive at the final
form of (4) as
(S3)

f . fr /

	

of c = b + ja based upon some simple
logic and simple arithmetic using values
a and b. The offset can be computed
using the signs of a + b and a - b, as
shown in Table 1. The determination of
the two samples needed for the computation of fr can also be performed using
the signs of a + b and a - b. It should
IEEE SIGNAL PROCESSING MAGAZINE

|

November 2017

|

-pl(0)
p (1)
=
.(S9)
p m(0)
2p (0)

be pointed out that the variable used in
the denominator in both expressions for
fr is always the largest absolute value
between a and b.
Note that, in [5], Shima used an
approximation for the one-variable
atan(x) (derived from a first-order
Lagrange polynomial interpolation of
Table of Contents for the Digital Edition of Signal Processing - November 2017
