Sampling—50 Years After Shannon
MICHAEL UNSER, FELLOW, IEEE
This paper presents an account of the current state of sampling,
50 years after Shannon’s formulation of the sampling theorem. The
emphasis is on regular sampling, where the grid is uniform. This
topic has benefited from a strong research revival during the past
few years, thanks in part to the mathematical connections that were
made with wavelet theory. To introduce the reader to the modern,
Hilbert-space formulation, we reinterpret Shannon’s sampling pro-
cedure as an orthogonal projection onto the subspace of band-lim-
ited functions. We then extend the standard sampling paradigm for
a representation of functions in the more general class of “shift-in-
variant” functions spaces, including splines and wavelets. Practi-
cally, this allows for simpler—and possibly more realistic—inter-
polation models, which can be used in conjunction with a much
wider class of (anti-aliasing) prefilters that are not necessarily ideal
low-pass. We summarize and discuss the results available for the
determination of the approximation error and of the sampling rate
when the input of the system is essentially arbitrary; e.g., nonban-
dlimited. We also review variations of sampling that can be under-
stood from the same unifying perspective. These include wavelets,
multiwavelets, Papoulis generalized sampling, finite elements, and
frames. Irregular sampling and radial basis functions are briefly
mentioned.
Keywords—Band-limited functions, Hilbert spaces, interpola-
tion, least squares approximation, projection operators, sampling,
sampling theorem, Shannon, splines, wavelets.
I. INTRODUCTION
In 1949, Shannon published the paper “Communication in
the Presence of Noise,” which set the foundation of informa-
tion theory [105], [106]. This paper is a masterpiece; both in
terms of achievement and conciseness. It is undoubtedly one
of the theoretical works that has had the greatest impact on
modern electrical engineering [145]. In order to formulate his
rate/distortion theory, Shannon needed a general mechanism
for converting an analog signal into a sequence of numbers.
This led him to state the classical sampling theorem at the
very beginning of his paper in the following terms.
Theorem 1 [Shannon]: If a function
contains no fre-
quencies higher than
(in radians per second), it is com-
pletely determined by giving its ordinates at a series of points
spaced
seconds apart.
Manuscript received September 17, 1999; revised January 4, 2000.
The author is with the Biomedical Imaging Group, Swiss Federal Institute
of Technology Lausanne CH-1015 Lausanne EPFL, Switzerland (e-mail:
Michael.Unser@epfl.ch).
Publisher Item Identifier S 0018-9219(00)02874-7.
While Shannon must get full credit for formalizing this
result and for realizing its potential for communication
theory and for signal processing, he did not claim it as his
own. In fact, just below the theorem, he wrote: “this is a
fact which is common knowledge in the communication
art.” He was also well aware of equivalent forms of the
theorem that had appeared in the mathematical literature;
in particular, the work of Whittaker [144]. In the Russian
literature, this theorem was introduced to communication
theory by Kotel’nikov [67], [68].
The reconstruction formula that complements the sam-
pling theorem is
(1)
in which the equidistant samples of
may be interpreted
as coefficients of some basis functions obtained by appro-
priate shifting and rescaling of the sinc-function: sinc
. Formula (1) is exact if is bandlimited to
; this upper limit is the Nyquist frequency, a
term that was coined by Shannon in recognition of Nyquist’s
important contributions in communication theory [88]. In the
mathematical literature, (1) is known as the cardinal series
expansion; it is often attributed to Whittaker in 1915 [26],
[143] but has also been traced back much further [14], [58].
Shannon’s sampling theorem and its corresponding recon-
struction formula are best understood in the frequency do-
main, as illustrated in Fig. 1. A short reminder of the key
sampling formulas is provided in Appendix A to make the
presentation self-contained.
Nowadays the sampling theorem plays a crucial role in
signal processing and communications: it tells us how to
convert an analog signal into a sequence of numbers, which
can then be processed digitally—or coded—on a computer.
While Shannon’s result is very elegant and has proven to
be extremely fruitful, there are several problems associated
with it. First, it is an idealization: real world signals or
images are never exactly bandlimited [108]. Second, there is
no such device as an ideal (anti-aliasing or reconstruction)
low-pass filter. Third, Shannon’s reconstruction formula is
rarely used in practice (especially with images) because of
the slow decay of the sinc function [91]. Instead, practi-
tioners typically rely on much simpler techniques such as
0018–9219/00$10.00 © 2000 IEEE
PROCEEDINGS OF THE IEEE, VOL. 88, NO. 4, APRIL 2000 569
Fig. 1. Frequency interpretation of the sampling theorem: (a)
Fourier transform of the analog input signal
f
(
x
)
, (b) the sampling
process results in a periodization of the Fourier transform, and (c)
the analog signal is reconstructed by ideal low-pass filtering; a
perfect recovery is possible provided that
!
=T
.
linear interpolation. Despite these apparent mismatches with
the physical world, we will show that a reconciliation is
possible and that Shannon’s sampling theory, in its modern
and extended versions, can perfectly handle such “nonideal”
situations.
Ten to 15 years ago, the subject of sampling had reached
what seemed to be a very mature state [26], [62]. The re-
search in this area had become very mathematically oriented,
with less and less immediate relevance to signal processing
and communications. Recently, there has been strong revival
of the subject, which was motivated by the intense activity
taking place around wavelets (see [7], [35], [80], and [85]).
It soon became clear that the mathematics of wavelets were
also applicable to sampling, but with more freedom because
no multiresolution is required. This led researchers to reex-
amine some of the foundations of Shannon’s theory and de-
velop more general formulations, many of which turn out to
be quite practical from the point of view of implementation.
Our goal in this paper is to give an up-to-date account of
the recent advances that have occurred in regular sampling.
Here, the term “regular” refers to the fact that the samples
are taken on a uniform grid—the situation most commonly
encountered in practice. While the paper is primarily con-
ceived as a tutorial, it contains a fair amount of review mate-
rial—mostly recent work: This should make it a useful com-
plement to the excellent survey article of Jerri, which gives
a comprehensive overview of sampling up to the mid-1970’s
[62].
The outline of this paper is as follows. In Section II, we
will argue that the requirement of a perfect reconstruction
is an unnecessarily strong constraint. We will reinterpret the
standard sampling system, which includes an anti-aliasing
prefilter, as an orthogonal projection operator that computes
the minimum error band-limited approximation of a not-nec-
essarily band-limited input signal. This is a crucial obser-
vation that changes our perspective: instead of insisting that
the reconstruction be exact, we want it to be as close as pos-
sible to the original; the global system, however, remains un-
changed, except that the input can now be arbitrary. (We can
obviously not force it to be bandlimited.)
In Section III, we will show that the concept extends nicely
to the whole class of spline-like (or wavelet-like) spaces gen-
erated from the integer shifts of a generating function. We
will describe several approximation algorithms, all based on
the standard three-step paradigm: prefiltering, sampling, and
postfiltering—the only difference being that the filters are
not necessarily ideal. Mathematically, these algorithms can
all be described as projectors. A direct consequence is that
they reconstruct all signals included within the reconstruc-
tion space perfectly—this is the more abstract formulation
of Shannon’s theorem.
In Section IV, we will investigate the issue of approxima-
tion error, which becomes relevant once we have given up
the goal of a perfect reconstruction. We will present recent
results in approximation theory, making them accessible to
an engineering audience. This will give us the tools to select
the appropriate sampling rate and to understand the effect of
different approximation or sampling procedures.
Last, in Section V, we will review additional extensions
and variations of sampling such as (multi)wavelets, finite el-
ements, derivative and interlaced sampling, and frames. Ir-
regular sampling will also be mentioned, but only briefly, be-
cause it is not the main focus of this paper. Our list of sam-
pling topics is not exhaustive—for instance, we have com-
pletely left out the sampling of discrete sequences and of
stochastic processes—but we believe that the present paper
covers a good portion of the current state of research on
regular sampling. We apologize in advance to those authors
whose work was left out of the discussion.
II. S
HANNON’S SAMPLING THEOREM REVISITED
Shannon’s sampling theory is applicable whenever the
input function is bandlimited. When this is not the case, the
standard signal-processing practice is to apply a low-pass
filter prior to sampling in order to suppress aliasing. The
optimal choice is the ideal filter
, which
suppresses aliasing completely without introducing any
distortion in the bandpass region. Its impulse response is
. The corresponding block diagram is shown
in Fig. 2. In this section, we provide a geometrical Hilbert
space interpretation of the standard sampling paradigm. For
notational simplicity, we will set
and rescale the time
dimension accordingly.
In 1941, the English mathematician Hardy, who was re-
ferring to the basis functions in Whittaker’s cardinal series
(1), wrote: “It is odd that, although these functions occur re-
peatedly in analysis, especially in the theory of interpolation,
it does not seem to have been remarked explicitly that they
form an orthogonal system” [55]. Orthonormality is a fun-
damental property of the sinc-function that has been revived
recently.
To understand the modern point of view, we have to con-
sider the Hilbert space
, which consists of all functions that
570 PROCEEDINGS OF THE IEEE, VOL. 88, NO. 4, APRIL 2000
Fig. 2. Schematic representation of the standard three-step
sampling paradigm with
T
=1
: 1) the analog input signal is
prefiltered with
h
(
x
)
(anti-aliasing step), 2) the sampling process
yields the sampled representation
c
(
x
)=
c
(
k
)
(
x
0
k
)
,
and 3) the reconstructed output
~
f
(
x
)=
c
(
k
)
'
(
x
0
k
)
is obtained by analog filtering of
c
with
'
. In the traditional
approach, the pre- and postfilters are both ideal low-pass:
h
(
x
)=
'
(
x
)=sinc(
x
)
. In the more modern schemes, the filters
can be selected more freely under the constraint that they remain
biorthogonal:
h
'
(
x
0
k
)
;
~
'
(
x
0
l
)
i
=
.
are squareintegrable in Lebesgue’s sense. The corresponding
-norm is
(2)
It is induced by the conventional
-inner product
(3)
We now assume that the input function that we want to
sample is in
, a space that is considerably larger than the
usual subspace of band-limited functions, which we will
refer to as
. This means that we will need to make an
approximation if we want to represent a non-band-limited
signal in the band-limited space
. To make an analogy,
is to what (the real numbers) is to (the integers).
The countable nature of
is apparent if we rewrite the
normalized form of (1) with
as
with (4)
where the
’s are some coefficients, and where the ’s
are the basis functions to which Hardy was referring. It is not
difficult to show that they are orthonormal
(5)
where the autocorrelation function is evaluated as follows:
This orthonormality property greatly simplifies the imple-
mentation of the approximation process by which a function
is projected onto . Specifically, the orthogonal pro-
jection operator
can be written as
(6)
where the inner product
represents the signal
contribution along the direction specified by
—the ap-
proximation problem is decoupled component-wise because
of the orthogonality of the basis functions. The projection
theorem (see [69]) ensures that this projection operation is
well defined and that it yields the minimum-error approxi-
mation of
into .
(7)
By a lucky coincidence, this inner product computation is
equivalent to first filtering the input function
with the ideal
low-pass filter and sampling thereafter. More generally, we
observe that any combination of prefiltering and sampling
can be rewritten in terms of inner products
with (8)
That is, the underlying analysis functions correspond to the
integer shifts of
, the time-reversed impulse
response of the prefilter
(which can be arbitrary). In the
present case,
and is the
ideal low-pass filtered version of
.
The conclusion of this section is that the traditional sam-
pling paradigm with ideal prefiltering yields an approxima-
tion
, which is the orthogonal projection of the
input function onto
(the space of band-limited functions).
In other words,
is the approximation of in with min-
imum error. In light of this geometrical interpretation, it is
obvious that
(since is a projection
operator), which is a more concise statement of Shannon’s
theorem.
III. S
AMPLING IN SHIFT-INVARIANT (OR SPLINE-LIKE)
S
PACES
Having reinterpreted the sampling theorem from the more
abstract perspective of Hilbert spaces and of projection op-
erators, we can take the next logical step and generalize the
approach to other classes of functions.
A. Extending Shannon’s Model
While it would be possible to consider arbitrary basis func-
tions, we want a sampling scheme that is practical and retains
the basic, shift-invariant flavor of the classical theory. This is
achieved by simply replacing
by a more general tem-
plate: the generating function
. Accordingly, we specify
our basic approximation space
as
(9)
This means that any function
, which is con-
tinuously defined, is characterized by a sequence of coeffi-
cients
; this is the discrete signal representation that will
be used to do signal processing calculations or to perform
coding. Note that the
’s are not necessarily the samples
UNSER: SAMPLING—50 YEARS AFTER SHANNON 571
of the signal, and that can be something quite different
from
. Indeed, one of our motivations is to discover
functions that are simpler to handle numerically and have a
much faster decay.
For such a continuous/discrete model to make sense, we
need to put a few mathematical safeguards. First, the se-
quence of coefficients must be square-summable:
.
Second, the representation should be stable
1
and unambigu-
ously defined. In other words, the family of functions
should form a Riesz basis of , which
is the next best thing after an orthogonal one [35]. The def-
inition of a Riesz basis is that there must exist two strictly
positive constants
and such that
(10)
where
is the squared -norm (or en-
ergy) of
. A direct consequence of the lower inequality
is that
implies . Thus, the basis
functions are linearly independent, which also means that
every signal
is uniquely specified by its co-
efficients
. The upper bound in (10) implies that the
-norm of the signal is finite, so that is a valid sub-
space of
. Note that the basis is orthonormal if and only if
, in which case we have a perfect norm equiva-
lence between the continuous and the discrete domains (Par-
seval’s relation). Because of the translation-invariant struc-
ture of the construction, the Riesz basis requirement has an
equivalent expression in the Fourier domain [9]
(11)
where
is the Fourier transform of
. Note that the central term in (11) is the Fourier trans-
form of the sampled autocorrelation function
(12)
It can therefore also be written
[see (A.7) in Appendix A].
The final requirement is that the model should have the
capability of approximating any input function as closely as
desired by selecting a sampling step that is sufficiently small
(similar to the Nyquist criterion
). As shown
in Appendix B, this is equivalent to the partition of unity
condition
(13)
In practice, it is this last condition that puts the strongest
constraint of the selection on an admissible generating func-
tion
.
Let usnow lookat some examples. The firstone, which has
already been discussed in great length, is the classical choice
1
By stable, we mean that a small variation of the coefficients must result
in a small variation of the function. Here, the upper bound of the Riesz con-
dition (10) ensures
L
-stability.
. It is easy to verify that the corresponding
Riesz bounds in (11) are
, which is consistent
with the orthonormality property (5). We now show that the
sinc function satisfies the partition of unity: using Poisson’s
summation formula (cf. Appendix A), we derive an equiva-
lent formulation of (13) in the Fourier domain
2
(14)
a relation that is obviously satisfied by
, the
Fourier transform of
.
The sinc functionis well localized in the frequencydomain
but has very poor time decay. At the other extreme, we can
look for the simplest and shortest function that satisfies (13).
It is the box function (or B-spline of degree 0)
.
(15)
The corresponding basis functions are clearly orthogonal for
they do not overlap.
By convolving this function with itself repeatedly, we con-
struct the B-splines of degree
, which are defined recur-
sively as
(16)
These functions are known to generate the polynomial
splines with equally spaced knots [98], [99]. Specifically,
if
, then the signals defined by (9) are
polynomials of degree
within each interval (
odd), respectively, (1/2) (1/2) when is even,
with pieces that are patched together such as to guarantee
the continuity of the function and its derivatives up to order
(i.e., ). The B-spline basis functions up
to degree 4 are shown in Fig. 3. They are symmetric and
well localized, but not orthogonal—except for
. Yet,
they all satisfy the Riesz basis condition and the partition
of unity. This last property is easily verified in the Fourier
domain using (14). The B-splines are frequently used in
practice (especially for image processing) because of their
short support and excellent approximation properties [126].
The B-spline of degree 1 is the tent function, and the
corresponding signal model is piecewise linear. This rep-
resentation is quite relevant because linear interpolation is
one of the most commonly used algorithm for interpolating
signal values.
As additional examples of admissible generating func-
tions, we may consider any scaling function (to be defined
below) that appears in the theory of the wavelet transform
[35], [79], [115], [139]. It is important to note, however,
that scaling functions, which are often also denoted by
, satisfy an additional two-scale relation, which is
not required here but not detrimental either. Specifically,
a scaling function is valid (in the sense defined by Mallat
2
The equality on the right hand side of (14) holds in the distributional
sense provided that
'
(
x
)
2
L
, or by extension, when
'
(
x
)+
'
(
x
+1)
2
L
, which happens to be the case for
'
(
x
)=sinc(
x
)
572 PROCEEDINGS OF THE IEEE, VOL. 88, NO. 4, APRIL 2000
Fig. 3. The centered B-splines for
n
=0
to
4
. The B-splines of
degree
n
are supported in the interval
[
0
((
n
+1)
=
2)
;
((
n
+1)
=
2)]
;
as
n
increases, they flatten out andget more and more Gaussian-like.
[81]) if and only if 1) it is an admissible generating function
(Riesz basis condition
partition of unity) and 2) it satisfies
the two-scale relation
(17)
where
is the so-called refinement filter. In other words,
the dilated version of
must live in the space , a prop-
erty that is much more constraining than the conditions im-
posed here.
B. Minimum Error Sampling
Having defined our signal space, the next natural question
is how to obtain the
’s in (9) such that the signal model
is a faithful approximation of some input function
. The optimal solution in the least squares sense is the
orthogonal projection, which can be specified as
(18)
where the
’s are the dual basis functions of .
This is very similar to (6), except that the analysis and syn-
thesis functions (
and , respectively) are not identical—in
general, the approximation problem is not decoupled. The
dual basis
with is unique and is deter-
mined by the biorthogonality condition
(19)
It also inherits the translation-invariant structure of the basis
functions:
.
Since
, it is a linear conbination of the ’s.
Thus, we may represent it as
where is a suitable sequence to be determined next. Let
us therefore evaluate the inner product
where is the autocorrelation sequence in (12). By
imposing the biorthogonality constraint
,
and by solving this equation in the Fourier domain [i.e.,
], we find
that
(20)
Note that the Riesz basis condition (11) guarantees that this
solution is always well defined [i.e., the numerator on the
right-hand side of (20) is bounded and nonvanishing].
Similar to what hasbeensaid for the band-limited case [see
(4)], the algorithm described by (18) has a straightforward
signal-processing interpretation (see Fig. 2). The procedure
is exactly the same as the one dictated by Shannon’s theory
(with anti-aliasing filter), except that the filters are not nec-
essarily ideal anymore. Note that the optimal analysis filter
is entirely specified by the choice of the generating func-
tion (reconstruction filter); its frequency response is given
by (20). If
is orthonormal, then it is its own analysis func-
tion (i.e.,
) and the prefilter is simply a flipped
version of the reconstruction filter. For example, this implies
that the optimal prefilter for a piecewise constant signal ap-
proximation is a box function.
C. Consistent Sampling
We have just seen how to design an optimal sampling
system. In practice however, the analog prefilter is often
specified a priori (acquisition device), and not necessarily
optimal or ideal. We will assume that the measurements of
a function
are obtained by sampling its prefiltered
version,
, which is equivalent to computing the series
of inner products
(21)
with analysis function
[see (8)]. Here, it is
important to specify the input space
such that the mea-
surements are square-summable:
. In the
most usual cases (typically,
), we will be able to
consider
; otherwise ( is a Delta Dirac or a dif-
ferential operator), we may simply switch to a slightly more
constrained Sobolev space.
3
Now, given the measurements
in (21), we want to construct a meaningful approximation of
the form (9) with synthesis function
. The solution is
to apply a suitable digital correction filter
, as shown in the
block diagram in Fig. 4.
Here, we consider a design based on the idea of consis-
tency [127]. Specifically, one seeks a signal approximation
that is such that it would yield exactly the same measure-
ments if it was reinjected into the system. This is a reason-
able requirement, especially when we have no other way of
probing the input signal: if it “looks” the same, we may as
well say that it is the same for all practical purposes.
3
The Sobolev space
W
specifies the class of functions that
are
r
times differentiable in the
L
-sense. Specifically,
W
=
f
f
:
(1+
!
)
j
^
f
(
!
)
j
d! <
+
1g
where
^
f
(
!
)
is the Fourier transform
of
f
(
x
)
.
UNSER: SAMPLING—50 YEARS AFTER SHANNON 573
![Fig. 3. The centered B-splines forn = 0 to 4. The B-splines of degreen are supported in the interval[ ((n+1)=2); ((n+1)=2)]; asn increases, they flatten out and get more and more Gaussian-like.](/figures/fig-3-the-centered-b-splines-forn-0-to-4-the-b-splines-of-1u3hpmf0.webp)