Nonuniform Sampling and Reconstruction in Shift-Invariant Spaces

TL;DR: A unified framework for uniform and nonuniform sampling and reconstruction in shift-invariant subspaces is provided by bringing together wavelet theory, frame theory, reproducing kernel Hilbert spaces, approximation theory, amalgam spaces, and sampling.

Abstract: This article discusses modern techniques for nonuniform sampling and reconstruction of functions in shift-invariant spaces. It is a survey as well as a research paper and provides a unified framework for uniform and nonuniform sampling and reconstruction in shift-invariant subspaces by bringing together wavelet theory, frame theory, reproducing kernel Hilbert spaces, approximation theory, amalgam spaces, and sampling. Inspired by applications taken from communication, astronomy, and medicine, the following aspects will be emphasized: (a) The sampling problem is well defined within the setting of shift-invariant spaces. (b) The general theory works in arbitrary dimension and for a broad class of generators. (c) The reconstruction of a function from any sufficiently dense nonuniform sampling set is obtained by efficient iterative algorithms. These algorithms converge geometrically and are robust in the presence of noise. (d) To model the natural decay conditions of real signals and images, the sampling theory is developed in weighted L p-spaces.

Summary

Introduction

• Nonuniform sampling, irregular sampling, sampling, reconstruction, wavelets, shift-invariant spaces, frame, reproducing kernel Hilbert space, weighted Lp-spaces, amalgam spaces AMS subject classifications.
• More information about modern techniques for nonuniform sampling and applications can be found in [16].
• These sequences can then be processed digitally and converted back to analog signals via (1.1).

In particular, if φ is continuous and has compact support, then the conclusions (i)–(iii) hold.

• This implies that the sampling is stable or, equivalently, that the reconstruction of f from its samples is continuous.
• The authors need to design reconstruction procedures that are useful and efficient in practical applications.
• Remark 3.1. (i) The hypothesis that X be separated is for convenience only and is not essen- tial.
• For arbitrary sampling sets, the authors can use adaptive weights to compensate for the local variations of the sampling density [48, 49].

Figures

Fig. 6.3 Missing data reconstruction. Top left: Original digital MRI image with 128×128 samples. Top right: MRI image with 50% randomly missing samples. Bottom left: Reconstruction using the iterative reconstruction algorithm (6.2) of Theorem 6.1. The corresponding shiftinvariant space is generated by φ(x, y) = β3(x)× β3(y), where β3 = χ[0,1] ∗ χ[0,1] ∗ χ[0,1] ∗ χ[0,1] is the B-spline function of degree 3.

Fig. 1.1 The sampling problem. Top: A function f defined on R has been sampled on a uniform grid. Bottom: The same function f has been sampled on a nonuniformly spaced set. The sampling locations xj are marked by the symbol ×, and the sampled values f(xj) by a circle o.

Fig. 6.2 Reconstruction of a function f with ‖f‖2 ≈ 3.5 using the iteration algorithm (6.2) of Theorem 6.1. Top left: Function f belonging to the shift-invariant space generated by the Gaussian function e−x 2/2σ2 , σ ≈ 0.81, and its sample values {f(xj) : xj ∈ X} marked by + (density γ ≈ 0.8). Top right: Error ‖f − fn‖L2 against the number of iterations. Bottom left: Final error f − fn after 10 iterations. Bottom right: Reconstructed function f10 (continuous line) and original samples {f(xj) : xj ∈ X}.

Fig. 6.1 The iterative reconstruction algorithm of Theorem 6.1.

Fig. 1.3 The missing data problem. Left: Original digital MRI image with 128×128 samples. Right: MRI image with 50% randomly missing samples.

Fig. 1.2 Sampling grids. Top left: Because of its simplicity the uniform Cartesian sampling grid is used in signal and image processing whenever possible. Top right: A polar sampling grid used in computerized tomography (see [90]). In this case, the two-dimensional Fourier transform f̂ is sampled with the goal of reconstructing f . Bottom left: Spiral sampling used for fast MRI by direct signal reconstruction from spectral data on spirals [21]. Bottom right: A typical nonuniform sampling set as encountered in spectroscopy, astronomy, geophysics, and other signal and image processing applications.

