Shannon-Type Sampling Theory

on Unions of Equally Spaced and

Noncommensurate Grids

ARL-TR-2305 February 2001

Terrence J. Moore 

Approved for public release; distribution unlimited.
The findings in this report are not to be construed as an

official Department of the Army position unless so

designated by other authorized documents.

Citation of manufacturer’s or trade names does not

constitute an official endorsement or approval of the use

thereof.

Destroy this report when it is no longer needed. Do not

return it to the originator.
ARL-TR-2305 February 2001

Army Research Laboratory

Adelphi, MD 20783-1197

Shannon-Type Sampling Theory

on Unions of Equally Spaced

and Noncommensurate Grids

Terrence J. Moore 

Computational and Information Sciences Directorate

Approved for public release; distribution unlimited.
Abstract

Sampling theory is that branch of mathematics that seeks to reconstruct

functions from its values at a discrete set of points. The fundamental result

in sampling theory known as “Shannon’s sampling theorem” has many applications to signal processing and communications engineering. I demonstrate Shannon’s result via complex interpolation methods. I then quote a

result that uses these methods to solve interpolation problems on unions of

noncommensurate lattices, which are created via a specific number of theoretic guidelines. These interpolations give Shannon-type reconstructions

on these lattices. I close by doing simulations in MATLAB of the sampling

reconstructions on these noncommensurate grids.

ii
Contents

1 Introduction 1

2 Preliminaries in Harmonic Analysis 2

3 Classical Sampling Theory 5

3.1 Classical Sampling Theorem ................... 5

3.2 A Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.3 An Example ............................ 7

3.4 Deriving the Cardinal Series ................... 9

3.4.1 Delta Method ....................... 9

3.4.2 Complex Method . . . . . . . . . . . . . . . . . . . . . 10

3.5 A Bit of History . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Multirate Sampling Theory 16

4.1 Multichannel Deconvolution . . . . . . . . . . . . . . . . . . 16

4.2 Noncommensurate Sampling Lattices . . . . . . . . . . . . . 18

4.2.1 Two Sampling Lattices . . . . . . . . . . . . . . . . . . 18

4.2.2 An Arbitrary Number of Lattices . . . . . . . . . . . . 19

4.3 Interpolators at Origin . . . . . . . . . . . . . . . . . . . . . . 20

4.3.1 Two Interpolators . . . . . . . . . . . . . . . . . . . . . 20

4.3.2 Arbitrary Number of Interpolators . . . . . . . . . . . 22

4.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

References 28

Appendixes 31

A Proof of the Polynomial Expansion 31

iii
B MATLAB Code 33

B-1 Shannon.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

B-2 MultirateShannon.m . . . . . . . . . . . . . . . . . . . . . . . 36

Distribution 39

Report Documentation Page 41

Figures

1 Original signal f(x) and its Fourier transform f(ω) ........ 7

2 Original signal f and reconstruction g ............... 8

3 Aliasing caused by undersampling ................. 8

4 Original signal f(t) and bisample reconstruction g(t) . . . . . . 25

5 Interpolator and its derivative . . . . . . . . . . . . . . . . . . . . 26

6 Original signal f(t) and trisample reconstruction g(t) . . . . . . 26

iv
1. Introduction

Sampling theory is the distinctive branch of mathematics that attempts

to solve the interpolation problem of a function with fixed growth from

known sampled values by reconstructing the original function. The theory is a subset of the general theory of interpolation, which constructs

functions that satisfy the known values. Sampling theory uses interpolation with other knowledge of a function to find methods that reconstruct

the original function. This involves restrictions on the samples of known

values.

The most significant contribution to this theory is called the “classical sampling theorem.” This theorem is known by a myriad of names—Whittaker,

Kotel’nikov, and Shannon—but the result is essentially the same. The theorem tells us that if a signal f is a function of finite energy on R (f ∈ L2(R))

with Fourier transform f

(ω)=0 for all |ω| ≥ Ω, and if T ≤ 1/2Ω, then for

all t ∈ R,

f(t) = T ∞

n=−∞

f(nT)

sin( π

T (t − nT))

π(t − nT) . (1)

In other words, if f is sampled uniformly at a sufficiently high rate, f can be

reconstructed from these samples by the formula in equation (1). This formula has a rich history extending back as far as the 1600s. However, much

of the crucial research in developing and establishing the theory has occurred in the twentieth century, primarily because of the research in the

various engineering fields. With the motivation of obtaining useful and

practical results, the classical sampling theorem has been generalized and

extended to a number of applications.

1
2. Preliminaries in Harmonic Analysis

Before I discuss the main results of this report, I first introduce a number

of preliminary definitions and theorems. Much of what follows is standard

notation in the literature and can be found in any text covering Fourier

analysis. (For example, see Benedetto [1], Dym and McKean [2], or Korner  ̈

[3]. Also see Apostol [4], Higgins [5], Marks [6], and Zayed [7].) First, I

define in the following what are absolutely and square integrable functions:

Function f. A function f is called absolutely integrable, i.e., f ∈ L1(R), if



R

|f(x)| dx ≡ ||f||1 < +∞ .

If f is in L1, one can say that ||f||1 is the L1 norm of f. Similarly, a function

is called square integrable, i.e., f ∈ L2(R), if



R

|f(x)|

2 dx ≡ ||f||2 < +∞ .

If f is in L2, one can say that ||f||2 is the L2 norm of f.

In this section, all functions are considered absolutely and square integrable

functions on the real line, unless otherwise noted. Likewise, all integrals

are assumed to be over the whole domain (either R or C, depending on the

context) unless noted otherwise. The Fourier series and Fourier transform

play integral parts in the sampling theory and Fourier analysis. Their definitions, from Benedetto [1] and Dym and McKean [2], are defined in the

following:

Fourier Series. Let f be a periodic, integrable function on R, with period 2Ω.

Then the Fourier series of f is represented as

f(x) = ∞

n=−∞

cn exp−iπnx/Ω , (2)

where the Fourier coefficients cn are defined by

cn = 1

2Ω  Ω

−Ω

f(x) expiπnx/Ω dx . (3)

Fourier Transform and Inversion Formula. The Fourier transform of f ∈ L1(R)

is defined as

f

(ω) = 

R

f(t) exp−2πitω dt , (4)

2
and its inversion formula for f

∈ L1(R) is

f(t) = 

R

f

(ω) exp2πiωt dω . (5)

The factor 2π in the exponential simplifies Plancheral’s formula, i.e., ||f||2 =

||f

||2. Often, the Fourier transform of f is also denoted F. The Fourier transform can be extended to square integrable functions via a density argument

on C∞

(R), the space of infinitely differentiable, rapidly decreasing functions on R. The transform is an isometry in L2(R). By itself, the Fourier

transform has useful applications in signal and image processing. If I consider the function f to be a signal, then I refer to the domain of f as the

time domain and the domain of the transform as representing the signal

frequency. The following definitions are useful when referring to the time

or frequency domain of a signal or function:

Support. The support of a function f : R → R, denoted supp(f), is the

closure of the set on which f is nonzero, i.e., supp(f) = {x ∈ R : f(x) 
= 0}.

The function f is said to have compact support if supp(f) is a compact set.

Band-limited. The function f : R→ Ris called Ω band-limited if the support

of its Fourier transform f

is contained in the bounded interval [−Ω, Ω] in R,

i.e., supp(f

) ⊂ [−Ω, Ω].

The Fourier transform has a multitude of mathematically convenient and

useful algebraic properties: linearity, symmetry, conjugation, translation

(time shifting), modulation (frequency shifting), and time dilation [1]. Also,

the transform has a number of important analytic properties: boundedness,

continuity, Riemann-Lebesgue Lemma, time differentiation, and frequency

differentiation [1–3,5]. One of the more relevant properties of the Fourier

transform involves a mathematical operation called a convolution.

Convolution. For f,g ∈ L1(R), the convolution of f and g is defined by

h(x) = f ∗ g(x) = f(y)g(x − y) dy (6)

for x ∈ R.

Convolution is a mathematical operation that consists of the filtering of

one signal (the source) passing it through a filter, which is simply another

signal.

Convolution is commutative, i.e., f ∗ g = g ∗ f, and if f, g, and h are of compact support, convolution is associative. The subject of convolution is important to signal reconstruction. The operation is the filtering of one function by passing it through a filter, which is simply another function. More

on this topic will be discussed in section 4.1. For now, I note the following

two significant properties from the Fourier perspective:

f

∗ g = f

· g, and

f

· g = f∗ g . 

3
These properties prove to be extremely useful in developing the theory on

deconvolution.

Finally, I conclude with the Paley-Wiener theorem:

Paley-Wiener. The Fourier-Laplace transform of an infinitely differentiable

function f with compact support contained in {|t| ≤ A} is an entire function f

(ζ) in C, which satisfies the following property:

For every integer N ≥ 0, a positive constant C = C(N) exists, such that

|f

(ζ)| ≤ C(1 + |ζ|)

−N exp2πA|Imζ| for all ζ ∈ C . (7)

Conversely, every entire function in C satisfying this property is the

Fourier-Laplace transform of a C∞ function, with compact support within

{|t| ≤ A}.

