1

Optimal High-Resolution Adaptive Sampling
of Deterministic Signals

arXiv:1611.01850v3 [] 20 Jun 2017

Yehuda Dar and Alfred M. Bruckstein

Abstract—In this work we study the topic of high-resolution
adaptive sampling of a given deterministic signal and establish
a connection with classic approaches to high-rate quantization.
Specifically, we formulate solutions for the task of optimal highresolution sampling, counterparts of well-known results for highrate quantization. Our results reveal that the optimal highresolution sampling structure is determined by the density of the
signal-gradient energy, just as the probability-density-function
defines the optimal high-rate quantization form. This paper has
three main contributions: the first is establishing a fundamental
paradigm bridging the topics of sampling and quantization.
The second is a theoretical analysis of sampling, for arbitrary
signal-dimension, relevant to the emerging field of high-resolution
signal processing. The third is a new approach to nonuniform
sampling of one-dimensional signals that is experimentally shown
to outperform an optimized tree-structured sampling technique.
Index Terms—High-resolution sampling, adaptive sampling,
high-rate quantization.

I. I NTRODUCTION
AMPLING and quantization are fundamental processes
in signal digitization and coding techniques. Each of
them is a field of research rich in theoretical and practical
studies. Quantization addresses the problem of discretizing
a range of values by a mapping function, usually based on
decomposition of the range into a finite set of non-intersecting
regions. Sampling a given signal, defined over a continuous
and bounded domain, is the task of discretizing the signal
representation (for example, representing a finite-length onedimensional signal as a vector). In this paper we consider
nonuniform sampling that relies on segmentation of the signal
domain into non-overlapping regions, each represented by a
single scalar coefficient (we do not consider here generalized
sampling that uses projections of the signal onto a discrete set
of orthonormal functions). The signals considered in this work
are deterministic in the sense that their continuous forms are
accessible to the sampler for its operation. Our main settings
address a given one-dimensional signal defined over the interval [0, 1), and its representation using a piecewise-constant
approximation that relies on a high-resolution nonuniform
segmentation (see Fig. 1).
The discretizations in quantization and sampling were first
implemented in their simplest forms relying on uniform divisions of the respective domains. Then, the quantizer designs progressed to utilize nonuniform structures, exploiting
input-data statistics to improve rate-distortion performance. In
his fundamental work, Bennett [1] suggested to implement

S

The authors are with the Department of Computer Science, Technion, Israel.
E-mail addresses: {ydar, freddy}@cs.technion.ac.il.

Fig. 1. The general procedure for nonuniform sampling and reconstruction
of a one-dimensional deterministic signal ϕ (t) defined for t ∈ [0, 1).

nonuniform scalar quantization based on a companding model
– where the input value goes through a nonlinear mapping
(compressor), the obtained value being uniformly quantized
and then mapped back via the inverse of the nonlinearity
(expander). Moreover, under high-rate assumptions, Bennett
derived a formula for approximating the quantizer distortion based on the source probability-density-function and
the derivative of the nonlinear compressor function. This
important formula is often referred to as Bennett’s integral.
Bennett’s work was followed by a long line of theoretic and
algorithmic studies of the nonuniform quantization problem.
A prominent branch of research addressed the scenario of
nonuniform quantization at high-rates (for example, see [2],
[3]), where the quantizer has a large number of representationvalues to be wisely located at the quantizer-design stage. The
popularity of high-rate studies is due not only to their relevance
in addressing high-quality coding applications, but also to the
possibility to gain useful theoretical perspectives. Specifically,
reasonable assumptions made for high-rate quantization often
led to convenient closed-form mathematical solutions that, in
turn, provided deep insights into rate-distortion trade-offs.
Sampling has been prevalently studied for the purpose
of acquiring signals based on coarse characterizations such
as their bandwidth. The classical uniform sampling theorems (see the detailed review in [4]) were also extended to
the nonuniform settings where the variable sampling-rate is
adapted to, e.g., the signal’s local-bandwidth estimate [5]–
[8]. In [6]–[8], a nonuniform sampling grid was designed
by mapping the signal’s time axis to expand portions corresponding to high local-bandwidth at the expense of segments
containing content of lower local-bandwidth. After applying
this transformation, signal acquisition could subsequently be
carried out using uniform sampling. While these works refer to acquiring signals of unknown content, we consider
here sampling of fully-accessible deterministic signals for
representation-oriented tasks such as compression (as will be
described below). Despite the above discussed essential difference in the sampling and reconstruction procedures, some
results presented in this paper conceptually resemble sampling
methods for acquisition of unknown signals using coordinate

2

transformation. The results presented here for sampling given
deterministic signals may also be interpreted as complements
of the adaptive sampling paradigms previously proposed for
signal acquisition. In our deterministic settings, the analysis
and sampling-rate design directly rely on properties obtained
from the continuous signal instead of coarse, local averagedbased, spectral characterizations.
Practical signal compression (see, e.g., [9], [10]) usually
considers a digital high-resolution input function that is going
through a nonuniform adaptive partitioning of the domain
followed by a quantization procedure that suits the partition
size and shape. Analytically, the digital input can be regarded
as a signal defined over a continuous domain that needs to
be nonuniformly sampled and quantized. In practical coding
procedures, the nonuniform domain segmentation often relies
on structured partitions of the domain in order to reduce
the bit-cost and computational complexity (examples for the
common use of tree-structured partitioning, e.g., defined based
on binary or quad trees, are available in [10], [11]). This
application is one of the main motivations to the study
of adaptive sampling of deterministic signals. Nevertheless,
we consider here the nonuniform sampling problem in its
most general form in order to provide insights to the basic
sampling problem, that may be useful to problems beyond
compression. For example, adaptive sampling is often used
in computer graphics in the task of Halftoning (e.g., [12],
[13]) and in rendering an image from a 2D/3D-graphical
model, where the nonuniform sampling pattern is described
by point distributions (e.g., [14]–[18]). It should be noted
that the signal processing and computer graphics contexts
of the sampling problem are quite different. Interestingly,
it was suggested in [16] to practically employ the Lloyd
algorithm [19], originally intended for quantizer design, for
improving the sampling point distribution. Recent studies
[17], [18] extended the point-distribution computation to rely
on kernel functions, so that the resulting task resembles a
generalized signal-sampling procedure. While the conceptual
relation between vector-quantization and data-representation
has been understood and used in graphics applications, we
argue that a clear mathematical analysis that demonstrates the
connection between quantization and high-resolution signalsampling has not been provided yet.
Some recent papers [20]–[22] explore the sampling task
from a stochastic, information theoretic rate-distortion tradeoff
perspective, considering lossy compression of the sampled
values. Their focus is on uniform [20], [22] and nonuniform
[21] sampling where the sampling design is based on spectral characteristics of stationary signals. Another interesting
direction of nonuniform sampling was explored in [23], [24]
where the sampling intervals are recursively determined based
on previously recorded data. This process does not require the
coding of the sampling intervals. The sampling design in [23],
[24] is based on either the local properties of the stochastic
process whose realizations are sampled, or the local behavior
of deterministic signals exhibited in their Taylor expansion.
In this paper we theoretically explore the task of highresolution adaptive sampling of a given deterministic signal.
The sampling analysis provided emerges from ideas similar