Fig. 1.4 Sampling and boundary reconstruction from ultrasonic images. Left: Detected edge points of the left ventricle of a heart from a two-dimensional ultrasound image constitute a nonuniform sampling of the left ventricle’s contour. Right: Boundary of the left ventricle reconstructed from the detected edge sample points (see [59]).

Fig. 7.1 Reconstruction of a function f with additive noise using the iterative algorithm (6.2) of Theorem 6.1. Top left: Function f belonging to the shift-invariant space generated by the Gaussian function e−x 2/2σ2 , σ ≈ 0.81. Top right: Function f with an additive white noise (SNR ≈ 0db). Bottom left: Noisy signal sampled on a nonuniform grid with maximal gap ≈ 0.51. Bottom right: Reconstructed function f10 after 10 iterations (continuous line) and original signal f (dotted line).

Full text

SIAM REVIEW
c
 2001 Society for Industrial and Applied Mathematics
Vol. 43, No. 4, pp. 585–620
Nonuniform Sampling and
Reconstruction in Shift-Invariant
Spaces
∗
Akram Aldroubi
†
Karlheinz Gr
¨
ochenig
‡
Abstract. This article discusses modern techniques for nonuniform sampling and reconstruction of
functions in shift-invariant spaces. It is a survey as well as a research paper and provides
a unified framework for uniform and nonuniform sampling and reconstruction in shift-
invariant spaces by bringing together wavelet theory, frame theory, reproducing kernel
Hilbert spaces, approximation theory, amalgam spaces, and sampling. Inspired by appli-
cations taken from communication, astronomy, and medicine, the following aspects will be
emphasized: (a) The sampling problem is well defined within the setting of shift-invariant
spaces. (b) The general theory works in arbitrary dimension and for a broad class of gener-
ators. (c) The reconstruction of a function from any sufficiently dense nonuniform sampling
set is obtained by efficient iterative algorithms. These algorithms converge geometrically
and are robust in the presence of noise. (d) To model the natural decay conditions of real
signals and images, the sampling theory is developed in weighted L
p
-spaces.
Key words. nonuniform sampling, irregular sampling, sampling, reconstruction, wavelets, shift-invar-
iant spaces, frame, reproducing kernel Hilbert space, weighted L
p
-spaces, amalgam spaces
AMS subject classifications. 41A15,42C15, 46A35, 46E15, 46N99, 47B37
PII. S0036144501386986
1. Introduction. Modern digital data processing of functions (or signals or im-
ages) always uses a discretized version of the original signal f that is obtained by
sampling f on a discrete set. The question then arises whether and how f can be
recovered from its samples. Therefore, the objective of research on the sampling prob-
lem is twofold. The first goal is to quantify the conditions under which it is possible to
recover particular classes of functions from different sets of discrete samples. The sec-
ond goal is to use these analytical results to develop explicit reconstruction schemes
for the analysis and processing of digital data. Specifically, the sampling problem
consists of two main parts:
(a) Given a class of functions V on R
d
, find conditions on sampling sets X =
{x
j
∈ R
d
: j ∈ J}, where J is a countable index set, under which a function
f ∈ V can be reconstructed uniquely and stably from its samples {f(x
j
) :
x
j
∈ X}.
∗
Received by the editors January 24, 2001; accepted for publication (in revised form) April 2, 2001;
published electronically October 31, 2001. The U.S. Government retains a nonexclusive, royalty-free
license to publish or reproduce the published form of this contribution, or allow others to do so, for
U.S. Government purposes. Copyright is owned by SIAM to the extent not limited by these rights.
http://www.siam.org/journals/sirev/43-4/38698.html
†
Department of Mathematics, Vanderbilt University, Nashville, TN 37240 (aldroubi@math.
vanderbilt.edu). This author’s research was supported in part by NSF grant DMS-9805483.
‡
Department of Mathematics U-3009, University of Connecticut, Storrs, CT 06269-3009 (groch@
math.uconn.edu).
585

