In our case assume ||b(n; ) ’ b(n) ||6 ||b(n) ||, where b(n; ) denotes a noisy sample. We terminate the

j

(n; )

CG iterations when the iterates (c )k satisfy for the ÿrst time [22]

||b(n) ’ (c(n; ) )k ||6 ||b(n) || (22)

for some ÿxed ¿ 1.

It should be noted that one can construct “academic” examples where this stopping rule does not

prevent CG from diverging, see [22], “most of the time” however it gives satisfactory results. We

refer the reader to [27,22] for a detailed discussion of various stopping criteria.

There is a variety of reasons, besides the ones we have already mentioned, that make the conjugate

gradient method and the nonuniform sampling problem a “perfect couple”. See Sections 3, 4:1, and

4:2 for more details.

By combining the truncated frame approach with the conjugate gradient method (with appropriate

stopping rule) we ÿnally arrive at a reconstruction method that is of some practical relevance.

However, the only existing method at the moment that can handle large scale reconstruction problems

seems to be the one proposed in the next section.

T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 305

3. Trigonometric polynomials and e cient signal reconstruction

In the previous section we have seen that the naive ÿnite-dimensional approach via truncated

frames is not satisfactory, it already leads to severe stability problems in the ideal case of regular

oversampling. In this section we propose a di erent ÿnite-dimensional model, which resembles much

better the structural properties of the sampling problem, as can be seen below.

The idea is simple. In practice, only a ÿnite number of samples {f(tj )}r is given, where without

j=1

loss of generality we assume ’M 6t1 ¡ · · · ¡ tr 6M (otherwise we can always re-normalize the

data). Since no data of f are available from outside this region we focus on a local approximation

of f on [ ’ M; M ]. We extend the sampling set periodically across the boundaries, and identify this

interval with the (properly normalized) torus T. To avoid technical problems at the boundaries in

the sequel we will choose the interval somewhat larger and consider either [ ’ M ’ 1=2; M + 1=2] or

[ ’ N; N ] with N = M + M=(r ’ 1). For theoretical considerations the choice [ ’ M ’ 1=2; M + 1=2]

is more convenient.

Since the dual group of the torus T is Z, periodic band-limited functions on T reduce to trigono-

metric polynomials (of course, technically f does then no longer belong to B since it is no longer in

L2 (R)). This suggests to use trigonometric polynomials as a realistic ÿnite-dimensional model for a

numerical solution of the nonuniform sampling problem. We consider the space PM of trigonometric

polynomials of degree M of the form

M

’1

ak e2 ikt=(2M +1)

p(t) = (2M + 1) : (23)

k=’M

The norm of p ∈ PM is

M

N

2 2

|ak |2 :

||p|| = |p(t)| dt =

’N k=’M

Since the distributional Fourier transform of p is p = (2M + 1)’1 M

ˆ k=’M ak k=(2M +1) we have

11

supp p ⊆{k=(2M + 1); |k|6M } ⊆ [ ’ 2 ; 2 ]. Hence PM is indeed a natural ÿnite-dimensional model

ˆ

for B.

In general the f(tj ) are not the samples of a trigonometric polynomial in PM , moreover the samples

are usually perturbed by noise, hence we may not ÿnd a p ∈ PM such that p(tj ) = bj = f(tj ). We

therefore consider the least-squares problem

r

|p(tj ) ’ bj |2 wj :

min (24)

p∈PM

j=1

Here the wj ¿ 0 are user-deÿned weights, which can be chosen for instance to compensate for

irregularities in the sampling geometry [14].

By increasing M so that r62M + 1 we can certainly ÿnd a trigonometric polynomial that inter-

polates the given data exactly. However in the presence of noise, such a solution is usually rough

and highly oscillating and may poorly resemble the original signal. We will discuss the question of

the optimal choice of M if the original bandwidth is not known and in presence of noisy data in

Section 4.2.

306 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

The following theorem provides an e cient numerical reconstruction algorithm. It is also the key

for the analysis of the relation between the ÿnite-dimensional approximation in PM and the solution

of the original inÿnite-dimensional sampling problem in B.

Theorem 3.1 (and Algorithm [19,14]). Given the sampling points ’M 6t1 ¡ : : : ; tr 6M; samples

{bj }r ; positive weights {wj }r with 2M + 16r.

j=1 j=1

Step 1: Compute the (2M + 1) — (2M + 1) Toeplitz matrix TM with entries

r

1

wj e’2 i(k’l)tj =(2M +1)

(TM )k; l = for |k|; |l|6M (25)

2M + 1 j=1

and yM ∈ C(2M +1) by

r

1

bj wj e’2 iktj =(2M +1)

