SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

arXiv:1709.02086v1 [cs.IT] 7 Sep 2017

CAN EVREN YARMAN

Abstract. R-limited functions are multivariate generalization of band-limited
functions whose Fourier transforms are supported within a compact region
R ⊂ Rn . In this work, we generalize sampling and interpolation theorems
for band-limited functions to R-limited functions. More precisely, we investigated the following question: “For a function compactly supported within a
region similar to R, does there exist an R-limited function that agrees with the
function over its support for a desired accuracy?”. Starting with the Fourier
domain definition of an R-limited function, we write the equivalent convolution and a discrete Fourier transform representations for R-limited functions
through approximation of the convolution kernel using a discrete subset of
Fourier basis. The accuracy of the approximation of the convolution kernel
determines the accuracy of the discrete Fourier representation. Construction
of the discretization can be achieved using the tools from approximation theory
as demonstrated in the appendices. The main contribution of this work is proving the equivalence between the discretization of the Fourier and convolution
representations of R-limited functions. Here discrete convolution representation is restricted to shifts over a compactly supported region similar to R. We
show that discrete shifts for the convolution representation are equivalent to
the spectral parameters used in discretization of the Fourier representation of
the convolution kernel. This result is a generalization of the cardinal theorem
of interpolation of band-limited functions. The error corresponding to discrete convolution representation is also bounded by the approximation of the
convolution kernel using discretized Fourier basis.

1. Introduction
R-limited functions are functions whose Fourier transforms are supported within
a region R ⊂ Rn . They are multivariate generalization of band-limited functions.
The terminology was coined by Slepian in [22]. In this work, we generalize sampling and interpolation theorems for band-limited functions to R-limited functions.
Specifically, we explore answers to the following questions: “For a function compactly supported within a region similar to R, does there exist an R-limited function
that agrees with the function over its support? If so, how shall we sample the function to construct an R-limited function that approximates the original function
within a desired accuracy?”. The first question has been answered in [22]. Combining these results with methods from approximation theory, we answer the second
question. Answering how to sample also provides a guide for where to sample.
In our exposition we choose an approximation theory perspective which provides
an alternative insight to understanding of band-limited functions through discretization of the sine cardinal function as well as a deterministic framework for constructing sampling schemes for R-limited functions. Starting with the Fourier domain
definition of a R-limited function, we write the equivalent convolution representation and write a discrete Fourier transform representation for R-limited functions
1

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

2

through approximation of the convolution kernel using a discrete subset of Fourier
basis. The accuracy of the approximation of the kernel determines the accuracy
of the discrete Fourier representation. Construction of the discretization can be
achieved using the tools from approximation theory such as generalization of Padé
approximation, which is summarized in Appendix A.
Our main results, Theorems 10 and 12, prove the equivalence between the discretization of the Fourier and convolution representations that approximate a compactly supported function within a region similar to R. We show that discrete shifts
for the convolution representation are equivalent to the spectral parameters used in
discretization of the Fourier representation of the convolution kernel. Discretization
of the convolution representation is also referred to as sampling and interpolation
theorem.
Theorem 12 is a generalization of the sampling and interpolation theorem for
band-limited functions summarized in Theorem 7. It also provides a way to analyze and approximate the resulting error. We show that the error corresponding to
discrete convolution representation is bounded by error obtained from discretization
of the Fourier transform of the convolution kernel. In single dimension, it provides
a new way to prove truncation of discrete Fourier series as well as sinc interpolation
formula. Furthermore, in single dimension, our result indicates that for a support
of interest, instead of uniform sampling, improvement in discrete the Fourier representation of band-limited functions can be obtained using Gauss-Legendre type
quadratures (see Appendix A). This also raises the questions on what is the most
cost efficient way to implement fast Fourier transforms using Gaussian quadratures
which may be addressed using the ideas from [6, 5, 7, 12] and left for a future
discussion. Similar discretizations are obtained for special cases of R-limited functions in Appendix C and D where we make use of cascaded quadratures that are
equivalent to Gauss-Legendre or Clenshaw-Curtis quadratures. While the body of
the manuscript contains our main results, the Appendices also provide as valuable
information by providing a constructive way for computing quadratures to discretize convolution kernels which can be utilized in the sampling and interpolation
theorems.
The outline of the manuscript is as follows. In Section 2, we present the conventions used in the rest of our discussion. To motivate the multivariate case,
in Section 3, we study discretization of Fourier transform and sinc interpolation
formula for one-dimensional (univariate) band-limited functions, or band-limited
projection of compactly supported functions. Both discrete Fourier transform and
sinc interpolation formulas have been studied in the literature with many books
devoted to this topic. We refer the reader to [13, 27, 26, 14] for a comprehensive
list of references on these topics. We give an alternative exposition, which leads
to proof of the equivalence of sampling in the domain of the function (Theorem
7) and its Fourier transform (Theorem 1). The necessary background material
for Section 3 is provided in Appendices A and B which discuss Generalization of
Pade approximation and approximations to sine cardinal function. Compared to
band-limited function, sampling and representation of multivariate functions whose
Fourier transforms’ support are not similar to a hypercube is studied and understood less. In Section 4, we extend our results for band-limited function (Theorems
1 and 7) to R-limited functions (Theorems 10 and 12). Examples of special cases
of convolution kernels for R-limited functions are presented in Appendices C and

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

3

D. In Appendix C, we provide a method to construct quadratures for isosceles triangle and trirectangular tetrahedron which are used to construct quadratures for
equilateral triangle and regular tetrahedron. Similar method is used in Appendix
D to construct quadratures for a finite cone and a ball in three dimensions which
have practical importance in multidimensional signal processing seismic data, image
processing and video processing.
2. Conventions
We employ the following conventions of Fourier transform, inverse Fourier transform and convolution.
The Fourier transform F [f ] (k) of f (x), an absolutely integrable function for
x, k ∈ RN , which we denote by fˆ (k), is defined by
Z
ˆ
(2.1)
F [f ] (k) = f (k) =
f (x) e−i2πk·x dx
RN

The inverse Fourier transform is defined by
Z
h i
(2.2)
F −1 fˆ (x) = f (x) =

fˆ (k) ei2πk·x dk

RN

Denoting the convolution operator by ∗, convolution of two functions is defined
by
Z
(2.3)

(f ∗ g) (x) =

f (y) g (x − y) dy
RN

The Fourier transform of the convolutions is the product of the Fourier transforms:
(2.4)

F [f ∗ g] (k) = fˆ (k) ĝ (k)

also referred to as convolution theorem.
3. Band-limited functions
We say that fB (t), for t ∈ R, is a band-limited function with band-limit B if
there exists an fˆB (ω) such that
Z B
(3.1)
fB (t) =
fˆB (ω) ei2πtω dω
−B

Z

1

fˆB (Bω) ei2πBtω dω

=B
−1

Given a function f (t), its band-limited projection PB [f ] (t), denoted by fB (t)
for short, is defined by
Z B
PB [f ] (t) = fB (t) =
fˆ (ω) ei2πtω dω
−B

Z
(3.2)

1

=B

fˆ (Bω) ei2πBtω dω

−1

or, equivalently, in the convolution representation using Parseval’s theorem
Z ∞
(3.3)
PB [f ] (t) = fB (t) =
f (τ ) 2B sinc (2πB [t − τ ]) dτ
−∞

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

4

where
(3.4)

sinc (Bt) =

1
2B

Z

B

eiωt dω =

−B

Z

1

cos (Bωt) dω =
0

sin(Bt)
Bt

is the sinc function normalized with band-limit B. Note that
(3.5)
fˆB (ω) = fˆ (ω) , ω ∈ [−B, B]
3.1. Discrete Fourier representation of band-limited approximation of
compactly supported functions: In this section we derive discrete Fourier approximations of band-limited projection of compactly supported functions starting
from their convolution representation.
Consider a discretization of the integral representation of sinc (see Figures B.1
and B.3 for two examples. Another example in Section 8 of [1].)
Z

B

2B sinc (Bt) = 2

cos (ωt) dω
0

=

(3.6)

M
X

2αm cos (Bωm t) + α0 + B (t)

m=1

for a given B ∈ R+ , with αm ∈ R+ and ωm ∈ [0, 1]. Equivalently, using exponentials
instead of cosines, we write
Z B
M
X
(3.7)
2B sinc (Bt) =
eitω dω =
αm eiBωm t + B (t)
−B

m=−M

where −ωm = ω−m and α−m = αm . Now we can prove:
Theorem 1. Given a compactly supported function f (t) over [−T, T ], restriction
of its band-limited projection, fB (t), onto interval [−T, T ] can be approximated as
a discrete sum of Fourier basis by
(3.8)

M
X

fB (t) =

αm fˆ (Bωm ) ei2πBωm t + f (t)

m=−M

using the approximation (3.7) with the error bound
(3.9)

max |f (t)| ≤ 2T

t∈[−T,T ]

max |f (t)|
t∈[−T,T ]

max
t∈[−2T,2T ]

|B (2πt)| .

Proof. For a function f (τ ) compactly supported on τ ∈ [−T, T ], substituting (3.7)
into (3.3), its band-limited projection can be approximated by (3.8) where
Z T
(3.10)
f (t) =
f (τ ) B (2π [t − τ ]) dτ
−T

(3.8) provides a discretization of (3.2) through approximation of the sinc function
as a sum of cosines.


Example 2. Choosing
(3.11)

M

(αm , ωm )m=−M =



2m
2B
,
2M + 1 2M + 1

M
m=−M

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

5

for some M ≥ 0, (3.8) becomes the discrete Fourier transform representation of
fB (t):


M
X
2m
2B
2m
ˆ
fB (t) =
(3.12)
f B
ei2πB 2M +1 t + f (t)
2M + 1
M +1
m=−M

The summation term is referred to as the discrete inverse Fourier transform of fˆ.
(3.12) is a Riemann sum approximation of the integral (3.1) for uniform sampling
of the interval [−B, B].
In practice measurements are performed over a finite duration. Thus it is desirable
to have the band-limited projection of a compactly supported function approximately agree with the function at least over t ∈ [−T, T ]. Thus, by (3.9), for fB (t)
to approximate accurately f (t) over [−T, T ], one needs to build up an approximation to sinc (Bt) that is accurate over the interval [−4πT, 4πT ]. In Appendix B,
we present two different approximations in the form of (3.7) (see Figures B.1 and
B.3), one using Gauss-Legendre quadratures (see Figure B.2) and the other using
uniform sampling. We show that, for a desired interval and bandwidth, a discrete
representation of sinc that is accurate upto machine precision can be achieved using
Gauss-Legendre quadratures without requiring as many uniform samples.
Example 3. From a finite duration measurement, only finite number of samples are
utilized for digital signal processing. This raises a natural question: “What should
be the sampling rate for a band-limited measurement such that band-limited projection of the the discrete measurement agree with with the discrete measurement?”.
In this regard, consider the following model for a discrete measurement


K
X
2k
(3.13)
T
f (t) =
fk δ t −
2K + 1
k=−K

where 2T /(2K + 1) is the sampling period. Then
fˆ (ω) =

(3.14)

K
X

2k

fk e−i2πω 2K+1 T

k=−K

Assuming that the measurement has band-limit B, let ωm = B 2m (2M + 1)−1 . By
(3.14), we rewrite (3.12) in terms of fˆ (ωm ) and obtain
!
K
M
X
X
2m
2k
2B
i2πB 2M
t−
T
(
)
+1
2K+1
fk
e
(3.15) fB (t) =
2M + 1
k=−K

m=−M

+ f (t)
which, for t = 2l/ (2K + 1) T , l = −K, . . . , K, becomes

(3.16) fB

2l
T
2K + 1



K
X
2B
=
fk
2M + 1
k=−K

M
X

!
e

2m
i2πB 2M
+1

2T
2K+1 (l−k)

m=−M


+ f

2l
T
2K + 1

K

(3.14) is known as the discrete Fourier transform of the vector {fk }k=−K .



SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

6

For T = (2K + 1) /(4B), (3.13), (3.14) for ω = B 2m (2M + 1)−1 , and (3.16)
become

