. 69
( 83 .)


= [0; ]; = exp ; = ∞:

Lemma 8.1. 1: If S is positive; then there exists a matrix of positive Borel measures with common
support on :
k; j
d = (d ); k = 1; : : : ; p; j = 1; : : : ; q
such that
k; j
k; j
k; j
z f (z) = :
2: The matrix is bounded; if and only if; the measure d has a compact support. In this case;
the measure is uniquely deÿned.
3: If matrix A is bounded; then the generalization of Markov theorem holds: the sequence of
approximants n converges to F on compact sets of the domain D = C ’ ; for some constant .
k; j
Proof. In Section 6, we have proved that for each pair of indices k; j, the rational functions n
can be decomposed into elementary elements:
k; j
d n (x)
k; j
k; j
z n (z) = ; (21)
where d n j are discrete positive measures with common support in the set of the zeros of the

polynomial qn . In matrix form
d n (x)
z n (z) = :
It means that d n is the solution of the ÿnite moment problem

x(p+q) d n (x) = S ; = 0; : : : ; Nn

with Nn ’ ∞. By the ÿrst Helly theorem, it is possible to choose a weak converging subsequence,
i.e., (d n )n∈ .
By the second Helly theorem, the limit measure d is a solution of the full moment problem

x(p+q) d (x) = S ; ¿0;

that is the same as
d (x)
z F(z) = :
Let us now prove the second point of the lemma. From the genetic sums, it follows that A is
bounded if and only if the moments Sn have a geometric estimation. This fact is, as well known,
294 V. Sorokin, J. Van Iseghem / Journal of Computational and Applied Mathematics 122 (2000) 275“295

equivalent to the compactness of the spectrum. The uniqueness of the solution follows from the
uniqueness theorem for holomorphic functions.
For the last point of the theorem: let be the minimal radius of the disk including the support
of d . The proof of the Markov theorem is standard. From (21), follows the uniform boundness of
the sequence ( n ) on compact sets of the domain D. By Montel™s theorem, the sequence ( n ) is
compact. Because it converges to F for the archimedian norm, and ∞ is in D, F is the unique limit
point of ( n ). Hence ( n ) converges to F on compact sets of D.

Remark 8.2. The measure d is the spectral measure of the operator A. If we know d , we can
construct the operator A by decomposition of F into a continued fraction. Note that the arbitrary
measure d does not have positive moments. For example, it is necessary that all elements of the
matrix d have a common support. But this condition is not su cient. There exist some su cients
conditions and examples.

Remark 8.3. The support of d is only a subset of the spectrum of the operator A, and the constant
is less than the spectrum radius of the operator.

Let us consider the case p = 3; q = 2. By initial data bn+3 (0) we construct the operator A and solve
the direct spectral problem, i.e., we look for the measure d by summing the continued fraction.
Consider the special dynamics of the spectral measure
d ˜ (x; t)
˜ 1; 1 d ˜ 1; 1 (x; t);
d ˜ (x; t) = exp(x5 t) d (x); S 0 (t) = d (x; t) = ;
˜ 1; 1
S 0 (t)

d( ˜ ) = x5 d ˜ :
Hence the power moments S n of the measure d ˜ satisfy the following di erential equation:
˜ ˜
(S n ) = S n+1 :

˜ ˜ 1; 1
Because the power moments of the measure d are Sn = S n = S 0 , then

˜ 1; 1 ˜ ˜ 1; 1
˜ ˜
Sn (S 0 ) S n+1 Sn S1
Sn = 1; 1 ’ S n 1; 1 = 1; 1 ’ 1; 1 1; 1
˜ ˜ ˜ ˜˜
S0 (S 0 ) 2 S0 S0 S0

= Sn+1 ’ Sn S1 1 :

We have proved that this last equation is equivalent to the dynamical system. Thus, the reconstruction
parameters bn+3 (t) is the inverse spectral problem and its solution is the decomposition of
1 d (x; t)
z z’x
in a continued fraction.
V. Sorokin, J. Van Iseghem / Journal of Computational and Applied Mathematics 122 (2000) 275“295 295