The classical Paley-Wiener theorem says that a square-integrable complexvalued function, defined over the real line, can be extended off the real line

as an entire function of exponential type if and only if its Fourier transform

f(x) is identically zero for |x| > a, i.e., if and only if f is band-limited to

[−a, a]. (For a derivation, see literature by Dym and McKean [2]. )

The Shannon series is a means of extending to the entire complex plane C.

The extension of the Paley-Wiener theorem to generalized functions (to

tempered distributions) is called the Paley-Wiener-Schwartz theorem. Additional background material can be found in literature by Benedetto [1],

Dym and McKean [2], and Korner [3].  ̈

4
3. Classical Sampling Theory

3.1 Classical Sampling Theorem

Sampling theory is an area of high interest and interacts with numerous

other mathematical fields. The idea behind sampling is to convert continuous time signals into appropriate discrete signals that represent the original

signal and vice versa. Reconstruction reverses the process to retrieve the

original signal from the samples. This theory of sampling was created via

an amalgam of mathematical, scientific, and engineering research. Without

dispute, the fundamental result from the conglomeration of these sources

has been what is called the classical sampling theorem: Let f be a function

of finite energy on R (f ∈ L2(R)) with Fourier transform f

(ω)=0 for all

|ω| ≥ Ω, i.e., f(t) is Ω band-limited.

If T ≤ 1/2Ω, then for all t ∈ R,

f(t) = T ∞

n=−∞

f(nT)

sin( π

T (t − nT))

π(t − nT) . (8)

If T ≤ 1/2Ω and f(nT)=0 for all n ∈ Z, then f ≡ 0. The rate 1

2Ω is known

as the Nyquist rate.

This theorem, though it is simple to state, has far-reaching consequences.

The theorem basically states that some signals, or functions, can be perfectly represented by sample values of the function taken at regular, or equispaced, intervals. Knowledge of the frequency bound is required, since

this determines the minimum rate at which the signal must be sampled

to reconstruct the signal completely. Because Nyquist discovered this minimum rate, it is commonly referred to as the Nyquist rate [8]. No information is lost if the signal is sampled at this rate, and also no additional

information is acquired by sampling faster than this rate. That is, the signal

is already perfectly reconstructed. Essentially, the band-limited condition

guarantees that sufficiently close sampling reveals a nicely behaved function between the sample points. Since f is band-limited, one can gather

data sufficiently faster than the signal’s oscillation in time and then recover

f using the “sinc” function sin(2πx)

x as an interpolator. Thus, one can think

of classical sampling as a uniqueness theorem from the theory of analytic

splines.

There is also a duality in the theorem. Not only can one reconstruct the

function in time, assuming the function is band-limited, but also in frequency, assuming that the function is time-limited. In this case, the function

and its transform switch roles.

There is an intriguing history involving the development of this theorem.

A dispute prevailed over who initially discovered this result and who is

5
credited. Consequently, the theorem has several names associated with it:

the cardinal series, the Whittaker sampling theorem, the Kotel’nikov theorem, and the Shannon sampling theorem. Often, the names are combined

in various forms. For simplicity, it is referred to in this report simply as the

classical sampling theorem.

3.2 A Proof

There are numerous proofs of the classical sampling theorem [1,5–7,9,10],

some more rigorous and satisfying than others. A simple proof of the theorem follows:

Proof. Let f be Ω band-limited, i.e., supp(f

) ⊂ [−Ω, Ω]. Since this condition

holds, f

can also be represented as the restriction of a 2Ω periodic function

to the support of f

as

f

(u) = ∞

n=−∞

cn exp−iπnu/Ω ·X[−Ω,Ω](u) , (9)

where the Fourier coefficients are simplified with the use of the Fourier

inversion formula and the nature of X . Note, XA is the characteristic function defined to be unity in the set A and zero elsewhere. Then the Fourier

coefficients are found to be

cn = 1

2Ω  Ω

−Ω

f

(u) expiπnu/Ω du

= 1

2Ω  ∞

−∞

f

(u) exp2iπ(n/2Ω)u du (10)

= 1

2Ωf( n2Ω) .

Substituting equation (10) back into equation (9) and solving for f with

the use of the inverse Fourier transform in equation (5) and the uniform

convergence of Fourier series, the classical sampling theorem is the result

f(x) =  ∞

−∞

f

(u) exp2πiux du

=

 ∞

−∞

1

2Ω

∞

n=−∞

f( n

2Ω) exp−iπnu/Ω ·X[−Ω,Ω] exp2πiux du

= 1

2Ω

∞

n=−∞

f( n

2Ω)

 Ω

−Ω

exp(2πix−πin/Ω)u du (11)

= 1

2Ω

∞

n=−∞

f( n

2Ω)2Ωsin π(2Ωx − n)π(2Ωx − n)

= ∞

n=−∞

f( n

2Ω)

sin π(2Ωx − n)

π(2Ωx − n) .

The third equality, where the integral and summation are exchanged, is

permitted by uniform convergence. So, if T = 1/2Ω, then equation (11) is

identical to equation (1).

6
3.3 An Example

The sampling theory with this band-limited condition, at first, appears too

restrictive for practical use. The formula

f(t) = T ∞

n=−∞

f(nT)

sin( π

T (t − nT))

π(t − nT)

requires that the signal’s frequency be bounded for the reconstruction to be

stable. However, in practice, this condition is not restrictive. The physical

interpretation of an Ω band-limited function is that it is a signal with a

frequency no greater than Ω cycles per second. In countless examples, this

reflects the reality of the physical world. For instance, many sounds are

band-limited naturally; i.e., the frequency of orchestral music is less than

20 kH and the human voice is less than 8 kH. Even our hearing is limited

to a certain range of frequencies, and outside of these frequencies, we miss

the signal. For instance, not many people can hear the frequency of a dog

whistle.

Suppose, one had a signal that was Ω band-limited. Then its Fourier transform would be equivalent to the product of S(ω) and F(ω) = X[−Ω,Ω](ω),

the characteristic function of [−Ω, Ω]. Then the signal is defined by the

convolution of the inverse Fourier transform of S(ω), perhaps s(x), and

f(x) = sin(2πx)

πx . For simplicity, let Ω=1 and the function S(ω)=1, then the

signal is represented just by f(x). (See fig. 1.)

According to the sampling theorem, since f is band-limited, it can be reconstructed from its values at the points nT, where T ≤ 1/2Ω = 1/2 and

n ∈ Z:

g(x) = T ∞

n=−∞

f(nT)

sin( π

T (x − nT))

π(x − nT) ,

where x ∈ R.

In the reconstruction in figure 2, I set T = 1/2Ω = 1/2. Because of the bandlimited condition, there is no gain to oversampling. To further illustrate

the significance of the band-limited condition, i.e., the need to sample at

or above the Nyquist rate, the signal function is undersampled. The next

sample, in figure 3, is taken at T = 1

2∗0.8 ≥ 12Ω, i.e., below Nyquist.

Knowledge of the frequency bound allows one to periodically extend the

transform to the real line. But when the sample is taken below the Nyquist

Figure 1. Original signal

f(x) and its Fourier

transform f(ω).

–5 0 5 –1

0

1

2

Time (x)

f(x)

–2 –1 0 1 2 –1

0

1

2

Frequency (ω)

f(ω)

7
Figure 2. Original signal

f and reconstruction g.

–5 –4 –3 –2 –10 1 2 3 4 5

–0.5

0

0.5

1

1.5

2

2.5

Time (x)

Functions signal values

f(x)

g(x)

Samples

Figure 3. Aliasing

caused by

undersampling.

–5 –4 –3 –2 –10 1 2 3 4 5

–1

–0.5

0

0.5

1

1.5

2

2.5

Time (x)

Functions signal values

f(x)

g(x)

Samples

8
rate, those periodic extensions contain overlapping sections and the information is not weighted appropriately. This effect, demonstrated in figure 3,

is known as aliasing.

3.4 Deriving the Cardinal Series

There are numerous ways in which to derive the cardinal series, or sampling formula, as attested in the literature concerning this. The historical development of its construction follows in section 3.5. In this report,

I display several methods that successfully derive the cardinal series. The

first, which is called the delta method [5], is natural in its approach but is

not rigorous! The second derivation uses the Cauchy integral formula and

shows this to be equivalent to the classical sampling theorem under certain

conditions.

3.4.1 Delta Method

For this derivation, the function is sampled equidistantly. Mathematically,

this is using a series of generalized functions with point support, i.e., delta

functions. Thus, a sampled version of f is

fs(t) = f(t)



n∈Z

δ(t − nτ ) = 

n∈Z

f(nτ )δ(t − nτ ) , (12)

where τ controls the rate of the sampling. This hints at the following reconstruction formula that uses some function g:

f(t) = 

n∈Z

f(nτ )g(t − nτ ) . (13)

The task is now to find an appropriate interpolating function g and rate τ .

Then equation (13) can be rewritten as

f(t) = 

R

f(u)g(t − u)



n∈Z

δ(u − nτ ) du . (14)

The Fourier series expansion of the delta sum is 

n∈Z δ(u − nτ ) =

1/τ 

n∈Z exp−2πinu/τ with the use of equation (12) [1]. So, again equation

(14) is rewritten as