PK
k
f (t) = k=−K fk δ t − 2B
(3.17)

fˆ B

(3.18)


(3.19)

fB

l
2B

2m
2M + 1


=

k=−K

K
X
2B
=
fk
2M + 1



m

PK

k=−K

fk e−i2π 2M +1 k

!

M
X

i2π 2Mm+1 (l−k)

e


+ f

m=−M

l
2B



which, for |l − k| ≤ 2M , is

(3.20)

fB

l
2B


= 2B

K
X


fk δkl + f

k=−K

l
2B




= 2Bfl + f

l
2B


,

where δkl is the Kronecker delta function. 1 The inequality |l − k| ≤ 2M imposes
that M ≥ K. For T = (2K + 1)/(4B), by (3.13) and (3.7), (3.10) becomes

 Z T

 
l
l
f
=
f (τ ) B 2π t −
dτ
2B
2B
−T

 
K
X
l
k
fk B 2π
−
=
2B
2B
k=−K
"
#
K
X
2Bsinc
(π [k − l])


PM
=
fk
i2π 2Mm+1 (l−k)
2B
− 2M
m=−M e
+1
k=−K

=

K
X


fk 2Bδk,l −

k=−K

(3.21)

2B
(2M + 1) δk,l
2M + 1



= 0.

The special sampling rate 1/(2B) that gave rise to the band-limited function fB (t)
whose values are equal to the original function at t = l/(2B). This sampling rate
is referred to as the Nyquist rate. By sinc interpolation, also known as WhittakerShannon interpolation formula [26],





∞
X
k
k
(3.22)
fB (t) =
fB
sinc 2πB t −
,
2B
2B
k=−∞

1Consider (3.13) for t = lT (K + 1)−1 :


 P
K+1
T
f l K+1
= K
k=−K fk T δ (l − k) .
For T = (K + 1)/(2B), f



l
2B



= 2B

PK

k=−K

fk δ (l − k).

Thus the factor 2B in front of

the sum in (3.20)
 is due to difference between continious and discrete nature of Dirac delta,
δ [l − k](2B)−1 = 2B δ (l − k), and Kronecker delta, δlk , functions.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

7

a band-limited function can be exactly determined from samples obtained using
Nyquist rate. In digital signal processing, because finitely many samples are measured, which are modelled by (3.17), (3.22) is approximated by





K
X
k
k
fB (t) ≈
(3.23)
f
sinc 2πB t −
,
2B
2B
k=−K

which is band-limited projection of (3.13). This approximation agrees with the
discrete measurements. However, there is an approximation error between the
sample locations as a result of implicit imposition f k(2B)−1 equal to zero for
|k| > K. This imposition is eliminated when, instead of discrete Fourier basis,
prolate spheroidal wave functions (PSWFs) are used as a basis to represent bandlimited functions. Expressing band-limited function in terms of PSWFs doesn’t
directly answer how a band-limited function should be sampled but provides the
necessary foundation to answer “How accurately can we approximate a band-limited
function from its samples given over a compact support?”, which is addressed in
Theorem (7).
3.2. Band-limited projections of compactly supported function and prolate spheroidal wave functions. In this section we present prolate spheroidal
wave functions, their properties and two methods on how we can numerically compute them. For the rest of our discussion we will assume that T = 1. This can be
compensated by choosing the band-limit to be T times more.
3.2.1. Prolate spheroidal wave functions (PSWF). Prolate spheroidal wave functions (PSWF), ϕn (t) can be defined as the eigenfunctions of the band-limited projection operator restricted to a compact support, which, without loss of generality,
is given by [24, 20]
Z 1
(3.24)
2B sinc (2πB (t − τ )) ϕn (τ ) dτ = µn ϕn (t) , t ∈ R
−1

PSWF form a basis for band-limited functions as well as L2 ([−1, 1]), and satisfy
the following properties [24, 20]:
(1) PSWF are real valued and corresponding eigenvalues µn are positive: ϕn (t) ∈
R, µn ∈ R+ .
(2) PSWF are orthogonal within the interval t ∈ [−1, 1] as well as over the real
line:
Z ∞
Z 1
(3.25)
µn
ϕn (t) ϕm (t) dt =
ϕn (t) ϕm (t) dt = δm,n
−∞

−1

where δm,n is the Kronocker delta function equal to 1 for m = n and zero
otherwise.
(3) PSWF are eigenfunctions of Fourier operator restricted to the interval
[−1, 1]:
Z 1
(3.26)
ϕn (t) ei2πBωt dt = λn ϕn (ω) , ω ∈ [−1, 1] , λn ∈ C
−1

The second property implies that if a band-limited function is known within an
interval then it is known
The eigenvalues satisfy the following properties:

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

8

(1) Multiplying both sides of (3.26) with e−i2πBωt , integrating with respect to
2
ω and comparing the result with (3.24) one obtains µn = B |λn | (see 3.48
in [20])
(2) [[17],Theorem 3.14 in [20]] Let N (α) denote the number of eigenvalues
µn > α for some 0 < α < 1. Then



1−α
1
log
(3.27)
N (α) = 4B +
log (2πB) + O (log (2πB))
π2
α
Thus there are about 4B eigenvalues µn that are close to one, on the order
of log (2πB) eigenvalues that decay rapidly, and the rest of them are very
close to zero.
For a comprehensive review on PSWF, we refer the reader to [20].
Lemma 4. Given a compactly supported function f (t) on t ∈ [−1, 1], it can be
expressed in terms of PSWFs by
X
f (t) =
(3.28)
fB,n ϕn (t)
n

where
Z

1

fB,n =

f (t) ϕn (t) dt
−1

Z 1
1
fB (t) ϕn (t) dt
µn −1
Z ∞
=
fB (t) ϕn (t) dt
=

(3.29)

−∞

Proof. Because 2B sinc (2πB (t − τ )) is band-limited, it can be expanded as a sum
of PSWF,
X
(3.30)
2B sinc (2πB (t − τ )) =
an (τ ) ϕn (t) ,
n

where
Z
(3.31)

1

2B sinc (2πB (t − τ )) ϕn (t) dt = µn ϕn (τ )

an (τ ) =
−1

leading to the decomposition of sinc in terms of PSWF:
X
(3.32)
2B sinc (2πB (t − τ )) =
µn ϕn (τ ) ϕn (t)
n

Given a compactly supported function f (t) on t ∈ [−1, 1], it’s band-limited projection fB (t) can be expanded in term of PSWF by substituting (3.32) in (3.3)
X
(3.33)
fB (t) =
µn fB,n ϕn (t) .
n

Then, by (3.25), the coefficients can be computed by either of the three ways in
(3.29).


SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

9

3.2.2. Approximating PSWF as eigenvectors of e−i2πBωk ωm . Substituting (3.7) into
2
(3.24), and recalling µn = B |λn | , by (3.26), we obtain
Z 1
1
ϕn (t) =
2B sinc (2πB (t − τ )) ϕn (τ ) dτ
µn −1
(3.34)

=

1
Bλn

M
X

αm ei2πBωm t ϕn (ωm ) + ϕn (t)

m=−M

where
(3.35)
with

ϕn (t) =

1
µn

Z

1

B (2π [t − τ ]) ϕn (τ ) dτ
−1

2

max |ϕn (t)| ≤

(3.36)

t∈[−1,1]

1
µn

max |B (2πt)|

t∈[−2,2]

In [1] (see equation (8.19) and (8.20) in [1]), (3.34) was used to build approximate
PSWFs by first solving the eigensystem
(3.37)

M
X

1
ϕn (ωk ) =
Bλn

αm ei2πBωm ωk ϕn (ωm )

m=−M

for the eigenvector ϕn (ωk ), followed by substituting ϕn (ωk ) back in (3.34):
(3.38)

ϕn (t) =

1
Bλn

M
X

αm ei2πBωm t ϕn (ωm ) + ϕn (t)

m=−M

where ϕn (ωk ) = 0, for k = −M, . . . , M .
Thus, the eigensystem (3.37) provides an approximation to PSWF over the interval [−1, 1] bounded by (3.36). Because (3.37) is a system of 2M + 1 equations,
it has 2M + 1 eigenvalues, which we will denote by µn=0,...,2M . By (3.27), in order to capture the dominant eigenvalues, i.e. eigenvalues around 1, one shall have
M ≥ d2B − 1/2e.
For sake of simplicity of the discussion, we will assume that M  d2B − 1/2e,
µ2M  1. Thus the corresponding 2M + 1 approximate PSWFs provides a sufficiently accurately approximate band-limited functions over the interval [−1, 1] and,
by the same token, sinc function over [−2, 2]. Thus we treat, (3.32) is equivalent to
its truncated version:
M
X
(3.39)
2B sinc (2πB (t − τ )) ≈
µn ϕn (τ ) ϕn (t) , t, τ ∈ [−1, 1] ,
n=−M

and similarly all the infinite sums over PSWFs as finite sums.
2By Hölder’s inequality,

|ϕn (t)|2 ≤

1
µ2n

≤

1
µ2n

Z

1

|B (2π [t − τ ])|2 |ϕn (τ )|2 dτ

−1

max

τ ∈[t−1,t+1]

|B (2πτ )|2

Z

1
−1

|

|ϕn (τ )|2 dτ
{z
}
=1

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

10

Example 5. Consider (3.26) and the quadratures of discrete inverse Fourier transform (3.11) for B = 2M4+1

M
1
2m
M
(αm , ωm )m=−M =
,
(3.40)
2 2M + 1 m=−M
for some positive integer M . Then (3.34) becomes


M
X
2
2m
1
(3.41)
eiπmt ϕn
+ ϕn (t)
ϕn (t) =
λn 2M + 1
2M + 1
m=−M

Consequently, the eigensystem for approximate PSWF is




M
X
mk
1
2
2k
2m
(3.42)
=
ei2π 2M +1 ϕn
ϕn
2M + 1
λn 2M + 1
2M + 1
m=−M
√
for k = √
−M, . . . , M , which for M ≥ 2 has four distinct eigenvalues, ±2 2M + 1
and ±i2 2M + 1 with multiplicities (see page 32 of [2]).
3.2.3. Approximating PSWF as eigenvectors of sinc (2πB (ωm − ωk )). An alternative to the method presented in Section 3.2.2, PSWF can be approximated through
discretization of (3.24).
Starting with (3.34) and using (3.26), we have
Z 1
1
ϕn (ω) =
(3.43)
e−i2πBωt ϕn (t) dt
λn −1
M
X
1
=
αm 2B sinc (2πB (ωm − ω)) ϕn (ωm ) + εϕn (ω)
(3.44)
Bµn
m=−M

where
εϕn (ω) =

(3.45)

1
Bλn

Z

1

e−i2πBωt ϕn (t) dt

−1

2

Because µn = B |λn | , by (3.36),
2
max |ϕ (t)|
|λn | t∈[−1,1] n
√
2 B
≤ 3/2 max |B (2πt)|
µn t∈[−2,2]

max |εϕn (t)| ≤

(3.46)

t∈[−1,1]

(3.47)

Similar to the method of [1], (3.43) can be used to build approximate PSWF by
first solving the eigensystem
(3.48)

ϕn (ωm ) =

M
X
1
αk 2B sinc (2πB (ωm − ωk )) ϕn (ωk )
Bµn
k=−M

for the eigenvector ϕn (ωk ), followed by substituting ϕn (ωk ) back in (3.43):
(3.49)

ϕn (t) =

M
X
1
αk 2B sinc (2πB (t − ωk )) ϕn (ωk ) + εϕn (t) ,
Bµn
k=−M

where εϕn (ωm ) = 0, for m = −M, . . . , M . The eigenvectors φn (ωk ) are generalizations of discrete prolate spheroidal sequences (DPSS) [23]. When ωk are uniformly

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

11

