# Multi-Dimensional Unlimited Sampling
and Robust Reconstruction

Dorian Florescu and Ayush Bhandari

###### Abstract

In this paper we introduce a new sampling and reconstruction approach for
multi-dimensional analog signals. Building on top of the Unlimited Sensing
Framework (USF), we present a new folded sampling operator called the multi-
dimensional modulo-hysteresis that is also backwards compatible with the
existing one-dimensional modulo operator. Unlike previous approaches, the
proposed model is specifically tailored to multi-dimensional signals. In
particular, the model uses certain redundancy in dimensions $2$ and above,
which is exploited for input recovery with robustness. We prove that the new
operator is well-defined and its outputs have a bounded dynamic range. For the
noiseless case, we derive a theoretically guaranteed input reconstruction
approach. When the input is corrupted by Gaussian noise, we exploit redundancy
in higher dimensions to provide a bound on the error probability and show this
drops to $0$ for high enough sampling rates leading to new theoretical
guarantees for the noisy case. Our numerical examples corroborate the
theoretical results and show that the proposed approach can handle a
significantly larger amount of noise compared to USF.

###### keywords:

Analog-to-digital conversion, approximation, bandlimited functions, modulo
sampling, Shannon sampling.

††journal: Journal of LaTeX Templatesmytitlenotemytitlenotefootnotetext: D.
Florescu and A. Bhandari are with the Department of Electrical and Electronic
Engineering, Imperial College London, SW72AZ, UK. This work was supported by
the UK Research and Innovation council’s Future Leaders Fellowship program
“Sensing Beyond Barriers” (MRC Fellowship award no. MR/S034897/1). E-mails:
{D.Florescu<EMAIL_ADDRESS>or<EMAIL_ADDRESS>Project page
for (future) release of hardware design, code and data: https://bit.ly/USF-
Link.

###### Contents

  1. 1 Introduction
    1. 1.1 Unlimited Sensing Framework (USF)
  2. 2 Modulo acquisition and recovery
    1. 2.1 Recovery from one-dimensional modulo data
  3. 3 Multi-dimensional modulo sampling
    1. 3.1 Multi-dimensional lattice sampling preliminaries
    2. 3.2 Previous approaches for recovery from multi-dimensional data
    3. 3.3 Towards multi-dimensional modulo-hysteresis acquisition
    4. 3.4 Properties of the proposed operator
    5. 3.5 Problem formulation
  4. 4 Detecting modulo-hysteresis discontinuities
    1. 4.1 Detecting folding times and signs
    2. 4.2 Detecting constants $M_{\bar{\mathbf{b}}}$
  5. 5 Input reconstruction
    1. 5.1 Recovery with the proposed operator
    2. 5.2 Comparison to ideal modulo recovery with Gaussian noise measurements
  6. 6 Numerical study
  7. 7 Proofs
    1. 7.1 Multi-Dimensional Modulo Properties
    2. 7.2 Multi-Dimensional Modulo Recovery
  8. 8 Conclusion

##  1 Introduction

Shannon’s sampling theory is the workhorse of almost all modern-world digital
systems. Its practical implementation is carried out via electronic hardware,
namely, the analog-to-digital converter (ADC). However, there is a gap between
theory and practice which leads to a few fundamental deviations from the ideal
sampling model, including, among others, quantization (see the extensive
survey by Gray & Neuhoff [1]) non-pointwise sampling [2] and ADC saturation.
The latter deviation arises from the fact that the ADC is a physical device
and hence, one can only record a fixed range of amplitudes (typically, a
prescribed voltage range). This input amplitude range defines the _dynamic
range_ (or DR) of the ADC, say $\lambda>0$. Any signal exceeding (in absolute
value) $\lambda$ would result in permanent loss of information due to
saturation or clipping. Mathematically, this is synonymous to _hard
thresholding_ [3, 4], but the difference is that it occurs in hardware and is
highly undesirable. Clipped sample values lead to high frequency components,
which in turn leads to aliasing. Typical solutions to the saturation problem
rely on:

  1. (a)

hardware approaches such us _companding_ [5] or adaptively matching the
dynamic range to the input signal range via _automatic gain control_. There
are also techniques that re-think ADC design (see, for example, [6]).

  2. (b)

algorithmic approaches that aim to solve the inverse problem of _de-clipping_
[7, 8] or _inpainting_ [9].

Clipping or saturation is also highly relevant in the context of digital
imaging, so much so that almost all modern smartphones are equipped with the
_High Dynamic Range_ or “HDR” mode, based on multiple captures that are
combined numerically [10].

The progress in the last many decades has led to deepened understanding of the
nuances involved with the quantization and limited DR aspects. Clearly, ADCs
need to be matched to the DR of the input signal to avoid saturation or
clipping. Beyond this calibration step—typically addressed by the
engineers—there is an additional challenge: higher dynamic range requires a
higher number of bits to achieve a given resolution; this in turn leads to a
higher power consumption in the ADC, thus highlighting the integral role of DR
in digital acquisition.

Figure 1: Signal acquisition and recovery pipeline for Unlimited Sensing
Framework [11, 12]. (a) Modulo sampling hardware [13]. (b) Continuous-time
signal waveforms on an oscilloscope. (c) Signal reconstruction using recovery
algorithms [11, 12, 13, 14].

###  1.1 Unlimited Sensing Framework (USF)

Recently, the Unlimited Sensing Framework (USF) [11, 15, 16, 17, 18, 12, 13]
has been proposed in the literature that serves as an alternative digital
acquisition protocol for avoiding the DR limitation in conventional ADCs. The
USF is based on a joint design of hardware and mathematical algorithms.

  * •

In hardware, the modulo non-linearity ensures that HDR inputs are folded back
in to the ADC’s DR; this is because the modulo threshold is chosen such that
the modulo ADC’s range is bounded by $\lambda$. Consequently, the modulo ADC
results in _folded samples_.

  * •

To recover the HDR input from folded, modulo samples, mathematically
guaranteed recovery algorithms are deployed.

Similar to the Shannon–Nyquist sampling criterion where a higher input
bandwidth can be traded off for higher sampling rates, it was shown that HDR
signals can also be tackled by sampling more densely. This is made precise by
the following theorem.

######  Theorem 1 (Unlimited Sampling Theorem [11]).

Let $f(t)$ be a continuous-time function with maximum frequency $\Omega$
(rads/s). Then, a sufficient condition for recovery of $f(t)$ from its modulo
samples (up to an additive constant) taken every $T$ seconds apart is
$T\leqslant 1/\left(2\Omega\mathrm{e}\right)$ where $\mathrm{e}$ is Euler’s
number.

Thus, the USF addresses a major bottleneck in physical sensors by allowing the
recovery of inputs beyond the sensor dynamic range. A first validation of the
USF with experiments based on a modulo ADC were presented in [13]. In
particular, it was shown that signals as large as $\approx 25\lambda$ can be
recovered in a laboratory setup. The full sampling and reconstruction pipeline
for the USF is shown in Fig. 1.

The initial works based on USF tackled signals supported on the real line
spanned in a bandlimited [11, 12] or spline spaces [19]. There are also
methods to recover signals with compressive priors [20, 21] and using wavelet
filters [22]. Further extensions of the USF include compactly supported inputs
[13], sparse signals [23] and also new acquisition models [14].

A new acquisition model called modulo-hysteresis was introduced in [24] and
further discussed in [14], which considers hardware non-idealities and enables
new recovery guarantees. The modulo-hysteresis was also implemented in a
hardware prototype [14]. This line of work also paved the path to novel and
exciting low-power acquisition neuromorphic applications [24, 25].

The methods discussed so far assume that the input is one-dimensional .
However, in many applications, such as photography [19], X-ray imaging or
Computed Tomography [26], the input signal is multi-dimensional.

#### Motivation for a multi-dimensional model

There have been attempts to address multi-dimensional inputs with modulo
architectures by rasterizing and processing the signal line-by-line in imaging
[27, 19, 26] or for lattice sampling [28]. However, the existing modulo
sampling approaches for multi-dimensional data are based on a one-dimensional
modulo operator that is applied sequentially on the slices of a multi-
dimensional input. This considers each slice a distinct signal, and does not
exploit that they are all part of a multi-dimensional input. In other words,
this modulo operation represents a separable transformation that does not
exploit the multi-dimensional nature of the input. Furthermore, it is known
that non-separable transformations represent much more powerful tools in
analysing multi-dimensional data [29, 30].

We consider only the problem of input recovery for noisy inputs, distinct from
that of input denoising, which was addressed before for modulo sampling [31,
32]. Methods such as USF recover the noise corrupted input samples while
keeping the noise sequence intact. However, USF (Theorem 1) is fundamentally
restricted to work with noise amplitudes smaller than the modulo threshold.
When this requirement is not satisfied, the input recovery is heavily
distorted. This limitation is carried over to the existing attempts to apply
USF to multi-dimensional data.

#### Contributions

Here we present a modulo model that exploits the multi-dimensional structure
of the data in the encoding process. Specifically, via multi-dimensional
sampling in $D$ dimensions, we are able to dedicate a $D-1$-dimensional
subspace to deal with noise reduction, leaving dimension $d=1$ for estimating
the modulo folds. Specifically, our contributions are below:

  1. $\textrm{C}_{1})$

We introduce a generalized $D$-dimensional modulo operator for sampling on a
lattice.

  2. $\textrm{C}_{2})$

We prove that the operator is well-defined and the folding discontinuities are
located along directions given by the lattice vectors.

  3. $\textrm{C}_{3})$

We provide recovery guarantees under noiseless assumption.

  4. $\textrm{C}_{4})$

Under Gaussian noise assumption, we provide an upper bound on the recovery
error probability that drops to $0$ for high enough sampling rates.

  5. $\textrm{C}_{5})$

Using a numerical study we show that the proposed model offers significantly
better noise robustness than USF.

#### Notation

For $x\in\mathbb{R}$, $[\\![x]\\!]=x-\lfloor x\rfloor$ denotes the fractional
part of $x$ and $\left\lfloor x\right\rfloor$ is the floor function. For a set
$\mathbb{S}$, $\mathbbmtt{1}_{\mathbb{S}}$ is the indicator function and
$\mathrm{cl}\left(\mathbb{S}\right)$ is the set closure. The set of real and
integer numbers are $\mathbb{R}$ and $\mathbb{Z}$, respectively. Let
$\mathbb{R}^{*}=\mathbb{R}\setminus\left\\{0\right\\}$ and
$\mathbb{Z}^{*}=\mathbb{Z}\setminus\left\\{0\right\\}$, and let the sets
restricted to positive numbers $x\geqslant 0$ be $\mathbb{R}_{+}$ and
$\mathbb{Z}_{+}$. We denote by $\emptyset$ the empty set.

We use bold lowercase for vectors such as
$\mathbf{x}=\left[x_{1},\dots,x_{D}\right]^{\top}$, assumed to be column
vectors unless otherwise specified. Matrices are denoted by bold uppercase,
e.g., $\mathbf{M}\in\mathbb{R}^{D\times D}$. The element on line $k_{1}$ and
column $k_{2}$ in matrix $\mathbf{M}$ is denoted by
$\left[\mathbf{M}\right]_{k_{1},k_{2}}$. Unless specified otherwise, we use
notation $\mathbf{x}$ to denote a vector $\mathbf{x}\in\mathbb{Z}^{D}$ and
$\bar{\mathbf{x}}$ to denote a vector $\bar{\mathbf{x}}\in\mathbb{Z}^{D-1}$.
When used in the same context, $\bar{\mathbf{x}}$ denotes the last $D-1$
coordinates of $\mathbf{x}$ such that
$\mathbf{x}=\left[x_{1},\bar{\mathbf{x}}\right]$. Similarly, we denote by
$\bar{\mathbf{V}}\in\mathbb{R}^{D\times\left(D-1\right)}$ a matrix containing
the last $D-1$ columns of matrix $\mathbf{V}\in\mathbb{R}^{D\times D}$. For
two vectors $\mathbf{v}_{1},\mathbf{v}_{2}\in\mathbb{R}^{D}$ we denote their
inner product by
$\left\langle\mathbf{v}_{1},\mathbf{v}_{2}\right\rangle=\sum_{d=1}^{D}v_{1,d}\cdot\mathtt{conj}({v}_{2,d})=\mathbf{v}_{1}^{\top}\cdot\mathtt{conj}(\mathbf{v}_{1})$,
where $\mathtt{conj}$ is the complex-conjugate. Norm $\|\mathbf{v}\|_{2}$
denotes the Euclidean norm for a vector $\mathbf{v}\in\mathbb{R}^{D}$ and norm
$\|\mathbf{v}\|_{\infty}$ is defined as
$\|\mathbf{v}\|_{\infty}=\max_{d=1,\dots,D}\left\lvert\left[\mathbf{v}\right]_{d}\right\rvert$.
We denote by $\det\left(\mathbf{V}\right)$ the determinant of matrix
$\mathbf{V}$. We denote by
$\mathbf{T}=\mathrm{diag}\left\\{T_{1},\dots,T_{D}\right\\}$ a matrix with
$T_{1},\dots,T_{D}$ on the main diagonal and $0$ otherwise.

For a function $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$, $\|f\|_{2}$ and
$\|f\|_{\infty}$ represent the $L^{2}\left(\mathbb{R}^{D}\right)$ and
$L^{\infty}\left(\mathbb{R}^{D}\right)$ norms, respectively. We denote by
$\mathcal{F}f$ the Fourier transform applied to $f$, defined as

$\mathcal{F}f\left(\boldsymbol{\omega}\right)=\int_{\mathbb{R}^{D}}f\left(\mathbf{x}\right)e^{-\jmath\left\langle\boldsymbol{\omega},\mathbf{x}\right\rangle}d\mathbf{x},$
(1)

where $\mathbf{x}=\left[x_{1},\cdots,x_{D}\right]^{\top}$ and
$\boldsymbol{\omega}=\left[\omega_{1},\cdots,\omega_{D}\right]^{\top}$. The
inverse Fourier transform $\mathcal{F}^{-1}F$ is defined as

$\mathcal{F}^{-1}F\left(\mathbf{x}\right)=\frac{1}{\left(2\pi\right)^{D}}\int_{\mathbb{R}^{D}}F\left(\boldsymbol{\omega}\right)e^{\jmath\left\langle\boldsymbol{\omega},\mathbf{x}\right\rangle}d\boldsymbol{\omega}.$
(2)

The support of sequence $\psi$ is denoted by $\mathrm{supp}\left(\psi\right)$
and the support of a function $f\left(\mathbf{x}\right)$ is
$\mathrm{supp}\left(f\right)$. For two multi-dimensional sequences
$\gamma_{1},\gamma_{2}:\mathbb{Z}^{D}\rightarrow\mathbb{R}$,
$\left\langle\gamma_{1},\gamma_{2}\right\rangle$ denotes the multi-dimensional
inner product defined as
$\left\langle\gamma_{1},\gamma_{2}\right\rangle=\sum_{\mathbf{k}\in\mathbb{Z}^{D}}\gamma_{1}\left[\mathbf{k}\right]\mathtt{conj}(\gamma_{2}\left[\mathbf{k}\right])$.
The coefficients for the forward finite difference of order $N$ are denoted by
$\Delta^{N}\left[k\right]$. Specifically, it is defined as
$\Delta^{N}\left[k\right]=\Delta_{-}^{N}\left[-k\right]$, where
$\Delta^{N}_{-}$ is defined recursively as
$\Delta^{N+1}_{-}\left[k\right]=\Delta^{N}_{-}\ast\Delta^{1}_{-}\left[k\right]$,
where $\Delta^{1}_{-}\left[-1\right]=1,\Delta^{1}_{-}\left[0\right]=-1$,
$\Delta^{1}_{-}\left[k\right]=0,k\in\mathbb{Z}\setminus\left\\{-1,1\right\\}$.

We denote by $\mathbf{v}_{1},\dots,\mathbf{v}_{D}$ the set of vectors defining
a lattice
$\boldsymbol{\Lambda}=\left\\{\mathbf{VTk},\mathbf{k}\in\mathbb{Z}^{D}\right\\}$
where
$\mathbf{V}=\left[\mathbf{v}_{1},\dots,\mathbf{v}_{D}\right],\mathbf{v}_{d}\in\mathbb{R}^{D},\forall
d\in\left\\{1,\dots,D\right\\}$ and
$\mathbf{T}=\mathrm{diag}\left\\{T_{1},\dots,T_{D}\right\\}$, where $T_{d}$
denotes the sampling period across dimension $d\in\left\\{1,\dots,D\right\\}$.
Without loss of generality, we assume that $\mathbf{v}_{d}$ are _versors_ ,
i.e., $\|\mathbf{v}_{d}\|_{2}=1,\forall d\in\left\\{1,\dots,D\right\\}$. The
vectors are assumed linearly independent, and thus $\mathbf{V}$ is invertible.
Therefore, a function $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ can be
equivalently evaluated using Cartesian coordinates as
$f\left(\mathbf{x}_{c}\right)$ or lattice coordinates as
$f\left(\mathbf{Vx}_{v}\right)$ such that
$\mathbf{x}_{v}=\mathbf{V}^{-1}\mathbf{x}_{c}$. Unless specified otherwise,
$\gamma\left[\mathbf{k}\right],\mathbf{k}\in\mathbb{Z}^{D}$, denote the
samples of the input function $f$ on lattice $\boldsymbol{\Lambda}$ such that
$\gamma\left[\mathbf{k}\right]=f\left(\mathbf{VTk}\right)$. The dual lattice
$\widehat{\boldsymbol{\Lambda}}$ is defined as
$\widehat{\boldsymbol{\Lambda}}=\left\\{\widehat{\mathbf{V}}\mathbf{T}^{-1}\mathbf{k}|\mathbf{k}\in\mathbb{Z}^{D}\right\\}$,
where
$\widehat{\mathbf{V}}=\mathbf{V}^{-\top}=\left(\mathbf{V}^{-1}\right)^{\top}$
and $\mathbf{T}^{-1}=\mathrm{diag}\left\\{1/T_{1},\dots,1/T_{D}\right\\}$.