to the ones that were applied to the study of high-rate
quantization [25], thereby linking sampling and quantization
in a new and enlightening way. We analytically formulate the
optimal high-resolution sampling of one-dimensional signals,
based on the Mean-Squared-Error (MSE) criterion, showing
that the optimal partitioning is determined by the cube-root
of the signal-derivative energy. This result corresponds to
the work by Panter and Dite [2], where the optimal onedimensional quantizer is designed based on the cube-root of
the probability-density-function. We also connect this result
to the fundamental analysis given by Bennett [1] for high-rate
nonuniform quantization based on companding. We continue
our analysis by addressing the problem of sampling Kdimensional signals. We show that the optimal sampling-point
K
-power of the
density is determined by the density of the K+2
signal-gradient energy, a generalization of the one-dimensional
result. We obtain this result based on assumptions that parallel
a famous conjecture given by Gersho [3] in his analysis of
high-rate quantization.
Gersho’s conjecture states that, for asymptotically high rate,
the optimal K-dimensional quantizer is formed by regions
that are approximately congruent and scaled versions of a Kdimensional convex polytope that optimally tessellates the Kdimensional space (where the polytope optimality is in the
sense of minimum normalized moment of inertia [3]). The
latter holds for a given K only if the optimal tessellation of
the K-dimensional space is a lattice, and therefore constructed
based on a single optimal polytope. This assumption significantly simplifies the explicit calculation of the quantizer’s
distortion. This conjecture draws its credibility from two
prominent sources. First, it is known that the best tessellation
is a lattice for K = 1 (based on equal-sized intervals) and for
K = 2 (based on the hexagon shape [26]). While not proven
yet for K = 3, it is also believed that the optimal threedimensional tessellation is the body-centered cubic lattice [27],
[28]. Second, Gersho’s distortion formula conforms to the
structure of the expression rigorously obtained by Zador [29].
Moreover, Gersho’s conjecture determines the value of the
multiplicative-constant (left unspecified) in Zador’s formula.
Hence, the possible inaccuracy in Gersho’s conjecture will
affect only the multiplicative constant, and the deviation is
assumed to be moderate [25]. Due to this, the conjecture,
still unproved for K ≥ 3, is widely considered as a valuable
tool for analysis of high-rate quantization (see the thorough
discussion in [25]).
The analysis provided in this paper to high-resolution multidimensional sampling is based on two main assumptions.
First, the signal is assumed to be approximately linear within
each of the sampling regions. Second, sampling regions are
assumed to be approximately congruent and scaled forms of
the optimal K-dimensional tessellating convex polytope – just
as in Gersho’s conjecture for high-rate quantization. These
high-resolution assumptions yield our main result that the
optimal sampling-point density is determined by the density of
K
-power of the signal-gradient energy. We emphasize
the K+2
the importance of the signal’s local-linearity assumption as
a prerequisite stage that mathematically connects the signalsampling problem to the conjecture on the high-resolution cell

3

arrangement.
This paper is organized as follows. In section II we mathematically analyze the optimal sampling of one-dimensional
signals. In section III we provide an experimental evaluation of
the proposed sampling method for one-dimensional signals. In
section IV we generalize our study by theoretically addressing
the optimal sampling of multidimensional signals. Section V
concludes this paper.
II. A NALYSIS FOR O NE -D IMENSIONAL S IGNALS
A. Optimal High-Resolution Sampling
Let us consider a one-dimensional signal ϕ (t) defined as a
differentiable function
ϕ : [0, 1) → [ϕL , ϕH ] ,

(1)

defined for t in the interval [0, 1) and having values from
a bounded range [ϕL , ϕH ]. This signal is sampled based
on its partition to N ∈ N non-overlapping variable-length
segments, where the ith subinterval [ai−1 , ai ) is associated
with the sample ϕi for i = 1, ..., N . We assume a segmentation
structure satisfying a0 = 0, aN = 1, and ai−1 < ai for
i = 1, ..., N . The sampling procedure is coupled with a
reconstruction that provides the continuous-time piecewiseconstant signal
ϕ̂ (t) = ϕi

for

t ∈ [ai−1 , ai ) .

(2)

The sampling is optimized to minimize the mean-squared-error
(MSE), expressed as


N −1

N



M SE {ai }i=1 , {ϕi }i=1 =

N Zai
X

where the remainder term o (|t − ti |) corresponds to our highresolution assumption in describing the approximation error
for t → ti .
Using the linear approximation, and by (4), the optimal
sample in the ith subinterval is given by
= ϕ (ti ) + o (∆i )
ϕopt
i

(6)

due to the fact that the average of a linear function over an
interval is its value at the interval’s center. The term o (∆i )
describes the error in the sample value due to the linear
approximation. The size of the error term as ∆i → 0 is
provided in Appendix A and relies on assuming that the second
derivative of the signal exists and bounded over the sampling
interval. Then, the MSE of the ith subinterval is expressed as
Zai

1
M SE i (ai−1 , ai ) =
∆i

2
ϕ (t) − ϕopt
dt
i

ai−1


1
2
(8)
(ϕ0 (ti )) ∆2i + o ∆2i
=
12
revealing the effect of signal-derivative energy on the sampling
MSE. The derivation of (8) is detailed in Appendix A.
Returning to the total sampling MSE, corresponding to optimal coefficients, and relying on its relation to the subintervals
MSE yields (recall the implicit normalization to unit-interval
length)
N

 X
N −1
M SE {ai }i=1 =
∆i M SE i (ai−1 , ai )
i=1
N

=

2

(ϕ (t) − ϕi ) dt


1 X 0
2
(ϕ (ti )) ∆3i + o ∆3max
12 i=1
(9)

i=1a

i−1

(3)
exhibiting the roles of the signal partitioning and sample
values. Note that in (3), as also in this entire section, averaging
over the unit-interval length is implicit.
N
Optimizing the sampling coefficients, {ϕi }i=1 , given a parN
titioning {ai }i=0 is a convex problem that can be analytically
solved to show that the optimal ith sample is the signal average
over the corresponding subinterval, namely,
ϕopt
i

1
=
∆i

(7)

Zai
ϕ (t) dt

(4)

ai−1

where ∆i , ai − ai−1 is the length of the ith subinterval.
We continue the analysis by assuming high-rate sampling,
meaning that N is large enough to result in small samplingintervals that, however, may still have different lengths. Furthermore, the sampling intervals are assumed to be sufficiently
small such that, within each of them, the signal is well approximated via a local linear form – an argument that is analyzed
next. Accordingly, for t ∈ [ai−1 , ai ) (i = 1, ..., N ), we consider the first-order Taylor approximation of the signal about
the center of the ith sampling interval, ti , 12 (ai−1 + ai ),
ϕ(t) = ϕ (ti ) + ϕ0 (ti ) · (t − ti ) + o (|t − ti |)

(5)

where ∆max , max{∆1 , ..., ∆N } is the largest subinterval.
Accordingly, the inaccuracies in evaluating the MSE
using the

signal linearity assumption are of size o ∆3max .
Let us connect our discussion to the classical approach
of studying high-rate quantization based on the reproductionvalue density function (see examples in [3], [19], [25], [30]).
Following our scenario of high-resolution sampling we assume
that the sampling-point layout can be described via a samplingpoint density function, λ (t), such that a small interval of
¯ around t̄ approximately contains ∆
¯ · λ (t̄) sampling
length ∆
points. Moreover, the sampling-point density is related to the
sampling intervals via
λ (t) ≈

1
,
N · ∆i

for t ∈ [ai−1 , ai ) .

(10)

As we consider arbitrarily large values of N , the density λ (t)
is assumed to be a smooth function.
Then, plugging the relation ∆i ≈ 1/ (N · λ (ti )) into the
sampling-error expression (9), in addition to approximating
the sum as an integral and omitting the explicit inaccuracy
term, yields
1
M SE (λ) ≈
12N 2

Z1
0

2

(ϕ0 (t))
dt.
λ2 (t)

(11)

4

Here the sampling structure and the resulting MSE are determined by the sampling-point density λ(t). The last MSE
expression can be interpreted as the sampling equivalent of
Bennett’s integral for nonuniform quantization [1]. Commonly,
the expression form in (11) is minimized via Hölder’s inequality (see examples in [3], [25], [30]). For our problem, see
details in Appendix B, we have the optimal sampling-point
density given by
q
2
3
(ϕ0 (t))
opt
,
(12)
λ (t) = 1 q
R 3
2
(ϕ0 (z)) dz
0

and the optimal sampling MSE is
 1
3
Z q

1
3
2

M SE λopt ≈
(ϕ0 (t)) dt .
12N 2

(13)

0

