p+q

k=1

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

d

(x)

k; j

k; j

z f (z) = :

z’x

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)

z’x

where d n j are discrete positive measures with common support in the set of the zeros of the

k;

polynomial qn . In matrix form

d n (x)

z n (z) = :

z’x

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) = :

z’x

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)

then

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 :

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

References

[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)

245“272.

[7] V. Kaliaguine, Hermite“PadÃ approximants and spectral analysis of nonsymmetric operators (Engl. transl.), Russian

e

Acad., Sci. Sb. Math. 82 (1995) 199“216.

[8] J.K. Moser, Three integrable Hamiltonian systems connected with isospectral deformations, Adv. Math. 16 (1975)

197“210.

[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)

287“296.

[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,

1998.

[15] H.S. Wall, Analytic Theory of Continued Fractions, Van Nostrand, Princeton, NJ, 1948.

Journal of Computational and Applied Mathematics 122 (2000) 297“316

www.elsevier.nl/locate/cam

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

Abstract

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

method

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

2

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)

n∈Z

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

j=1

ˆ 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

j

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)

j

j

We deÿne the analysis operator T by

T : f ∈ H ’ Ff = { f; f }j∈Z (4)

j