f(t) = 

R

f(u)g(t − u)

1

τ



n∈Z

exp−2πinu/τ du

= f(·)



n∈Z

exp−2πin·/τ ∗

1

τ

g(·)(t) , (15)

where the convolution operator ∗ is defined in section 2. Then, the Fourier

transform of both sides results in

f

(x) = 1

τ

g(x)



n∈Z

f

(x +

n

τ

) . (16)

9
Thus, the transform of f is equivalent to some weighted spectral repetition of itself. To make sense, g must be such that it chooses only one copy

of the spectrum. This is inherently impossible if the repetitions overlap.

Supposing that f is ω band-limited, then 1/τ ≥ 2ω is needed to prevent

overlapping repetitions. This reveals the importance of the rate τ ≤ 1/2ω

in sampling the function. Now, g(x)/τ is a window in which to keep the

desired copy of the spectrum, where n = 0. Thus, define g so that

g(x) = τ · X[− 1

2τ , 12τ ]

(x) , (17)

and note that gneeds only to be equivalent to the characteristic function

over the interval [−ω, ω] ⊂ [−1/2τ, 1/2τ ]. This is equivalent to equation

(17) and keeps the consistency in the formula dependent on the sampling

rate exclusively. This reveals the interpolator to be

g(t) = τ ·

sin 2π 1

2τ t

πt = sin π

τ t

π

τ t , (18)

and therefore, the final rewriting of equation (13) becomes the sampling

formula

f(t) = 

n∈Z

f(nτ )g(t − nτ )

= 

n∈Z

f(nτ )

sin π

τ (t − nτ )

π

τ (t − nτ ) . (19)

Again, the classical sampling formula is derived. This is not the only

method of deriving the classical sampling theorem. A method that uses

complex variables is demonstrated in the next section (3.4.2).

3.4.2 Complex Method

The sampling theorem is intimately connected to a wide variety of mathematical theories. Under certain conditions, it has been shown to be equivalent to a number of other important formulas: the Poisson summation formula, the Cauchy integral formula, etc [5]. Its connection to Cauchy integral theory is demonstrated here.

The relation between complex integration theory and the sampling formula

is inherently obvious. The contour integral around a set of simple poles

is simply a linear combination of some function’s values at the poles (of

course, with the poles removed). This is a kind of sampling relationship,

made obvious when one considers that in the classical sampling theorem,

the “sinc” function contains the appropriate zero values at the sampled

data points of the signal function. This technique fits under the umbrella

of general complex interpolation. The splitting of the integral is called the

Jacobi interpolation formula. This technique was used by Berenstein, Casey,

Yger, and Walnut to solve deconvolution problems.

To show this connection mathematically, I will start with Cauchy’s integral

formula. See Marsden and Hoffman [11] for the necessary background in

complex analysis.

10
Cauchy Integral Formula Theorem [11]. Let f be analytic on a region D and let

C be a closed curve in D that is homotopic to a point, i.e., let C be a simple,

closed, and rectifiable curve in C and also let z ∈ D be a point in D inside

C. Then, for r ∈ N0,

η(C, z) · f(r)

(z) = r!

2πi 

C

f(γ)

(γ − z)r+1 dγ , (20)

where η(C, z) = 1/2πi 

C 1/(γ−z) dγ is the winding number. Note, η(C, z) =

1 if C is positively oriented and loops around z ∈ int C only once.

Given this theorem, I can derive the classical sampling formula. For simplicity, assume that Ω ≤ 1/2T and place the sample values at the integers,

n ∈ Z. Thus, the sampling rate T = 1 and Ω ≤ 1/2. Let f be a function

that satisfies the conditions in the theorem. By the Paley-Weiner theorem

in section 2, f is an analytic function in t ∈ R, which analytically continues to the entire complex plane. Moreover, this entire function satisfies the

Paley-Weiner growth bound, i.e., for Ω > 0,

|f(z)| ≤ A exp−2πΩ|Imz| = A1

(1 + |ζ|)N exp−2πΩ|Imz| . (21)

Now, define μ(ζ) = sin πζ. This has zeros at the integer values, i.e., the

sampling rate. So the contour will need to avoid the zeros of μ(ζ). I define

Cm to be a square contour with corners on

(±1 ± i)(m +

1

2

) .

Thus, η(C, z)=1. Now, first use the Jacobi interpolation formula (see Casey

and Walnut [12]) that gives

f(z) = 1

2πi 

Cm

f(ζ)

ζ − z

dζ

= 1

2πi 

Cm

f(ζ)

ζ − z ·

μ(ζ)

μ(ζ)

dζ (22)

= 1

2πi 

Cm

f(ζ)

ζ − z ·

μ(ζ) − μ(z)

μ(ζ) dζ +

1

2πi 

Cm

f(ζ)

ζ − z ·

μ(z)

μ(ζ)

dζ

= (1) + (2) .

Next, I shall show that (2) → 0 as m → ∞, where

(2) = μ(z)

2πi 

Cm

f(ζ)

(ζ − z)μ(ζ)

dζ . (23)

Recall, that Cm is a square contour with corners (±1±i)(m+1/2). Let R be

a closed, compact disc centered at 0. Then when m + 1/2 > R and ζ ∈ Cm,

|ζ − z| ≥ m + 1/2 − R, such that

|(2)| = |

μ(z)

2πi 

Cm

f(ζ)

(ζ − z)μ(ζ)

dζ|

≤ |μ(z)|

2π



Cm

|f(ζ)|

|(ζ − z)|·|μ(ζ)|

dζ (24)

≤ |μ(z)|

2π(m + 1

2 − R)



Cm

|f(ζ)|

|μ(ζ)|

dζ .

11
Now, let ζ = u + iv. Then, define C(1)

m as the side of Cm that joins (1 −

i)(m + 1/2) and (1 + i)(m + 1/2), C(2)

m that joins (1 + i)(m + 1/2) and (−1 +

i)(m + 1/2), C(3)

m that joins (−1 + i)(m + 1/2) and (−1 − i)(m + 1/2), and

C(4)

m that completes the square contour that joins (−1 − i)(m + 1/2) and

(1 − i)(m + 1/2).

On the side C(1)

m , ζ = m + 1/2 + iv, thus

|sin(πζ)|

2 = sin2(πu) + sinh2(πv) = 1 + sinh2(πv) = cosh2(πv) ≥ exp2π|v| .

Also |f(ζ)| ≤ A exp−2πΩ|v| in equation (21), so that if Ω < 1

2 ,



C(1)

m

|f(ζ)|

|μ(ζ)|

dζ ≤ A

 m+ 1

2

−m− 1

2

exp−π|v|(2Ω+1) dv

= Aexp−π(2Ω+1)(m+ 1

2 ) −1

π(2Ω + 1) , (25)

which disappears as m → ∞, since Ω > 0. However, if Ω=1/2, then



C(1)

m

|f(ζ)|

|μ(ζ)|

dζ ≤

 m+ 1

2

−m− 1

2

A1

(1 + |ζ|)N dζ

=

 m+ 1

2

−m− 1

2

A1

(1 + (m + 1

2 )2 + v2)N dv (26)

= A1 · (2m + 1)

(1 + (m + 1

2 )2)N .

Thus, the integral goes to 0 as m → ∞. A similar approach can show

the same for C(3)

m . On the side C(2)m , ζ = u + i(m + 1/2), and |sin(πζ)| ≥

sinh(π(m + 1/2)). Thus, if Ω ≤ 1/2, then



C(2)

m

|f(ζ)|

|μ(ζ)|

dζ ≤

 m+ 1

2

−m− 1

2

A1

(1 + u2 + (m + 1

2 )2)N

exp−2π|m+ 1

2 |Ω

sinh π(m + 1

2 )

du

= A1 · (2m + 1)

(1 + (m + 1

2 )2)N

exp−2πΩ(m+ 1

2 )

sinh π(m + 1

2 ) , (27)

which also disappears as m → ∞, since Ω > 0. Likewise, C(4)

m is integrated

in the same manner. So, equation (22) becomes

f(z) = (1) = 1

2πi 

Cm

f(ζ)

ζ − z ·

μ(ζ) − μ(z)

μ(ζ) dζ

= μ(z)

2πi 

Cm

f(ζ)

ζ − z ·

	

1

μ(z) − 1μ(ζ)




dζ . (28)

At this point, I need to note a couple items before proceeding. First,

sin πζ = lim

N→∞ sN (ζ) = lim N→∞ πζ 

N

j=1 

1 + ζ2

j2




(29)

12
is the product representation of a sine function. Second, note the following

partial fractions decomposition:

1

(ζ − z)sN (ζ) = 1(ζ − z)sN (z) + 

|n|≤N

1

(ζ − n)(n − z)s


N (n) . (30)

To show this, I need the following: Let p(ζ) be a polynomial with only simple zeros at {rj}s

j=1, then (shown in app A)

1

p(ζ) = s

j=1

1

(ζ − rj )p
(rj ) . (31)

Then, if m<N and defining Rm = (2) on Cm,

f(z) = μ(z)

2πi 

Cm

f(ζ)

ζ − z ·

	

1

μ(z) − 1μ(ζ)




dζ + Rm

= lim

N→∞

sN (z)