sampled they are equivalent to DPSS, which asymptotically approximate PSWF
[23]. Similar to the discussion in Section 3.2.2, we say the eigensystem (3.48) provides an approximation to PSWF. Because it can only capture 2M + 1 of the eigenvalues, which we denote by µn=0,...,2M , by (3.27), one shall choose M ≥ d2B − 1/2e
in order to capture all the eigenvalues close to one and some of the eigenvalues in
the transition zone from one to zero, depending on the desired accuracy of the
approximation.
Example 6. Consider the quadratures of discrete inverse Fourier transform for
B = 2M4+1


M

(αm , ωm )m=−M =

(3.50)

1
2m
,
2 2M + 1

M
m=−M

for some positive integer M . Then (3.43) becomes
(3.51)

ϕn (t) =

 



M
2M + 1
1 X
2k
sinc π
t − k ϕn
+ εϕn (t)
µn
2
2M + 1
k=−M

Consequently, the eigensystem for approximate PSWF is

ϕn

2m
2M + 1


=



M
1 X
2k
sinc (π [m − k]) ϕn
µn
2M + 1
k=−M



M
1 X
2k
δm,k ϕn
µn
2M + 1
k=−M


1
2m
=
ϕn
µn
2M + 1
=

(3.52)

for m = −M,
 . . ., M which implies that µn = 1 for n = 0, . . . , 2M . As mentioned
above, ϕn M2m
+1 are referred to as discrete prolate spheroidal sequences and were
studied in [23] along with their relationship to periodic discrete prolate spheroidal
sequences (P-DPSS). This example shows that, similar to P-DPSS [28], eigenvalues
of DPSS are not necessarily simple and therefore definition of DPSS can be non
unique.
3.3. Discrete convolution representation of band-limited approximation
of compactly supported functions.
Theorem 7. Consider a function f (t) compactly supported on t ∈ [−1, 1]. Its
band-limited projection fB (t) can be computed by
(3.53)

fB (t) =

M
X

2B sinc (2πB (t − ωk )) fk + fB (t)

k=−M

where
(3.54)

fk =

M
X
m=−M

f (ωm ) αm Rm (ωk )

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

12

and (αm , ωm ) satisfy (3.7) and
Rm (t) =

(3.55)

2M
X

µ−1
n ϕn (ωm ) ϕn (t)

n=0

Proof. Substituting (3.43) in (3.33), we obtain
(3.56)

M
X

fB (t) =

αk 2B sinc (2πB (t − ωk )) f˜B (ωk ) + fB (t)

k=−M

where
2M
1 X
fB,n ϕn (t)
f˜B (t) =
B n=0

(3.57)

fB (t) =

(3.58)

2M
X

µn fB,n εϕn (t)

n=0

with fB,n defined in (3.29). Similarly, substituting (3.43) in (3.29), we obtain
Z 1
M
X
1
1
fB,n =
αk ϕn (ωk ) fB (ωk ) +
(3.59)
fB (t) εϕn (t) dt,
Bµn
µn −1
k=−M

and substituting (3.59) in (3.57), we obtain
(3.60)
1
f˜B (t) = 2
B

M
X

αm fB (ωm ) Rm (t) +

2M
X

Z


fB (τ ) εϕn (τ ) dτ

.

−∞

n=0

m=−M

∞

ϕn (t)


Corollary 8. The error fB (t), for t ∈ [−1, 1] is bounded by
(3.61)

max |fB (t)| ≤ C max |fB (t)| max |B (2πt)|

t∈[−1,1]

t∈[−1,1]

t∈[−2,2]

for some constant C.
Proof. By (3.58), for t ∈ [−1, 1]
max |fB (t)| ≤

t∈[−1,1]

(3.62)

2M
X

µn |fB,n | |εϕn (t)|

n=0

≤

2M
X
n=0

√
µn |fB,n | 2 Bµn−3/2 max |B (2πt)|
t∈[−2,2]

where we used (3.47) to write the second inequality. By (3.59) and (3.47), because
αk ∈ R+ , we have
#
"
PM
B −1 k=−M αk
−1
√ −1/2
|fB,n | ≤ µn max |fB (t)|
t∈[−1,1]
+4 Bµn
maxt∈[−2,2] |B (2πt)|
 −1

B (2B − B (0))
√
(3.63)
= µ−1
max
|f
(t)|
B
−1/2
n
t∈[−1,1]
+4 Bµn
maxt∈[−2,2] |B (2πt)|

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

13

PM
where we used the fact that k=−M αk = 2B − B (0) by (3.7), ϕn (ωm ) are eigenP
2
vector obtained by solving (3.48), hence have unit norm, i.e.
m |ϕn (ωm )| = 1,
and maxm |ϕn (ωm )| ≤ 1. Then
max |fB (t)| ≤ C max |B (2πt)| max |fB (t)|

(3.64)

t∈[−1,1]

t∈[−2,2]

t∈[−1,1]

where the constant C is given by

2M 
−3/2
X
2B −1/2 µn
(2B − B (0))
C=
(3.65)
.
+8Bµ−2
n maxt∈[−2,2] |B (2πt)|
n=0

Example 9. Considering the quadratures (αm , ωm ) =



2m
2B
2M +1 , 2M +1

M

of the

m=−M
3
, B (0)

discrete inverse Fourier transform for fixed B, µn = 1 and, by (B.26)
= 0.
Then, using (3.56) and Corollary 8, the error between the nodes is bounded by


maxt∈[−1,1] |fB (t)| maxt∈[−2,2] |B (2πt)|

P2M  1/2
|fB (t)| ≤
× n=0 4B
+ 8B maxt∈[−2,2] |B (2πt)|


maxt∈[−1,1] |f
maxt∈[−2,2] |B (2πt)|
 B (t)|

(3.66)
=
× (2M + 1) 4B 1/2 + 8B maxt∈[−2,2] |B (2πt)|


−2
for t ∈ [−1, 1]. By (B.27), one can achieve maxt∈[−2,2] |B (2πt)| = O (2M + 1)


−1
consequently |fB (t)| ≤ maxt∈[−1,1] |fB (t)| O (2M + 1)
which is in the order
of the truncation errors presented in Section VI of [13].
4. R-limited functions
In this section we introduce an equivalent of R-limited functions with respect to a
linear transformation and R-Slepian functions which are multivariate generalization
of band-limited functions and prolate spheroidal wave functions, respectively. Then
we prove the generalizations of Theorems 1 and 7 to R-limited functions.
Let GL (R, N ) denote the general linear group, the set of invertible matrices in
RN ×N , and

RA = k = Ax | x ∈ R ⊂ RN , A ∈ GL (R, N )
for some compact R ⊂ RN . Employing the terminology introduced in [22], we
define RB -limited functions by
Z
(4.1)
fB (x) =
fˆB (k) ei2πk·x dk
RB

where B ∈ GL (R, N ) is a real symmetric matrix and
Z
ˆ
(4.2)
fB (k) =
fB (x) e−i2πk·x dx.
RN

Here B is a multidimensional analogue of band-limit. When B is the identity
matrix, RI = R, one obtains definition of R-limited functions of [22].
3lim

x→0

sinc (Bx) = limx→0

sin(Bx)(2M +1)−1
sin(Bx(2M +1)−1 )

=1

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

14

Alternatively, we can write
Z Z
fB (x) =
RB

Z
=
ZR

N

R

fB (y) |det (B)| K (2πB (x − y)) dy

=

(4.3)

fB (y) ei2πk·(x−y) dydk


Z
i2πBk·(x−y)
fB (y) |det (B)|
e
dk dy
RN

RN

where det (B) denotes the determinant of B and
Z
(4.4)
K (x) =
eik·x dk.
R

Given a function f (x), its RB -limited projection PB [f ] (x), denoted by fB (x) for
short, is defined by
Z
(4.5)
PB [f ] (x) = fB (x) =
fˆ (k) ei2πk·x dk
RB

or, equivalently, in the convolution representation using convolution theorem
Z
(4.6)
PB [f ] (x) =
f (y) |det (B)| K (2πB (x − y)) dy.
RN

Note that fˆB (k) = fˆ (k), for k ∈ RB .
4.1. Discrete Fourier representation of R-limited approximation of compactly supported functions.
Theorem 10. Consider a discretization of the integral representation (4.4) of
K (x)
X
(4.7)
|det (B)| K (Bx) =
αm eikm ·Bx + K (x)
m

for km ∈ R. For a function f (x) compactly supported within a region X ⊂ Rn , its
R-limited projection can be approximated by
X
fB (x) =
(4.8)
αm ei2πBkm ·x fˆ (Bkm ) + f (x)
m

and (4.8) provides a discretization of (4.1) with
(4.9)

max |f (x)| ≤ |X| max |f (x)| max |K (2πx)| ,
x∈X

x∈X

x∈X+X

where X + X = {x|x = x1 + x2 , x1 , x2 ∈ X} and |X| =

R
X

dx.

Proof. Substituting (4.7) into (4.6), we obtain (4.8) where
Z
(4.10)
f (x) =
f (y) K (2π [x − y]) dy.
X


For fB (x) to accurately approximate the compactly supported function f (x)
over x ∈ X, by (4.9), one needs to build up an approximation of K (Bx) that is
accurate over the set 2π (X + X) = {y|y = 2πx, x ∈ X + X}.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

15

4.2. R-Slepian functions: A multivariate generalization of PSWF. We can
generalize the PSWF and their approximations presented in Sections 3.2.2 and 3.2.3
to multiple-variables to define R-Slepian functions and construct their approximations.
4.2.1. R-Slepian functions. Consider the RB -limited projection of a compactly supported function with support S ∈ RN :
Z
PB [f ] (x) = fB (x) =
(4.11)
f (y) |det (B)| K (2πB (x − y)) dy
S

restricted to x ∈ S. PB [f ] (x) is a positive definite operator. Furthermore, if R is
symmetric, i.e. R = −R = {−x | x ∈ R}, then K (x) = K (−x) is real,
Z
(4.12)
cos (k · x) dk,
K (x) =
R

PB is a positive definite real symmetric operator,
Z Z
f (y) |det (B)| K (2πB (x − y)) dy f (x) dx
S S
Z Z
=
f (y) |det (B)| K (2πB (y − x)) dy f (x) dx
S S
Z
2
(4.13)
=
fˆ (k) dk ≥ 0,
RB

and, consequently, by Mercer’s theorem (see page 245 [21]), accepts a discrete
eigendecomposition
Z
(4.14)
µn ϕn (x) =
ϕn (y) |det (B)| K (2πB (x − y)) dy
S

with positive eigenvalues µn and real eigenfunctions ϕn (x) ∈ R for x ∈ S. We refer
to eigenfunctions ϕn (x) as the R-Slepian functions. For the sake of simplicity of
the discussion we will consider symmetric R. The case of non-symmetric R can be
reduced to the symmetric case by translation of R away from the origin to exclude
origin and consider R ∪ −R.
Consider, S = R and solutions ψ (x) of the equation
Z
(4.15)
λψ (x) =
ψ (k) ei2πBk·x dk, x ∈ R
R

Define ψe (x) = [ψ (x) + ψ (−x)] /2 and ψo (x) = [ψ (x) − ψ (−x)] /2 as the even
and odd parts of ψ (x). For symmetric R, i.e. R = −R, we have
Z
λ [ψe (x) + ψo (x)] =
ψ (k) ei2πBk·x dk
R
 R

ψe (k) cos (2πBk · x) dk
R
R
(4.16)
=
+i R ψo (k) sin (2πBk · x) dk
Considering the equations
Z
(4.17)

ψe (k) cos (2πBk · x) dk

βe ψe (x) =
ZR

(4.18)

ψo (k) sin (2πBk · x) dk,

βo ψo (x) =
R

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

16

which have real symmetric kernel with real eigenvalues and eigenfunctions, the real
and imaginary eigenvalues of (4.15) are associated with eigenfunctions of (4.17)
and (4.18), respectively. Completeness follow from Fourier theory using the same
arguments as in [22, 21].
Eigenfunctions ϕn (x) of (4.15),
Z
λn ϕn (x) =
(4.19)
ϕn (k) ei2πBk·x dk, x ∈ R,
R

are also eigenfunctions of (4.14):
Z
ϕn (y) |det (B)| K (2πB (x − y)) dy
R
Z

