# Time Encoding via Unlimited Sampling:
Theory, Algorithms and Hardware Validation

Dorian Florescu and Ayush Bhandari This work is supported by the UK Research
and Innovation (UKRI) council’s _Future Leaders Fellowship_ program “Sensing
Beyond Barriers” (MRC Fellowship award no. MR/S034897/1). Project page for
(future) release of hardware design, code and data: https://bit.ly/USF-Link.D.
Florescu and A. Bhandari are with the Department of Electrical and Electronic
Engineering, Imperial College London, SW72AZ, UK. E-mails: {D.Florescu,
<EMAIL_ADDRESS>or<EMAIL_ADDRESS>received on Feb.
XX, 2022; accepted with minor revisions on August XX, 2022.

###### Abstract

An alternative to conventional uniform sampling is that of time encoding,
which converts continuous-time signals into streams of trigger times. This
gives rise to Event-Driven Sampling (EDS) models. The data-driven nature of
EDS acquisition is advantageous in terms of power consumption and time
resolution and is inspired by the information representation in biological
nervous systems. If an analog signal is outside a predefined dynamic range,
then EDS generates a low density of trigger times, which in turn leads to
recovery distortion due to aliasing. In this paper, inspired by the Unlimited
Sensing Framework (USF), we propose a new EDS architecture that incorporates a
modulo nonlinearity prior to acquisition that we refer to as the modulo EDS or
MEDS. In MEDS, the modulo nonlinearity folds high dynamic range inputs into
low dynamic range amplitudes, thus avoiding recovery distortion. In
particular, we consider the asynchronous sigma-delta modulator (ASDM),
previously used for low power analog-to-digital conversion. This novel MEDS
based acquisition is enabled by a recent generalization of the modulo
nonlinearity called modulo-hysteresis. We design a mathematically guaranteed
recovery algorithm for bandlimited inputs based on a sampling rate criterion
and provide reconstruction error bounds. We go beyond numerical experiments
and also provide a first hardware validation of our approach, thus bridging
the gap between theory and practice, while corroborating the conceptual
underpinnings of our work.

###### Index Terms:

Event-driven, nonuniform sampling, modulo sampling, analog-to-digital
conversion (ADC).

###### Contents

  1. I Introduction
  2. II The Asynchronous Sigma-Delta Modulator
    1. II-A Model Description
    2. II-B Reconstruction from Asynchronous Sigma-Delta Samples
  3. III The Proposed Modulo Event-Driven Encoder
    1. III-A The Ideal Modulo
    2. III-B The Modulo-Hysteresis Model
    3. III-C The Proposed Modulo Event-Driven Sampling Model
  4. IV Recovery from Modulo Event-Driven Data
    1. IV-A The Nonuniform Difference Operator
    2. IV-B The Case of a Single Fold
    3. IV-C The Case of Multiple Folds
    4. IV-D Input Reconstruction
  5. V Numerical and Hardware Experiments
    1. V-A Reconstruction using Synthetic Data
    2. V-B Reconstruction using Experimental Data
  6. VI Proofs
  7. VII Conclusions and Future Work

##  I Introduction

In Shannon’s sampling paradigm, acquisition is performed by recording the
amplitude of a signal at predefined, uniform time locations. Alternatively, in
_time encoding_ the signal is converted into a nonuniform sequence of time
events, leading to the acquisition paradigm called Event-Driven Sampling (EDS)
[1, 2, 3, 4].

Such time events are induced by _significant_ changes in the input signal
values, thus enabling a data-driven approach to sampling. EDS has been adopted
in engineering fields such as control engineering and signal processing [5]
and has a wide range of applications, including neuromorphic vision [6],
machine learning [7, 8] and brain-machine interfaces (BMIs) [9]. In biology,
EDS is used to model the signal transmission in the nervous system of
vertebrates [10]. When compared to uniform sampling, event-driven sampling
schemes are not controlled by an external clock signal and are characterized
by low power consumption [11, 12] and increased time resolution [6].

Figure 1: The output of the proposed modulo event-driven sampling (MEDS) model
and the standalone EDS implemented as an asynchronous sigma-delta modulator
(ASDM) for a high dynamic range input. (a) The MEDS model folds the input
amplitude thus preventing saturation. (b) The standalone ASDM has a suppressed
output when the input amplitude is large. Our MEDS hardware prototype is
depicted in Fig. 5.

Bottleneck in the Current EDS Acquisition. In EDS, the input is encoded into
an asynchronous sequence of time events or trigger times. The information
content in such a sequence can be characterised by the sampling density, a
positive quantity measuring the number of trigger instants per unit time. This
leads to an encoding barrier: no information is being transmitted for inputs
whose amplitudes are outside the range where the sampling density is positive.
The barrier is illustrated for the particular case of the asynchronous sigma-
delta modulator based EDS in Fig. 1. Therefore, the accurate recovery of the
input is affected, which requires a minimal sampling density to prevent
aliasing [13]. Thus, in EDS, the input is typically constrained to a
predefined dynamic range.

The problem of limited dynamic range is common in traditional encoding schemes
such as analog-to-digital converters (ADCs), which clip or saturate the input
signals exceeding a predefined threshold. To overcome this fundamental
bottleneck, the Unlimited Sensing Framework (USF) [14, 15, 16, 17, 18] was
introduced and recently validated in hardware via modulo ADCs [18]—it was
shown that, in practice, signals up to $24$ times larger than ADC’s threshold
can be recovered using modulo ADCs [18], despite non-idealities and
quantization effects. In USF, the introduction of modulo nonlinearity converts
a signal from high dynamic range (HDR) to low dynamic range. This is done by
_folding_ the HDR, continuous-time signal back into the dynamic range of the
ADC, whenever the input amplitude is outside a predefined range [14, 15, 16,
17, 18]. The follow-up work on USF includes,

  * •

Different signal classes, such as sparse signals [19, 20], sinusoidal mixtures
[21] and spline spaces [22].

  * •

New sampling architectures [23, 24, 25, 26].

  * •

Application frontiers such as array signal processing [27, 28] and
computational imaging [22, 29, 30].

A more general model for modulo nonlinearity was introduced recently, namely,
_modulo-hysteresis_[25, 23], which, compared to ideal modulo, comprises two
additional parameters that can tackle non-idealities present in experimental
observations. Moreover, modulo-hysteresis provided guarantees for a new class
of reconstruction methods based on thresholding.

To address the dynamic range limitation in event-driven sampling schemes, we
propose a new encoding model consisting of a modulo-hysteresis nonlinearity in
series with an EDS model, which we call modulo EDS or MEDS. In our setting the
EDS model is fixed, and the EDS input (modulo output) is guaranteed to satisfy
the EDS dynamic range by selecting an appropriate modulo nonlinearity (Fig.
1(a)). None of the existing recovery methods for EDS models can address the
reconstruction of the modulo output, which is a folded signal. Therefore, the
use of advanced recovery algorithms enables high dynamic range acquisition via
MEDS that is conventionally not possible.

The EDS model in this paper is an asynchronous sigma-delta modulator (ASDM)
[1]. Thanks to its low power consumption [31] and modular design [32], the
ASDM was first implemented as an alternative to conventional ADCs in [33] and
subsequently included in applications such as brain-machine interfaces [9].
Assuming that the input amplitude is in a predefined range of values, the ASDM
transmits trigger times measuring changes in the input integral [1]. For
inputs of amplitude outside this range, the ASDM transmits no information, as
shown in Fig. 1(b).

Related Work. The ASDM dynamic range can be extended by changing the model
parameters or architecture. This includes increasing the sampling density by
adjusting the triggering threshold [1], or introducing an adaptation mechanism
via a new ASDM architecture [34]. Another popular EDS model is the integrate-
and-fire (IF) neuron, inspired from the biological neuron, for which input
reconstruction was studied for a wide range of input classes [3, 35, 36]
including recovery methods for EDS leveraging uniform sampling theory [37].
The recovery methods for IF also assume that the input is in a predefined
dynamic range. It is possible to encode high dynamic range signals with an IF
by generating bidirectional event-driven samples [35, 38]. In this case,
however, recovery guarantees do not exist for low amplitude inputs, for which
the model does not generate output. Event based cameras, which use models
closely related to the IF or ASDM, capture high dynamic range data with
logarithmic pixels. This, in turn, decreases the pixel resolution for high
intensities[6]. None of the methods above can recover modulo folded inputs
from the output of an EDS. A MEDS architecture that extends the IF dynamic
range was first numerically validated in [25].

Contributions. Our contributions are as follows,

  1. $\bullet$

We introduce a new event-driven modulo sampling scheme, comprising a modulo-
hysteresis model in conjunction with an ASDM, which addresses the dynamic
range restriction for an ASDM model without any parameter or architecture
alteration. The new scheme is fully compatible with the existing ASDM
methodology, but can also accommodate inputs that do not satisfy the dynamic
range requirement.

  2. $\bullet$

By exploiting the modulo-hysteresis architecture, we provide theoretical
conditions for which the input can be recovered from the output trigger times.

  3. $\bullet$

Extensive numerical experiments and a validation of the MEDS-based hardware
prototype (cf. Fig. 5) show the validity of our approach.

Notation. We use $\mathbb{Z}$ and $\mathbb{R}$ to denote the set of integers
and real numbers, respectively, and $\mathbb{N}$ to denote the set of positive
integers. Continuous functions are denoted as $f\left(t\right)$ and discrete
sequences as $f\left[k\right]$. The Hilbert space of square-integrable
functions is denoted as $L^{2}\left(\mathbb{R}\right)$. The norm in any space
$H$ is denoted as $\|f\|_{H}$. The derivative of order $N$ is denoted as
$f^{\left(N\right)}\left(t\right)$ and, for sequences, the finite difference
of order $N$ is $\Delta^{N}f\left[k\right]$, which is computed by applying
recursively
$\Delta^{N+1}f\left[k\right]=\Delta^{N}f\left[k+1\right]-\Delta^{N}f\left[k\right]$,
where $\Delta^{0}f=f$. The space of square integrable functions bandlimited to
$\Omega$ is the Paley-Wiener space denoted $\mathsf{PW}_{\Omega}$. The
indicator function $\mathbbmtt{1}_{S}\left(t\right)$ is $1$ for $t\in S$ and
$0$ otherwise. The floor and ceiling functions are
$\left\lfloor\cdot\right\rfloor$ and $\left\lceil\cdot\right\rceil$,
respectively, and $[\\![x]\\!]=x-\lfloor x\rfloor$ denotes the fractional part
of $x$. Furthermore, let $\mathrm{supp}\left(f\right)$ denote the support of
sequence $f\left[k\right]$, let $\left\lvert S\right\rvert$ denote the
cardinality of a set $S$, and let $\emptyset$ denote the empty set. By
$f_{\infty}$ or $\|f\|_{\infty}$ we denote the absolute sequence norm
$\|f\|_{\ell^{\infty}}$ or the absolute function norm $\|f\|_{L^{\infty}}$,
respectively. Let $\mathrm{sinc}_{\Omega}$ be the sinc function defined here
as $\mathrm{sinc}_{\Omega}\left(t\right)\triangleq\frac{\sin\left(\Omega
t\right)}{\pi t}$.

Scope and Organization. The ASDM model and associated reconstruction methods
are summarised in Section II. The MEDS architecture is introduced in Section
III. The new reconstruction method from MEDS samples and associated
theoretical guarantees are given in Section IV. The numerical demonstration
and MEDS hardware prototype are in Section V. Section VI comprises proofs for
some of the theoretical results in Section IV, and the conclusions are in
Section VII.

##  II The Asynchronous Sigma-Delta Modulator

The ASDM is an event-driven sampling model consisting of a feedback loop
containing an adder, an integrator, and a noninverting Schmitt trigger, as
illustrated in Fig. 2. We introduce the formal definition of the ASDM model in
Section II-A and summarise the recovery method for bandlimited inputs and its
shortcomings in Section II-B.

Figure 2: The asynchronous sigma-delta modulator (ASDM) adds input
$g\left(t\right)$ with $z\left(t\right)$, the current state of the ASDM, and
the result is processed with an integrator with initial condition $y_{0}$ in
series with a Schmitt trigger with parameters $\delta,b$.

Figure 3: The effect of ASDM saturation observed in hardware measurements (cf.
Fig. 5(b)). Input $g(t)$ goes below the ASDM dynamic range, thus the
integrator output $y(t)$ is no longer monotonic. The next trigger time is only
generated when $g(t)$ is back in the dynamic range which pushes $y(t)$ to
cross threshold $-b$. This illustration is based on the experiment reported in
[39].

###  II-A Model Description

At the initial time point $t{=t_{0}}=0$, the output of the Schmitt trigger is
set as $z\left(t\right)=-b$. The signal $z\left(t\right)$ retains this state
until the first trigger time occurs, as will be explained next. In a vicinity
of the $t=0$ the integrator output $y(t)$ satisfies:

$y\left(t\right)=y_{0}+\int_{0}^{t}\left(g\left(s\right)+b\right)ds,$ (1)

where $g\left(t\right)$ is the input, $y_{0}$ is a constant denoting the
integrator’s initial condition, and $b,\delta$ are the Schmitt trigger
parameters. The ASDM generates output only if $|g\left(t\right)|\leqslant
c<b$, which implies that $y\left(t\right)$ is strictly increasing (1). When
$|g\left(t\right)|>c$, $y\left(t\right)$ is no longer monotonic, and thus may
not cross the threshold $\delta$. This saturates the ASDM, preventing any
triggers until the input is back in the ASDM dynamic range (see Fig. 3). The
trigger time $t_{1}$ corresponds to $y\left(t_{1}\right)=\delta$, which flips
the output of the Schmitt trigger to $b$. This enforces $y\left(t\right)$ to
be strictly decreasing, and the next trigger time satisfies
$y\left(t_{2}\right)=-\delta$. The sequence of trigger times
$\left\\{t_{k}\right\\}$ is directly linked with the input and ASDM parameters
via the following recursive equation, also known as the _t-transform_ [1]

$\int_{t_{k}}^{t_{k+1}}\left(g\left(s\right)+\left(-1\right)^{k}b\right)ds=\left(-1\right)^{k}2\delta,\quad{k\in\mathbb{Z}_{+}^{*}}.$
(2)

We model the link between the trigger times and input via the operator
$\left\\{t_{k}\right\\}=\mathrm{ASDM}_{\delta,b}\left(g\right)$. The
$t$-transform was also studied in the context of other EDS models such as the
IF neuron [2]. By defining
$G\left(t\right)\triangleq\int_{0}^{t}g\left(s\right)ds$, it can be shown via
induction that $G\left(t_{k}\right)$ can be derived directly from (2). This
creates a strong link between ASDM sampling and nonuniform sampling, and
allows using similar analysis and recovery approaches [13]. A hardware
realisation of an ASDM is in Fig. 5(b).

###  II-B Reconstruction from Asynchronous Sigma-Delta Samples

The ASDM input $g(t)$ can be perfectly reconstructed [2] if it is bandlimited
$g\in\mathsf{PW}_{\Omega}$ and the maximum distance between the ASDM trigger
times satisfies the Nyquist rate condition, i.e.,

$\left\lvert t_{k+1}-t_{k}\right\rvert<\frac{\pi}{\Omega}.$

This condition is sufficient to guarantee recovery from nonuniform samples
$G\left(t_{k}\right)$ [13].

In the case of the ASDM, if $\left\lvert g(t)\right\rvert>b$ the encoder does
not generate any output, meaning that $\left\lvert
t_{k+1}-t_{k}\right\rvert<\frac{\pi}{\Omega}$ does not always hold. To
overcome this problem it is assumed that $g_{\infty}<b$, which , via (2),
guarantees that $t_{k+1}-t_{k}<\frac{2\delta}{b-g_{\infty}}$. Therefore, the
Nyquist rate condition for recovery is guaranteed if

$g_{\infty}<b-\frac{2\delta\Omega}{\pi}.$ (3)

The input is recovered as [1]

$g\left(t\right)=\sum\limits_{m\in\mathbb{Z}}c_{m}\mathrm{sinc}_{\Omega}\left(t-\bar{t}_{m}\right),$
(4)

where $\bar{t}_{m}=\left(t_{m}+t_{m+1}\right)/{2}$ are the midpoints of the
intervals between the trigger times. The coefficients $c_{m}$ are computed via
least squares from

$\int_{t_{k}}^{t_{k+1}}g\left(s\right)ds\stackrel{{\scriptstyle(\ref{eq:ttransform})}}{{=}}\left(-1\right)^{k}2\delta-\left(-1\right)^{k}b\left(t_{k+1}-t_{k}\right).$

In the context of nonuniform sampling, this procedure is also known as
reconstruction from local averages [13].

Intuitively, whenever $\left\lvert g\left(t\right)\right\rvert>b$, the input
does not trigger the ASDM (Fig. 1(b)). This is because the integrator output
$y\left(t\right)$ is no longer monotonic and might not reach the threshold
$\pm b$. This locks in the ASDM output until the threshold is reached. This
effect is also observed in hardware, as illustrated in Fig. 3. As will be
shown in Section V, the ASDM saturation leads to unstable reconstructions.

Figure 4: The proposed encoding model.

##  III The Proposed Modulo Event-Driven Encoder

Here we introduce a new event-driven modulo encoder model, comprising a modulo
encoder with hysteresis and an ASDM model, depicted in Fig. 4. We first
present the ideal modulo encoder in III-A and explain some of its limitations.
The modulo-hysteresis is then described in III-B. The proposed event-driven
modulo is introduced in III-C.

###  III-A The Ideal Modulo

Let $g\in\mathsf{PW}_{\Omega}$ where $\mathsf{PW}_{\Omega}$ is the Paley-
Wiener space of $\Omega$–bandlimited functions. As in [19], the modulo
nonlinearity is defined by function
$\mathscr{M}_{\lambda}:\mathsf{PW}_{\Omega}\rightarrow
L^{2}\left(\mathbb{R}\right)$

$\mathscr{M}_{\lambda}\left(g\left(t\right)\right)=2\lambda\left({\left[\kern-3.00003pt\left[{\frac{g\left(t\right)}{{2\lambda}}+\frac{1}{2}}\right]\kern-3.00003pt\right]-\frac{1}{2}}\right)$
(5)

where $[\\![x]\\!]=x-\lfloor x\rfloor$ is the fractional part of $x$, and
$L^{2}\left(\mathbb{R}\right)$ denotes the space of functions with finite
energy. The instantaneous times when $[\\![x]\\!]$ switches from $1$ to $0$
are called are called _folding times_ and denoted $\tau_{r},r\in\mathbb{Z}$.
The output of the model
$x\left(t\right)=\mathscr{M}_{\lambda}\left(g\left(t\right)\right)$ can be
expanded as,

$x\left(t\right)=g\left(t\right)-\varepsilon_{g}\left(t\right)\mbox{ where
}\varepsilon_{g}\left(t\right)=2\lambda\sum\limits_{r\in\mathbb{Z}}s_{r}\mathbbmtt{1}_{[\tau_{r},\infty)}\left(t\right),t\in\mathbb{R},\quad
s_{r}=\pm 1.$

In Unlimited Sampling the data is sampled uniformly with period $T$, which
yields

$x\left(kT\right)=g\left(kT\right)-\varepsilon_{g}\left(kT\right),\quad
k\in\mathbb{Z}.$ (6)

The input $g\left(kT\right)$ is recovered by applying a finite difference
filter $\Delta^{N}$ to (6), which annihilates
$\Delta^{N}\varepsilon_{g}\left(kT\right)$ given $\Delta^{N}x\left(kT\right)$.
This is possible because $\Delta^{N}$ vanishes gradually
$\Delta^{N}g\left(kT\right)$ but always leaves
$\Delta^{N}\varepsilon_{g}\left(kT\right)$ on an equally spaced grid that is
then annihilated via a modulo operation. However this approach does not work
for experimental modulo data due to various non-idealities [18, 23]. This
approach also would not work when sampling the output with an ASDM, which
integrates the input on intervals of various durations. (2).

Figure 5: Hardware experiment. (a) The US-ADC hardware prototype with variable
parameters controlling the threshold $\lambda$ and hysteresis $h$. (b) The
ASDM hardware realisation. (c) Oscilloscope screenshot showing the input and
outputs of the modulo-hysteresis and ASDM hardware prototypes.

Figure 6: Encoding a random bandlimited waveform with two types of modulo
nonlinearities. (a) The ideal modulo resets the input to the opposite
threshold value, and can generate arbitrarily close folding times [14, 19].
(b) The modulo-hysteresis nonlinearity separates the thresholds and post-reset
values to guarantee a minimal separation between the folding times.

A different reconstruction approach from uniform samples for the ideal modulo
consists of filtering the modulo output with a kernel that amplifies the
folding times, which are subsequently detected via _thresholding_ [40, 24].
The thresholding approach is particularly desirable when the data is affected
by non-idealities, such as transient folding samples, where the Unlimited
Sampling approach fails [23]. However, the input reconstruction is not
guaranteed when the folding times are very close, which can happen for the
ideal modulo, as depicted in Fig. 6(a). This problem can be solved using the
modulo-hysteresis model [23], presented in the next subsection, which will be
used in the present work.

###  III-B The Modulo-Hysteresis Model

The phenomenon of _hysteresis_ , present in ADC circuits [41], results from a
difference between the quantization threshold values when the value is
approached from above or below [42]. This effect is also commonly exploited in
circuits such as the Schmitt trigger, which avoids noise induced spurious
triggering during ADC conversion. Recently, that effect was exploited in a
modulo-hysteresis model and associated hardware implementation [23] see Fig.
5(a), which limits the minimum distance between two consecutive folding times.
We give below the formal definition of modulo-hysteresis.

######  Definition 1 (Modulo-Hysteresis)

The modulo-hysteresis with threshold $\lambda$ and hysteresis parameter
$h\in\left[0,\lambda\right)$, denoted $\mathscr{M}_{\boldsymbol{\mathsf{H}}}$,
$\boldsymbol{\mathsf{H}}=\left[\lambda\ h\right]$, generates an nonuniform
sequence $\left\\{\tau_{r}\right\\}$ and an analog function
$x\left(t\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\left(t\right),$
$t\geqslant\tau_{0}$ in response to input $g\in\mathsf{PW}_{\Omega}$. The
sequence $\left\\{\tau_{r}\right\\}$ is computed recursively as
$\displaystyle\tau_{1}$
$\displaystyle=\min\left\\{t>\tau_{0}|\mathscr{M}_{\lambda}\left(g\left(t\right)+\lambda\right)=0\right\\}$
$\displaystyle\tau_{r+1}$
$\displaystyle=\min\left\\{t>\tau_{r}|\mathscr{M}_{\lambda}\left(g\left(t\right)-g\left(\tau_{r}\right)+hs_{r}\right)=0\right\\},$
where
$s_{r}=\mathrm{sign}\left(g\left(\tau_{r}\right)-g\left(\tau_{r-1}\right)\right),n\geqslant
1$. The function $x\left(t\right)$ returned by the encoder is then given by
$x\left(t\right)=g\left(t\right)-\varepsilon_{g}\left(t\right),$ (7) where
$\varepsilon_{g}\left(t\right)=2\lambda_{h}\sum\nolimits_{r\in\mathbb{Z}}s_{r}\mathbbmtt{1}_{[\tau_{r},\infty)}\left(t\right)$
and $\lambda_{h}\triangleq\lambda-h/2$.

Modulo-hysteresis has the following important property.

######  Lemma 1 (Minimum Separation between Folds)

Let $\left\\{\tau_{r}\right\\}_{r\in\mathbb{Z}}$ be the sequence of folding
times corresponding to encoder $\mathscr{M}_{\boldsymbol{\mathsf{H}}}$,
$\boldsymbol{\mathsf{H}}=\left[\lambda\ h\right]$, for an input
$g\in\mathsf{PW}_{\Omega}$. Then

$\tau_{r+1}-\tau_{r}\geqslant\frac{h^{*}}{\Omega g_{\infty}},\quad
h^{*}=\min\left\\{h,2\lambda_{h}\right\\}.$ (8)

###### Proof.

From Bernšteĭn’s inequality
$\|g^{\left(n\right)}\|_{p}\leqslant\Omega^{n}\|g\|_{p}$ (cf. pg. 116, [43])
and Definition 7, for $p=\infty$, $n=1$ and $r\geqslant 1$,

$\tau_{r+1}-\tau_{r}\geqslant\frac{\left\lvert
g\left(\tau_{r+1}\right)-g\left(\tau_{r}\right)\right\rvert}{\max\left\lvert
g^{\prime}\left(t\right)\right\rvert}\geqslant\frac{\min\left\\{h,2\lambda-h\right\\}}{\Omega\max\left\lvert
g\left(t\right)\right\rvert}=\frac{h^{*}}{\Omega g_{\infty}}.$

This proves the result. ∎

When the data is sampled with period $T$, it was shown that here
$g\left(kT\right)$ can still be recovered from
$x\left(kT\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\left(kT\right)$ with
USF, given that $\Delta^{N}\varepsilon_{g}\left(kT\right)$ is on a uniform
grid with spacing $2\lambda-h$ [23]. This property no longer holds when the
data is sampled with an asynchronous sigma-delta modulator, as will be
explained next.

###  III-C The Proposed Modulo Event-Driven Sampling Model

The proposed model consists of a modulo-hysteresis in series with an ASDM
model. Then, for $k\geqslant 0$, the trigger times $\left\\{t_{k}\right\\}$
generated by the model satisfy (2)

$\begin{gathered}\int_{t_{k}}^{t_{k+1}}\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\left(s\right)ds=q_{k},\\\
q_{k}\triangleq\left(-1\right)^{k}2\delta-\left(-1\right)^{k}b\left(t_{k+1}-t_{k}\right).\end{gathered}$
(9)

We assume $g\in\mathsf{PW}_{\Omega}$ and $\lambda\leqslant
b-\frac{2\delta\Omega}{\pi}$ such that if $\left\lvert
g\left(t\right)\right\rvert<\lambda$ then (3) is satisfied, and thus $g$ is
left unchanged by the modulo nonlinearity. This ensures that (9) is equivalent
to (2), and thus the proposed encoder is backwards compatible with the
existing results on event-driven sampling. When $\left\lvert
g\left(t\right)\right\rvert>b-\frac{2\delta\Omega}{\pi}$, the asynchronous
sigma-delta modulator leads to information loss, which is prevented by the
proposed circuit by folding the input back into the range
$\left[-\lambda,\lambda\right]$. According to Definition 7,
$x\left(t\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\left(t\right)$
satisfies

$x\left(t\right)\stackrel{{\scriptstyle(\ref{eq:zt})}}{{=}}g\left(t\right)-2\lambda_{h}\sum\limits_{r\in\mathbb{Z}}s_{r}\mathbbmtt{1}_{\left[\tau_{r},\infty\right)}\left(t\right).$
(10)

Let $X\left(t\right)=\int_{{0}}^{t}x\left(s\right)ds$. Then
${X\left(t_{0}\right)=X\left(0\right)=0}$ and (9) implies that

$X\left(t_{k}\right)=\sum\limits_{m={0}}^{k-1}q_{m},\quad k\geqslant{1}.$ (11)

Furthermore, let
$E_{g}\left(t\right)\triangleq\int_{0}^{t}\varepsilon_{g}\left(s\right)ds$ and
$G\left(t\right)\triangleq\int_{0}^{t}g\left(s\right)ds$, and thus, by
integrating (10) we get

$X\left(t_{k}\right)=G\left(t_{k}\right)-E_{g}\left(t_{k}\right).$ (12)

We can then compute samples $X\left(t_{k}\right)$ using just the model
parameters and output trigger times via (11). Further on, the input and
residue will be separated using (12). Equation (12) has a similar form as (6),
but here $E_{g}\left(t_{k}\right)$ no longer lie on an uniform grid as it
involves integration on intervals of variable length. Therefore, a different
recovery strategy than in the case of USF is required, which is explained in
the following.

##  IV Recovery from Modulo Event-Driven Data

Our goal is to recover $g\left(t\right)$ from (9) given time instances
$\left\\{t_{k}\right\\}_{k\in\mathbb{Z}}$ and design parameters
$\left\\{\lambda,h,b,\delta\right\\}$. The proposed encoding circuit
guarantees the asynchronous sigma-delta modulator input is in a predefined
range. Generally, the EDS input $x\left(t\right)$ is not bandlimited, and does
not even belong to a class of shift-invariant spaces, given that
${\tau_{r}\geqslant 0}$ do not lie on a uniform grid. Thus $x\left(t\right)$
cannot be recovered with the traditional methods for EDS models, which assume
that the input is bandlimited, or at least smooth [2, 44, 3, 23].

Given samples $X\left(t_{k}\right)=G\left(t_{k}\right)-E\left(t_{k}\right)$,
the main challenge is to recover both $G\left(t_{k}\right)$ and
$E\left(t_{k}\right)$. By applying a generic filter $\mathscr{D}$ to the data
$X\left(t_{k}\right)$ for $k\geqslant 2$ we get (7, 11)

$\mathscr{D}X\left(t_{k}\right)=\mathscr{D}G\left(t_{k}\right)-\mathscr{D}E_{g}\left(t_{k}\right)=\mathscr{D}\left(\sum\limits_{m={0}}^{k-1}q_{m}\right).$
(13)

The input reconstruction is organised in two stages. In the first stage,
$\mathscr{D}G\left(t_{k}\right)$ is interpreted as noise and an estimation
$\widetilde{E}_{g}\left(t_{k}\right)$ of $E_{g}\left(t_{k}\right)$ is computed
by thresholding $\mathscr{D}X\left(t_{k}\right)$. This approach was used
before to reconstruct modulo folded data in the context of uniform samples
[40, 24, 23]. In the second recovery stage, an estimation
$\widetilde{G}\left(t_{k}\right)$ of $G\left(t_{k}\right)$ is computed from
the integral form of the modulo-hysteresis input $X\left(t_{k}\right)$ and
$\widetilde{E}_{g}\left(t_{k}\right)$ through (12), and the respective
estimate for the continuous input $\widetilde{g}\left(t\right)$ is
reconstructed from the estimated samples $\widetilde{G}\left(t_{k}\right)$,
such that the error norm $\|g-\widetilde{g}\|_{L^{2}}$ is bounded.

The section is organised as follows. To accommodate nonuniform samples, we
define the nonuniform difference operator $\mathscr{D}^{N}$ in Section IV-A.
We then describe the first recovery stage, and provide guarantees for
computing $E_{g}\left(t_{k}\right)$ in the particular case of one folding time
in Section IV-B, and extend it to multiple folding times in Section IV-C.
Section IV-D describes the second recovery stage and gives the input
reconstruction algorithm and associated error bound.

###  IV-A The Nonuniform Difference Operator

Here we formulate a thresholding based recovery for the case of the ASDM
model. Given samples
$X\left(t_{k}\right)=G\left(t_{k}\right)-E\left(t_{k}\right)$, the main
challenge is to recover both $G\left(t_{k}\right)$ and $E\left(t_{k}\right)$.
We threshold samples $\mathscr{D}X\left(t_{k}\right)$ to recover the residual
$E_{g}\left(t\right)$, which requires identifying the folding times $\tau_{r}$
and signs $s_{r}$ (10). Operator $\mathscr{D}$ is designed to enhance the
effect of $E_{g}$ in (13). In [23] the data was processed with the uniform
finite difference filter $\Delta^{N}$. Here we recover the input from
nonuniform samples and, to this end, we choose $\mathscr{D}=\mathscr{D}^{N}$
defined recursively the nonuniform difference operator of order $N$,

$\mathscr{D}^{N}f\left[k\right]=\frac{\mathscr{D}^{N-1}f\left[k+1\right]-\mathscr{D}^{N-1}f\left[k\right]}{t_{k+N}-t_{k}},\quad
N>0,$ (14)

where $\mathscr{D}^{0}f\left[k\right]=f\left[k\right]$ and $f\left[k\right]$
is a generic sequence. Here, $\mathscr{D}^{N}$ is an extension of the
difference operator $\Delta^{N}$ used for uniform samples. In (13), the values
$\mathscr{D}^{N}\sum_{m={0}}^{k-1}q_{m}$ can be computed directly from the
measurements. As is the case for USF, here the measurements lead to a mixture
of the samples $\mathscr{D}^{N}G\left(t_{k}\right)$ and
$\mathscr{D}^{N}E_{g}\left(t_{k}\right)$. However, unlike the uniform sampling
scenario, values $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$ do not lie on an
uniform grid. We will show that the separation of
$\mathscr{D}^{N}G\left(t_{k}\right)$ and
$\mathscr{D}^{N}E_{g}\left(t_{k}\right)$ is still possible via a thresholding
approach. To this end, we provide the following bound, which will be used in
our recovery approach (proof in Section VI).

###### Proposition 1

Let $\mathscr{D}^{N}$ be the nonuniform finite difference operator defined in
(14). Then
$\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert\leqslant\frac{1}{\Omega
N!}\left(\frac{T_{\mathsf{max}}}{T_{\mathsf{min}}}\Omega
e\right)^{N}g_{\infty},$ (15) where $T_{\mathsf{min}},T_{\mathsf{max}}>0$ are
two constants satisfying $T_{\mathsf{min}}\leqslant t_{k+1}-t_{k}\leqslant
T_{\mathsf{max}}$.

This bound is an extension to nonuniform samples of the bound derived in [14]
for $\left\lvert\Delta^{N}g\left(kT\right)\right\rvert$. Specifically, we note
that in the limit case where $T_{\mathsf{min}}=T_{\mathsf{max}}=T$, and thus
$t_{k}=kT$, we have that

$\mathscr{D}^{N}G\left(t_{k}\right)\stackrel{{\scriptstyle(\ref{eq:NFD})}}{{=}}\frac{\Delta^{N}G\left(kT\right)}{T^{N}N!}.$

According to Bernšteĭn’s we have $\|G^{\prime}\|_{\infty}\leqslant\Omega
G_{\infty}\Rightarrow g_{\infty}\leqslant\Omega G_{\infty}$. Then, in this
particular case, (15) becomes
$\left\lvert\Delta^{N}G\left(kT\right)\right\rvert\leqslant\left(T\Omega
e\right)^{N}G_{\infty}$, which is in accordance to the bound derived in [14,
16].

###  IV-B The Case of a Single Fold

For simplicity, let us first assume we have only one folding time $\tau_{1}$.
Let $K_{1}$ be the index of the trigger time preceding the folding time, such
that $\tau_{1}\in\left[t_{K_{1}},t_{K_{1}+1}\right)$. The following result
calculates the expression of $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$.

###### Proposition 2

For all $k\in\mathbb{Z}$, sequence $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$
satisfies
$\frac{\mathscr{D}^{N}E_{g}\left(t_{k}\right)}{s_{1}\left(2\lambda-h\right)}=\mu_{K_{1}}^{N}\left[k\right]p+\beta_{K_{1}}^{N}\left[k\right]\left(1-p\right),\forall
N\geqslant 2,$ (16) where $k\in\mathbb{Z}$,
$p=\frac{t_{K_{1}+1}-\tau_{1}}{t_{K_{1}+1}-t_{K_{1}}}\in\left[0,1\right]$ and
$\mu_{K_{1}}^{N}\left[k\right]$, $\beta_{K_{1}}^{N}\left[k\right]$ satisfy
$\mu_{K_{1}}^{N+1}\left[k\right]=\frac{\mu_{K_{1}}^{N}\left[{k+1}\right]-\mu_{K_{1}}^{N}\left[{k}\right]}{t_{k+N+1}-t_{k}}\mbox{
and
}\beta_{K_{1}}^{N+1}\left[k\right]=\frac{\beta_{K_{1}}^{N}\left[{k+1}\right]-\beta_{K_{1}}^{N}\left[{k}\right]}{t_{k+N+1}-t_{k}},$
(17) for $N\geqslant 2$, where
$\mu_{K_{1}}^{2}\left[k\right]=\frac{\mathbbmtt{1}_{\left\\{K_{1}-1\right\\}}\left(k\right)}{t_{k+2}-t_{k}}\mbox{
and
}\beta_{K_{1}}^{2}\left[k\right]=\frac{\mathbbmtt{1}_{\left\\{K_{1}\right\\}}\left(k\right)}{t_{k+2}-t_{k}}.$

###### Proof.

We compute $E_{g}\left(t\right)$ by integrating
$\varepsilon_{g}\left(t\right)$ defined in (7)

$E_{g}\left(t\right)=\left(2\lambda-h\right)s_{1}\mathbbmtt{1}_{\left[\tau_{1},\infty\right)}\left(t\right)\cdot\left(t-\tau_{1}\right).$
(18)

We evaluate $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$ for $N=1,2$ as below

$\displaystyle\frac{\mathscr{D}^{1}E_{g}\left(t_{k}\right)}{s_{1}\left(2\lambda-h\right)}$
$\displaystyle=p\cdot\mathbbmtt{1}_{\left\\{K_{1}\right\\}}\left(k\right)+\mathbbmtt{1}_{\left\\{K_{1}+1,\dots,\infty\right\\}}\left(k\right),$
$\displaystyle\frac{\mathscr{D}^{2}E_{g}\left(t_{k}\right)}{s_{1}\left(2\lambda-h\right)}$
$\displaystyle=\frac{p}{t_{K_{1}+1}-t_{K_{1}-1}}\mathbbmtt{1}_{\left\\{K_{1}-1\right\\}}\left(k\right)+\frac{1-p}{t_{K_{1}+2}-t_{K_{1}}}\mathbbmtt{1}_{\left\\{K_{1}\right\\}}\left(k\right)=\mu_{K_{1}}^{2}\left[k\right]p+\beta_{K_{1}}^{2}\left[k\right]\left(1-p\right),$

and the result follows via (14) and (17). ∎

Using (16), the non-zero values of $\mathscr{D}^{N}E_{g}$ therefore satisfy

$\mathrm{supp}\left(\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right)\subseteq\mathbb{S}_{N},$
(19)

where $\mathbb{S}_{N}=\left\\{{K_{1}-N+1},\dots,{K_{1}}\right\\}$. We only
have access to $\mathscr{D}^{N}X\left(t_{k}\right)$, which also includes the
effect of $\mathscr{D}^{N}G\left(t_{k}\right)$. In the first recovery stage,
we compute the support of $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$ via
thresholding. In the case of uniform sampling, a similar strategy was proposed
with a fixed threshold $\Psi^{N}=\frac{\lambda_{h}}{2N}$ [23]. In the present
case of nonuniform sampling, we can achieve input reconstruction in a more
general scenario by allowing a time-varying threshold
$\Psi^{N}\left[k\right]$, which depends on the nonuniform trigger times
$\left\\{t_{k}\right\\}$. In order to constrain the effect of
$\mathscr{D}^{N}G\left(t_{k}\right)$ on the measurements
$\mathscr{D}^{N}X\left(t_{k}\right)$, the following lemma (proof in Section
VI) derives a bound for $\mathscr{D}^{N}G\left(t_{k}\right)$ and subsequently
computes support of $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$.

###### Lemma 2

Assume that

$\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert<\Psi^{N}\left[k\right],$
(20)

where $\Psi^{N}\left[k\right]$ is defined as

$\Psi^{N}\left[k\right]=\lambda_{h}\cdot\min\left\\{\phi_{0,k+N-1},\phi_{0,k+N-2},\phi_{1,k-1},\phi_{1,k}\right\\}$

and sequences $\phi_{i,l},i=0,1,$ $l\in\mathbb{N}$, satisfy

$\displaystyle\begin{split}\phi_{0,l}&=\frac{\left\lvert\mu_{l}^{N}\left[l-N+1\right]\right\rvert\cdot\left\lvert\beta_{l}^{N}\left[l-N+2\right]\right\rvert}{\left\lvert\mu_{l}^{N}\left[l-N+1\right]\right\rvert+\left\lvert\beta_{l}^{N}\left[l-N+2\right]\right\rvert+\left\lvert\mu_{l}^{N}\left[l-N+2\right]\right\rvert},\\\
\phi_{1,l}&=\frac{\left\lvert\mu_{l}^{N}\left[l-1\right]\right\rvert\cdot\left\lvert\beta_{l}^{N}\left[l\right]\right\rvert}{\left\lvert\mu_{l}^{N}\left[l-1\right]\right\rvert+\left\lvert\beta_{l}^{N}\left[l\right]\right\rvert+\left\lvert\beta_{l}^{N}\left[l-1\right]\right\rvert},\end{split}$
(21)

where $\mu_{l}^{N}\left[k\right]$ and $\beta_{l}^{N}\left[k\right]$ are given
in (17). Furthermore, let $\mathcal{K}_{N}=\left\\{k_{m},k_{M}\right\\}$ be
such that

$\displaystyle\begin{split}k_{m}=\min\mathbb{M}_{N},\ \
k_{M}=\max\mathbb{M}_{N},\\\ \mathbb{M}_{N}\triangleq\left\\{k\in\mathbb{Z}\
|\
\left\lvert\mathscr{D}^{N}Z\left(t_{k}\right)\right\rvert\geqslant\Psi^{N}\left[k\right]\right\\}.\end{split}$
(22)

Then $k_{m}\in\left\\{K_{1}-N+1,K_{1}-N+2\right\\}$ and
$k_{M}\in\left\\{K_{1}-1,K_{1}\right\\}$.

The next theorem is our main result, showing that $s_{1}\in\left\\{\pm
1\right\\}$ can be correctly identified and $\tau_{1}$ can be computed with an
accuracy dependent on the sampling rate.

######  Theorem 1 (Estimation of Folding Times)

Let $g\in\mathsf{PW}_{\Omega}$,
$\left\\{t_{k}\right\\}_{k\in\mathbb{Z}}=\mathrm{ASDM}_{\delta,b}\left(\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\right)$
and assume
$\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert<\Psi^{N}\left[k\right],\forall
k\in\mathbb{Z}$. Let $\mathcal{K}_{N}=\\{k_{m},k_{M}\\}$ defined as in (22),
and let $\widetilde{s}_{1}$ be the estimation of the folding sign given by
$\widetilde{s}_{1}=-\mathrm{sign}\left[\mathscr{D}^{N}X\left(t_{k_{m}}\right)\right].$
(23) Let $T_{\mathsf{min}},T_{\mathsf{max}}>0$ be two constants satisfying
$T_{\mathsf{min}}\leqslant t_{k+1}-t_{k}\leqslant T_{\mathsf{max}}$. The
folding time estimation $\widetilde{\tau}_{1}$ is defined as follows. 1. If
$k_{M}-k_{m}=N-1$, let $\widetilde{\tau}_{1}=\frac{t_{k_{M}}+t_{k_{M}+1}}{2}$
and thus
$\widetilde{s}_{1}=s_{1},\quad\left\lvert\tau_{1}-\widetilde{\tau}_{1}\right\rvert\leqslant
T_{\mathsf{max}}/2.$ (24) 2. If $k_{M}-k_{m}=N-2$, let
$\widetilde{\tau}_{1}=t_{k_{M}+1}$. Then
$\widetilde{s}_{1}=s_{1},\quad\left\lvert\tau_{1}-\widetilde{\tau}_{1}\right\rvert\leqslant
T_{\mathsf{max}}.$ (25) 3. If $k_{M}-k_{m}=N-3$, then
$\widetilde{\tau}_{1}=\frac{t_{k_{M}+1}+t_{k_{M}+2}}{2}$ and (24) holds.

Given that $t_{k+1}-t_{k}\leqslant\frac{2\delta}{b-g_{\infty}},\forall
k\in\mathbb{Z}$, by choosing $T_{\mathsf{max}}=\frac{2\delta}{b-g_{\infty}}$,
then $\tau_{1}$ can be recovered with any accuracy for a small enough $\delta$
given by (24,25). We note that in USF the folding times
$\left\\{\tau_{r}\right\\}$ are estimated as $\tau_{r}\approx K_{r}T$, leading
to a similar error as here, bounded by the sampling period $T$.

###  IV-C The Case of Multiple Folds

In the general case there are $R$ folding times
$\left\\{\tau_{r}\right\\}_{r=1}^{R}$, and each has an associated set
$\mathbb{S}_{N}^{r}=\left\\{t_{K_{r}-N+1},\dots,t_{K_{r}}\right\\}$, such that
the non-zero values of $\mathscr{D}^{N}E_{g}\left(t_{k}\right)$ satisfy (26)

$\mathrm{supp}\left(\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right)\subseteq{\bigcup}_{r=1,\dots,R}\mathbb{S}_{N}^{r}.$
(26)

Figure 7: Illustration of the thresholding principle. (a) Each sample from the
filtered data $\mathscr{D}^{N}X\left(t_{k}\right)$ is computed from $N+1$
trigger times. For $N=2$, each folding time $\tau_{r}$ displaces up to two
values of $\mathscr{D}^{N}X\left(t_{k}\right)$ above threshold
$\Psi\left[k\right]$. One can have a minimum of one value of
$\mathscr{D}^{N}X\left(t_{k}\right)$ above the threshold for each folding
time. The minimum is achieved in this example because $\tau_{2}$ is close to
$t_{7}$. (b) The separation of folding times $\tau_{1}$ and $\tau_{2}$ is
enabled by the modulo-hysteresis parameter $h$. If at least $N+1=3$ trigger
times are located between each two folds, the estimation can be performed
sequentially to determine all folding times.

Here the results in Lemma 2 and Theorem 1 can be applied sequentially provided
that $\mathbb{S}_{N}^{r_{1}}\cap\mathbb{S}_{N}^{r_{2}}=\emptyset,\forall
r_{1},r_{2}=1\dots,R,r_{1}\neq r_{2}$. In other words, this means there are at
least $N$ trigger times between each two consecutive folds
$\tau_{r},\tau_{r+1}$, which is guaranteed if $\tau_{r+1}-\tau_{r}\geqslant
NT_{\mathsf{max}}$. A sufficient condition can be obtained via Lemma 1

$\frac{h^{*}}{\Omega g_{\infty}}\geqslant NT_{\mathsf{max}}.$ (27)

Then if (20,27) are true then $\left\\{\tau_{r}\right\\}_{r=1}^{R}$ and
$\left\\{s_{r}\right\\}_{r=1}^{R}$ can be computed as follows. Let
$\mathcal{K}_{N}=\bigcup\nolimits_{r\in\left\\{1,R\right\\}}{\left\\{k_{m}^{r},k_{M}^{r}\right\\}}$,
where $k_{m}^{r},k_{M}^{r}$ are computed as below. For $r=1$,

$\displaystyle\begin{split}k_{m}^{1}&=\min\left\\{k\ |\ \left\lvert
X_{k}^{N}\right\rvert\geqslant\Psi^{N}\left[k\right]\right\\},\\\
k_{M}^{1}&=\max\left\\{k\leqslant k_{m}^{1}+N-1\ |\ \left\lvert
X_{k}^{N}\right\rvert\geqslant\Psi^{N}\left[k\right]\right\\}.\end{split}$
(28)

For $r=2,\dots,R$

$\displaystyle\begin{split}k_{m}^{r}&=\min\left\\{k>k_{M}^{r-1}\ |\
\left\lvert X_{k}^{N}\right\rvert\geqslant\Psi^{N}\left[k\right]\right\\},\\\
k_{M}^{r}&=\max\left\\{k\leqslant k_{m}^{r}+N-1\ |\ \left\lvert
X_{k}^{N}\right\rvert\geqslant\Psi^{N}\left[k\right]\right\\},\end{split}$
(29)

where $X_{k}^{N}\triangleq\mathscr{D}^{N}X\left(t_{k}\right)$. Then
$\widetilde{s}_{r},\widetilde{\tau}_{r},r=1,\dots,R$ can be computed
sequentially with Theorem 1 from $\mathcal{K}_{N}$. To better explain the
proposed technique, we illustrated the thresholding strategy for multiple
folds in Fig. 7. The diagram shows the role of $h$ in separating the folds
which, in turn, enables an independent detection of the folding times. The
following lemma provides sufficient recovery conditions.

######  Lemma 3 (Sufficient Recovery Conditions)

The recovery is possible when either of the conditions 1) and 2) below is
true. 1. The values $\widetilde{s}_{r},\widetilde{\tau}_{r},r=1,\dots,R$ can
be recovered from $\mathscr{D}^{N}X\left(t_{k}\right)$ with an error given by
(24,25) if $\mathrm{S}_{1}$) $\left(\frac{C^{2}}{b-\lambda}2\delta\Omega
e\right)^{N-1}g_{\infty}<\frac{\lambda_{h}}{CNe}$. $\mathrm{S}_{2}$)
$N\frac{2\delta\Omega g_{\infty}}{b-\lambda}\leqslant h^{*}$, where
$h^{\ast}=\min\left\\{h,2\lambda-h\right\\}$ and
$C=\frac{b+\lambda}{b-\lambda}$. 2. Furthermore, $\mathrm{S}_{1})$ and
$\mathrm{S}_{2})$ are satisfied if
$\delta<\frac{\left(b-\lambda\right)h^{*}\kappa}{2N\Omega g_{\infty}}\mbox{
where
}\kappa=\min\left\\{{1,\frac{{{\lambda_{h}}}}{e^{2}h^{*}C^{2+\frac{1}{N-1}}}}\right\\}.$
(30)