(yM )k = √ for |k|6M: (26)

2M + 1 j=1

Step 2: Solve the system

TM aM = yM : (27)

Step 3: Then the polynomial pM ∈ PM that solves (24) is given by

M

1

(aM )k e2 ikt=(2M +1)

pM (t) = √ : (28)

2M + 1 k=’M

Numerical Implementation of Theorem/Algorithm 3.1.

Step 1: The entries of TM and yM of Eqs. (25) and (26) can be computed in O(M log M +

r log(1= )) operations (where is the required accuracy) using Beylkin™s unequally spaced FFT

algorithm [4].

Step 2: We solve TM aM = yM by the conjugate gradient (CG) algorithm [18]. The matrix-vector

multiplication in each iteration of CG can be carried out in O(M log M ) operations via FFT [8].

Thus the solution of (27) takes O(kM log M ) operations, where k is the number of iterations.

Step 3: Usually, the signal is reconstructed on regularly space nodes {ui }N . In this case pM (ui ) in

i=1

(28) can be computed by FFT. For non-uniformly spaced nodes ui we can again resort to Beylkin™s

USFFT algorithm.

There exists a large number of fast algorithms for the solution of Toeplitz systems. Probably the

most e cient algorithm in our case is CG. We have already mentioned that the Toeplitz system

(27) can be solved in O(kM log M ) via CG. The number of iterations k depends essentially on the

clustering of the eigenvalues of TM , cf. [8]. It follows from equation (31) below and perturbation

theory [10] that, if the sampling points stem from a perturbed regular sampling set, the eigenvalues

of TM will be clustered around ÿ, where ÿ is the oversampling rate. In such cases we can expect

a very fast rate of convergence. The simple frame iteration [26,1] is not able to take advantage of

such a situation.

For the analysis of the relation between the solution pM of Theorem 3.1 and the solution f of the

original inÿnite-dimensional problem we follow Grochenig [20]. Assume that the samples {f(tj )}j∈Z

of f ∈ B are given. For the ÿnite-dimensional approximation we consider only those samples f(tj )

T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 307

for which tj is contained in the interval [’M ’ 1 ; M + 1 ] and compute the least-squares approximation

2 2

pM with degree M and period 2M + 1 as in Theorem 3.1. It is shown in [20] that if (TM ) ⊆ [ ; ÿ]

for all M with ¿ 0 then

|f(t) ’ pM (t)|2 dt = 0

lim (29)

M ’∞ [’M; M ]

and also limpM (t) = f(t) uniformly on compact sets.

Under the Nyquist condition sup (tj+1 ’tj ):= ¡ 1 and using weights wj =(tj+1 ’tj’1 )=2 Grochenig

has shown that

(TM ) ⊆ [(1 ’ )2 ; 6]; (30)

independently of M , see [20]. These results validate the usage of trigonometric polynomials as

ÿnite-dimensional model for nonuniform sampling.

Example 1 (Reconsidered). Recall that in Example 1 of Section 2 we have considered the recon-

struction of a regularly oversampled signal f ∈ B. What does the reconstruction method of Theorem

3.1 yield in this case? Let us check the entries of the matrix TM when we take only those samples

n

in the interval [ ’ n; n]. The period of the polynomial becomes 2N with N = n + r’1 where r is the

number of given samples. Then

r nm

1 2 i(k’l)tj =(2N )

e2 i(k’l) j=(2nm+1)

(TM )k; l = e = =m (31)

k; l

2N j=’nm

j=1

for k; l = ’M; : : : ; M , where is Kronecker™s symbol with the usual meaning = 1 if k = l and

k; l k; l

0 else. Hence we get

TM = mI;

where I is the identity matrix on C2M +1 , thus TM resembles the structure of the inÿnite-dimensional

frame operator S in this case (including exact approximation of the frame bounds). Recall that the

truncated frame approach leads to an “artiÿcial” ill-posed problem even in such a simple situation.

The advantages of the trigonometric polynomial approach compared to the truncated frame ap-

proach are manifold. In the one case we have to deal with an ill-posed problem which has no

speciÿc structure, hence its solution is numerically very expensive. In the other case we have to

solve a problem with rich mathematical structure, whose stability depends only on the sampling

density, a situation that resembles the original inÿnite-dimensional sampling problem.

In principle, the coe cients aM = {(aM )k }M

k=’M of the polynomial pM that minimizes (24) could

also be computed by directly solving the Vandermonde-type system

W VaM = Wb; (32)

√

where Vj; k = (1=( 2M + 1))e’2 iktj =(2M +1) for j = 1; : : : ; r; k = ’M; : : : ; M and W is a diagonal matrix