The error expression in (13) is a product composed of two
parts: sampling error for a simple linear signal (with a unit
slope), and a term expressing the nonlinearity of the given
signal based on its derivative energy. Moreover, Hölder’s
inequality also shows that (13) is the global minimum. The
MSE expression in (13) is the sampling counterpart of the famous Panter-Dite formula for high-rate quantization MSE [2].
Evidently, while the quantization derivations were determined
by the probability-density-function of the source, the sampling
analysis provided here depends on the signal-derivative energy.
Using the sampling-point density, λ(t), we can implement
our nonuniform sampling via the companding design. Companding [1] is a widely-known technique for implementing
nonuniform quantization based on a uniform quantizer. This
is achieved by applying a nonlinear compressor function on the
input value, then applying uniform quantization and mapping
the result back via an expander function (the inverse of the
compressor). Since in sampling we address the problem of
discretizing the signal domain, the corresponding compressor
and expander functions operate on the signal domain (i.e., as
a nonlinear scaling of the time axis). The optimal compressor
function is defined based on the optimal density (12) as
Rt q
2
3
(ϕ0 (z)) dz
Zt
0
u (t) = λopt (z)dz = 1 q
.
(14)
R 3
2
(ϕ0 (z)) dz
0
0

The corresponding expander, v (τ ), is the inverse function of
the compressor, hence, it can be defined via the relation
v(τ
Z )

λopt (z)dz = τ.

(15)

1) Using a Strictly-Positive Extension of the Signal’s
Derivative-Energy:
Consider the expander function in (15) to construct the boundaries of the nonuniform segmentation of [0, 1) via
 
i
opt
, i = 1, ..., N,
(16)
ai = v
N
i.e., evaluating the inverse of the compressor function at N
equally-spaced points. However, for non-monotonic signals
(that also may have constant-valued segments), it is likely to
need expander values at points where the compressor function
is not invertible – an issue that can be solved as follows. Since
the problem occurs where the signal derivative is zero, we
define the following extension of the derivative energy:
(
2
2
(ϕ0 (t)) for (ϕ0 (t)) > ε
2
(17)
gε (t) =
ε
otherwise
where ε > 0 is an arbitrarily small constant. The corresponding extension of the optimal sampling-point density (12) is
p
3
gε2 (t)
opt
λε (t) = 1
(18)
R p
3
2
gε (z)dz
0

λopt
ε (t)

enables treatment of nonAccordingly, the density
monotonic signals, while closely approximating the density
λopt (t) in (12). Replacing λopt (t) with λopt
ε (t) in (14)-(15)
provides a practically useful compressor-expander pair, in the
sense that the compressor function is invertible everywhere in
the domain, assuring the computations in (16).
2) Sequential Solution via Integration Thresholds:
Eq. (15) defines the segment level aopt
(i = 1, ..., N ) as
i
opt
a
Zi

λopt (z)dz =

i
N

(19)

0

implying that aopt
can be evaluated given the former partii
tioning level aopt
via
i−1
opt
a
Zi

q
3

(ϕ0

1
(t)) dt =
N
2

aopt
i−1

Z1 q
3

2

(ϕ0 (z)) dz

(20)

0

q
2
3
i.e., we just need to continuously integrate (ϕ0 (t)) starting
opt
at t = ai−1 until we reach the threshold level defined by the
right side of the last equation. Then, the time of achieving the
threshold level defines aopt
i . Since a partitioning level is placed
at the (first) time the threshold level is obtained, intervals
of zero signal derivative-energy do not cause ambiguity in
segmentation definitions. Importantly, one should note that the
(optimal) threshold, defined according to Eq. (20) as

0

The last equation has a unique solution for strictly monotonic
signals where the signal-derivative energy is positive over
the entire domain, leading to a strictly-monotonic increasing
compressor function (14) that is invertible and, thus, defines
the expander function. Two suggestions for the treatment of
non-monotonic signals will be described next.

Topt

1
=
N

Z1 q
3
2
(ϕ0 (z)) dz,

(21)

0

requires knowing the signal derivative over the entire [0, 1)
interval. Note however that the segmentation levels are here
well-defined for any given nonzero signal.

5

B. The Sampling-Quantization Duality at High Resolution
The above developments for optimal high-resolution sampling show a clear correspondence to well-known results
from high-rate quantization studies. The connection between
sampling and quantization, in the high-resolution settings,
emerges from local linear approximation of the signal within
each of the sampling intervals, implying a constant signalderivative within each interval. This mirrors the assumption
of locally constant probability density in the high-rate quantization problem. In turn, our construction yields an optimal
sampling procedure determined by the signal’s derivativeenergy density. Specifically, the MSE forms in (11) and (13)
and the corresponding companding approach are as in the
classical high-rate quantization results, however, here they
correspond to the signal’s derivative-energy density instead of
the probability density function of a source to be quantized.
In the sequel, let us consider the scalar quantizer design
for a random variable X, defined by the probability density
function pX (x) that may be positive only for x ∈ [xL , xH ).
Then, recall Bennett’s integral [1] for the MSE of high-rate
quantization using N reproduction values:
1
M SE Q (λQ ; pX ) ≈
12N 2

ZxH
xL

pX (x)
dx
λ2Q (x)

(22)

where, as interpreted in [19], λQ (x) is the reproduction-point
density that defines the quantizer structure to optimize, and
pX (x) is the given input probability density function.
In this subsection we distinguish between optimizations
originating in sampling and quantization procedures by the
following notations: the quantizer design procedure addresses
Bennett’s integral (22) to minimize the quantization distortion,
M SEQ , as a function of the reproduction-point density λQ .
Similarly, the sampling procedure considers the sampling
counterpart of Bennett’s integral presented in (11), that using
this subsection’s notations is formulated as
1
M SE S (λS ; ϕ) ≈
12N 2

Z1

2

(ϕ0 (t))
dt,
λ2S (t)

0

(23)

where M SES , the distortion of sampling the signal ϕ(t), is
to be minimized by optimizing the sampling-point density λS .
1) Optimal Sampling via Optimal Quantizer Design:
For a given differentiable signal, ϕ (t) defined for t ∈ [0, 1),
we can define a probability density function, pϕ (x), defined
for x ∈ [0, 1) via
0

pϕ (x) =

(ϕ (x))
R1

(ϕ0

2

,

(24)

2

(z)) dz

of pϕ (x) from (24) into (22), resulting in
1
M SE Q (λQ ; pϕ ) ≈
12N 2 Eϕ0

Z1
0

2

(ϕ0 (x))
dx.
λ2Q (x)

where λQ (·) is the reproduction-point density defining the
high-rate quantizer structure, and the total derivative-energy
of the signal is denoted here as
Z1
2
Eϕ0 , (ϕ0 (z)) dz.
(26)
0

Then, the high-resolution sampling MSE, M SES , formulated
in (23) is related to the high-rate quantization MSE in (25)
via
1
M SE S (λQ ; ϕ) ,
(27)
M SE Q (λQ ; pϕ ) ≈
Eϕ0
i.e., the MSE of quantizing a source defined by pϕ using
a quantizer structured according to the partitioning induced
from λQ is equivalent, up to a normalization in the total
derivative-energy of the signal, to the MSE of sampling the
signal ϕ (t) according to the segmentation defined by λQ .
The MSE relation (27) shows that the optimal λQ for the
above quantization problem, is also optimal for sampling of
the signal ϕ (t). Indeed, Bennett’s integral for quantization of
pϕ (x) is minimized for λQ formulated exactly as the optimal
sampling-point density in (12).
2) Optimal Quantizer Design via Optimal Sampling:
Given a random variable X corresponding to the probability
density function pX (x) defined for x ∈ [xL , xH ), we can
construct a signal ϕX (t) defined for t ∈ [0, 1) via
1
ϕX (t) =
xH − xL

xL +t·(x
Z H −xL )

p

pX (x)dx.

(28)

xL

We take the positive square-root of pX (x), hence, the signal
ϕX (t) is monotonically non-decreasing. The derivative of
ϕX (t) (with respect to t) is
q

ϕ0X (t) = pX xL + t · (xH − xL ) .
(29)
Then, according to (23), the sampling MSE of ϕX (t) for a
sampling-point density λS is
Z1 0
2
1
(ϕX (t))
M SE S (λS ; ϕX ) ≈
dt
(30)
12N 2
λ2S (t)
0