2πi 

Cm

f(ζ)

ζ − z ·

	

1

sN (z) − 1sN (ζ)




dζ + Rm

= lim

N→∞

sN (z)

2πi 

Cm

·



|n|≤N

f(ζ)

(ζ − n)(n − z)s


N (n)

dζ + Rm

= lim

N→∞

sN (z)

2πi



|n|≤N

·



Cm

f(ζ)

(ζ − n)(n − z)s


N (n)

dζ + Rm (32)

= lim

N→∞ sN (z)



|n|≤m

f(n)

(z − n)s


N (n) + Rm

= sin πz 

|n|≤m

f(n) · (−1)n

π(z − n) + Rm

= 

|n|≤m

f(n)sinc(z − n) + Rm .

Thus, as m → ∞, then Rm → 0, and this is the classical sampling theorem. After the limit, the contour integral covers over all the zeros of the

sampling function and thus integrates over the whole function. Note, the

switch between limit and integration in the second equality and between

summation and integration in the fourth is permitted because of uniform

convergence. A similar approach is undertaken by Higgins [5], and he also

demonstrates other and more general equivalence relations involving the

classical sampling formula.

3.5 A Bit of History

Sampling theory grew out of a plethora of other mathematical fields. It is

intimately connected with interpolation theory and approximation theory.

In fact, for certain cases, it can be shown that the cardinal series (the classical sampling formula), the Cauchy integral formula, and the Poisson summation formula are equivalent [5]. Sampling theory also has applications

13
in numerous other fields of systems theory, prediction theory, information

theory, stochastic processes, signal and image processing, etc. The theory

of interpolation grew out of the need to calculate intermediate values of

functions with some known values. In 1624, Briggs introduced methods

for achieving this by using successive and modified differences. Wallis also

had some influence in this area in the 1650s and is credited with coining

the term “interpolation.” However, it was not until the work of Gregory

and Newton that the polynomial nature of these methods began to be recognized. In 1670, Gregory introduced the first interpolation series, and later

Newton introduced divided and adjusted divided differences. Both methods derive from taking the polynomial as a close approximation to the (continuous) function to be interpolated [13].

Although these methods were somewhat successful in providing approximations, it became desirable to construct a polynomial interpolant without

the need of finding successive differences. This was achieved with the identity below:

f(λ0)+(t − λ0)f(λ0, λ1) + ··· + (t − λ0)...(t − λM−1)f(λ0, ... ,λM) (33)

= 

M

j=0

f(λj ) GM(t)

G


M(λj )(t − λj ) ,

where

f(λ0, λ1) = f(λ1) − f(λ0)

λ1 − λ0

f(λ0, λ1, λ2) = f(λ1, λ2) − f(λ0, λ1)

λ2 − λ0 ...

and

GM(t) = 

M

j=0



1 − t

λj



.

Both sides are identical because both are polynomials of degree M taking

on the same values on {λj}. The right side of the equality in equation (33)

is known as Lagrange’s formula. Essentially, the Lagrange interpolator acts

as a delta function over the known values. This formula first appeared in

Lagrange’s lectures given in 1795. However, it was previously discovered

by Waring in 1779 [5].

Parallel to the above approach, it became desirable to have a periodic interpolant rather that a polynomial one. Given a periodic function of the form

p(t) = M

j=−M cj expijt, then appropriate formulas were found for these

cases. One such formula is

p(t) = 1

2M + 1



2M

j=0

p



2πj

2M + 1
sin[( 2M+1

2 )(t − 2πj2M+1 )]

sin 1

2 (t − 2πj2M+1 ) .

14
This finite sampling series is credited to Cauchy (1841), although Gauss introduced primitive forms as far back as 1805 [5]. Cauchy’s other significant

contribution was to be the first to note the importance of the rate of the

sampling.

The research conducted in the nineteenth century extended Lagrange’s formula from interpolation on a finite number of points to the more general

case of infinitely many points {λn}. Assuming f(λn) = an for all n ∈ N,

f(z) = 

n

an

φ(z)

φ
(λn)(z − λn) ,

where φ is a function with simple zeros at {λn}. A more general formula

was also constructed:



n

an( z

λn

)

sn φ(z)

φ
(λn)(z − λn) ,

where (sn) is a sequence of integers appropriately chosen to ensure convergence [14]. These formulas were found useful by Cazzaniga in 1882 for constructive purposes. The century ended with the first explicit statement of

the cardinal series by Borel. In 1899, Borel set λn = n and noted that sin(πz)

has the appropriate zeros, returning the classical sampling formula [5].

During the first half of the twentieth century, this special case of the Lagrange formula with knowledge of equidistantly spaced points, i.e., the

cardinal series, was rediscovered in this form several times [5,14]. Whipple

did so in 1910, and E. T. Whittaker in 1915 [10]. Whittaker would also point

out the band-limited nature of the sum. His son, J. M. Whittaker, would

coin the phrase ”cardinal series” in 1929. Ogura, in 1920, was the first to

provide a fully rigorous proof of the cardinal series by using the calculus of

residues. In 1928, Nyquist established the time-bandwidth component of

a signal, the importance of rates, in connection with telegraphy. His work

showed that stability and reconstruction were not possible below a certain sampling rate. In 1941, Hardy noticed that the cardinal series is an

orthogonal expansion, an important development contributing to Hilbert

and Banach spaces. In 1933, Kotel’nikov introduced sampling theory into

communication theory. However, this fact was unknown in the west until

the 1960s. For the Western world, Shannon introduced it in 1949 [14]. Shannon’s original statement of the theorem is ”If a function f(t) contains no

frequencies higher than W cps, it is completely determined by giving its

ordinates at a series of points spaced (1/2 W) s apart.”

In fact, he claims that this fact is common knowledge in communication

and credits Whittaker for an earlier form of the theorem. He also notes

Nyquist’s and Gabor’s use of this fact. However, his name remains attached

to the theorem.

Since Shannon’s use of the sampling theorem, much has been developed in

this field of research. Parzen [15] in 1956 and Petersen and Middleton [16] in

1962 extended the sampling theorem to multiple dimensions, and Kramer

developed a generalization of the classical sampling theorem in 1959 [17].

15
4. Multirate Sampling Theory

The classical sampling theorem shows that one can effectively recover a signal using a regular sampling rate determined by knowledge of the signal

frequency’s bandwidth. However, if the sampling rate is slower than the

Nyquist rate, the function cannot be reconstructed via the classical theorem. Research in this direction has proven effective via other methods. The

procedure here relies on both the classical theorem and deconvolution theory to show how a band-limited function that is sampled below Nyquist

can still be completely recovered.

4.1 Multichannel Deconvolution

Deconvolution has proven itself to be a useful mathematical tool in the field

of signal and image processing. For images, it acts as an enhancing filter to

correct blurs in a picture. For signals, it can be used to correct distorted

line shapes without a loss of signal-noise ratio. In essence, deconvolution

is useful in signal and image processing, where the incoming data contains

a high degree of information. The reason this works so naturally is because

the convolution equation models a number of linear systems.

Before I progress further, it is appropriate to understand what convolution is defined as mathematically. Recall the definition of convolution, also

given earlier (sect. 2), for f,g ∈ L1(R), the convolution of f and g is defined

by

f ∗ g(x) = f(y)g(x − y) dy

for x ∈ R.

For the purpose of this report, the convolution equation s = f ∗ μ models linear, translation invariant systems (e.g., sensors, linear filters). In this

model, f is the input signal function, μ is the system impulse response distribution, and s is the output (received) signal. However, in many physical

applications, s is often a poor approximation of the signal f. This motivates one to deconvolve f from μ to attain the original signal. Results have

shown that if the convolver μ is time-limited (i.e., compactly supported)

and nonsingular (i.e., not a delta generalized function), then this problem

is ill-posed, in the sense of Hadamard [12]. It has been shown to be ill-posed

for all realizable convolvers—all convolvers that can be built. For the circumvention of this scenario, a theory of multichannel deconvolution has

been developed to solve these equations. A multichannel system preserves

information about the signal that would otherwise be lost. Thus, data lost

16
by one convolver can still be retained by another convolver. The signal is

now overdetermined, i.e.,

si = f ∗ μi, i = 1 ...n .

Then, if the convolvers {μi} satisfy the condition of being strongly coprime,

then deconvolving for f is now a well-posed problem.

Strongly coprime. Strongly coprime is defined as a set of convolvers {μi} that

satisfy the inequality

n

i=1

μi|(ζ)|

2


1

2

≥ A exp−B|Imz| ·(1 + |ζ|)

−N (34)

for every ζ ∈ C, where A and B are positive constants and N is a positive

integer, is said to be strongly coprime.

Note, this tells us that 1/n

i=1 |μi(ζ)|

2 satisfies the Paley-Weiner growth

bound. If this is so, then there exists a set of time-limited deconvolvers {νi}

such that

μ1 ∗ ν1 + ··· + μn ∗ νn = δ ,

and consequently,

μ1 · ν1 + ··· + μn · νn = 1 ,

where δ is the Dirac delta function. The second of these two equations is the

analytic Bezout equation. The existence of such deconvolvers is guaranteed

by the following theorem:

H ̈ormander [12]. For compactly supported distributions {μi}n