###### Proof.

  1. 1.

We use the upper bound for
$\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert$ in Proposition 1
to derive a sufficient condition for (20) as

$\left(\frac{T_{\mathsf{max}}^{3}}{T_{\mathsf{min}}^{2}}\Omega
e\right)^{N}g_{\infty}<\frac{\lambda_{h}}{Ne}\frac{T_{\mathsf{min}}}{T_{\mathsf{max}}}.$
(31)

To prove (31), note that, in the limit case,
$T_{\mathsf{max}}=T_{\mathsf{min}}=T$ is always satisfied for large enough $N$
if

$\frac{T_{\mathsf{max}}^{3}}{T_{\mathsf{min}}^{2}}\Omega e=T\Omega e<1.$

The latter is the sampling rate condition for uniform sampling [14]. Then
$\left.{\mathrm{S}}_{1}\right)$ follows from (31) and
$\left.{\mathrm{S}}_{2}\right)$ from (27) where $T_{\mathsf{min}}$ and
$T_{\mathsf{max}}$ are selected as the following bounds for the intervals
between the trigger times [1]

$T_{\mathsf{min}}=\frac{2\delta}{b+\lambda}\leqslant\left\lvert
t_{k+1}-t_{k}\right\rvert\leqslant\frac{2\delta}{b-\lambda}=T_{\mathsf{max}}.$
(32)

  2. 2.