586 AKRAM ALDROUBI AND KARLHEINZ GR
¨
OCHENIG
0 100 200 300 400 500 600 700 800 900 1000
0
0.5
1
1.5
2
0 100 200 300 400 500 600 700 800 900 1000
0
0.5
1
1.5
2
Fig. 1.1 The sampling problem. Top: A function f defined on R has been sampled on a uniform
grid. Bottom: The same function f has been sampled on a nonuniformly spaced set. The
sampling locations x
j
are marked by the symbol ×, and the sampled values f (x
j
) by a
circle o.
(b) Find efficient and fast numerical algorithms that recover any function f ∈ V
from its samples on X.
In some applications, it is justified to assume that the sampling set X = {x
j
: j ∈ J}
is uniform, i.e., that X forms a regular n-dimensional Cartesian grid; see Figures
1.1 and 1.2. For example, a digital image is often acquired by sampling light inten-
sities on a uniform grid. Data acquisition requirements and the ability to process
and reconstruct the data simply and efficiently often justify this type of uniform data
collection. However, in many realistic situations the data are known only on a nonuni-
formly spaced sampling set. This nonuniformity is a fact of life and prevents the use of
the standard methods from Fourier analysis. The following examples are typical and
indicate that nonuniform sampling problems are pervasive in science and engineering.
• Communication theory: When data from a uniformly sampled signal (func-
tion) are lost, the result is generally a sequence of nonuniform samples. This
scenario is usually referred to as a missing data problem. Often, missing sam-
ples are due to the partial destruction of storage devices, e.g., scratches on
a CD. As an illustration, in Figure 1.3 we simulate a missing data problem
by randomly removing samples from a slice of a three-dimensional magnetic
resonance (MR) digital image.
• Astronomical measurements: The measurement of star luminosity gives rise
to extremely nonuniformly sampled time series. Daylight periods and adverse
nighttime weather conditions prevent regular data collection (see, e.g., [111]
and the references therein).
• Medical imaging: Computerized tomography (CT) and magnetic resonance
imaging (MRI) frequently use the nonuniform polar and spiral sampling sets
(see Figure 1.2 and [21, 90]).

NONUNIFORM SAMPLING AND RECONSTRUCTION 587
0 2 4 6 8
0
1
2
3
4
5
Cartesian uniform sampling grid
0 2 4 6 8
0
1
2
3
4
5
Nonuniform sampling grid
-5 0 5
-5
0
5
Polar sampling grid
-5 0 5
-5
0
5
Spiral sampling grid
Fig. 1.2 Sampling grids. Top left: Because of its simplicity the uniform Cartesian sampling grid
is used in signal and image processing whenever possible. Top right: A polar sampling
grid used in computerized tomography (see [90]). In this case, the two-dimensional Fourier
transform
ˆ
f is sampled with the goal of reconstructing f . Bottom left: Spiral sampling used
for fast MRI by direct signal reconstruction from spectral data on spirals [21]. Bottom right:
A typical nonuniform sampling set as encountered in spectroscopy, astronomy, geophysics,
and other signal and image processing applications.
Original digital image Digital image with missing data
Fig. 1.3 The missing data problem. Left: Original digital MRI image with 128×128 samples. Right:
MRI image with 50% randomly missing samples.
Other applications using nonuniform sampling sets occur in geophysics [92], spec-
troscopy [101], general signal/image processing [13, 22, 103, 106], and biomedical
imaging [20, 59, 90, 101] (see Figures 1.2 and 1.4). More information about modern
techniques for nonuniform sampling and applications can be found in [16].