The $D$-dimensional Paley-Wiener space
$\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ of bandwidth
$\boldsymbol{\Omega}$ relative to lattice $\boldsymbol{\Lambda}$ consists of
functions $f\in L^{2}\left(\mathbb{R}^{D}\right)$ such that

$\mathrm{supp}\left(\mathcal{F}f\right)\subseteq\left\\{{\widehat{\mathbf{V}}}\boldsymbol{\omega}=\sum_{d=1}^{D}\omega_{d}\widehat{v}_{d}\in\mathbb{R}^{D}|\left\lvert\omega_{d}\right\rvert<\Omega_{d},\forall
d\in\left\\{1,\dots,D\right\\}\right\\}.$ (3)

For the one-dimensional case $D=1$, the lattice matrices
$\widehat{\mathbf{V}}$ reduce to the trivial case
$\widehat{\mathbf{V}}=\mathbf{V}=\left[1\right]$ and
$\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ reduces to the
classical one-dimensional Paley-Wiener space
$\mathsf{PW}_{\Omega}\left(\mathbb{R}\right)$.

For a random variable $\eta$ we denote by p.d.f. its probability density
function. We denote by $\mathcal{N}\left(\mu,\sigma^{2}\right)$ the normal
distribution with mean $\mu$ and standard deviation $\sigma$. A random
variable $\eta$ drawn from the Gaussian normal distribution is denoted as
$\eta\sim\mathcal{N}\left(\mu,\sigma^{2}\right)$.

##  2 Modulo acquisition and recovery

###  2.1 Recovery from one-dimensional modulo data

The centered modulo with threshold $\lambda$ is a function
$\mathscr{M}_{\lambda}:\mathbb{R}\rightarrow\mathbb{R}$ satisfying [12]

$\mathscr{M}_{\lambda}\left(x\right)=2\lambda\left({\left[\kern-3.00003pt\left[{\frac{x}{{2\lambda}}+\frac{1}{2}}\right]\kern-3.00003pt\right]-\frac{1}{2}}\right).$
(4)

When applied to a one-dimensional function $g$ the modulo non-linearity
generates values
$\mathscr{M}_{\lambda}\left(g\left(t\right)\right)\in\left[-\lambda,\lambda\right]$.

In analogy to Shannon-Nyquist sampling theory, the first recovery result in
the Unlimited Sensing Framework (USF) utilized bandlimited inputs, namely,
$g\in\mathsf{PW}_{\Omega}$. In the noiseless scenario, the _unlimited sampling
theorem_ [11, 12] guarantees that the input of the ideal modulo encoder can be
recovered from the output samples provided that the sampling period satisfies
$T<\frac{1}{2\Omega e}$. Furthermore, reconstruction is also possible in the
case of data corrupted by bounded noise if the following is true [12]

$\left(T\Omega e\right)^{N}g_{\infty}+2^{N}\eta_{\infty}<\lambda.$ (5)

The recovery approach used in [11, 12] aims to reconstruct the residual
function $\varepsilon_{g}\left(t\right)$ defined as
$\varepsilon_{g}\left(t\right)\triangleq
g\left(t\right)-\mathscr{M}_{\lambda}\left(g\left(t\right)\right)\in
2\lambda\cdot\mathbb{Z}$. In other words, for the ideal modulo encoder the
values of $\varepsilon_{g}\left(t\right)$ lie on an equally spaced grid with
step $2\lambda$. However, this is not true for non-ideal modulo encoders
exhibiting phenomena such as hysteresis, leading to reconstruction
distortions.

A generalized model of the modulo operator, called _modulo-hysteresis_ , was
introduced for the one-dimensional scenario [33, 34, 14]. Here, we generalize
this model to multi-dimensional sampling. As in the one-dimensional case, we
will show that modulo-hysteresis enables the separation of the folding times,
which will be used in the recovery in Section 3, 4, and 5. We begin with the
definition of the one-dimensional modulo-hysteresis.

######  Definition 1 (One-dimensional modulo-hysteresis).

The operator $\mathscr{M}_{\boldsymbol{\mathsf{H}}}$ with threshold $\lambda$
and hysteresis $h\in\left[0,2\lambda/3\right)$, where
$\mathbf{\mathsf{H}}=\left[\lambda\ h\right]$, generates a function
$z\left(t\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\left(t\right)$ for
input $g\in\mathsf{PW}_{\Omega}\left(\mathbb{R}\right)$, such that, for
$t\geqslant 0$

$z\left(t\right)=g\left(t\right)-\varepsilon_{g}\left(t\right),$ (6)

where

  * •

$\varepsilon_{g}\left(t\right)=2\lambda_{h}\sum\nolimits_{r=1}^{R}s_{r}\mathbbmtt{1}_{[\tau_{r},\infty)}\left(t\right)+hM,\quad$
$\lambda_{h}\triangleq\lambda-h/2$,

  * •

$M=\left\lfloor\tfrac{g\left(0\right)+\lambda}{h}\right\rfloor-1$,

  * •

$\tau_{r}$ and $s_{r}$ are the folding time and sign respectively, satisfying
$\tau_{0}=s_{0}=0$ and

$\displaystyle\tau_{1}$
$\displaystyle=\inf\left\\{t>\tau_{0}|\mathscr{M}_{\lambda}\left(g\left(t\right)+\lambda\right)=0\right\\},$
$\displaystyle s_{r}$
$\displaystyle=\mathrm{sign}\left(g\left(\tau_{r}\right)-g\left(\tau_{r-1}\right)\right),$
(7) $\displaystyle\tau_{r+1}$
$\displaystyle=\inf\left\\{t>\tau_{r}|\mathscr{M}_{\lambda}\left(g\left(t\right)-g\left(\tau_{r}\right)+hs_{r}\right)=0\right\\},\quad
r\geqslant 1.$

Furthermore, for $t<0$ we have
$z\left(t\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}\left[g\left(-\cdot\right)\right]\left(-t\right)$.
Let $\tau^{-}_{r},s^{-}_{r},r\geqslant 1$ be the sequence of folding times and
signs computed via (7) for
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}\left[g\left(-\cdot\right)\right]$. Then
we define $\tau_{r}\triangleq-\tau^{-}_{-r},s_{r}\triangleq
s^{-}_{-r},r\in\mathbb{Z},r<0$.

A key property of the 1D modulo-hysteresis is the folding time separation [14,
33]

$\tau_{r+1}-\tau_{r}\geqslant\frac{h}{\Omega\|g\|_{\infty}}.$ (8)

We note that the ideal modulo, which satisfies
$\mathscr{M}_{\lambda}\left(g\left(t\right)\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}g\left(t\right)$
for $h=0$ does not guarantee any separation via (8). The reconstruction
problem proposed aims to recover $g\left(kT\right)$ from
$y\left[k\right]=z\left(kT\right)$. Furthermore, it was shown that this
approach enables handling a number of modulo non-idealities [14, 35, 33, 24].

##  3 Multi-dimensional modulo sampling

###  3.1 Multi-dimensional lattice sampling preliminaries

Let $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ be a $D$-dimensional scalar
function. The data is then sampled on a lattice
$\boldsymbol{\Lambda}=\\{\mathbf{V}\mathbf{T}\mathbf{k}|\mathbf{k}\in\mathbb{Z}^{D}\\}$.
We denote the resulting samples by
$\gamma\left[\mathbf{k}\right]=f\left(\mathbf{VTk}\right)$. We assume that $f$
has a Fourier transform satisfying

$\mathrm{supp}\left(\mathcal{F}f\right)\subseteq\mathbb{D}.$

Upon sampling, the spectrum of $f$ is copied periodically, to produce the
multi-dimensional discrete-time Fourier transform
$F_{\boldsymbol{\Lambda}}\left(\boldsymbol{\omega}\right)=\sum_{\mathbf{k}\in\mathbb{Z}^{D}}f\left(\mathbf{VTk}\right)e^{-\jmath\left\langle\boldsymbol{\omega},\mathbf{k}\right\rangle}$,
whose support satisfies

$\mathrm{supp}\left(F_{\boldsymbol{\Lambda}}\right)=\mathbf{T}\mathbf{V}^{\top}\cdot\underset{\mathbf{n}\in\widehat{\boldsymbol{\Lambda}}}{\bigcup}\left(\mathbb{D}+2\pi\mathbf{n}\right).$
(9)

It was shown that $f$ can be recovered from its lattice samples
$\gamma\left[\mathbf{k}\right]$ if [36, 37]

$\mathbf{T}\mathbf{V}^{\top}\left(\mathbb{D}+2\pi\mathbf{n}_{1}\right)\cap\mathbf{T}\mathbf{V}^{\top}\left(\mathbb{D}+2\pi\mathbf{n}_{2}\right)=\emptyset,\quad\forall\mathbf{n_{1}},\mathbf{n_{2}}\in\widehat{\boldsymbol{\Lambda}},\mathbf{n_{1}}\neq\mathbf{n_{2}}.$
(10)

To ensure that recovery is possible, we assume that the spectrum of $f$ has a
compact support satisfying
$\mathbb{D}\subseteq\widehat{\mathbf{V}}\prod_{d=1}^{D}\left(-\Omega_{d},\Omega_{d}\right)$.
Formally, our assumption is
$f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$. Then $f$
can be reconstructed from samples $\gamma\left[\mathbf{k}\right]$ if (10) is
satisfied, which is sufficiently guaranteed if we replace $\mathbb{D}$ by
$\widehat{\mathbf{V}}\cdot\prod_{d=1}^{D}\left(-\Omega_{d},\Omega_{d}\right)$,
for which the terms on the left-hand-side of (10) are computed as

$\mathbf{T}\mathbf{V}^{\top}\left(\mathbb{D}+2\pi\mathbf{n}\right)=\mathbf{T}\mathbf{V}^{\top}\left[\mathbf{V}^{-\top}\prod_{d=1}^{D}\left(-\Omega_{d},\Omega_{d}\right)+2\pi\mathbf{n}\right],\quad\mathbf{n}\in\widehat{\boldsymbol{\Lambda}}.$
(11)

We use that
$\mathbf{n}=\mathbf{V}^{-\top}\mathbf{T}^{-1}\mathbf{k},\quad\mathbf{k}\in\mathbb{Z}^{D}$,
yielding

$\displaystyle\mathbf{T}\mathbf{V}^{\top}\left(\mathbb{D}+2\pi\mathbf{n}\right)=\mathbf{T}$
$\displaystyle\prod_{d=1}^{D}\left(-\Omega_{d},\Omega_{d}\right)+2\pi\mathbf{k},\quad\mathbf{k}\in\mathbb{Z}^{D}$
(12) $\displaystyle=$
$\displaystyle\prod_{d=1}^{D}\left(-T_{d}\Omega_{d},T_{d}\Omega_{d}\right)+2\pi\mathbf{k}.$
(13)

Therefore, (10) is true if

$T_{d}<\frac{\pi}{\Omega_{d}},\quad\forall d\in\left\\{1,\dots,D\right\\}.$
(14)

Just as in the one-dimensional case, considering the problem of sensor
saturation motivates using the concept of modulo folding also in the multi-
dimensional case. Next we go through some of the attempts to apply modulo for
multi-dimensional data.

###  3.2 Previous approaches for recovery from multi-dimensional data

As in the one-dimensional case, the problem with computing directly
$\gamma\left[\mathbf{k}\right]$ is that the sample values may be very large
which would saturate an analog-to-digital (ADC) acquisition device, which has
a restricted dynamic range [12]. The modulo operator was applied previously
for multi-dimensional inputs by processing a 1D slice at a time [19, 28].
However, these methods are not truly multi-dimensional because they don’t
exploit the multi-dimensional structure of the data. Furthermore, when dealing
with noise in 1D, the modulo samples
$y\left[\mathbf{k}\right]=\gamma\left[\mathbf{k}\right]+\eta\left[\mathbf{k}\right]-\varepsilon_{\gamma}\left[\mathbf{k}\right]$
require the separation of both residual $\varepsilon_{\gamma}$ and $\eta$
within the same dimension. This turns out to be contradictory, as detecting
$\varepsilon_{\gamma}$ requires a high-pass filter (such as $\Delta^{N}$ in
the case of USF), while denoising is typically done with low-pass filters
[34]. Furthermore, a denoising approach on modulo data was tested for multi-
dimensional signals [31]. However, we center our analysis on purely modulo
inversion techniques, where the noise sequence remains unaltered.

We define the following functions, representing slices of function
$f(\mathbf{V}\mathbf{x})$. Let
$f_{\mathbf{V}}\left(\mathbf{x}\right)\triangleq f(\mathbf{V}\mathbf{x})$. Let
$f_{\bar{\mathbf{x}}}:\mathbb{R}\rightarrow\mathbb{R}$ denote the slice along
dimension $x_{1}$ defined as $f_{\bar{\mathbf{x}}}\left(x\right)\triangleq
f\left(x\mathbf{v}_{1}+\sum_{d=2}^{D}x_{d}\mathbf{v}_{d}\right),\forall\bar{\mathbf{x}}\in\mathbb{R}^{D-1},\bar{\mathbf{x}}=\left[x_{2},\dots,x_{D}\right]$.
In the next proposition we also use a generic slices defined as the one-
dimensional function $g_{d}\left(x_{d}\right)\triangleq
f\left(\sum_{n=1}^{D}x_{n}\mathbf{v}_{n}\right)$ by fixing dimensions
$\left\\{x_{1},\dots,x_{d-1},x_{d+1},\dots,x_{D}\right\\},d\neq 1$. The
following proposition was proven in [28].

######  Proposition 1 (Bandlimited slices).

The function $f_{\bar{\mathbf{x}}}$ satisfies
$f_{\bar{\mathbf{x}}}\in\mathsf{PW}_{\Omega_{1}}\left(\mathbb{R}\right)$.
Furthermore, $g_{d}\in\mathsf{PW}_{\Omega_{d}}\left(\mathbb{R}\right)$,
$\forall d\in\left\\{1,\dots,D\right\\},d\neq 1$.

###### Proof.

The proof is in Section 7.1. ∎

The proposition above proves that any _slice_ of the multi-dimensional
function $f$ along lattice dimension $d$ has the spectrum compactly supported
within $\left(-\Omega_{d},\Omega_{d}\right)$. Furthermore, this implies that a
Bernšteĭn bound can be applied for each variable such that

$\displaystyle\begin{split}\left\lvert\frac{\partial}{\partial
x_{d}}f_{\mathbf{V}}\left(\mathbf{x}\right)\right\rvert&\leqslant\Omega_{d}\max_{x_{d}}\left\lvert
f_{\mathbf{V}}\left(\mathbf{x}\right)\right\rvert\leqslant\Omega_{d}\|f_{\mathbf{V}}\|_{\infty}\\\
&=\Omega_{d}\|f\|_{\infty},\quad\forall
d\in\left\\{1,\dots,D\right\\},\forall\bar{\mathbf{x}}\in\mathbb{R}^{D-1}.\end{split}$
(15)

The works in [19] and [28] apply one-dimensional ideal modulo to
$f_{\bar{\mathbf{x}}}$:

$\mathscr{M}_{\lambda}f_{\bar{\mathbf{x}}}\left(x\right)=f_{\bar{\mathbf{x}}}\left(x\right)-2\lambda_{h}\sum_{r\in\mathbb{Z}}s_{\bar{\mathbf{x}},r}\mathbbmtt{1}_{\left[\tau_{\bar{\mathbf{x}},r},\infty\right)}\left(x\right),x\in\mathbb{R}.$

Therefore we can define the "folded" multi-dimensional function
$z\left(\mathbf{x}\right)$ as

$z\left(\mathbf{Vx}\right)=\mathscr{M}_{\lambda}\left(f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right),$
(16)

where
$\mathbf{x}=\left[x_{1},\dots,x_{D}\right],\bar{\mathbf{x}}=\left[x_{2},\dots,x_{D}\right],x_{d}\in\mathbb{R},\forall
d\in\left\\{1,\dots,D\right\\}$. Subsequently, the output samples are
$z\left(\mathbf{VTk}\right)=\mathscr{M}_{\lambda}f_{\bar{\mathbf{T}}\bar{\mathbf{k}}}\left(k_{1}T_{1}\right)$,
where $\bar{\mathbf{T}}=\mathrm{diag}\left\\{T_{2},\dots,T_{D}\right\\}$.

Then, recovering $f_{\bar{\mathbf{T}}\bar{\mathbf{k}}}\left(k_{1}T_{1}\right)$
for all $\bar{\mathbf{T}}\bar{\mathbf{k}}$ from $z\left(\mathbf{VTk}\right)$
represents a line-by-line approach used in [19, 26, 28], which is guaranteed
to work if (5) holds true. However, as explained previously, this approach
does not exploit the multi-dimensional structure of the input data. This is
further motivated by the accepted knowledge in image processing that non-
separability in multiple dimensions has a lot more to offer than separability
[29, 30]. This motivates introducing a multi-dimensional modulo-hysteresis
model in the next section. We will show that the new model allows a
significantly large amount of noise, which is not possible with USF that
processes the data line-by-line.

###  3.3 Towards multi-dimensional modulo-hysteresis acquisition

Inspired from the 1D modulo-hysteresis operator that showed improvements for
noise robustness [34], we define in the following a new operator called multi-
dimensional modulo-hysteresis that addresses the issues discussed in the
previous subsection. The idea is to split the domain $\mathbb{R}^{D-1}$ in
disjoint sets confined in polytopes
$\mathcal{P}_{\bar{\mathbf{b}}},\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$ defined
as

$\mathcal{P}_{\bar{\mathbf{b}}}=\left\\{\bar{\mathbf{V}}B\left(\bar{\mathbf{b}}+\bar{\boldsymbol{\alpha}}\right)|\forall\bar{\boldsymbol{\alpha}}\in\left[0,1\right)^{D-1}\right\\},$
(17)