that here equals to
1
M SE S (λS ; ϕX ) ≈
12N 2

Z1
0


pX xL + t · (xH − xL )
dt.
λ2S (t)
(31)

0

i.e., the probability density is the signal’s local squaredderivative normalized by the total derivative energy. Evidently,
pϕ (x) is non-negative valued and integrates to one, hence, it
is a valid probability density function.
Bennett’s integral for high-rate quantization of a source distributed according to pϕ is obtained by plugging the definition

(25)

By changing the integration variable to x = xL +t·(xH − xL )
we get
ZxH
1
pX (x)
M SE S (λS ; ϕX ) ≈
dx.
L
12N 2 (xH − xL )
λ2S ( xx−x
)
H −xL
xL

(32)

6

By referring to (22), the last form shows the following relation
between the sampling MSE and the quantization MSE
M SE S (λS ; ϕX ) ≈

1
M SE Q (λQ ; pX ) (33)
(xH − xL )

where, for x ∈ [xL , xH ),

λQ (x) = λS

x − xL
xH − xL


.

(34)

The last result means that for a random variable X, with a
given probability density function pX (x), one can design the
optimal high-rate quantizer by considering the high-resolution
sampling problem of the signal ϕX (t) defined in (28). Equations (33)-(34) show that the MSEs of the two procedures
are equal up to a normalization by the width of the valuerange of X, and a linear transformation of the coordinates of
the partitioning structure (see Eq. (34)). Hence, implementing
the optimal sampling-point density (12) for the signal ϕX (t),
together with the appropriate linear-transformation of the
coordinates, provides the structure that minimizes the highrate quantization MSE for pX (x).
C. Numerical Demonstrations
We now turn to study our theoretical results by applying
them for analytic signals.
1) Exponential Signals: We consider an exponential signal
of the form
ϕ(t) = eαt

, t ∈ [0, 1)

(35)

where α > 0 is a real-valued parameter determining the growing rate (see examples in Fig. 2a). The signal-derivative energy
is expressed as (ϕ0 (t))2 = α2 e2αt , which is positive valued
for t ∈ [0, 1) (Fig. 2b). Therefore, the optimal sampling-point
density for this monotonic signal is expressed as
2

λopt (t) =

e 3 αt
2α
· 2α
3 e3 − 1

(36)

Plugging (36) into (15) defines the optimal expander, v (τ ),
via

We evaluate the gain of our approach with respect to the uniform sampling. (for i = 1, ..., N ). The MSE of uniform highorm
resolution sampling is calculated by setting aunif
= Ni in
i
the MSE expression in (9), yielding
n
oN −1 

α
unif orm
e2α − 1 . (41)
=
M SE
ai
2
24N
i=1
The MSE of the two sampling methods are compared in
Fig. 2c and Fig. 2d for various exponential signals (via the
α parameter) and sampling resolutions (the parameter N ),
respectively. Note that the MSE values are normalized by the

R1 2
1
e2α − 1 .
signal energy, here expressed as
ϕ (t) = 2α
t=0

Figure 3 presents an example for the compressor function
and the corresponding nonuniform sampling of an exponential
signal.
2) Cosine Signal: Let us demonstrate our approach for a
non-monotonic signal of the form:
ϕ(t) = cos (2παt)

, t ∈ [0, 1)

(42)

where α > 0 is an integer determining the number of periods
contained in the [0, 1) interval (see example for α = 5 in Fig.
4a). The derivative energy of the cosine signal in (42) is
(ϕ0 (t))2 = 4π 2 α2 sin2 (2παt)

, t ∈ [0, 1).

(43)

As demonstrated in Fig. 4b, the signal-derivative energy is
j
for j = 0, 1, ..., 2α − 1.
zero only at the points t = 2α
The compressor function does not lend itself here to a simple
analytic form, nevertheless, it can be constructed numerically
via its definition as a cumulative-density-function (see Eq.
(14)) providing the compressor curve in Fig. 4c. The 2α points
of zero signal-derivative energy are not a real obstacle in
our numerical construction and, anyway, their corresponding
values can be replaced by an arbitrarily small ε value1 as suggested above. The expander function is numerically formed as
the inverse of the compressor curve. The resulting nonuniform
sampling structure (Fig. 4d) shows its adaptation to the local
signal derivative and to the periodic nature of the signal.
3) Chirp Signal: The cosine signal in (42) can be extended
to the following chirp signal with a linearly increasing frequency:

2

e 3 αv(τ ) − 1
2

e3α − 1

=τ

(37)

that yields
v (τ ) =

 2


3
log e 3 α − 1 · τ + 1 ,
2α

(38)

and the corresponding optimal nonuniform partitioning is
determined via (16) as


 i
2
3
3α − 1
log
e
·
+
1
, i = 1, ..., N.(39)
aopt
=
i
2α
N
The MSE corresponding to the optimal nonuniform sampling
is calculated using (13) and expressed as


 2
3
9
N −1
M SE aopt
=
e3α − 1 .
(40)
i
2
i=1
32αN

ϕ(t) = cos (2πt (1 + αt))

, t ∈ [0, 1).

(44)

Here the α > 0 parameter determines the linear growth-rate
of the frequency (Fig. 5a exemplifies this for α = 5). The
signal-derivative energy of the chirp (44) is
2

(ϕ0 (t))2 = 4π 2 (1 + 2αt) sin2 (2πt (1 + αt))

, t ∈ [0, 1).
(45)

The nonuniform sampling of the chirp is demonstrated in
Fig. 5d. Comparison between the nonuniform sampling of the
cosine signal (Fig. 4) and the chirp signal (Fig. 5) reveals the
influence of the varying frequency embodied in the chirp.
1 In the numerical demonstrations provided here we set the value of ε to
the smallest positive floating-point value available in Matlab via the command
eps, which returns the value 2−52 .

7

(a)

(b)

(c)

(d)

Fig. 2. Demonstration for an exponential signal ϕ(t) = exp(αt) for α > 0. (a) The exponential signal for several α values. (b) The signal-derivative energy
for several α values. In (c)-(d), theoretical reconstruction-MSE obtained via nonuniform and uniform sampling procedures are compared. (c) evaluated for
various α values and N = 50. (d) evaluated for a range of sampling resolutions (N ) and α = 3.

(46)
and the corresponding


25
(ϕ0 (t))2 = 1


0

signal-derivative energy is
for t ∈ [0, α) ∪ [1 − α, 1)
for t ∈ [α, 2α) ∪ [1 − 2α, 1 − α)
for t ∈ [2α, 1 − 2α) .
(47)

(a) Compressor Curve

Figure 6 demonstrates the sampling procedure for a piecewiselinear signal constructed for α = 81 . Fig. 6d shows that using
the extended sampling-point density function from (18), where
small derivative-energy values are clipped to ε, provides a
good solution to the constant-valued region of this signal.
III. E XPERIMENTAL R ESULTS : S AMPLING OF
O NE -D IMENSIONAL S IGNALS
In this section we present experimental evaluations of the
procedure proposed in section II for nonuniform sampling of
one-dimensional signals.
A. The Main Competing Method: Tree-Structured Nonuniform
Sampling

(b) Nonuniform Sampling
Fig. 3. Optimal nonuniform sampling (N = 50) of an exponential signal
ϕ(t) = exp(3t). (a) the mapping between uniform to nonuniform samplingspacing. (b) shows the original signal (magenta), the reconstructed signal from
nonuniform sampling (blue), and the partitioning to sampling intervals (red).

4) Piecewise-Linear Signal: Now we turn to evaluate our
sampling method on a piecewise-linear signal containing a
significant portion of a constant value. This will show the
suitability of the proposed method also to signals with a
zero derivative-energy over a large continuous sub-interval.
We define here a specific signal form as


5t
for t ∈ [0, α)





(t
−
α)
+
5α
for t ∈ [α, 2α)

ϕ(t) = 6α
for t ∈ [2α, 1 − 2α)