From $\left.{\mathrm{S}}_{1}\right)$ and $\left.{\mathrm{S}}_{2}\right)$ it
follows that
$2\delta<\min\left\\{B_{1}\left(N\right),B_{2}\left(N\right)\right\\}$, where

$B_{1}\left(N\right)=\frac{b-\lambda}{C^{2}\Omega
e}\left(\frac{\lambda_{h}}{g_{\infty}CNe}\right)^{\frac{1}{N-1}}\ \text{and}\
B_{2}\left(N\right)=\frac{\left({b-\lambda}\right)h^{*}}{N\Omega g_{\infty}}.$

A direct calculation yields
$B_{1}\left(N\right)\frac{e^{2}h^{*}C^{2}\sqrt[N-1]{C}}{\lambda_{h}}\geqslant
B_{2}\left(N\right)$, which completes the proof.

∎

If $g_{\infty}$ is unknown yet an upper bound is known, the results in Lemma 3
still hold true by replacing $g_{\infty}$ with the upper bound. Furthermore,
we note that the conditions in Lemma 3 are only sufficient, and can be
conservative. Therefore, in practice, one can achieve good reconstructions
with higher values of $N$, some of which that might not satisfy the lemma.

###  IV-D Input Reconstruction

Assuming $\left\\{\widetilde{\tau}_{r},\widetilde{s}_{r}\right\\}$ known, the
estimation of $E_{g}\left(t_{k}\right)$ is