where $\bar{\mathbf{b}}=\left[b_{2},\dots,b_{D}\right]$ and
$B\in\mathbb{R}_{+}^{*}$ is the polytope edge length along directions parallel
with versors $\mathbf{v}_{2},\dots,\mathbf{v}_{D}$. The set
$\mathcal{P}_{\bar{\mathbf{b}}}$ is created via the last $D-1$ vectors of the
lattice basis $\mathbf{v}_{2},\dots,\mathbf{v}_{D}$ where the basis
coordinates lie in a set of rectangular polytopes
$\mathcal{R}_{\bar{\mathbf{b}}}$ such that
$\mathcal{P}_{\bar{\mathbf{b}}}=\left\\{\bar{\mathbf{V}}\bar{\mathbf{x}}|\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}\right\\}$,
where

$\mathcal{R}_{\bar{\mathbf{b}}}=\left\\{B\left(\bar{\mathbf{b}}+\bar{\boldsymbol{\alpha}}\right)|\bar{\boldsymbol{\alpha}}\in\left[0,1\right)^{D-1}\right\\}=\prod_{d=2}^{D}\left[b_{d}B,\left(b_{d}+1\right)B\right).$
(18)

Note that
$\mathbb{R}^{D-1}=\bigcup_{\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}}\mathcal{P}_{\bar{\mathbf{b}}}$.
Sets $\mathcal{P}_{\bar{\mathbf{b}}}$ thus can be used to split the domain of
function $f$ in disjoint _bands_ of width $B$ given by
$\mathcal{B}_{\bar{\mathbf{b}}}=\left(\mathbb{R}\mathbf{v}_{1}\right)\times\mathcal{P}_{\bar{\mathbf{b}}}$
identified using the indices in $\bar{\mathbf{b}}$, such that
$\mathbb{R}^{D}=\bigcup_{\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}}\mathcal{B}_{\bar{\mathbf{b}}}$.
By exploiting the smoothness of $f$ we can derive that, for a fixed
$x_{1}\in\mathbb{R}$, $f$ has bounded variation within each band, and thus the
folding can occur simultaneously on all coordinates $x_{2},\dots,x_{D}$, which
gives the folded signal a particular structure to be exploited in recovery.
The definition of the new operator is given as follows.

######  Definition 2 (Multi-dimensional modulo-hysteresis).

The operator $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$ with threshold
$\lambda$ and hysteresis $h\in\left[0,2\lambda/3\right)$, where
$\boldsymbol{\mathsf{H}}=\left[\lambda\ h\right]$, generates a function $z$
for input $f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$
such that, for $\left[\mathbf{x}\right]_{1}=x_{1}\geqslant 0$

$z\left(\mathbf{Vx}\right)=f\left(\mathbf{Vx}\right)-\varepsilon_{f}\left(\mathbf{Vx}\right),$
(19)

where $\varepsilon_{f}$ denotes the modulo-hysteresis residual defined as

$\varepsilon_{f}\left(\mathbf{Vx}\right)=h\left[M_{\bar{\mathbf{b}}}+\sum_{r=0}^{R_{\bar{\mathbf{b}}}^{+}}s_{\bar{\mathbf{b}},r}\mathbbmtt{1}_{\left[\tau_{\bar{\mathbf{b}},r},\infty\right)}\left(x_{1}\right)\right],\forall\bar{\mathbf{x}}=\left[x_{2},\dots,x_{D}\right]\in\mathcal{R}_{\bar{\mathbf{b}}},$
(20)

where $\mathcal{R}_{\bar{\mathbf{b}}}$ satisfies (18),
$\mathbf{x}=\left[x_{1},\dots,x_{D}\right]$,
$\bar{\mathbf{x}}=\left[x_{2},\dots,x_{D}\right]$, and
$M_{\bar{\mathbf{b}}}\in\mathbb{Z}$ satisfies

$M_{\bar{\mathbf{b}}}=\left\lfloor\frac{\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(0\right)+\lambda}{h}\right\rfloor-1,$
(21)

and
$\tau_{\bar{\mathbf{b}},r}\in\left\\{0,\dots,R_{\bar{\mathbf{b}}}^{+}\right\\},s_{\bar{\mathbf{b}},r}\in\left\\{-1,0,1\right\\}$
are the folding times and signs in band $\mathcal{B}_{\bar{\mathbf{b}}}$,
respectively, defined as

$\displaystyle\tau_{\bar{\mathbf{b}},r+1}=\inf\left\\{x_{1}>\tau_{\bar{\mathbf{b}},r}|\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\varepsilon_{\bar{\mathbf{b}},r}\left(x_{1}\right)\right\rvert=\lambda\right\\},$
(22) $\displaystyle
s_{\bar{\mathbf{b}},r+1}=\mathrm{sign}\left[f_{B\bar{\mathbf{b}}}\left(\tau_{\bar{\mathbf{b}},r+1}\right)-\varepsilon_{\bar{\mathbf{b}},r}\left(\tau_{\bar{\mathbf{b}},r+1}\right)\right],$
(23)
$\displaystyle\varepsilon_{\bar{\mathbf{b}},r+1}\left(x_{1}\right)=\varepsilon_{\bar{\mathbf{b}},r}\left(x_{1}\right)+hs_{\bar{\mathbf{b}},r+1}\mathbbmtt{1}_{\left[\tau_{\bar{\mathbf{b}},r+1},\infty\right)}\left(x_{1}\right),$
(24)

where $r\in\left\\{0,\dots,R_{\bar{\mathbf{b}}}^{+}-1\right\\}$,
$\varepsilon_{\bar{\mathbf{b}},r}$ is a recursive sequence of functions for
computing the residual $\varepsilon_{f}$ such that
$\varepsilon_{\bar{\mathbf{b}},0}\left(x_{1}\right)=hM_{\bar{\mathbf{b}}}$,
$\varepsilon_{\bar{\mathbf{b}},R_{\bar{\mathbf{b}}}^{+}}\left(x_{1}\right)=\varepsilon_{f}\left(\mathbf{Vx}\right),\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}$
and $\tau_{\bar{\mathbf{b}},0}=s_{\bar{\mathbf{b}},0}=0$. Furthermore
$R_{\bar{\mathbf{b}}}^{+}\in\mathbb{Z}_{+}$ satisfies

$\left\\{x_{1}>\tau_{\bar{\mathbf{b}},R_{\bar{\mathbf{b}}}^{+}}|\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\varepsilon_{\bar{\mathbf{b}},R_{\bar{\mathbf{b}}}^{+}}\left(x_{1}\right)\right\rvert=\lambda\right\\}=\emptyset.$
(25)

Furthermore, for $x_{1}<0$ we have

$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{Vx}\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}\left[f^{-}\right]\left(-x_{1}\mathbf{v}_{1}+x_{2}\mathbf{v}_{2}+\dots+x_{D}\mathbf{v}_{D}\right),$

where $f^{-}\left(\mathbf{x}\right)=f_{\bar{\mathbf{x}}}\left(-x_{1}\right)$.
Let $\tau_{\bar{\mathbf{b}},r}^{-},s_{\bar{\mathbf{b}},r}^{-},r\geqslant 1$ be
the folding times and signs computed via (22), (23) for
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}\left[f^{-}\right]$. Then we define
$\tau_{r}\triangleq-\tau^{-}_{-r},s_{r}\triangleq
s^{-}_{-r},r\in\mathbb{Z},r<0$.

We note that the operator in Definition 2 is backwards compatible with the
one-dimensional operator in Definition 1. Specifically, if one chooses
$\mathcal{R}_{\bar{\mathbf{b}}}$ to be a single point
$\mathcal{R}_{\bar{\mathbf{b}}}=B\bar{\mathbf{b}}\in\mathbb{R}^{D-1}$ instead
of a hypercube, then $\tau_{\bar{\mathbf{b}},r}$ and $s_{\bar{\mathbf{b}},r}$
in (22) are the same as $\tau_{r}$ in Definition 1. Furthermore, it was shown
that, for $h=0$, the one-dimensional modulo-hysteresis operator in Definition
1 is identical to an ideal modulo operator (4) [14].

---
|  |
---|---
|
Figure 2: A random two-dimensional bandlimited input
$f\left(x_{1},x_{2}\right)$ was generated (a). The ideal modulo output
$\mathscr{M}_{\lambda}f$ is in (b) and the generalized modulo output
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f$ in (c). The corresponding
residual functions are depicted in (d) for ideal modulo and in (e) for
generalized modulo. For $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$, the
lattice $\boldsymbol{\Lambda}$ consists of vectors
$\mathbf{v}_{1}=\left[0.97,0.25\right]^{\top},\mathbf{v}_{2}=\left[0.32,0.95\right]^{\top}$.
The folding curves of the ideal modulo – the contours in (b) & (d) – are
unknown _a priori_. In the case of the modulo-hysteresis, folding occurs along
straight lines with directions dictated by the lattice $\boldsymbol{\Lambda}$
which is known _a priori_. This property will be exploited in recovery.

###  3.4 Properties of the proposed operator

In the following we give a number of properties of the multi-dimensional
modulo-hysteresis operator for a bandlimited input.

######  Proposition 2 (Folding time separation).

Assume that $\tau_{\bar{\mathbf{b}},r}$ are well-defined in (22) for
$f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ and
$r\in\left\\{1,\dots,R\right\\}$ where $R\geqslant 1$ and
$\tau_{\bar{\mathbf{b}},0}=0$. Furthermore, assume that
$\mathcal{D}f<\min\left\\{h/2,2\lambda-3h\right\\}$, where

$\mathcal{D}f\triangleq\mathop{\sup_{\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}}}_{x_{1}\in\mathbb{R}}\left[\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right].$
(26)

Then

$\tau_{\bar{\mathbf{b}},r+1}-\tau_{\bar{\mathbf{b}},r}\geqslant\frac{h}{\Omega_{1}\|f\|_{\infty}},\quad\forall\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}.$
(27)

###### Proof.

The proof is in Section 7.1. ∎

######  Proposition 3 (Well-defined operator).

For input $f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$
operator $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$ in Definition 1 is _well-
defined_ if $\mathcal{D}f<\min\left\\{h/2,2\lambda-3h\right\\}$, where
$\mathcal{D}f$ satisfies (26).

###### Proof.

The proof is in Section 7.1. ∎

######  Proposition 4 (Modulo output dynamic range).

Let $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$ be the multi-dimensional
modulo-hysteresis operator in Definition 2 and
$f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ such that
$\mathcal{D}f$ (26) satisfies
$\mathcal{D}f<\min\left\\{h/2,2\lambda-3h\right\\}$. Then
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{x}\right)\in\left[-\lambda,\lambda\right],\forall\mathbf{x}\in\mathbb{R}^{D}$.

###### Proof.

The proof is in Section 7.1. ∎

Proposition 4 shows that operator $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$
has a similar effect as the one-dimensional modulo nonlinearity, in that it
keeps a signal within a fixed dynamic range $\left[-\lambda,\lambda\right]$.

######  Corollary 1 (Bound for intra-band variation).

The quantity $\mathcal{D}f$ defined in Proposition 27 can be bounded as

$\mathcal{D}f\leqslant\|f\|_{\infty}\cdot
B\sqrt{D}\cdot\|\boldsymbol{\Omega}\|_{2}.$ (28)

###### Proof.

The proof is in Section 7.1. ∎

In Fig. 2 the variation of $M_{\bar{\mathbf{b}}}$, which can be seen as folds
along coordinates $x_{2},\dots,x_{D}$ for $x_{1}=0$, is _gradual_ , meaning
that $M_{\bar{\mathbf{b}}}$ changes by $1$ between neighboring bands.
Formally, we define by $\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$ the
set comprising all neighboring bands of band $\bar{\mathbf{b}}$ below.

######  Definition 3 (Neighboring bands).

The set $\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$ of vectors
neighboring $\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$ is defined as the set of all
$\bar{\mathbf{b}}^{*}\in\mathbb{Z}^{D-1}$ for which $\exists
d^{*}\in\left\\{2,\dots,D\right\\}$ such that
$\left\lvert\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}-\left[\bar{\mathbf{b}}\right]_{d^{*}}\right\rvert=1$
and
$\left[\bar{\mathbf{b}}\right]_{d}=\left[\bar{\mathbf{b}}^{*}\right]_{d},\forall
d\in\left\\{2,\dots,D\right\\}\setminus d^{*}$.

In the following, we provide conditions for which
$M_{\bar{\mathbf{b}}^{*}}-M_{\bar{\mathbf{b}}}\in\left\\{-1,0,1\right\\}$
where
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$.

######  Proposition 5 (Variation of $M_{\bar{\mathbf{b}}}$).

For $\forall\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$, let
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$. Let
$M_{\bar{\mathbf{b}}^{*}}$ and $M_{\bar{\mathbf{b}}}$ be the modulo-hysteresis
constants for an input
$f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ satisfying
the following condition as per Definition 2

$\|f\|_{\infty}B\sqrt{D}\cdot\sqrt{\sum_{d=1}^{D}\Omega_{d}^{2}}<\min\left\\{\frac{h}{2},2\lambda-3h\right\\}.$
(29)

Then
$M_{\bar{\mathbf{b}}^{*}}-M_{\bar{\mathbf{b}}}\in\left\\{-1,0,1\right\\}$.

###### Proof.

We first note that sets $\mathcal{R}_{\bar{\mathbf{b}}}$ and
$\mathcal{R}_{\bar{\mathbf{b}}^{*}}$ are neighboring polytopes which, via
their definition, satisfy
$\mathrm{cl}\left(\mathcal{R}_{\bar{\mathbf{b}}}\right)\cap\mathrm{cl}\left(\mathcal{R}_{\bar{\mathbf{b}}^{*}}\right)\neq\emptyset$.
Using this in conjunction with the properties of supremum and infimum, we get

$\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(0\right)=\max_{\bar{\mathbf{x}}\in{\mathrm{cl}\left(\mathcal{R}_{\bar{\mathbf{b}}}\right)}}f_{\bar{\mathbf{x}}}\left(0\right)\geqslant\min_{\bar{\mathbf{x}}\in{\mathrm{cl}\left(\mathcal{R}_{\bar{\mathbf{b}}^{*}}\right)}}f_{\bar{\mathbf{x}}}\left(0\right)=\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}^{*}}}f_{\bar{\mathbf{x}}}\left(0\right).$
(30)

Condition (29) implies $\mathcal{D}f<h/2$ due to Corollary 1. Therefore,

$\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(0\right)-\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(0\right)<h/2.$
(31)

Using (30) and (31)

$\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}^{*}}}f_{\bar{\mathbf{x}}}\left(0\right)-\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(0\right)<h/2.$
(32)

Given that, by definition,
$M_{\bar{\mathbf{b}}^{*}}=\left\lfloor\frac{\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}^{*}}}f_{\bar{\mathbf{x}}}\left(0\right)+\lambda}{h}\right\rfloor-1$,
it can be shown by direct derivation that $M_{\bar{\mathbf{b}}^{*}}\leqslant
M_{\bar{\mathbf{b}}}+1$. Furthermore, by swapping $\bar{\mathbf{b}}$ and
$\bar{\mathbf{b}}^{*}$ in the derivation above we get

$M_{\bar{\mathbf{b}}^{*}}\in\left\\{M_{\bar{\mathbf{b}}}-1,M_{\bar{\mathbf{b}}},M_{\bar{\mathbf{b}}}+1\right\\}.$
(33)

∎

Therefore, using Definition 2, for $x_{1}=0$, in neighboring bands
$\bar{\mathbf{b}}$ and $\bar{\mathbf{b}}^{*}$ the residual $\varepsilon_{f}$
is either the same, or differs by $h$. This is similar to the behavior of the
residual around a folding time $\tau_{\bar{\mathbf{b}},r}$ along dimension
$x_{1}$.

###  3.5 Problem formulation

The measurements $y\left[\mathbf{k}\right]$ are assumed to be samples on a
multi-dimensional lattice
$\boldsymbol{\Lambda}=\left\\{\mathbf{VTk}|\mathbf{k}\in\mathbb{Z}^{D}\right\\}$,
such that

$y\left[\mathbf{k}\right]=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f(\mathbf{VTk})+\eta\left[\mathbf{k}\right]=\gamma\left[\mathbf{k}\right]-\varepsilon_{\gamma}\left[\mathbf{k}\right]+\eta\left[\mathbf{k}\right],$
(34)

where $\gamma\left[\mathbf{k}\right]=f\left(\mathbf{VTk}\right)$,
$\varepsilon_{\gamma}\left[\mathbf{k}\right]=\varepsilon_{f}\left(\mathbf{VTk}\right)$
and $\eta\left[\mathbf{k}\right]\sim\mathcal{N}\left(0,\sigma^{2}\right)$. The
known variables are the input bandwidth $\boldsymbol{\Omega}$, number of
dimensions $D$, the lattice $\boldsymbol{\Lambda}$, modulo-hysteresis
parameters $B,\lambda,h$ and output samples $y\left[\mathbf{k}\right]$. The
proposed reconstruction problem is to compute the input lattice samples
$\widetilde{\gamma}\left[\mathbf{k}\right]$ defined as

$\widetilde{\gamma}\left[\mathbf{k}\right]=\gamma\left[\mathbf{k}\right]+\eta\left[\mathbf{k}\right]+Mh,$
(35)

where $M\in\mathbb{Z}$ is an unknown integer. The input
$\gamma\left[\mathbf{k}\right]$ can only be reconstructed up to an integer
multiple of $h$ given that
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{x}\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}\left[f+Mh\right]\left(\mathbf{x}\right),$
$\forall M\in\mathbb{Z}$ (19,20).

##  4 Detecting modulo-hysteresis discontinuities

We define $N_{d}^{B}=\frac{B}{T_{d}}$. For simplicity, we assume that
$N_{d}^{B}\in\mathbb{Z},N_{d}^{B}\geqslant 1,\forall d\in\left\\{2,\dots
D\right\\}$, i.e., for a polytope $\mathcal{P}_{\bar{\mathbf{b}}}$ we have an
integer number of sampling periods $T_{d}$ along each of its edges of length
$B$ and direction $\mathbf{v}_{d}$. The recovery is achieved in several steps

  1. 1.

Compute $M_{\bar{\mathbf{b}}}$, the folding times $\tau_{\bar{\mathbf{b}},r}$
and signs $s_{\bar{\mathbf{b}},r}$ in each band
$\mathcal{B}_{\bar{\mathbf{b}}}$.

  2. 2.

Compute residual $\varepsilon_{f}\left(\mathbf{Vx}\right)$ (20).

  3. 3.