−
(t
−
(1
−
2α))
+
6α
for t ∈ [1 − 2α, 1 − α)



−5 (t − (1 − α)) + 5α for t ∈ [1 − α, 1)

Let us consider sampling based on a nonuniform partitioning that is represented using a binary tree. The approach
examined here is inspired by the general framework given in
[31] for optimizing tree-structures for various tasks, and is also
influenced by the discrete Lagrangian optimization approach
[32] and its application in coding [33], [34].
The suggested approach relies on an initial tree, which is a
full d-depth binary-tree, representing a uniform segmentation
of the interval [0, 1) into 2d sampling intervals of 2−d length
(see example in Fig. 7a). The segmentation of the interval
[0, 1) is described by the leaves of the binary tree: the interval
location and length are defined by the leaf position in the tree,
specifically, the tree-level where the leaf resides in determines
the interval length. The examined nonuniform segmentations
are induced by all the trees obtained by repeatedly pruning
neighboring-leaves having the same parent node (examples are
given in Figures 7b-7c). The initial d-depth full-tree together
with all its pruned subtrees form the set of relevant trees,
denoted here as Td .
The leaves of a tree T ∈ Td form a set denoted as
L(T ), where the number of leaves is referred to as |L(T )|.

8

(a)

(b)

(c)

(d)

Fig. 4. Demonstration for a cosine signal ϕ(t) = cos(10πt). (a) The signal. (b) The signal-derivative energy. (c) The optimal compressor curve. (d) Optimal
nonuniform sampling using N = 100 samples (the partitioning of the signal domain is in red).

(a)

(b)

(c)

(d)

Fig. 5. Demonstration for a chirp signal ϕ(t) = cos (2πt (1 + 5t)). (a) The signal. (b) The signal-derivative energy. (c) The optimal compressor curve. (d)
Optimal nonuniform sampling using N = 100 samples (the partitioning of the signal domain is in red).

(a)

(b)

(c)

(d)

1
.
8

Fig. 6. Demonstration for a piecewise-linear signal as in Eq. (46) with α =
(a) The signal. (b) The signal-derivative energy. (c) The optimal compressor
curve. (d) Optimal nonuniform sampling via the extended sampling-point density using N = 20 samples (the partitioning of the signal domain is in red).

Accordingly, the tree T represents a (possibly) nonuniform
partitioning of the [0, 1) interval into |L(T )| sampling intervals. A leaf l ∈ L(T ) resides in
of the
h the h(l) level

t right
tree and corresponds to the interval alef
,
a
of
length
(l)
(l)
∆ (l) = 2−h(l) . Following the analysis in section II, the optimal sample corresponding to the leaf l ∈ L(T ) is expressed
via (4) as
aright
(l)

ϕopt
(l) =

Z

1
∆ (l)

ϕ (t) dt.

(48)

t
alef
(l)

Consequently, as the tree leaves correspond to a segmentation
of the [0, 1) interval, the sampling MSE induced by the tree
T ∈ Td is calculated based on the leaves, L(T ), via
aright
(l)

M SE (T ) =

X
l∈L(T )

Z
t
alef
(l)



ϕ (t) −

ϕopt
(l)

2

dt.

(49)

For a signal ϕ(t) and a budget of N samples, one can
formulate the optimization of a tree-structured nonuniform
sampling as
minimize
T ∈Td

M SE (T )
(50)

subject to |L(T )| = N,
i.e., the optimization searches for the tree with N leaves
that provides minimal sampling MSE. The unconstrained
Lagrangian form of (50) is
min {M SE (T ) + µ|L(T )|} ,

T ∈Td

(51)

where µ ≥ 0 is a Lagrange multiplier that reflects the
constraint |L(T )| = N . However, it should be noted that due to
the discrete nature of the problem such µ does not necessarily
exist for any N value (for details see, e.g., [31], [32]). The

9

problem (51) can also be written as


aright


(l)


Z

 X
2

opt
min
ϕ (t) − ϕ(l) dt + µ|L(T )| . (52)
T ∈Td 




l∈L(T ) alef t
(l)

Note that, due the non-intersecting sampling intervals, the
contribution of a leaf, l ∈ L(T ), to the Lagrangian cost is
aright
(l)

Z
C (l) =



ϕ (t) − ϕopt
(l)

2

(a)

dt + µ,

(53)

t
alef
(l)

evaluated for the corresponding sampling interval.
The discrete optimization problem (52) of finding the optimal tree for a given signal and a Lagrange multiplier µ is
addressed by the following procedure. Start from the full ddepth tree and determine the corresponding sampling intervals
and their optimal samples, squared errors, and contributions
to the Lagrangian cost (53). Go through the tree levels from
bottom and up, in each tree level find the pairs of neighboring
leaves having the same parent node and evaluate the pruning
condition: if
C (left child) + C (right child) > C (parent)

(54)

is true, then prune the two leaves – implying that two sampling
intervals are merged to form a single interval of double length
(thus, the total samples in the partitioning is reduced by one).
If the condition (54) is false, then the two leaves (and the
corresponding sampling intervals) are kept. This procedure is
continued until reaching a level where no pruning is done, or
when getting to the tree root.
B. Evaluation and Comparison
The proposed sampling method is compared to two
other sampling approaches. The first is the trivial, however
commonly-used, uniform sampling, where the signal domain is
partitioned into equal-size sampling-intervals. Specifically, for
a budget of N samples the signal domain [0, 1) is segmented
according to ai = Ni for i = 0, 1, ..., N . The samples are
determined as the averages of the corresponding sampling
intervals. The second competing method is a nonuniform
sampling based on a binary-tree structure that is adapted to
the signal via a Lagrange optimization, as presented in the
previous subsection.
We examined sampling of several signals defined analytically. A grid of µ values was set for the tree-structured
Lagrange optimization, this determined the number of samples to consider2 . First, the sampling of the cosine signal
ϕ (t) = 255 · cos (10πt) was examined, and showed that our
method consistently outperforms the uniform and the treestructured techniques for various amounts of samples (Fig.
8a). These observations were further established by examining
2 Note that the proposed method and the uniform sampling do not rely
on µ and operate directly based on a given number of samples, however,
we defined the examined sample budgets based on the µ grid of the treestructured sampling in order to maintain an accurate comparison between all
the sampling methods.

(b)

(c)

Fig. 7. Segmentations of the [0, 1) interval produced by binary trees. The
leaves, which determine the partitioning, are colored in green. (a) A full
binary tree of depth 3 and the corresponding uniform segmentation into 8
sub-intervals. (b) A tree obtained by a single pruning of the full tree, thus the
partitioning includes 7 sub-intervals. (c) A tree obtained by several prunings
resulting in 4 leaves/segments.

sampling of the chirp signal, ϕ (t) = 255 · cos (2πt (1 + 5t)),
and the piecewise-linear signal from Eq. (46) that was scaled
by a factor of 255 (Figures 8b-8c).
IV. A NALYSIS FOR M ULTIDIMENSIONAL S IGNALS
In this section we discuss the theoretic analysis of optimal high-resolution sampling by studying the problem for
sampling multidimensional signals. The analysis for the onedimensional case is relatively simple since the partition is
characterized by sampling-interval lengths (see section II).
However, when considering the multidimensional problem
the required analysis becomes more intricate as the sampling regions are, in general, of arbitrary shape and size.
Accordingly, our analysis in this section is conceptually and
mathematically similar to the study of multidimensional highrate quantization provided by Gersho [3], which generalized
the one-dimensional theory of Bennett [1] and the distortion
formula of Panter and Dite [2]. Thus, we provide a theoretic
framework for optimal multidimensional signal sampling at
high resolution.
The given differentiable signal, ϕ (x), is defined over a
K-dimensional unit cube, C K , [0, 1]K , and is scalar realvalued from a bounded range, i.e., ϕ (x) ∈ [ϕL , ϕH ] for any
x ∈ C K . The signal goes through a sampling procedure in
order to provide a discrete representation using N ∈ N (scalar
N
valued) samples, {ϕi }i=1 , corresponding to a partitioning of
N
C K to N distinct multidimensional regions, {Ai }i=1 , such that
K
∪N
i=1 Ai = C . Again, we consider a reconstruction procedure
providing the continuous-domain piecewise-constant signal
ϕ̂ (x) = ϕi