√

with entries Wj; j = wj , cf. [31]. Several algorithms are known for a relatively e cient solution

of Vandermonde systems [5,31]. However this is one of the rare cases, where, instead of directly

solving (32), it is advisable to explicitly establish the system of normal equations

TM aM = yM ; (33)

where T = V — W 2 V and y = V — W 2 b.

308 T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316

The advantages of considering the system TM aM = yM instead of the Vandermonde system (32)

are manifold:

• The matrix TM plays a key role in the analysis of the relation of the solution of (24) and the

solution of the inÿnite-dimensional sampling problem (9), see (29) and (30) above.

• TM is of size (2M + 1) — (2M + 1), independently of the number of sampling points. Moreover,

since (TM )k; l = r wj e2 i(k’l)tj , it is of Toeplitz type. These facts give rise to fast and robust

j=1

reconstruction algorithms.

• The resulting reconstruction algorithms can be easily generalized to higher dimensions, see

Section 3.1. Such a generalization to higher dimensions seems not to be straightforward for

fast solvers of Vandermonde systems such as the algorithm proposed in [31].

We point out that other ÿnite-dimensional approaches are proposed in [16,7]. These approaches

may provide interesting alternatives in the few cases where the algorithm outlined in Section 3 does

not lead to good results. These cases occur when only a few samples of the signal f are given in

an interval [a; b] say, and at the same time we have |f(a) ’ f(b)|0 and |f (a) ’ f (b)|0, i.e.,

if f is “strongly nonperiodic” on [a; b]. However the computational complexity of the methods in

[16,7] is signiÿcantly larger.

3.1. Multi-dimensional nonuniform sampling

The approach presented above can be easily generalized to higher dimensions by a diligent

book-keeping of the notation. We consider the space of d-dimensional trigonometric polynomials

PM as ÿnite-dimensional model for B d . For given samples f(tj ) of f ∈ B d , where tj ∈ Rd , we

d

compute the least-squares approximation pM similar to Theorem 3.1 by solving the corresponding

system of equations TM aM = yM .

In 2-D for instance the matrix TM becomes a block Toeplitz matrix with Toeplitz blocks [37]. For

a fast computation of the entries of T we can again make use of Beylkin™s USFFT algorithm [4].

And similar to 1-D, multiplication of a vector by TM can be carried out by 2-D FFT.

d

Also the relation between the ÿnite-dimensional approximation in PM and the inÿnite-dimensional

solution in B d is similar as in 1-D. The only mathematical di culty is to give conditions under which

the matrix TM is invertible. Since the fundamental theorem of algebra does not hold in dimensions

larger than one, the condition (2M + 1)d 6r is necessary but no longer su cient for the invertibility

of TM . Su cient conditions for the invertibility, depending on the sampling density, are presented

in [21].

4. Bandwidth estimation and regularization

In this section we discuss several numerical aspects of nonuniform sampling that are very important

from a practical viewpoint, however only few answers to these problems can be found in the

literature.

T. Strohmer / Journal of Computational and Applied Mathematics 122 (2000) 297“316 309

4.1. A multilevel signal reconstruction algorithm

In almost all theoretical results and numerical algorithms for reconstructing a band-limited signal

from nonuniform samples it is assumed that the bandwidth is known a priori. This information

however is often not available in practice.

A good choice of the bandwidth for the reconstruction algorithm becomes crucial in case of noisy

data. It is intuitively clear that choosing a too large bandwidth leads to over-ÿt of the noise in the

data, while a too small bandwidth yields a smooth solution but also to under-ÿt of the data. And

of course we want to avoid the determination of the “correct” by trial-and-error methods. Hence

the problem is to design a method that can reconstruct a signal from nonuniformly spaced, noisy

samples without requiring a priori information about the bandwidth of the signal.

The multilevel approach derived in [34] provides an answer to this problem. The approach applies

to an inÿnite-dimensional as well as to a ÿnite-dimensional setting. We describe the method directly

for the trigonometric polynomial model, where the determination of the bandwidth translates into

the determination of the polynomial degree M of the reconstruction. The idea of the multilevel

algorithm is as follows.

Let the noisy samples {bj }r ={f (tj )}r of f ∈ B be given with r |f(tj )’b (tj )|2 6 2 ||b ||2

j=1 j=1 j=1

and let QM denote the orthogonal projection from B into PM . We start with initial degree M = 1

and run Algorithm 3:1 until the iterates p0; k satisfy for the ÿrst time the inner stopping criterion

r

|p1; k (tj ) ’ bj |2 62 ( ||b || + ||Q0 f ’ f||)||b ||

j=1