[1] A. Aptekarev, V. Kaliaguine, J. Van Iseghem, Genetic sums representation for the moments of system of Stieltjes
and applications, Constr. Approx., 16 (2000) 487“524.
[2] D. Barrios, G. Lopez, A. Martinez-Finkelshtein, E. Torrano, Finite-dimensional approximations of the resolvent of
an inÿnite band matrix and continued fractions, Sb. Math. 190 (4) (1999) 501“519.
[3] B. Beckermann, Complex jacobi matrices, Pub ANO403, 1999, submitted for publication.
[4] O.I. Bogoyavlensky, Some constructions of integrable dynamical systems, Math. USSR Izv. 31 (1) (1988) 47“75.
[5] O.I. Bogoyavlensky, Integrable dynamical systems associated with the K de V equation, Math. USSR Izv. 31 (3)
(1988) 435“454.
[6] M.G. de Bruin, The interruption phenomenon for generalized continued fractions, Bull. Austral. Math. Soc. 19 (1978)
[7] V. Kaliaguine, Hermite“Padà approximants and spectral analysis of nonsymmetric operators (Engl. transl.), Russian
Acad., Sci. Sb. Math. 82 (1995) 199“216.
[8] J.K. Moser, Three integrable Hamiltonian systems connected with isospectral deformations, Adv. Math. 16 (1975)
[9] E.M. Nikishin, V.N. Sorokin, Rational approximations and orthogonality, Trans. AMS 92 (1991).
[10] V.I. Parusnikov, The Jacobi“Perron algorithm and simultaneous approximation of functions, Math. Sb. 42 (1982)
[11] V.N. Sorokin, Integrable non linear dynamical systems of Langmuir lattice type, Math zametki (11) (1997).
[12] V.N. Sorokin, J. Van Iseghem, Algebraic aspects of matrix orthogonality for vector polynomials, J. Approx. Theory
90 (1997) 97“116.
[13] V.N. Sorokin, J. Van Iseghem, Matrix continued fraction, J. Approx. Theory 96 (1999) 237“257.
[14] J. Van Iseghem, Matrix continued fraction for the resolvent function of the band operator, Pub ANO 390, Lille,
[15] H.S. Wall, Analytic Theory of Continued Fractions, Van Nostrand, Princeton, NJ, 1948.
Journal of Computational and Applied Mathematics 122 (2000) 297“316

Numerical analysis of the non-uniform sampling problem
Thomas Strohmer
Department of Mathematics, University of California, Davis, CA-95616, USA

Received 23 September 1999; received in revised form 26 November 1999

We give an overview of recent developments in the problem of reconstructing a band-limited signal from nonuniform
sampling from a numerical analysis view point. It is shown that the appropriate design of the ÿnite-dimensional model
plays a key role in the numerical solution of the nonuniform sampling problem. In the one approach (often proposed in the
literature) the ÿnite-dimensional model leads to an ill-posed problem even in very simple situations. The other approach that
we consider leads to a well-posed problem that preserves important structural properties of the original inÿnite-dimensional
problem and gives rise to e cient numerical algorithms. Furthermore, a fast multilevel algorithm is presented that can
reconstruct signals of unknown bandwidth from noisy nonuniformly spaced samples. We also discuss the design of e cient
regularization methods for ill-conditioned reconstruction problems. Numerical examples from spectroscopy and exploration
geophysics demonstrate the performance of the proposed methods. c 2000 Elsevier Science B.V. All rights reserved.

MSC: 65T40; 65F22; 42A10; 94A12

Keywords: Nonuniform sampling; Band-limited functions; Frames; Regularization; Signal reconstruction; Multi-level

1. Introduction

The problem of reconstructing a signal f from nonuniformly spaced measurements f(tj ) arises
in areas as diverse as geophysics, medical imaging, communication engineering, and astronomy. A
successful reconstruction of f from its samples f(tj ) requires a priori information about the signal,
otherwise the reconstruction problem is ill-posed. This a priori information can often be obtained from
physical properties of the process generating the signal. In many of the aforementioned applications
the signal can be assumed to be (essentially) band-limited.

The author was supported by NSF DMS grant 9973373.
E-mail address: strohmer@math.ucdavis.edu (T. Strohmer)

0377-0427/00/$ - see front matter c 2000 Elsevier Science B.V. All rights reserved.
PII: S 0 3 7 7 - 0 4 2 7 ( 0 0 ) 0 0 3 6 1 - 7
298 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