$\widetilde{E}_{g}\left(t_{k}\right)=2\lambda_{h}\sum\limits_{r=1}^{R}\widetilde{s}_{r}\mathbbmtt{1}_{\left[\widetilde{\tau}_{r},\infty\right)}\left(t_{k}\right)\cdot\left(t_{k}-\widetilde{\tau}_{r}\right).$
(33)

Finally, $\widetilde{g}\left(t\right)$ is recovered recursively using
reconstruction from local averages,

$\Delta\widetilde{G}\left(t_{k}\right)=\widetilde{G}\left(t_{k+1}\right)-\widetilde{G}\left(t_{k}\right)\stackrel{{\scriptstyle(\ref{eq:rec_local_av})}}{{=}}\int_{t_{k}}^{t_{k+1}}g\left(s\right)ds,$

where
$\widetilde{G}\left(t_{k}\right)=X\left(t_{k}\right)+\widetilde{E}_{g}\left(t_{k}\right)$.
This recovery approach was proposed in the context of irregular sampling in
[13] and adapted for event-driven ASDM sampling in [1]. In our case, sequence
$\widetilde{g}_{n}\left(t\right)$ is computed as

$\displaystyle\begin{split}\widetilde{g}_{0}\left(t\right)&=\mathcal{S}_{\Omega}\left[{\Delta\widetilde{G}\left(t_{k}\right)}\right]\left(t\right),\\\
\widetilde{g}_{n+1}\left(t\right)&=\widetilde{g}_{n}\left(t\right)+\widetilde{g}_{0}\left(t\right)-\mathcal{S}_{\Omega}\left[{\mathcal{L}\widetilde{g}_{n}}\right]\left(t\right),\end{split}$
(34)

where,

  * •

$\mathcal{S}_{\Omega}:\ell^{2}\rightarrow L^{2}\left(\mathbb{R}\right)$.

  * •

$\mathcal{S}_{\Omega}\left[{a_{k}}\right]\triangleq\sum\limits_{k\in\mathbb{Z}}a_{k}\cdot\mathrm{sinc}_{\Omega}\left(t-s_{k}\right)$,
$s_{k}=\frac{t_{k}+t_{k+1}}{2}$, and,

  * •

$\mathcal{L}:L^{2}\left(\mathbb{R}\right)\rightarrow\ell^{2},\left(\mathcal{L}f\right)_{k}=\int\limits_{t_{k}}^{t_{k+1}}f\left(s\right)ds$.

We outline the recovery procedure in Algorithm 1.

Data: $\lambda,h,\delta,b\in\mathbb{R}_{+},\\{t_{k}\\},\Omega>0,N\geqslant 2$.

Result: $\widetilde{g}\left(t\right)$

  1. 1.

Compute
$X_{k}^{N}=\mathscr{D}^{N}\left(\sum\limits_{m={0}}^{k-1}q_{m}\right)$, where
$q_{m}=\left(-1\right)^{m}2\delta-\left(-1\right)^{m}b\left(t_{m+1}-t_{m}\right)$.

  2. 2.

Compute $R,\ \left\\{k_{m}^{r},k_{M}^{r}\right\\}_{r=1}^{R}$ with (28), (29).

  3. 3.

For $r=1,\dots,R$, compute $\widetilde{s}_{r}$ with (23) and,

$\widetilde{\tau}_{r}=\frac{1}{2}\begin{cases}t_{k_{M}^{r}}+t_{k_{M}^{r}+1}&\mbox{
if }\left(k_{M}^{r}-k_{m}^{r}\right)=N-1\\\
2\left(t_{k_{M}^{r}+1}\right)&\mbox{ if
}\left(k_{M}^{r}-k_{m}^{r}\right)=N-2\\\
t_{k_{M}^{r}+1}+t_{k_{M}^{r}+2}&\mbox{ if
}\left(k_{M}^{r}-k_{m}^{r}\right)=N-3\end{cases}.$

  4. 4.

Compute $\widetilde{E}_{g}\left(t_{k}\right)$ using (33) and
$\widetilde{G}\left(t_{k}\right)$ using (12).

  5. 5.

Compute $\widetilde{g}_{n}\left(t\right)$ from samples
$\Delta\widetilde{G}\left(t_{k}\right)$ with (34).

Algorithm 1 Recovery Algorithm.

The following result proven in Section VI gives a bound for the input
reconstruction error.

###### Proposition 3

Assuming that the conditions in Lemma 3 are satisfied and that the data
consists of $R$ folding times, then the input can be recovered recursively as
$\widetilde{g}_{n}\left(t\right)$, where $n$ is the iteration number. The
reconstruction error satisfies:
$\|g-\widetilde{g}_{n}\|_{L^{2}}\leqslant\frac{4\lambda_{h}2\delta
R\sqrt{\Omega\pi}}{\pi\left(b-\lambda\right)-2\delta\Omega}+\left(\frac{2\delta\Omega}{\pi\left(b-\lambda\right)}\right)^{n+1}\|g\|_{L^{2}}.$
(35)

We note that the second term in the error bound (35) is guaranteed to vanish
for $n\rightarrow\infty$ if
$\frac{2\delta\Omega}{\pi\left(b-\lambda\right)}<1$, which is satisfied if
condition $\mathrm{S}_{2})$ in Lemma 3 is true. Furthermore, the error in (35)
can be made arbitrarily small for a small enough ASDM threshold $\delta$.
Furthermore, we note that Algorithm 1 has an intrinsic method to indicate that
the data does not verify the recovery conditions, if
$k_{M}^{r}-k_{m}^{r}\not\in\left\\{N-3,N-2,N-1\right\\}$ for any
$r\in\left\\{1,\dots,R\right\\}$. Although theoretically there could be
situations where the algorithm runs with no error for an input not satisfying
the recovery conditions, we found this to be very unlikely in practice.

##  V Numerical and Hardware Experiments

We validate the proposed event-driven model using both synthetic and real data
measurements. We generate a bandlimited function $g\in\mathsf{PW}_{\Omega}$,
which is encoded directly using the standalone ASDM into trigger times
$t_{k}^{\mathsf{ASDM}}$. The same input is then processed with the proposed
encoding model, which triggers times $t_{k}^{\mathsf{MEDS}}$. The direct
reconstruction from ASDM samples is denoted
$\widetilde{g}_{\mathsf{ASDM}}\left(t\right)$, and the proposed reconstruction
from the MEDS samples is denoted
$\widetilde{g}_{\mathsf{MEDS}}\left(t\right)$. We measure the relative
recovery error via the ratio $\mathsf{Err}\left(\widetilde{g},g\right)$,
defined as

$\mathsf{Err}\left(\widetilde{g},g\right)=100\cdot\frac{\|g-\widetilde{g}\|_{L^{2}}}{\|g\|_{L^{2}}}\left(\%\right).$
(36)

We denote the recovery errors with each method by
$\mathsf{Err}_{\mathsf{ASDM}}=\mathsf{Err}\left(\widetilde{g}_{\mathsf{ASDM}},g\right)$
and
$\mathsf{Err}_{\mathsf{MEDS}}=\mathsf{Err}\left(\widetilde{g}_{\mathsf{MEDS}},g\right)$,
respectively.

The recovery results from synthetic data are in V-A, and the ones from
experimental data are presented in V-B.

###  V-A Reconstruction using Synthetic Data

Here we compare the proposed model with the ASDM encoder. The input is
$g\left(t\right)=Ag_{0}\left(t\right)$, where

$g_{0}\left(t\right)=\sum\nolimits_{n\in\mathbb{Z}}c_{n}\mathrm{sinc}_{\Omega}\left(t-n\pi/\Omega\right),$

where $\Omega=150\ \mathrm{rad/s}$ and $c_{n}$ are drawn from the uniform
distribution $U\left(\left[-1,1\right]\right)$. The signal is truncated such
that $t\in\left[0,0.13\ \mathrm{s}\right]$, and normalised such that
$\mathrm{max}_{t}\left\lvert g_{0}\left(t\right)\right\rvert=1$. The ASDM
parameters are $\delta=2.5\times 10^{-3},b=9$. The Nyquist rate condition is
$\left\lvert g\left(t\right)\right\rvert<g_{\textsf{MAX}}$, where
$g_{\textsf{MAX}}$ denotes the ASDM dynamic range $g_{\textsf{MAX}}=\frac{\pi
b-2\delta\Omega}{\pi}=8.76$.

The proposed encoder consists of a modulo model with threshold
$\lambda=g_{\textsf{MAX}}/2$ and hysteresis parameter $h=\lambda/2$ in series
with an ASDM whose parameters are the same as the standalone ASDM encoder.
Therefore the models are equivalent for $\left\lvert
g\left(t\right)\right\rvert<\lambda$. For $A=34.6$, the ASDM generates trigger
times $\left\\{t_{k}^{\textsf{ASDM}}\right\\}_{k=1}^{14}$ in response to input
$g\left(t\right)$ when its amplitude is within the ASDM dynamic range
$\left[-8.76,8.76\right]$, and sends no information when the input amplitude
is outside the range. However, the proposed model generates trigger times
$\left\\{t_{k}^{\textsf{MEDS}}\right\\}_{k=1}^{220}$, which cover the whole
duration of $g\left(t\right)$, as depicted in Fig. 8 (a). The input is first
recovered from the ASDM samples
$\left\\{t_{k}^{\textsf{ASDM}}\right\\}_{k=1}^{14}$ [1], which leads to a high
reconstruction error especially where the input has a large amplitude, as
depicted in Fig. 8(d). For the proposed model, we used the $\mathscr{D}^{3}$
to identify correctly the $21$ folding times. The filtered data
$\mathscr{D}^{3}X\left(t_{k}^{\textsf{MEDS}}\right)$, the threshold
$\Psi_{N}\left[k\right]$ and the estimated folding times are in Fig. 8 (b,c).
We then recovered $\widetilde{g}_{\mathsf{MEDS}}\left(t\right)$ using
Algorithm 1. The resulted errors for the two methods are
$\mathsf{Err}_{\textsf{ASDM}}=6.11\cdot 10^{3}\%$ and
$\mathsf{Err}_{\textsf{MEDS}}=0.25\%$.