588 AKRAM ALDROUBI AND KARLHEINZ GR
¨
OCHENIG
Fig. 1.4 Sampling and boundary reconstruction from ultrasonic images. Left: Detected edge points
of the left ventricle of a heart from a two-dimensional ultrasound image constitute a nonuni-
form sampling of the left ventricle’s contour. Right: Boundary of the left ventricle recon-
structed from the detected edge sample points (see [59]).
1.1. Sampling in Paley–Wiener Spaces: Bandlimited Functions. Since infi-
nitely many functions can have the same sampled values on X = {x
j
}
j∈J
⊂ R
d
, the
sampling problem becomes meaningful only after imposing some a priori conditions
on f . The standard assumption is that the function f on R
d
belongs to the space
of bandlimited functions B
Ω
; i.e., the Fourier transform
ˆ
f(ξ) =
R
R
d
f(x)e
−2πihξ,xi
dx
of f is such that
ˆ
f(ξ) = 0 for all ξ /∈ Ω = [−ω, ω]
d
for some ω < ∞ (see, e.g.,
[15, 44, 47, 55, 62, 72, 78, 88, 51, 112] and the review papers [27, 61, 65]). The reason
for this assumption is a classical result of Whittaker [114] in complex analysis which
states that, for dimension d = 1, a function f ∈ L
2
(R) ∩ B
[−1/2,1/2]
can be recovered
exactly from its samples {f(k) : k ∈ Z} by the interpolation formula
(1.1) f(x) =
X
k∈Z
f(k) sinc(x − k),
where sinc(x) =
sin πx
πx
. This series gave rise to the uniform sampling theory of Shan-
non [96], which is fundamental in engineering and digital signal processing because
it gives a framework for converting analog signals into sequences of numbers. These
sequences can then be processed digitally and converted back to analog signals via
(1.1).
Taking the Fourier transform of (1.1) and using the fact that the Fourier transform
of the sinc function is the characteristic function χ
[−1/2,1/2]
shows that for any ξ ∈
[−1/2, 1/2]
ˆ
f(ξ) =
X
k
f(k)e
2πikξ
=
X
k
h
ˆ
f, e
i2πk·
i
L
2
(−1/2,1/2)
e
i2πkξ
.
Thus, reconstruction by means of the formula (1.1) is equivalent to the fact that
the set {e
i2πkξ
, k ∈ Z} forms an orthonormal basis of L
2
(−1/2, 1/2) called the har-
monic Fourier basis. This equivalence between the harmonic Fourier basis and the
reconstruction of a uniformly sampled bandlimited function has been extended to
treat some special cases of nonuniformly sampled data. In particular, the results by
Paley and Wiener [87], Kadec [71], and others on the nonharmonic Fourier bases
{e
i2πx
k
ξ
, k ∈ Z} can be translated into results about nonuniform sampling and re-
construction of bandlimited functions [15, 62, 89, 94]. For example, Kadec’s theorem

NONUNIFORM SAMPLING AND RECONSTRUCTION 589
[71] states that if X = {x
k
∈ R : |x
k
− k| ≤ L < 1/4} for all k ∈ Z, then the set
{e
i2πx
k
ξ
, k ∈ Z} is a Riesz basis of L
2
(−1/2, 1/2); i.e., {e
i2πx
k
ξ
, k ∈ Z} is the image
of an orthonormal basis of L
2
(−1/2, 1/2) under a bounded and invertible operator
from L
2
(−1/2, 1/2) onto L
2
(−1/2, 1/2). Using Fourier transform methods, this re-
sult implies that any bandlimited function f ∈ L
2
∩ B
[−1/2,1/2]
can be completely
recovered from its samples f (x
k
), k ∈ Z, as long as the sampling set is of the form
X = {x
k
∈ R : |x
k
− k| < 1/4}
k∈Z
.
The sampling set X = {x
k
∈ R : |x
k
− k| < 1/4}
k∈Z
in Kadec’s theorem is just
a perturbation of Z. For more general sampling sets, the work of Beurling [23, 24],
Landau [74], and others [18, 58] provides a deep understanding of the one-dimensional
theory of nonuniform sampling of bandlimited functions. Specifically, for the exact and
stable reconstruction of a bandlimited function f from its samples {f(x
j
) : x
j
∈ X},
it is sufficient that the Beurling density
(1.2) D(X) = lim
r→∞
inf
y∈R
#X ∩ (y + [0, r])
r
satisfies D(X) > 1. Conversely, if f is uniquely and stably determined by its samples
on X ⊂ R, then D(X) ≥ 1 [74]. The marginal case D(X) = 1 is very complicated
and is treated in [79, 89, 94].
It should be emphasized that these results deal with stable reconstructions. This
means that an inequality of the form
kfk
p
≤ C


X
x
j
∈X
|f(x
j
)|
p