Compute the samples $\widetilde{\gamma}\left[\mathbf{k}\right]$.

For step 1, we define a filter $\psi$ and compute $\left\langle
y,\psi\right\rangle$. For detecting the folding times and signs, we choose
$\psi=\psi_{\bar{\mathbf{b}},m}$, which is centered in sample $m$ along
dimension $x_{1}$ in band $\bar{\mathbf{b}}$ as

$\psi_{\bar{\mathbf{b}},m}\left[\mathbf{k}\right]=\left\\{\begin{array}[]{cc}\Delta^{N}\left[k_{1}-m\right]\cdot\frac{1}{N^{B}},&\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}},\\\
0,&\bar{\mathbf{T}}\bar{\mathbf{k}}\notin\mathcal{R}_{\bar{\mathbf{b}}},\end{array}\right.$
(36)

where $\bar{\mathbf{k}}\in\mathbb{Z}^{D-1}$ and
$N^{B}=\prod_{d=2}^{D}N_{d}^{B}$ is the number of samples in each band
$\bar{\mathbf{b}}$ for $k_{1}$ fixed. For fixed
$\bar{\mathbf{k}}\in\mathbb{Z}^{D-1}$,
$\psi_{\bar{\mathbf{b}},m}\left[\mathbf{k}\right]$ is a finite difference
filter along dimension $k_{1}$.

For detecting the $M_{\bar{\mathbf{b}}}$, we need a filter detecting the
change in $M_{\bar{\mathbf{b}}}$ for two bands $\bar{\mathbf{b}}$,
$\bar{\mathbf{b}}^{*}$, such that
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$ and
$\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}=\left[\bar{\mathbf{b}}\right]_{d^{*}}+1$.
We define
$\psi=\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\left[\mathbf{k}\right]$ as

$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\left[\mathbf{k}\right]=\left\\{\begin{array}[]{cc}\Delta^{N}\left[k_{d^{*}}-k_{\bar{\mathbf{b}}}\right]\cdot\frac{1}{N^{B*}},&\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}\cup\mathcal{R}_{\bar{\mathbf{b}}^{*}},k_{1}=0,\\\
0,&\quad\mathrm{otherwise},\end{array}\right.$ (37)

where
$k_{\bar{\mathbf{b}}}=N_{d^{*}}^{B}\cdot\left[\bar{\mathbf{b}}\right]_{d^{*}}-1$
and $N^{B*}=\prod_{d=2,d\neq d^{*}}^{D}N_{d}^{B}$ is the number of samples in
each band $\bar{\mathbf{b}}$ for $k_{1}$, $k_{d^{*}}$ fixed. Therefore,
similar to $\psi_{\bar{\mathbf{b}},m}$, filter
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ is a finite difference filter
along dimension $x_{d^{*}}$ which is perpendicular to the hyperplane
separating bands $\bar{\mathbf{b}}$ and $\bar{\mathbf{b}}^{*}$. Furthermore,
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ is constant within each band
$\bar{\mathbf{b}}$ and $\bar{\mathbf{b}}^{*}$. This means that the finite
difference filter is repeated $N^{B*}$ times across dimensions
$d\in\left\\{2,\dots,D\right\\}\setminus d^{*}$, which has a noise averaging
effect.

###  4.1 Detecting folding times and signs

To detect $\tau_{\bar{\mathbf{b}},r}$ and $s_{\bar{\mathbf{b}},r}$, we use
filter $\psi_{\bar{\mathbf{b}},m}$ to compute sequence
$y_{\bar{\mathbf{b}},m}$

$y_{\bar{\mathbf{b}},m}=\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle=\underbrace{\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle}_{\textsf{Input}}-\underbrace{\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle}_{\textsf{Residual}}+\underbrace{\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle}_{\textsf{Noise}}.$
(38)

In (38) the filtered samples $y_{\bar{\mathbf{b}},m}$ are composed of three
terms: the _input term_
$\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle$, the _residual
term_
$\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle$ and
the _noise term_ $\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle$.
Given that all three are unknown, the general recovery strategy is to separate
them via thresholding; as will be shown later, thresholding samples
$y_{\bar{\mathbf{b}},m}$ allows to compute the folding times and signs. While
this will be derived rigorously later in propositions 6 and 7, here we give a
brief intuitive explanation of the recovery method, by explaining the effect
that filter $\psi_{\bar{\mathbf{b}},m}$ has on all $3$ terms in the right-
hand-side of (38). A similar analysis applies to filter
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ which will be described in
Section 4.2.

As noted before, $\psi_{\bar{\mathbf{b}},m}$ is a finite difference filter
along dimension $x_{1}$. It was shown for the one-dimensional case that this
causes $\gamma\left[\mathbf{k}\right]$ to vanish for large $N$. Furthermore,
it generates peaks at the folding times in residual $\varepsilon_{\gamma}$ and
also amplifies the noise $\eta$ [12, 14]. This latter effect is undesirable
for recovery. To decrease the effect of the noise the finite difference filter
was convolved with a spline in the one-dimensional case, but this also makes
the detection of $\varepsilon_{\gamma}$ more difficult [34]. Here we can
address noise filtering without affecting the folding time detection by
exploiting the multi-dimensional structure of
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$.

For fixed $k_{1}\in\mathbb{Z}$, the filter
$\psi_{\bar{\mathbf{b}},m}\left[\mathbf{k}\right]$ is constant along
dimensions $k_{2},\dots,k_{D}$ as long as
$\bar{\mathbf{V}}\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}$
and thus the inner product has an averaging effect. However, we know that the
modulo residual corresponding to
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{Vx}\right),\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}$
is also constant within $\mathcal{R}_{\bar{\mathbf{b}}}$ by definition, and
therefore the averaging effect does not affect the residual edges, which are
along dimension $k_{1}$. Moreover, given that $f$ is smooth and changes slowly
within a band $\mathcal{B}_{\bar{\mathbf{b}}}$, the filter averaging along
dimensions $\bar{\mathbf{k}}$ has very little effect on $\gamma$. Therefore,
along dimensions $\bar{\mathbf{k}}$, the filter acts mainly on the noise
sequence $\eta\left[\bar{\mathbf{k}}\right]$ by narrowing its p.d.f. around
the origin such that its effect gradually vanishes.

As in the one-dimensional case, the modulo output $z$ is smooth in-between the
folds, and has discontinuities at the folding times. The filter
$\psi_{\bar{\mathbf{b}},m}$ responds with pulses of non-zero support to the
input discontinuities. To account for this, we define by $\mathbb{S}_{N}$ the
support of the filtered residual such that (see [14] for details)

$\mathbb{S}_{N}=\mathrm{supp}\left[\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle\right]=\underset{r\in\mathbb{Z}^{*}}{\bigcup}\left\\{\left\lfloor\frac{\tau_{\bar{\mathbf{b}},r}}{T_{1}}\right\rfloor-N,\dots,\left\lfloor\frac{\tau_{\bar{\mathbf{b}},r}}{T_{1}}\right\rfloor\right\\}.$

The following theorem shows how $\mathbb{S}_{N}$ can be recovered by
thresholding sequence $\left\langle y,\psi_{\bar{\mathbf{b}},m}\right\rangle$,
which is an important step in computing folding times
$\tau_{\bar{\mathbf{b}},r}$.

######  Proposition 6 (Detection of folding times).

Let $f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ and let
$z\left(\mathbf{x}\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{x}\right)$
be the output of a multi-dimensional modulo-hysteresis model with parameters
$\lambda,h,B$. Furthermore, let
$y\left[\mathbf{k}\right]=z\left(\mathbf{VTk}\right)+\eta\left[\mathbf{k}\right]$
be the samples of the modulo output computed on lattice $\boldsymbol{\Lambda}$
corrupted by a noise sequence
$\eta\left[\mathbf{k}\right]\sim\mathcal{N}\left(0,\sigma^{2}\right)$.
Furthermore, assume that $T_{d}<B,\forall d\in\left\\{2,\dots,D\right\\}$ and
that

$\displaystyle\begin{split}\|f\|_{\infty}B\sqrt{D}\cdot\sqrt{\sum_{d=1}^{D}\Omega_{d}^{2}}&<\min\left\\{\frac{h}{2},2\lambda-3h\right\\},\\\
\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}&<h/2.\end{split}$ (39)

If $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant h/2$ then
$m\in\mathbb{S}_{N}$ with probability $p>1-p_{\mathsf{err}}$ where

$p_{\mathsf{err}}\leqslant e^{-C^{2}},\quad\mathrm{where}\
C=\frac{h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}}{\sigma\sqrt{2^{N+1}}}\prod_{d=2}^{D}\sqrt{\frac{B}{T_{d}}}.$
(40)

###### Proof.

The proof is in Section 7.2. ∎

Due to Proposition 6, for
$\forall\sigma,\lambda,h,\boldsymbol{\Omega}\in\mathbb{R}^{D}_{+}$ satisfying
(39) and a fixed $m\in\mathbb{Z}$, one can choose $T_{1},\dots,T_{D}>0$ such
that the truth value of $m\in\mathbb{S}_{N}$ is evaluated correctly with an
arbitrarily large probability. We note that $p_{\mathsf{err}}$ measures the
probability when a recovery error is possible, but not guaranteed, therefore
the error probability is smaller in a real scenario. A small error in
Proposition 6 means a large $C$, which can be achieved by decreasing the
sampling periods $T_{1},\dots,T_{d}$ or increasing $B,h$ or number of
dimensions $D$.

The residual $\varepsilon_{\gamma}$, used for reconstructing $\gamma$,
requires detecting constants $M_{\bar{\mathbf{b}}}$ in addition to the folding
times ${\tau}_{\bar{\mathbf{b}},r}$ and signs ${s}_{\bar{\mathbf{b}},r}$, as
will be explained in the next subsection.

###  4.2 Detecting constants $M_{\bar{\mathbf{b}}}$

Given that we can only recover $\gamma$ up to an integer multiple of $h$ (35),
we define $\widetilde{M}_{\bar{\mathbf{b}}}\triangleq
M_{\bar{\mathbf{b}}}-M_{\bar{\mathbf{0}}}$, where $\bar{\mathbf{0}}$ is the
null vector of $\mathbb{Z}^{D-1}$, and recover
$\widetilde{M}_{\bar{\mathbf{b}}}$. Just as the folding times, different
values of $M_{\bar{\mathbf{b}}}$ for adjacent bands cause discontinuities.
However, unlike the detection of the folding times, here we have additional
information. Specifically, we know that the discontinuities may only be
located at the neighboring sides of polytopes
$\mathcal{P}_{\bar{\mathbf{b}}}$. We use the filter
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ (37) to detect the
discontinuities in a similar fashion to detecting the folding times via
$\psi_{\bar{\mathbf{b}},m}$. This time, however, the finite differences
$\Delta^{N}$ computed via $\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$
evaluate variations across dimensions $x_{2},\dots,x_{D}$.

######  Proposition 7 (Detection of constants $M_{\bar{\mathbf{b}}}$).

Let
$y\left[\mathbf{k}\right]=z\left(\mathbf{VTk}\right)+\eta\left[\mathbf{k}\right]$,
where $\eta\left[\mathbf{k}\right]\sim\mathcal{N}\left(0,\sigma^{2}\right)$
and $z\left(\mathbf{x}\right)$ is the output of a modulo-hysteresis operator
$z\left(\mathbf{x}\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{x}\right)$
with parameters $\lambda,h$ and
$f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ satisfying

$\|f\|_{\infty}B\sqrt{D}\cdot\sqrt{\sum_{d=1}^{D}\Omega_{d}^{2}}<\min\left\\{h/2,2\lambda-3h\right\\},$
(41)

Furthermore, let $\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$,
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$ and
$d^{*}\in\left\\{2,\dots,D\right\\}$ such that
$\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}=\left[\bar{\mathbf{b}}\right]_{d^{*}}+1$.
Assume that

$\displaystyle\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}<h/2,$
(42) $\displaystyle\left(N+1\right)T_{d^{*}}<B.$ (43)

Then the following is true with probability $p>1-p_{\mathsf{err}}$

$\begin{cases}M_{\bar{\mathbf{b}}^{*}}=M_{\bar{\mathbf{b}}}+\mathrm{sign}\left(\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)&\mathrm{if}\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\geqslant
h/2,\\\
M_{\bar{\mathbf{b}}^{*}}=M_{\bar{\mathbf{b}}}&\mathrm{otherwise}.\end{cases}$
(44)

where

$p_{\mathsf{err}}\leqslant e^{-\kappa^{2}},\quad\mathrm{where}\
\kappa=\frac{h/2-\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}}{\sigma\sqrt{2^{N+1}}}\prod_{d=2,d\neq
d^{*}}^{D}\sqrt{\frac{B}{T_{d}}}.$ (45)

###### Proof.

The proof is in Section 7.2. ∎

##  5 Input reconstruction

###  5.1 Recovery with the proposed operator

We begin with the noiseless input recovery scenario $\sigma=0$ where, via
Proposition 6, the set $\mathbb{S}_{N}$ is perfectly identified with
probability $1$. Furthermore, a constant
$\widetilde{M}_{\bar{\mathbf{b}}}=M_{\bar{\mathbf{b}}}-M_{\bar{\mathbf{0}}}$
can be perfectly recovered from $\widetilde{M}_{\bar{\mathbf{b}}^{*}}$ where
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$
according to Proposition 7. The following theorem proves the input recovery
conditions in the case $\sigma=0$.

######  Theorem 2 (Noiseless input reconstruction).

Let $f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ and let
$z\left(\mathbf{x}\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{x}\right)$
be the output of a multi-dimensional modulo-hysteresis model with parameters
$\lambda,h,B$. Furthermore, let
$y\left[\mathbf{k}\right]=z\left(\mathbf{VTk}\right)$ be the samples of the
modulo output computed on lattice $\boldsymbol{\Lambda}$. Furthermore, for
$N\in\mathbb{Z},N\geqslant 1$, assume that

$\displaystyle\|f\|_{\infty}B\sqrt{D}\cdot\|\boldsymbol{\Omega}\|_{2}<\min\left\\{h/2,2\lambda-3h\right\\},$
(46)
$\displaystyle\left(T_{d}\Omega_{d}e\right)^{N}\|f\|_{\infty}<h/2,\quad\forall
d\in\left\\{1,\dots,D\right\\},$ (47)
$\displaystyle\left(N+1\right)T_{1}<\frac{h}{\Omega_{1}\|f\|_{\infty}},$ (48)
$\displaystyle\left(N+1\right)T_{d}<B,\quad\forall
d\in\left\\{2,\dots,D\right\\}.$ (49)

Then samples
$\widetilde{\gamma}\left[\mathbf{k}\right]=\gamma\left[\mathbf{k}\right]-hM_{\bar{\mathbf{0}}}$
can be perfectly reconstructed from $y\left[\mathbf{k}\right]$.

###### Proof.

The proof is in Section 7.2. ∎

The interpretation of the sufficient conditions in Theorem 2 is as follows.
The modulo-hysteresis is well-defined due to a bounded intra-band variation
guaranteed by (46). Condition (47) bounds the $N$-th order difference of the
input along all of the dimensions, ensuring that the filter has enough
shrinking effect on the input. Finally, (48) and (49) guarantee enough samples
in between the folds (48) and within each band (49) so that the supports of
the filters detecting consecutive discontinuities don’t overlap.

In the general case where $\sigma>0$ the following result holds true.

######  Theorem 3 (Noisy input reconstruction).

Let $f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ and
$z\left(\mathbf{x}\right)=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}f\left(\mathbf{x}\right)$.
Assume that
$y\left[\mathbf{k}\right]=z\left(\mathbf{VTk}\right)+\eta\left[\mathbf{k}\right]$
are known for $k_{d}\in\left\\{1,\dots,K^{d}_{\mathsf{max}}\right\\}$ where
$\eta\left[\mathbf{k}\right]\sim\mathcal{N}\left(0,\sigma^{2}\right)$, such
that
$B_{\mathsf{max}}^{d}\triangleq\frac{K^{d}_{\mathsf{max}}}{N_{d}^{B}}\in\mathbb{Z}$
for $d\in\left\\{2,\dots,D\right\\}$. Then, if (46-49) are true, then the
input samples $\widetilde{\gamma}\left[\mathbf{k}\right]$ can be recovered
from $y\left[\mathbf{k}\right]$ with a probability $p>p_{\mathsf{acc}}$ such
that

$p_{\mathsf{acc}}\geqslant\left(1-e^{-C^{2}}\right)^{K_{\mathsf{max}}^{1}\cdot\prod_{d=2}^{D}B_{\mathsf{max}}^{d}}\cdot\left(1-e^{-\kappa_{\mathsf{min}}^{2}}\right)^{\prod_{d=2}^{D}B_{\mathsf{max}}^{d}},$

where

$\displaystyle C$
$\displaystyle=\frac{h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}}{\sigma\sqrt{2^{N+1}}}\prod_{d=2}^{D}\sqrt{\frac{B}{T_{d}}},$
$\displaystyle\kappa_{\mathsf{min}}$
$\displaystyle=\frac{h/2-\left[\left(\max_{d}T_{d}\Omega_{d}\right)\cdot
e\right]^{N}\|f\|_{\infty}}{\sigma\sqrt{2^{N+1}}}\sqrt{\frac{B^{D-2}}{T_{\mathsf{max}}^{D-2}}},$

where $T_{\mathsf{max}}=\max_{d\in\left\\{2,\dots,D\right\\}}T_{d}$.

###### Proof.

Theorem 2 assumes that Proposition 6 and 7 hold with $p_{\mathsf{err}}=0$ for
all filters $\psi_{\bar{\mathbf{b}},m}$ and
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$. To calculate the overall error
probability when this assumption is not true, we count the filters above, when
used in reconstruction, as follows. There are a total of
$\prod_{d=2}^{D}B_{\mathsf{max}}^{d}$ bands, and $K_{\mathsf{max}}^{1}$
samples along dimension $x_{1}$. Then the probability that Proposition 6 holds
for all filters $\psi_{\bar{\mathbf{b}},m}$ is
$\left(1-e^{-C^{2}}\right)^{K_{\mathsf{max}}^{1}\cdot\prod_{d=2}^{D}B_{\mathsf{max}}^{d}}$.