To see the effect of changing the ASDM threshold $\delta$ on the recovery
performance, we repeated the recovery for $10$ uniformly spaced values of
$\delta$ in interval $\left[10^{-3},3\cdot 10^{-3}\right]$. We evaluated the
input error as before, the number of spikes generated, and the folding time
error evaluated as
$\mathsf{Err}_{\tau}=100\cdot\tfrac{\|\tau-\widetilde{\tau}\|_{\ell^{2}}}{\|\tau\|_{\ell^{2}}}\left(\%\right)$.
The results, depicted in Fig. 9, show that when $\delta>2.6\cdot 10^{-3}$ the
recovery algorithm does not identify the folds correctly leading to a large
error. Choosing $\delta<1.8$ leads to an exponential increase in sample size
without a significant decrease in error. Moreover, we remark that, in this
example, the input error is around $6$ times larger than the folding time
error. This is in line with the theoretical derivation that the error for
computing $\Delta G\left(t_{k}\right)$ is $2\lambda_{h}=6.57$ times larger
than the folding time error as shown in (57).

|
---|---
|
Figure 8: Simulations with synthetic data. (a) Input, output of the modulo
nonlinearity and corresponding output trigger times. (b) Estimating the
folding times by thresholding the filtered data. (c) Reconstruction of the
residual function. (d) Input, proposed reconstruction and conventional IF
based recovery.

Figure 9: The effect of the ASDM threshold on recovery accuracy and number of
samples.

###  V-B Reconstruction using Experimental Data

In this example, the input is
$g\left(t\right)=4.51\cdot\sin\left(\Omega\left(t-\tau_{0}\right)\right)$
where $\Omega=125\ \mathrm{rad/s}$ and $\tau_{0}=1.6\cdot 10^{-2}$. The output
data was generated using a MEDS hardware prototype shown in Fig. 5 that is
coarsely calibrated such that the modulo output is within the ASDM dynamic
range. The corresponding acquisition pipeline is shown in Fig. 10. The MEDS
parameters, estimated using line search based optimization [39], are
$\lambda=1.53$, $h=1.51$, $\delta=2.07\times 10^{-4}$, and $b=2.22$. Therefore
the ASDM dynamic range is $g_{\mathsf{MAX}}=\frac{\pi
b-2\delta\Omega}{\pi}=2.206$. The output of the prototype in response to
$g\left(t\right)$ is $\left\\{t_{k}^{\textsf{MEDS}}\right\\}_{k=1}^{214}$
depicted in Fig. 11 (a).

Figure 10: The diagram of the hardware setup used for experiments. The input
$g\left(t\right)$ is generated via a digital-to-analog converter (DAC) from
digital samples. The MEDS hardware prototype is implemented using an US-ADC in
series with an ASDM circuit. The continuous MEDS output $z\left(t\right)$ and
input $g\left(t\right)$ are then sampled with an oscilloscope to generate
uniform samples $z\left(kT\right)$ and $g\left(kT\right)$. The ASDM trigger
times $\left\\{t_{k}\right\\}$ are subsequently computed from the zero-
crossings of $z\left(kT\right)$.

Figure 11: Input reconstruction from measurements based on our MEDS hardware
(Fig. 5). (a) Estimating the folding times by thresholding the filtered data.
(b) Reconstruction of the residual function. (c) The proposed reconstruction
superimposed on the input data.

We used the $\mathscr{D}^{2}$ to identify correctly the $12$ folding times.
The filtered data $\mathscr{D}^{2}X\left(t_{k}^{\textsf{MEDS}}\right)$, the
threshold $\Psi_{N}\left[k\right]$ and the estimated folding times are in Fig.
11 (b,c). The input is recovered with Algorithm 1. The resulted error is
$\mathsf{Err}_{\textsf{MEDS}}=0.68\%$.

##  VI Proofs

Here we prove Proposition 1, Lemma 2, Theorem 1, and Proposition 3.

###### Proof of Proposition 1.

We fix an arbitrary value $k\in\mathbb{Z}$ and expand $G(t_{k+l}),\quad l\in$
$\left\\{0,\dots,N\right\\}$ in Taylor series around
$\tau=\frac{t_{k+N}+t_{k}}{2}$

$\displaystyle G\left(t_{k+l}\right)=G(\tau)+$
$\displaystyle\sum\limits_{n=1}^{N-1}\left(t_{k+l}-\tau\right)^{n}g^{\left(n-1\right)}\left(\tau\right)+\frac{\left(t_{k+l}-\tau\right)^{N}}{N!}g^{\left(N-1\right)}\left(\zeta_{k+l}\right),$
(37)

where $\zeta_{k+l}$ is between $t_{k+l}$ and $\tau$. The equation above uses
that $G^{\left(1\right)}(t)=g(t).$ Let $r\left(k+l\right)$ denote the Taylor
series remainder
$r\left(k+l\right)=\frac{\left(t_{k+l}-\tau\right)^{N}}{N!}g^{\left(N-1\right)}\left(\zeta_{k+l}\right)$.
We use the fact that the operator $\mathscr{D}^{N}$ cancels out polynomials of
degree up to $N-1$, and thus

$\displaystyle\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert=\left\lvert\mathscr{D}^{N}r\left(k\right)\right\rvert$
$\displaystyle=\left\lvert\frac{\mathscr{D}^{N-1}r\left(k+1\right)-\mathscr{D}^{N-1}r\left(k\right)}{t_{k+N}-t_{k}}\right\rvert$
$\displaystyle\leqslant\frac{\left\lvert\mathscr{D}^{N-1}r\left(k+1\right)\right\rvert+\left\lvert\mathscr{D}^{N-1}r\left(k\right)\right\rvert}{t_{k+N}-t_{k}}$
$\displaystyle\leqslant\frac{2\cdot\underset{l=0,1}{\max}\left\lvert\mathscr{D}^{N-1}r\left(k+l\right)\right\rvert}{NT_{\mathsf{min}}}.$

Similarly, via induction, it can be shown that

$\displaystyle\begin{split}\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert&\leqslant\frac{2^{N}}{T_{\mathsf{min}}^{N}\cdot
N!}\cdot\underset{l=0,\dots,N}{\max}\left\lvert\mathscr{D}^{0}r\left(k+l\right)\right\rvert\\\
&=\frac{2^{N}}{T_{\mathsf{min}}^{N}\cdot
N!}\cdot\underset{l=0,\dots,N}{\max}\left\lvert\frac{\left(t_{k+l}-\tau\right)^{N}}{N!}g^{\left(N-1\right)}\left(\zeta_{k+l}\right)\right\rvert\\\
&=\frac{2^{N}}{T_{\mathsf{min}}^{N}\cdot
N!}\cdot\|g^{\left(N-1\right)}\|_{\infty}\underset{l=0,\dots,N}{\max}\left\lvert\frac{\left(t_{k+l}-\tau\right)^{N}}{N!}\right\rvert\\\
&\leqslant\frac{2^{N}}{T_{\mathsf{min}}^{N}\cdot
N!}\cdot\frac{\|g^{\left(N-1\right)}\|_{\infty}}{N!}\cdot\left(\frac{t_{k+N}-t_{k}}{2}\right)^{N}\\\
&\leqslant\left(\frac{N\cdot
T_{\mathsf{max}}}{T_{\mathsf{min}}}\right)^{N}\cdot\frac{\|g^{\left(N-1\right)}\|_{\infty}}{\left(N!\right)^{2}}.\end{split}$

The following two inequalities follow using the Stirling and Bernšteĭn (cf.
pg. 116, [43]) inequalities, respectively.

$\displaystyle\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert$
$\displaystyle\leqslant\left(\frac{e\cdot
T_{\mathsf{max}}}{T_{\mathsf{min}}}\right)^{N}\cdot\frac{\|g^{\left(N-1\right)}\|_{\infty}}{N!}$
$\displaystyle\leqslant\frac{1}{\Omega
N!}\left(\frac{T_{\mathsf{max}}}{T_{\mathsf{min}}}\Omega
e\right)^{N}\|g\|_{\infty}.$

∎

To prove Lemma 2, we first provide a few properties of
$\mu_{l}^{N}\left[k\right]$ and $\beta_{l}^{N}\left[k\right]$ in the following
proposition.

######  Proposition 4 (Properties of $\mu_{l}^{N}$ and $\beta_{l}^{N}$)

The following hold, $\mathrm{supp}\
\mu_{l}^{N}\left[k\right]=\left\\{l-N+1,\dots,l-1\right\\}\quad\mbox{and}\quad\mathrm{supp}\
\beta_{l}^{N}\left[k\right]=\left\\{l-N+2,\dots,l\right\\}.$ (38) For all
$k\in\mathbb{Z},\forall N\geqslant 2$ it follows that
$\displaystyle\begin{split}\mathrm{sign}\left(\mu_{l}^{N}\left[k\right]\right)\cdot\mathrm{sign}\left(\mu_{l}^{N}\left[{k+1}\right]\right)&\leqslant
0,\\\
\mathrm{sign}\left(\beta_{l}^{N}\left[k\right]\right)\cdot\mathrm{sign}\left(\beta_{l}^{N}\left[{k+1}\right]\right)&\leqslant
0,\\\
\mathrm{sign}\left(\mu_{l}^{N}\left[k\right]\right)\cdot\mathrm{sign}\left(\beta_{l}^{N}\left[{k}\right]\right)&\leqslant
0.\end{split}$ (39) The following bounds hold true.
$\displaystyle\begin{split}\frac{1}{N!\hskip
1.00006ptT_{\mathsf{max}}^{N-1}}&\leqslant\left\lvert\beta_{l}^{N}\left[k\right]\right\rvert\leqslant\frac{1}{N!\hskip
1.00006ptT_{\mathsf{min}}^{N-1}},\ k\in\left\\{l-N+2,l\right\\},\\\
\frac{1}{N!\hskip
1.00006ptT_{\mathsf{max}}^{N-1}}&\leqslant\left\lvert\mu_{l}^{N}\left[k\right]\right\rvert\leqslant\frac{1}{N!\hskip
1.00006ptT_{\mathsf{min}}^{N-1}},\ k\in\left\\{l-N+1,l-1\right\\},\end{split}$
(40)
$\left\lvert\mu_{l}^{N}\left[{l-N+2}\right]\right\rvert\leqslant\frac{N-2}{N!\
T_{\mathsf{min}}^{N-1}},$ (41)
$\left\lvert\beta_{l}^{N}\left[{l-1}\right]\right\rvert\leqslant\frac{N-2}{N!\
T_{\mathsf{min}}^{N-1}}.$ (42)

###### Proof.

Equation (38) is derived directly from the definitions of $\mu_{l}^{N}$ and
$\beta_{l}^{N}$ in Proposition 2. Moreover, (39) can be shown directly via
induction for $N\geqslant 2$. Bounds (40) are shown recursively by using
$nT_{\mathsf{min}}\leqslant t_{k+n}-t_{k}\leqslant nT_{\mathsf{max}}$.

For (41), we compute recursively an upper bound for
$m_{N}\triangleq\left\lvert\mu_{l}^{N}\left[{l-N+2}\right]\right\rvert$, where
$m_{2}=0$ (38) and, for $N>2$,

$\displaystyle
m_{N+1}=\frac{\left\lvert\mu_{l}^{N}\left[{l-N+2}\right]-\mu_{l}^{N}\left[{l-N+3}\right]\right\rvert}{t_{l+2}-t_{l-N+1}}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:signs_mu_beta})}}{{=}}\frac{m_{N}+\left\lvert\mu_{l}^{N}\left[{l-N+1}\right]\right\rvert}{t_{l+2}-t_{l-N+1}}$
$\displaystyle\leqslant\frac{m_{N}+\left\lvert\mu_{l}^{N}\left[{l-N+1}\right]\right\rvert}{\left(N+1\right)T_{\mathsf{min}}}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:edge_bound_mu_beta})}}{{=}}{\leqslant}\frac{m_{N}}{\left(N+1\right)T_{\mathsf{min}}}+\frac{1}{\left(N+1\right)!\hskip
1.00006ptT_{\mathsf{min}}^{N}}.$

Via induction, we derive (41). The same procedure is used for the derivation
of (42). ∎

###### Proof of Lemma 2.

Due to (19, 20) it directly follows that

$\left\lvert\mathscr{D}^{N}X\left(t_{k}\right)\right\rvert=\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)-\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right\rvert<\Psi^{N}\left[k\right],$
(43)

for all
$k\notin\mathbb{S}_{N}\Rightarrow\mathcal{K}_{N}\subseteq\mathbb{M}_{N}\subseteq\mathbb{S}_{N}$.
Furthermore, we note that

$\displaystyle\begin{split}\left\lvert\mathscr{D}^{N}X\left(t_{k}\right)\right\rvert&\geqslant\left\lvert\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right\rvert-\left\lvert\mathscr{D}^{N}G\left(t_{k}\right)\right\rvert\\\
&>\left\lvert\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right\rvert-\Psi^{N}\left[k\right],\forall
k\in{\mathbb{Z}_{+}}.\end{split}$ (44)

Specifically, we will identify the values of $k$ for which
$\left\lvert\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right\rvert\geqslant
2\Psi^{N}\left[k\right]$, which are contained in $\mathbb{M}_{N}$ according to
(44). This yields upper bounds for $k_{m}$ and lower bounds for $k_{M}$. With

$E_{k}^{N}\triangleq\frac{\left\lvert\mathscr{D}^{N}E_{g}\left(t_{k}\right)\right\rvert}{2\lambda_{h}},$

the values $k$ of interest satisfy

$E_{k}^{N}\geqslant\frac{\Psi^{N}\left[k\right]}{\lambda_{h}}=\min\left\\{\phi_{0,k+N-1},\phi_{0,k+N-2},\phi_{1,k-1},\phi_{1,k}\right\\}.$
(45)

An important step in the proof is to show that

$\displaystyle\max\left\\{E_{l-N+1}^{N},E_{K_{1}-N+2}^{N}\right\\}\geqslant\phi_{0,K_{1}},$
(46)
$\displaystyle\max\left\\{E_{K_{1}-1}^{N},E_{K_{1}}^{N}\right\\}\geqslant\phi_{1,K_{1}}.$
(47)

To show (46), we denote by $f_{1},f_{2}:\left[0,1\right]\rightarrow\mathbb{R}$
two functions

$\displaystyle f_{1}\left(x\right)$
$\displaystyle=\left\lvert\bar{\mu}_{1}x+\bar{\beta}_{1}\left(1-x\right)\right\rvert\overset{\eqref{eq:supp_mu_beta}}{=}\left\lvert\bar{\mu}_{1}\right\rvert
x,$ $\displaystyle f_{2}\left(x\right)$
$\displaystyle=\left\lvert\bar{\mu}_{2}x+\bar{\beta}_{2}\left(1-x\right)\right\rvert\overset{\eqref{eq:signs_mu_beta}}{=}\left\lvert\left\lvert\bar{\mu}_{2}\right\rvert
x-\left\lvert\bar{\beta}_{2}\right\rvert\left(1-x\right)\right\rvert,$