Z
i2πBx·k
−i2πBy·k
= |det (B)|
e
ϕn (y) e
dy dk
R
R
Z
ei2πBx·k ϕn (k) dk
= λn |det (B)|
R

(4.20)

2

= |λn | |det (B)| ϕn (x)

with
2

µn = |det (B)| |λn | .

(4.21)

Similar to the case of PSWF, R-Slepian functions satisfy double orthogonality relation [22]
Z
Z
(4.22)
µn
ϕn (x) ϕm (x) dx =
ϕn (x) ϕm (x) dx = δm,n
RN

R

4.2.2. Approximating R-Slepian functions as eigenfunction of ei2πBkm ·xl . Similar
to the single dimensional case, if one can build up an approximation to (4.12) or,
equivalently, (4.4),
(4.23)

K (x) =

M
X

αm exp (ikm · x) + K (x)

m=1

for some km ∈ R, then, substituting (4.23) in (4.14), we obtain
(4.24)

M
X

λn ϕn (x) =

αm ei2πBkm ·x ϕn (km ) + ϕn (x)

m=1

where
Z
(4.25)

ϕn (y) K (2πB [x − y]) .

ϕn (x) =
R

We can approximate R-Slepian functions by substituting the eigenvectors ϕn (xl )
of the equation
(4.26)

λn ϕn (xl ) =

M
X

αm ei2πBkm ·xl ϕn (km ) ,

m=1

xl ∈

[

{km }

m

into (4.24) for n, l = 1, . . . , M .
Because there are |det (B)| |S| |R| number of eigenvalues µn close to one (see
Theorem 3 in [16]4), in order to capture all the eigenvalues close to one, one shall
4Recently this theorem is rediscovered in [8].

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

17

have M ≥ d|det (B)| |S| |R|e. A detailed analysis of the characterization of the
eigenvalues of the projection operator defined in (4.11) around one, zero and the
transition zone can be found in [25].
4.2.3. Approximating R-Slepian functions as eigenfunction of K (2πB [km − kl ]) .
Multiplying both sides of (4.24) by |det (B)| e−i2πBk·x , integrating over x and using
(4.21), we obtain
2

|λn | |det (B)| ϕn (k) = µn ϕn (k)
=

(4.27)

M
X

αm |det (B)| K (2πB [k − km ]) ϕn (km ) + εϕn (k)

m=1

where
Z
εϕn (k) =

(4.28)

ϕn (x) e−i2πBk·x dx

R

Consequently, by (4.27), we can approximate R-Slepian functions by substituting
the eigenvectors ϕn (km ) of the equation
(4.29)

µn ϕn (kl ) =

M
X

αm |det (B)| K (2πB [kl − km ]) ϕn (km )

m=1

for n = 1, . . . M .
4.3. Discrete convolution representation of R-limited approximation of
compactly supported functions. Using R-Slepian functions, we can prove a
generalization of the sampling and interpolation theorem, Theorem 7, for RA limited functions. To do this, we first show the sampling theorem of RB -limited
functions for symmetric B and then generalize it to RA -functions for an arbitrary
A ∈ GL (N, R).
Lemma 11. Given a symmetric B ∈ GL (N, R), i.e. B = B T and a function
f (x), RB -limited projection fB (x) of f (x) can be approximated by
(4.30)

fB (x) ≈

M
X

fm |det (B)| K (2πB (x − km ))

m=1

where
(4.31)

fm =

M
X

fB (kn ) αn |det (B)| K (2πB (km − kn )) αm .

n=1

Proof. Consider a symmetric B ∈ GL (N, R). Because ϕn (x) are complete for
compactly supported functions as well as RB -limited functions, we can expand any
RB -limited function fB using (4.27) as follows:
Z
fB (x) =
f (x) |det (B)| K (2πB (x − y))
R

=

M
X

µn fB,n ϕn (x)

n=1

(4.32)

(P
)
M
αm |det (B)| K (2πB (x − km )) f˜B (km )
m=1
PM
=
+ n=1 µn fB,n εϕn (x)

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

18

where
f˜B (k) =

(4.33)

M
X

fB,n ϕn (k)

n=1

(4.33) is a multivariate generalization of (3.33)
We rewrite fB,n as
(4.34)

fB,n =

Z
M
1 X
αm ϕn (km ) fB (km ) +
fB (x) εϕn (x) dx.
µn m=1
RN

Substituting (4.34) in (4.33), we have
(4.35)

f˜B (k) =

M
X

αm fB (km ) Rm (k) +

m=1

M
X

!

Z
ϕn (k)

n=1

fB (x) εϕn (x) dx
RN

where
M
X
1
ϕn (km ) ϕn (k) .
Rm (k) =
µ
n=1 n

(4.36)

Similar to the single dimensional case, a regularized approximation to Rm (k) is
given by |det (B)| K (2πB (k − km )), leading to the approximate interpolation formula
(
)
M
X
αmh |det (B)| K (2πB (x − km ))
i
PM
(4.37)
fB (x) ≈
×
n=1 αn fB (kn ) |det (B)| K (2πB (km − kn ))
m=1


Theorem 12. Consider the RA -limited function
Z

(4.38)
fA (x) =
f (y) |det (A)| K 2πAT (x − y) dy
AR

for some, not necessarily symmetric, A ∈ GL (N, R). Then
X

(4.39) fA (x) =
αm |det (B)| K 2πAT (x − Akm ) f˜A (Akm )
m

+

X


µn gB,n εϕn A−1 x ,

x ∈ RN , Akm ∈ RA

n

where f˜A (Akm ) and gB,n are defined by (4.44) and (4.43), respectively.
Proof. Let B = AT A. Then gB (x) = fA (Ax) is an RB -limited projection of
g (x) = f (Ax):
Z

2
fA (Ax) =
f (Ay) |det (A)| K 2πAT (Ax − Ay) dy
ZR
=
f (Ay) |det (B)| K (2πB (x − y)) dy
ZR
(4.40)
gB (x) =
g (y) |det (B)| K (2πB (x − y)) dy
R

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

19

which, by (4.32), can be approximated by
(4.41) gB (x) =

X

αm |det (B)| K (2πB (x − km )) g̃B (km )

m

+

X

µn gB,n εϕn (x) ,

x ∈ RN , km ∈ R.

n

or equivalently
(4.42) fA (x) =

X


αm |det (B)| K 2πAT (x − Akm ) f˜A (Akm )

m

+

X


µn gB,n εϕn A−1 x ,

x ∈ RN , Akm ∈ RA .

n

Here
Z
gB,n =

fA (Ax) ϕn (x) dx
Z

1
fA (x) ϕn A−1 x dx
=
det (A) RN
Z
=
fA (Ax) ϕn (x) dx
R
 PM

1
Rm=1 αm fA (Akm ) ϕn (km ) ,
=
+ R fA (Ax) εϕn (x) dx
µn
RN

(4.43)
and
(4.44)

f˜A (Ak) = g̃B (k) =

X

gB,n ϕn (k)

n

with ϕn (x) being the eigenvector of the projection operator PB with respect to the
symmetric matrix B.

4.4. Construction of discrete Fourier approximation of the kernel (4.7). It
is important to note that both Theorems 10 and 12 rely on finding an approximation
of the convolution kernel in the form of (4.23). Although there is no unique way
of finding an approximation in the form of (4.23), it can be constructed using tools
from approximation theory. Considering that two and three dimensional domains
can be approximated using tetrahedral and triangular meshes along with their
multidimensional extensions [4, 18], it is necessary to build quadrature (αm , km ) to
approximate triangle-limited and tetrahedral-limited (shortly T -limited) functions.
Because our results do not require RA to be
they can be generalized
 connected, and
N
to RA = ∪L
and Rl ⊂ RN such that intersection of
l=1 Al Rl where Al ∈ GL R
L
{Al Rl }l=1 has measure zero, i.e. ∩L
l=1 Al Rl = 0, quadratures obtained for T limited functions can be patched together to construct an approximation of the
form (4.23). In this case, the convolution kernel becomes
(4.45)

KΣ (x) =

M
X
m=1

|det (Am )| Km 2πATm x



SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

20

R
where Km (x) = Rm eik·x dk and, consequently, the RA -limited projection of a
function is given by
Z
f (y) KΣ (x − y) dy =

fA (x) =

(4.46)

RN

L
X

fAl (x)

l=1

where
Z
(4.47)

fAl (x) =
RN


f (y) |det (Al )| K 2πATl (x − y) dy.

By Corollary 12, because each fAl (x) can be approximated by
(4.48) fAl (x) ≈

X


αml |det (Bl )| K 2πATl (x − Akml ,l ) g̃Bl (kml ,l )

ml

+

X


µml gB,ml εϕml A−1
l x ,

x ∈ RN , Al kml ,l ∈ Al Rl

ml

then fA (x) can be approximated using samples of f (x) for x ∈ {Al kml ,l }

l=1,...,L .
ml =1,...,Ml

In Appendix C, we provide a method to construct quadratures for T -limited
functions, specifically for isosceles triangle and trirectangular tetrahedron which
are used to construct quadratures for equilateral triangle and regular tetrahedron,
respectively, that satisfy the corresponding symmetry properties. We present two
ways to construct the quadrature for equilateral triangle, one capturing the symmetries of the triangle and the other doesn’t. Although approximations are constructed
to capture the behavior of the kernel and its derivatives at zero, the quadrature
that satisfy the symmetry properties provide a more accurate approximation within
a larger vicinity of zero with fewer number of nodes.
Another special case of R-limited functions that have practical importance in
multidimensional signal processing seismic data is considered in Appendix D which
can also be extended to image processing in 2D and video processing 3D. We present
the convolution kernels and construction of corresponding quadratures that can be
used in sampling and interpolation Theorems 10 and 12 in Appendix D.

5. Conclusion
In this manuscript, we proved duality between the discretization of Fourier and
convolution representations of R-limited functions which lead to the sampling and
interpolation theorem, Theorem 12, where the interpolation is to be understood
as an approximation within a desired accuracy over a compact region. Because
discretization of the Fourier representation is over a compact support, so is the discretization of the convolution representation. Thus, an R-limited function can be
approximated from samples over a compact support that is similar to R. We provided examples of convolution kernels for some special cases of R-limited functions,
namely T -limited and C-limited functions whose Fourier transforms are supported
in a triangle, or tetrahedron, and cone, respectively. We constructed discretization
of the Fourier representation of these kernels which can be used along with sampling
and inteprolation theorems.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

21

Acknowledgments
I would like to thank Lucas Monzon who introduced moment problems and
quadrature methods to me. He has been a mentor, a colleague and most importantly
a dear friend. Initial sketches of results in Sections 3.2.3 and 4.2.3 were obtained
with him in 2012. He also reviewed the manuscript in detail which made the content
more accurate and clearer. Our collaboration wouldn’t have been possible without
support and trust of Konstantin Osypov. I would like to thank Kemal Özdemir
for extensive discussions during the final preparation of this manuscript and his
invaluable signal processing perspective. I would like to thank Ozan Öktem and
Daan Huybrechs for giving me the opportunity to present parts of this work at KTH
Royal Institute of Technology and University of Leuven where follow up questions,
comments and discussions have improved the flow and content of the manuscript.
Finally I thank Garret Flagg and Vladimir Druskin for reading the initial draft of
the manuscript and providing constructive feedback.
References
[1] G. Beylkin and L. Monzón. On generalized Gaussian quadratures for exponentials and their
applications. Applied and Computational Harmonic Analysis, 12(3):332–373, 2002.
[2] William L Briggs et al. The DFT: An Owners’ Manual for the Discrete Fourier Transform.
SIAM, 1995.
[3] Emmanuel J Candes. Multiscale chirplets and near-optimal recovery of chirps. Technical
report, Technical Report, Stanford University, 2002.
[4] Siu-Wing Cheng, Tamal K Dey, and Jonathan Shewchuk. Delaunay mesh generation. CRC
Press, 2012.
[5] AJW Duijndam and MA Schonewille. Nonuniform fast fourier transform. Geophysics,
64(2):539–551, 1999.
[6] Alok Dutt and Vladimir Rokhlin. Fast fourier transforms for nonequispaced data. SIAM
Journal on Scientific computing, 14(6):1368–1393, 1993.
[7] Alan Edelman, Peter McCorquodale, and Sivan Toledo. The future fast fourier transform?
SIAM Journal on Scientific Computing, 20(3):1094–1114, 1998.
[8] Massimo Franceschetti. On Landau’s eigenvalue theorem and information cut-sets. Information Theory, IEEE Transactions on, 61(9):5042–5051, 2015.
[9] F.G. Friedlander and M.S. Joshi. Introduction to the Theory of Distributions. Cambridge
University Press, 1998.
[10] Walter Gautschi. Moments in quadrature problems. Computers & Mathematics with Applications, 33(1):105–118, 1997.
[11] I. S. Gradshteyn and I. M. Ryzhik. Table of integrals, series, and products. Elsevier/Academic
Press, Amsterdam, seventh edition, 2007. Translated from the Russian, Translation edited
and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows,
Macintosh and UNIX).
[12] Leslie Greengard and June-Yub Lee. Accelerating the nonuniform fast fourier transform.
SIAM review, 46(3):443–454, 2004.
[13] Abdul J Jerri. The Shannon sampling theorem – Its various extensions and applications: A
tutorial review. Proceedings of the IEEE, 65(11):1565–1596, 1977.
[14] Abdul J Jerri. The Gibbs phenomenon in Fourier analysis, splines and wavelet approximations, volume 446. Springer Science & Business Media, 2013.
[15] Sun-Yuan Kung. A new identification and model reduction algorithm via singular value decomposition. In Proc. 12th Asilomar Conf. Circuits, Syst. Comput., Pacific Grove, CA, pages
705–714, 1978.
[16] HJ Landau. On Szegö’s eingenvalue distribution theorem and non-Hermitian kernels. Journal
d’Analyse Mathématique, 28(1):335–357, 1975.
[17] HJ Landau and Harold Widom. Eigenvalue distribution of time and frequency limiting. Journal of Mathematical Analysis and Applications, 77(2):469–481, 1980.
[18] Daniel SH Lo. Finite Element Mesh Generation. CRC Press, 2014.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