Next, in the case of Proposition 7, we bound the error probability as follows

$\displaystyle p_{\mathsf{err}}\leqslant e^{-\kappa^{2}}\leqslant
e^{-\kappa_{\mathsf{min}}^{2}},\quad\kappa_{\mathsf{min}}=\tfrac{h/2-\left[\left(\max_{d\in\left\\{2,\dots,D\right\\}}T_{d}\Omega_{d}\right)e\right]^{N}\|f\|_{\infty}}{\sigma\sqrt{2^{N+1}}}\cdot\sqrt{\tfrac{B^{D-2}}{T_{\mathsf{max}}^{D-2}}}.$
(50)

We note that we do not use all filters
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$. Given that we use a set of
samples that is contiguous along all dimensions, any band $\bar{\mathbf{b}}$
containing samples has at least one neighboring band
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$ that
contains samples. Then, each constant
$\widetilde{M}_{\bar{\mathbf{b}}}=M_{\bar{\mathbf{b}}}-M_{\bar{\mathbf{0}}}$
can be computed using a single evaluation of
$\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$, which, in total, is evaluated
$\prod_{d=2}^{D}B_{\mathsf{max}}^{d}-1$ times, and the theorem follows. ∎

###  5.2 Comparison to ideal modulo recovery with Gaussian noise measurements

The USF was not analysed in the presence of Gaussian noise, but rather on
bounded noise [12, 28]. However, USF can still be applied for recovery in the
context of Gaussian noise, and the reconstruction would still be accurate in
the instances when the noise sample with maximum amplitude satisfies the USF
conditions. The modulo operator with Gaussian noise was considered before, but
the main objective was denoising, rather than input reconstruction [32, 31].
In order to assess the advantage of the new
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$ operator, we provide some insight
on reconstruction via USF for multi-dimensional inputs. Specifically, we note
that, for an input
$f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ and a
lattice $\boldsymbol{\Lambda}$, the ideal modulo output is decomposed as (see
also Section 2.1)

$y\left[\mathbf{k}\right]=\gamma\left[\mathbf{k}\right]-\varepsilon_{\gamma}\left[\mathbf{k}\right]+\eta\left[\mathbf{k}\right].$

The recovery method from [28] involves a _line-by-line_ approach, meaning that
the recovery is performed along dimension $k_{1}$,
$\forall\bar{\mathbf{k}}=\left[k_{2},\dots,k_{D}\right]$. By defining
$y_{\bar{\mathbf{k}}}\left[k_{1}\right]=y\left[\mathbf{k}\right]$,
$\gamma_{\bar{\mathbf{k}}}\left[k_{1}\right]=\gamma\left[\mathbf{k}\right]$,
$\gamma_{\bar{\mathbf{k}}}\left[k_{1}\right]=\gamma\left[\mathbf{k}\right]$,
$\varepsilon_{\gamma,\bar{\mathbf{k}}}\left[{k_{1}}\right]=\varepsilon_{\gamma}\left[\mathbf{k}\right]$,
and $\eta_{\bar{\mathbf{k}}}\left[{k_{1}}\right]=\eta\left[\mathbf{k}\right]$,
the recovery is performed by computing

$\displaystyle\begin{split}\left\langle
y_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle&=\left\langle\gamma_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle-\left\langle\varepsilon_{\gamma,\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\\\
&+\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle.\end{split}$
(51)

We remark that the processing in (51) is equivalent to applying filter
$\psi_{\bar{\mathbf{b}},m}$ (36) in the case of the multidimensional modulo-
hysteresis operator when there is only one sample per band in all dimensions,
i.e., $N_{2}^{B}=N_{3}^{B}=\dots=N_{d}^{B}=N^{B}=1$. Even though the noise
here is not bounded, we can derive the condition when the USF would work for a
specific noise instance, which is [12]

$\left\lvert\left\langle\gamma_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle+\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert\leqslant\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}+\left\lvert\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert<\lambda.$

for all $m\in\mathbb{Z}$. Thus, recovery is only guaranteed if the noise
instance is bounded by

$\left\lvert\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert<\lambda-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}.$
(52)

In a similar fashion to the derivation of (98) it can be shown that
$\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\sim\mathcal{N}\left(0,\sigma^{2}\cdot
2^{N}\right)$. We remark that the standard deviation of
$\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle$
is always at least $\sigma\sqrt{2}$. Therefore, depending on the values of
$T_{1},\Omega_{1}$ and $\lambda$, the probability that (52) holds may be very
low. Conversely, in the recovery with the proposed operator
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$,
$\left\langle\eta_{\bar{\mathbf{k}}},\psi_{\bar{\mathbf{b}},m}\right\rangle$
satisfies
$\left\langle\eta_{\bar{\mathbf{k}}},\psi_{\bar{\mathbf{b}},m}\right\rangle\sim\mathcal{N}\left(0,\sigma^{2}\cdot\frac{2^{N}}{N^{B}}\right)$
(98), where the standard deviation can be made arbitrarily small by increasing
the number of samples $N^{B}$ within each band $\bar{\mathbf{b}}$. In Section
6 this fact will be exploited to achieve significantly higher recovery
performance for the proposed operator
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$ compared to ideal modulo
$\mathscr{M}_{\lambda}$.

##  6 Numerical study

Let $\mathbf{V}=\left[\mathbf{v}_{1},\mathbf{v}_{2}\right]$ be a randomly
generated matrix such that $\|\mathbf{v}_{1}\|_{2}=\|\mathbf{v}_{1}\|_{2}=1$.
The input $f:\mathbb{R}^{2}\rightarrow\mathbb{R}$ was restricted to two
variables $x_{1},x_{2}$ for visualisation purposes, and was generated as
$f\left(\mathbf{x}\right)=f_{0}\left(\mathbf{Vx}\right)$, where

$f_{0}\left(x_{1},x_{2}\right)=\sum_{\mathbf{k}\in\mathbb{Z}^{2}}c_{\mathbf{k}}\frac{\sin\left(\Omega_{1}\left(x_{1}-k_{1}\frac{\pi}{\Omega_{1}}\right)\right)}{{\Omega_{1}\left(x_{1}-k_{1}\frac{\pi}{\Omega_{1}}\right)}}\cdot\frac{\sin\left(\Omega_{2}\left(x_{2}-k_{2}\frac{\pi}{\Omega_{2}}\right)\right)}{{\Omega_{2}\left(x_{2}-k_{2}\frac{\pi}{\Omega_{2}}\right)}}.$
(53)

Figure 3: Reconstruction accuracy comparison for (a) the ideal modulo
$\mathscr{M}_{\lambda}$ and (b) modulo-hysteresis
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{D}$ for a two-dimensional input.

We selected $\Omega_{1}=\Omega_{2}=1$, and computed $f$ for
$\mathbf{x}\in\left[-5,5\right]^{2}$. The coefficients $c_{\mathbf{k}}$ were
randomly generated for $\left\lvert k_{1}\right\rvert\leqslant 1,\left\lvert
k_{2}\right\rvert\leqslant 1$, drawn from the uniform distribution on
$\left[-1,1\right]$. The dynamic range of $f$ is $\left[-1,1\right]$. We
encoded $f$ using the ideal modulo $\mathscr{M}_{\lambda}$ with threshold
$\lambda=0.3$ and the proposed modulo-hysteresis
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{2}$ with $\lambda=0.3,h=0.19,B=0.32$.
The output samples, computed on lattice $\boldsymbol{\Lambda}$ with basis
vectors $\mathbf{v}_{1},\mathbf{v}_{2}$ and sampling periods $T_{1},T_{2}$. We
kept $T_{1}=0.02\ \mathrm{s}$ constant because its choice affects in a similar
fashion recovery from $\mathscr{M}_{\lambda}$ and
$\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{2}$ (see [14]). We varied $T_{2}$ in
the range $\left[0.005\ s,0.08\ s\right]$. The output samples are

$\displaystyle
y_{1}\left[\mathbf{k}\right]=\mathscr{M}_{\lambda}\left(f\left(\mathbf{VTk}\right)\right)+\eta\left[\mathbf{k}\right],$
(54) $\displaystyle
y_{2}\left[\mathbf{k}\right]=\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{2}f\left(\mathbf{VTk}\right)+\eta\left[\mathbf{k}\right],$
(55)

where $\eta\left[\mathbf{k}\right]\sim\mathcal{N}\left(0,\sigma^{2}\right)$.
We varied $\sigma$ in the range $\left[0.04,0.08\right]$, and recovered the
input
$\widetilde{\gamma}\left[\mathbf{k}\right]=f\left(\mathbf{VTk}\right)+\eta\left[\mathbf{k}\right]$
up to a constant multiple of $2\lambda$ for $\mathscr{M}_{\lambda}$ and $2h$
for $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{2}$. We generated $100$ random
inputs $f$ and $100$ noise sequences $\eta\left[\mathbf{k}\right]$ and counted
the number of inputs correctly reconstructed using each method. In our
context, correctly reconstructed means that the recovery conditions hold true.
The results are depicted in Fig. 3. We note that, while the accuracy increases
significantly for $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{2}$ for small
$T_{2}$, as proven by Theorem 3, for $\mathscr{M}_{\lambda}$ the reverse
happens. This is because $\mathscr{M}_{\lambda}$ processes the input in a
line-by-line fashion, and does not exploit in any way the higher resolution
along dimension $x_{2}$. In fact, here a higher resolution simply adds more
noise samples from sequence $\eta\left[\mathbf{k}\right]$, which are not
filtered and thus increase the probability that (52) does not hold. We also
note that, for larger sampling periods $T_{2}\sim 0.08\ \mathrm{s}$,
$\mathscr{M}_{\lambda}$ performs slightly better for low noise, i.e.,
$\sigma<0.05$. This is explained by the fact that the modulo-hysteresis
requirement
$\left\lvert\left\langle\eta_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert<h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}$
is more strict than (52) given that $h/2<\lambda$. This is a small trade-off
that enables $\mathscr{M}_{\boldsymbol{\mathsf{H}}}^{2}$ to handle arbitrarily
large values of $\sigma$ for small enough sampling periods $T_{2}$.

##  7 Proofs

###  7.1 Multi-Dimensional Modulo Properties

###### Proof for Proposition 1 (Bandlimited slices).

We consider the slice $f_{\mathbf{V}}=f\left(\mathbf{Vx}\right)$ along
dimension $x_{d}$, with $d\in\left\\{1,\dots,D\right\\}$ fixed. To this end,
we apply the $(D-1)$ – dimensional inverse Fourier transform to
$f_{\mathbf{V}}$ corresponding to all variables apart from $x_{d}$, such that

$\displaystyle\frac{1}{\left(2\pi\right)^{D-1}}$
$\displaystyle\int_{\mathbb{R}^{D-1}}\mathcal{F}f_{\mathbf{V}}\left(\mathbf{\omega}\right)e^{\jmath\left\langle\mathbf{\omega},\mathbf{x}\right\rangle}d\omega_{1}\dots
d\omega_{d-1}d\omega_{d+1}\dots d\omega_{D}$ (56) $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}f_{\mathbf{V}}\left(\mathbf{x}\right)e^{-\jmath\omega_{d}x_{d}}dx_{d}$
(57) $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}f\left(\sum\nolimits_{d=1}^{D}x_{d}\mathbf{v}_{d}\right)e^{-\jmath\omega_{d}x_{d}}dx_{d}=\mathcal{F}g_{d}\left(\omega_{d}\right).$
(58)

Therefore, the spectrum of $g_{d}$ depends on the spectrum of
$f_{\mathbf{V}}\left(\mathbf{x}\right)$, which will be evaluated in the
following. To this end, via the change of variable
$\mathbf{x}=\mathbf{Vx}^{*}$, we get that

$\displaystyle\begin{split}\mathcal{F}f\left(\mathbf{\omega}\right)=\int_{\mathbb{R}^{D}}f\left(\mathbf{x}\right)e^{-\jmath\left\langle\mathbf{\omega},\mathbf{x}\right\rangle}d\mathbf{x}&=\int_{\mathbb{R}^{D}}f\left(\mathbf{Vx}^{*}\right)e^{-\jmath\left\langle\mathbf{\omega},\mathbf{Vx}^{*}\right\rangle}\left\lvert\det\left(\mathbf{V}\right)\right\rvert
d\mathbf{x}^{*}\\\
&=\int_{\mathbb{R}^{D}}f\left(\mathbf{Vx}^{*}\right)e^{-\jmath\left\langle\mathbf{\mathbf{V}^{\top}\omega},\mathbf{x}^{*}\right\rangle}\left\lvert\det\left(\mathbf{V}\right)\right\rvert
d\mathbf{x}^{*}\\\
&=\left\lvert\det\left(\mathbf{V}\right)\right\rvert\cdot\mathcal{F}f_{\mathbf{V}}\left(\mathbf{V}^{\top}\mathbf{\omega}\right).\end{split}$
(59)

Using $f\in\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ we
get that
$\mathrm{supp}\left(\mathcal{F}f\right)\subseteq\mathbf{V}^{-\top}\cdot\prod_{d=1}^{D}\left(-\Omega_{d},\Omega_{d}\right)$.
Then, via (59), we have that
$\mathrm{supp}\left(\mathcal{F}f_{\mathbf{V}}\right)\subseteq\prod_{d=1}^{D}\left(-\Omega_{d},\Omega_{d}\right)$,
and therefore
$\mathcal{F}f_{\mathbf{V}}\left(\mathbf{\omega}\right)=0,\forall\omega_{d}\in\mathbb{R},\left\lvert\omega_{d}\right\rvert>\Omega_{d}$
in (56). It follows that
$\mathcal{F}g_{d}\left(\omega_{d}\right)=0,\forall\omega_{d}\in\mathbb{R},\left\lvert\omega_{d}\right\rvert>\Omega_{d}$,
$d\neq 1$. Furthermore, choosing $d=1$ gives us
$\mathcal{F}f_{\bar{\mathbf{x}}}\left(\omega_{1}\right)=0,\forall\omega_{1}\in\mathbb{R},\left\lvert\omega_{1}\right\rvert>\Omega_{1}$,
which finalizes the proof. ∎

###### Proof for Proposition 27 (Folding time separation).

We first show an intermediate result in the following lemma.

######  Lemma 1.

For a fixed vector $\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$ indicating the modulo
band and $r\in\left\\{1,\dots,R\right\\}$, let
$g_{\bar{\mathbf{x}},r}:\mathbb{R}\rightarrow\mathbb{R}$ and
$g_{r}:\mathbb{R}\rightarrow\mathbb{R}$ be two functions defined as

$\displaystyle g_{\bar{\mathbf{x}},r}\left(x_{1}\right)$
$\displaystyle\triangleq
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\varepsilon_{\bar{\mathbf{b}},r-1}\left(x_{1}\right),\quad\forall\mathbf{x}\in\mathbb{R}^{D},$
(60) $\displaystyle g_{r}\left(x_{1}\right)$
$\displaystyle\triangleq\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(x_{1}\right)\right\rvert,\quad\forall
x_{1}\in\mathbb{R}.$ (61)

Then $g_{\bar{\mathbf{x}},r}$ and $\left\lvert
g_{\bar{\mathbf{x}},r}\right\rvert$ are Lipschitz-continuous as functions of
$\mathbf{x}=\left[x_{1},\bar{\mathbf{x}}\right]\in\mathbb{Z}^{D}$ and $g_{r}$
is Lipschitz-continuous as function of $x_{1}\in\mathbb{R}$. Furthermore, the
Lipschitz constant in all cases is
$\|\boldsymbol{\Omega}\|_{2}\cdot\|f\|_{\infty}$.

###### Proof.

We begin by deriving a bound for the Lischitz constant of function
$f_{\mathbf{V}}=f\left(\mathbf{V}\cdot\right)$. Given that $f$ is
differentiable, according to the _mean value theorem_

$f_{\mathbf{V}}\left(\mathbf{x}\right)-f_{\mathbf{V}}\left(\boldsymbol{\raisebox{1.5pt}{$\chi$}}\right)=\left\langle\nabla
f_{\mathbf{V}}\left(\mathbf{c}\right),\mathbf{x}-\boldsymbol{\raisebox{1.5pt}{$\chi$}}\right\rangle,\quad\forall\mathbf{x},\boldsymbol{\raisebox{1.5pt}{$\chi$}}\in\mathbb{R}^{D},$
(62)

where
$\mathbf{c}=\alpha\mathbf{x}+\left(1-\alpha\right)\boldsymbol{\raisebox{1.5pt}{$\chi$}},\alpha\in\left[0,1\right]$
is an intermediate point on the segment joining $\mathbf{x}$ and $\chi$. Due
to the _Cauchy-Schwartz_ inequality

$\left\lvert
f_{\mathbf{V}}\left(\mathbf{x}\right)-f_{\mathbf{V}}\left(\boldsymbol{\raisebox{1.5pt}{$\chi$}}\right)\right\rvert\leqslant\|\nabla
f_{\mathbf{V}}\left(\mathbf{c}\right)\|_{2}\cdot\|\mathbf{x}-\boldsymbol{\raisebox{1.5pt}{$\chi$}}\|_{2}.$
(63)

Next, we use the Bernšteĭn bounds in (15) to derive a bound for $\|\nabla
f_{\mathbf{V}}\|$ as follows

$\displaystyle\|\nabla
f_{\mathbf{V}}\|_{2}^{2}=\sum_{d=1}^{D}\left\lvert\frac{\partial}{\partial
x_{d}}f_{\mathbf{V}}\left(\mathbf{c}\right)\right\rvert^{2}\leqslant\sum_{d=1}^{D}\Omega_{d}^{2}\|f\|^{2}_{\infty}$
(64)

which means $f_{\mathbf{V}}$ is a Lipschitz-continuous function satisfying

$\left\lvert
f_{\mathbf{V}}\left(\mathbf{x}\right)-f_{\mathbf{V}}\left(\boldsymbol{\raisebox{1.5pt}{$\chi$}}\right)\right\rvert\leqslant\|\mathbf{\Omega}\|_{2}\cdot\|f\|_{\infty}\cdot\|\mathbf{x}-\boldsymbol{\raisebox{1.5pt}{$\chi$}}\|_{2},\quad\forall\mathbf{x},\boldsymbol{\raisebox{1.5pt}{$\chi$}}\in\mathbb{R}^{D}.$
(65)