for

x ∈ Ai ,

(55)

10

(a) Cosine Signal

(b) Chirp Signal

(c) Piecewise-Linear Signal

Fig. 8. Performance comparison of the proposed sampling method, the uniform sampling approach, and optimized tree-structured sampling. The curves
present the sampling MSE obtained for various sample budgets.

and optimization in the sense of minimizing the overall MSE,
here formulated as (note the implicit normalization in the unitcube volume)
N Z
 X

2
N
N
(ϕ (x) − ϕi ) dx.
M SE {Ai }i=1 , {ϕi }i=1 =
i=1 A

i

(56)
Moreover, as before, the optimal sampling coefficients given
N
the partitioning {Ai }i=1 can be analytically determined, showing that
Z
1
ϕopt
=
ϕ (x) dx,
i = 1, ..., N
(57)
i
V (Ai )
Ai

where V (Ai ) is the volume of the region Ai .
We consider the case of optimal high-resolution sampling
(i.e., N is very large), and assume that the optimal sampling
regions are all small enough and appropriately shaped such
that we can further presume that the signal ϕ (x) is well approximated within each region by a linear form that is locally
determined. The last assumption emerges as a generalization
of the assumptions for the one-dimensional case that were
presented and analyzed in detail in Section II-A. Specifically,
here
ϕ(x) ≈ βiT x + γi

for

x ∈ Ai ,

(58)

where βi is a K-dimensional column vector and γi is a scalar
value. Moreover, the above linear form can be determined in
the ith region by the first-order Taylor approximation
of the
R
1
signal around the region center, xi , V (A
xdx,
namely,
i)
Ai

ϕ(x) ≈ ϕ (xi ) + ∇ϕ(xi ) (x − xi )

(59)

where ∇ϕ(xi ) is the signal gradient, evaluated at the region
center, here having the form of a K-length row vector consisting of the K partial derivatives, i.e.,


∂ϕ (x)
∂ϕ (x)
∂ϕ (x)
∇ϕ(xi ) =
,
,...,
,
∂x1 x=xi ∂x2 x=xi
∂xK x=xi
(60)
where
j

th

∂ϕ(x)
∂xj

central point xi . Accordingly, the linear-form parameters of
the ith region are set to
βiT = ∇ϕ(xi )

(61)

γi = ϕ (xi ) − ∇ϕ(xi ) xi .

(62)

The local linear approximation within each region (58)
yields an optimal sampling coefficient that is approximately
the signal value at the region center, i.e.,
Z

1
βiT x + γi dx = βiT xi + γi ≈ ϕ (xi ) .
ϕopt
≈
i
V (Ai )
Ai

(63)
th

Then, the sampling MSE of the i region is
Z
2
1
M SE i (Ai ) =
ϕ (x) − ϕopt
dx
(64)
i
V (Ai )
Ai
Z
 T
2
1
βi x + γi − βiT xi + γi dx
≈
V (Ai )
Ai
Z
1
2
2
=
kβi k2 kx − xi k2 dx.
V (Ai )
Ai

Inspired by the analysis given by Gersho [3] to high-rate
quantization, we turn to interpret the last error expression using
the normalized moment of inertia of the region Ai around its
center xi , defined as
R
2
kx − xi k2 dx
M (Ai ) ,

Ai
2

K · V (Ai )1+ K

,

(65)

where this quantity is invariant to proportional scaling of the
region. Now the region MSE becomes
2

2

M SE i (Ai ) ≈ K kβi k2 M (Ai )V (Ai ) K ,

(66)

and the total MSE is expressed as
N

 X
N
M SE {Ai }i=1 =
V (Ai ) · M SE i (Ai )
i=1

≈K

N
X

2

2

kβi k2 M (Ai )V (Ai )1+ K .

(67)

i=1
x=xi

is the partial derivative of the signal in the

standard direction (j = 1, ..., K) measured at the region

We now assume that in the high-resolution scenario there
is a sampling-point density function, λ (x), such that for any

11

small volume A that contains x, the fraction of sampling points
contained in it is approximately λ (x) V (A). Furthermore, the
density function satisfies the approximate relation
λ (x) ≈

1
,
N · V (Ai )

for x ∈ Ai .

(68)

approximately congruent and scaled versions of the optimal
polytope used for sampling a linear signal. Then, Eq. (72) is
satisfied also in our case of sampling a nonlinear signal, since
scaling does not change the normalized moment of inertia. The
latter hypothesis mirrors the conjecture made by Gersho in [3]
for high-rate quantization. Consequently, the MSE expression
(70) is reformed to
Z
K · M (A∗K )
1
M SE (λ) ≈
β 2 (x)
2
2 dx.
K
K
N
λ
(x)
K

Hence, the above assumption implies that adjacent sampling
regions have similar density values.
Na and Neuhoff [30] introduced (in the context of vector
quantization) the important notion of the inertial profile, denoted here as m (x). The function m (x) is assumed to be
smooth and to approximate the normalized moment of inertia
of the cells (around their mass centers) in the neighborhood
of x. The smoothness of the inertial profile is based on the
assumption that neighboring regions have similar values of
normalized moment of inertia in a high-resolution segmentation.
The definitions of the sampling-point density function in
(68) and the inertial profile let us to express the MSE in (67)
as
N


K X
2 m (xi )
N
kβi k2
M SE {Ai }i=1 ≈
2
2 V (Ai ).
K
N i=1
λ (xi ) K
(69)

We optimize the sampling procedure by characterizing the
best sampling-point density, λopt (x), that minimizes the MSE
as expressed in (73). Similar to the optimization in the onedimensional case (see Appendix B), we rely on Hölder’s
inequality that provides us a lower bound to the MSE in (73)
in the form of

Furthermore, due to the high-resolution assumption we approximate the previous sum by the following integral,
Z
K
m (x)
M SE (λ, m) ≈
β 2 (x)
2 dx.
2
NK
λ (x) K
K

density is

x∈C

(70)

x∈C

(73)

M SE (λ) ≥

(74)
1+ K2



K · M (A∗K ) 
N

Z

 K
β 2 (x) K+2 dx

2
K

x∈C K

Here, following the application of Hölder’s inequality, the
MSE lower bound is attained when λ(x) is proportional
to β 2 (x) 1 2 implying that the optimal sampling point
(λ(x)) K

λ

opt

(x) =

β 2 (x)
R

K
 K+2
K

.

(75)

(β 2 (z)) K+2 dz

z∈C K

where
β 2 (x) , ∇ϕ(x)

2
2

.

(71)

The error expression in (70) can be interpreted as Bennett’s
integral for multidimensional sampling. This formula shows
that the sampling MSE for a given signal, which is represented
by its gradient-energy density, is determined by the sampling
point density and the inertial profile of the sampling structure.
We now argue that optimal high-resolution sampling of a
multidimensional linear signal with a uniform gradient-energy
density evaluated as 1 everywhere in C K is obtained by
partitioning the signal domain, C K , based on a tessellation
generated by a single optimal (in the sense of minimum
normalized moment of inertia) convex polytope, A∗K , such
that all the regions in the segmentation are congruent to
it3 . Accordingly, this optimal division results in a constant
normalized moment of inertia to all of the sampling regions,
namely,
M (Ai ) = M (A∗K )

for

i = 1, ..., N,

(72)

or in inertial profile terms: m (x) = M (A∗K ), i.e., a constant
function. We further assume that optimal high-resolution sampling of a nonlinear signal is obtained by regions that are
3 Based on the high-resolution assumption, we neglect cells that are intersected by the boundary of the signal domain, C K , and thus may not be
congruent to the optimal tessellating polytope.