i=1 on R, compactly supported distributions {νi}n

i=1 exist such that

δ = μ1 ∗ ν1 + ··· + μn ∗ νn .

That is, they satisfy the analytic Bezout equation if and only if the set of

distributions {μi}n

i=1 is strongly coprime.

By Hormander  ̈ ’s theorem, a strongly coprime set has a solution in the analytic Bezout equation. Thus, given the deconvolvers and output signals, f

is naturally produced:



i

si ∗ νi = 

i

(f ∗ μi) ∗ νi = 

i

f ∗ (μi ∗ νi) = f ∗



i

(μi ∗ νi) = f ∗ δ = f .

The system is such that no information is lost in this process. This occurs

because the condition of being strongly coprime guarantees that the zeros

in the analytic Bezout equation do not cluster quickly as |ζ|→∞. If the

{μi} did have a common zero, then si(ζ)=0 at that zero and information

about f would be lost. Because the system is engineered toward eliminating any common zero, no information about f is lost (and the problem is

17
well-posed). Thus, the signal f is gathered by this strongly coprime system,

and these received signals are then filtered by the deconvolvers to reconstruct f.

These methods are linear and realizable; thus deconvolution at a time sample only depends on the information near that time sample. Unfortunately,

Hormander  ̈ ’s theorem is an existence theorem; it does not reveal what the

deconvolvers might be. One obvious solution is

νi = 

μi

i |μi|

2 .

The deconvolvers have been shown not to be unique, and so, for certain

scenarios, one set of deconvolvers may be better than another. More in this

direction is given by Casey and Walnut [18].

4.2 Noncommensurate Sampling Lattices

The theory developed here merges together the two ideas of multichannel deconvolution and classical sampling theory. It has been developed by

Casey [19], Casey and Sadler [20], and Casey and Walnut [12,18].

By the results in multichannel deconvolution, a set of strongly coprime convolvers {μi} needs to be created for the problem to be well-posed. First, a

definition is needed, then the following theorem will be useful in creating

this set.

Definition. A real number α is said to be poorly approximated by rationals

if an integer N ≥ 2 and a constant C = C(α) exist such that for all integers

p, q with q ≥ 2,









α −

p

q









≥ C|q|

−N . (35)

This class of numbers will be denoted by P.

Theorem. Let 0 < r1 < ··· < rm, m ≥ d = 1 satisfy the condition that for

all i 
= j, ri/rj are in P, then {X[−ri,ri]d } is a strongly coprime set. Then to

create a set of appropriate convolvers, one needs only to find such a set

whose ratios are inadequately approximated by rationals.

4.2.1 Two Sampling Lattices

Suppose f ∈ L2(R). Let α be an irrational in the class P. Now, the pair μ1 =

X[−1,1], μ2 = X[−α,α] are strongly coprime. Thus, the problem of solving the

analytic Bezout equation for ν1 and ν2

μ1 · ν1 + μ2 · ν2 = 1

is well-posed. Now, μ1(ω) = sin(2πω)

πω and μ2(ω) = sin(2παω)πω have zeros

Zμ1 =

n

2



, Zμ2 =

 n

2α



18
for n ∈ Z\{0}. Thus, one solution to the analytic Bezout equation is ν1(ω) =

1/2μ1(ω) for ω ∈ Zμ1 and, likewise, ν2(ω)=1/2μ2(ω) for ω ∈ Zμ2 . Then, the

problem simply becomes an interpolation problem on the set of zeros Γ =

Zμ1

 Zμ2 . So, one has a new reconstruction formula for f from its values

on Γ. However, this can be done for any poorly approximated rational. The

following theorem results:

Theorem. Let

Γ = ±k

2

±k

2α



(36)

for n ∈ N, and let λi ∈ Γ. Let f be a (1 + α)-band-limited function, then f

can be conditionally reconstructed by the formula

f(t) = 

λi∈Γ

f(λi) G(t)

G
(λi)(t − λi) + [f(0)K1(t) + f


(0)K2(t)] , (37)

where

G(t) = sin(2πt) · sin(2παt) (38)

K1(t) = G(t)

G(0)

2! t2 (39)

K2(t) = G(t)

G(0)

2! t . (40)

Note that the interpolators at the origin appropriately model the function at

the origin, i.e., K1(0) = 1 while K


1(0) = 0 and K2(0) = 0 while K
2(0) = 1.

How these interpolators are calculated will be demonstrated later in section 4.3.1. Also, note that the information in the signal can be reconstructed

uniquely by its samples on Γ∪{0}. Here, it is also important to note that

the sampling rates correspond to 1-band-limited functions and α-bandlimited functions. These add up to the band-limit of the function to be sampled. Thus, a (1+α)-band-limited function is reconstructed via significantly

lower sampling rates.

4.2.2 An Arbitrary Number of Lattices

This result generalizes for an arbitrary number of sampling rates. Let

{ri}n

i=1 be rates such that all ri/rj , i 
= j, are in P, the class of irrationals poorly approximated by rationals. Then, {ri} is a strongly coprime

set. Thus, our convolvers have the form μi(t) = X[−ri,ri](t). These {μi}

model the impulse response of a multichannel system. The convolvers have

Fourier transforms

μi(ζ) = sin(2πriζ)

πζ

with zero sets Zi = { ±k

2ri

}, both respectively, for k ∈ N. Note that exclusive

of the origin, the zeros sets are nonrepetitive. Let Γi = Zi and Γ = n

i=1 Γi =

19
{λi}. Thus, f can be reconstructed on its values at Γ. This can also be done

via techniques in complex interpolation theory.

Thus, the following theorem results: Let {ri}n

i=1 be a set whose ratios are

poorly approximated by rationals, and let f be a (

ri)-band-limited function. Let

Γi =

±k

2ri



for k ∈ N and i = 1,... ,n. Also, let

n

i=1

Γi =Γ= {λi} .

Then f is uniquely determined by

{f(λi)}  

f(0),... ,f(n−1)(0).

Furthermore, f can be reconstructed from its values on Γ {0} by the following formula:

f(t) ≈ 

λi∈Γ

f(λi) G(t)

G
(λi)(t − λi) +n

i=1

f(i−1)(0) · Ki(0) , (41)

where

G(t) = n

i=1

sin(2πrit) (42)

and the interpolating functions Ki(t) at the origin are a linear combination

of G(t)/tj , j = 1, . . . , n, chosen so that K(l−1)

i (0) = δi,l, i, l = 1,... ,n.

4.3 Interpolators at Origin

What remains for this interpolating function in equation (41) is to find appropriate interpolators at the origin. To remove this ambiguity, I construct

interpolators here. Under the assumption that there are two interpolators,

then G(z) = sin(2πz) · sin(2παz), where α is the square root of any prime

number. What is desired, as stated in the theorem in section 4.2.2, is for

the interpolators to satisfy the following condition: {Ki(t)} is a linear combination of G(t)/ti

, i = 1,... ,n, chosen so that K(j−1)

i (0) = δi,j , where

j = 1,... ,n.

4.3.1 Two Interpolators

Assume that there are two interpolators. As stated earlier in equation (42),

G(t) = sin(2πt) · sin(2παt) .

20
To find the appropriate interpolators at the origin, I need to manipulate this

function. Remember that the Taylor expansion formula of a sine function is

sin(t) = t − t

3

3! + O(t

5) . (43)

Thus, inserting equation (43) into G(t) obtains

G(t) = sin(2πt) · sin(2παt)

=

	

2πt − (2πt)3

3! + O(t

5)




·

	

2παt − (2παt)3

3! + O(t

5)




= 4π2αt2 −

	

16π4α3

3! +

16π4α

3! 


t

4 + O(t6) .

And so, the following are representations for the linear combinations of

G(t)/ti

, i = 1,... ,n:

G(t)

t

= 4π2αt −

	

16π4α3

3! +

16π4α

3! 


t

3 + O(t5)

G(t)

t2 = 4π2α −

	

16π4α3

3! +

16π4α

3! 


t

2 + O(t4) .

Their respective derivatives follow:

d

dt

G(t)

t

= 4π2α − 3

	

16π4α3

3! +

16π4α

3! 


t

2 + O(t4)

d

dt

G(t)

t2 = 2	

16π4α3

3! +

16π4α

3! 


t + O(t

3) .

Next, note the limits of the functions and their derivatives at t = 0 —

lim

t→0

G(t)

t

= 0

lim

t→0

G(t)

t2 = 4π2α

lim

t→0

d

dt

G(t)

t

= 4π2α

lim

t→0

d

dt

G(t)

t2 = 0 .

The first interpolant K1 is desired to be such that K1(0) = 1 and K


1(0) = 0.

Also, K2 is desired such that K2(0) = 0 and K


2(0) = 1. Setting the following relations satisfy these requirements:

K1(t) = G(t)

G(0)

2! t2 , K2(t) = G(t)

G(0)

2! t , (44)

where, conveniently, one can see that G

(0)/2! = 4π2α.

21
4.3.2 Arbitrary Number of Interpolators

It is convenient to further construct this idea to an arbitrary multiple of interpolators. In this scenario, assume that there are n interpolators involved.

Thus, as the theorem in section 4.2.2 dictates, there are n sampling rates