For all
$\bar{\mathbf{x}},\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}\in\mathbb{R}^{D-1}$
and $\forall x_{1},\raisebox{1.5pt}{$\chi$}_{1}>\tau_{\bar{\mathbf{b}},r-1}$
we have
$\varepsilon_{\bar{\mathbf{b}},r-1}\left(x_{1}\right)=\varepsilon_{\bar{\mathbf{b}},r-1}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)$,
which implies

$\displaystyle\begin{split}\left\lvert\left\lvert
g_{\bar{\mathbf{x}},r}\left(x_{1}\right)\right\rvert-\lvert{g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}},r}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)}\rvert\right\rvert&\leqslant\left\lvert
g_{\bar{\mathbf{x}},r}\left(x_{1}\right)-g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}},r}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert=\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert\\\
&=\left\lvert
f_{\mathbf{V}}\left(\mathbf{x}\right)-f_{\mathbf{V}}\left(\boldsymbol{\raisebox{1.5pt}{$\chi$}}\right)\right\rvert\leqslant\|\mathbf{\Omega}\|_{2}\cdot\|f\|_{\infty}\cdot\|\mathbf{x}-\boldsymbol{\raisebox{1.5pt}{$\chi$}}\|_{2}.\end{split}$
(66)

The final step is to show that $g_{r}\left(x_{1}\right)$ is also Lipschitz,
i.e., that the supremum does not change the Lipschitz constant. To this end,
we note the following two properties of the supremum. Specifically, for
$\forall\epsilon>0$ and $\forall
x_{1},\raisebox{1.5pt}{$\chi$}_{1}\in\mathbb{R}$,
$\exists\bar{\mathbf{x}}_{\epsilon},\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}\in\mathcal{R}_{\bar{\mathbf{b}}}$
such that

$\displaystyle\begin{split}\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(x_{1}\right)\right\rvert\leqslant
g_{r}\left(x_{1}\right)<\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(x_{1}\right)\right\rvert+\epsilon,\\\
\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert\leqslant
g_{r}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)<\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert+\epsilon.\end{split}$
(67)

We derive that
$g_{r}\left(x_{1}\right)-g_{r}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\leqslant\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(x_{1}\right)\right\rvert+\epsilon-\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert$.
Similarly, we get
$g_{r}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)-g_{r}\left(x_{1}\right)\leqslant\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert+\epsilon-\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(x_{1}\right)\right\rvert$.
Finally, we restrict $x_{1},\raisebox{1.5pt}{$\chi$}_{1}$ in (67) to satisfy
$x_{1},\raisebox{1.5pt}{$\chi$}_{1}>\tau_{\bar{\mathbf{b}},r-1}$ and select
$\bar{\mathbf{x}},\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}$ in (66) as
$\bar{\mathbf{x}}=\bar{\mathbf{x}}_{\epsilon},\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}=\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}$,
which yields

$\displaystyle\left\lvert
g_{r}\left(x_{1}\right)-g_{r}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert$
$\displaystyle\leqslant\max\left\\{\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(x_{1}\right)\right\rvert-\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert,\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert-\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(x_{1}\right)\right\rvert\right\\}+\epsilon$
$\displaystyle\leqslant\max\left\\{\Big{\lvert}{\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(x_{1}\right)\right\rvert-\left\lvert
g_{\bar{\mathbf{x}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert}\Big{\rvert},\left\lvert\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(\raisebox{1.5pt}{$\chi$}_{1}\right)\right\rvert-\left\lvert
g_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{\epsilon}}\left(x_{1}\right)\right\rvert\right\rvert\right\\}+\epsilon$
$\displaystyle\leqslant\left\lvert{\Omega}_{1}\right\rvert\cdot\|f\|_{\infty}\cdot\left\lvert
x_{1}-\raisebox{1.5pt}{$\chi$}_{1}\right\rvert+\epsilon,\quad\forall\epsilon>0.$

Taking $\epsilon\rightarrow 0$ above proves the required result. ∎

We begin by evaluating
$g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)$

$g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)=f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},0}\right)-hM_{\bar{\mathbf{b}}}=f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},0}\right)-h\left\lfloor\frac{\inf_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}}\left(\tau_{\bar{\mathbf{b}},0}\right)+\lambda}{h}\right\rfloor+h.$
(68)

Using that $x-1<\left\lfloor x\right\rfloor\leqslant x$ we derive

$\displaystyle\begin{split}\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)&\geqslant\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left[f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},0}\right)-h\tfrac{\inf_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}}\left(\tau_{\bar{\mathbf{b}},0}\right)+\lambda}{h}+h\right]\\\
&=-\lambda+h\\\
\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)&\leqslant\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left[f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},0}\right)-h\tfrac{\inf_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}}\left(\tau_{\bar{\mathbf{b}},0}\right)+\lambda}{h}+2h\right]\\\
&\leqslant-\lambda+2h+\mathcal{D}f\leqslant\lambda-h.\end{split}$ (69)

Then, given that

$g_{1}\left(\tau_{\bar{\mathbf{b}},0}\right)=\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert
g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)\right\rvert=\max\left\\{\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}{g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)},-\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}{g_{\bar{\mathbf{x}},1}\left(\tau_{\bar{\mathbf{b}},0}\right)}\right\\},$

we get $g_{1}\left(\tau_{\bar{\mathbf{b}},0}\right)\leqslant\lambda-h$.
Furthermore, $g_{1}$ is Lipschitz-continuous due to Lemma 1, and therefore
continuous. Thus, given that
$\tau_{\bar{\mathbf{b}},1}=\inf\left\\{x_{1}>\tau_{\bar{\mathbf{b}},0}|g_{1}\left(x_{1}\right)=\lambda\right\\}$
is true by definition, we get
$g_{1}\left(\tau_{\bar{\mathbf{b}},1}\right)=\lambda$.

We then investigate the variation of function $g_{1}$ between folding times
$\tau_{\bar{\mathbf{b}},0}$ and $\tau_{\bar{\mathbf{b}},1}$ as follows. Using
Lemma 1, we get

$g_{1}\left(\tau_{\bar{\mathbf{b}},0}\right)-g_{1}\left(\tau_{\bar{\mathbf{b}},1}\right)\leqslant\Omega_{1}\|f\|_{\infty}\cdot\left\lvert\tau_{\bar{\mathbf{b}},1}-\tau_{\bar{\mathbf{b}},0}\right\rvert,$

and then we use that $g_{1}(\tau_{\bar{\mathbf{b}},1})=\lambda$ and
$g_{1}(\tau_{\bar{\mathbf{b}},0})<\lambda-h$ to derive

$h\leqslant\left\lvert
g_{1}\left(\tau_{\bar{\mathbf{b}},0}\right)-g_{1}\left(\tau_{\bar{\mathbf{b}},1}\right)\right\rvert.$

Therefore the folding times satisfy the following bound

$\left\lvert\tau_{\bar{\mathbf{b}},1}-\tau_{\bar{\mathbf{b}},0}\right\rvert\geqslant\frac{h}{\Omega_{1}\|f\|_{\infty}}.$
(70)

We note the resemblance between (70) the one-dimensional case (8). Next, we
will show by induction that
$g_{r}\left(\tau_{\bar{\mathbf{b}},r-1}\right)\leqslant\lambda-h,\
g_{r}\left(\tau_{\bar{\mathbf{b}},r}\right)=\lambda$ and then
$\left\lvert\tau_{\bar{\mathbf{b}},r}-\tau_{\bar{\mathbf{b}},r-1}\right\rvert\geqslant\frac{h}{\Omega_{1}\|f\|_{\infty}}$
holds where sequence
$\left\\{\tau_{\bar{\mathbf{b}},1},\dots,\tau_{\bar{\mathbf{b}},r}\right\\}$
is computed according to (22). The _base case_ $r=1$ is shown in the
derivation to (70). For the _induction step_ we proceed with computing
$g_{r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)$.

$\displaystyle
g_{r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)\right\rvert=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r}\right)-\varepsilon_{\bar{\mathbf{b}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right\rvert$
(71)
$\displaystyle=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r}\right)-\varepsilon_{\bar{\mathbf{b}},r-1}\left(\tau_{\bar{\mathbf{b}},r}\right)+hs_{\bar{\mathbf{b}},r}\mathbbmtt{1}_{\left[\tau_{\bar{\mathbf{b}},r}\infty\right)}\left(\tau_{\bar{\mathbf{b}},r}\right)\right\rvert$
(72)
$\displaystyle=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-hs_{\bar{\mathbf{b}},r}\right\rvert,\
\mathrm{where}\
s_{\bar{\mathbf{b}},r}=\mathrm{sign}\left[g_{\bar{\mathbf{b}}B,r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right].$
(73)

We have that
$\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right\rvert=\lambda$,
which leads to two possible cases

  1. 1.

$\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right\rvert=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)}=\lambda$

Here we use that

$\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)}-\inf_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)}=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r}\right)}-\inf_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r}\right)}<\mathcal{D}f,$

which then implies
$\inf_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)}>\lambda-\mathcal{D}f>\lambda-h/2>0$.
Then

$\mathrm{sign}\left[g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right]=s_{\bar{\mathbf{b}},r}=\mathrm{sign}\left[g_{\bar{\mathbf{b}}B,r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right]=1.$

Then, continuing the derivation in (73), we get

$g_{r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-hs_{\bar{\mathbf{b}},r}\right\rvert=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h\right\rvert.$
(74)

Next, we evaluate the sign of
$g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h$ as follows

$\displaystyle g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h$
$\displaystyle\leqslant\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h=\lambda-h,$
(75) $\displaystyle
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h$
$\displaystyle\geqslant\inf_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h=\lambda-\mathcal{D}f-h$
(76) $\displaystyle\geqslant\lambda-\frac{3h}{2}>0.$ (77)

The final inequality follows from the assumption $h<2\lambda/3$. It follows
that (74)

$g_{r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)-h}\leqslant\lambda-h.$
(78)

  2. 2.

$\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right\rvert=-\inf_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}{g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)}=\lambda$.
As before, here we prove that
$\mathrm{sign}\left[g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)\right]=-1$,
$g_{r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)=\sup_{\bar{\boldsymbol{x}}\in\mathcal{R}_{\bar{\boldsymbol{b}}}}\left\lvert
g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)+h\right\rvert$,
and finally
$g_{\bar{\mathbf{x}},r}\left(\tau_{\bar{\mathbf{b}},r}\right)+h\in[-\lambda+h,0)$
which leads to (78).

As in the _base case_ $r=1$, given that
$\tau_{\bar{\mathbf{b}},r+1}=\inf\left\\{x_{1}>\tau_{\bar{\mathbf{b}},r}|g_{r+1}\left(x_{1}\right)=\lambda\right\\}$
is true by definition and using the continuity of $g_{r+1}$, we get
$g_{r+1}\left(\tau_{\bar{\mathbf{b}},r+1}\right)=\lambda$. Using (78) and the
same reasoning as in the _base case_ leading to (70) the proposition follows.
∎

###### Proof for Proposition 3 (Well-defined operator).

For this to be true we need to show the existence of
$R_{\bar{\mathbf{b}}}^{+}$ in (25). First, we derive some preliminary
properties of $f_{\bar{\mathbf{x}}}$. Due to Proposition 1 we have that
$f_{\bar{\mathbf{x}}}\in\mathsf{PW}_{\Omega_{1}}\left(\mathbb{R}\right)$ and
implicitly $f_{\bar{\mathbf{x}}}\in L^{2}\left(\mathbb{R}\right)$. Using the
properties of the
$\mathsf{PW}_{\boldsymbol{\Omega}}\left(\mathbb{R}^{D}\right)$ space and
$\|f\|_{\infty}<\infty$, it follows that

$\lim_{x_{1}\rightarrow\infty}f_{\bar{\mathbf{x}}}\left(x_{1}\right)=0,\quad\forall\bar{\mathbf{x}}\in\mathbb{R}^{D-1}.$
(79)

We require to show that the limit in (79) is uniform for all
$\bar{\mathbf{x}}$, which is done in the following lemma.

######  Lemma 2 (Uniform convergence).

The following holds true

$\forall\epsilon>0,\exists x_{1}^{*}>0,\forall
x_{1}>x_{1}^{*},\forall\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}},\quad\mathrm{s.t.}\quad\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right\rvert<\epsilon.$

###### Proof.

We prove by contradiction. Thus we assume

$\exists\epsilon>0,\forall x_{1}^{*}>0,\exists
x_{1}>x_{1}^{*},\exists\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}},\quad\mathrm{s.t.}\quad\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right\rvert\geqslant\epsilon.$

We select
$x_{1}^{*}=n,x_{1}=x_{1,n},\bar{\mathbf{x}}=\bar{\mathbf{x}}_{n}\in\mathcal{R}_{\bar{\mathbf{b}}}$
which leads to

$\exists\epsilon>0,\forall n\in\mathbb{N}^{*},\exists
x_{1,n}>n,\exists\bar{\mathbf{x}}_{n}\in\mathcal{R}_{\bar{\mathbf{b}}},\quad\mathrm{s.t.}\quad\left\lvert
f_{\bar{\mathbf{x}}_{n}}\left(x_{1,n}\right)\right\rvert\geqslant\epsilon.$
(80)

We use that, because $\mathcal{R}_{\bar{\mathbf{b}}}$ is bounded, then any
sequence $\bar{\mathbf{x}}_{n}\in\mathcal{R}_{\bar{\mathbf{b}}}$ has a
subsequence $\bar{\mathbf{x}}_{\kappa_{m}}\in\mathcal{R}_{\bar{\mathbf{b}}}$
that converges to a point in the closure of $\mathcal{R}_{\bar{\mathbf{b}}}$,
i.e.,

$\lim_{m\rightarrow\infty}\bar{\mathbf{x}}_{\kappa_{m}}=\bar{\mathbf{x}}_{\infty}\in\mathrm{cl}\left(\mathcal{R}_{\bar{\mathbf{b}}}\right).$

We then select $n=\kappa_{m}$ in (80) and get

$\exists\epsilon>0,\forall m\in\mathbb{N}^{*},\exists
x_{1,\kappa_{m}}>\kappa_{m},\exists\bar{\mathbf{x}}_{\kappa_{m}}\in\mathcal{R}_{\bar{\mathbf{b}}},\quad\mathrm{s.t.}\quad\left\lvert
f_{\bar{\mathbf{x}}_{\kappa_{m}}}\left(x_{1,{\kappa_{m}}}\right)\right\rvert\geqslant\epsilon.$

Furthermore, $f$ is a Lipschitz-continuous function with constant
$\|\boldsymbol{\Omega}\|_{2}\cdot\|f\|_{\infty}$ as shown in (65), which
implies that

$\displaystyle\left\lvert
f_{\bar{\mathbf{x}}_{\kappa_{m}}}\left(x_{1,\kappa_{m}}\right)\right\rvert-\left\lvert
f_{\bar{\mathbf{x}}_{\infty}}\left(x_{1,\kappa_{m}}\right)\right\rvert$
$\displaystyle\leqslant\left\lvert
f_{\bar{\mathbf{x}}_{\kappa_{m}}}\left(x_{1,\kappa_{m}}\right)-f_{\bar{\mathbf{x}}_{\infty}}\left(x_{1,\kappa_{m}}\right)\right\rvert$
$\displaystyle\leqslant\|\boldsymbol{\Omega}\|_{2}\cdot\|f\|_{\infty}\cdot\|\bar{\mathbf{x}}_{\kappa_{m}}-\bar{\mathbf{x}}_{\infty}\|_{2}.$

This allows defining the following lower bound on $\left\lvert
f_{\bar{\mathbf{x}}_{\infty}}\left(x_{1,\kappa_{m}}\right)\right\rvert$

$\left\lvert
f_{\bar{\mathbf{x}}_{\infty}}\left(x_{1,\kappa_{m}}\right)\right\rvert\geqslant\left\lvert
f_{\bar{\mathbf{x}}_{\kappa_{m}}}\left(x_{1,\kappa_{m}}\right)\right\rvert-\|\boldsymbol{\Omega}\|_{2}\cdot\|f\|_{\infty}\cdot\|\bar{\mathbf{x}}_{\kappa_{m}}-\bar{\mathbf{x}}_{\infty}\|_{2}.$

We know that $\left\lvert
f_{\bar{\mathbf{x}}_{\kappa_{m}}}\left(x_{1,\kappa_{m}}\right)\right\rvert\geqslant\epsilon$
and
$\lim_{m\rightarrow\infty}\|\bar{\mathbf{x}}_{\kappa_{m}}-\bar{\mathbf{x}}_{\infty}\|_{2}=0$.
Then it follows that $\left\lvert
f_{\bar{\mathbf{x}}_{\infty}}\left(x_{1,\kappa_{m}}\right)\right\rvert$ cannot
converge to $0$ for $m\rightarrow\infty$, which directly contradicts (79).
Then the starting assumption is wrong, and the Lemma follows. ∎

We approach this proof by contradiction. If we assume that
$R_{\bar{\mathbf{b}}}^{+}\geqslant 0$ does not exist it follows that
$\tau_{\bar{\mathbf{b}},r}$ in (22) is well-defined for $r\in\mathbb{Z}_{+}$.
Given that, by definition, $\tau_{\bar{\mathbf{b}},r}$ is increasing as a
function of $r$, then due to Proposition 27 it follows that
$\lim_{r\rightarrow\infty}\tau_{\bar{\mathbf{b}},r}=\infty$. We use Lemma 2
for $\epsilon=\frac{h}{4}$, which yields $\exists x_{1}^{*}>0$ such that
$\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right\rvert<\frac{h}{4},\forall\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}},\forall
x_{1}>x_{1}^{*}$. Given our assumption on $\tau_{\bar{\mathbf{b}},r}$ then
$\exists r^{*}\in\mathbb{Z}_{+}$ such that
$\tau_{\bar{\mathbf{b}},r^{*}}>x_{1}^{*}$. Thus, by definition,
$\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert
f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)-\varepsilon_{\bar{\mathbf{b}},r^{*}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)\right\rvert=\lambda$.
Using (73)–(78) where $r$ is replaced by $r^{*}$ and functions
$g_{r},g_{\bar{\mathbf{x}},r}$ are defined in (60),(61), one can show that
$g_{r^{*}+1}\left(\tau_{\bar{\mathbf{b}},r}\right)\leqslant\lambda-h$. The
next folding time $\tau_{\bar{\mathbf{b}},r^{*}+1}$ satisfies
$g_{r^{*}+1}\left(\tau_{\bar{\mathbf{b}},r^{*}+1}\right)=\lambda$. In the
following we will show this is not possible. Specifically, for
$x_{1}>\tau_{\bar{\mathbf{b}},r^{*}}$,

