846 IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-34, NO.

4, AUGUST 1986

Sampling Theorem for Polynomial Interpolation

FRANCESCO PALMIER1

Abstract-The classic polynomial interpolation approach is used to sidered in this paper includes superposition of exponen-
derive a sampling theorem for the class of signals that are the response tially damped sinusoids that, in general, cannot be inter-
ofsystemsdescribed by differentialequationswithconstant coeffi-
cients. In particular, the polynomial that interpolates the signal be- polated at theNyquist rate using the cardinal interpolation
tween the sampling points for increasing order will give increasing ac- formula.Our attempt is to look at the polynomial ap-
curacy for stable one-sided sequences, if the sampling rate is at least proach that preceded and resulted in the discovery of the
six times the highest pole frequency (Bolgiano sampling rate [7]).The cardinal interpolation formula. Polynomials, in fact, rep-
convergence is ensured also for nonstable poles that lie in a certain
resent a wider class of functions [4] and can be useful for
region of the complex plane.
If, instead,we use the symmetrical polynomial approach, it i s enough the reconstruction, from evenly spaced samples, of sig-
to sample at a rate that is just two times the highest pole frequency nals that are the response of linear systems described by
(Nyquist sampling rate), with some constraintson the real part of the differential equations with constant coefficients. Basi-
poles. cally, we will address the problem of selecting the sam-
A bound for the erroris derived for both cases and a comparison to pling rate according to the position of the poles in the
the Shannon-Whittaker sampling theoremis presented.
complex plane. In many practical cases, in fact, the use
of higher order polynomials may not improve the quality
I. INTRODUCTION of the reconstruction.
Thus, we will derive, using a very simple test, a con-
T HE Shannon-Whittaker sampling theorem has been
an essential contribution to the development of the
techniques for digital signal processing even if often the
dition for the convergence of the interpolation process.
Furthermore, since in many situations the Lagrange re-
basic hypothesis of dealing with band-limited functions is mainder gives an unreasonably large estimation for the
violated. Sharp edge signals, finite duration signals, and error, we will use the condition for uniform convergence
decaying exponentials have spectra that are far from band- to derive some tighter bounds.
limited.
The origin of the cardinal interpolation formula goes 11. POLYNOMIAL INTERPOLATION
back to the work of Whittaker [21], in 1915, that posed Given a real functionf(t), t E R , and (n + 1) values at
- -
the problem of finding a ‘‘cotabular’’ function that passes the points (to, tl, * , t n > ,the unique polynomial p,(t)
through a set of points. A few years later, his son J. M. of degree n that interpolates the given values can be ex-
Whittaker,in about 1920, reevaluated the problem [4]; pressed in the Lagrange form [l]
and, only in 1449, Shannon, in his famous paper [24],
pointed out that band-limited functions are exactly rep-
resented by the cardinal interpolation expansion.
We also have to mention that the same result was pub- where
lished in the Russian communication literature in 1933 by
Kotel’nikov [23]. W n ( t ) = (t - to)(t - t l ) * * (t - tn), (2)
Ever since, many investigations have been carried out and
to evaluate the performances of the cardinal interpolation
formula applied to real data. In particular, the purpose has w;(tk) = ( t k - t o ) * ( t k - t k - 1)

been to predict the degree of accuracy produced by its

* ( t k - t k + ~ )* . * ( t k - t n ) , (3)
truncated version applied to signals that have a band-lim-
ited Fourier transform [32]-[36] or band-limited power or in the Newton form [l]
spectral densities [23]. n
Unfortunately, when we are faced with the problem of
Pn(t> = kC f ~ t o ,. * tk>(t - to)
3
reconstruction of signals that are not band-limited, it is =O
not always clear if the cardinal interpolation formula can (t - tl) ‘ * (t - t k - l ) , (4)
give a reasonable approximation.
The classof nonband-limited functions that will be con- wheref(to, -- * , tk) denotes the kth divided difference of

Manuscript received July 10, 1985; revised December 13, 1985.

f on the points (to, - - , t k } . The natural questions that
The author is with the Electrical Engineering Department, University of
arise are: does the polynomial p,(t) accurately represent
Delaware, Newark, DE 19716. the original function between the sampling points; what
IEEE Log Number 8608129. is the distribution of the knots that produces the best ap-

0096-3518/86/0800-0846$01.OO 0 1986 IEEE

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 PALMIERI: SAMPLING THEOREM
INTERPOLATION
FOR POLYNOMIAL 847

proximation; and more specifically, does increasing the Basically, an integral function is a function that has only
number of samples lead to a better approximation between one pole at infinity.
the sampling points. These questions face the so-called It can be shown [ l , p. 1231 that for integral functions,
problem of “convergence of the interpolation process.” the sequence of interpolating polynomials corresponding
To explain this, let us referto the triangular matrix to a matrix of the kind (5) converges uniformly tof(t), if
t belongs to a finite interval [a, b].
t:
However, the problem we are usually interested in for
t p tp digital processing of signals deals with uniformly spaced
t p t{a t p samples that span an infinite range, and the above result
cannot be applied.
... Since, as we said before, to ensure convergence of the
t$4 $1 . . . $0 (5) interpolation processwe need somehypothesison the
structure of the signal, wewill restrict our attention to the
where the values t i j ) belong to a certain subset Z of the class of functions very common in the engineering liter-,
real axis that could be finite or infinite. If, at each row, ature
we associate a polynomial p,(t), we say that the interpo-
lation process is convergent if N

lim p,(t) = f(t) Vt E I. (6)

f(t) = C altkleSit, t E R, ( 10)
1=1
,403

To explain the above statement, letus consider, as an ex- where we will refer to the parameterssl = ( x 1 + j w , as the
ample, the finite interval [a, b]. At first glance, it would “poles” of the Laplace transform of f(t). We will use
seem that, for any function, sampling the interval [a, b] such a denomination even if it is, in general, arbitrary
more,densely would lead to a better approximation. This because it requires the definition of the Laplace transform
is, in general, not true! In fact we can state thefollowing of f(t) and its domain of convergence. However, if we
theorem due to Farber [l]. assume that f(t) = 0 for t < to, (10) corresponds to the
For any defined matrix of points ( 5 ) , there is a Con- impulse response of a system that has a rational transfer
tinuous function f(t) such that the Lagrange inter- function with N poles s j , each with multiplicity kj - 1 .
polation polynomials, constructed foritonthese Such a system is usually referred to as a “natural sys-
points, do not converge uniformly on [a, b] to f(t). tem.” I

We will consider the interpolation for the single term

A classic example shown by Bernstein is the interpo- with s = Q! + j w
lation of the function f(t) = (tl for t E [I, 11. The poly-
nomial p&) for evenly spaced samples for increasing n f ( t ) = atk e’‘, (1 1)
does not converge to f ( t ) at any point other than - 1, 0,
since the results for the sum of them can be obtained by
1 . As another example, Runge [2] considered the function
- -
superposition.
1 The function (1 1), in general, is a complex function
in
f(0 = 1 + 100t2’ t E [-1, 11. (7) the real variable t. If we extendits definition inallthe
complex plane, we can say that f ( t ) is an integral func-
Theinterpolationprocessfor any knot distribution,for tion, and can be expanded in Taylor series around any
increasing n, does not ‘Onverge tof(t)* This could be ex- finite point. In particular, it will have a polynomial rep-
plained by expanding in series and observing that resentation. Now if its polynomial representation con-
03
verges V t for increasing order, it will converge to the cor-
1 +
1
loot2
= c (-100t2)k,
k=O
(8) rect valuef(t) V t .

only for It( < &. SEQUENCES

111. ONE-SIDED
The problem of the convergence has been largely stud-
ied in the literature, and since many of the results are re- Let US consider the following triangular matrix corre-
ferred to integral functions, it isprobably useful to remind sponding to evenly spaced samples:
the reader of the following simple definition. 0
An integral (or entire) function f(z) is an analytic O T
function which is regular for all finite z. In particu-
lar, it can be represented in the form of a Taylor 0 T 2T
series ...
f(zj = a. + al(z - z0) + az(z - z ~ ) ~ 0 T 2T * - nT,
(12)
+ - + a,(z - zo)” + - ,
* * (9) the corresponding interpolating polynomial in the
La-
which converges
allfor values of z. grange form is

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 848 IEEE TRANSACTIONS ON ACOUSTICS,
SPEECH, AND SIGNAL
PROCESSING, VOL. ASSP-34, NO. 4, AUGUST 1986

Sufficient condition for uniform convergence of (20) can

be found using the ratio criterion

lim gk
k-+m
___
1 I
-!-l(t)

gk(t) <"
or in the Newton form where gk(t) is the kth term of (20). In particular,
n
p,(t) = A")j
t(t - T) * . (t - (k - 1) T )
3 (14)
k=O k! Tk
where A(k)fdenotes the kth finite descending difference on and for any finite t,
the set of points (0, T, ---
, nT) . (The series { p n ( t ))
written in the form (14) for T = 1 was discovered in 1670
by James Gregory, and isusually referred as the Gregory-
Newton series [ 131.) Therefore, the condition to satisfy is
Let us restrict our attention to the case of the pole with
multiplicity one. lesT - 11 < 1,
f(t) = ae". (15) that can be rewritten as
The finite descending difference for this function on the -a
samples ( f ( O ) , f V ) , * * ,f(nT)) can be written as
aT < In (2cos UT),
2
-
2
< UT < -.a
The condition (24) not only ensures uniform conver-
gence, but also a series with terms decreasing in absolute
value. In fact (22) for t E [0, kT] can be bounded

= a(eST- ilk. (16)

It is interesting to derive also the divided difference using and if lesT - 11 < 1,
the Peano theorem [5] Igk+l(t)l < Igk(t)l vt E [o, kT1
f(t0, 11, - * + ,t"> Vk = 0, 1, 2, * * -. (27)
-
- E
1
izo zk
w x ; to, tl, * * - , t k ) f ( k ) ( X ) h, (17) The condition (24) was found by Bolgiano [6]-[12] by a
different means and is depicted graphically in Fig. 1 as
where f k ) ( x )is the kth derivative of f ( x ) and M ( x ; to, t l , the area in the plane esTincluded in the circle with radius
- * * , tk) is the kth-order B-spline defined as the divided
1 and center (0, 1). Sufficient condition on the sampling
rate for stable poles will reduce to [UTI < a/3, that is,
difference on the points {to, tl, * * , tk)of thePeano
kernel 1
- > 6f Bolgianosamplingrate (28)
M(x;t) k(x - t):. (18) T
Therefore, in our case where f = d 2 a . Notice that the Bolgiano sampling rate
depends only on the natural frequency of the poles.
f ( t 0 , tl, -- * tk) = -

- esxdx
ST
k ! a o M(x; 0, T, * * , kT)sk Practical use of the condition stated above requires
some measureof the error whenthe interpolating series is
truncated to a finite number of terms. An expression for
the residual in polynomial interpolation is the well-known
ask Lagrange remainder
-
- - 6: [M(t; 0, T,
k!
. * , kT)]

f ( t ) - P A 4 = f'"'''(f)
(n l)! t(t - T )
+ e . (t - nT), (29)

where 4 is some point in the interval [O, nT]. In our case

The Newton interpolation formula will be (29) would become
1 8E
.
f ( t ) - P&> = a
(n + l ) ! t(t - T ) . (t - nT),
* (30)

t(t - T) * (t - (k - 1) T). (20) that for t E [0, nT] can be bounded

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
AMPLING PALMIERI: THEOREM
INTERPOLATION
FOR POLYNOMIAL 849

The condition (24) leads to a boundfor the error that is

a geometric succession approaching zero for increasing n
as (esT- 11".
If instead of considering the simple pole ( 1 3 , we refer
to a pole with multiplicity q > 0, the condition (21) be-
comes

and applying 1'Hopital rule backward for t E [0, kT] we

get

I F

Therefore, (24) is still a sufficient condition for uniform

convergence for the case of multiple poles. However, it
-?r I is easy to show that the inequality (27) will be satisfied
only for q = 1, and consequently, the bound (35) will be
(b)
Fig. 1. Region of convergence for one-sided polynomial interpolation (a)
valid only for simple and double poles.
in the plane esT,(b) in the plane (aT,UT). The aboveresults could be summarized in the following
theorem.
Theorem I .- Given a signalof the form (lo), a sufficient
condition for uniform convergence of the interpolation
process on the points (12) is that all the poles s = CY jo +
satisfy the condition (24) that corresponds to the region
This expression, unfortunately in 'our case, gives an un- inside a circle withradius 1 centered at (1, 0) in the com-
reasonably large estimation of the error, and it is probablyplex plane esT. This condition will ensure, for poles of
of little practical interest. Nevertheless, due to the con- multiplicity less or equal to two, convergence with a se-
dition (24), the inequality (26) allowsustoderivea ries that has terms decreasing in absolute value corre-
tighter bound for the error in the following simple way. sponding to an error bounded by the monotonically de-
From (24) and (26), creasing sequence (35).
The above theorem can be very useful if we aredealing
Igk+l(t)l < lesT - 1llgk$!)l V t E [o, nT] with segments of a signal that can be reconstructed be-
v k = 0, 1, e -, * (32) tween the sampling points via the unique polynomialp,(t),
namely, in therange [0, n T ] . The reco.nstructed signal
Igk+l(t)l < lesT - Ilkgo = lesT - IIa, (33) could be,forinstance, used to comparethe alias-free
0) short-time Fourier transform [101.
At this point, it is worthwhile to make further com-
ments onthe range for thevariable t in which we consider
therefore, the convergence of the polynomial interpolation process.
0)
Since the results of theorem 1 are valid V t E [0, k T ] ,
they will hold in particular for t E [0, T I . This means that
the samplingcriterion given by theorem 1 can also beused
n \ if we consider the polynomial of degree n on the points
{iT,(i 1)T, *+ - +
, (i n ) T ) andwe are interested in
the reconstruction of the signal in the interval [iT,(i +
and 1)T], V i = 0 , 1 , - -
. Thisistheapproach usually re-
ferred as the one-sided piecewise polynomial interpola-
tion, since in each interval between sampling points the
signal will be approximated by a different polynomial.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 850 IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-34, NO. 4, AUGUST 1986

IV. SYMMETRICAL PIECEWISE POLYNOMIAL

INTERPOLATION
If, instead of considering polynomials that span the
range only at theright of the interval inwhich we consider
the reconstruction, we can place the polynomial interpo-
lant symmetrically with respect to, say, [iT, (i 1)T], +
degree n will spantherange [(i - (n - 1)/2)T, (i +
(n +1)/2)T].We will refer to this approach as central
piecewice polynomial interpolation. Let us consider the
triangular matrix
0
O T
= /e-Ts(esT -
-T 0 T

' )(t + 2k(2k

-T 0 T 2T (k l)T)(t kT)
- -
-2T -T 0 T 2T 1)T2 -
...
(42)
The above expression can be bounded for t E [O, TI
The Newton polynomial,
that in this form is usually re- k
ferred to as the Gauss-Newton series, can be written ex-
plicitly in terms of the finite central differences

(43)
therefore, a sufficient condition for the two series to be
t(t - T )
+ 2! T2
A("( - T, 0 , T ) convergent and with terms decreasing in absolute value
for t E [0, T ] is
+ ...
(44)
+ t(t - T ) ( t + T) . (t - (k - l ) T ) ( t
* + (k - l ) T )
(2k - 1 ) ! T2k-1 or
. A(2k- 1)f ( - ( k - 1)T, . * * , 0 , * * . , k T )
+ t(t - T)(t + T) * *

2kT2k
+
(t (k - l ) T ) ( t - k T )
that can be rewritten as
I I:
sinh - < 1, (45)

. Ac2ky(-kT . . . , 0, . . . , kT) cosh a T < COS UT + 2, (46)

+ t(t - T)(t +(2kT )+ 1 ) !(tTZk+l

- kT)(t + kT)
*
aT < cosh-' (cos U T + 2) for aT >0
(47)
a T > cosh-' (cos U T + 2) for aT < 0 ,
A(2k+11f(-kT,. ,0, * * , (k + 1)T),(39)
* *
or
where n = 2k + 1. If we find convergence conditions for
the central interval [0, T I , the result will apply for shift
(YT < In [COS U T + 2 + J(COS + 212 - ij UT

invariance to any interval centered with respect to the for aT >0

knots considered.
The series ( p , ( t ) } can be seen as the sum of two series CXT> In [COS UT + 2 - J(COS + 212 - 13 UT

formed, respectively, with the even and odd terms. Ap- for a T < 0. (48)
plying the ratio test to both of them, the condition for
uniform convergence will be The corresponding regions of the complex planes ( a T ,
U T ) and esTfor which the condition (48) is satisfied, are
depicted in Figs. 2 and 3.
It can be proved, in a way similar to the one used for

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 PALMIERI:SAMPLINGTHEOREM FOR POLYNOMIALINTERPOLATION 851

z tWT 10-10
0 . 0 0 2 . 0 0 Y . 0 0 6 . 0 0 8.00 10.0 12.0 111.0 16.0 18.0 2 0 . 0 2
(a)
B=1.170 alpha=-0.50 om=1.57 Tx1.00 C0,nTI

'i 10-2

f 10"
10-5
I 0 - l
10-6

2.00
(b)
Fig. 2. Region of convergence for symmetrical polynomial interpolation
(a) in the plane esT, (b) in the plane (aT,UT).

I
::::
10-9 i
10 -100 . 0 0 2 . 0 0 9 . 0 0 6 . 0 0 8 00 10.0 12.0

(b)
n
lY.O 18.0 2 016.0
. 0 2i

Fig. 4. Maximum error for the nonconvergent case for the one-sided poly-
nomial interpolation,with (a,w ) = (-112, d 2 ) . (a) The reconstruction
error js computed in the interval IO, T I . (b) The reconstruction error is
computed in [0, n T ] .

wT/?r
Fig. 3. Particular region of convergence for symmetrical polynomial in-
terpolation in the plane (aT,wT).

one-sided sequences, that the uniform convergence is ver-

ified in the same region for multiple poles.
Also, a bound for the error can be derived in a similar and for t E [ O , TI the error can be bounded in the follow-
way. We will restrict our attention to the case of simple ing way:
poles. The bound will be the sum of the two bounds cor-
responding, respectively, to the series with even terms and
the series with odd terms. Namely, defining A as

(49)

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 852 TRANSACTIONS
IEEE ON ACOUSTICS,
SPEECH, AND SIGNAL
PROCESSING, VOL. ASSP-34,
NO, 4, AUGUST 1986

B=0.714 alpha=-0.50 om=0.79 T=1.00 CO,T3 B=O 765 alpha= 0 om=0.79 T=l.OO CO,Tl

I:
101
I ::: ,
100

10-1 -
#? j
10-2

10-3

IO+
z 10-~
10-6

10-9 4
10-100.00 2.00 't.00 6 00 8 0 0 10.0 12.0 l't.O 16 0 18.0 20.0 2 2 . 0

B=0.714 alpha=-0 50
(a)
om=0.79 T = 1 . 0 0 C0,nTl
0.00 2 00 't.00 6 00 8.00 10.0 12.0 19.0 16.0 18.0 20.0 2 2 . 0

B=0.765 alpha= 0 .
(a)
om-0.79T=1.00 [O,nTJ

IO'
100 2
10-1 -
'j 10-2
k 10-3
at
10-9
10-5
10-6

10-10
0.00 2.00 't.00 6 00 8 00 10.0 12.0 111.0 16.0 18.0 20.0 22.(I 0.00 2 00 9.00 6.00 8.00 10.0 12.0 l't.0 16.0 16.0 20.0 2

n n

(b)
Fig. 5 . Maximum error for the convergent case for the one-sided polyno- Fig. 6. Maximum error for the convergent case for the one-sided polyno-
mialinterpolationwith (a,w ) = (-1/2, xi4). (a) Thereconstruction mial interpolation with (a,w ) = (0,x/4). (a) The reconstruction error
error is computed in the interval [0, T I . (b) The reconstruction error is is computed in the interval [ 0 , TI. (b) The reconstruction error is com-
computed in [0, n T ] . 0-bound and 0-error. puted in [0, a T ] .

Theorem 2: Given a signal of the form (lo), a sufficient

condition for uniform convergence of the symmetrical
polynomial approximation is that the poles and the sam-
But gdt> = a and gl(t) can be bounded for t E [0, TI pling rate are such that the condition (45) is satisfied. The
central piecewise approximation forsimple poles will
converge with terms decreasing in absolute value and a
bound for the error is given by (55).
< If(T) -f(O)I < aleST- 11. (54) Theorems 1 and 2 will allow us to select, according to
the specific application, the sampling rate for polynomial
Finally, since n = 2k +
1 we have
reconstruction if the signal is assumed to have the form
A'" + 1)/2 (10). In particular, for piecewise approximation, as dis-
E , < a(1 + lesT - 11) ~
(55) cussed by Schoemberg [30] and Hopcroft and Steiglitz
1-A'
[ 151, it is possible to express the interpolated signal as the
Notice that the sufficient condition still holds, even if the convolution
bound (55) would not be the same, in general, if we con-
+m
sider any finite value of t instead of the only interval [0,
TI. This implies that (45) is a sufficient condition for uni-
form convergence of the symmetrical polynomial inter-
polation also for t E [ -iT, (i 1)Tl. + where a,(t) is a pulse function such that a,(O) = 1 and
The above results can be summarized in the following a,(kT) = 0 vk = * -2, -1,1, -
2, * * , selected
theorem. according to our approximation criterion.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 PALMIERI: SAMPLING THEOREM FOR POLYNOMIALINTERPOLATION

'
853

B=0.71+alpha=-0.50om=0.79 T=l 00 CO,Tl

10' 1 I

b 10-2

l"
10-3
10'
X
;10"

:I,
10-9
10-10

10'
, I , ,
0 . O O 2 . 0 0 Y . 0 06 . 0 0
(a)
B=0.87+ alpha= 0.20 om=0.79 T=l 00
, I , , ,
8.OD 10.0 12.0 19.016.0
I .

CO,nTl
I , ,
18.0 20.0 2 2 . 0
"

10"
10 -10
~ ~

0 0 02 . 0 0

B10.765alpha=
, ~ , 1 , , , 1

Lt.00 6 . 0 0 8 . 0 0 10.012.0

0.
(a)
,

om=0.79T=1.00
1 , 1 . , ' ,
~

19.0 16.018.020.02

C0,Tl
~ ~ , 1

10-6 -
- :L-
110-8
10-10
0-~

10-10Lt.00 6 . 0 0 8 . 0 0 1 0 . 0 12.0 lr.0 1 6 . 0

10" 0 . 0 0 2.00 18.0 , 20.0 22.0
0 . 0 0 2.00 9 . 0 06 . 0 0 8 . 0 0 10.0
12.0
lY.O 16.0 18.0 2 0 . 0 2i
n n
(b) (b)
Fig. 7. Maximum error for the convergentcase for the one-sided polyno- Fig. 8. Maximumerror for double poles in one-sided polynomialinter-
mial interpolation with (a,o) = (1/5, d 4 ) . (a) The reconstruction error polation with (a,a) = (a) (-112, d4); (b) (0, d4).
is computed in the interval [O, T I . (b) The reconstmction error is com-
puted in [O, FIT].
where
It is well known [27], [31] that the above formula (56)
can be used to design a digital filter that simulates H ( j w )
for input of the form (10) as
~ ( z )= Z[F-'[H(jw) A ~ W J ) I L = ~ ~ , and
where A,( j w ) is the Fouriertransform of the interpolation H(t) = lim Hn(t).
n-+m
pulse a,(t) that can be easily computed in the way sug-
gested in [l5]. Defining
V . COMPARISON TO THE SHANNON-WHITTAKER
SAMPLING THEOREM
As already pointed out in the Introduction, thecardinal
series was originally derived from the Lagrange interpo- (60) can be written as
lation formula. In fact, we can rewrite (1) on the set of
-
points {--tn, * , to, - * , tn>inthe following way [4]:

(57) If the knots are equally spaced, that is, tk = kT,

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 854 TRANSACTIONS
IEEE ACOUSTICS,
ON SPEECH,
SIGNAL
ANDPROCESSING, VOL. ASSP-34, NO. 4, AUGUST 1986

A=1.055 alpha=-0 60 om=2 75 T=1.00 A=0.962 alpha= 0 . om=2.75T=1.00

10'
100
10-1

:
LA
10-100.00 2 00 11 00 6.00 8 00 10 0 12.0 111.0

A = l 355
(a)
alpha=-1.90 om=0.79
0 18.0
16

T=1.00
22.0
20.0 10-100 00 2.00 11 00 6 00 8.00

A=O 500
(a)
alpha= 0 .
10.0111 0 16.0

om=l 57
12.0

T=l . O O
18.0
0 22 0 20

:::
101

I
:f
10 -100 . 0 0 2.00 11 00 6 00 8.00 1 0 . 0

(b)
n
12 0 111.0 16.0 18 0 20.0 22.0

Fig. 9. Maximumerror for nonconvergentpiecewisesymmetricalpoly-

nomial interpolation with ( a , o) = (a) (-315, 7n/8);(b) (-19110, n14).
10-100.00 2 P O 11.00 6 00 8.00 1 0 . 0 12.0 111.0 16.0 18.0 20.0 2i

(b)
n

Fig. 10. Maximumerrorforconvergentpiecewisesymmetricalpolyno-

mial interpolation with ( a , w ) = (a) (0,7 ~ 1 8 )(b)
; ( 0 , n12).
0

t we can state the following theorem reported by Whittaker

a, sin x -
H(t) = t
i= 1 (1 - &) = T-, x
T
c62)
in [4].
If the cardinal series

s)
and R
+m sin - (t - nT)
01 T
n!n! P(t) = f(nT>
+(n, k) = (1 - = . (63) x
k) ! (n + k) !
n= -01
i=n+l (n - - (t - nT)
T
If a series
R
n

converges, we say that the series h k is convergent by the

method of De la Valee' Poussin ( l l n 2 ) [4, p. 641.
Now since (64) can be rewritten as

is convergent,the Gauss-Newton series (64) con-

~ 2 n ( t )= 7+ Wn7 k)(-1lk
verges to the same sum.If the Gauss-Newton series
(64) is convergent, the cardinal interpolation series
(65) is convergent by the method of De la Vallee'
(@) Poussin to thesamesum.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 PALMIERI:SAMPLING THEOREM FOR POLYNOMIALINTERPOLATION 855

A=O 187 alpha= 0 .40 om=0.79 T=l . O O A=O 500 alpha= O . om=l.57 T=l . O O
103
102
10' 101
100 100
10-1
tj 10-2
10-3
b lo*
X
2 -

:
10+
10-6 -
10-~
.
10-8

l o
0
- ~
bo 2 . 0 09 . 0 0
~ ~ ~ l ~ l . l
6 00 8 00 1 0 . 0 1 2 . 01 9 . 01 6 . 0
~ l ' l .
18.0 2 0 . 02 2 . 0
l ' l ' i:::
I '
1 0 - ~0.00
10-10 I2 00~ 9 I00 ~
6 . 0 0 l 8 00 10 0 12.0 19.0 0 12 60 . 00 2
2.
18

(a) 4
A=O ,315 alpha=-0.80 om=O .79 T=l . O O A=O 500 alpha= 0 . om=l.57 T=l . D O
lo3
102
101 101
100

lo-'
tj 10-2
10-3
b 10''
X
10-5
10-6
10-~

10-~
l o
1 - ~ ~ ~ ~ ~ ~ l ~
0.00 2 . 0 0 9.00 6 . 0 0 8 00 1 0 . 0 1 2 . 0 19.0 1 6 . 0 18.0 2 0 . 02 2 . 0
l ~ l * l ~ l ~ l
0.00 2 . 0 09 . 0 0
~ l ' l ' l ' l
6 . 0 0 8 . 0 01 0 . 01 2 . 01 9 . 01 6 . 01 8 . 02 0 . 0 2:
n n
(b) (b)
Fig. 11. Maximumerrorforconvergent piecewise symmetricalpolyno- Fig. 12. Maximumerrorforconvergent piecewise symmetricalpolyno-
mial interpolation with ( a , w ) = (a) ( 2 / 5 , ~ 1 4 ) ;(b) (-4/5,d4). mial interpolation with( a , w ) = (0,~ 1 2 )(a)
. Double pole. (b) Pole with
multiplicity 5 .

In our case forthe signals that are responses of natural VI. NUMERICAL
EXAMPLES
systems, we can say the following. We consider the signal
If the poles of the systemare inside the region of the f(t) = tieaTsin (wt + 4) Vt E R, (66)
complex plane (aT, U T ) depicted in Fig. 2, the car-
dinal interpolation (Shannon-Whittaker formula) and T = 1 .
converges withthe method of De la Vallee'Poussin. One-sided Sequences: Defining B = (esT- 11, the
Vice versa if the signal is represented by the cardinal condition for convergence for one-sided sequences is B
interpolation formula, the interpolating series (64) < 1.
will convergeuniformly to the same sum. There- Fig. 4 shows the maximum error between the polyno-
fore, the Newton-Gauss expansion represents a class mial that goesthroughthesamples { f(O), f(T), , --
of functions larger than the cardinal interpolation f ( n T ) ) . The signalf(t) has a simple pole (CY,w) = (i,$),
formula. As a particular case, the band-limited func- which corresponds to B > 1 . Fig. 4(a) shows the error
tions of the Shannon-Whittaker sampling theorem computed only in the interval [0, TI and Fig. 4(b) shows
will be located on the unit circle that is contained inthe error on all the range [0, n T ] V n = 0 , 1 , * , 21. --
the regions of convergence for the polynomialinter- The oscillating behavior is evident. Figs.5 , 6 , and 7 show
polation. the same calculation, respectively, for. (a,'3) = (- i, $) ,
(0, %),(3, 5) that correspond in all the cases to B < 1 .
Notice that the above result cannot be, in general, ap- The bound (35) is obviously decreasing and follows fairly
plied to one-sided sequences since the cardinal expansion well the behavior of the maximum error. Fig.8(a) and (b)
has a symmetrical kernel. refers to double poles for the error in the interval [0, 11.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
 856 TRANSACTIONS
IEEE ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-34, NO. 4, AUGUST 1986

Central PiecewiseApproximation: For the central [3] M. L. Cartwright, IntegralFunctions. NewYork:CambridgeUni-

piecewise approximation, we consider the polynomial of versity Press, 1956.
degree n that goes through the points { i - (n - 1)/2, , -- J. M. Whittaker, lnterpolafor))Function Theor). NewYork:Cam-
bridge University Press, 1935.
0, i + +
(n 1)/2}, and we compute the maximum error I. J.Schoenberg, CardinalSplineInterpolation. Philadelphia,PA:
+
in the interval [iT,( i 1)T] i = 0 , 1 , 2, 3 . Fig. 9(a) and SIAM,1973.
K. L. Kabir and L. P. Bolgiano, Jr., “The interpolation of signals
(b) shows two divergent cases in which the parameter A sampled at less than the Nyquist rate,” in Trans. Micro Delcon ’82,
in (49) is greater than one. Figs. 10 and 11 show conver- Newark, DE, 1982.
gent cases with the bound (55). Fig. 12(a) and (b) shows L. P. Bolgiano, Jr., “Interpolation of the natural response of physical
systems,” Numer. Function. Anal., Optimiz., vol. 2 , pp.323-329,
the maximum error for poles of multiplicity 2 and 5 . It Apr.1980.
. can be noticed that in the latter cases there is still conver- -, “A sampling theorem for system identification,” presented at
gence, but it is slower. the Int. Symp. Ill-PosedProblems,Univ.Delaware,Newark, Oct.
1979.
VII. CONCLUSIONS L. P. Bolgiano, Jr., and K. L. Kabir, “Computation of Fourier in-
tegralusingpolynomial interpolation,” in Proc. IEEE Inr. Con$
We have given in this paper a practical criterion to se- Acoust., Speech, Signal Processing, Washington DC, Apr. 1979, pp.
lect the sampling rate for polynomial interpolation of sig- 494-497.
K. L. Kabir, “A spectral computational technique using finite differ-
nals of the form (IO) that can be considered as the re- encesandthe Newton-Gregoryinterpolating series,” M.S. thesis,
sponse of systems described by linear differential equa- Dep. Elec. Eng., Univ. Delaware, Newark, June 1979.
tions with constant coefficients. A condition is derived for K. L. Kabirand L. P. Bolgiano, Jr., “A new approach for digital
processing-A new sampling theorem,” in Con$ Rec. Micro Delcon
uniform convergence of the one-sided interpolation pro- ’79, 2nd Annu. Conf., Comput. Soc. Chapter Delaware Bay, Mar.
cess that leads to a useful bound for the truncation error. 20,1979.
The same reasoning is used to derive a very simple con- L. P. Bolgiano, Jr., “Convolution preserving transformations, A new
interpretation,” presented at the Septieme Colloque Sur Le Traite-
dition on the sampling rate for symmetrical polynomial ment Du Signal Et Ses Applications, Nice, France, June 1979.
approximation. In practice, for both cases, the choice of L. P. Bolgiano,Jr. et al., “A nondifferentiablepulsewithamore
the sampling rate can be made by computing only one band-limited spectrum than a differentiable one,” Proc. IEEE, vol.
73, pp. 155-160, Jan. 1985.
parameter that defines also the rate of convergence for C. Lanczos, LinearDifferentialOperations. NewYork:Van Nos-
increasing order. Some numerical examples confirm the trand,1961.
simple theory. J. E. Hopcroft and K. Steiglitz, “A class of finite memory interpo-
lation filters,” IEEE Trans.CircuitTheory, vol.CT-15, pp, 105-
Finally, a direct comparison to the Shannon-Whittaker 111, June 1968.
sampling theorem based on the results reported by Whit- G. J. Joseph and L. P. Bolgiano, Jr., “Cubicsplinereconstruction
taker in [4] is presented. of a signal smoothed by a continuous physical system,”in Proc. 15th
Annu. Con$ Inform. Sci. Syst., John Hopkins Univ., Baltimore, MD,
It is the intention of the author to pursue or stimulate Mar. 1981, pp. 441-445.
further investigations for evaluating the effect that exter- R . W. Schafer and L. R . Rabiner, “A digital signal processing ap-
nal or modeling noise would produce on the interpolating proach tointerpolation,” Proc.IEEE, vol.61, pp.692-702,June
1973.
process. J. R. Higgins, “Five short stories about the cardinal series,” Bull.
The same convergence considerations could perhaps be Amer. Math. Soc., vol. 12, no. 1, pp. 45-89, 1985.
used to evaluate theeffectiveness of extrapolating the sig- J . McNamee, F. Stenger, and E. L. Whitney, “Whittaker’s cardinal
function in retrospect,” Math. Comput., vol. 25, pp. 141-154, Jan.
nal outside the range in which the samples are contained. 1971.
Interesting further investigations could also be aimed at D. L. Jagerman and L. J. Fogel, “Some general aspects of the sam-
deriving a practical sampling criterion forinterpolation by plingtheorem,” IRETrans.Inform.Theory, pp.139-146, vol. IT-
3, Dec.1956.
splines. E. T. Whittaker, “On the functions which are represented by the ex-
pansion in the interpolatory theory,” Proc. Roy. SOC.Edinburgh, vol.
ACKNOWLEDGMENT 35, pp. 181-194,1915.
This work is the result of the profound inspiration I re- F. Stenger, “Approximations via Whittaker’s cardinal function,” J.
Approx. Theory, vol. 17, pp. 222-240, 1976.
ceived from my advisor at the University of Delaware, L. A. J. Jerri, “The Shannon sampling theorem-Its various extensions
P. Bolgiano, Jr., who gave somecontributions on the topic and applications: A tutorial review,” Proc. IEEE, vol. 65, pp. 1565-
a few years ago. Therefore, I would like to thank him for 1596,Nov.1977.
C. E. Shannon, “Communication in presence of noise,” Proc. IRE,
his constant scientific and moral support during the time vol. 37. pp. 10-21, Jan. 1949.
of this work. A Papoulis, “Error analysis in sampling theory,” Proc. IEEE, vol.
54, pp. 947-955, July 1966.
I would like to thank Dott. R. Magnanini, a visiting K. F. Cheung and R . J. Marks, 11, ‘‘Ill-posed sampling theorems,”
scholar at the University of Delaware from the Diparti- IEEE Trans. Circuits Syst., vol. CAS-32, pp. 481-484, May 1985.
mento di Matematica dell’universita’ di Firenze in Italy, S. D. Stearns, Digital Signal Analysis. Rochelle Park, NJ: Hayden,
for some useful discussions on the topic. 1975.
DelaVallCe Poussin, “Quasi analytic functions,”. Rice Inst. Pam-
REFERENCES phlet, vol. XII, pp. 105- 172, Apr. 1925.
E. De Castro, Complementi di Analisi Matematica. Bologna, Italy:
[ l j I. S. Berezin and N. P. Zhidkov, Computing Methods, Vol. 1 . Elms- Zanichelli, 1972.
ford, NY: Pergamon 1965. I. J. Schoenberg, “Contribution to the problem of approximation of
[2] P. M.Prenter, Splines and Variational Methods. New York: Wiley? equidistant data byanalytical functions,” Quart. Appl. Math., vol.
1975. IV, part A, pp. 45-49, part B, pp. 112-141, 1946.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.
SAMPLING PALMIERI: FOR POLYNOMIAL INTERPOLATION 857

[31] J . R. Ragazzini and G.F. Franklin, Sampled-Data Control Systems. Francesco Palmieri was born in S. Maria C. V.,
New York: McGraw-Hill, 1958. Caserta, Italy, on August 28, 1956. He received
[32] H. D. Helms and J. B. Thomas, “Truncation error of sampling-theo- the Laurea in Ingegneria Elettronica degree from
rem expansions,” Proc. IRE, vol. 50, pp. 179-184, Feb. 1962. the Universita’ degli Studi di Napoli, Naples, It-
[33] K. Yao and J. B. Thomas, “On truncation error for sampling repre- aly,in1980.In1983hereceivedaFullbright
sentations of band-limited signals,” ZEEE Trans. Aerosp. Electron. scholarship as part of the program of cultural ex-
SySt. V O ~AES-2,
. pp. 640-647, NOV. 1966. change between Italy and the United States, and
[34] A. Papoulis, “Truncated sampling expansions,” ZEEE Trans. Auto- he joined the University of Delaware, Newark,
mat. Contr., vol. AC-12, pp. 604-605, Oct. 1967. where he received the Master’s degree in applied
1351 J. L. Brown, “Bounds for truncation error in sampling expansions of sciences in electrical engineering in June 1985. He
band-limited signals,” ZEEE Trans. Inform. Theory, vol. IT-15, pp. is currently completing the Ph.D. degree at the
440-444, July 1969. same University.
[36] A. J. Lee, “Approximate interpolation and the sampling theorem,” From 1982 to 1983 he worked as a Designer for digital communication
SZAM J. Appl. Math., vol. 32, June 1977. systems with the ITT firms: Face Sud Selettronica, in Salerno, Italy, and
Bell Telephone Man. Co. in Antwerpen, Belgium. His research interests
are in the fieldof communications and signal processing. In particular his
experience includes investigations on: speech processing for speaker rec-
ognition, spectral estimation, digital filter design, image processing with
emphasis on problems of reconstruction from projections, sampling and
reconstruction techniques, and robust filtering.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on January 03,2025 at 04:53:40 UTC from IEEE Xplore. Restrictions apply.