where
$\bar{\mu}_{i}=\mu_{K_{1}}^{N}\left[K_{1}-N+i\right],\bar{\beta}_{i}=\beta_{K_{1}}^{N}\left[K_{1}-N+i\right],i=1,2$.
Then we have that $E_{K_{1}-N+1}^{N}=f_{1}\left(p\right)$ and
$E_{K_{1}-N+2}^{N}=f_{2}\left(p\right)$ (16), which gives

$\max\left\\{E_{K_{1}-N+1}^{N},E_{K_{1}-N+2}^{N}\right\\}\geqslant\min_{x\in\left[0,1\right]}\left(\max_{i=1,2}f_{i}\left(x\right)\right).$
(48)

Function $f_{1}\left(x\right)$ is strictly increasing, whereas
$f_{2}\left(x\right)$ has one zero $f\left(x^{\ast}\right)=0$, where

$x^{\ast}=\frac{\left\lvert\bar{\beta}_{2}\right\rvert}{\left\lvert\bar{\mu}_{2}\right\rvert+\left\lvert\bar{\beta}_{2}\right\rvert}.$

Then $f_{2}\left(x\right)$ is strictly decreasing for
$x\in\left[0,x^{\ast}\right]$, and thus, given that $f_{1}(0)=0$ and
$f_{2}(0)>0$, it follows that there exists a unique point

$x^{\ast\ast}\in\left(0,x^{\ast}\right)$

such that $f_{1}\left(x^{\ast\ast}\right)=f_{2}\left(x^{\ast\ast}\right)$.
Indeed, when solving $f_{1}\left(x\right)=f_{2}\left(x\right)$ for $x$, there
are two solutions

$x_{k}=\frac{\left\lvert\bar{\beta}_{2}\right\rvert}{\left\lvert\bar{\beta}_{2}\right\rvert+\left\lvert\bar{\mu}_{2}\right\rvert+\left(-1\right)^{k}\left\lvert\bar{\mu}_{1}\right\rvert},\quad
k=1,2,$

such that $x_{2}<x^{\ast}<x_{1}$, and thus

$x^{\ast\ast}=\frac{\left\lvert\bar{\beta}_{2}\right\rvert}{\left\lvert\bar{\beta}_{2}\right\rvert+\left\lvert\bar{\mu}_{2}\right\rvert+\left\lvert\bar{\mu}_{1}\right\rvert}.$

It follows that
$\max_{i=1,2}f_{i}\left(x\right)=f_{2}\left(x\right),x\in\left[0,x^{\ast\ast}\right]$
and
$\max_{i=1,2}f_{i}\left(x\right)=f_{1}\left(x\right),x\in\left[x^{\ast\ast},x^{\ast}\right]$.
Given that for $x>x^{\ast}$ $f_{2}\left(x\right)$ is strictly increasing, then
$f_{2}\left(x\right)>f_{2}\left(x^{\ast}\right),\forall x>x^{\ast}$, and thus

$\max_{i=1,2}f_{i}\left(x\right)\geqslant\max_{\underset{x\in\left[0,x^{\ast}\right]}{i=1,2}}f_{i}\left(x\right)\geqslant
f_{1}\left(x^{\ast\ast}\right)={\left\lvert\bar{\mu}_{1}\right\rvert}x^{\ast\ast}=\phi_{0,K_{1}},$

which yields (46) via (48).

The bound (47) is proven following the same steps as for (46). Specifically,
we denote by $f_{3},f_{4}:\left[0,1\right]\rightarrow\mathbb{R}$

$\displaystyle f_{3}\left(x\right)$
$\displaystyle=\left\lvert\bar{\mu}_{3}x+\bar{\beta}_{3}\left(1-x\right)\right\rvert\overset{\eqref{eq:signs_mu_beta}}{=}\left\lvert\left\lvert\bar{\mu}_{3}\right\rvert
x-\left\lvert\bar{\beta}_{3}\right\rvert\left(1-x\right)\right\rvert,$
$\displaystyle f_{4}\left(x\right)$
$\displaystyle=\left\lvert\bar{\mu}_{4}x+\bar{\beta}_{4}\left(1-x\right)\right\rvert\overset{\eqref{eq:supp_mu_beta}}{=}\left\lvert\bar{\beta}_{4}\right\rvert\left(1-x\right),$

where

$\bar{\mu}_{i}=\mu_{K_{1}}^{N}\left[K_{1}-4+i\right],\bar{\beta}_{i}=\beta_{K_{1}}^{N}\left[K_{1}-4+i\right],i=3,4.$

Using $y=1-x$ we get the same problem as before, which yields (47) via direct
calculation. Finally, (46,47) together with (44) imply that $k_{m}\leqslant
K_{1}-N+2$ and $k_{M}\geqslant K_{1}-1$. This completes the proof of the lemma
via $\mathcal{K}_{N}\subseteq\mathbb{S}_{N}$. ∎

###### Proof of Theorem 1.

To show that $\widetilde{s}_{1}=s_{1}$, we observe that
$\mathscr{D}^{N}G\left(t_{k_{m}}\right)<\Psi^{N}\left[k_{m}\right]$ and
$\mathscr{D}^{N}X\left(t_{k_{m}}\right)\geqslant\Psi^{N}\left[k_{m}\right]$
which implies that
$\mathrm{sign}\left[\mathscr{D}^{N}X\left(t_{k_{m}}\right)\right]=-\mathrm{sign}\left[\mathscr{D}^{N}E_{g}\left(t_{k_{m}}\right)\right]$
and thus (16)

$\widetilde{s}_{1}=\mathrm{sign}\left[\mu_{K_{1}}^{N}\left[{k_{m}}\right]p+\beta_{K_{1}}^{N}\left[{k_{m}}\right]\left(1-p\right)\right]\cdot
s_{1}.$ (49)

  1. 1.

If $k_{M}-k_{m}=N-1$, then, according to Lemma 2, we know that
$k_{m}=K_{1}-N+1,k_{M}=K_{1}$. Then, due to (38),

$\mu_{K_{1}}^{N}\left[{k_{m}}\right]p+\beta_{K_{1}}^{N}\left[{k_{m}}\right]\left(1-p\right)=\mu_{K_{1}}^{N}\left[{k_{m}}\right]p.$

We know that
$\mu_{K_{1}}^{2}\left[{K_{1}-1}\right]=\frac{1}{t_{K_{1}+1}-t_{K_{1}-1}}>0$.
It can be shown via induction that
$\mathrm{sign}\left(\mu_{K_{1}}^{N}\left[{K_{1}-N+1}\right]\right)=1,\forall
N\geqslant 2$, and thus $\widetilde{s}_{1}=s_{1}$ via (49). Moreover,
$\widetilde{\tau}_{1}=\frac{t_{K_{1}}+t_{K_{1}+1}}{2}$ which satisfies (24)
given that $\tau_{1}\in\left[t_{K_{1}},t_{K_{1}+1}\right]$.

  1. 2.

If $k_{M}-k_{m}=N-2$, via Lemma 2 there are two options

    1. (a)

$k_{m}=K_{1}-N+1,k_{M}=K_{1}-1$. Here (25) based on the same reasoning as in
1).

    2. (b)

$k_{m}=K_{1}-N+2,k_{M}=K_{1}$. Here we have that

$\widetilde{s}_{1}={\mathrm{sign}\left(f(p)\right)}\cdot s_{1},$ (50)

where,

       * •

$f\left(p\right)=\bar{\mu}_{2}p+\bar{\beta}_{2}\left(1-p\right)$.

       * •

$\bar{\mu}_{2}=\mu_{K_{1}}^{N}\left[K_{1}-N+2\right]$, and,

       * •

$\bar{\beta}_{2}=\beta_{K_{1}}^{N}\left[K_{1}-N+2\right]$.

We know that
$\left\lvert\mathscr{D}^{N}{Z}\left(t_{K_{1}-N+1}\right)\right\rvert<\Psi^{N}\left[K_{1}-N+1\right]$,
which implies that
$\left\lvert\mathscr{D}^{N}E_{g}\left(t_{K_{1}-N+1}\right)\right\rvert<2\Psi^{N}\left[K_{1}-N+1\right]$
via (44). This can be rewritten as

$\displaystyle\frac{\left\lvert\mathscr{D}^{N}E_{g}\left(t_{K_{1}-N+1}\right)\right\rvert}{2\lambda_{h}}=\left\lvert\bar{\mu}_{1}\right\rvert\cdot
p$
$\displaystyle<\frac{\Psi^{N}\left[K_{1}-N+1\right]}{\lambda_{h}}\leqslant\phi_{0,K_{1}},$

where $\bar{\mu}_{1}=\mu_{K_{1}}^{N}\left[K_{1}-N+1\right]$. Using (21), this
gives us $p<x^{\ast\ast}$, where
$x^{\ast\ast}=\tfrac{\left\lvert\bar{\beta}_{2}\right\rvert}{\left\lvert\bar{\beta}_{2}\right\rvert+\left\lvert\bar{\mu}_{2}\right\rvert+\left\lvert\bar{\mu}_{1}\right\rvert}$.
Moreover, using that $\mathrm{sign}\left(\bar{\mu}_{1}\right)=1$ from 1) and
(39), it can be derived that $f\left(p\right)$ is a strictly decreasing line
and $f\left(x^{\ast\ast}\right){>}0$. Given that $p<x^{\ast\ast}$ we get
$f(p)>0$ which implies $\widetilde{s}_{1}=s_{1}$ via (50). Moreover,
$\widetilde{\tau}_{1}=t_{k_{M}+1}=t_{K_{1}+1}$ which satisfies (25).

  2. 3.

If $k_{M}-k_{m}=N-3$, Lemma 2 implies $k_{m}=K_{1}-N+2,k_{M}=K_{1}-1$. Then
$\widetilde{s}_{1}=s_{1}$ via the same derivation as in 2) b). Moreover,
$\widetilde{\tau}_{1}=\frac{t_{K_{1}}+t_{K_{1}+1}}{2}$ which satisfies (24).

∎

###### Proof of Proposition 3.

Assuming the values $\Delta G\left(t_{k}\right)=\left(\mathcal{L}g\right)_{k}$
known, $g$ can be recovered recursively as [45, 46]

$\displaystyle\begin{split}{g}_{0}\left(t\right)&=\mathcal{S}_{\Omega}\left[{\Delta{G}\left(t_{k}\right)}\right]\left(t\right),\\\
{g}_{n+1}\left(t\right)&={g}_{n}\left(t\right)+{g}_{0}\left(t\right)-\mathcal{S}_{\Omega}\left[{\mathcal{L}{g}_{n}}\right]\left(t\right),n>0.\end{split}$
(51)

Then, provided that $T_{\mathsf{max}}<\frac{\pi}{\Omega}$ the following hold
[37]

$\displaystyle\|g-g_{n}\|_{L^{2}}\leqslant\left(\frac{T_{\mathsf{max}}\Omega}{\pi}\right)^{n+1}\|g\|_{L^{2}},$
(52)
$\displaystyle\|\mathcal{I}-\mathcal{S}_{\Omega}\mathcal{L}\|_{L^{2}}\leqslant\frac{T_{\mathsf{max}}\Omega}{\pi},$
(53)

where $\mathcal{I}$ is the identity operator, and thus
$\lim_{n\rightarrow\infty}g_{n}=g$.

When comparing (51) with (34), the main difference is in the initial condition
$g_{0}$, which is determined by samples $\Delta{G}\left(t_{k}\right)$. In the
following, we evaluate the error in computing these samples. We begin by
making two observations.

_Observation 1._ As derived at the end of the Proof for Theorem 1,

$\widetilde{\tau}_{r}\in\left\\{t_{K_{r}+1},\frac{t_{K_{r}}+t_{K_{r}+1}}{2}\right\\},$

and thus, via the definition of $K_{r}$, we get
${\tau}_{r},\widetilde{\tau}_{r}\in\left[t_{K_{r}},t_{K_{r}+1}\right]$, i.e.,
there are no spike times in between the true and estimated fold.

_Observation 2._ Because $N\frac{2\delta\Omega g_{\infty}}{b-\lambda}\leqslant
h^{*}$ (condition $S_{2}$ in Lemma 3) we are guaranteed to have at least two
spike times between each two consecutive folding times $\tau_{r},\tau_{r+1}$.

Using (12) and that $X\left(t_{k}\right)$ is known, we get that

$\left\lvert\Delta\widetilde{G}\left(t_{k}\right)-\Delta
G\left(t_{k}\right)\right\rvert=\left\lvert\Delta\widetilde{E}_{g}\left(t_{k}\right)-\Delta
E_{g}\left(t_{k}\right)\right\rvert.$

We will derive a bound for the latter by considering two cases.

  1. 1.

$k\in\mathbb{Z}_{+}^{*}\setminus\left\\{K_{1},\dots,K_{R}\right\\}$. We first
consider the trivial case where $k<K_{1}\Rightarrow
E_{g}\left(t_{k}\right)=\widetilde{E}_{g}\left(t_{k}\right)=0$ (Observation
1). If $k>K_{1}$ then $\exists r_{0}\in\left\\{1,\dots,R\right\\}$ s.t.
$k\in\left\\{K_{r_{0}}+1,\dots,K_{r_{0}+1}-1\right\\}$ (Observation 2). Using
the expression of $E_{g}$ we get that

$\displaystyle\Delta E_{g}\left(t_{k}\right)$
$\displaystyle=\sum\limits_{r=1}^{r_{0}}2\lambda_{h}s_{r}\Delta\left(\mathbbmtt{1}_{\left[\tau_{r},\infty\right)}\left(t_{k}\right)\left(t_{k}-\tau_{r}\right)\right)$
(54)
$\displaystyle=\sum\limits_{r=1}^{r_{0}}2\lambda_{h}s_{r}\left(t_{k+1}-t_{k}\right).$
(55)

The last equality uses that there are no folds in $\left[t_{k},t_{k+1}\right]$
given that $k\neq K_{r_{0}}$. Using the same derivation for
$\Delta\widetilde{E}_{g}$ and that $s_{r}=\widetilde{s}_{r},\forall
r\in\left\\{1,\dots,R\right\\}$, we get that
$\Delta\widetilde{E}_{g}\left(t_{k}\right)=\Delta E_{g}\left(t_{k}\right)$.

  2. 2.

$k\in\left\\{K_{1},\dots,K_{R}\right\\}$. Let
$k=K_{r_{0}},r_{0}\in\left\\{1,\dots,R\right\\}$. Due to Observation 1 we have
$\tau_{r_{0}},\widetilde{\tau}_{r_{0}}\in\left[t_{K_{r_{0}}},t_{K_{r_{0}}+1}\right]$.
Then, using that
$\mathbbmtt{1}_{\left[\tau_{r_{0}},\infty\right)}\left(t_{k}\right)=0$ in (54)
we get

$\Delta
E_{g}\left(t_{k}\right)=\sum_{r=1}^{r_{0}-1}2\lambda_{h}s_{r}\left(t_{k+1}-t_{k}\right)+2\lambda_{h}s_{r_{0}}\left(t_{k+1}-\tau_{r_{0}}\right)$
(56)

Therefore $\left\lvert\Delta
E_{g}\left(t_{k}\right)-\Delta\widetilde{E}_{g}\left(t_{k}\right)\right\rvert=2\lambda_{h}\left\lvert\tau_{r_{0}}-\widetilde{\tau}_{r_{0}}\right\rvert$.