{ri} = {

√p0 = 1,√p1,... , √pn−1}, where p1,... ,pn−1 are the first n − 1

primes, which generate the following reconstruction interpolant:

G(t) = n

i=1

sin(2πrit) . (45)

Now, the process of finding the interpolants at the origin continues as before. With the use of the Taylor expansion formula for the sine function,

then

G(t) = n

i=1 	

2πrit − 8π3r3

i

3! t

3 + O(t5)




= n

i=1 	

2πrit − 8π3pi−1 · ri

3! t

3 + O(t5)




= 2nπn

n

i=1

ri




t

n − 2n+2πn+2

3! n

i=1

ri


1 +

n

−1

i=1

pi




t

n+2 + O(tn+4) (46)

= 2nπn√p1p2 ··· pnt

n − 2n+2πn+2

3!

√p1p2 ··· pn



1 +

n

−1

i=1

pi




t

n+2 + O(tn+4) .

Now, the interpolants at the origin {Ki} again will be linear combinations

of G(t)/ti

, i = 1,... ,n, chosen appropriately so that K(j−1)

i (0) = δi,j ,where

j = 1,... ,n. The general form of G(t)/ti is

G(t)

tk = 2nπn√p1p2 ··· pnt

n−k − 2n+2πn+2

3!

√p1p2 ··· pn



1 +

n

−1

i=1

pi




t

n−k+2 + O(tn−k+4) . (47)

Let Hk(t) = G(t)/tk. Then, the formula for the jth derivative of Hk is

H(j)

k (t)=2nπn (n − k)!

(n − k − j)!

√p1p2 ··· pnt

n−k−j

−2n+2πn+2

3!

(n − k + 2)!

(n − k − j + 2)!

√p1p2 ··· pn



1 +

n

−1

i=1

pi




t

n−k−j+2 + O(tn−k−j+4). (48)

Allowing t to approach 0, the following results are obtained:

lim

t→0

Hk(t)=2nπn√p1p2 ··· pnδn,k

= G(n)

(0)

n! , (49)

22
and

lim

t→0

H(j)

k (t)=2nπn (n − k)!

(n − k − j)!

√p1p2 ··· pnδn,k+j (50)

−2n+2πn+2

3!

(n − k + 2)!

(n − k − j + 2)!

√p1p2 ··· pn



1 +

n

−1

i=1

pi




δn,k+j−2 .

Now, using these evaluations in equation (50) as before, I can construct a

formula of appropriate interpolants using the set {Hk}. First, note that each

interpolant is of the form

Kk(t) = 

k

i=1

ciHi(t)

for all k = 1,... ,n. Note that the ci’s vary with the k’s. As required of the

sampling theorem, coefficients ci are needed such that K(j)

k (0) = δk−1,j . By

linearity of differentiation, I have

K(j)

k (0) = 

k

i=1

ciH(j)

i (0)

= 

k

i=1

ci ·

	

2nπn (n − i)!

(n − i − j)!

√p1p2 ··· pnδn,i+j

−2n+2πn+2

3!

(n − i + 2)!

(n − i − j + 2)!

√p1p2 ··· pn



1 +

n

−1

i=1

pi




δn,i+j−2




= 

k

i=1 	

ci2nπn (n − i)!

(n − i − j)!

√p1p2 ··· pn

−ci+2

2n+2πn+2

3!

(n − i)!

(n − i − j)!

√p1p2 ··· pn



1 +

n

−1

i=1

pi



δn,i+j (51)

= 

k

i=1

2nπn (n − i)!

(n − i − j)!

√p1p2 ··· pn

	

ci − ci+2

22π2

3! 

1 +

n

−1

i=1

pi



δn,i+j

= 2nπn√p1p2 ··· pn · j!

	

cn−j − cn−j+2

22π2

3! 

1 +

n

−1

i=1

pi





= G(n)

(0)

(n)! · j!

	

cn−j − cn−j+2

22π2

3! 

1 +

n

−1

i=1

pi





= δk−1,j .

In equation (51), ck+1, ck+2 = 0 is defined so the summation evaluates correctly. Now, remember the following relation must hold: K(k−1)

k (0) = 1,

23
since

K(k−1)

k (0) = G(n)

(0)

(n)! · (k − 1)!	

cn−k+1 − cn−k+3

22π2

3! 

1 +

n

−1

i=1

pi



 . (52)

Since K(j)

k (0) = 0 if j 
= k − 1, then

	

cn−j − cn−j+2

22π2

3! 

1 +

n

−1

i=1

pi



 = 0 . (53)

Thus, via cancelation, cn−k+1 = (G(n)

(0)/(n)! · (k − 1)!)−1. Note, also from

the summations just mentioned, that there is the restriction n − j ≤ k, or

restated as j ≥ n − k. So, continuing in this direction, I find the following

recursive equations:

cn−k+1+2m = [

22π2

3! (1 + n−1

i=1 pi)]m

G(n)(0)

(n)! · (k − 1)!

(54)

cn−k+2+2m = 0 . (55)

This generates interpolants at the origin with the desired properties. Thus,

the interpolants can be formulated as

Kn−k+1(t) =

n

−k

i=1

ciHi(t)

=

(n−k)/2


m=0 

22π2

3! 

1 + n−1

i=1 pi

m

t

2mG(t)

G(n)(0)

(n)! · (k − 1)! · tn−k+1 . (56)

Now, with these interpolants, I can use this new reconstruction method

to recover a function of a larger bandwidth than our individual sampling

rates.

4.4 Simulations

Here I present a number of simulations of the reconstructions of a signal

using the multirate sampling formula. The same function (signal) used in

section 3.3, f(t) = sin(2πt)/πt. Recall that this is an Ω=1 band-limited

function. The classical sampling theory stipulates that the sampling rate be

T ≤ 1/2Ω = 1/2. The classic theory dictates that the Nyquist rate condition holds to successfully reconstruct the function. As shown earlier, it was

sufficient for the signal to be sampled at Nyquist. Now, I wish to sample

slower than the Nyquist rate and still reconstruct the function. First, the

signal will be sampled twice below Nyquist at the sampling values

Γ1 =







±k

2 1

(1+√2)





 , Γ2 =







±k

2

√2

(1+√2)





 .

24
Note that these sample rates are both slower than the Nyquist rate. Also,

note that these are appropriately chosen, i.e., 1/(1+√2)+√2/(1+√2) = 1,

and the ratio of the two samples is √2 ∈ P. Recall, that P is the class of

irrationals poorly approximated by rationals. Thus, by either theorem in

section 4.2.2 or 4.2.3, the signal can be reconstructed by

g(t) = 

λi∈Γ

f(λi) G(t)

G
(λi)(t − λi) + [f(0)K1(t) + f


(0)K2(t)] ,

with

G(t) = sin  2πt

(1 + √2)

· sin 

2π

√2t

(1 + √2)


K1(t) = G(t)

G(0)

2! t2

K2(t) = G(t)

G(0)

2! t .

Essentially, G(t) looks exactly as it should, i.e,. the product of two sine functions. The reconstruction of the signal is in figure (4). It gives the appearance of completely reconstructing the function signal; however, upon closer

analysis, minor variations can be seen between the orignal signal and the

reconstruction. These points of aliasing are due to locations in the sampling

grid where the samples are close. Methods of eliminating this error are discussed by Casey and Walnut [12,19]. The best way to improve this sampling

scheme is to back off the bandwidth slightly, i.e., increase the sampling rates

by about 10 percent or more. Then, the signal can be reconstructed, practically, within an acceptable margin of error.

Figure 4. Original signal

f(t) and bisample

reconstruction g(t).

–5 –4 –3 –2 –1 0 1 2 3 4 5

–1

–0.5

0

0.5

1

1.5

2

2.5

Time (x)

f(x)

g(x)

Samples

Functions f & 

g

25
Again, as in the classical formula, the rate is critically important in achieving an accurate reconstruction of the signal. Research on the bounds of

these rates is ongoing.

Two items of interest are the interpolator and its derivative and, most importantly, the reconstruction of the signal. The interpolator and its derivative are displayed in figure 5.

In figure 6, I demonstrate the reconstruction of the same signal but use three

noncommensurate sampling rates below Nyquist. Here, I use the following

sampling rates:

Γ1 =







±k

2 1

(1+√2+√3)





 , Γ2 =







±k

2

√2

(1+√2+√3)





 , Γ3 =







±k

2

√3

(1+√2+√3)





 .

Figure 5. (a) Interpolator

and (b) its derivative.

–5 –4 –3 –2 –1 0 1 2 3 4 5 –1

–0.5

0

0.5

1

–5 –4 –3 –2 –1 0 1 2 3 4 5 –4

–2

0

2

4

G′(t) G(t) (a)

(b)

Figure 6. Original signal

f(t) and trisample

reconstruction g(t).

–5 –4 –3 –2 –10 1 2 3 4 5

–1

–0.5

0

0.5

1

1.5

2

2.5

Time (x)

f(x)

g(x)

Samples

Signals f & 

g

26
Note that these rates are still below Nyquist. The reconstruction formula

used is the same as presented in section 4.2.2. The only real point of interest

to present is the interpolators at the origin. From the derivations in section