1/p
holds for all bandlimited functions f ∈ L
p
∩ B
Ω
. A sampling set for which the recon-
struction is stable in this sense is called a (stable) set of sampling. This terminology
is used to contrast a set of sampling with the weaker notion of a set of uniqueness. X
is a set of uniqueness for B
Ω
if f|
X
= 0 implies that f = 0. Whereas a set of sampling
for B
[−1/2,1/2]
has a density D ≥ 1, there are sets of uniqueness with arbitrarily small
density. See [73, 25] for examples and characterizations of sets of uniqueness.
While the theorems of Paley and Wiener and Kadec about Riesz bases consisting
of complex exponentials e
i2πx
k
ξ
are equivalent to statements about sampling sets that
are perturbations of Z, the results about arbitrary sets of sampling are connected to
the more general notion of frames introduced by Duffin and Schaeffer [40]. The concept
of frames generalizes the notion of orthogonal bases and Riesz bases in Hilbert spaces
and of unconditional bases in some Banach spaces [2, 5, 6, 12, 14, 15, 20, 28, 29, 46,
66, 97].
1.2. Sampling in Shift-Invariant Spaces. The series (1.1) shows that the space
of bandlimited functions B
[−1/2,1/2]
is identical with the space
(1.3) V
2
(sinc) =
(
X
k∈Z
c
k
sinc(x − k) : (c
k
) ∈ `
2
)
.
Since the sinc function has infinite support and slow decay, the space of band-
limited functions is often unsuitable for numerical implementations. For instance, the
pointwise evaluation
f 7→ f(x
0
) =
X
k∈Z
c
k
sinc(x
0
− k)

Frequently asked questions

Q1. What have the authors contributed in "Nonuniform sampling and reconstruction in shift-invariant spaces∗" ?

This article discusses modern techniques for nonuniform sampling and reconstruction of functions in shift-invariant spaces. Inspired by applications taken from communication, astronomy, and medicine, the following aspects will be emphasized: ( a ) The sampling problem is well defined within the setting of shift-invariant spaces. 

Q2. What is the subspace of continuous functions W0(Lp)?

The subspace of continuous functions W0(Lpν) = W (C,L p ν) ⊂ W (Lpν) is a closed subspace of W (Lpν) and thus also a Banach space [43, 45]. 

Q3. What is the result of Lemma 2.10?

Assume that the submultiplicative weight ω satisfies the so-calledBeurling–Domar condition (mentioned in section 2.1)(2.14) ∞∑n=1logω(nk) n2 <∞ ∀ k ∈ Zd . 

Q4. What is the standard condition in wavelet theory?

For ν = 1 and p = 2, the standard condition in wavelet theory is often stated in the Fourier domain as (2.6) 0 < m ≤ âφ(ξ) = ∑j∈Zd |φ̂(ξ + j)|2 ≤M <∞ for almost every ξ,for some constants m > 0 and M > 0 [80, 81]. 

Q5. What is the earliest result of the theory of frames?

The earliest results [31, 77] concentrate on perturbation of regular sampling in shift-invariant spaces and are therefore similar in spirit to Kadec’s result for bandlimited functions. 

Q6. What is the significance of the Shannon sampling theory?

These results imply that Shannon’s sampling theory can be viewed as a limiting case of polynomial spline interpolation when the order of the spline tends to infinity [11, 109]. 

Q7. What are the connections between reproducing kernel Hilbert spaces and sampling?

Extensions of frame theory and reproducing kernel Hilbert spaces to Banach spaces are discussed, and the connections between reproducing kernels in weighted Lp-spaces, Banach frames, and sampling are described. 

Q8. What is the difference between the pointwise evaluationf 7 f(x0) and?

For instance, the pointwise evaluationf 7→ f(x0) = ∑k∈Z ck sinc(x0 − k)is a nonlocal operation, because, as a consequence of the long-range behavior of sinc, many coefficients ck will contribute to the value f(x0). 

Q9. What is the common practice in engineering when using the bandlimited theory?

when using the bandlimited theory, the common practice in engineering is to force the function f to become bandlimited before sampling. 

Q10. How is the theory of frames used in the case of shift-invariant spaces?

This is accomplished by bringing together wavelet theory, frame theory, reproducing kernel Hilbert spaces, approximation theory, amalgam spaces, and sampling. 

Q11. What is the weighted sequence space p?

The authors also consider the weighted sequence spaces `pν(Zd) with weight ν: a sequence {(ck) : k ∈ Zd} belongs to `pν if ((cν)k) = (ckνk) belongs to `p with norm ‖c‖`pν = ‖νc‖`p , where (νk) is the restriction of ν to Zd.2.2.