By combining 1) and 2), we get:

$\left\lvert\Delta\widetilde{G}\left(t_{k}\right)-\Delta
G\left(t_{k}\right)\right\rvert\leqslant
2\lambda_{h}\max_{r}\left\lvert\tau_{r}-\widetilde{\tau}_{r}\right\rvert.$
(57)

To bound the estimation of $g_{n}$, we define
$e_{n}\triangleq\widetilde{g}_{n}-g_{n}\Rightarrow
e_{0}=\mathcal{S}_{\Omega}\left[{\Delta\widetilde{G}\left(t_{k}\right)-\Delta
G\left(t_{k}\right)}\right]$ and thus

$e_{0}\left(t\right)=\sum_{k\in\mathbb{K}_{r}}\left[\Delta\widetilde{G}\left(t_{k}\right)-\Delta
G\left(t_{k}\right)\right]\cdot\mathrm{sinc}_{\Omega}\left(t-s_{k}\right).$
(58)

Given that $\left\lvert\mathbb{K}_{r}\right\rvert\leqslant 2R$, via the
triangle inequality and (57) we have that

$\displaystyle\begin{split}\|e_{0}\|_{L^{2}}&\leqslant
2R\max_{k}\left\lvert\Delta\widetilde{G}\left(t_{k}\right)-\Delta
G\left(t_{k}\right)\right\rvert\cdot\|\mathrm{sinc}_{\Omega}\|_{L^{2}}\\\
&\leqslant 4\lambda_{h}R\cdot
T_{\mathsf{max}}\sqrt{\frac{\Omega}{\pi}}.\end{split}$ (59)

Furthermore, we have that (34, 51)

$\displaystyle\|e_{n}\|_{L^{2}}$
$\displaystyle=\|e_{n-1}+e_{0}-\mathcal{S}_{\Omega}\mathcal{L}e_{n-1}\|_{L^{2}}$
(60)
$\displaystyle\leqslant\|e_{0}\|_{L^{2}}+\|\left(\mathcal{I}-\mathcal{S}_{\Omega}\mathcal{L}\right)e_{n-1}\|_{L^{2}}$
(61)
$\displaystyle\overset{\eqref{eq:op_bound}}{\leqslant}\|e_{0}\|_{L^{2}}+\|e_{n-1}\|_{L^{2}}\cdot\frac{T_{\mathsf{max}}\Omega}{\pi}.$
(62)

Via induction, given that $T_{\mathsf{max}}<\frac{\pi}{\Omega}$, we get

$\displaystyle\|e_{n}\|_{L^{2}}$
$\displaystyle\leqslant\|e_{0}\|_{L^{2}}\sum_{l=0}^{n}\left(\frac{T_{\mathsf{max}}\Omega}{\pi}\right)^{l}=\|e_{0}\|_{L^{2}}\frac{1-\left(\frac{T_{\mathsf{max}}\Omega}{\pi}\right)^{n+1}}{1-\frac{T_{\mathsf{max}}\Omega}{\pi}}\leqslant\frac{\|e_{0}\|_{L^{2}}}{1-\frac{T_{\mathsf{max}}\Omega}{\pi}}.$
(63)

Finally, we bound the reconstruction error as

$\|\widetilde{g}_{n}-g\|_{L^{2}}\leqslant\|\widetilde{g}_{n}-g_{n}\|_{L^{2}}+\|g-g_{n}\|_{L^{2}},$
(64)

which yields the desired result via (59,63,52), where $T_{\mathsf{max}}$
satisfies (32). ∎

##  VII Conclusions and Future Work

Summary of Results. We proposed a new encoding model based on the modulo
event-driven sampling (MEDS) pipeline consisting of a modulo-hysteresis
nonlinearity in series with an asynchronous sigma-delta modulator (ASDM)
model. Our result contributes to the existing efforts of alleviating the
dynamic range restrictions of existing EDS sensors. We provide mathematical
guarantees for the input reconstruction, and validate the MEDS model with
synthetic and hardware experiments. Numerical simulations show that the
proposed method can recover inputs with a dynamic range up to $20$ times the
modulo threshold.

Future Work. We plan on following up the results in several directions,

  1. 1.

The method can be extended to modulo models in series with other event-driven
encoders, such as level-crossing encoders and biphasic integrate-and-fire
neurons.

  2. 2.

The nonuniform finite difference could be replaced by a more general operator
for achieving a better recovery.

  3. 3.

The current methodology can be extended to inputs in multiple dimensions.

New exciting hardware applications such as event-driven cameras show the
promise that this encoding scheme holds in the field of signal processing.
Despite this, the sensor dynamic range restriction remains a fundamental
bottleneck in the field. Recently, it was shown that a high dynamic range
image can be recovered from uniform samples in a single capture [22]. Given
that event-driven signals have an inherently low dynamic range, this line of
research represents a natural step in the field that could lead to promising
new encoding strategies.

## Acknowledgements

The authors thank AB’s undergraduate mentee Jerry Chen for help with the
experiments, specially during the pandemic and the anonymous reviewers whose
comments have improved the readability of our work.

## References

  * [1] A. A. Lazar and L. T. Tóth, “Perfect recovery and sensitivity analysis of time encoded bandlimited signals,” _IEEE Trans. Circuits Syst. I_ , vol. 51, no. 10, pp. 2060–2073, 2004.
  * [2] A. A. Lazar and E. A. Pnevmatikakis, “Faithful representation of stimuli with a population of integrate-and-fire neurons,” _Neural Comput._ , vol. 20, no. 11, pp. 2715–2744, 2008.
  * [3] D. Gontier and M. Vetterli, “Sampling based on timing: Time encoding machines on shift-invariant subspaces,” _Applied and Computational Harmonic Analysis_ , vol. 36, no. 1, pp. 63–78, 2014.
  * [4] P. Martínez-Nuevo, H.-Y. Lai, and A. V. Oppenheim, “Delta-ramp encoder for amplitude sampling and its interpretation as time encoding,” _IEEE Trans. Signal Process._ , vol. 67, no. 10, pp. 2516–2527, 2019.
  * [5] M. Miskowicz, _Event-based control and signal processing_. CRC press, 2018.
  * [6] G. Gallego, T. Delbruck, G. M. Orchard, C. Bartolozzi, B. Taba, A. Censi, S. Leutenegger, A. Davison, J. Conradt, K. Daniilidis, and D. Scaramuzza, “Event-based vision: A survey,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. In Press, pp. 1–26, 2020.
  * [7] W. Maass, “Networks of spiking neurons: the third generation of neural network models,” _Neural Networks_ , vol. 10, no. 9, pp. 1659–1671, 1997.
  * [8] D. Florescu and D. Coca, “Learning with precise spike times: A new decoding algorithm for liquid state machines,” _Neural Computation_ , vol. 31, no. 9, pp. 1825–1852, 2019.
  * [9] K. Ozols, M. Greitans, and R. Shavelis, “EEG data acquisition system based on asynchronous sigma-delta modulator,” in _Biennial Baltic Electronics Conf._ , 2012, pp. 183–186.
  * [10] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski, _Neuronal dynamics: From single neurons to networks and models of cognition_. Cambridge University Press, 2014.
  * [11] M. Renaudin, “Asynchronous circuits and systems: a promising design alternative,” _Microelectronic engineering_ , vol. 54, no. 1-2, pp. 133–149, 2000.
  * [12] D. Wang, S. J. Kim, M. Yang, A. A. Lazar, and M. Seok, “A background-noise and process-variation-tolerant 109nw acoustic feature extractor based on spike-domain divisive-energy normalization for an always-on keyword spotting device,” in _IEEE Intl. Solid-State Circuits Conference (ISSCC)_ , vol. 64, 2021, pp. 160–162.
  * [13] H. G. Feichtinger and K. Gröchenig, “Theory and practice of irregular sampling,” _Wavelets: Mathematics and Applications_ , vol. 1994, pp. 305–363, 1994.
  * [14] A. Bhandari, F. Krahmer, and R. Raskar, “On unlimited sampling,” in _Intl. Conf. on Sampling Theory and Applications (SampTA)_ , 2017, pp. 31–35.
  * [15] A. Bhandari and F. Krahmer, “On identifiability in unlimited sampling,” in _Proc. Int. Conf. Sampling Theory Appl._ IEEE, 2019, pp. 1–4.
  * [16] A. Bhandari, F. Krahmer, and R. Raskar, “On unlimited sampling and reconstruction,” _IEEE Trans. Signal Process._ , vol. 69, pp. 3827 – 3839, Dec. 2020.
  * [17] ——, “Methods and apparatus for modulo sampling and recovery,” US Patent US10 651 865B2, May, 2020.
  * [18] A. Bhandari, F. Krahmer, and T. Poskitt, “Unlimited sampling from theory to practice: Fourier-Prony recovery and prototype ADC,” _IEEE Trans. Signal Process._ , Sep. 2021.
  * [19] A. Bhandari, F. Krahmer, and R. Raskar, “Unlimited sampling of sparse signals,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc._ , 2018, pp. 4569–4573.
  * [20] A. Bhandari, “Back in the US-SR: Unlimited sampling and sparse super-resolution with its hardware validation,” _IEEE Signal Process. Lett._ , vol. 29, pp. 1047–1051, 2022.
  * [21] A. Bhandari, F. Krahmer, and R. Raskar, “Unlimited sampling of sparse sinusoidal mixtures,” in _IEEE Intl. Symp. on Inf. Theory (ISIT)_ , 2018, pp. 336–340.
  * [22] A. Bhandari and F. Krahmer, “HDR imaging from quantization noise,” in _IEEE Intl. Conf. on Image Proc. (ICIP)_ , 2020, pp. 101–105.
  * [23] D. Florescu, F. Krahmer, and A. Bhandari, “The surprising benefits of hysteresis in unlimited sampling: Theory, algorithms and experiments,” _IEEE Trans. Signal Process._ , vol. 70, pp. 616–630, 2022.
  * [24] O. Graf, A. Bhandari, and F. Krahmer, “One-bit unlimited sampling,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc._ , 2019, pp. 5102–5106.
  * [25] D. Florescu, F. Krahmer, and A. Bhandari, “Event-driven modulo sampling,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , 2021.
  * [26] D. Florescu and A. Bhandari, “Unlimited sampling with local averages,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , 2022.
  * [27] S. Fernández-Menduina, F. Krahmer, G. Leus, and A. Bhandari, “Doa estimation via unlimited sensing,” in _28th European Signal Processing Conference (EUSIPCO)_ , 2020, pp. 1866–1870.
  * [28] S. Fernandez-Menduina, F. Krahmer, G. Leus, and A. Bhandari, “Computational array signal processing via modulo non-linearities,” _IEEE Trans. Signal Process._ , Jul. 2021.
  * [29] A. Bhandari, M. Beckmann, and F. Krahmer, “The modulo Radon transform and its inversion,” in _28th European Sig. Proc. Conf. (EUSIPCO)_ , 2020, pp. 770–774.
  * [30] M. Beckmann, F. Krahmer, and A. Bhandari, “The modulo Radon transform: Theory, algorithms and applications,” _SIAM Journal of Imaging Sciences_ , vol. 15, no. 2, 2022.
  * [31] E. Roza, “Analog-to-digital conversion via duty-cycle modulation,” _IEEE Trans. Circuits Syst. II_ , vol. 44, no. 11, pp. 907–914, 1997\.
  * [32] K. Ozols, “Implementation of reception and real-time decoding of ASDM encoded and wirelessly transmitted signals,” in _25th International Conference Radioelektronika_ , 2015, pp. 236–239.
  * [33] D. Wei, V. Garg, and J. G. Harris, “An asynchronous delta-sigma converter implementation,” in _IEEE Intl. Symposium on Circuits and Systems_ , 2006, pp. 4903–4906.
  * [34] R. Shavelis, K. Ozols, and M. Greitans, “Amplitude adaptive ASDM circuit,” in _3rd International Conference on Event-Based Control, Communication and Signal Processing (EBCCSP)_ , 2017, pp. 1–8.
  * [35] R. Alexandru and P. L. Dragotti, “Reconstructing classes of non-bandlimited signals from time encoded information,” _IEEE Trans. Signal Process._ , vol. 68, pp. 747–763, 2019.
  * [36] S. Rudresh, A. J. Kamath, and C. S. Seelamantula, “A time-based sampling framework for finite-rate-of-innovation signals,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc._ , 2020, pp. 5585–5589.
  * [37] D. Florescu and D. Coca, “A novel reconstruction framework for time-encoded signals with integrate-and-fire neurons,” _Neural Comput._ , vol. 27, no. 9, pp. 1872–1898, 2015.
  * [38] H. G. Feichtinger, J. C. Príncipe, J. L. Romero, A. S. Alvarado, and G. A. Velasco, “Approximate reconstruction of bandlimited functions for the integrate and fire sampler,” _Adv. Comput. Math._ , vol. 36, no. 1, pp. 67–78, 2012.
  * [39] D. Florescu and A. Bhandari, “Modulo event-driven sampling: System identification and hardware experiments,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP)_ , 2022.
  * [40] S. Rudresh, A. Adiga, B. A. Shenoy, and C. S. Seelamantula, “Wavelet-based reconstruction for unlimited sampling,” in _IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc._ , 2018, pp. 4584–4588.
  * [41] J. D. Pereira, A. C. Serra, and P. S. Girao, “Dithered ADC systems in the presence of hysteresis errors,” in _Proc. of the IEEE Instrumentation and Meas. Tech. Conf._ , vol. 3, 1999, pp. 1648–1652.
  * [42] “IEEE standard for digitizing waveform recorders,” _IEEE Std 1057-2017 (Revision of IEEE Std 1057-2007)_ , pp. 1–179, 2018.
  * [43] S. M. Nikol’skii, _Approximation of functions of several variables and imbedding theorems_. Springer Science & Business Media, 2012, vol. 205.
  * [44] A. A. Lazar and E. A. Pnevmatikakis, “Consistent recovery of sensory stimuli encoded with mimo neural circuits,” _Computational Intelligence and Neuroscience_ , Feb 2010.
  * [45] D. Florescu, _Reconstruction, identification and implementation methods for spiking neural circuits_. Springer, 2017.
  * [46] A. A. Lazar, “Time encoding with an integrate-and-fire neuron with a refractory period,” _Neurocomputing_ , vol. 58, pp. 53–58, 2004.


# On the theory of ideal Bose-gas at a finite particle number

A.I. Bugrij 1111E-mail<EMAIL_ADDRESS>V.M. Loktev1,2777E-mail:
<EMAIL_ADDRESS>
1Bogolyubov Institute for Theoretical Physics of the NAS of Ukraine,
14-b, Metrolohichna Str., Kyiv 03143, Ukraine
2National Technical University of Ukraine “Igor Sikorsky Kyiv Politechnical
Institute“
37, Peremogy Ave., Kyiv 03056, Ukraine

The ideal Bose-gas with finite number $N$ of particles is investigated. The
exact expressions for the partition functions and occupation numbers in the
grand canonical, canonical and microcanonical ensembles are found. The
asymptotic expressions (in the case $N\gg 1$) for the partition functions and
occupation numbers in the canonical and microcanonical ensembles are
evaluated. It is shown that the chemical potential $\mu$ of the ideal Bose-gas
can lie in the range $-\infty<\mu<\infty$ oppositely to the widely adopted
opinion that the value of this potential is negative.

PACS: 05.30.Jp, 05.30.-d, 05.30.Ch