Recall that a signal (function) is band-limited with bandwidth if it belongs to the space B ,
given by
B = {f ∈ L2 (R): f(!) = 0 for |!| ¿ }; (1)
where f is the Fourier transform of f deÿned by
ˆ f(t)e’2 i!t
f(!) = dt:

For convenience and without loss of generality we restrict our attention to the case = 1 , since any
other bandwidth can be reduced to this case by a simple dilation. Therefore, we will henceforth use
the symbol B for the space of band-limited signals.
It is now more than 50 years ago that Shannon published his celebrated sampling theorem [35].
His theorem implies that any signal f ∈ B can be reconstructed from its regularly spaced samples
{f(n)}n∈Z by
sin (t ’ n)
f(t) = f(n) : (2)
(t ’ n)

In practice, however, we seldom enjoy the luxury of equally spaced samples. The solution of the
nonuniform sampling problem poses much more di culties, the crucial questions being:
• Under which conditions is a signal f ∈ B uniquely deÿned by its samples {f(tj )}j∈Z ?
• How can f be stably reconstructed from its samples f(tj )?
These questions have led to a vast literature on nonuniform sampling theory with deep mathemati-
cal contributions see [11,25,3,6,15] to mention only a few. There is also no lack of methods claiming
to e ciently reconstruct a function from its samples [42,41,1,14,40,26,15]. These numerical meth-
ods naturally have to operate in a ÿnite-dimensional model, whereas theoretical results are usually
derived for the inÿnite-dimensional space B. From a numerical point of view the “reconstruction”
of a band-limited signal f from a ÿnite number of samples {f(tj )}r amounts to computing an
ˆ at su ciently dense (regularly) spaced grid points in an interval (t1 ; tr ).
approximation to f (or f)
Hence in order to obtain a “complete” solution of the sampling problem following questions have
to be answered:
• Does the approximation computed within the ÿnite-dimensional model actually converge to the
original signal f, when the dimension of the model approaches inÿnity?
• Does the ÿnite-dimensional model give rise to fast and stable numerical algorithms?
These are the questions that we have in mind, when presenting an overview on recent advances
and new results in the nonuniform sampling problem from a numerical analysis view point.
In Section 2 it is demonstrated that the celebrated frame approach does only lead to fast and stable
numerical methods when the ÿnite-dimensional model is carefully designed. The approach usually
proposed in the literature leads to an ill-posed problem even in very simple situations. We discuss
several methods to stabilize the reconstruction algorithm in this case. In Section 3 we derive an
alternative ÿnite-dimensional model, based on trigonometric polynomials. This approach leads to a
well-posed problem that preserves important structural properties of the original inÿnite-dimensional
problem and gives rise to e cient numerical algorithms. Section 4 describes how this approach
T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 299

can be modiÿed in order to reconstruct band-limited signals for the in practice very important case
when the bandwidth of the signal is not known. Furthermore, we present regularization techniques
for ill-conditioned sampling problems. Finally Section 5 contains numerical experiments from spec-
troscopy and geophysics.
Before we proceed we introduce some notation that will be used throughout the paper. If not
otherwise mentioned ||h|| always denotes the L2 (R)-norm (˜2 (Z)-norm) of a function (vector). For
operators (matrices) ||T || is the standard operator (matrix) norm. The condition number of an in-
vertible operator T is deÿned by „(A) = ||A||||A’1 || and the spectrum of T is (T ). I denotes the
identity operator.

1.1. Nonuniform sampling, frames, and numerical algorithms

The concept of frames is an excellent tool to study nonuniform sampling problems [13,2,1,24,15,44].
The frame approach has the advantage that it gives rise to deep theoretical results and also to the
construction of e cient numerical algorithms “ if (and this point is often ignored in the literature)
the ÿnite-dimensional model is properly designed.
Following Du n and Schae er [11], a family {f }j∈Z in a separable Hilbert space H is said to
be a frame for H , if there exist constants (the frame bounds) A; B ¿ 0 such that
A||f||2 6 | f; f |2 6B||f||2 ; ∀f ∈ H : (3)

We deÿne the analysis operator T by
T : f ∈ H ’ Ff = { f; f }j∈Z (4)


. 69
( 83 .)