22

[19] Steve Mann and Simon Haykin. The chirplet transform: Physical considerations. Signal Processing, IEEE Transactions on, 43(11):2745–2761, 1995.
[20] Andrei Osipov, Vladimir Rokhlin, and Hong Xiao. Prolate spheroidal wave functions of order
zero. Springer Ser. Appl. Math. Sci, 187, 2013.
[21] Frigyes Riesz and Béla Sz Nagy. Functional analysis. Dover Publications, 1990.
[22] David Slepian. Prolate spheroidal wave functions, fourier analysis and uncertainty - IV: extensions to many dimensions; generalized prolate spheroidal functions. Bell System Technical
Journal, 43(6):3009–3057, 1964.
[23] David Slepian. Prolate spheroidal wave functions, fourier analysis, and uncertainty - V: The
discrete case. Bell System Technical Journal, 57(5):1371–1430, 1978.
[24] David Slepian and Henry O Pollak. Prolate spheroidal wave functions, fourier analysis and
uncertainty - I. Bell System Technical Journal, 40(1):43–63, 1961.
[25] A Sobolev. Pseudo-differential operators with discontinuous symbols: Widom’s conjecture,
volume 222. American Mathematical Society, 2013.
[26] Radomir S Stankovic, Jaakko T Astola, and Mark G Karpovsky. Some historical remarks on
sampling theorem. In Proceedings of the 2006 International TICSP Workshop on Spectral
Methods and Multirate Signal Processing, SMMSP2006, Florence, Italy, pages 2–3, 2006.
[27] Michael Unser. Sampling-50 years after Shannon. Proceedings of the IEEE, 88(4):569–587,
2000.
[28] Wen Yuan Xu and Christodoulos Chamzas. On the periodic discrete prolate spheroidal sequences. SIAM Journal on Applied Mathematics, 44(6):1210–1217, 1984.
[29] Can Evren Yarman and Garret Flagg. Generalization of Padé approximation from rational
functions to arbitrary analytic functions - Theory. Math. Comp., 84:1835–1860, 2015.

Appendix A. Generalization of Padé approximation
Let f (x) and g (x) be two analytic functions related to each other by the Cauchy
integral
Z
(A.1)
f (x) =
ρ (z) g (zx) dz
Γ

for some closed contour Γ ∈ C and a weighting function ρ (z). A generalization of
Padé approximation is achieved by finding a rational approximation to the weighting function
1 X αm
ρ (z) =
(A.2)
+ ρ (z) , z ∈ C
2πi m z − γm
for some distinct γm ∈ C and error ρ (z). Then we refer to
X
(A.3)
f (x) =
αm g (γm x) +  (x)
m

as the generalization of Padé approximation from rational function to analytic functions, for some error function  (x). Substituting the power series expansion of f
and g at zero into (A.1)
" ∞
#
Z
∞
X
X
n
n
(A.4)
f (x) =
fn x = ρ (z)
gn (zx) dz
m=0

m=0

and equating the terms of the series, one obtains that the moments of ρ (z) are
given by the ratio of the power series coefficients, which we denote by hn
Z
X
fn
n
(A.5)
=
αm γm
+ n
ρ (z) z n dz = hn =
gn
m
for some error n . Because (A.3) is a discrete approximation to the integral (A.1)
, (αm , γm ) are referred to as the quadratures. Individually, we refer to αm and γm

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

g (x)

f2n

g2n

Moment problem
P
2n
hn =
m α m γm

sinc (Bx)

cos (x)

J0 (Bx)

cos (x)

(−1)n B 2n
(2n+1)!
(−1)n B 2n

(−1)n
(2n)!
(−1)n

B 2n
2n+1
(2n)!B 2n

(2n n!)2
(−1)n B 2n
n!
(−1)n B 2n
(2n+1)!
(−1)n B 2n

(2n)!
(−1)n
(2n)!
(−1)n
n!
(−1)n

(2n n!)2
(2n)!B 2n
n!
(−1)n B 2n
(2n+1)!
(2n+1)!B 2n

(2n n!)2
(−1)n
22n+1 n!(n+1)!

(2n+1)!
(−1)n

(2n n!)2
(2n+2)!
22n+1 n! (n+1)!

f (Bx)

e−(Bx)

2

sinc (Bx)

cos (x)
exp −x2



J0 (Bx)

sinc (x)

(Bx)−1 J1 (Bx)

x−1 cosinc (x)

(2n+2)!

23

Quadrature
names
Gauss-Legendre
Clenshaw-Curtis
Gauss-Hermite

Table 1. Example of moment problems corresponding
to approxP
imation of some even functions f (Bx) ≈ m αm g (γm x) in terms
of other even functions g (x), along with the known quadrature
names (see [10]).

as weights and nodes, respectively. In [29], we presented the detailed theory of this
generalization of Padé approximation and a method to compute the quadratures
(αm , γm ) which is based on [15]. Some examples of moment problems are given in
Table 1.
Appendix B. On approximations of sinc
A good approximation to sinc (x) within the vicinity of zero can be achieved by
building up quadratures for the integral representation of sinc (Bx),

(B.1)

sinc (Bx) =

1
B

Z

B

cos (ωx) dω,
0

and then rescaling the approximation by 1/B. One way to do this is using the
method presented in Appendix A to obtain
X
sinc (Bx) =
(B.2)
αm cos (ωm x) + B (x) ,
m

where (αm , ωm ) satisfies the moment problem
X
fn
2n
(B.3)
=
αm ωm
+ n
hn =
gn
m
for some small |n | (see Figure B.2). Here
n

(B.4)

fn =

(−1) B 2n
(2n + 1)!

and
n

(B.5)

gn =

(−1)
(2n)!

are the Taylor series coefficients of sinc (x) and cos (x) at zero, respectively, and
(B.6)

B (x) =

∞
X
n=0

n x2n .

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

24

Solution to the moment problem is equivalent to computing Gauss-Legendre quadratures.
The approximation given in equation (B.2) yields a highly accurate approximation to the sinc (Bx) in a neighborhood of zero (see Figure B.1). However, due to
the rapid increase in the values of the moments hn for large band-limit B, the construction of this approximation suffers from numerical instabilities, and therefore
requires B to be in the range 0 < B ≤ 2. To overcome these challenges, approximation (B.2) can be coupled with a scaling property of the sinc, for example



1
2 cos 2 3n−1 Bx + 1 sinc 3n−1 Bx ,
3
to derive an error bound on the approximation of sinc (3n Bx) in terms of the error
in the approximation to the lower bandwidth sinc (Bx) as a sum of scaled cosines:
(B.7)

sinc (3n Bx) =

Lemma 13. Let
B (x) = sinc (Bx) −

(B.8)

X

αm cos (Bθm x) .

m

Then
n

(3 −1)/2
X
1 X
(B.9) sinc (3 Bx) − n
αm
cos ((θm + 2k) Bx) ≤ |B (x)| ,
3 m
n
n

for n ≥ 0.

k=−(3 −1)/2

Proof. We will prove by induction. For n = 0, this is trivial by assumption (B.8).
Let us define
n

(B.10) 3n B (x)

n

sinc (3 Bx) −

=

(3 −1)/2
X
1 X
α
cos ((θm + 2k) Bx) .
m
n
3 m
n
k=−(3 −1)/2

For n = 1, substituting (B.2) into the scaling property (B.7), we obtain
(B.11)
sinc (3Bx) =

1
X
1X
1
αm
cos (B (θm − 2k) x) + [2 cos (2Bx) + 1] B (x) ,
3 m
3
k=−1

implying
(B.12)

|3B (x)| =

Multiplying 3n B (x) by
(B.13)

1
3

1
[2 cos (2Bx) + 1] B (x) ≤ |B (x)| .
3

[2 cos (2 3n Bx) + 1], we have

1
[2 cos (2 3n Bx) + 1] 3n B (x) = 3n+1 B (x)
3

which implies
(B.14)

|3n+1 B (x)| ≤ |3n B (x)| ≤ . . . ≤ |3B (x)| ≤ |B (x)| .


Corollary 14. The error is given by



∞
X
πl
n+1
n+1
(B.15)
3 B (x) =
sinc 3
B x−
B (x)
B
l=−∞

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

25

which for x = mπ/B becomes
3n+1 B

(B.16)

 mπ 
B

= B

 mπ 
B

.

Furthermore,
(B.17)

lim 3

n+1

n→∞

∞
X

2B 3n+1 B (x) =



l=−∞

πl
δ x−
B


B (x) .

Proof. At the end of the proof of Lemma 13 we showed that the error satisfies the
scaling property (B.13). Using this, we can write
1
[2 cos (2 3n Bx) + 1] 3n B (x)
3


1 
= 2 2 cos 2 3n−1 Bx + 1 [2 cos (2 3n Bx) + 1] 3n−1 B (x)
3
n


1 Y
2 cos 2 3l Bx + 1 B (x)
= n+1
3

3n+1 B (x) =

(B.18)

l=1

which in the Fourier domain can be written as




1
δ k − 2 3l B 
n
ˆ3n+1 B (k) = n+1 ∗l=0
∗ ˆB (k)
+δ k + 2 3l B + δ (k)
3
!
∞
X
χ[−3n+1 B,3n+1 B] (k)
(B.19)
2B
δ (k − 2Bl) ∗ ˆB (k)
=
3n+1 2B
l=−∞

where
−1

Z

ˆB (k) = (2π)

(B.20)

B (x) eikx dx

is the inverse Fourier transform of B (x) and ∗nl=0 fl (k) = (f0 ∗ f1 ∗ · · · ∗ fn ) (k)
denotes a cascaded convolution operator. Taking the inverse Fourier transform, we
obtain

(B.21)

∞
X

!
2πl
3n+1 B (x) = sinc 3
Bx ∗
δ x−
B (x)
2B
l=−∞



∞
X
πl
=
sinc 3n+1 B x −
B (x)
B
n+1





l=−∞