4.3.2, the set of interpolators is shown to be

K1(t) = G(t) + 4π2

6 · (1 + √2 + √3) · t

2 · G(t)

G(3)(0)

3! t3

K2(t) = G(t)

G(3)(0)

3! t2

K3(t) = G(t)

G(3)(0)

3! t

.

A number of other simulations were run with an encouraging success rate,

demonstrating the viability of this multisampling technique. The success

rate is somewhat improved when the sampling rates are more equally

spaced. A brief discussion of applications follows in the next section.

4.5 Applications

A benefit to using this technique is that the individual sampling rates are

lower than the traditional Nyquist sampling rate in the classical sampling

theorem. For example, for a (1+α)-band-limited function, the Nyquist rate

is 1/2(1+α), but when sampling with two lattices, the sampling rates are

only 1/2 and 1/2α. Both these results are below Nyquist. The generalization shows that the comparisons between the sampling rates and the frequency band limit are more outstanding as the number of lattices increase.

Suppose that I have an Ω band-limited signal, but the sampling rate is restricted to 1/2β, then choose the first n primes so that

(1 + p1 + p2 + ··· + pn) · β

pn

≥ Ω . (57)

Then, n + 1 sampling lattices can be constructed via the theorem in section

4.2.3 that can effectively reconstruct f with sampling rates no faster than

1/2β.

Unfortunately, this structure has some difficulties. The first is that the sampling grid is extremely rigid; perturbations on the lattices result in losing

information. Also, because the sampling points might become arbitrarily

close, ripples occur at points near each other. If the bandwidth is backed

off, the simulations results improve, but exact bounds are unknown.

These results also extend to higher dimensions, where one can think of

the grids and the interpolating function as Cartesian products. Casey and

Sadler [20] have used these techniques to develop new analog-digital transformers for signal processing. In their signal processing literature, these are

known as A-D converters.

27
References

1. J. J. Benedetto, Harmonic Analysis and Applications, CRC Press, Boca

Raton, Florida (1997).

2. H. Dym and H. P. McKean, Fourier Series and Integrals, Academic

Press, New York (1972).

3. T. W. Korner,  ̈ Fourier Analysis, Cambridge University Press, Cambridge, Great Britain (1988).

4. T. M. Apostol, Mathematical Analysis, Addison–Wesley, Reading, Massachusetts (1974).

5. J. R. Higgins, Sampling Theory in Fourier and Signal Analysis: Foundations, Clarendon Press, Oxford (1996).

6. Robert J. Marks II, Introduction to Shannon Sampling and Interpolation

Theory, Springer-Verlag, New York (1991).

7. A. I. Zayed, Advances in Shannon’s Sampling Theory, CRC Press, Boca

Raton, Florida (1993).

8. H. Nyquist, “Certain topics in telegraph transmission theory,” AIEE

Trans. 47 (1928), pp 617–644.

9. C. Shannon, “A mathematical theory of communication,” Bell System

Technical Journal 27 (1948), pp 379, 623.

10. E. T. Whittaker, “On the functions which are represented by the expansions of the interpolation theory,” Proc. Royal Society of Edinburg

35 (1915), pp 181–194.

11. J. E. Marsden and M. J. Hoffman, Basic Complex Analysis, W. H. Freeman, New York (1987).

12. S. Casey and D. Walnut, “Systems of convolution equations, deconvolution, shannon sampling and the gabor and wavelet transform,”

SIAM Review 36, No. 4 (1994), pp 537–577.

13. D. T. Whiteside, “Patterns of mathematical thought in the later seventeenth century,” Arch. Hist. Exact Sci. 1 (1961), pp 178–388.

14. A. J. Jerri, “The Shannon sampling theorem—Its various extensions

and applications: A tutorial review,” Proc. IEEE 65, No. 11 (1977),

pp 1565–1596.

28
15. E. Parzen, A Simple Proof and Some Extensions of Sampling Theorems,

Stanford, CA, Technical Report 7 (1956).

16. D. P. Petersen and D. Middleton, “Sampling and reconstruction of

wave number–limited function in N-dimensional Euclidean spaces,”

Informat. Contr. 5 (1962), pp 279–323.

17. H. P. Kramer, “A generalized sampling theorem,” J. Math. Phys. 38

(1959), pp 68–72.

18. S. Casey and D. Walnut, “Residue and sampling techniques in deconvolution,” 1997 International Workshop on Sampling Theory and Applications (SampTA ’97) (1997), pp 279–284.

19. S. Casey, “Sampling on unions of noncommensurate lattices via complex interpolation theory,” 1999 International Workshop on Sampling

Theory and Applications (1999), pp 48–53.

20. S. Casey and B. Sadler, New Directions in Sampling Multi-Rate A-D

Conversion via Number Theoretic Methods, IEEE International Conference on Acoustics, Speech, and Signal Processing 1 (2000), pp 336–339.

29

Appendix A. Proof of the Polynomial Expansion

I shall show the result in section 3.4.2 of the main report via induction on

the degree of the polynomial.

Proof. Suppose p(ζ) is a one degree polynomial with leading coefficient 1

and one simple zero, then p(ζ) = ζ − r1. Thus, clearly,

1

p(ζ) = 1ζ − r1

. (A-1)

Now, assuming that this statement is true for any k degree polynomial with

leading coefficient 1 and k distinct simple zeros, where k = 1,... ,s, I want

to show the same result for an s + 1 degree polynomial with s simple zeros.

Suppose p(ζ) is such a polynomial. Then, I have a representation for p using

the set of its zeros:

p(ζ) = s

j=0

(ζ − rj ) . (A-2)

I can equivalently state this as p(ζ)=(ζ − r0) · q(z), where q is an s degree polynomial with only s simple zeros. Thus, I have a formula by the

assumption in the induction for q. Using this, I obtain the following:

1

p(ζ) = 1(ζ − r0) · 1q(z)

= 1

(ζ − r0) ·

s

j=1

1

(ζ − rj ) · q
(rj ) (A-3)

= s

j=1

1

(ζ − r0)(ζ − rj )q
(rj ) .

31
Now, let m(ζ)=(ζ − r0)(ζ − rj )q


(rj ), then m
(ζ)=(ζ − rj )q
(rj )+(ζ −

r0)q


(rj ). Also, note that p
(ζ) = q(ζ)+(ζ − r0)q
(ζ). So continuing, I have

1

p(ζ) = s

j=1

1

(ζ − r0)(ζ − rj )q
(rj )

= s

j=1

1

(ζ − r0)m
(r0) +

1

(ζ − rj )m
(rj )

= s

j=1

1

(ζ − r0)(r0 − rj )q
(rj ) +

1

(ζ − rj )(rj − r0)q
(rj ) (A-4)

= 1

(ζ − r0)

s

j=1

1

(r0 − rj )q
(rj ) +s

j=1

1

(ζ − rj )p
(rj )

= 1

(ζ − r0) · 1q(r0) +s

j=1

1

(ζ − rj )p
(rj )

= 1

(ζ − r0) · 1p
(r0) +s

j=1

1

(ζ − rj )p
(rj ) .

Equation A-4 gives the desired result.

32
Appendix B. MATLAB Code

The following MATLAB *.m files were used to perform the simulations. I

used MATLAB version 5.3 by The Mathworks, Inc.

B-1 Shannon.m

%SHANNON EXAMPLE % This m-file generates plots for the demonstration of

Shannon's formula. Three plots will be created. The first will show the signal and its

Fourier transform, the second will demonstrate the reconstruction procedure, and

the last plot will show the effects of aliasing in the sampling.

% Terrence Moore, 3/6/2000

% Revision: 1.3 Date: 7/5/2000 18:59:33

%begin of code

%generate the domain

x = linspace(-5,5,1000);

%the signal function

y = sin(2*pi*x)./(pi*x);

%plot of signal

figure;

subplot(3,2,3);

plot(x,y);

title('original signal f(x)');

xlabel('x');

ylabel('f(x)');

%the transform of the signal

x1 = [-2 -1 -1 1 1 2];

y1 = [0 0 1 1 0 0];

%plot of the transform

subplot(3,2,4);

plot(x1,y1);

axis([-2 2 -1 2]);

title('Fourier transform f(\ omega)');

xlabel('\ omega');

ylabel('f(\ omega)');

33
%our recontruction formula only for this signal

%omega = 1

s=1;

z=0;

count=0;

for n=-20:20

if n=0

z=z+sin(pi*n/s)./(pi*n/(2*s)).*sin(pi*(2*s*x-n))./(pi*(2*s*x-n));

count=count+1;

end

end

z=z+2*sin(pi*(2*s*x))./(pi*(2*s*x));

figure;

plot(x,y,'-',x,z,':');

hold;

%the sampled data

N = -20:20;

N = N/(2*s);

yN = sin(2*pi*N)./(pi*N);

plot(N,yN,'o');plot(0,2,'o');

axis([-5 5 -0.5 2.5])

xlabel('x');

legend('f(x)','g(x)','samples');

%omega = 0.8 (undersampling - aliasing)

s=0.8;

z=0;

count=0;

for n=-20:20

if n =0

z=z+sin(pi*n/s)./(pi*n/(2*s)).*sin(pi*(2*s*x-n))./(pi*(2*s*x-n));