Note that we used the fact that integrating a density function
over the entire domain should be 1.
The optimal sampling-point density demonstrates that in
regions where the derivative energy is higher, the sampling
should be denser by reducing the volumes of the relevant sampling regions. Moreover, returning to the discrete formulation
for high-resolution sampling MSE in (67) and by utilizing (68)
together with (75) and (72) shows that in the optimal solution
all the sampling regions contribute the same amount of MSE.
In addition, the results in this section are generalization of
those obtained in the analysis of one-dimensional signals in the
previous section, this can be observed by setting K = 1 and
1
M (A∗1 ) = 12
, which is the normalized moment of inertia for
one-dimensional intervals (or any other K-dimensional cube)
around their center.
The high-resolution analysis provides a theoretic evaluation
of a sampler based on its sampling-point density function.
As we described above, the suggested framework lets to
determine the optimal sampling procedure in terms of the
best sampling-point density function. In the one-dimensional
case the sampling-point density can be directly translated to
a practical sampling procedure via the companding model
(see sections II and III-B). However, in the multidimensional
case, in general, there are no direct ways to implement a
sampler based on a given sampling-point density function.
This conclusion is based on the following results from the

12

quantization field [3], [35], [36]. In [35], [36] it was shown
that optimal companding requires a compressor function that
is a conformal mapping. This result implies that a direct
implementation of multidimensional sampler based on a given
point density is limited to a minor class of signals with a
suitable gradient-energy density – thus, in general, optimal
multidimensional companding is impractical. This result was
followed by a treatment of multidimensional companding
for the vector quantization problem in limited settings that
consider suboptimal solutions and/or particular source distributions [37]–[40]. Consequently, as in the high-rate quantization
literature, the analysis provided here for the multidimensional
case is a theoretic framework for studying sampling of multidimensional signals. Specifically, it describes the optimal
sampler and allows to assess its theoretic performance. This
together with Bennett’s integral for multidimensional sampling
(Eq. (70)) can be used to evaluate the performance of practical
suboptimal sampling procedures (similar to the analysis of
practical vector quantizers in [30]).
V. C ONCLUSION
We analyze the topic of nonuniform sampling of deterministic signals as a mirror-image of nonuniform quantization. With
the advent of new technologies, adaptive sampling becomes
a viable alternative to be considered in data compression
applications. In all the above developments, the crucial localdensity controlling parameter turns out to be the local energy
of the signal gradient. A new adaptive sampling method
for one-dimensional signals is proposed and experimentally
established as a leading nonuniform-sampling approach.
A PPENDIX A
A NALYSIS OF I NACCURACIES D UE TO THE S IGNAL
L INEARITY A SSUMPTION

where the last inequality conforms with Ri (t) = o(|t − ti |)
for t → ti , as in Eq. (5).
As explained in Section II, the optimal sample value is
the signal average over the sampling interval (see Eq. (4)).
Therefore, averaging the signal in its form from (76) gives
ϕopt
i

Hence, the remainder average is the amount of inaccuracy in
the optimal sample value due to the signal linearity assumption. We analyze this quantity by bounding it
Zai

1
∆i

|Ri (t)| dt
ai−1
Zai

≤

M
2∆i

(80)

2

(t − ti ) dt =

M 2
∆
24 i

ai−1

where we used the remainder bound from (78). Using the
bound in (80) we can state that the inaccuracy in the sample
value is
Zai
1
Ri (t)dt = o(∆i )
(81)
∆i
ai−1

for ∆i → 0, as was presented in (6).
Now we proceed to analyzing the sampling MSE in the ith
interval. The basic expression given in (7) is equivalent to
1
M SEi =
∆i

Zai

ϕ2 (t)dt − ϕopt
i

2

(82)

Then, using the expressions from (76) and (79) the interval
MSE becomes
1
2
M SEi =
(ϕ0 (ti )) ∆i 2 + Di
(83)
12
where
Zai
1
Di =
Ri2 (t)dt
(84)
∆i
ai−1

1
+2ϕ (ti )
∆i
0


1
−
∆i

Considering the remainder in its integral form will be useful
for the analysis in this appendix. The absolute value of the
remainder is bounded for t ∈ [ai−1 , ai ) as follows
M
2
(t − ti )
2

2

Zai

Ri (t)dt

evaluates the deviation from the MSE obtained for the linearapproximation signal. Let us bound the MSE deviation term,
Di , as follows:
|Di | ≤
≤

(78)

(t − ti ) Ri (t)dt

ai−1

(77)

ti

Zai
ai−1

where the remainder (for the ith sampling interval) is

|Ri (t)| ≤

Zai

1
Ri (t)dt ≤
∆i

ai−1

(76)

(t − z) ϕ00 (z) dz

(79)

ai−1

ϕ(t) = ϕ (ti ) + ϕ0 (ti ) · (t − ti ) + Ri (t),

Ri (t) =

Ri (t)dt.

ai−1

Our main high-resolution assumption suggests to consider
the signal via its local linear approximation. Let us assume that
the second derivative of the signal, ϕ00 (t), is continuous and
bounded, i.e., |ϕ00 (t)| ≤ M for some positive constant M . Expressing the remainder of the first-order Taylor approximation
using its integral form, lets us to rewrite (5), for t ∈ [ai−1 , ai )
(recall that ti is the interval center), as

Zt

Zai

1
= ϕ (ti ) +
∆i

1
∆i

(85)
Zai

ai−1

Ri2 (t)dt +

2
|ϕ0 (ti )|
∆i

Zai
(t − ti ) Ri (t)dt .
ai−1

13

Then, using the bound (78) we get
Zai

1
∆i

Ri2 (t)dt ≤

M2
4∆i

ai−1

Zai

M2 4
∆ , (86)
320 i

4

(t − ti ) dt =
ai−1

and also
Zai

Zai

|t − ti | |Ri (t)| dt

(t − ti ) Ri (t)dt ≤

(87)

ai−1

ai−1

Zai

M
≤
2

3

|t − ti | dt =

M 4
∆ .
64 i

ai−1

Here, Hölder’s inequality is attained with equality when
2
λ (t) is proportional to (ϕ0 (t)) · λ21(t) , therefore, the optimal
sampling-point density is
q
2
3
(ϕ0 (t))
λopt (t) = R q
.
(94)
1 3
2
0
(ϕ (z)) dz
0
Setting the lower bound from (93), which is achieved by
λopt (t), in the MSE expression from (90) gives the following
optimal sampling MSE
 1
3
Z q

1
 3 (ϕ0 (t))2 dt .
M SE λopt ≈
(95)
12N 2
0

Plugging (86) and (87) into (85), together with assuming that
|ϕ0 (ti )| < M1 , yields
|Di | ≤

M 2 4 M · M1 3
∆ +
∆i
320 i
32

(88)

implying that the MSE deviation follows

Di = o ∆2i

(89)

A PPENDIX B
S AMPLING -P OINT D ENSITY O PTIMIZATION VIA H ÖLDER ’ S
I NEQUALITY FOR T HE O NE -D IMENSIONAL C ASE
We aim at minimizing the sampling MSE given in (11) as
Z1

2

(ϕ0 (t))
dt.
λ2 (t)

(90)

0

Using Hölder’s inequality we can write
 1
Z q
3
2

(ϕ0 (t)) ·
0

1
2

λ 3 (t)

3

 32
 13  1
Z 
 32
2
λ 3 (t) dt
dt 
0

Z1 q
3

>

2

(ϕ0 (t)) dt.

0

(91)
Simplifying the left side of (91) gives
 1
Z
 (ϕ0 (t))2 ·

 32
 13  1
Z
1
dt  λ (t) dt
λ2 (t)

0

(92)

0

Z1 q
3

>

2

(ϕ0 (t)) dt.

0

Then, since λ (t) is a density function its integration over the
entire domain equals to 1, reducing (92) into
 1
3
Z1
Z q
1
3
2
2
(ϕ0 (t)) · 2 dt > 
(ϕ0 (t)) dt . (93)
λ (t)
0

0

The authors thank the reviewers for very pertinent and
useful comments.
R EFERENCES

as ∆i → 0.

1
M SE (λ) ≈
12N 2

ACKNOWLEDGMENT