which for x = mπ/B is
3n+1 B
(B.22)

 mπ 
B

=

=

∞
X
l=−∞
∞
X
l=−∞

  mπ 
sinc 3n+1 π [m − l] B
B
δml B

 mπ 
B

.

By using the identity
(B.23)

lim 2a sinc (ax) = δ (x)

a→∞

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

26

Algorithm 1 Representation of sinc (B0 x) as a sum of scaled cosines
Given 0 ≤ B0 ∈ R
(1) Compute n = log3 bB0 c + 1
(2) Set B = B0 3−n .
(3) Solve the moment problem
X

−1
2 n
hn = B 2n (2n + 1) =
αm ωm
m

for (αm , ωm ) using the method of [29] (see Appendix A and Figure B.2)
(4) Set θm = ωm B −1
(5) Form the approximation (see Figure B.1)
1 X
αm
sinc (B0 x) ≈ n
3 m

(3n −1)/2

X

cos ((θm + 2k) Bx)

k=−(3n −1)/2

and the dominated convergence theorem (see page 14 of [9]), we have
n+1

lim 3

n→∞

(B.24)

∞
X

2B 3n+1 B (x) = lim

n→∞

=

3

n+1



2B sinc 3

l=−∞

∞
X
l=−∞

n+1



πl
B x−
B (x)
B



πl
B (x) .
δ x−
B


The practical implications
of Lemma 13 are as long as one has a good approxiP
mation sinc (x) ≈ m αm cos (θm x) over an interval around zero, the (αm , θm ) can
be used to build up an approximation of (i) sinc (Bx), for any B ∈ R, on the same
interval, (ii) sinc (x) on any interval around zero or, equivalently, (iii) sinc (Bx),
for any B ∈ R, on any interval around zero, as accurate as the initial approximation to sinc (x). (see Figures B.1 and B.3) Algorithm 1 outlines our approach to
approximating a sinc of arbitrary bandwidth as a sum of scaled cosines.
How about uniform sampling? Consider approximation of the integral
Z B
1
(B.25)
f (x) =
eixω dω = sinc (Bx)
2B −B
by discretization of the integral using uniform sampling over the interval [−B, B]:
N
X
2B
1
2B
f˜B,N (x) =
eix 2N +1 n
2B 2N + 1
n=−N

(B.26)

1
sin (Bx)


=
2N + 1 sin Bx
2N +1

Because f˜B,N (x) is periodic with period πB −1 (2N + 1), it is also referred to as
periodic sinc function.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

approximation to sinc(0.74074x)
1
approx
org

0
-1
-100

-80

-60

-40

-20

0

20

40

60

80

100

approximation to sinc(20x)
1
approx
org

0
-1
-100
0

-80

-60

-40

-20

0

20

40

60

80

100

logarithmic absolute errors: log10 |ǫ 20 (x)| and log 10 |ǫ 0.74074(x)|
ǫ20 (x)

-10

ǫ0.74074 (x)

-20
-100

-80

-60

-40

-20

0

20

40

60

80

100

Figure B.1. Approximation of sinc (Bx) by (B.2) using Algorithm 1. On top and middle plots, sinc (Bx) and sinc (B0 x)
(red dashed) along with their approximations (solid blue), for
B = 3−3 20 and B0 = 20, respectively. On the bottom plot, the
logarithmic absolute errors for B (red dashed) and B0 (blue solid).
As derived the error corresponding to B0 is less than that of B.
0.15

αm

0.1

0.05

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

ωm = B θ m
6

×10 -3

αm

4
2
0
0

2

4

6

8

10

12

14

16

18

20

B (θ m +2k)

Figure B.2. On top plot, solution (αm , ωm ) to the moment problem (B.3) for B = 3−3 B0 and B0 = 20 which is used to approximate sinc (Bx) (see top plot in Figure B.1). On the bottom plot,
(33 −1)/2
(αm , |θm + 2k| B)k=−(33 −1)/2 used to approximatesinc (B0 x) (see
middle plot in Figure B.1).

27

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

28

approximation to sinc(0.74074x)
1
approx
org

0
-1
-100

-80

-60

-40

-20

0

20

40

60

80

100

approximation to sinc(20x)
1
approx
org

0
-1
-100
0

-80

-60

-40

-20

0

20

40

60

80

100

logarithmic absolute errors: log10 |ǫ 20 (x)| and log 10 |ǫ 0.74074(x)|
ǫ20 (x)

-10
-20
-100

ǫ0.74074 (x)

-80

-60

-40

-20

0

20

40

60

80

100

Figure B.3. Approximation of sinc (Bx) by (B.26). On top
and middle plots, sinc (Bx) and sinc (B0 x) (red dashed) along
with their approximations f˜B,N (x) and f˜B0 ,N0 (x)
 (solid blue) us−3
ing uniform
sampling,
for
(B,
N
)
=
3
20,
13
and (B0 , N0 ) =

20, 33 13 , respectively. On the bottom plot, the logarithmic absolute errors for B (red dashed) and B0 (blue solid). As derived
the error corresponding to B0 is less than that of B.
−1

Using the Taylor series expansion of (1 − x) and sinc (x) around zero 5 series
representation of the error becomes
B (x) = f (x) − f˜B,N (x)


sin (Bx) 
1
h

i 
=
1−
Bx
Bx
1 − 1 − sinc 2N
+1
"∞
2n#m
∞
X
X (−1)n  B
m+1
= sinc (Bx)
(−1)
(B.27)
x
(2n + 1)! 2N + 1
m=1
n=1


−2
within the vicinity of zero and increases away
which decays like O (2N + 1)
−1

from zero for |x| < π (2B) (2N + 1). Consequently, maximum absolute error is
−1
obtained at |x| = π (2B) (2N + 1) which is




 π −1 
(2N + 1) π
(2N + 1) π
= sinc
1 − sinc
B
2B
2
2
π

2
(B.28)
=
−1
(2N + 1) π 2

and decays in the order of N . For (N, B) = 13, 20 × 3−3 and (N0 , B0 ) =

13 × 33 , 20 we present f˜B,N (x) and f˜B0 ,N0 (x) in Figure B.3.
(−1)n
5,(1 − x)−1 = P∞ xm and sinc (x) = P∞
2n
m=0
n=0 (2n+1)! x

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

29

Lemma 13 can be generalized to any arbitrary function:
Theorem 15. Let
B (x) = sinc (Bx) − f (x) .

(B.29)
Then

 n

(3 −1)/2
X
1
(B.30) sinc (3n Bx) − n 
2 cos (2Blx) + 1 f (x) ≤ |B (x)| , for n ≥ 0.
3
l=1

Proof. The proof follows from arguments that are similar to those in Lemma 13.



Corollary 16. Let
B (x) = sinc (Bx) −

(B.31)

X


αm exp −γm x2 ,

m

such that Re {γm } > 0. Then
n

(B.32)
where am,l

(B.33)

(3 −1)/2
X
1 X
sinc (3 Bx) − n
am,l gm,l (x) ≤ |B (x)| , for n ≥ 0.
3 m
l=−(3n −1)/2


(Bl)2
= αm exp i Im{γ
and
m}

gm,l (x) = exp −Re {γm } x2

2 !
Bl
× exp −iIm {γm } x −
.
Im {γm }
n

Proof. This is a direct consequence of Theorem 15 and identities
(B.34)
 n

−1 /2
(3n+1
(3 −1)/2
X)
X


2


2 cos (2Blx) + 1 exp −γm x =
exp −γm x2 + i2Blx
l=−(3n+1 −1)/2

l=1


exp −γm x2 + i2Blx
(B.35)
= exp −Re {γm } x

2





Bl
exp −iIm {γm } x −
Im {γm }

2 !

2

(Bl)
exp i
Im {γm }

!


Corollary 16 says that the sinc function can be approximated as a sum of shifted,
Gaussian tapered chirps. One can determine (αm , γm ) using the method in Appendix A by solving the appropriate moment problem (see Step 3 of Algorithm 2).
This type of approximations of sinc (x) can be used to construct a multiresolution
scheme for band-limited function as an alternative to existing multiscale approaches.
It is important to point out that unlike chirplet decomposition methods presented
in [19, 3], the moment problem provides an explicit solution for (αm , γm ) while
coupling the real and imaginary part of the complex Gaussian parameters γm . Algorithm 2 outlines approximating a sinc of arbitrary bandwidth as a sum of scaled
cosines based on the moment problem and Corollary 16. A corresponding example
is presented in Figure B.4.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

30

Algorithm 2 Representation of sinc (B0 x) as a sum of chirplets
Given 0 ≤ B0 ∈ R
(1) Compute n = log3 bB0 c + 1
(2) Set B = B0 3−n .
(3) Solve the moment problem
X
n!
n
hn = B 2n
=
αm γm
(2n + 1)!
m
for (αm , γm ) using the method of [29] (see Appendix A)
(4) Form the approximation





(Bl)2
(3n −1)/2
exp i Im{γ
exp −Re {γm } x2
X
X
}
m
1



2  
αm 
sinc (B0 x) ≈ n

Bl
3 m
× exp −i Im {γm } x − Im{γ }
l=−(3n −1)/2
m

approximation to sinc(0.74074x)
1
approx
org

0
-1
-100

-80

-60

-40

-20

0

20

40

60

80

100

approximation to sinc(20x)
1
approx
org

0
-1
-100
0

-80

-60

-40

-20

0

20

40

60

80

100

logarithmic absolute errors: log10 |ǫ 20 (x)| and log 10 |ǫ 0.74074(x)|
ǫ20 (x)

-10
-20
-100

ǫ0.74074 (x)

-80

-60

-40

-20

0

20

40

60

80

100

Figure B.4. Approximation of sinc (Bx) as a sum of chirplets
(see Corollary 16) using Algorithm 2. On top and middle plots,
sinc (B x) and sinc (B0 x) (red dashed) along with their approximations (solid blue), for B = 3−3 20 and B0 = 20, respectively. On
the bottom plot, the logarithmic absolute errors for B (red dashed)
and B0 (blue solid). As derived the error corresponding to B0 is
less than that of B.
Appendix C. T -limited functions
In this section we consider triangle and tetrahedral limited functions both of
whom are referred to as T-limited functions. The distinction of two class of Tlimited function should be clear from the dimensions of their variables.
C.1. Triangle-limited functions. We say a function f (x, y), (x, y) ∈ R2 is triangle limited if its Fourier transform fˆ (kx , ky ) is supported within a triangular region

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS
ℜ{K}

31

ℑ{K}

-1

-1

600

3000
2500

-0.5

400

-0.5

2000
1500

0

200
0

0

1000
500

0.5

-200

0.5

-400

0

1

-600

1
-1

-0.5

0

0.5

1

-1

ℜ{F [K]}

-0.5

0

0.5

1

×10 -15

ℑ{F [K]}

-50

-50

6
1
4
0.8
2
0.6

0

0

0

0.4

-2

0.2

-4

0

50
-50

0

50

-6

50
-50

0

50

Figure C.1. Real,
=
imaginary parts of K (x, y)
K4 (x, y) exp −i2π∆p 23 x and its Fourier transform for ∆p = 75
√ −1
and s = 3 . The horizontal axis are x and kx axis and the
vertical axis are y and ky , appropriately.
T ⊂ R2 . Without loss of generality let T ⊂ R2 be parametrize by
(C.1)

T = {(kx , ky ) | 0 ≤ kx ≤ ∆p, |ky | ≤ kx s, }

for some ∆p, s ∈ R+ . Define K4 (x, y) to be
Z ∆p Z kx s
K4 (x, y) =
ei2π(kx x+ky y) dky dkx
0
−kx s


∆p cosinc (2π∆p (x + sy)) − cosinc (2π∆p (x − sy))
(C.2)
=
2πy −i [sinc (2π∆p (x + sy)) − sinc (2π∆p (x − sy))]
We refer to K4 as the kernel for T -limited functions. Let us consider the kernel