$\displaystyle g_{r^{*}+1}\left(x_{1}\right)$ $\displaystyle=\sup\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\varepsilon_{\bar{\mathbf{b}},r^{*}}\left(x_{1}\right)\right\rvert=\sup\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\varepsilon_{\bar{\mathbf{b}},r^{*}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)\right\rvert$
$\displaystyle=\sup\left\lvert
f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)-\varepsilon_{\bar{\mathbf{b}},r^{*}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)+f_{\bar{\mathbf{x}}}\left(x_{1}\right)-f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)\right\rvert$
$\displaystyle\leqslant\sup\left\lvert
f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)-\varepsilon_{\bar{\mathbf{b}},r^{*}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)\right\rvert+\sup\left\lvert
f_{\bar{\mathbf{x}}}\left(x_{1}\right)-f_{\bar{\mathbf{x}}}\left(\tau_{\bar{\mathbf{b}},r^{*}}\right)\right\rvert$
$\displaystyle\leqslant\lambda-h+h/2=\lambda-h/2<\lambda.$

We conclude that the definition of $\tau_{\bar{\mathbf{b}},r^{*}+1}$ via (22)
is therefore not possible, and thus our assumption that
$\tau_{\bar{\mathbf{b}},r}$ is well-defined for $r\in\mathbb{Z}$ is false, and
the proposition follows. ∎

###### Proof for Proposition 4 (Modulo output dynamic range).

We first assume that $x_{1}$ satisfies
$x_{1}\in\left[\tau_{\bar{\mathbf{b}},r},\tau_{\bar{\mathbf{b}},r+1}\right)$
and then extrapolate to the whole real axis. Given the definition of function
$\varepsilon_{\bar{\mathbf{b}},r}$ (24) it follows that
$\varepsilon\left(\mathbf{Vx}\right)=\varepsilon_{\bar{\mathbf{b}},r}\left(x_{1}\right)$
and thus

$z\left(\mathbf{Vx}\right)=f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\varepsilon_{\bar{\mathbf{b}},r}\left(x_{1}\right)=g_{\bar{\mathbf{x}},r+1}\left(x_{1}\right).$

It was shown before that
$g_{\bar{\mathbf{x}},r+1}\left(\tau_{\bar{\mathbf{b}},r}\right)\leqslant\lambda-h$
(78). Furthermore, due to the definition of $\tau_{\bar{\mathbf{b}},r+1}$ (22)
we get that $g_{r+1}\left(\tau_{\bar{\mathbf{b}},r+1}\right)=\lambda$ and
$\left\lvert g_{\bar{\mathbf{x}},r+1}\left(x_{1}\right)\right\rvert<\lambda$
when our assumption
$x_{1}\in\left[\tau_{\bar{\mathbf{b}},r},\tau_{\bar{\mathbf{b}},r+1}\right)$
holds. By repeating the process above for $r\in\mathbb{Z}_{+}$ we get
$\left\lvert z\left(\mathbf{Vx}\right)\right\rvert<\lambda,\forall
x_{1}\geqslant 0$. For $x_{1}<0$, the process above is reproduced for
$x_{1}\in\left(\tau_{\bar{\mathbf{b}},r-1},\tau_{\bar{\mathbf{b}},r}\right],\forall
r\in\mathbb{Z}_{-}$. ∎

###### Proof for Corollary 1 (Bound for intra-band variation).

Let
$\mathbf{x},\mathbf{\raisebox{1.5pt}{$\chi$}}\in\mathcal{R}_{\bar{\mathbf{b}}}$.
As shown before (65), $f_{\mathbf{V}}=f\left(\mathbf{V}\cdot\right)$ is a
Lipschitz-continuous function with constant
$\|\mathbf{\Omega}\|_{2}\cdot\|f\|_{\infty}$. Then it follows that

$\displaystyle\left\lvert
f\left(\mathbf{Vx}\right)-f\left(\mathbf{V\raisebox{1.5pt}{$\chi$}}\right)\right\rvert$
$\displaystyle\leqslant\|\mathbf{\Omega}\|_{2}\cdot\|f\|_{\infty}\cdot\|\mathbf{x}-\mathbf{\raisebox{1.5pt}{$\chi$}}\|_{2}$
(81) $\displaystyle\leqslant\|\mathbf{\Omega}\|_{2}\cdot\|f\|_{\infty}\cdot
B\sqrt{D}.$ (82)

We recall that $\mathcal{D}f$ satisfies (26)

$\mathcal{D}f\triangleq\sup_{\bar{\mathbf{b}}\in\mathbb{Z}^{D-1},x_{1}\in\mathbb{R}}\left[\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right].$
(83)

We fix a $\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$, and define
$\mathcal{D}_{\bar{\mathbf{b}}}f$ as

$\mathcal{D}_{\bar{\mathbf{b}}}f\triangleq\sup_{x_{1}\in\mathbb{R}}\left[\sup_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(x_{1}\right)-\inf_{\bar{\mathbf{x}}\in\mathcal{R}_{\bar{\mathbf{b}}}}f_{\bar{\mathbf{x}}}\left(x_{1}\right)\right].$
(84)

We then use the property that there always exists a sequence in a set
converging to the infimum or supremum. Therefore,
$\exists\bar{\mathbf{x}}_{n},\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{n},x_{1,n}$
such that

$\lim_{n\rightarrow\infty}f_{\bar{\mathbf{x}}_{n}}\left(x_{1,n}\right)-f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{n}}\left(x_{1,n}\right)=\mathcal{D}_{\bar{\mathbf{b}}}f.$
(85)

Given that $f_{\bar{\mathbf{x}}_{n}}\left(x_{1,n}\right)$ converges to the
supremum and
$f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{n}}\left(x_{1,n}\right)$ to
the infimum, it follows that $\exists N\in\mathbb{Z}_{+}$ such that
$f_{\bar{\mathbf{x}}_{n}}\left(x_{1,n}\right)-f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{n}}\left(x_{1,n}\right)\geqslant
0,\forall n\in\mathbb{Z},n\geqslant N$. Then

$\displaystyle
f_{\bar{\mathbf{x}}_{n}}\left(x_{1,n}\right)-f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{n}}\left(x_{1,n}\right)$
$\displaystyle=\left\lvert
f_{\bar{\mathbf{x}}_{n}}\left(x_{1,n}\right)-f_{\bar{\boldsymbol{\raisebox{1.5pt}{$\chi$}}}_{n}}\left(x_{1,n}\right)\right\rvert,\quad\forall
n\geqslant N$ (86) $\displaystyle\leqslant\|f\|_{\infty}\cdot
B\sqrt{D}\cdot\sqrt{\sum_{d=1}^{D}\Omega_{d}^{2}}.$ (87)

By taking $n\rightarrow\infty$ above we have that
$\mathcal{D}_{\bar{\mathbf{b}}}f\leqslant\|f\|_{\infty}\cdot
B\sqrt{D}\cdot\sqrt{\sum_{d=1}^{D}\Omega_{d}^{2}},\quad\forall\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$,
and computing
$\mathcal{D}f=\sup_{\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}}\mathcal{D}_{\bar{\mathbf{b}}}f$
leads to the desired bound. ∎

###  7.2 Multi-Dimensional Modulo Recovery

###### Proof for Proposition 6 (Detection of folding times).

As in the continuous-time scenario, we denote the one-dimensional slices of
the samples in each band $\bar{\mathbf{b}}$ by

$\displaystyle\gamma_{\bar{\mathbf{k}}}\left[k_{1}\right]$
$\displaystyle=\gamma\left[k_{1}T_{1}\mathbf{v}_{1}+\sum_{d=2}^{D}k_{d}T_{d}\mathbf{v}_{d}\right],$
(88) $\displaystyle\varepsilon_{\gamma,\bar{\mathbf{k}}}\left[k_{1}\right]$
$\displaystyle=\varepsilon_{\gamma}\left[k_{1}T_{1}\mathbf{v}_{1}+\sum_{d=2}^{D}k_{d}T_{d}\mathbf{v}_{d}\right],$
(89) $\displaystyle\eta_{\bar{\mathbf{k}}}\left[k_{1}\right]$
$\displaystyle=\eta\left[k_{1}T_{1}\mathbf{v}_{1}+\sum_{d=2}^{D}k_{d}T_{d}\mathbf{v}_{d}\right].$
(90)

where $\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}$ and
$\bar{\mathbf{k}}=\left[k_{2},\dots,k_{D}\right]\in\mathbb{Z}^{D-1}$. We note
that, due to Definition 2,
$\varepsilon_{\gamma,\bar{\mathbf{k}}}\left[k_{1}\right]$ does not change with
$\bar{\mathbf{k}}$ as long as
$\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}$. The
filtered samples $y_{\bar{\mathbf{b}},m}$ satisfy

$y_{\bar{\mathbf{b}},m}=\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle=\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle-\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle+\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle$

We first exploit that
$\gamma_{\bar{\mathbf{k}}}\left[k_{1}\right]=f_{\bar{\mathbf{T}}\bar{\mathbf{k}}}\left(k_{1}T_{1}\right)$
where $f_{\bar{\mathbf{x}}}$ is bandlimited to $\Omega_{1}\ \mathrm{rad/s}$,
which yields [12, 14]

$\displaystyle\begin{split}\left\lvert\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert&=\bigg{\lvert}\tfrac{1}{N^{B}}\sum_{\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\langle\gamma_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\bigg{\rvert}\\\
&\leqslant\tfrac{1}{N^{B}}\sum_{\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert\left\langle\gamma_{\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert\\\
&\leqslant\tfrac{1}{N^{B}}\sum_{\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}=\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}.\end{split}$
(91)

Similarly, for the residual samples $\varepsilon_{\gamma}$ the following holds

$\displaystyle\begin{split}\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert&=\bigg{\lvert}\tfrac{1}{N^{B}}\sum_{\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\langle\varepsilon_{\gamma,\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\bigg{\rvert}\\\
&\leqslant\tfrac{1}{N^{B}}\sum_{\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}}\left\lvert\left\langle\varepsilon_{\gamma,\bar{\mathbf{k}}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert=\left\lvert\left\langle\varepsilon_{\gamma,\bar{\mathbf{k}}^{*}},\Delta^{N}\left[\cdot-m\right]\right\rangle\right\rvert,\end{split}$
(92)

for all $\bar{\mathbf{k}}^{*}$ satisfying
$\bar{\mathbf{T}}\bar{\mathbf{k}}^{*}\in\mathcal{R}_{\bar{\mathbf{b}}}$, given
that $\varepsilon_{\gamma,\bar{\mathbf{k}}}$ does not change within the band
$\bar{\mathbf{b}}$ as a function of $\bar{\mathbf{k}}$ (24). We note that
$\mathrm{supp}\left(\Delta^{N}\left[\cdot-m\right]\right)=\left\\{m,\dots,m+N\right\\}$
and derive that

$\left\\{\begin{array}[]{cc}\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant
h&\mathrm{if}\ m\in\mathbb{S}_{N},\\\
\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert=0&\mathrm{otherwise},\end{array}\right.$
(93)

Therefore, if $\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}<h/2$ and
$\left\lvert\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert<h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}$,
then, $\mathrm{if}\ m\in\mathbb{S}_{N}$,

$\displaystyle\begin{split}\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert&=\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle-\left(\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle+\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle\right)\right\rvert\\\
&\geqslant\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert-\left\lvert\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle+\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\\\
&\geqslant h-h/2=h/2.\end{split}$ (94)

Furthermore, for $m\not\in\mathbb{S_{N}}$, we get $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\leqslant\left\lvert\left\langle\gamma,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert+\left\lvert\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert<h/2$.
Then, we identify if $m\in\mathbb{S}_{N}$ by thresholding sequence
$\left\langle y,\psi_{\bar{\mathbf{b}},m}\right\rangle$ via the following
inequalities.

$\left\\{\begin{array}[]{cc}\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant h/2&\mathrm{if}\
m\not\in\mathbb{S}_{N},\\\ \left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert<h/2&\mathrm{otherwise},\end{array}\right.$
(95)

Sequence $\eta\left[\bar{\mathbf{k}}\right]$ is drawn from the normal
distribution and is not bounded, therefore we can only guarantee that
$\left\lvert\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert<h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}$
holds with a given probability. However, we will show that modulo-hysteresis
allows increasing this probability exponentially.

We use two properties of the p.d.f. of Gaussian distributions. First, given
$M$ random variables $\eta_{i}\sim\mathcal{N}\left(0,\sigma_{i}^{2}\right),$
$i=1,\dots,M$ their summation satisfies
$\sum_{i=1}^{M}\eta_{i}\sim\mathcal{N}\left(0,\sigma^{2}\right)$, where
$\sigma=\sqrt{\sigma_{1}^{2}+\dots+\sigma_{M}^{2}}$. Second, for a random
variable $\eta_{1}\sim\mathcal{N}\left(0,\sigma^{2}\right)$, the
multiplication with a constant $\alpha$ yields
$\alpha\eta_{1}\sim\mathcal{N}\left(0,\left(\alpha\sigma\right)^{2}\right)$.

Then the noise term $\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle$
satisfies

$\displaystyle\begin{split}\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle&=\sum_{\bar{\mathbf{k}}\in\mathbb{Z}^{D-1}}\sum_{k_{1}\in\mathbb{Z}}\frac{\Delta^{N}\left[k_{1}-m\right]\cdot\eta\left[\mathbf{k}\right]}{N^{B}}\\\
&=\sum_{k_{1}\in\mathbb{Z}}\Delta^{N}\left[k_{1}-m\right]\cdot\bar{\eta}\left[k_{1}\right],\quad\bar{\eta}\left[k_{1}\right]=\frac{1}{N^{B}}\sum_{\bar{\mathbf{k}}\in\mathbb{Z}^{D-1}}\eta\left[\mathbf{k}\right].\end{split}$
(96)

Using the two p.d.f. properties above,
$\bar{\eta}\sim\mathcal{N}\left(0,\left(\sigma\sqrt{N^{B}}/N^{B}\right)^{2}\right)=\mathcal{N}\left(0,\frac{\sigma^{2}}{N^{B}}\right)$.
As expected, averaging gradually narrows down the p.d.f. of the distribution
around the origin. Furthermore, we can write

$\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle=\Delta^{N}_{-}\ast\bar{\eta}\left[m\right],\quad\Delta^{N}_{-}\left[k_{1}\right]=\Delta^{N}\left[-k_{1}\right],\forall
k_{1}\in\mathbb{Z}.$ (97)

Given that $\Delta_{-}^{n}=\Delta_{-}^{1}\ast\dots\ast\Delta_{-}^{n-1},\forall
n\in\mathbb{Z}$, we will compute recursively the p.d.f. of
$\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle$ as follows. First,
$\Delta_{-}^{1}\ast\bar{\eta}\left[m\right]=\bar{\eta}\left[m+1\right]-\bar{\eta}\left[m\right]$.
Given that
$-\bar{\eta}\left[m\right]\sim\mathcal{N}\left(0,\frac{\sigma^{2}}{N^{B}}\right)$,
we get
$\Delta_{-}^{1}\ast\bar{\eta}\left[m\right]\sim\mathcal{N}\left(0,\sigma^{2}{\frac{2}{N^{B}}}\right)$.
Recursively,

$\Delta_{-}^{N}\ast\bar{\eta}\left[m\right]\sim\mathcal{N}\left(0,\sigma^{2}{\frac{2^{N}}{N^{B}}}\right).$
(98)

In the equation above one can notice that the finite difference degree $N$
leads to an exponential increase in the standard deviation of the noise term.
A very similar result was reported for bounded noise in the one-dimensional
case [12, 14, 34]. However, in this multi-dimensional case we have the option
to decrease the noise by increasing the number of samples
$N^{B}=\prod_{d=2}^{D}\frac{B}{T_{d}}$. This can be done either by increasing
the band size $B$ within the allowable range ensuring
$\mathcal{D}f<\min\left\\{h/2,2\lambda-3h\right\\}$, but also by decreasing
the sampling periods $T_{d},d=2,\dots,D$. Both of these act only on dimensions
$x_{2},\dots,x_{D}$ and are fully independent of dimension $x_{1}$.

Therefore, the noise term
$\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle$ in (38) represents a
random variable that allows to correctly evaluate if $m\in\mathbb{S}_{N}$ via
(95) when
$\left\langle\eta,\psi_{\bar{\mathbf{b}},m}\right\rangle<h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}$.
The probability that this doesn’t hold is denoted by $p_{\mathsf{err}}$ which
is calculated using the p.d.f. of the normal distribution as

$p_{\mathsf{err}}=2\int_{\eta_{\mathsf{max}}}^{\infty}\frac{1}{\sigma_{0}\sqrt{2\pi}}\cdot
e^{-\frac{1}{2}\left(\frac{x}{\sigma_{0}}\right)^{2}}dx,$ (99)

where $\sigma_{0}=\sigma\cdot\sqrt{\frac{2^{N}}{N^{B}}}$ and
$\eta_{\mathsf{max}}=h/2-\left(T_{1}\Omega_{1}e\right)^{N}\|f\|_{\infty}$. The
integral above can be bounded in terms of the _complementary error function_
as follows, which finalizes the proof [38, 39]

$p_{\mathsf{err}}\leqslant
e^{-\left(\frac{\eta_{\mathsf{max}}}{\sigma_{0}\sqrt{2}}\right)^{2}}.$ (100)

∎

###### Proof for Proposition 7 (Detection of constants
$M_{\bar{\mathbf{b}}}$).

We define by $y_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ the samples $y$
filtered with $\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ where
$\bar{\mathbf{b}}^{*}\in\mathrm{Neighbors}\left(\bar{\mathbf{b}}\right)$, such
that
$\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}=\left[\bar{\mathbf{b}}\right]_{d^{*}}+1$
and