[1] W. R. Bennett, “Spectra of quantized signals,” Bell System Technical
Journal, vol. 27, no. 3, pp. 446–472, 1948.
[2] P. Panter and W. Dite, “Quantization distortion in pulse-count modulation with nonuniform spacing of levels,” Proceedings of the IRE, vol. 39,
no. 1, pp. 44–48, 1951.
[3] A. Gersho, “Asymptotically optimal block quantization,” IEEE Transactions on information theory, vol. 25, no. 4, pp. 373–380, 1979.
[4] A. J. Jerri, “The Shannon sampling theorem – Its various extensions
and applications: A tutorial review,” Proceedings of the IEEE, vol. 65,
no. 11, pp. 1565–1596, 1977.
[5] K. Horiuchi, “Sampling principle for continuous signals with timevarying bands,” Information and Control, vol. 13, no. 1, pp. 53–61,
1968.
[6] J. Clark, M. Palmer, and P. Lawrence, “A transformation method for the
reconstruction of functions from nonuniformly spaced samples,” IEEE
Transactions on Acoustics, Speech, and Signal Processing, vol. 33, no. 5,
pp. 1151–1165, 1985.
[7] N. Brueller, N. Peterfreund, and M. Porat, “Non-stationary signals: optimal sampling and instantaneous bandwidth estimation,” in Proceedings
of the IEEE-SP International Symposium on Time-Frequency and TimeScale Analysis, 1998, pp. 113–115.
[8] D. Wei and A. V. Oppenheim, “Sampling based on local bandwidth,” in
Conference Record of the Forty-First Asilomar Conference on Signals,
Systems and Computers (ACSSC), 2007, pp. 1103–1107.
[9] G. J. Sullivan, J. Ohm, W.-J. Han, and T. Wiegand, “Overview of the
high efficiency video coding (HEVC) standard,” IEEE Transactions on
circuits and systems for video technology, vol. 22, no. 12, pp. 1649–
1668, 2012.
[10] R. Shukla, P. L. Dragotti, M. N. Do, and M. Vetterli, “Rate-distortion
optimized tree-structured compression algorithms for piecewise polynomial images,” IEEE Transactions on Image Processing, vol. 14, no. 3,
pp. 343–359, 2005.
[11] E. Shusterman and M. Feder, “Image compression via improved quadtree
decomposition algorithms,” IEEE Transactions on Image Processing,
vol. 3, no. 2, pp. 207–215, 1994.
[12] R. A. Ulichney, “Review of halftoning techniques,” in Electronic Imaging, 1999, pp. 378–391.
[13] D. L. Lau and R. Ulichney, “Blue-noise halftoning for hexagonal grids,”
IEEE Transactions on Image Processing, vol. 15, no. 5, pp. 1270–1284,
2006.
[14] D. P. Mitchell, “Generating antialiased images at low sampling densities,” in ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, 1987, pp.
65–72.
[15] ——, “Spectrally optimal sampling for distribution ray tracing,” in ACM
SIGGRAPH Computer Graphics, vol. 25, no. 4, 1991, pp. 157–164.
[16] M. McCool and E. Fiume, “Hierarchical poisson disk sampling distributions,” in Proceedings of the Conference on Graphics Interface, vol. 92,
1992, pp. 94–105.
[17] R. Fattal, “Blue-noise point sampling using kernel density model,” in
ACM Transactions on Graphics (TOG), vol. 30, no. 4, 2011, p. 48.

14

[18] Z. Zhong and J. Hua, “Kernel-based adaptive sampling for image reconstruction and meshing,” Computer Aided Geometric Design, vol. 43,
pp. 68–81, 2016.
[19] S. Lloyd, “Least squares quantization in PCM,” IEEE Transactions on
Information Theory, vol. 28, no. 2, pp. 129–137, 1982.
[20] D. L. Neuhoff and S. S. Pradhan, “Information rates of densely sampled
data: Distributed vector quantization and scalar quantization with transforms for Gaussian sources,” IEEE Transactions on Information Theory,
vol. 59, no. 9, pp. 5641–5664, 2013.
[21] A. Kipnis, A. J. Goldsmith, and Y. C. Eldar, “Gaussian distortionrate function under sub-Nyquist nonuniform sampling,” in 52nd Annual
Allerton Conference on Communication, Control, and Computing, 2014,
pp. 874–880.
[22] A. Kipnis, A. J. Goldsmith, Y. C. Eldar, and T. Weissman, “Distortion
rate function of sub-Nyquist sampled Gaussian sources,” IEEE Transactions on Information Theory, vol. 62, no. 1, pp. 401–429, 2016.
[23] S. Feizi, V. K. Goyal, and M. Médard, “Locally adaptive sampling,”
in 48th Annual Allerton Conference on Communication, Control, and
Computing, 2010, pp. 152–159.
[24] ——, “Time-stampless adaptive nonuniform sampling for stochastic
signals,” IEEE Transactions on Signal Processing, vol. 60, no. 10, pp.
5440–5450, 2012.
[25] R. M. Gray and D. L. Neuhoff, “Quantization,” IEEE Transactions on
Information Theory, vol. 44, no. 6, pp. 2325–2383, 1998.
[26] D. Newman, “The hexagon theorem,” IEEE Transactions on Information
Theory, vol. 28, no. 2, pp. 137–139, 1982.
[27] E. Barnes and N. Sloane, “The optimal lattice quantizer in three
dimensions,” SIAM Journal on Algebraic Discrete Methods, vol. 4, no. 1,
pp. 30–41, 1983.
[28] Q. Du and D. Wang, “The optimal centroidal voronoi tessellations and
the gersho’s conjecture in the three-dimensional space,” Computers &
Mathematics with Applications, vol. 49, no. 9, pp. 1355–1373, 2005.
[29] P. Zador, “Asymptotic quantization error of continuous signals and the
quantization dimension,” IEEE Transactions on Information Theory,
vol. 28, no. 2, pp. 139–149, 1982.
[30] S. Na and D. L. Neuhoff, “Bennett’s integral for vector quantizers,”
IEEE Transactions on Information Theory, vol. 41, no. 4, pp. 886–900,
1995.
[31] P. A. Chou, T. Lookabaugh, and R. M. Gray, “Optimal pruning with
applications to tree-structured source coding and modeling,” IEEE
Transactions on Information Theory, vol. 35, no. 2, pp. 299–315, 1989.
[32] H. Everett III, “Generalized lagrange multiplier method for solving
problems of optimum allocation of resources,” Operations Research,
vol. 11, no. 3, pp. 399–417, 1963.
[33] Y. Shoham and A. Gersho, “Efficient bit allocation for an arbitrary
set of quantizers,” IEEE Transactions on Acoustics, Speech, and Signal
Processing, vol. 36, no. 9, pp. 1445–1453, 1988.
[34] A. Ortega and K. Ramchandran, “Rate-distortion methods for image and
video compression,” IEEE Signal Processing Magazine, vol. 15, no. 6,
pp. 23–50, 1998.
[35] J. Bucklew, “Companding and random quantization in several dimensions,” IEEE Transactions on Information Theory, vol. 27, no. 2, pp.
207–211, 1981.
[36] ——, “A note on optimal multidimensional companders (corresp.),”
IEEE Transactions on Information Theory, vol. 29, no. 2, pp. 279–279,
1983.
[37] P. Moo and D. Neuhoff, “Optimal compressor functions for multidimensional companding,” in IEEE International Symposium on Information
Theory, 1997, p. 515.
[38] S. F. Simon, “On suboptimal multidimensional companding,” in Data
Compression Conference, 1998. DCC’98. Proceedings, 1998, pp. 438–
447.
[39] J. Samuelsson, “Multidimensional companding quantization of the gaussian source,” IEEE Transactions on Information Theory, vol. 49, no. 5,
pp. 1343–1351, 2003.
[40] Z. H. Perić, O. D. Milanović, and A. Ž. Jovanović, “Optimal companding vector quantization for circularly symmetric sources,” Information
Sciences, vol. 178, no. 22, pp. 4375–4381, 2008.