2
K (x, y) = K4 (x, y) exp −i2π∆p x
3
√ −1
for s = 3 . Fourier transform of K (x, y) is the characteristic function over an
equilateral triangle whose center of mass is at the origin. For ∆p = 75, we present
the kernel K (x, y) and its Fourier transform in Figure C.1. Next we will give two
ways to construct quadratures for discrete Fourier approximation of K4 (x, y) which
can be generalized to construct quadratures for simplexes in higher dimensions, too.
Proposition 17. K4 (x, y) satisfies the following scaling property

 

1 K4 x2 , y2 1 + 2eiπ∆px
cos (π∆psy)

(C.3)
K4 (x, y) =
y
4 +ei2π∆px K4 − x2 , 2
Consequently, let K4,m (x, y) = K4 (2m x, 2m y). Then
 

m
1 K4,m−1 (x, y) 1 + 2eiπ2 ∆px cos (π2m ∆psy)
(C.4)
K4,m (x, y) =
m+1
∆px
4 +eiπ2
K4,m−1 (−x, y)

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

32

Proof. This is a direct consequence self similarity of isosceles triangle which is used
to decompose integral representation of K4 (x, y) using the identity
Z

∆p

Z

kx s

dky dkx

(C.5)
−kx s

0

" R ∆p/2 R kx s
=

and 14 K4

x y
2, 2



=

R ∆p/2 R kx s
0

−kx s

R ∆p R (∆p−k )s
+ ∆p/2 −(∆p−kx x )s
−kx s
R ∆p R kx s
R ∆p R (k −∆p)s
+ ∆p/2 (∆p−k
+ ∆p/2 −kxx s
x )s
0

#
dky dkx

ei2π(kx x+ky y) dky dkx



Corollary 18. Let K4,m (x, y) = K4 (2m x, 2m y). Then
 
m
K4,m (x, y) 1 + 2e−iπ2 ∆px cos (π2m ∆psy)
m+1

∆px
−eiπ2
K4,m (−x, y)


cos2 (π2m ∆psy) + cos (π2m ∆px) cos (π2m ∆psy) 


(C.6)



K4,m−1 (x, y) = 


Proof. Considering K4,m (x, y) and K4,m (−x, y), computation of
(C.7)

K4,m (x, y)
K4,m (−x, y)
m ∆px
m+1 ∆px −
−iπ2
iπ2
(1 + 2e
cos (π2m ∆psy))
e


K
(x, y)

1  4,m−1
m
2

=
|1+2eiπ2 ∆px cos(π2m ∆psy)| −1
4 × eiπ2m+1 ∆px (1+2e−iπ2m ∆px cos(π2m ∆psy))


leads to (C.6).

C.1.1. Discrete Fourier approximation of K4 (x, y) . For a discrete representation
of the kernel let us consider a bounded region (x, y) ∈ S. T -limited projection
operator P4 and T -limited projection f4 of f restricted to the region S are defined
by

(C.8)

P4 [f ] (x, y) = f4 (x, y)
Z
=
f (x0 , y 0 ) K4 (x − x0 , y − y 0 ) dx0 dy 0
S

Noticing that argument of K4 ranges over
(C.9)

n
o
S + S = (x, y) | (x, y) = (x1 , y1 ) + (x2 , y2 ) , (xn , yn )n=1,2 ∈ S ,

in order to compute f4 (x, y) accurately over the region S, one should have an
accurate representation of K4 inside S + S .
Recalling (C.2),

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

Z

∆p

Z

kx s

ei2π(kx x+ky y) dky dkx

K4 (x, y) =
0

(C.10)

=

∆p
s
2πi

∆p
≈
s
2πi
(C.11)

−kx s
Z 1
−1

= ∆p2 s

1


ei2π∆p(x+sky y)kx dkx dky

0

"

1

∂x

M
X

#
αm e

i2π∆p(x+sky y)kx [m]

−1

m=1

M
X

αmhkx [m] ei2π∆p x kx [m]
i
R1
× −1 ei2π∆p s kx [m] y ky dky

m=1
2

Z

∂x
Z

≈ 2∆p s

33

M
X
m=1

dky
!

αmhkx [m] ei2π∆p x kx [m] e−i2π2∆p s kx [m] iy
PN (m)
i2π∆p s kx [m] y 2ky [m,n]
×
n=1 βm,n e

!

where (αm , kx [m]) and (βm,n , ky [m, n]) are quadratures for approximating sinc as a
sum of cosines (See (B.2) and consider B equal to 2π∆p (X + sY ) and 4π∆p s kx [m] Y ,
respectively, where X = max(x,y)∈S+S |x| and Y = max(x,y)∈S+S |y|.) for a desired
accuracy, which determines the accuracy of approximating K4 inside S + S.
Corollary 19. Let K̃4,m (x, y) be an approximation of K4,m (x, y) and 4,m (x, y) =
K4,m (x, y) − K̃4,m (x, y) be the associated error. Then
1
[3 |4,m−n (x, y)| + |4,m−n (−x, y)|] .
4
Furthermore, if |4,m0 (x, y)| = |4,m0 (−x, y)| for an m0 ∈ Z, then |4,m (x, y)| ≤
|4,m0 (x, y)| for all m ≥ m0 . Consequently,


K4 (x, y) − K̃4,m 2−m x, 2−m y = 4,m 2−m x, 2−m y

(C.13)
≤ 4 2−m x, 2−m y
(C.12)

|4,m (x, y)| ≤

Proof. Proof by induction using the identity


1 4,m−1 (x, y) 1 + 2ei2π∆px cos (2π∆psy)
4,m (x, y) =
(C.14)
4 +ei4π∆px 4,m−1 (−x, y)
obtained from (C.4).



Corollary 19 says that for approximating K4,m (x, y) with a desired error bound
 over a desired region centered around zero, it is sufficient to find an approximation
to K4,m0 (x, y) for any integer (including negative integers) m0 < m whose error is
less than or equal to  within the vicinity of zero. We present K̃4,m (2−m x, 2−m y)
and |4,m (2−m x, 2−m y)| for m0 = 0 and m = 0, 1, . . . , 7 in Figures C.4 and C.5.
C.1.2. Nodes capturing rotational invariance of equilateral triangle. Let us consider
the equilateral triangle TE with each side equal to 1 and center of mass at the origin.
√ −1
√
TE is equivalent to the triangle T of (C.1) with ∆p = 3/2 and s = 3 . While
we can construct discrete Fourier approximation of the kernel for the equilateral
triangle using the method in Appendix C.1.1, constructed nodes (kx [m] , ky [m, n])
do not necessarily satisfy the rotational invariance of equilateral triangle (see Figure
C.2).

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

m

RNI

RI

0

1

2

3

4

5

6

7

Figure C.2. Rotationally not invariant (RNI column) and rotationally invariant (RI column) Fourier quadratures for K4 (x, y).
The horizontal axis is kx and the vertical axis is ky .

34

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

35

Figure C.3. Construction of nodes that satisfy the rotational invariance of equilateral triangle can be achieved by first constructing nodes for the lightest gray triangle followed by rotating and
accumulating the constructed nodes.
In order to preserve the rotational invariance of the equilateral triangle among
the nodes, one can construct nodes for the isosceles triangle TI which is triangle T
√
√ −1
of (C.1) with ∆p = 3/6 and s = 3 , then rotate these nodes by 2π/3 and 4π/3
to construct nodes satisfying the rotational invariance of the equilateral triangle
(see Figure C.2) as illustrated in Figure C.3. Our observation is, for the same or
less number of nodes, the nodes with rotational symmetry provide a more accurate
discretization of the Fourier approximation of the kernel K4 (x, y) compared to
nodes without rotational symmetry (see Figures C.4 and C.5).
C.2. Tetrahedron limited functions. We say a function f (x, y, z), (x, y, z) ∈ R3
is tetrahedron limited if its Fourier transform fˆ (kx , ky , kz ) is supported within a
tetrahedral region T ⊂ R3 :
Z
(C.15)
f (x, y) =
fˆ (kx , ky , kz ) e2πi(kx x+ky y+kz z) dkz dky dkx
T

Without loss of generality, let T be parametrized by
(C.16)

T = {(kx , ky , kz ) | 0 ≤ kz ≤ h, 0 ≤ ky ≤ ∆pkz , |kx | ≤ ky s, }

for some ∆p, s, h ∈ R+ . Define K/ (x, y, z) to be
K/ (x, y, z)
Z h Z kz ∆p Z
=
0

0
2

(C.17)

h ∆p
=−
2πx

ky s

ei2π(kx x+ky y+kz z) dkx dky dkz

−ky s

( expc(i2πh[z+∆p(y+sx)])−expc(i2πhz)
2πh∆p(y+sx)

−

expc(i2πh[z+∆p(y−sx)])−expc(i2πhz)
2πh∆p(y−sx)

)

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS
K̃m 2−m x, 2−m y

m
ℜ{K̂}

-1



ℑ{K̂}

-1

-K − K̂ -

-1

3000

2000

2000

36

2000

0

0

0

0

0
1000

-2000

-2000

1

0

1

1
-1

0

1

-1

ℜ{K̂}

-1

3000

0

-1

1

ℑ{K̂}

-1

2000

0

1

-K − K̂ -

-1

2500
2000
1500
1000
500

2000
1000

0

0

0

0

0
-1000

1

1

-2000

1
-1

0

1

-1

ℜ{K̂}

-1

3000

0

1000

-1

0

0

2

0

1

-1

ℜ{K̂}

-1

3000

2000
1500

0

1000
500

1

1
-1

1

-1000

0

1

0

-K − K̂ -

1000

2000

0

-1

ℑ{K̂}

-1

1

1

0

-1

1

ℑ{K̂}

-1

2000

0

1

-K − K̂ -

-1

2000

2000

0

1000

1500

0

0

0

1

-2000

1

1000
500

0

1

3

-1

0

1

-1

ℜ{K̂}

-1

3000

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1
500

1000

2000

0

1000

0

0

1

-1

ℜ{K̂}

-1

3000

600
200

1

1
-1

0

-500

0

1

4

0

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1
500

600

2000

0

1000

0

0

400

0

200

1

5

-500

0

1

1
-1

0

1

-1

ℜ{K̂}

-1

3000

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1

400

500

300

2000

0

1000

0

-1

0

1

-1

ℜ{K̂}
3000

200
100

1

1
-1

0

-500

0

1

6

0

-1

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

500

2

2000

0

1000

0

0

1

1
0.5

1

1
-1

1.5

0

-500

0

1

7

0

-1

0

1

-1

0

1

Figure C.4. Approximation of the kernel K4 (x, y) by
K̃4,m (2−m x, 2−m y) using the RNI quadratures in Figure C.2. The
horizontal axis is x and the vertical axis is y.

where expc (x) = (exp (x) − 1) x−1 , implying expc (ix) = sinc (x) + icosinc (x), with
q
cosinc (x) = (1 − cos (x)) x−1 . For an equilateral tetrahedron, choose h =

2
3,

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS
K̃4,m 2−m x, 2−m y

m
ℜ{K̂}

-1

3000

ℑ{K̂}

-1

1000

-1000

0

-1

0

0

-1

1000
500

1

1
-1

3000

1500

0

-1000

1

ℜ{K̂}

-K − K̂ -

1000

0

0

1



-1

2000

0

0

-1

1

ℑ{K̂}

-1

1000

0

1

-K − K̂ -

-1

1500

2000

0

1000

1000

0

0

1

-1000

1
-1

0

1

-1

ℜ{K̂}

-1

3000

0
500

0

1

0

-1

1

ℑ{K̂}

-1

1
0

1

-K − K̂ -

-1

800

500
2000

0

1000

600

0

0

0

1

2

0

1

3000

400
200

1
-1

ℜ{K̂}

-1

0

-500

1
-1

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1

800

500
2000

0

1000

600

0

0

3

0

1

-1

ℜ{K̂}

-1

3000

400
200

1

1
-1

0

-500

0

1

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1

500

500

400

2000

0

1000

0

0

1
0

1

-1

ℜ{K̂}

-1

3000

200
100

1

1
-1

300

0

-500

0

4

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1

400

500

300

2000

0

1000

0

0

0

1

-1

ℜ{K̂}