Keywords: ideal Bose-gas, Bose-distribution, canonical ensembles.

##  1 Introduction

The theory of ideal Bose-gas is referred to the old and rather well studied
field of statistical physics of macroscopic systems [1-4]. However, the real
experiments aimed at its verification are performed with systems consisting of
a finite number $N$ of the particles. For example, in the experiments with
atomic Bose-gases [5-8] whose results are usually interpreted as the
experimental confirmation of the phenomenon of Bose-condensation, the number
of particles is at most $\sim 10^{4}$, or does not attain the Avogadro number.
At the same time, no corresponding calculations for the ideal Bose-gas
consisting of a finite number of particles are available in the literature.
Therefore, the goal of the our paper is to fill this gap. Moreover, we will
calculate the observable quantities for all basic statistical ensembles,
namely, the grand canonical (GCE), canonical (CE), and microcanonical (MCE)
ones.

##  2 Ideal Bose-gas

A quantum particle located in the vessel with impermeable walls with volume
$V$ has the discrete energy spectrum $\varepsilon_{k}$, where the index $k$
runs the values $0,1,\ldots,\infty,$ as the energy $\varepsilon_{k}$ of the
state increases.

Consider the system consisting of $N$ bosons non-interacting with one another
that are placed into the vessel. The system is characterized by a
configuration $[n]$, i.e., by a set of the occupation numbers
$\mathbf{n}_{k}$. Each of them indicates the number of particles in the $k$-th
state with the energy $\varepsilon_{k}$. It is obvious that $\mathbf{n}_{k}$
_cannot exceed_ the number $N$ of particles in the system. In this case, their
total number and total energy of the system are given by the equalities:

$N=\sum_{k=0}^{\infty}{\mathbf{n}}_{k},\quad
E=\sum_{k=0}^{\infty}{\mathbf{n}}_{k}\varepsilon_{k}.$ (1)

Let us place the vessel into a thermostat with the temperature $T$. Under the
action of thermal fluctuations on the vessel walls, the configuration $[n]$
will vary (stochastically) with the time. As is known, from experiments one
obtain the time-averaged values of the observable quantities. In statistical
physics for their calculations the ergodic hypothesis is as a rule accepted.
According to it, the mean over the time coincides with the mean over the
ensemble.

By definition, the ensemble is a collection of systems with all possible
configurations or, in other words, with different distributions of particles
over states. The ensembles can be very different. Each specific ensemble is
determined by the own distribution function $f[n]$, and the most popular among
them are, as mentioned above, GCE, CE, and MCE.

##  3 Grand canonical ensemble

The distribution function for GCE takes the form

$f[n]={\mathrm{e}}^{-\beta(E-\mu
N)}=\prod_{k}{\mathrm{e}}^{-\mathbf{n}_{k}\beta(\varepsilon_{k}-\mu)}\,,$ (2)

where $\beta=(k_{\scriptstyle{B}}T)^{-1}$ stands for the reciprocal
temperature, and $\mu$ is as above the chemical potential.

Since $f[n]$ (2) is factorized, the partition function of such system can be
easily calculated [9]:

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\sum_{[n]}f[n]=\sum_{\mathbf{n}_{0}=0}^{N}{\mathrm{e}}^{-\mathbf{n}_{0}\beta(\varepsilon_{0}-\mu)}\sum_{\mathbf{n}_{1}=0}^{N}{\mathrm{e}}^{-\mathbf{n}_{1}\beta(\varepsilon_{1}-\mu)}\ldots\sum_{\mathbf{n}_{k}=0}^{N}{\mathrm{e}}^{-\mathbf{n}_{k}\beta(\varepsilon_{k}-\mu)}\ldots$
(3) $\displaystyle=$
$\displaystyle\prod_{k}\frac{1-{\mathrm{e}}^{-\beta(\varepsilon_{k}-\mu)(N+1)}}{1-{\mathrm{e}}^{-\beta(\varepsilon_{k}-\mu)}}.$

The mean value of the number of occupation can be also determined and takes
the form

$\overline{n}_{k}=Z^{-1}\sum_{[\mathbf{n}]}\mathbf{n}_{k}f[\mathbf{n}]=n_{k}-(N+1){m}_{k},$
(4)

where $n_{k}$ and ${m}_{k}$ are defined as follows:

$n_{k}=\frac{1}{{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)}-1},\quad{m}_{k}=\frac{1}{{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)(N+1)}-1}.$

It is seen that the first term on the right-hand side of (4) is the ordinary
Bose-distribution, and the second one introduces the dependence on the number
of particles of the system into the average over the ensemble. This point is
worth to be noted, since the mean number of occupation following from
relations (3) and (4) coincides formally with that for the parastatistics [3]
under the condition $p=N$, though it is quite obvious that the order $p$ of
the parastatistics and the number $N$ are different physical quantities
without any connection with each other. In addition, the derivation of the
corresponding formula in [3] seems not to be quite proper, because the
Stirling factorial formula is used in it for the quantities less than 1.
Nevertheless, the final result turns out proper and coincides with formula (4)
which is obtained with the help of exact calculations.

It is worth to note that the number $\overline{n}_{k}$ has a finite value for

$-\infty<\mu<\infty,$

rather than only for $\mu<0$. So, the latter is the assertion which is not
absolutely true. Indeed, the number of occupation for the lowest state

$\overline{n}_{0}=\frac{1}{{\mathrm{e}}^{\beta(\varepsilon_{0}-\mu)}-1}-\frac{N+1}{{\mathrm{e}}^{\beta(\varepsilon_{0}-\mu)(N+1)}-1}.$
(5)

From whence, we get $\overline{n}_{0}=N$ for $\mu\rightarrow\infty$,
$\overline{n}_{0}={N/2}$ for $\mu\rightarrow\varepsilon_{0}$, and
$\overline{n}_{0}=0$ for $\mu\rightarrow-\infty$. In this case, the parameters
$\beta$ and $\mu$ can be presented in terms of $N$ and $E$ by means of the
system of equations (cf. (1))

$N=\sum_{k=0}^{\infty}\overline{n}_{k},\quad
E=\sum_{k=0}^{\infty}\overline{n}_{k}\varepsilon_{k}.$ (6)

##  4 Canonical ensemble

Unlike the previous case, the distribution function for CE takes the form

$f_{CE}[n]=f[n]\delta\biggl{(}N-\sum_{k}\mathbf{n}_{k}\biggr{)}.$ (7)

The calculation of the partition function becomes harder due to the presence
of the $\delta$-function in (7), but the exact analytic calculation can be
carried out.

We represent $\delta(N-\sum_{k}\mathbf{n}_{k})$ as the integral:

$\delta(N-\sum_{k}\mathbf{n}_{k})=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}dx\,{\mathrm{e}}^{ixN}\prod_{k}{\mathrm{e}}^{-ix\mathbf{n}_{k}}.$
(8)

Then the partition function of CE takes the form

$Z_{{CE}}=\sum_{[n]}f_{{CE}}[n]=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}dx\,{\mathrm{e}}^{-w(x)}\,$
(9)

where

$w(x)=-ixN-\sum_{k}\ln\frac{1-{\mathrm{e}}^{-[\beta(\varepsilon_{k}-\mu)+ix](N+1)}}{1-{\mathrm{e}}^{-[\beta(\varepsilon_{k}-\mu)+ix]}},$
(10)
$\overline{n}_{k}^{{CE}}=\frac{Z_{{CE}}^{-1}}{2\pi}\int\limits_{-\pi}^{\pi}dx\,{\mathrm{e}}^{-w(x)}\left(\frac{1}{{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)+ix}-1}-\frac{N+1}{{\mathrm{e}}^{[\beta(\varepsilon_{k}-\mu)+ix](N+1)}-1}\right).$
(11)

For $N\gg 1,$ integrals (9) and (11) can be calculated with the help of the
saddle-point method. Corresponding saddle point is determined from the
equation

$\frac{\partial w(x)}{\partial
x}=0=-i\left(N-\sum_{k}\biggl{(}\frac{1}{{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)+ix}-1}-\frac{N+1}{{\mathrm{e}}^{[\beta(\varepsilon_{k}-\mu)+ix](N+1)}-1}\biggr{)}\right)\,$

which holds, due to equalities (6), for $x=0$. The second derivative at the
saddle point takes the form

$w_{xx}=\left.\frac{\partial^{2}w(x)}{\partial
x^{2}}\right|_{x=0}=\sum_{k}\left(\frac{{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)}}{({\mathrm{e}}^{\beta}(\varepsilon_{k}-\mu)-1)^{2}}-\frac{(N+1)^{2}{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)(N+1)}}{({\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)(N+1)}-1)^{2}}\right).$
(12)

As a result, if the condition $N\gg 1$ holds, we obtain

$Z_{CE}=\frac{Z}{\sqrt{2\pi w_{xx}}}\,.$ (13)

The mean value of the number of occupation in CE for $N\gg 1$ is

$\overline{n}_{k}^{CE}=-\beta^{-1}\frac{\partial\ln
Z_{CE}}{\partial\varepsilon_{k}}=\overline{n}_{k}-m_{k}^{CE},$ (14)

where

$m_{k}^{CE}=\frac{1}{2w_{xx}}\left(n_{k}(n_{k}+1)(2n_{k}+1)-(N+1)^{3}{m}_{k}({m}_{k}+1)(2{m}_{k}+1)\right),$
(15)

and the quantities $n_{k}$ and $m_{k}$ are defined in (4). The direct
verification gives that the numbers of occupation in the ground state and in
the excited ones in CE are, respectively, larger and less, than in GCE.

##  5 Microcanonical ensemble

The distribution function for MCE is easily set and takes the form

$f_{MCE}[n]=f[n]\delta(N-\sum_{k}\mathbf{n}_{k})\delta(E-\sum_{k}\mathbf{n}_{k}\varepsilon_{k})\,.$
(16)

For the second $\delta$-function, we also use the integral representation
analogous to (8), namely:

$\delta\bigl{(}E-\sum_{k}\mathbf{n}_{k}\varepsilon_{k}\bigr{)}=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}dy{\mathrm{e}}^{iyE}\prod_{k}{\mathrm{e}}^{-iy\mathbf{n}_{k}\varepsilon_{k}}\,,$
(17)

and can easily sum over the configurations. Then, for the partition function
in MCE, we get (cf. (9))

$Z_{MCE}=\sum_{[n]}f_{MCE}[n]=\frac{1}{(2\pi)^{2}}\int\limits_{-\pi}^{\pi}dxdy\,{\mathrm{e}}^{-w(x,y)},$
(18)

where

$w(x,y)=-iNx-
iEy-\sum_{k}\ln\frac{1-{\mathrm{e}}^{-[\beta(\varepsilon_{k}-\mu)+ix+iy\varepsilon_{k}](N+1)}}{1-{\mathrm{e}}^{-[\beta(\varepsilon_{k}-\mu)+ix+iy\varepsilon_{k}]}}\,.$
(19)

In this case, the mean value of the number of occupation is presented in terms
of the integral

$\displaystyle\overline{n}_{k}^{MCE}$ $\displaystyle=$
$\displaystyle\frac{Z^{-1}_{MCE}}{(2\pi)^{2}}\int\limits_{-\pi}^{\pi}dxdy\,{\mathrm{e}}^{-w(x,y)}\times$
(20) $\displaystyle\times$
$\displaystyle\left(\frac{1}{{\mathrm{e}}^{\beta(\varepsilon_{k}-\mu)+ix+iy\varepsilon_{k}}-1}-\frac{N+1}{{\mathrm{e}}^{[\beta(\varepsilon_{k}-\mu)+ix+iy\varepsilon_{k}](N+1)}-1}\right),$

which can be considered as a generalization of integral (11). In the
calculation of the partition function (18), we take $N\gg 1$ and use the
saddle-point method whose point has the coordinates $x=y=0$ . As a result, we
obtain

$Z_{MCE}=\frac{1}{2\pi}\frac{Z}{\sqrt{d}},$ (21)

where

$d=w_{xx}w_{yy}-w_{xy}^{2},$ (22)
$w_{xx}=\left.\frac{\partial^{2}w(x,y)}{\partial x^{2}}\right|_{x=y=0},\quad
w_{xy}=\left.\frac{\partial^{2}w(x,y)}{\partial x\partial
y}\right|_{x=y=0},\quad w_{yy}=\left.\frac{\partial^{2}w(x,y)}{\partial
y^{2}}\right|_{x=y=0},$

$\overline{n}_{k}^{MCE}=\overline{n}_{k}-m_{k}^{MCE},$ (23)

and

$m_{k}^{MCE}=\frac{1}{2d}\left(\bigl{(}w_{yy}+\varepsilon_{k}^{2}w_{xx}\bigr{)}b_{k}-\bigl{(}2\varepsilon_{k}w_{xx}-w_{xy}\bigr{)}\frac{a_{k}}{\beta}\right),$
(24)

$a_{k}=n_{k}(n_{k}+1)-(N+1)^{2}{m}_{k}({m}_{k}+1),$

$b_{k}=n_{k}(n_{k}+1)(2n_{k}+1)-(N+1)^{3}{m}_{k}({m}_{k}+1)(2{m}_{k}+1).$

As is seen from formulas (22) and (23), the mean values of the numbers of
occupation of the ground state and of the excited ones in MCE are,
respectively, larger and less, than in the above-considered GCE and CE.

##  6 Conclusions

The above results are different from the commonly known ones (see, e.g.,
[1-4]) and supplement them. In our opinion, the cause lies in that the
previous studies (at least, the studies of systems in which the number of
particles (bosons, in this case) is given) neglected some obvious fact,
namely, the physical limitation imposed on the numbers of occupation
$\mathbf{n}_{k}\leq N$, which requires a special consideration.

Nevertheless, even for a sufficiently large number of bosons that usually
corresponds to the experiments on the Bose–Einstein condensation of cold
atomic gases, their thermodynamic quantities (in particular, the condensate
density) depend, as is shown above, on the number $N$. Therefore, we believe
that the presented results can be useful for the thermodynamics and
statistical physics of Bose-systems with finite numbers of particles.

We are grateful to A.S. Kovalev for his interest in the problem and the useful
discussions.

The work was performed in the frame of the programs 0117U00236 and 0116U003191
(State CPCEC 6541210 and 6541230) and the Special scientific program
0117U00240 of the Department of Physics and Astronomy of the NAS of Ukraine.

## References

  * [1] L.D. Landau and E.M. Lifshitz. Statistical Physics. Pergamon Press, New York (1980).
  * [2] K. Huang. Statistical Mechanics. Wiley, New York (1987).
  * [3] A. Isihara. Statistical Physics. Acad. Press, New York (1971).
  * [4] N.N.Bogolyubov. Selected Works, Vol. 2, Naukova Dumka, Kiev, 1970, p. 351.
  * [5] M.H. Anderson, J.R.Ensher, M.R.Vatthews, C.E.Wieman, and E.A.Cornell. Science 269, 198-201 (1995).
  * [6] K.Burnett, M.Edwards, and C.W.Clark. Physics Today 52, 37-42 (1999).
  * [7] Y.Shin, M.Saba, A.Schirotzek, T.A.Pasquini, A.E.Leanhardt, D.E.Prithard, and W.Kttterle. Phys. Rev. Lett. 92, 150401 (2004).
  * [8] G.K.Chaudhary, A.Chattopadhyay, and R.Ramakumar, Int. J. Mod. Phys. B 25, 3927-3940 (2012).
  * [9] A.I.Bugrij and V.M.Loktev. LTP 47, No. 2 (2021).