. 71
( 83 .)


which remains to be controlled by an appropriate stopping rule [27,22].
In our case assume ||b(n; ) ’ b(n) ||6 ||b(n) ||, where b(n; ) denotes a noisy sample. We terminate the
(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
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
ak e2 ikt=(2M +1)
p(t) = (2M + 1) : (23)

The norm of p ∈ PM is
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
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
|p(tj ) ’ bj |2 wj :
min (24)

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
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
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
(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
(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
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

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

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.
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
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
|p1; k (tj ) ’ bj |2 62 ( ||b || + ||Q0 f ’ f||)||b ||


. 71
( 83 .)