-1

3000

200
100

1

1
-1

0

-500

0

1

5

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

-1

150

500

2000

0

1000

6

0

0

50

1

1
-1

-1

0

1

-1

ℜ{K̂}
3000

100

0

-500

0

1

-1

0

-1

1

ℑ{K̂}

-1

0

1

-K − K̂ -

500

0.8

2000

0

1000

0

0

0

1

0.4
0.2

1

1
-1

0.6

0

-500

0

1

7

37

-1

0

1

-1

0

1

Figure C.5. Approximation of the kernel K4 (x, y) by
K̃4,m (2−m x, 2−m y) using the RI quadratures in Figure C.2. The
horizontal axis is x and the vertical axis is y.

√
√
∆p = 63 h, s = 3 and add the resulting kernel with its
versions around the z-axis (see Figure C.6).

2
3π

and

4
3π

rotated

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

38

Figure C.6. Real, imaginary  partsh of K
z)
=
 (x, y,
√ i
P2
B −3 n=0 K/ (Bxn , Byn , Bzn ) exp i2πB − 21 z + 12 + 33 y
q
√
√
3
2
and its Fourier transform for h =
3 and
3 , ∆p = 6 h, s =

T
T
B = 20. Here [xn , yn , zn ] = R[0,0,1] 32 πn [x, y, z] for n = 0, 1, 2.

Similar to the triangle-limited case, construction of the nodes is equivalent to discretization of the Fourier integral using cascaded Gauss-Legendre quadratures that
can accurately approximate the representation kernel within a region of interest:
K/ (x, y, z)
Z hZ kz ∆pZ
=
0

0

ky s

ei2π(kx x+ky y+kz z) dkx dky dkz

−ky s

Z 1 Z 1Z 1
h∆ps
∂
∂
ei2πh(∆p[s kx x+y]ky +z)kz dkz dky dkx
y z
4π 2
−1 0 0

X αm βm,n γm,n,l k 2 [m] ky [m, n]
3
2
z
= 2h ∆p s
×ei2πh(∆p[s (2kx [m,n,l]−1)x+y]ky [m,n]+z)kz [m]

=−
(C.18)

m,n,l

where (αm , kz [m]), (βm,n , ky [m, n]) and (γm,n,l , kx [m, n, l]) are quadratures for approximating sinc (Bx) as a sum of cosines (see (B.2)) for B equal to 2πh (∆p (Y + sX) + Z)
and 2πh∆p kz [m] (Y + sX) and 4πh∆p s kz [m] ky [m, n] X, respectively, where X =
max(x,y,z)∈S+S |x|, Y = max(x,y,z)∈S+S |y| and Z = max(x,y,z)∈S+S |z|.
Let us consider a unit tetrahedron, i.e. a tetrahedron with all sides equal to one.
In order to construct nodes that satisfy the symmetries of the unit tetrahedron, first
√
√ −1
√
construct nodes for the sub-tetrahedron with s = 3, ∆p = 2, h = 24 (see
Figure C.7) and then use the symmetry
 group of the regular tetrahedron.
 Namely,
apply rotations Rv̂4 32 πn Rv̂1 32 πm , for n, m = 0, 1, 2, and Rv̂4 23 πn Rv̂2 43 π
to the quadrature of sub-tetrahedron. Here v̂ = v/ |v| is the unit vector pointing

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

39

Figure C.7. Construction of the nodes that satisfies the isometries of the tetrahedron can be obtained by construct nodes for
the shaded tetrahedron (top) and apply symmetry group of the
tetrahedon (see Figure C.8) to the constructed nodes (bottom).
along vector v with vi , for i = 1, 2, 3, 4, given by
"
r #
√
1
3 1 3
(C.19)
v1 = − , −
,−
2
6
6 2
"
r #
√
1
3 1 3
v2 =
(C.20)
,−
,−
2
6
6 2
" √
r #
3 1 3
(C.21)
,−
v3 = 0,
3
6 2
"
r
r #
2 1 3
(C.22)
v4 = 0, 0,
−
.
3 6 2
and


(C.23)

r1 (θ)
Ru (θ) =  r123 (θ)
r132−(θ)


r123 (−θ) r123 (θ)
r2 (θ)
r231 (−θ) 
r231 (θ)
r3 (θ)

is the matrix for a rotation by an angle θ around the unit vector u = [u1 , u2 , u3 ]
with
(C.24)
(C.25)

ri (θ) = cos θ + u2i (1 − cos θ)
rijk (θ) = ui uj (1 − cos θ) + uk sin θ.

The quadrature generated for the regular tetrahedron using the discussed steps is
presented in Figure C.8 .
Appendix D. Cone-limited functions
In seismic or electromagnetic signal processing the signal is modeled through the
wave equation. For an acoustic homogeneous medium with wave speed c = p−1 ,

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

40

Figure C.8. A regular tetrahedron (top left) is symmetric under
rotation with respect to one of its faces (top middle). Its rotational
symmetries with respect to 3 out of 4 faces (top right). There are
a total of 4 × 3 × 2 = 24 isometries of the tetrahedron which form
the symmetry group of the tetrahedron. Out of these 12 of them
preserves orientation which are used to generate the corresponding
quadrature (bottom) from the quadrature of the sub-tetrahedron
in Figure C.7.

the wave equation provides a dispersion relationship between the frequency ω and
wave number k, |k| = ωp. For a heterogeneous medium, the dispersion relationship
becomes an inequality |k| ≤ ωpmax where the maximum slowness pmax = c−1
min is the
one over the minimum speed cmin of the heterogeneous medium. Given the maximum frequency ω0 of the recording system, the Fourier transform of the measurement is supported inside the cone C = {(ω, k) ∈ R × Rn |ω ∈ [−ω0 , ω0 ] , |k| ≤ ωpmax }
for t ∈ R and x ∈ Rn , which is referred to as the signal cone. Temporal and spatial Fourier transform of video images also have their Fourier transforms supported
effectively in a similar cone.
We say a function f (t, x) is cone-limited, C-limited for short, if its Fourier
transform fˆ (ω, k) is supported within the cone C. C-limited functions are invariant
under convolution with the kernel K (t, x), whose Fourier transform, K̂ (ω, k), is
equal to one within C:
Z
(D.1)

K (t, x) =
C

ei2π(ωt−k·x) dω dk.

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

41

For n odd, K (t, x) can be represented in terms of elementary functions. For example, for n = 1 and n = 3, we have


ω0
cosinc (2πω0 [t + pmax x])
(D.2)
K (t, x) =
−cosinc (2πω0 [t − pmax x])
πx
and
(D.3)

ω0
K (t, x) = − ∂r
πr

 

1 cosinc (2πω0 [t + pmax r])
,
r −cosinc (2πω0 [t − pmax r])

respectively. On the other hand, for n even, K (t, x) is a multivariate special function. For example, for n = 2 , we have
2π

Z

Z

pmax

Z

ω0

ei2πω(t−pr cos θ) ω 2 p dω dp dθ

K (t, x) =
0

(D.4)

0

2ω02 pmax
=
πr

Z

−ω0
1

J1 (2πωrω0 pmax ) cos(2πωω0 t) ω dω
0

where |x| = r and Jn (t) is the nth order Bessel function of the first kind. Becuase
all the cases for n even, requires evaluation of an integral of the form (D.4), for the
ease of our discussion we will focus on the case n = 2.
PM
Following [29], J1 (x) can be approximate by J1 (x) ≈
m=1 αm cosinc(γm x),
where (αm , γm ) satisfies a moment problem. (see Table 1 in Appendix A) . Consequently, we have
K (t, x) ≈ K̃ (t, x)
(D.5)



M
sinc(2πω
0 t)


ω0 X αm 

sinc (2πω0 [γm pmax r − t])
= 2
− 21
π m=1 γm
+sinc (2πω0 [γm pmax r + t])

The least square error is given by
(D.6)
Z

2

K(t, x) − K̃(t, x) dt dx =

4ω03 p2max 1 PM
P
α w
α
,
( 2 − m=1 αm γm + M
m,m0 =1 m m,m0 m0 )
π2 3

where
Z
(D.7)

wm,m0 =
R+

sin2 (γm r/2) sin2 (γm0 r/2) dr
,
γm r/2
γm0 r/2
r

which can be explicitly computed using integration by parts and the identities
3.827-3.828 on pages 462-463 of [11].
Because

(D.8)

1
K̃p (t, x) = 2
r





1 sinc (2πω0 [pr − t])
sinc (2πω0 t) −
2 +sinc (2πω0 [pr + t])

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

42

has its Fourier transform
Z
(D.9)
K̃p (t, x) ei2π(ωt−k·x) dt dx

−1

= 2π 2 |ωp| χ[−1,1] ωω0


q
2
2
(ωp) − |k|
,
arcsinh 
|k|
0 < |k| ≤ |ωp| ,

supported over

Cp = (ω, k) ∈ R × R2 |ω ∈ [−ω0 , ω0 ] , |k| ≤ ωp ,

K̃ (t, x) is C-limited within C̃ = (ω, k) ∈ R × R2 |ω ∈ [−ω0 , ω0 ] , |k| ≤ ωpmax maxm {γm } .
If maxm {γm } ≈ 1, then the cone-limit C of K (t, x) is approximated by the conelimit C̃ of K̃ (t, x) which is the case in practice. For ω0 = 50 and pmax = 1, we
present K̃ (t, x) and its Fourier transform in Figure D.1.
For n = 2, discretizaton of the integral representation of K (t, x) can be obtained
by
Z 1 Z 1Z 1
ω2 p
2
3
K (t, x) = 2pmax ω0
dω dp dτ
ei2πφ(ω,p,τ ;t,x,y) √
1 − τ2
−1 0
−1
X
(D.11)
= 2p2max ω03
am,n,l ei2πφ[m,n,l](t,x,y)
(D.10)

m,n,l

where
(D.12)
(D.13)
(D.14)

ω[m]2 p[m,n]

am,n,l =αm βm,n γm,n,l √

1−(τ [m,n,l])2

√
φ(ω,p,τ ;t,x,y)=ω0 ω (t−pmax p[xτ +y 1−τ 2 ])

h
i
√
φ[m,n,l](t,x,y)=ω0 ω[m] t−pmax p[m,n] xτ [m,n,l]+y 1−(τ [m,n,l])2

with (αm , ω [m]) and (βm,n , p [m, n]) are quadratures for approximating sinc (Bx) as
a sum of cosines (see (B.2)) for B equal to 4πω0 (T + pmax R) and 2πω0 ω [m] pmax R,
respectively, and (γm,n,l , τ [m, n, l]) is the quadrature for approximating J0 (Bx)p
as a
sum of cosines for B equal to ω0 ω [m] pmax p [m, n] R, where R = max(t,x,y)∈S+S x2 + y 2 ,
and T = max(t,x,y)∈S+S |t|, for some region of interest S ⊂ R × R2 . While
(αm , ω [m]) and (βm,n , p [m, n]) are equivalent to Gauss-Legendre quadrature, and
computation of (γm,n,l , τ [m, n, l]) requires solving the following moment problem
PM
related to the approximation J0 (x) ≈ m=1 αm cos(γm x), which is equivalent to
finding the Clenshaw-Curtis quadrature (see Table 1 in Appendix A).
Note that quadrature for B-limited functions, whose Fourier transforms are supported within a ball B = {k ∈ Rn | |k| ≤ kmax }, can be generated in a similar fashion
(see Figure D.3).
Schlumberger, High Cross, Madingley Road, Cambridge,CB3 0EL, United Kingdom
(cyarman@slb.com)

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

Figure D.1. Real, imaginary parts of K̃ (t, x) and its Fourier
transform for ω0 = 50 and pmax = 1.

Figure D.2. Quadrature (amnl , ω0 (ω [m] , pmax p [m, n, l])) for Climited functions for ω0 = 1 and pmax = 1.

43

SAMPLING FOR APPROXIMATING R-LIMITED FUNCTIONS

Figure D.3. Quadrature (amnl , kmax km,n,l ) for C-limited functions for kmax = 1.

44