count=count+1;

end

end

z=z+2*sin(pi*(2*s*x))./(pi*(2*s*x));

figure;

plot(x,y,'-',x,z,':');

hold;

%the sampled data

N = -20:20;

N = N/(2*s);

yN = sin(2*pi*N)./(pi*N);

34
plot(N,yN,'o');plot(0,2,'o');

axis([-5 5 -1 2.5])

xlabel('x');

legend('f(x)','g(x)','samples');

%omega = 1.2 (oversampling -- aliasing)

s=1.2;

z=0;

count=0;

for n=-20:20

if n =0

z=z+sin(pi*n/s)./(pi*n/(2*s)).*sin(pi*(2*s*x-n))./(pi*(2*s*x-n));

count=count+1;

end

end

z=z+2*sin(pi*(2*s*x))./(pi*(2*s*x));

figure;

plot(x,y,'-',x,z,':');

hold;

%the sampled data

N = -20:20;

N = N/(2*s);

yN = sin(2*pi*N)./(pi*N);

plot(N,yN,'o');plot(0,2,'o');

axis([-5 5 -0.5 2.5])

title('reconstruction function g(x) and signal f(x)');

xlabel('x');

35
B-2 MultirateShannon.m

t = linspace(-5,5,1000);

y = sin(2*pi*t)./(pi*t);

r0 = 1 + sqrt(2);

r = [1 sqrt(2)];

r = r./r0;

N = 20;

rates = 1./(2*r);

Ns = round(N*(1./rates) +1);

tsamp = zeros(length(r),2*max(Ns)+1);

for i = 1:length(r)

tsamp(i,:) = [(-1*Ns(i):Ns(i)) zeros(1,2*max(Ns)-2*Ns(i))];

end

samp = diag(rates)*tsamp;

samp = samp(:);

indz = find(samp == 0);

samp(indz) = [];

lambda = sort([samp.']);

lambda = lambda;% + .001;

ST = sin(2*pi*r.'*t);

St = prod(ST,1);

CT = 2*pi*diag(r)*cos(2*pi*r.'*t);

CT = flipud(CT);

Sp = sum(ST.*CT,1);

figure;

subplot(2,1,1);

plot(t,St);

subplot(2,1,2);

plot(t,Sp);

Sl = sin(2*pi*r.'*lambda);

Cl = 2*pi*diag(r)*cos(2*pi*r.'*lambda);

Cl = flipud(Cl);

bp = Sl.*Cl;

bp = sum(bp,1);

%sample values to be interpolated

tp = sin(2*pi*lambda)./(pi*lambda);

fp = tp./bp;

36
temp = zeros(size(t));

for i = 1:length(lambda)

temp = temp + fp(i)./(t-lambda(i));

end

f = temp.*St;

ST2 = 2*4*pi*pi*prod(r,2);

K1 = St./(t.*t);

K1 = K1/ST2*2;

K2 = St./(t);

K2 = K2/ST2*2;

f=f+ 2*K1 + 0*K2;

figure;

plot(t,y,'-',t,f,':');

hold;

%the sampled data

N = lambda;

yN = sin(2*pi*N)./(pi*N);

plot(N,yN,'o');plot(0,2,'o');

axis([-5 5 -1 2.5])

xlabel('x');

legend('f(x)','g(x)','samples');

37

39

Distribution

Admnstr

Defns Techl Info Ctr

ATTN DTIC-OCP

8725 John J Kingman Rd Ste 0944

FT Belvoir VA 22060-6218

DARPA

ATTN S Welby

3701 N Fairfax Dr

Arlington VA 22203-1714

Ofc of the Secy of Defns

ATTN ODDRE (R&AT)

The Pentagon

Washington DC 20301-3080

Ofc of the Secy of Defns

ATTN OUSD(A&T)/ODDR&E(R) R J Trew

3080 Defense Pentagon

Washington DC 20301-7100

AMCOM MRDEC

ATTN AMSMI-RD W C McCorkle

Redstone Arsenal AL 35898-5240

Army Rsrch Physics Div

ATTN AMSRL-RO-EM Assoc Dir Math &

Computer Sci Div R Launer

PO Box 12211

Research Triangle Park NC 27709-2211

US Army TRADOC

Battle Lab Integration & Techl Dirctrt

ATTN ATCD-B

FT Monroe VA 23651-5850

US Military Acdmy

Mathematical Sci Ctr of Excellence

ATTN MADN-MATH MAJ M Huber

Thayer Hall

West Point NY 10996-1786

Dir for MANPRINT

Ofc of the Deputy Chief of Staff for Prsnnl

ATTN J Hiller

The Pentagon Rm 2C733

Washington DC 20301-0300

SMC/CZA

2435 Vela Way Ste 1613

El Segundo CA 90245-5500

TECOM

ATTN AMSTE-CL

Aberdeen Proving Ground MD 21005-5057

US Army ARDEC

ATTN AMSTA-AR-TD

Bldg 1

Picatinny Arsenal NJ 07806-5000

US Army Info Sys Engrg Cmnd

ATTN AMSEL-IE-TD F Jenia

FT Huachuca AZ 85613-5300

US Army Natick RDEC Acting Techl Dir

ATTN SBCN-T P Brandler

Natick MA 01760-5002

US Army Simulation Train & Instrmntn

Cmnd

ATTN AMSTI-CG M Macedonia

ATTN J Stahl

12350 Research Parkway

Orlando FL 32826-3726

US Army Tank-Automtv Cmnd RDEC

ATTN AMSTA-TR J Chapin

Warren MI 48397-5000

Nav Surfc Warfare Ctr

ATTN Code B07 J Pennella

17320 Dahlgren Rd Bldg 1470 Rm 1101

Dahlgren VA 22448-5100

The American Univ

Dept of Mathematical Statistics

ATTN Nolan

ATTN S D Casey

4400 Massachusets Ave NW

Washington DC 20016-8050

Hicks & Assoc Inc

ATTN G Singley III

1710 Goodrich Dr Ste 1300

McLean VA 22102
Distribution (cont’d)

40

Director

US Army Rsrch Lab

ATTN AMSRL-RO-D JCI Chang

ATTN AMSRL-RO-EN W D Bach

PO Box 12211

Research Triangle Park NC 27709

US Army Rsrch Ofc

ATTN AMSRL-RO-EM J Lavery

ATTN AMSRL-RO-EM J Wu

ATTN AMSRL-RO-EM S Davis

ATTN AMSRL-RO-EM B J West

PO Box 12211

Research Triangle Park NC 27709-2211

US Army Rsrch Lab

ATTN AMSRL-D D R Smith

ATTN AMSRL-DD J M Miller

ATTN AMSRL-CI-AI-R Mail & Records

Mgmt

ATTN AMSRL-CI-AP Techl Pub (2 copies)

ATTN AMSRL-CI-LL Techl Lib (2 copies)

ATTN AMSRL-CI-C J Gowens

ATTN AMSRL-CI-CN A Swami

ATTN AMSRL-CI-CN B Sadler

ATTN AMSRL-CI-CN T Moore (15 copies)

Adelphi MD 20783-1197
1. AGENCY USE ONLY

8. PERFORMING ORGANIZATION

 REPORT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

12a. DISTRIBUTION/AVAILABILITY STATEMENT

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHOR(S)

REPORT DOCUMENTATION PAGE

2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

11. SUPPLEMENTARY NOTES

14. SUBJECT TERMS

13. ABSTRACT (Maximum 200 words)

Form Approved

OMB No. 0704-0188

(Leave blank)

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

12b. DISTRIBUTION CODE

15. NUMBER OF PAGES

16. PRICE CODE

17. SECURITY CLASSIFICATION

 OF REPORT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

19. SECURITY CLASSIFICATION

 OF ABSTRACT 20. LIMITATION OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Shannon-Type Sampling Theory on Unions of Equally

Spaced and Noncommensurate Grids

February 2001 Final, 10/99–9/00

Sampling theory is that branch of mathematics that seeks to reconstruct functions from its values

at a discrete set of points. The fundamental result in sampling theory known as “Shannon’s

sampling theorem” has many applications to signal processing and communications engineering. I

demonstrate Shannon’s result via complex interpolation methods. I then quote a result that uses

these methods to solve interpolation problems on unions of noncommensurate lattices, which are

created via a specific number of theoretic guidelines. These interpolations give Shannon-type

reconstructions on these lattices. I close by doing simulations in MATLAB of the sampling

reconstructions on these noncommensurate grids.

Signal processing, Shannon’s Sampling Theorem, multichannel

deconvolution sampling lattices, interpolation, Cauchy Integral Formula,

signal reconstruction, multirate sampling

Unclassified

ARL-TR-2305

1FEPC1

611102H4811

AH48

611102H48

ARL PR:

AMS code:

DA PR:

PE:

Approved for public release; distribution

unlimited.

45

Unclassified Unclassified SAR

2800 Powder Mill Road

Adelphi, MD 20783-1197

U.S. Army Research Laboratory

Attn: AMSRL-CI-CN tm1917a@american.edu

U.S. Army Research Laboratory

41

Terrence J. Moore

email:

2800 Powder Mill Road

Adelphi, MD 20783-1197