$y_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}=\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle=\left\langle\gamma,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle-\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle+\left\langle\eta,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle$

Along the same lines as (91), we derive

$\left\lvert\left\langle\gamma,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\leqslant\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}.$
(101)

Along dimension $x_{d^{*}}$ and for
$\bar{\mathbf{T}}\bar{\mathbf{k}}\in\mathcal{R}_{\bar{\mathbf{b}}}\cup\mathcal{R}_{\bar{\mathbf{b}}^{*}},k_{1}=0$,
the support of $\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ is

$\left\\{N_{d^{*}}^{B}\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}-1,\dots,N_{d^{*}}^{B}\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}-1+N\right\\}.$

Furthermore, the discontinuity between the bands would be located in between
the samples
$\left(N_{d^{*}}^{B}\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}-1\right)T_{d^{*}},N_{d^{*}}^{B}\left[\bar{\mathbf{b}}^{*}\right]_{d^{*}}T_{d^{*}}$.
Using the same reasoning as before (92-93)

$\left\\{\begin{array}[]{cc}\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\geqslant
h&\mathrm{if}\ M_{\bar{\mathbf{b}}}\neq M_{\bar{\mathbf{b}}^{*}},\\\
\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert=0&\mathrm{otherwise},\end{array}\right.$
(102)

Therefore, if $\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}<h/2$
and
$\left\lvert\left\langle\eta,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert<h/2-\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}$,
then, $\mathrm{if}\ M_{\bar{\mathbf{b}}}\neq M_{\bar{\mathbf{b}}^{*}}$, as
before, we get (94)

$\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\geqslant
h/2.$

Therefore, as before (95), we identify if $M_{\bar{\mathbf{b}}}\neq
M_{\bar{\mathbf{b}}^{*}}$ by thresholding sequence $\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle$ via the following
inequalities:

$\left\\{\begin{array}[]{cc}\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\geqslant
h/2&\mathrm{if}\ M_{\bar{\mathbf{b}}}\neq M_{\bar{\mathbf{b}}^{*}},\\\
\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert<h/2&\mathrm{otherwise}.\end{array}\right.$
(103)

Assuming that $M_{\bar{\mathbf{b}}}\neq M_{\bar{\mathbf{b}}^{*}}$, we compute
$\mathrm{sign}\left(\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)$ as
follows. We first show that

$\mathrm{sign}\left(\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)=-\mathrm{sign}\left(\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)$
(104)

If we assume by contradiction that $\mathrm{sign}\left(\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)=\mathrm{sign}\left(\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)$
we get $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle+\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert=\left\lvert\left\langle\gamma,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert<h/2$.
However, from (102)-(103) we have that $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\geqslant
h/2$ and
$\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert>0$.
Given our assumption that the two quantities have the same sign, we get a
contradiction and thus (104) is true.

Using the expression of $\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}$ in
(37), the definition of the residual in (20) and Proposition 5, it can be
shown directly that

$\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle=h\left(M_{\bar{\mathbf{b}}^{*}}-M_{\bar{\mathbf{b}}}\right)\left(-1\right)^{N+1}\Rightarrow\mathrm{sign}\left(\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)=\left(M_{\bar{\mathbf{b}}^{*}}-M_{\bar{\mathbf{b}}}\right)\left(-1\right)^{N+1},$

and thus $\mathrm{sign}\left(\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)=\left(M_{\bar{\mathbf{b}}^{*}}-M_{\bar{\mathbf{b}}}\right)\left(-1\right)^{N}$.
Finally, we get that

$\begin{cases}M_{\bar{\mathbf{b}}^{*}}=M_{\bar{\mathbf{b}}}+\mathrm{sign}\left(\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right)\left(-1\right)^{N}&\mathrm{if}\
\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert\geqslant
h/2,\\\
M_{\bar{\mathbf{b}}^{*}}=M_{\bar{\mathbf{b}}}&\mathrm{otherwise}.\end{cases}$
(105)

Furthermore, the filtered noise satisfies (96-98)

$\left\langle\eta,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\sim\mathcal{N}\left(0,\sigma^{2}\frac{2^{N}}{N^{B*}}\right).$
(106)

Then, the probability that
$\left\lvert\left\langle\eta,\psi_{\bar{\mathbf{b}},\bar{\mathbf{b}}^{*}}\right\rangle\right\rvert<h/2-\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}$
doesn’t hold, denoted by $p_{\mathsf{err}}$, satisfies

$p_{\mathsf{err}}\leqslant
e^{-\left(\frac{\eta_{\mathsf{max}}}{\sigma_{0}\sqrt{2}}\right)^{2}}.$ (107)

where $\sigma_{0}=\sigma\cdot\sqrt{\frac{2^{N}}{N^{B*}}}$ and
$\eta_{\mathsf{max}}=h/2-\left(T_{d^{*}}\Omega_{d^{*}}e\right)^{N}\|f\|_{\infty}$.
∎

###### Proof for Theorem 2 (Noiseless input recovery).

We note that, for $\sigma\rightarrow 0$, the results in Proposition 6 and 7
hold true with probability $1$. We then compute
$\widetilde{M}_{\bar{\mathbf{b}}}$ using Proposition 7 as follows. Given that
$\widetilde{M}_{\bar{\mathbf{0}}}=0$ by definition, one can compute
successively $\widetilde{M}_{\left[b_{2},0,\dots,0\right]}$ from
$\widetilde{M}_{\left[b_{2}-1,0,\dots,0\right]}$ for $\forall
b_{2}\in\mathbb{Z}_{+}$ and subsequently
$\widetilde{M}_{\left[b_{2},0,\dots,0\right]}$ from
$\widetilde{M}_{\left[b_{2}+1,0,\dots,0\right]}$ for $\forall
b_{2}\in\mathbb{Z}_{-}$. Repeating the process for
$b_{d},d\in\left\\{3,\dots,D\right\\}$ yields
$\widetilde{M}_{\bar{\mathbf{b}}},\forall\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$.

To compute the folding times via Proposition 6, we require that the sets
characterized by each folding time in $\mathbb{S}_{N}$ do not overlap.
Specifically we require that

$\left\\{\left\lceil\frac{\tau_{\bar{\mathbf{b}},r_{1}}}{T_{1}}\right\rceil-N,\dots,\left\lceil\frac{\tau_{\bar{\mathbf{b}},r_{1}}}{T_{1}}\right\rceil\right\\}\cap\left\\{\left\lceil\frac{\tau_{\bar{\mathbf{b}},r_{2}}}{T_{1}}\right\rceil-N,\dots,\left\lceil\frac{\tau_{\bar{\mathbf{b}},r_{2}}}{T_{1}}\right\rceil\right\\}=\emptyset,$

for $\forall r_{1},r_{2}\in\mathbb{Z},r_{1}\neq r_{2}$. A sufficient condition
for this is

$\left\lceil\frac{\tau_{\bar{\mathbf{b}},r_{2}}}{T_{1}}\right\rceil-N-\left\lceil\frac{\tau_{\bar{\mathbf{b}},r_{1}}}{T_{1}}\right\rceil>{\frac{\tau_{\bar{\mathbf{b}},r_{2}}}{T_{1}}}-{\frac{\tau_{,\bar{\mathbf{b}},r_{1}}}{T_{1}}}-\left(N+1\right)>0,$
(108)

which can be guaranteed via Proposition 27 if

$\left(N+1\right)T_{1}<\frac{h}{\Omega_{1}\|f\|_{\infty}}.$ (109)

Without reducing the generality we first assume that $x_{1}\geqslant 0$ and
thus $r,k_{1}\geqslant 0$. As before, the case $x_{1}\leqslant 0$ is treated
as a mirrored version of $x_{1}\geqslant 0$. An immediate consequence of (108)
is that for
$r_{1}=0,r_{2}=1\Rightarrow\left\lceil\frac{\tau_{\bar{\mathbf{b}},1}}{T_{1}}\right\rceil\geqslant
N+1$. Because there is no actual jump taking place at
$\tau_{\bar{\mathbf{b}},0}=0\Rightarrow\left\lvert\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},0}\right\rangle\right\rvert=0$
and $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},0}\right\rangle\right\rvert<h/2$ via Proposition 6.
The smallest $m$ for which filtered output satisfies $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant h/2$ is
$m=m_{\mathsf{min}}^{1}\triangleq\left\lceil\frac{\tau_{\bar{\mathbf{b}},1}}{T_{1}}\right\rceil-N$.
The last index $m$ corresponding to folding time $\tau_{\bar{\mathbf{b}},1}$
detected via $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant h/2$ is
$m=m_{\mathsf{max}}^{1}\triangleq\left\lceil\frac{\tau_{\bar{\mathbf{b}},1}}{T_{1}}\right\rceil$.
We can compute $m_{\mathsf{min}}^{1}$ and $m_{\mathsf{max}}^{1}$ as

$\displaystyle m_{\mathsf{min}}^{1}$
$\displaystyle=\min\left\\{m>0|\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant h/2\right\\},$
(110) $\displaystyle m_{\mathsf{max}}^{1}$
$\displaystyle=m_{\mathsf{min}}^{1}+N.$ (111)

Assuming (109) to be true, one can then compute recursively sequences
$m_{\mathsf{min}}^{r},m_{\mathsf{max}}^{r}$ corresponding to folding time
$\tau_{\bar{\mathbf{b}},r}$ as follows

$\displaystyle m_{\mathsf{min}}^{r}$
$\displaystyle=\min\left\\{m>m_{\mathsf{min}}^{r-1}+N|\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m}\right\rangle\right\rvert\geqslant h/2\right\\},$
(112) $\displaystyle m_{\mathsf{max}}^{r}$
$\displaystyle=m_{\mathsf{min}}^{r}+N.$ (113)

The folding time is estimated as
$\widetilde{\tau}_{\bar{\mathbf{b}},r}=\left[m_{\mathsf{min}}^{r}+N\right]\cdot
T_{1}$. As in the case $r=1$ we can show that
$m_{\mathsf{min}}^{r}=\left\lceil\frac{\tau_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil-N$.
Therefore
$\widetilde{\tau}_{\bar{\mathbf{b}},r}=\left\lceil\frac{\tau_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil\cdot
T_{1}$. Even though the folding time is not perfectly computed, this has no
effect on the input recovery because
$\left\lceil\frac{\widetilde{\tau}_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil=\left\lceil\frac{{\tau}_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil$
and we only evaluate the residual at the sampling locations $\mathbf{kT}$.
This means that replacing ${\tau}_{\bar{\mathbf{b}},r}$ by
$\left\lceil\frac{{\tau}_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil\cdot T_{1}$
in the expression of $\varepsilon_{f}\left(\mathbf{VTk}\right)$ yields the
same values (see Definition 2)

$\displaystyle\begin{split}\varepsilon_{\gamma}\left[\mathbf{k}\right]&=h\left[M_{\bar{\mathbf{b}}}+\sum_{i=0}^{r}{s}_{\bar{\mathbf{b}},r}\mathbbmtt{1}_{\left[{\tau}_{\bar{\mathbf{b}},r},\infty\right)}\left(k_{1}T_{1}\right)\right]\\\
&=h\left[M_{\bar{\mathbf{b}}}+\sum_{i=0}^{r}{s}_{\bar{\mathbf{b}},r}\mathbbmtt{1}_{\left[\left\lceil\frac{{\tau}_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil\cdot
T_{1},\infty\right)}\left(k_{1}T_{1}\right)\right].\end{split}$ (114)

We note that, as explained before, we do not recover $M_{\bar{\mathbf{b}}}$
but
$\widetilde{M}_{\bar{\mathbf{b}}}=M_{\bar{\mathbf{b}}}-M_{\bar{\mathbf{0}}}$.
This will be accounted for at the final input reconstruction stage.

Furthermore, we estimate the sign as
$\widetilde{s}_{\bar{\mathbf{b}},r}=-\mathrm{sign}\left\langle
y,\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle$. We will show
that $\widetilde{s}_{\bar{\mathbf{b}},r}={s}_{\bar{\mathbf{b}},r}$ as follows.
Given that
$\left\lvert\left\langle\gamma,\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle\right\rvert<h/2$
and $\left\lvert\left\langle
y,\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle\right\rvert\geqslant
h/2$, then, via (38), it follows that

$\mathrm{sign}\left\langle
y,\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle=-\mathrm{sign}\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle.$

We use the fact that $\varepsilon_{\gamma}\left[\mathbf{k}\right]$ does not
change for
$\bar{\mathbf{k}}\bar{\mathbf{T}}\in\mathcal{R}_{\bar{\mathbf{b}}}$. For
$N\geqslant 1$, using the expression of $\Delta^{N}$ and (114),

$\displaystyle\begin{split}\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle&=\left\langle\varepsilon_{\gamma},\Delta^{N}\left[\cdot-
m_{\mathsf{min}}^{r}\right]\right\rangle\\\
&=h{s}_{\bar{\mathbf{b}},r}\sum_{k_{1}\in\mathbb{Z}}\mathbbmtt{1}_{\left[\left\lceil\frac{{\tau}_{\bar{\mathbf{b}},r}}{T_{1}}\right\rceil\cdot
T_{1},\infty\right)}\left(k_{1}T_{1}\right)\cdot\Delta^{N}\left[k_{1}-m_{\mathsf{min}}^{r}\right]\\\
&=h{s}_{\bar{\mathbf{b}},r}\sum_{k_{1}\in\mathbb{Z}}\mathbbmtt{1}_{\left[\left(m_{\mathsf{min}}^{r}+N\right)\cdot
T_{1},\infty\right)}\left(k_{1}T_{1}\right)\cdot\Delta^{N}\left[k_{1}-m_{\mathsf{min}}^{r}\right]\end{split}$
(115)

By applying the change of variable $k_{1}^{*}=k_{1}-m_{\mathsf{min}}^{r}-N$

$\displaystyle\left\langle\varepsilon_{\gamma},\psi_{\bar{\mathbf{b}},m_{\mathsf{min}}^{r}}\right\rangle$
$\displaystyle=h{s}_{\bar{\mathbf{b}},r}\sum_{k_{1}^{*}\in\mathbb{Z}}\mathbbmtt{1}_{\left[0,\infty\right)}\left(k_{1}^{*}T_{1}\right)\cdot\Delta^{N}\left[k_{1}^{*}+N\right]$
(116)
$\displaystyle=h{s}_{\bar{\mathbf{b}},r}\sum_{k_{1}^{*}\in\mathbb{Z}_{+}}\Delta^{N}\left[k_{1}^{*}+N\right]=h{s}_{\bar{\mathbf{b}},r}.$
(117)

The last equality can be shown recursively via direct calculation for
$N\geqslant 1$, given that $k_{1}\geqslant 0$, which proves that
${s}_{\bar{\mathbf{b}},r}=\widetilde{s}_{\bar{\mathbf{b}},r}$.

After the folding times and signs are computed as above for all
$\bar{\mathbf{b}}\in\mathbb{Z}^{D-1}$, the input samples are reconstructed as

$\widetilde{\gamma}\left[\mathbf{k}\right]=y\left[\mathbf{k}\right]+\widetilde{\varepsilon}_{\gamma}\left[\mathbf{k}\right],$
(118)

where
$\widetilde{\varepsilon}_{\gamma}\left[\mathbf{k}\right]=h\left[\widetilde{M}_{\bar{\mathbf{b}}}+\sum_{i=0}^{r}{s}_{\bar{\mathbf{b}},r}\mathbbmtt{1}_{\left[{\tau}_{\bar{\mathbf{b}},r},\infty\right)}\left(k_{1}T_{1}\right)\right]=\varepsilon_{\gamma}\left[\mathbf{k}\right]-h{M}_{\bar{\mathbf{0}}}$,
which leads to
$\widetilde{\gamma}\left[\mathbf{k}\right]=\gamma\left[\mathbf{k}\right]-h{M}_{\bar{\mathbf{0}}}$.
∎

##  8 Conclusion

The Unlimited Sampling Framework (USF) provides sampling rate guarantees that
allow tackling high dynamic range signals in the one-dimensional case. For
multi-dimensional signals, USF is typically applied sequentially, thus not
exploiting the multi-dimensional structure of the input. In this paper, we

  * •

introduced the first multi-dimensional modulo operator and associated input
reconstruction method from lattice samples,

  * •

derived sampling rate conditions under which the reconstruction is perfect in
the noiseless scenario,

  * •

provided probability error bounds under Gaussian noise assumption,

  * •

showed numerically that, while USF does not allow noise amplitudes larger than
the modulo threshold, the proposed approach allows arbitrarily high noise for
sufficiently small sampling times.

This work can be extended in a number of ways

  1. 1.

It can be coupled with modulo denoising approaches such as [31] to yield
enhanced reconstruction algorithms.

  2. 2.

While it is assumed that the input is bandlimited, this work can be extended
for inputs generated with B-splines or sparse inputs.

  3. 3.

The model in our work can be extended to a wider range of models such as
modulo neuromorphic architectures, that could exploit multi-dimensional inputs
in a similar way to the biological systems [40] [24].

  4. 4.

Alternative sampling mechanisms that would benefit from a multi-dimensional
modulo operation include one-bit sampling [41] or average sampling [33].

  5. 5.

The current line of work can lead to a the implementation of a new multi-
dimensional hardware prototype.

## References

  * [1] R. M. Gray, D. L. Neuhoff, Quantization, IEEE Trans. Inf. Theory 44 (6) (1998) 2325–2383.  doi:10.1109/18.720541.
  * [2] W. Sun, X. Zhou, Reconstruction of band-limited functions from local averages, Constructive Approximation 18 (2) (2002) 205–222.  doi:10.1007/s00365-001-0011-y.