# On the Foundation of Sparse Sensing (Part I): Necessary and Sufficient
Sampling Theory and Robust Remaindering Problem

Hanshen Xiao, Yaowen Zhang, and Guoqiang Xiao Hanshen Xiao is with CSAIL and
the EECS Department, MIT, Cambridge, USA. E-mail<EMAIL_ADDRESS>Zhang
and Guoqiang Xiao are with the College of Computer and Information Science,
Southwest University, Chongqing, China. E-mail<EMAIL_ADDRESS>
###### Abstract

In the first part of the series papers, we set out to answer the following
question: given specific restrictions on a set of samplers, what kind of
signal can be uniquely represented by the corresponding samples attained, as
the foundation of sparse sensing. It is different from compressed sensing,
which exploits the sparse representation of a signal to reduce sample
complexity (compressed sampling or acquisition). We use sparse sensing to
denote a board concept of methods whose main focus is to improve the
efficiency and cost of sampling implementation itself. The “sparse” here is
referred to sampling at a low temporal or spatial rate (sparsity constrained
sampling or acquisition), which in practice models cheaper hardware such as
lower power, less memory and throughput.

We take frequency and direction of arrival (DoA) estimation as concrete
examples and give the necessary and sufficient requirements of the sampling
strategy. Interestingly, we prove that these problems can be reduced to some
(multiple) remainder model. As a straightforward corollary, we supplement and
complete the theory of co-prime sampling, which receives considerable
attention over last decade.

On the other hand, we advance the understanding of the robust multiple
remainder problem, which models the case when sampling with noise. A sharpened
tradeoff between the parameter dynamic range and the error bound is derived.
We prove that, for $N$-frequency estimation in either complex or real
waveforms, once the least common multiple (lcm) of the sampling rates selected
is sufficiently large, one may approach an error tolerance bound independent
of $N$.

##  I introduction

In 1928, Nyquist pointed out that for an arbitrary band-limited signal, once a
system uniformly samples it at a rate that exceeds the signal’s highest
frequency by at least a factor of two, termed as the Nyquist rate, the
discrete sample sequence acquired can perfectly recover the original signal
[1]. In 1948, Shannon rigorously proved Nyquist’s claim and pointed out that
Nyquist rate is indeed also the necessary condition that a uniform sampling
can uniquely recover the original band-limited signal in general [2]. The
Nyquist-Shannon Sampling Theorem serves as the foundation of modern signal
processing, which essentially characterizes the relationship between discrete
samples and a continuous signal. As for the non-uniform sampling, Landau in
1967 gave a generic conclusion, known as the Landau necessary condition, that
the average sampling rate must be twice the occupied bandwidth [3]. However,
those strong necessary conditions are not the end of the study on sampling
theory: Nyquist rate can be unnecessary given more assumptions on signals and
undersampling becomes possible provided stronger prior knowledge.

One representative is compressed sensing [4, 5, 6, 7], rooted in an under-
determined linear system, $y=\Phi x$, with $x\in\mathbb{R}^{d}$ and
$\Phi\in\mathbb{R}^{d^{\prime}\times d}$. In general, finding the unique
solution $x$ is impossible when the number of measurements $d^{\prime}$ is
less than the dimensionality $d$. Nonetheless, with sparsity restriction, say
$\|x\|_{0}\ll d$, and proper randomness in the encoding matrix $\Phi$, a
sparse $x$ can be reconstructed with high probability by solving an $L_{1}$
convex programming $\min_{z\in\mathbb{R}^{d}}\|y-\Phi z\|_{1}$ [8]. The above
statements provide a generic framework for signal reconstruction if a signal
can be sparsely represented in some domain. Theoretically, the number of
measurements can be almost only dependent on the sparsity parameter
$\|x\|_{0}$ for a high-probability reconstruction.

However, the implementation of compressed sampling is not always trivial in
practice, which is especially true when handling the analog signal, i.e., how
to transform the signal to a domain, where it can be represented by a sparse
vector $x$ and acquire samples through random linear measurements. In general,
before analyzing, communication or even storage of an analog signal, one has
to first convert it into a bit sequence digitally. Different from simply a
data processing, the digital representation of a signal is restricted by the
power, memory and other hardware constraints of the sampler applied. Though
compressed sensing presents a theoretical framework which allows fewer
measurements to reconstruct a sparse vector, the cost of each measurement may
not be always low. 111We have to stress that the point we want to make here is
that not all (sparse) signal processing tasks perfectly fit the compressed
sensing framework due to the sampler restriction. However, in some scenarios,
for example in Magnetic Resonance Imaging [9, 10, 11] and some optical imaging
systems [12], the signal of interest has already been represented by a sparse
vector and sampling operator itself can be modeled by a linear sensing, and
thus compressed sensing techniques definitely strengthen the sampling
efficiency.

This raises a fundamental question that, given specific restrictions, what
kind of signals can be represented by the corresponding sample sequence. Low
energy-cost hardware with loose architecture requirement is always desired in
many application scenarios. Recently, the series works [13, 14, 15], by
Kipnis, Eldar and Goldsmith, give a complete characterization of the
fundamental trade-off among sampling rate, compression bitrate (quantification
rate) and reconstruction error under the analog-to-digital framework, which
generalizes the classic sampling theorem to more practical situations.

On the other hand, to enable undersampling beyond the sampling theorem
restriction, roughly speaking, there are two main relaxation directions. One
is to make proper assumptions on the signal form, say standard complex
waveforms, where reconstructing the signal is instance-based and equivalent to
determining the unknown parameters, such as frequency and phase. The other is
to use multiple undersampling samplers with a finer grained coordination. On
the algorithmic level, there are two successful examples emerged during last
two decades, the Chinese Remainder Theorem (CRT) based parameter estimation,
short for CRT method in the following, pioneering by Xia [16, 17], and the co-
prime sampling (array), pioneering by Pal and Vaidyanathan [18, 19].

CRT Based Reconstruction: CRT method is spectrum based. A complex waveform
$x(t)=e^{j2\pi ft}$, with an unknown frequency $f$, is sampled uniformly at a
rate $m$, $m<f$, and Fourier transform is applied to the discrete sequence. It
is well known that undersampling will result in a spectrum aliasing: the
location of the spectrum peak is indeed the $f$ modulo $m$, denoted by
$\langle f\rangle_{m}$ in the following. Said another way, when $m<f$, the
residue, $\langle f\rangle_{m}=f-m\lfloor f/m\rfloor$, is of ambiguity to true
$f$, unless the folding number $\lfloor f/m\rfloor$ can be determined. On the
other hand, CRT characterizes the relationship between a number and its
residues modulo multiple smaller numbers, termed as moduli. Mathematically,
CRT states that given $L$ moduli, the set of residues can uniquely determine a
positive number $f$ once the least common multiple (lcm) of the moduli is
larger than $f$. Therefore, with a careful selection of sampling rates, say
pairwise co-prime, a frequency $f$ can be reconstructed via exponentially
small sampling rates $O(f^{\frac{1}{L}})$. This idea is further generalized to
other parameter estimations, such as phase unwrapping [20] with applications
in synthetic aperture radar (SAR) [21, 22].

However, in practice, CRT method encounters two main challenges. One is
robustness, where the residue system is non-weighted and conventional CRT is
very sensitive to small error. The follow-up works [23, 24, 25] addressed the
problem by adding redundancy to the moduli (sampling rate) by setting all
moduli to share a sufficiently large common divisor. The other challenge is
the correspondence relation ambiguity between multiple numbers and their
residues. Though a single objective can be elegantly estimated in the above
model, things becomes trickier when handling multiple objectives
simultaneously. In general, for a signal $x(t)=\sum_{i=1}^{N}e^{j2\pi f_{i}t}$
of multiple frequencies $\\{f_{1},f_{2},...,f_{N}\\}$, the correspondence
between the spectrum peaks and the frequencies is indistinguishable without
any further assumptions. Therefore, this problem cannot be simplified to $N$
independent single-frequency estimations. One strategy is to sacrifice the
encoding dynamic range, i.e., the range of $f_{i}$ can be uniquely determined,
to produce sufficient redundancy to enable decoding. Such tradeoff between the
dynamic range and $N$ is studied in [26], [27, 28], whereas the tight bound is
only known when $N=2$ [29]. To tackle both robustness and correspondence
ambiguity simultaneously, the polynomial-time (statistical) robust multi-
objective estimation is only known recently in [30, 31, 32, 33].

Co-prime Sampling (Array): Different from CRT method, which handles the
spectrum directly, co-prime sampling estimates the auto-correlation on the
temporal domain. Technically, the foundation of co-prime sampling is Bezout
Theorem, another fundamental number theory theorem, which states that for
arbitrary two numbers $P$ and $Q$, there exist $\alpha$ and $\beta$ such that

$P\beta-Q\alpha=gcd(P,Q),$

where $gcd(P,Q)$ denotes the greatest common divisor of $P$ and $Q$.
Especially, if $P$ and $Q$ are co-prime, i.e., $gcd(P,Q)=1$, as a
straightforward corollary, there exists $\beta_{1}$ and $\alpha_{1}$ such that

$P\beta_{1}-Q\alpha_{1}=1,$

and multiplying $l$ on both sides, we have

$P\cdot(l\beta_{1})-Q\cdot(l\alpha_{1})=l,$

for arbitrary $l$.

With the above preparation, now we can give a high-level picture of co-prime
sampling. Assume that the time interval corresponding to Nyquist rate is $T$
and two samplers with undersampling rates $P$ and $Q$, respectively, produce
two discrete sample sequences $\\{x[PTn],n\in\mathcal{N}\\}$ and
$\\{x[QTn],n\in\mathcal{N}\\}$ of a continuous signal $x(t)$. From Bezout
theorem, once $P$ and $Q$ are co-prime, there always exist $\alpha_{l}$ and
$\beta_{l}$ such that

$P\cdot\beta_{l}-Q\cdot\alpha_{l}=l,$

and correspondingly the pair $(x[P\beta_{l}T],x[Q\alpha_{l}T])$ renders an
estimation of the auto-correlation $\mathbb{E}_{n}\big{\\{}x[nT]\cdot
x^{*}[(n-l)T]\big{\\}}$ at lag $l$. Since auto-correlation and signal spectrum
are a pair of Fourier transform, the frequencies contained in a signal can be
estimated from the signal spectrum.

Switching to the spatial case, where $T$ is replaced by the half of
wavelength, the spatial sampling becomes to determine the physical location of
sensors and such construction is termed as co-prime array [19].

Though both CRT method and co-prime sampling are rooted in number theory, or
more precisely, the Diophantine equation theory, they were separately studied
over last decade. As far as we know, our paper is the first work to connect
them together and take them as the foundation of the study on sparse sensing.
In short, our conclusion is somewhat surprising: The underlying model in CRT
method is closer to characterize the necessary and sufficient condition that a
signal can be uniquely represented in sparse sensing, while co-prime sampling
indeed gives a more efficient “high-probability” reconstruction. 222The high
probability here is either with respect to natural distribution of the
parameter to be estimated, or a random undersampling rate selection for any
given parameters to be estimated, which will be specified later.

We start with these two representatives of sparse sensing to give some
intuition that the two fundamental questions we set out to address.

  1. 1.

First, theoretically, how “sparse”, regarding both the temporal and spatial
sampling rate, can we expect to be? Or equivalently, given the sparsity
restriction on the sampling, what kind of signal can be uniquely represented
by the attained samples.

  2. 2.

Second, under sparsity restriction, algorithmically, how “efficient”,
regarding both the freedom degree of parameters to be estimated and the time
complexity with associated performance, can the reconstruction scheme be?
Also, can the efficiency be significantly improved by resorting to a small
failure probability, relaxed from a deterministic reconstruction?

This paper mainly focuses on the first “encoding limitation” question to give
a systematical study on the necessary and sufficient conditions of sparse
sampling theory. Concrete sampling algorithms (“decoding”) will be presented
in the second paper, where interestingly we show given sufficiently many
samplers/senors, it could be sparse enough, i.e., arbitrarily low sampling
frequency in temporal domain and arbitrarily large spacing between sensors,
with a linear sampling time.

Organization and Contributions: We summarize the organization of rest contents
as follows. In Section II, we take frequency and DoA estimation as two
concrete examples and give the necessary and sufficient conditions of the
signal that can be represented given sampling sparsity restrictions. In
Section III, we give a complete analysis of co-prime sampling in Theorem 3,
which indicates that co-prime sampling is indeed a high-probability
reconstruction. In Section IV, we consider the more general case where
sampling with noise and we study this problem from a view of error-correcting
codes. We provide a characterization of encoding dynamic range and error
tolerance of complex waveforms in Theorem 4. In Section V, we provide further
generalization to the scenario of real waveforms. We conclude in Section VI.

##  II Necessary and Sufficient Conditions

In general, sampling can be viewed as a encoding procedure, where the
objective continuous signal is transformed to a discrete sample sequence.
Therefore, the necessary and sufficient conditions of sampling are equivalent
to determining if there is a bijection between the objective signal and the
encoded sample sequence. With the specification on the assumptions of signals
and restrictions on the samplers, such conditions are deterministic, at least
theoretically in the noiseless scenario. In the following, we take frequency
and DoA estimation as two concrete examples in the temporal and spatial
domain, respectively.

###  II-A Frequency and Phase Estimation

Without loss of generality, we consider a complex waveform

$x(t)=\sum_{i=1}^{N}A_{i}e^{j(2\pi f_{i}t+\theta_{i})},$ (1)

where $A_{i}$ is the amplitude, $f_{i}$ is the frequency and $\theta_{i}$ is
the phase of the $i$-th source. The continuous signal $x(t)$ is sampled by a
group of samples with sampling rates $\\{m_{1},m_{2},...,m_{L}\\}$,
respectively. Accordingly, $L$ sequences are obtained in a form,

$x_{l}[n]=\sum_{i=1}^{N}A_{i}e^{j\big{(}2\pi(f_{i}/m_{l})n+\theta^{l}_{i}\big{)}},$
(2)

where $\theta^{l}_{i}$ is the phase of the the $i$-th source at the $l$-th
sampler. Since Fourier Transform (FT) is an invertible transformation,
applying FT to (2), we have,

$\mathcal{F}_{l}(k)=FT(x_{l}[n])=\sum_{i=1}^{n}A_{i}\delta(k-\langle
f_{i}\rangle_{m_{l}})e^{j\theta^{l}_{i}},$ (3)

where $\delta(k)$ denotes the indicator, which equals $0$ iff $k=0$. Without
any prior assumptions on the rest parameters $\theta^{l}_{i}$ and $A_{i}$, the
sequence $x_{l}[n]$ at the $l$-th sampler is uniquely characterized by an
unordered set $\mathcal{R}_{l}=\\{r_{il}=\langle
f_{i}\rangle_{m_{l}},i=1,2,...,N\\}$. The necessary and sufficient conditions
that the unknown frequency set $\bm{f}=\\{f_{1},f_{2},...,f_{n}\\}$ can be
uniquely reconstructed with the $L$ sequences together are that, $N$
frequencies $\bm{f}\in[0,D)^{N}$, where $[0,D)$ represents the dynamic range,
are bijective to the encoded $L$ unordered residue sets $\mathcal{R}_{[1:L]}$.
It is worth mentioning that the above statement relies on the fact that no
additional prior knowledge is assumed. If $A_{i}$ or $\theta^{l}_{i}$ can be
distinguished, $\mathcal{R}_{l}$ becomes an “ordered” set, where the
correspondence between each $r_{il}$ and $f_{i}$ is determined. Given the
correspondence, the necessary and sufficient conditions become
$D<lcm(m_{1},m_{2},...,m_{L})$, characterized by regular CRT. We summarize the
conclusion as the following theorem.

######  Theorem 1.

Without additional prior assumptions, the necessary and sufficient conditions
that $N$ frequencies $\bm{f}=\\{f_{1},f_{2},...,f_{N}\\}$ within $[0,D)^{N}$
defined in (1) can be reconstructed given a set of samplers with rates
$\mathcal{M}=\\{m_{1},m_{2},...,m_{L}\\}$ are that there exists bijective
correspondence between $\bm{f}$ and $\mathcal{R}_{[1:L]}$, where
$\mathcal{R}_{l}=\\{r_{il}=\langle f_{i}\rangle_{m_{l}},i=1,2,...,N\\}$ is an
unordered residue set. Especially, if $A_{i}$ or $\theta^{l}_{i}$ are assumed
to be distinct and distinguishable, the necessary and sufficient conditions
become $D<lcm(m_{1},m_{2},...,m_{L})$.

Analogous to the frequency estimation, it is not hard to see that the
necessary and sufficient conditions of the phase-difference based distance
measurement share a similar form, where $\bm{f}$ will correspond to the
distance(s) between the target and the sensors to be estimated. We will
specify the dynamic range $D$ in Section IV and generalize the above analysis
to real waveforms in Section V.

###  II-B DoA Estimation

With a similar reasoning, we can also formulate the DoA estimation as an
encoding and decoding problem. Analogous to undersampling a signal in time
domain, spatial undersampling indeed also exists when the distance between
adjacent array sensors is larger than a half of wavelength. In general, assume
there are $N$ independent signals $\\{s_{i}(t),i=1,2,...,N\\}$, whose DoAs are
$\bm{\theta}_{[1:N]}=\\{\theta_{1},\theta_{2},...,\theta_{N}\\}$,
respectively. Consider an linear array of $L$ sensors, which are placed at
$\\{p_{1},p_{2},...,p_{L}\\}$, where without loss of generality, $p_{i}\geq 0$
and $p_{1}$ is set to be $0$ as the reference. Thus, the sample sequence at
the $l$-th sensor can be written as

$x_{l}[n]=\sum_{i=1}^{N}e^{-j\frac{2\pi
p_{l}}{\lambda}\sin(\theta_{i})}s_{i}[n].$ (4)

Therefore, let $\bm{x}[n]=(x_{1}[n],x_{2}[n],...,x_{L}[n])^{T}$, and (4) can
be rewritten in a matrix form,

$\bm{x}[n]=\bm{A}\bm{s}[n],$ (5)

where
$\bm{A}=(\bm{a}(\theta_{1}),\bm{a}(\theta_{2}),\dots,\bm{a}(\theta_{N}))$ and
$\bm{s}[n]=(s_{1}[n],\dots,s_{N}[n])^{T}$. $\bm{A}$ is usually called array
manifold matrix, where $\bm{a}(\theta_{i})=[1,e^{-j\frac{2\pi
p_{2}}{\lambda}\sin(\theta_{i})},...,e^{-j\frac{2\pi
p_{L}}{\lambda}\sin(\theta_{i})}]^{T}$, usually called the steering vector.
Here, $\lambda$ represents the wavelength. Moreover, if each sensor takes $K$
snapshots for $n=1,2,...,K$, the samples set $\bm{x}[1:K]=\bm{A}\bm{S}$, where
$\bm{S}=(\bm{s}[1],\bm{s}[2],...,\bm{s}[K])$. We assume $K>N$ in the following
to guarantee sufficient samples. Without any restriction,
$\bm{\theta}_{[1:N]}$ are distributed within $[-\pi/2,\pi/2)^{N}$. Generally,
the distance between any adjacent sensors is set to be $\lambda/2$ for most
uniform arrays. When the distance between adjacent sensors is larger than
$\lambda/2$, there is an ambiguous range for the $l$-th sensor [34]

$\Psi_{\theta_{i}}=[-\frac{\pi}{2},-\arcsin(\frac{\lambda}{p_{l}})]\cup[\arcsin(\frac{\lambda}{p_{l}}),\frac{\pi}{2})$

regarding $\theta_{i}$. In the following, we present a generic
characterization on the condition that $N$ DoAs can be uniquely encoded given
an arbitrary array.

######  Theorem 2.

When $\bm{S}$ is a full-rank matrix, the necessary and sufficient conditions
that $N$ DoAs, $\bm{\theta}_{[1:N]}$, can be uniquely reconstructed by sample
sequences from an array of $L$ sensors with locations specified as
$\\{p_{1}=0,p_{2},...,p_{L}\\}$ are that $\mathcal{C}\geq 2$, where

$\mathcal{C}=\min_{c_{2},...,c_{L}\in\mathbb{Z}^{+}}\frac{c_{2}\lambda}{p_{2}}=\frac{c_{3}\lambda}{p_{3}}=...=\frac{c_{L}\lambda}{p_{L}},$
(6)

is the least common multiple (lcm) of the numbers
$\\{\frac{\lambda}{p_{2}},\frac{\lambda}{p_{3}},...,\frac{\lambda}{p_{L}}\\}.$

###### Proof.

Since $\bm{S}$ is assumed to be of full-rank, the encoding procedure between
$\bm{A}$ and $K$ snapshots $\bm{x}[1:K]$ is invertible. Thus, it is equivalent
to considering the condition that the mapping between $\bm{\theta}_{[1:N]}$
and $\bm{A}$ is bijection, which means for arbitrary
$\theta\not=\theta^{\prime}$, $\bm{a}(\theta)\not=\bm{a}(\theta^{\prime})$. On
the other hand, if $\bm{a}(\theta)=\bm{a}(\theta^{\prime})$, it is equivalent
to that, there exist $k_{l}\in\mathbb{Z}$, $l=1,2,...,L$ and

$2\pi\frac{p_{l}}{\lambda}\sin(\theta_{i})=2\pi(\frac{p_{l}}{\lambda}\sin(\theta^{\prime}_{i})+k_{l}).$
(7)

Rewriting (7) into a modulo form, we have

$\sin(\theta_{i})\equiv\sin(\theta^{\prime}_{i})\mod\frac{\lambda}{p_{l}},l=1,2,...,L.$
(8)

Here, we use the generalized modulo operation where for any two real numbers,
$\langle a\rangle_{b}=a-\lfloor\frac{a}{b}\rfloor b$. Thus, a bijection is
equivalent to that arbitrary $\sin(\theta_{i})\in[-1,1)$ has a unique residue
representation modulo $\frac{\lambda}{p_{l}},l=2,3,...,L.$ Hence, the rest of
the proof is indeed to describe a generalized CRT with real moduli: Indeed,
CRT for real numbers is analogous to the integer case, where the largest
representability is still determined by a generalized least common multiple of
$\\{\frac{\lambda}{p_{l}},l=2,3,...,L\\}$, defined by

$\mathcal{C}=\min_{c_{2},...,c_{L}\in\mathbb{Z}^{+}}\frac{c_{2}\lambda}{p_{2}}=\frac{c_{3}\lambda}{p_{3}}=...=\frac{c_{L}\lambda}{p_{L}}.$

The proof is straightforward. On one hand, it is trivial to verify that the
residues of $a=0$ and $b=\mathcal{C}$ modulo
$\\{\frac{\lambda}{p_{l}},l=2,3,...,L\\}$ are identical. On the other hand,
for any two real numbers $a,b\in[0,\mathcal{C})$, if they share exactly the
same residues, which implies

$a-b\equiv 0\mod\frac{\lambda}{p_{l}},$

for $l=2,3,...,L$, a contradiction to the minimization assumption of
$\mathcal{C}$. Finally, it is noted that $\sin(\theta)\in[-1,1)$, therefore we
need $\mathcal{C}$ to be at least $2$ and the claim follows. ∎

Analogous to the DoA estimation, in the case of phase unwrapping (Doppler
shift) in multi-frequency antenna array synthetic aperture radar (SAR), the
condition is also similar to Theorem 2.

Throughout the two examples, we have shown that the fundamental problems to
characterize the frequency and DoA estimation under sparse sampling are
closely related to the remainder model. This is not surprising due to the
cycle property of waveform, where sampling itself can be viewed as somewhat
modulo operation on the signal. Though the necessary and sufficient conditions
describe a preliminary picture that when a signal can be perfectly
reconstructed given samplers with specific sparsity restriction, there are two
important generalizations worth further exploring:

  1. 1.

How could restrictions in Theorem 1 and 2 be further relaxed resorting to some
failure probability? Beyond perfect reconstruction on all signals, compressed
sensing provides a successful example to produce high probability recovery but
much sharpened sample complexity via random sampling.

  2. 2.

Theorem 1 and 2 only characterize the noiseless case while the noise and
quantification error are almost inevitable in practice. With further
robustness consideration, how should Theorem 1 and 2 be modified.

We address the above two questions in the following.

##  III A Complete Analysis of Co-prime Sampling

In the Introduction, we briefly mentioned the key idea of co-prime sampling
[18, 35]. The elegant construction by selecting two co-prime under-sampling
rates, $P$ and $Q$, gives a unified approach to estimate the auto-correlation
of the sample sequence attained by Nyqiust rate. Due to Bezout theorem, one
can properly select infinitely many pairs from two undersampled sequences such
that their sampling-time difference to be any multiple of the time interval
sampled by Nyquist rate. Once one can exactly recover the auto-correlation of
the sequence sampled at Nyquist rate, according to Wiener–Khinchin theorem,
the power spectrum can be reconstructed by applying Fourier transform on the
auto-correlation.

Co-prime sampling theory seems to violate the necessary and sufficient
condition described in Theorem 1, since $P$ and $Q$ can be arbitrarily large,
i.e., one can use two arbitrarily small sampling rates to recover the spectrum
of a signal with a arbitrarily large bandwidth. As shown below, we will
rigorously prove co-prime sampling is indeed an efficient “high-probability”
reconstruction but a failure could happen.

Assuming that two samplers whose sampling rates are $m_{1}=1/{PT}$ and
$m_{2}=1/{QT}$, respectively, where the Nyquist rate is $1/{T}$. With the same
setup described in Section II-A, $x(t)=\sum_{i=1}^{N}e^{j2\pi f_{i}t}$, where
we normalize the phase and amplitude without loss of generality, is sampled
and two sequences are attained as

$x_{1}[n_{1}]=\sum_{i=1}^{N}e^{j2\pi f_{i}TPn_{1}},\leavevmode\nobreak\
\leavevmode\nobreak\ x_{2}[n_{2}]=\sum_{i=1}^{N}e^{j2\pi f_{i}TQn_{2}}.$ (9)

Since $P$ and $Q$ are co-prime integers, the two sequences indeed form a cycle
of time $PQT$. In the $k$-th cycle, i.e., we consider the sampling time at
$n_{1}=Qk+r_{1}$ and $n_{2}=Pk+r_{2}$, and denote the sample vectors by
$\bm{x}^{1}_{k}=(x_{1}[Qk+1],x_{1}[Qk+2],...,x_{1}[Q(k+1)])$ and
$\bm{x}^{2}_{k}=(x_{2}[Pk+1],x_{2}[Pk+2],...,x_{2}[P(k+1)])$. Following [18,
35], let $\tilde{\bm{x}}_{k}=(\bm{x}^{1}_{k},\bm{x}^{2}_{k})$, be the
Cartesian product of the two vectors $\bm{x}^{1}_{k}$ and $\bm{x}^{2}_{k}$,
and consider
$\bm{R}=\mathbb{E}_{l}\big{\\{}(\tilde{\bm{x}}^{l})^{T}\cdot(\tilde{\bm{x}}^{l})^{*}\big{\\}}$,
which ideally is expected to be the auto-correlation estimation. Here,
$\bm{a}^{T}$ and $\bm{a}^{*}$ denotes the transpose and conjugate of a vector
$\bm{a}$, respectively.

When $N=1$, the co-prime sampling theory is indeed a special case of Theorem
1. It is noted that

$lcm(\frac{1}{PT},\frac{1}{QT})=\big{(}\min_{c_{1},c_{2}\in\mathbb{Z}^{+}}\frac{c_{1}}{PT}=\frac{c_{2}}{QT}\big{)}=\frac{1}{T},$
(10)

where the quality holds when $c_{1}=P$ and $c_{2}=Q$ due to the co-prime
assumption on $P$ and $Q$.

However, when $N>1$, co-prime sampling comes with a failure, where the
estimation $\bm{R}$ could be biased. Consider
$\mathbb{E}_{k}\\{x_{1}[Qk+u]\cdot x_{2}[Pk+v]\\}$, which is took as an
estimation of auto-correlation at lag $(uP-vQ)$, i.e.,
$\mathbb{E}_{n}\big{\\{}x[nT]\cdot x^{*}[(n-(uP-vQ))T]\big{\\}}$. With some
simple calculation,

$\displaystyle\mathbb{E}_{k}$ $\displaystyle\\{x_{1}[Qk+u]\cdot
x^{*}_{1}[Pk+v]\\}$ (11) $\displaystyle=\mathbb{E}_{n}\\{x[nT]\cdot
x^{*}[(n-(uP-vQ))T]\\}$
$\displaystyle+\underline{\sum_{i\not=l}\mathbb{E}_{k}\big{[}e^{j2\pi
f_{i}P(Qk+u)T}\cdot e^{-j2\pi f_{l}Q(Pk+v)T}\big{]}}$

Similarly, if the estimation is from the cross product from one sample
sequence, say $x_{1}[\cdot]$, where we expect
$\mathbb{E}_{k}\big{\\{}x_{1}[kQ+u]\cdot x_{1}[kQ+v]\big{\\}}$ equals the
auto-correlation at lag $(u-v)P$, i.e., $\mathbb{E}_{n}\big{\\{}x[nT]\cdot
x^{*}[(n-(u-v)P)T]\big{\\}}$, it becomes

$\displaystyle\mathbb{E}_{k}$ $\displaystyle\\{x_{1}[Qk+u]\cdot
x^{*}_{1}[Pk+v]\\}$ (12) $\displaystyle=\mathbb{E}_{n}\\{x[nT]\cdot
x^{*}[(n-(u-v)P)T]\\}$
$\displaystyle+\underline{\sum_{i\not=l}\mathbb{E}_{k}\big{[}e^{j2\pi
f_{i}P(Qk+u)T}\cdot e^{-j2\pi f_{l}P(Qk+v)T}\big{]}}$

Thus, co-prime sampling forms an unbiased estimation if and only if the
underlying terms in (11) and (12) should equal 0. With a bit more calculation,
a sufficient condition is that,

$\mathbb{E}_{k}e^{j2\pi k\frac{PQ(f_{i}-f_{l})}{f_{s}}}=0$ (13)

which requires that $PQ(f_{i}-f_{l})T$ cannot be an integer.

######  Theorem 3.

When $N>2$, the co-prime sampler defined in (9) can perfectly reconstruct the
$N$ frequencies $\\{f_{1},f_{2},...,f_{N}\\}$, when for any $i\not=j\in[1:N]$,
$PQ(f_{i}-f_{j})$ cannot be divided by the Nyquist rate $1/{T}$. 333Another
takeaway from the above analysis is that, in general by applying the co-prime
sampling idea on generic continuous signal, the assumption that the components
are independent is not sufficient since after sampling (processing), the
independence may not hold.

Theorem 3 characterizes a failure set, within which the signal cannot be
uniquely reconstructed. However, it is worth noting that such set is very
sparse (indeed has zero measure) across the whole domain $[0,1/T]^{N}$.
Compared to the best known dynamic range in CRT decoding method [36], which
states that given $L$ samplers with co-prime rates
$\\{m_{1},m_{2},...,m_{L}\\}$, one can uniquely recover $\bm{f}\in[0,D]^{N}$
for $D=\Omega\big{(}(\prod_{l=1}^{L}m_{l})^{1/N}\big{)}$. Theorem 3 is
inspiring and gives a positive example that, resorting to a negligible failure
probability (when $\bm{f}$ is assumed to be uniformly distributed, or a random
selection of $P$ and $Q$), it is possible that one can use much sharpened
undersampling rate to almost reconstruct the signal, compared to the
deterministic decoding bound.

##  IV Robust Reconstruction for Frequency Estimation of Complex Waveform

Throughout the previous sections, we have derived the necessary and sufficient
reconstruction condition in the noiseless case. In the rest of this paper, we
set out to advance the understanding of robust reconstruction and stick to the
more complicated (multiple) frequency estimation model. We adopt the noise
model considered in [20, 23, 24, 37, 31] and study the robustness from an
error-correcting decoding viewpoint.

As described in Theorem 1, under $L$ sampling rates $m_{[1:L]}$, the frequency
set $\bm{f}=(f_{1},f_{2},...,f_{N})$ to be estimated is essentially encoded
into $L$ unordered residue sets $\mathcal{R}_{[1:L]}$, where
$\mathcal{R}_{l}=\\{\langle f_{i}\rangle_{m_{l}},i=1,2,...,N\\}$. To formally
study the robust decoding condition, we introduce the robust remaindering
problem as follows.

######  Definition 1 (Robust Remaindering Problem).

Given $m_{[1:L]}$ and $N$, what is the relationship between the dynamic range
$D$ and the maximal error bound $\delta$ such that for arbitrary
$\bm{X}=(X_{1},X_{2},...,X_{N})\in[0,D)^{N}$, provided $L$ noisy unordered
residue sets $\mathcal{R}_{[1:L]}$, where $\mathcal{R}_{l}=\\{\langle
X_{i}\rangle_{m_{l}}+\Delta_{il},i=1,2,...,N\\}$ for arbitrary noise
$|\Delta_{il}|<\delta$, $\bm{X}$ can be robustly estimated by
$\hat{\bm{X}}=\\{\hat{X}_{1},\hat{X}_{2},...,\hat{X}_{N}\\}$ and the
estimation error $|X_{i}-\hat{X}_{i}|=O(\delta)$ for any $i$.

Before we proceed to study the tradeoff between $D$ and $\delta$ in general,
we first provide the following results on the sampling rate selection. We
first move our attention to a single number case, i.e., $N=1$. In the
following, we prove the optimal modulus selection in terms of the error
tolerance capacity.

###  IV-A Sampling Rate Selection

In [36], it is shown that when $N=1$, the maximal error bound is,

$\delta=\min\max_{\min m_{l}\leq X<D}\frac{\|{\textbf{x}}\|_{\infty}}{4}$ (14)

Here, x denotes the residue representation of a number $X$ and
$\|{\textbf{x}}\|_{\infty}$ is the largest coordinate. In the following, we
show a corollary concerning the lower bound of $\delta$.

######  Corollary 1.

When $D=lcm(m_{1},m_{2},...,m_{L})$,

$\delta\leq\min_{S\subset\\{1,2,...,L\\}}gcd(lcm(m_{l},l\in
S),lcm(m_{l},l\in\bar{S}))$ (15)

###### Proof.

Let us consider a class of residue vectors, in which the residues modulo
$m_{l}$, for $l\in S$, are set to be 0 and the others are set to be some
$a<\min m_{l}$. If such a residue vector exists, $X$ can be represented by
$X=k_{1}lcm(m_{l},l\in S)=k_{2}lcm(m_{l},l\in\bar{S})+a$. Since
$D=lcm(lcm(m_{l},l\in S),lcm(m_{l},l\in\bar{S}))$, the existence of such $X$
means that the following Diophantine equation has a solution,

$k_{1}lcm(m_{l},l\in S)-k_{2}lcm(m_{l},l\in\bar{S})=a.$ (16)

It is further equivalent to $gcd(lcm(m_{l},l\in S),lcm(m_{l},l\in\bar{S}))|a$.
Therefore, our claims hold. ∎

Clearly, when $M_{l}={m_{l}}/{\Gamma}$ are relatively co-prime to each other,
the error bound is $\delta={\Gamma}/{4}$. Correspondingly,
$D=\Gamma\prod_{l=1}^{L}M_{l}$. On the other hand, consider that the moduli
are not in such form, while we still expect to achieve a same error bound
$\delta={\Gamma}/{4}$. From Corollary 1, for any two disjoint partitions of
$\mathcal{M}=m_{[1:L]}$, the lcm of each part must share a factor no less than
$\Gamma=4\delta$. Now, we select $S$ to contain only a single modulus, which
from Corollary 1 indicates that each modulus $m_{l}$ must share a common
factor no less than $\Gamma$ with the lcm of the rest moduli. With the
consideration that $D=lcm(m_{1},m_{2},...,m_{L})$, if such a common factor is
not the same, more redundancy is needed compared with the case
$M_{l}={m_{l}}/{\Gamma}$ and $D=\Gamma\prod_{l=1}^{L}M_{l}$.

###  IV-B Error Tolerance Bound and Reconstruction Scheme

In this section, we study the relationship between robustness and the encoding
dynamic range via a reconstruction of a robust decoding algorithm. Assume that
$\bm{X}=\\{X_{1},X_{2},...,X_{N}\\}$ are $N$ distinct real numbers within the
dynamic range $[0,D)$ to be determined. Suppose a modulus set
$\mathbb{\mathcal{M}}=\\{m_{1},m_{2},...,m_{L}\\}$ is given, where
$m_{l}=\Gamma M_{l}$ with pairwise coprime $M_{l}$, and without loss of
generality $M_{[1:L]}$ are assumed to be in an ascending order, i.e.,
$M_{1}<M_{2}<...<M_{L}$. We use $\mathcal{R}_{l}=\\{\widetilde{r}_{il}=\langle
X_{i}+\Delta_{il}\rangle_{m_{l}},i=1,2,...,N\\}$ to denote the noisy
observations (encoding) in each modulus $m_{l}\in\mathcal{M}$. Moreover,
$\widetilde{r}^{c}_{il}=\langle\widetilde{r}_{il}\rangle_{\Gamma}$ denotes a
noisy observation of the common residue $r^{c}_{i}=\langle
X_{i}\rangle_{\Gamma}.$

######  Theorem 4.

For any
$\bm{X}\in\big{[}0,\Gamma(\prod_{l=1}^{\lceil\frac{L}{N}\rceil}M_{l}-1)\big{)}^{N}$,
if the errors underlying satisfy
$\max_{il}|\Delta_{il}|<\delta=\frac{\Gamma}{4}$, there exists a robust
decoding algorithm returning an estimation $\hat{\bm{X}}=\hat{X}_{[1:N]}$ such
that $|\hat{X}_{i}-X_{i}|<\frac{3\Gamma}{4}$ for each $i$.

We present a sketch of the proof, where the details can be found in Appendix
A. The decoding algorithm can be found in Algorithm 1. Roughly speaking, the
main procedure can be summarized into two steps: hypothesis and test. As
mentioned before, one key challenge in decoding is the unknown correspondence.
Therefore, given an arbitrary clustering proposal by dividing the residues
into $N$ subsets, a verification should exist to check whether any residue is
misclassified 444It is possible that different $X_{i}$ and $X_{j}$, or their
folding numbers $\lfloor X_{i}/\Gamma\rfloor$ and $\lfloor
X_{j}/\Gamma\rfloor$, share the same residue. We will take those residues
interchangeable in our analysis.. To complete the whole reconstruction,
besides a correct guess of clustering, another statistic needs to be recovered
is the order of errors: the ascending order of $\Delta_{il}$. It is noted that
modulo operation results in non-weighted residues with ambiguity. For example,
under noise perturbation, when we observe a noisy residue 11 modulo 12, it can
either be $13-2$ or $9+2$. This raises another key challenge in the robust
decoding that the empirical average of residues is not robust. However, the
ascending order of noises $\Delta_{il}$ enables us to remove the modulo
operation and simply represent the noisy residue over the real axis instead.
Continued to the earlier example, for some number $X$ with a true residue $1$
modulo 12, where we have two observations $11=\langle 1-2\rangle_{12}$ and
$2=\langle 1+1\rangle_{12}$. A simple average of $2$ and $11$, which gives an
estimate $6.5$, produces an estimation error larger than a quarter of the
modulus $12$ and we will show such trivial estimation cannot give a satisfied
error control at the end. Nonetheless, given the order of errors, where we
represent the noisy observation in the following equivalent form $(11-1\times
12)$ and $(2-0\times 12)$, then a trivial CRT is indeed robust already. In
Algorithm 1, the proposed verification is formed by two subroutines and we
show that if the proposal passes through the two validity tests, then it can
produce desired error control.

With the above understanding, when $\delta=\Gamma/4$, a correct proposal on
both clustering and the order of noise can be essentially modeled by
$S_{i}=\\{\tilde{r}_{(il)},l=1,2,...,L\\}$ and $\tau_{(il)}\in\\{0,1\\}$,
where $(il)$ denotes the index given the clustering, such that $S_{i}$ are the
residues of $X_{i}$ and meanwhile
$\\{\langle\tilde{r}_{(il)}\rangle_{\Gamma}-\tau_{(il)}\Gamma,l=1,2,...,L\\}$
in an ascending order corresponds to that of
$\\{\Delta_{(il)},l=1,2,...,L\\}$. The proof heavily relies on the following
two facts:

  1. 1.

Given the correct proposal of clustering and $\tau_{(il)}$, the shifted common
residues satisfy

$|(\tilde{r}^{c}_{(il)}-\tau_{(il)})-(\tilde{r}^{c}_{(il^{\prime})}-\tau_{(il^{\prime})})|<{\Gamma}/{2}.$

  2. 2.

Based on the pigeonhole principle, for any proposal of clustering, there exist
at least $\lceil\frac{L}{N}\rceil$ from the $L$ residues are all from one
integer $X_{i_{0}}$.

The first criterion is stemmed from the assumption that the magnitude of error
is smaller than $\Gamma/4$, and thus the distance between any two noisy
observations of the common residue should not be larger than $\Gamma/2$, after
they have been sorted according to the order of errors.

Algorithm 1 Robust Remaindering Decoding of Complex Waveform

Input: modulus set: $\mathcal{M}=\\{m_{l}=M_{l}\Gamma|l=1,2,\dots,L\\}$, where
$M_{[1:L]}$ are sorted in ascending order;
Residue Sets: $\mathcal{R}_{l}=\\{\widetilde{r}_{il}|i=1,2,\dots,N\\}$,
$l=1,2,\dots,L$.

1: Repeat: Propose a clustering assignment

2: Following the proposed clustering by selecting corresponding residue from
each $\mathcal{R}_{l}$ to obtain a $L$-residue clustering
$S_{i}=\\{\widetilde{r}_{(il)}|l=1,2,\dots,L\\}$, $i=1,2,...,N$. Here, $(il)$
denotes the index of the residues assigned to $S_{i}$.

3: Assign a binary parameter $\tau_{(il)}\in\\{0,1\\}$ to each residue in
$S_{i}$ randomly.

4: Calculate the shifted common residues
$\hat{r}^{c}_{(il)}=\langle\widetilde{r}_{(il)}\rangle_{\Gamma}-\tau_{(il)}\Gamma$
corresponding to each $\widetilde{r}_{(il)}$ in $S_{i}$, for $i=1,2,...,N$.

5: For each $S_{i}$, calculate
$q_{i}\equiv\frac{\widetilde{r}_{(il)}-\hat{r}^{c}_{(il)}}{\Gamma}\mod M_{l}$
via CRT.

6: Until: Each $q_{i}$ satisfies that $q_{i}\in[0,D_{q})$ where
$D_{q}=\prod_{l=1}^{\lceil\frac{L}{N}\rceil}M_{l}\Gamma$ and each residue in
$S_{i}$ satisfies equation (23) for $i=1,2,\dots,N$.

Output: $\hat{X_{i}}=q_{i}\Gamma+{\sum^{L}_{l=1}\hat{r}^{c}_{(il)}}/{L}$.

Before the end of this section, we give an interesting corollary that,
excluding a negligible failure set, the reconstruction error bound can be
further improved from $3\Gamma/4$ to $\Gamma/4$.

######  Corollary 2.

For any
$X_{i}\in[0,\Gamma(\prod^{\lceil\frac{L}{N}\rceil}_{l=1}M_{l}-1))^{N}$, with
the errors introduced in residues satisfy
$\max_{il}|\Delta_{il}|<\delta=\frac{\Gamma}{4}$, if $\min_{d=0,\pm
1}|r_{i_{0}l}-r_{i_{1}l}+dM_{l}\Gamma|>3\Gamma$ holds for $l=1,2,\dots L$ and
$1\leq i_{0}<i_{1}\leq N$, $X_{i}$ can be recovered robustly with error
bounded by $\frac{\Gamma}{4}$.

##  V Generalization to Real Waveforms

###  V-A Necessary and Sufficient Condition

Slightly different from the complex waveform, a sinusoidal real waveform has
two complex components,

$x(t)=A_{i}\cos(2\pi f_{i}t)=A_{i}e^{2\pi jf_{i}t}+A_{i}e^{-2\pi jf_{i}t},$
(17)

which results in two symmetric peaks over the spectrum domain, i.e., with an
additional negative duplicate compared to the complex waveform case. With a
same reasoning, the sufficient and necessary condition to robustly recover the
frequency(ies) in a real waveform can be similarly modeled as the encoding of
$N$ real numbers $\\{X_{1},X_{2},\dots,X_{N}\\}$ by $L$ residue sets
$\mathcal{R}_{l}=\\{\widetilde{r}^{+}_{i,l}=\langle
X_{i}+\Delta^{+}_{i,l}\rangle_{m_{l}},\widetilde{r}^{-}_{i,l}=\langle-
X_{i}+\Delta^{-}_{i,l}\rangle_{m_{l}}|i=1,2,\dots,N\\},$ for $l=1,2,\dots,L$,
where $\Delta^{\pm}_{i,l}$ represents the error. We still assume that
$m_{[1:L]}$ are in an ascending order, where $m_{l}=\Gamma M_{l}$ with co-
prime $M_{[1:L]}$. To make the problem nontrivial, it is assumed that
$m_{l}<2\max X_{i}$.

We first consider the noiseless real waveform model with a single frequency,
i.e., ${R}_{l}(X)=\\{r^{+}_{l}=\langle
X\rangle_{m_{l}},r^{-}_{l}=\langle-X\rangle_{m_{l}}\\}$, is given for each
sampling rate $m_{l}$. For a set of moduli, $\\{m_{1},m_{2},...,m_{L}\\}$, its
maximum dynamic range is the maximum $D$ such that any $X\in[0,D)$ has a
unique residue set representation $\mathcal{R}_{l}(X)$. In this section, we
provide a closed-form equation for the maximum dynamic range under the
errorless real waveform model.

Let us suppose that $D$ is the maximum dynamic range given
$\mathcal{M}=\\{m_{1},m_{2},...,m_{L}\\}$. By definition, there must exist
another $Y\in[0,D)$ such that $\mathcal{R}_{l}(D)=\mathcal{R}_{l}(Y)$ for
$l\in[1:L]$. Therefore, for each $m_{l}$, $D$ must satisfy at least one of the
following equations,

$D+Y\equiv 0\mod m_{l}$ (18)

or

$D-Y\equiv 0\mod m_{l}$ (19)

Let $\mathbb{U}_{1}$ and $\mathbb{U}_{2}$ be two subsets of
$\\{m_{1},m_{2},...,m_{L}\\}$, where $\mathbb{U}_{1}$ and $\mathbb{U}_{2}$
denote the set of the moduli satisfying (18) and (19), respectively. The two
sets might have an intersection and their union is
$\\{m_{1},m_{2},...,m_{L}\\}$, i.e.,
$\mathbb{U}_{1}\cup\mathbb{U}_{2}=\\{m_{1},m_{2},...,m_{L}\\}$. This yields
that

$\left\\{\begin{array}[]{lr}D+Y\equiv 0\mod lcm(m_{l}\in\mathbb{U}_{1})\\\
D-Y\equiv 0\mod lcm(m_{l}\in\mathbb{U}_{2})\end{array}\right..$

Since $D>Y>0$, $D+Y$ and $D-Y$ are both non-zero. Thus, $D+Y\geq
lcm(m_{l}\in\mathbb{U}_{1})$ and $D-Y\geq lcm(m_{l}\in\mathbb{U}_{2})$, which
implies that

$D\geq\frac{lcm(m_{l}\in\mathbb{U}_{1})+lcm(m_{l}\in\mathbb{U}_{2})}{2}.$ (20)

If we let
$D=\frac{lcm(m_{l}\in\mathbb{U}_{1})+lcm(m_{l}\in\mathbb{U}_{2})}{2}$ and
$Y=|\frac{lcm(m_{l}\in\mathbb{U}_{1})-lcm(m_{l}\in\mathbb{U}_{2})}{2}|$,
$\mathcal{R}_{l}(D)$ is exactly the same as $\mathcal{R}_{l}(Y)$. It means
that equation (20) is satisfied with a tight lower bound of $D$. Notice that
$\mathbb{U}_{1}$ and $\mathbb{U}_{2}$ can be arbitrary two sets that the union
of both is $\\{m_{1},m_{2},...,m_{L}\\}$. Thus, finding the maximum dynamic
range is equivalent to finding a proper subset
$\mathbb{U}\subset\\{m_{1},m_{2},...,m_{L}\\}$ to minimize

$D=\min_{\mathbb{U}}\frac{lcm(m_{l}\in\mathbb{U})+lcm(\\{m_{1},m_{2},...,m_{L}\\}/\mathbb{U})}{2}.$
(21)

In the following, we give the robustness bound for the single frequency
estimation.

######  Theorem 5.

If the error introduced in residues such that
$\max|\Delta^{\pm}_{l}|<\delta=\frac{\Gamma}{4}$, $X$ can be recovered
robustly error bounded by $\frac{3\Gamma}{4}$ for any $X\in[0,D)$, where 555If
we assume that the common residue $r_{c}=\langle
X\rangle_{\Gamma}\not=0\leavevmode\nobreak\ \text{ or }\leavevmode\nobreak\
\frac{\Gamma}{2}$, the maximum dynamic range is equivalent to
$lcm(m_{1},m_{2},...,m_{L})$, where $\Gamma$ represents
$gcd(m_{1},m_{2},...,m_{L})$. Such claim is based on the fact that for all $X$
satisfying our assumptions, its positive and negative residues $r^{+}_{c,l}$
and $r^{-}_{c,l}$ are not the same, i.e., $\langle
r^{+}_{c,l}\rangle_{\Gamma}=\langle
X\rangle_{\Gamma}\not=\langle-X\rangle_{\Gamma}=\langle
r^{-}_{c,l}\rangle_{\Gamma}.$ This allows us to distinguish $r^{+}_{c,l}$ from
$r^{-}_{c,l}$ and easily recover $X$. When the errors exist in our model, if
the magnitudes of errors are smaller than $\frac{\min\\{\langle
2X\rangle_{\Gamma},\Gamma-\langle 2X\rangle_{\Gamma}\\}}{2},$ $r^{+}_{c,l}$
and $r^{-}_{c,l}$ can still be distinguished.

$D\leq[0,\min_{\mathbb{U}\subset\\{M_{1},M_{2},\dots,M_{L}\\}}(\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1)\Gamma).$

The robust decoding algorithm and proof can be found in Appendix C and D,
respectively. We have to emphasize that Algorithm 2 proposed can be further
improved to be of linear decoding time with the trick from Generalized CRT
[38, 36] with a slight compromise in dynamic range, where $D$ becomes

$D\leq\Gamma(\prod_{l=1}^{\lceil\frac{L}{2N}\rceil}M_{l}-1).$

###  V-B Multiple Frequencies Estimation from Real Waveforms

In this section, we will present a generic characterization of the robustness
and encoding dynamic range in the real waveform model with $N$ frequencies.
With a similar reasoning, the corresponding model can be described as follows.
Given the modulus set $\mathcal{M}=\\{m_{l}=M_{l}\Gamma|l=1,2,\dots,L\\}$,
where $M_{l}$ are pairwise coprime and sorted in ascending order, we aim to
recover $N$ distinct real numbers $\bm{X}=\\{X_{1},X_{2},...,X_{N}\\}$ within
the dynamic range $[0,D)$ with $L$ residue sets
$\mathcal{R}_{l}=\\{\widetilde{r}^{+}_{i,l}=\langle
X_{i}+\Delta^{+}_{i,l}\rangle_{m_{l}},\widetilde{r}^{-}_{i,l}=\langle-
X_{i}+\Delta^{-}_{i,l}\rangle_{m_{l}}|i=1,2,...,N\\}.$ Here,
$\Delta^{\pm}_{i,l}$ similarly captures the noises underlying.

Indeed, the robust decoding problem in the real waveform case can be
analogously addressed with the similar idea applied in the complex waveform
case. Once the noisy residues can be correctly clustered, the $N$-number
robust decoding problem can be reduced to $N$ independent single-number
decoding, which has been solved in Section V-A. The essential challenge is
that now each objective number is encoded by two symmetrically positive and
negative residues with noises in each sampler, which makes decoding more
complicated.

We still apply the hypothesis-then-testing framework proposed in Section IV,
but in the real waveform case we view the $N$ real numbers with positive and
negative copies as $2N$ numbers. Thus, the residues are clustered into $2N$
subsets, still denoted by $S_{i}=\\{\widetilde{r}_{(i,l)}|l=1,2,\dots,L\\}$,
for $i=1,2,...,2N$. The two criteria can be similarly derived.

######  Theorem 6.

In the real waveform case, if the errors introduced in residues satisfy
$\max_{i,l}|\Delta^{\pm}_{i,l}|<\delta=\frac{\Gamma}{4}$, $X_{i}$ can be
reconstructed error bounded by $\frac{3\Gamma}{4}$ for any
$X_{i}\in[0,D)^{N}$, where
$D=\min(\frac{lcm(M_{l}\in\mathbb{U})+lcm(M_{l}\in\overline{\mathbb{U}})}{2}-1)\Gamma$
and
$\mathbb{U}\cup\overline{\mathbb{U}}=\\{M_{1},M_{2},\dots,M_{\lceil\frac{L}{N}\rceil}\\}$.

##  VI Conclusion

From a (robust) encoding-then-decoding viewpoint, this paper advances the
understanding of sampling theory from a new perspective: given the sparsity
restriction of a set of samplers temporally or spatially, what is the
necessary and sufficient condition of the signal reconstruction from acquired
samples? We characterize and prove such conditions in several common parameter
estimation tasks such as frequency, phase difference, DoA and Doppler shift.
One main takeaway is that, with a random selection of sampling parameter and a
negligible failure probability, it is possible that one can asymptotically
sharpen the deterministic necessary sampling constraint. We leave a
systematical study on such improvement to our future work.

## References

  * [1] H. Nyquist, “Abridgment of certain topics in telegraph transmission theory,” Journal of the A.I.E.E., vol. 47, no. 3, pp. 214–217, 1928.
  * [2] C. Shannon, “Communication in the presence of noise,” Proceedings of the IEEE, vol. 72, no. 9, pp. 1192–1201, 1984.
  * [3] H. Landau, “Necessary density conditions for sampling and interpolation of certain entire functions,” Acta Mathematica, vol. 117, no. 1, pp. 37–52, 1967.
  * [4] D. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
  * [5] E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006.
  * [6] E. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Communications on Pure and Applied Mathematics, vol. 59, no. 8, pp. 1207–1223, 2010.
  * [7] E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 21–30, 2008.
  * [8] E. J. Candes and T. Tao, “Decoding by linear programming,” IEEE transactions on information theory, vol. 51, no. 12, pp. 4203–4215, 2005.
  * [9] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse mri: The application of compressed sensing for rapid mr imaging,” Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007.
  * [10] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing mri,” IEEE signal processing magazine, vol. 25, no. 2, pp. 72–82, 2008\.
  * [11] J. P. Haldar, D. Hernando, and Z.-P. Liang, “Compressed-sensing mri with random encoding,” IEEE transactions on Medical Imaging, vol. 30, no. 4, pp. 893–903, 2010.
  * [12] R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Optical Engineering, vol. 50, no. 7, p. 072601, 2011.
  * [13] A. Kipnis, A. J. Goldsmith, Y. C. Eldar, and T. Weissman, “Distortion rate function of sub-nyquist sampled gaussian sources,” IEEE transactions on information theory, vol. 62, no. 1, pp. 401–429, 2015.
  * [14] A. Kipnis, Y. C. Eldar, and A. J. Goldsmith, “Fundamental distortion limits of analog-to-digital compression,” IEEE Transactions on Information Theory, vol. 64, no. 9, pp. 6013–6033, 2018.
  * [15] A. Kipnis, Y. C. Eldar, and A. J. Goldsmith, “Analog-to-digital compression: A new paradigm for converting signals to bits,” IEEE Signal Processing Magazine, vol. 35, no. 3, pp. 16–39, 2018.
  * [16] X.-G. Xia, “On estimation of multiple frequencies in undersampled complex valued waveforms,” IEEE Transactions on Signal Processing, vol. 47, no. 12, pp. 3417–3419, 1999.
  * [17] X.-G. Xia, “An efficient frequency-determination algorithm from multiple undersampled waveforms,” IEEE Signal Processing Letters, vol. 7, no. 2, pp. 34–37, 2000.
  * [18] P. Pal and P. P. Vaidyanathan, “Nested arrays: A novel approach to array processing with enhanced degrees of freedom,” IEEE Transactions on Signal Processing, vol. 58, no. 8, pp. 4167–4181, 2010.
  * [19] P. P. Vaidyanathan and P. Pal, “Sparse sensing with co-prime samplers and arrays,” IEEE Transactions on Signal Processing, vol. 59, no. 2, pp. 573–586, 2011.
  * [20] X.-G. Xia and G. Wang, “Phase unwrapping and a robust chinese remainder theorem,” IEEE Signal Processing Letters, vol. 14, no. 4, pp. 247–250, 2007.
  * [21] G. Wang, X.-G. Xia, V. C. Chen, and R. Fielder, “Detection, location, and imaging of fast moving targets using multifrequency antenna array sar,” IEEE Transactions on Aerospace and Electronic Systems, vol. 40, no. 1, pp. 345–355, 2004.
  * [22] J. Xu, Z.-Z. Huang, Z.-R. Wang, L. Xiao, X.-G. Xia, and T. Long, “Radial velocity retrieval for multichannel sar moving targets with time–space doppler deambiguity,” IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 1, pp. 35–48, 2017.
  * [23] W. Wang and X.-G. Xia, “A closed-form robust chinese remainder theorem and its performance analysis,” IEEE Transactions on Signal Processing, vol. 58, no. 11, pp. 5655–5666, 2010.
  * [24] L. Xiao, X.-G. Xia, and W. Wang, “Multi-stage robust chinese remainder theorem,” IEEE Transactions on Signal Processing, vol. 62, no. 18, pp. 4772–4785, 2014.
  * [25] L. Xiao, X.-G. Xia, and Y.-P. Wang, “Exact and robust reconstructions of integer vectors based on multidimensional chinese remainder theorem (md-crt),” IEEE Transactions on Signal Processing, vol. 68, pp. 5349–5364, 2020.
  * [26] H. Liao and X.-G. Xia, “A sharpened dynamic range of a generalized chinese remainder theorem for multiple integers,” IEEE Transactions on Information Theory, vol. 53, no. 1, pp. 428–433, 2007.
  * [27] L. Xiao, X.-G. Xia, and H. Huo, “Towards robustness in residue number systems,” IEEE Transactions on Signal Processing, vol. 65, no. 6, pp. 1497–1510, 2016.
  * [28] L. Xiao, X.-G. Xia, and H. Huo, “New conditions on achieving the maximal possible dynamic range for a generalized chinese remainder theorem of multiple integers,” IEEE Signal Processing Letters, vol. 22, no. 12, pp. 2199–2203, 2015.
  * [29] W. Wang, X. Li, X.-G. Xia, and W. Wang, “The largest dynamic range of a generalized chinese remainder theorem for two integers,” IEEE Signal Processing Letters, vol. 22, no. 2, pp. 254–258, 2015.
  * [30] H. Xiao, C. Cremers, and H. K. Garg, “Symmetric polynomial & crt based algorithms for multiple frequency determination from undersampled waveforms,” in 2016 IEEE Global Conference on Signal and Information Processing (GlobalSIP), 2016.
  * [31] H. Xiao, Y. Huang, Y. Ye, and G. Xiao, “Robustness in chinese remainder theorem for multiple numbers and remainder coding,” IEEE Transactions on Signal Processing, vol. 66, no. 16, pp. 4347–4361, 2018.
  * [32] H. Xiao and G. Xiao, “On solving ambiguity resolution with robust chinese remainder theorem for multiple numbers,” IEEE Transactions on Vehicular Technology, vol. 68, no. 5, pp. 5179–5184, 2019.
  * [33] H. Xiao, N. Du, Z. Wang, and G. Xiao, “Wrapped ambiguity gaussian mixed model with applications in sparse sampling based multiple parameter estimation,” Signal Processing, vol. 179, p. 107825, 2021.
  * [34] J. Benesty, I. Cohen, and J. Chen, Fundamentals of signal enhancement and array signal processing.  John Wiley & Sons, 2017.
  * [35] P. Pal and P. P. Vaidyanathan, “Coprime sampling and the music algorithm,” in 2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE), pp. 289–294, IEEE, 2011.
  * [36] H. Xiao, Y. Huang, Y. Ye, and G. Xiao, “Robustness in chinese remainder theorem for multiple numbers and remainder coding,” IEEE Transactions on Signal Processing, vol. 66, no. 16, pp. 4347–4361, 2018.
  * [37] X. Li, X.-G. Xia, W. Wang, and W. Wang, “A robust generalized chinese remainder theorem for two integers,” IEEE Transactions on Information Theory, vol. 62, no. 12, pp. 7491–7504, 2016.
  * [38] L. Xiao and X.-G. Xia, “A generalized chinese remainder theorem for two integers,” IEEE Signal Processing Letters, vol. 21, no. 1, pp. 55–59, 2013\.

##  Appendix A Proof of Theorem 4

###### Proof.

With the fact that
$\lfloor\frac{X_{i}}{\Gamma}\rfloor\equiv\lfloor\frac{\langle
X_{i}\rangle_{m_{l}}}{\Gamma}\rfloor\equiv\lfloor\frac{\langle
X_{i}\rangle_{m_{l}}-\langle
X_{i}\rangle_{\Gamma}}{\Gamma}\rfloor\equiv\frac{r_{i,l}-r^{c}_{i,l}}{\Gamma}\mod
M_{l}$, to estimate the folding number $\lfloor\frac{X_{i}}{\Gamma}\rfloor$,
we obtain

$\langle\lfloor\frac{\widetilde{r}_{il}}{\Gamma}\rfloor\rangle_{M_{l}}=\left\\{\begin{array}[]{lr}1).\langle\lfloor\frac{X_{i}}{\Gamma}\rfloor\rangle_{M_{l}},\leavevmode\nobreak\
when\leavevmode\nobreak\ r^{c}_{i}+\Delta_{il}\in[0,\Gamma)\\\
2).\langle\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\rangle_{M_{l}},\leavevmode\nobreak\
when\leavevmode\nobreak\ r^{c}_{i}+\Delta_{il}\in(-\Gamma,0)\\\
3).\langle\lfloor\frac{X_{i}}{\Gamma}\rfloor+1\rangle_{M_{l}},\leavevmode\nobreak\
when\leavevmode\nobreak\
r^{c}_{i}+\Delta_{il}\in[\Gamma,2\Gamma)\end{array}\right.$ (22)

Since $2\delta=\frac{\Gamma}{2}$, the three mentioned cases cannot happen
simultaneously. But case 1) and 2), or case 1) and 3), can happen at the same
time, which means there could exist $l_{1},l_{2}\in\\{1,2,...,L\\}$ such that
$r^{c}_{i}+\Delta_{il_{1}}\in[0,\Gamma)$ and
$r^{c}_{i}+\Delta_{il_{2}}\in(-\Gamma,0)$. Therefore,
$\langle\frac{\widetilde{r}^{c}_{il_{1}}}{\Gamma}\rangle_{M_{l_{1}}}$ and
$\langle\frac{\widetilde{r}^{c}_{il_{2}}}{\Gamma}\rangle_{M_{l_{2}}}$ can be
residues of either one of
$\\{\lfloor\frac{X_{i}}{\Gamma}\rfloor,\lfloor\frac{X_{i}}{\Gamma}\rfloor+1,\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\\}$.
Simply aggregating them via CRT will bring unpredictable reconstruction
errors. Without loss of generality, we only consider that case 1) and 2) occur
in the following.

To this end, following [31], we consider figuring out the order of $l$ such
that $\Delta_{il}$ are sorted in an ascending order. Said another way, it is
equivalent to determining a binary parameter $\tau_{il}\in\\{0,1\\}$ such that
the order of $l$, where
$\hat{r}^{c}_{il}=\widetilde{r}^{c}_{il}-\tau_{il}\Gamma$ are in ascending
order, is the same as that where $\Delta_{il}$ are sorted non-decreasingly. If
so,
$\langle\lfloor\frac{\widetilde{r}_{il}+\tau_{il}\Gamma}{\Gamma}\rfloor\rangle_{M_{l}}=\langle\frac{\widetilde{r}_{il}-\widetilde{r}^{c}_{il}}{\Gamma}+\tau_{il}\rangle_{M_{l}}$
are the residues of one integer, one of
$\\{\lfloor\frac{X_{i}}{\Gamma}\rfloor,\lfloor\frac{X_{i}}{\Gamma}\rfloor+1,\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\\}$.
Clearly, with proper $\tau_{il}$, the problem can be addressed by the
generalized CRT [26]. Furthermore,
$|(\widetilde{r}^{c}_{il_{1}}-\tau_{il_{1}}\Gamma)-(\widetilde{r}^{c}_{il_{2}}-\tau_{il_{2}}\Gamma)|$
becomes the minimum distance between $\widetilde{r}^{c}_{il_{1}}$ and
$\widetilde{r}^{c}_{il_{2}}$ in the ring of length $\Gamma$. Thus, for any
$l_{1},l_{2}$,

$|\widetilde{r}^{c}_{il_{1}}-\tau_{il_{1}}\Gamma-(\widetilde{r}^{c}_{il_{2}}-\tau_{il_{2}}\Gamma)|<2\delta=\frac{\Gamma}{2}$
(23)

Equation (23), called the first criterion in the following, is one of the keys
in the rest of the proof, which is also overlooked in previous works.

Now, we prove the following fact: Given a selection of $\\{\tau_{il}\\}$, if
they satisfy the first criterion, the shifted residues
$\\{\langle\lfloor\frac{\widetilde{r}_{il}+\tau_{il}\Gamma}{\Gamma}\rfloor\rangle_{M_{l}},l=1,2,...,L\\}$
are all residues of one of
$\\{\lfloor\frac{X_{i}}{\Gamma}\rfloor,\lfloor\frac{X_{i}}{\Gamma}\rfloor+1,\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\\}$.

For an arbitrary selection of $\tau_{il}$ assigned, let
$\mathbb{U}_{1}\in\\{1,2,...,L\\}$, where case 1) occurs for
$l_{1}\in\mathbb{U}_{1}$. Correspondingly,
$\mathbb{U}_{2}=\overline{\mathbb{U}}_{1}$, where case 2) happens for
$l_{2}\in\mathbb{U}_{2}$.

Analogously, we also specify the indices of $M_{l}$, where $\tau_{il}$ are set
to $0$ or $1$, respectively. We use $\mathbb{U}_{11}\subset\mathbb{U}_{1}$ to
denote those $l_{11}\in\mathbb{U}_{11}$ such that $\tau_{il_{11}}=0$.
Similarly, $l_{12}\in\mathbb{U}_{12}\subset\mathbb{U}_{1},$ where
$\tau_{il_{12}}=1$; $l_{21}\in\mathbb{U}_{21}\subset\mathbb{U}_{2}$, where
$\tau_{il_{22}}=0$, and $l_{22}\in\mathbb{U}_{22}\subset\mathbb{U}_{2}$,
$\tau_{il_{22}}=1$. Then, we have the following observations for the four
subcases, respectively: (1).
$\langle\frac{\widetilde{r}_{il_{11}}-\widetilde{r}^{c}_{il_{11}}}{\Gamma}+\tau_{il_{11}}\rangle_{M_{l}}=\langle\frac{X_{i}}{\Gamma}\rangle_{M_{l_{11}}}$;
(2).
$\langle\frac{\widetilde{r}_{il_{12}}-\widetilde{r}^{c}_{il_{12}}}{\Gamma}+\tau_{il_{12}}\rangle_{M_{l}}=\langle\frac{X_{i}}{\Gamma}+1\rangle_{M_{l_{12}}}$;
(3).
$\langle\frac{\widetilde{r}_{il_{21}}-\widetilde{r}^{c}_{il_{21}}}{\Gamma}+\tau_{il_{21}}\rangle_{M_{l}}=\langle\frac{X_{i}}{\Gamma}-1\rangle_{M_{l_{21}}}$;
(4).
$\langle\frac{\widetilde{r}_{il_{22}}-\widetilde{r}^{c}_{il_{22}}}{\Gamma}+\tau_{il_{22}}\rangle_{M_{l}}=\langle\frac{X_{i}}{\Gamma}\rangle_{M_{l_{22}}}$.

For those
$\langle\frac{\widetilde{r}_{il}-\widetilde{r}^{c}_{il}}{\Gamma}+\tau_{il}\rangle_{M_{l}}$,
we assume that there exists a solution of $N$ integers
$\\{q_{1},q_{2},...,q_{N}\\}$ within the dynamic range, i.e.,
$q_{i}\equiv\frac{\widetilde{r}_{k_{l}(i)l}-\widetilde{r}^{c}_{k_{l}(i)l}}{\Gamma}+\tau_{k_{l}(i)l}\mod
M_{l}$. Here, $k_{l}(i)$ is a permutation on $\\{1,2,...,N\\}$, indicating the
correspondence between $\widetilde{r}_{k_{l}(i)l}$ and $q_{i}$ for each
$M_{l}$. We elaborate on the wrong residue classifications below.

  * •

There exists $i\in\\{1,2,...,N\\}$ such that $k_{l_{11}}(i)=i_{0}$ and
$k_{l_{12}}(i)=i_{0}$. Then,
$|\widetilde{r}^{c}_{i_{0}l_{11}}-\tau_{i_{0}l_{11}}\Gamma-(\widetilde{r}^{c}_{i_{0}l_{12}}-\tau_{i_{0}l_{12}}\Gamma)|=|\widetilde{r}^{c}_{i_{0}l_{11}}-\widetilde{r}^{c}_{i_{0}l_{12}}+\Gamma|=|r^{c}_{i_{0}}+\Delta_{i_{0}l_{11}}-r^{c}_{i_{0}}-\Delta_{i_{0}l_{12}}+\Gamma|>\Gamma-2\delta=\frac{\Gamma}{2}$

  * •

There exists $i\in\\{1,2,...,N\\}$ such that $k_{l_{11}}(i)=i_{0}$ and
$k_{l_{21}}(i)=i_{0}$. Then
$|\widetilde{r}^{c}_{i_{0}l_{11}}-\tau_{il_{11}}\Gamma-(\widetilde{r}^{c}_{i_{0}l_{21}}-\tau_{il_{21}}\Gamma)|=|\widetilde{r}^{c}_{i_{0}l_{11}}-\widetilde{r}^{c}_{i_{0}l_{21}}|=|r^{c}_{i_{0}}+\Delta_{i_{0}l_{11}}-(r^{c}_{i_{0}}+\Delta_{i_{0}l_{21}}+\Gamma)|>\frac{\Gamma}{2}$

  * •

There exists $i\in\\{1,2,...,N\\}$ such that $k_{l_{21}}(i)=i_{0}$ and
$k_{l_{22}}(i)=i_{0}$. Then
$|\widetilde{r}^{c}_{i_{0}l_{21}}-\tau_{il_{21}}\Gamma-(\widetilde{r}^{c}_{i_{0}l_{22}}-\tau_{il_{22}}\Gamma)|=|r^{c}_{i_{0}}+\Delta_{i_{0}l_{21}}+\Gamma-(r^{c}_{i_{0}}+\Delta_{i_{0}l_{22}}+\Gamma-\Gamma)|>\frac{\Gamma}{2}$

  * •

There exists $i\in\\{1,2,...,N\\}$ such that $k_{l_{12}}(i)=i_{0}$ and
$k_{l_{22}}(i)=i_{0}$. Then
$|\widetilde{r}^{c}_{i_{0}l_{12}}-\tau_{il_{12}}\Gamma-(\widetilde{r}^{c}_{i_{0}l_{22}}-\tau_{il_{22}}\Gamma)|=|r^{c}_{i_{0}}+\Delta_{i_{0}l_{12}}-\Gamma-(r^{c}_{i_{0}}+\Delta_{i_{0}l_{22}}+\Gamma-\Gamma)|>\frac{\Gamma}{2}$

  * •

There exists $i\in\\{1,2,...,N\\}$ such that $k_{l_{12}}(i)=i_{0}$ and
$k_{l_{21}}(i)=i_{0}$. Then
$|\widetilde{r}^{c}_{i_{0}l_{12}}-\tau_{il_{12}}\Gamma-(\widetilde{r}^{c}_{i_{0}l_{21}}-\tau_{il_{21}}\Gamma)|=|r^{c}_{i_{0}}+\Delta_{i_{0}l_{12}}-\Gamma-(r^{c}_{i_{0}}+\Delta_{i_{0}l_{21}}+\Gamma)|>\frac{\Gamma}{2}$

Therefore, any incorrect selection of $\tau_{il}$ in $S_{i}$ will result in a
contradiction to the first criterion.

Next, we state the second criterion that the recovered
$\\{q_{1},q_{2},\dots,q_{N}\\}$ should be within the dynamic range
$[0,D_{q})$, where $D_{q}=\prod_{l=1}^{\lceil\frac{L}{N}\rceil}M_{l}.$ In the
following, we prove that a proper selection of $\tau_{il}$ and clustering
satisfying both criteria are sufficient for a robust reconstruction.

For $X_{i}<(\prod^{\lceil\frac{L}{N}\rceil}_{l=1}M_{l}-1)\Gamma$, we obtain
$\lfloor\frac{X_{i}}{\Gamma}\rfloor+1<\prod^{\lceil\frac{L}{N}\rceil}_{l=1}M_{l}$.
If the residues of $q_{i}$ are from one integer, we have
$q_{i}\in\\{\lfloor\frac{X_{i}}{\Gamma}\rfloor,\lfloor\frac{X_{i}}{\Gamma}\rfloor+1,\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\\}$,
i.e., $q_{i}<\prod^{\lceil\frac{L}{N}\rceil}_{l=1}M_{l}\Gamma$. We assume that
there exists another solution
$\\{q_{1}^{\prime},q_{2}^{\prime},\dots,q_{N}^{\prime}\\}$. Then, for the
residue classification of $q^{\prime}_{1}$, there must be at least
$\lceil\frac{L}{N}\rceil$ residues of $q^{\prime}_{1}$ are from one of
$\\{\lfloor\frac{X_{i}}{\Gamma}\rfloor,\lfloor\frac{X_{i}}{\Gamma}\rfloor+1,\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\\}$.
On the other hand, $q^{\prime}_{1}$ cannot share the residues of
$\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor$,
$\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor+1$ and
$\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor-1$ simultaneously. Otherwise, the
first criterion is violated. Thus, without loss of generality, we assume there
exist $\lceil\frac{L}{N}\rceil$ residues of $q^{\prime}_{1}$ are from
$\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor$, i.e.,
$q^{\prime}_{1}\equiv\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor\mod
lcm(M_{l}\in\mathbb{U}^{\prime})$, where
$|\mathbb{U}^{\prime}|=\lceil\frac{L}{N}\rceil$. Thus,
$q^{\prime}_{1}=klcm(M_{l}\in\mathbb{U}^{\prime})+\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor$,
where $k\geq 0$. With our assumption of $D_{q}$, if $k=0$,
$q^{\prime}_{1}=\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor$. Otherwise, $k\geq 1$,
$q^{\prime}_{1}>D_{q}$, a contradiction. Now, after determining the folding
numbers $q_{i}$, we obtain

$X_{i}=q_{i}\Gamma+\frac{\sum^{l=L}_{l=1}\hat{r}_{k_{l}(i)l}}{L}.$ (24)

Finally, we determine the worst-case reconstruction error. From the above
discussions, with a correct selection of $\tau_{il}$ and clustering such that
$S_{i}$ are all residues of $X_{i}$, which clearly satisfy the two criteria,
and then the reconstruction error is upper bounded by the maximal magnitude of
$\Delta_{il}$. However, there is a special case that the folding number of
different $X_{i}$ and $X_{j}$ may share the same residues: there exist some
$l\in\\{1,2,...,L\\}$ such that $i\not=j\in\\{1,2,...,N\\}$, $q_{i}\equiv
q_{j}\mod M_{l}$. They are essentially interchangeable, which will not lead to
the failure to recover $q_{i}$, but we have to quantify the reconstruction
error carefully.Since $\hat{r}^{c}_{jl}$ is assigned as a common residue to
$X_{i}$, i.e., the common residues assigned to $X_{i}$ are a combination of
$\hat{r}^{c}_{il_{1}}$ and $\hat{r}^{c}_{jl_{2}}$. Without loss of generality,
we assume that $q_{i_{0}}=\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor-1$, which
means there are at least $\lceil\frac{L}{N}\rceil$ common residues are from
$X_{i_{0}}$. It must fall into case 2) and $\tau_{i_{0}l}=0$. Since
$r^{c}_{i_{0}}+\Delta_{i_{0}l}<0$,
$\hat{r}^{c}_{i_{0}}=\langle\widetilde{r}^{c}_{i_{0}l}\rangle_{\Gamma}=r^{c}_{i_{0}}+\Delta_{i_{0}l}+\Gamma$.
The rest common residues satisfy
$\max|\hat{r}^{c}_{il}-\hat{r}^{c}_{i_{0}l}|<\frac{\Gamma}{2}$. Then,
$|\hat{X}_{i_{0}}-X_{i_{0}}|$ is equal to

$\displaystyle|(\lfloor\frac{X_{i_{0}}}{\Gamma}\rfloor-1)\Gamma+\frac{\sum^{l=\lceil\frac{L}{N}\rceil}_{l=1}\hat{r}^{c}_{i_{0}l}+\sum^{l=L}_{l=\lceil\frac{L}{N}\rceil+1}\hat{r}^{c}_{il}}{L}-X_{i_{0}}|$
(25)
$\displaystyle<|\frac{\sum^{l=\lceil\frac{L}{N}\rceil}_{l=1}\max\Delta_{i_{0}l}+\sum^{l=L}_{l=\lceil\frac{L}{N}\rceil+1}(\max\Delta_{i_{0}l}+\frac{\Gamma}{2})}{L}|$
$\displaystyle<|\frac{\lceil\frac{L}{N}\rceil\frac{\Gamma}{4}+(L-\lceil\frac{L}{N}\rceil)\frac{3\Gamma}{4}}{L}|<\frac{3\Gamma}{4}.$

Therefore, the reconstruction is error-bounded by $\frac{3\Gamma}{4}$. ∎

##  Appendix B Proof of Corollary 2

###### Proof.

Likewise, let $q_{i}$ denote the folding number estimations, where
$q_{i}\in\\{\lfloor\frac{X_{i}}{\Gamma}\rfloor,\lfloor\frac{X_{i}}{\Gamma}\rfloor+1,\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\\}$.
It is noted that $\min_{d\in\\{0,\pm
1\\}}|r_{i_{0}l}-r_{i_{1}l}+dM_{l}\Gamma|$ represents the minimum distance
between $r_{i_{0}l}$ and $r_{i_{1}l}$ on the circle of length $M_{l}\Gamma$.
Then, we have the following facts and at least one of the three cases happens:
1). $r_{il}=\langle q_{i}\Gamma\rangle_{M_{l}\Gamma}+r^{c}_{i}$, when
$q_{i}=\lfloor\frac{X_{i}}{\Gamma}\rfloor$; 2). $r_{il}=\langle
q_{i}\Gamma\rangle_{M_{l}\Gamma}+r^{c}_{i}-\Gamma$, when
$q_{i}=\lfloor\frac{X_{i}}{\Gamma}\rfloor+1$; 3). $r_{il}=\langle
q_{i}\Gamma\rangle_{M_{l}\Gamma}+r^{c}_{i}+\Gamma$, when
$q_{i}=\lfloor\frac{X_{i}}{\Gamma}\rfloor-1$.

When $r_{i_{0}l}$ and $r_{i_{1}l}$ satisfy 2) and 3) respectively, since
$\min|r_{i_{0}l}-r_{i_{1}l}+dM_{l}\Gamma|>3\Gamma$, replacing $r_{i_{0}l}$ and
$r_{i_{1}l}$ with the right hand of 2) and 3), we have

$\min|\langle q_{i_{0}}\Gamma\rangle_{M_{l}\Gamma}-\langle
q_{i_{1}}\Gamma\rangle_{M_{l}\Gamma}+dM_{l}\Gamma+(r^{c}_{i_{0}}-\Gamma-r^{c}_{i_{1}}-\Gamma)|>3\Gamma,$
(26)

and thus $|r^{c}_{i_{0}}-r^{c}_{i_{1}}-2\Gamma|<3\Gamma$, which results in
$q_{i_{0}}\not\equiv q_{i_{1}}\mod M_{l}$. The same conclusion can be derived
when $r_{i_{0}l}$ and $r_{i_{1}l}$ satisfy any two equations of 1), 2) and 3),
i.e., for each $l$, $q_{i_{0}}\not\equiv q_{i_{1}}\mod M_{l}$ holds. If
$\hat{r}^{c}_{i_{0}l}$ and $\hat{r}^{c}_{i_{1}l}$ are clustered into one set
$S_{i}$, the estimated folding number $q_{i}$ satisfies

$\left\\{\begin{array}[]{lr}q_{i}\equiv q_{i_{0}}\mod
lcm(M_{l}\in\mathbb{U}_{i_{0}})\\\ q_{i}\equiv q_{i_{1}}\mod
lcm(M_{l}\in\mathbb{U}_{i_{1}}).\end{array}\right.$ (27)

Based on the pigeonhole principle, at least $\lceil\frac{L}{N}\rceil$ residues
in $S_{i}$ are from one integer. Without loss of generality, we assume that
$|\mathbb{U}_{i_{0}}|\geq\lceil\frac{L}{N}\rceil$. Therefore,
$q_{i}=klcm(M_{l}\in\mathbb{U}_{i_{0}})+q_{i_{0}}$. Because
$q_{i_{0}}\not\equiv q_{i_{1}}\mod M_{l}$, we have $q_{i}\not=q_{i_{0}}$ and
$k>0$, which leads to a contradiction to the second criterion. That is to say,
if for each $l$, $q_{i_{0}}\not\equiv q_{i_{1}}\mod M_{l}$, clustering
residues from different integers violates the second criterion. So, for each
$q_{i}$, the residues must be all from $X_{i}$, which provides a sharpened
robust reconstruction error bounded by $\frac{\Gamma}{4}$. 666If $r_{il}$ is
assumed to be uniformly distributed across $[0,M_{l}\Gamma)$, we have
$\Pr(\min_{d}|r_{i_{0}l}-r_{i_{1}l}+dM_{l}\Gamma|>3\Gamma)=\prod_{l=1}^{L}\frac{M_{l}-6}{M_{l}}$.
∎

##  Appendix C Closed-form Robust Reconstruction for Single Tone Real
Waveforms

To be self-contained, we first introduce the setup and notations. Assume that
the moduli are in a form $\\{m_{l}=\Gamma M_{l}|l=1,2,\dots,L\\}$, where
$\\{M_{l}\\}$ are relatively co-prime, ranging in an ascending order. The
dynamic range of $X$ is represented by $D$, i.e., $X\in[0,D)$, where
$D\leq\min_{\mathbb{U}\subset\\{M_{1},M_{2},\dots,M_{L}\\}}(\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1)\Gamma$.
The residue set is $\mathcal{R}_{l}=\\{\widetilde{r}^{+}_{l}=\langle
X+\Delta^{+}_{l}\rangle_{m_{l}},\widetilde{r}^{-}_{l}=\langle-X+\Delta^{-}_{l}\rangle_{m_{l}}\\}$,
where $|\Delta^{\pm}_{l}|<\delta=\frac{\Gamma}{4}$.

With a similar idea, we consider selecting a residue from each
$\mathcal{R}_{l}$ randomly to obtain a $L$-residue clustering $S$, which
possibly contains residues from $X$ and $-X$ simultaneously. To recover the
folding number $\lfloor\frac{X}{\Gamma}\rfloor$ robustly, we apply (28) and
(29) on residues respectively. At least one of the operations leads to
residues sorted in ascending order, i.e., residues are sorted in the same
order as that of $\Delta^{\pm}_{l}$. Clearly, there exists an efficient
solution when all residues in $S$ are from $X$ and are properly sorted.
Otherwise, as proved in the next section, one can distinguish that, where at
least one of the two criteria, (30) and (31), will not hold.

Algorithm 2 Robust Remaindering Decoding of Single Tone Real Waveform

Input: Moduli set: $\mathcal{M}=\\{m_{l}=M_{l}\Gamma|l=1,2,\dots,L\\}$;
Residue Sets:
$\mathcal{R}_{l}=\\{\widetilde{r}^{+}_{l},\widetilde{r}^{-}_{l}\\}$,
$l=1,2,\dots,L$.

1: Repeat: Propose a clustering assignment

2: Following the proposed clustering by selecting one residue from each
$\mathcal{R}_{l}$ to obtain a $L$-residue clustering
$S=\\{\widetilde{r}_{(l)}|l=1,2,\dots,L\\}$, where $(l)$ denotes the index of
the residues assigned to $S$.

3: For each $\widetilde{r}_{(l)}$ in $S$, calculate the corresponding shifted
common residues $\hat{r}_{c,(l)}$ based on (28) and (29) respectively. We
obtain the shifted common residue set $S_{c1}$ and $S_{c2}$.

4: For $S_{c1}$ and $S_{c2}$, calculate
$q\equiv\frac{\widetilde{r}_{(l)}-\hat{r}_{c,(l)}}{\Gamma}\mod M_{l}$ via CRT
respectively.

5: Until: For $S_{c1}$ or $S_{c2}$, if each shifted common residue satisfies
(30) and the corresponding $q$ satisfies that $q\in[0,D_{q})$, where
$D_{q}=\min_{\mathbb{U}\subset\\{M_{1},M_{2},\dots,M_{L}\\}}\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1$,
we output the one that passes the two tests.

Output: $\hat{X}=q\Gamma+\frac{\sum^{L}_{l=1}\hat{r}_{c,(l)}}{L}$.

##  Appendix D Proof of Theorem 5

###### Proof.

With the fact that, $r^{-}_{c}=\Gamma-r^{+}_{c}$, it is clear that
$\\{r^{+}_{c},r^{-}_{c}\\}$ must fall into one of the following four cases,

  1. 1.

$r^{+}_{c}\in[{\Gamma}/{4},{\Gamma}/{2})$ and
$r^{-}_{c}\in({\Gamma}/{2},{3\Gamma}/{4}]$

  2. 2.

$r^{+}_{c}\in[{\Gamma}/{2},{3\Gamma}/{4})$ and
$r^{-}_{c}\in({\Gamma}/{4},{\Gamma}/{2}]$

  3. 3.

$r^{+}_{c}\in[0,{\Gamma}/{4})$ and $r^{-}_{c}\in({3\Gamma}/{4},\Gamma]$

  4. 4.

$r^{+}_{c}\in[{3\Gamma}/{4},\Gamma)$ and $r^{-}_{c}\in(0,{\Gamma}/{4}]$

Recall (22), to ensure robustness, we consider applying
$\tau^{+}_{l}\in\\{0,1\\}$ to determine the order of $l$ such that
$\Delta^{+}_{l}$ are in ascending order, i.e.,
$\hat{r}^{+}_{c,l}=\widetilde{r}^{+}_{c,l}-\tau^{+}_{l}\Gamma$ are sorted in
the order of $\Delta^{+}_{l}$. Here, we use
$\widetilde{r}^{+}_{c,l}=\langle\widetilde{r}^{+}_{l}\rangle_{\Gamma}$, the
noisy observation of common residue from the $l$-th sampler. Similarly,
$\hat{r}^{-}_{c,l}=\widetilde{r}^{-}_{c,l}-\tau^{-}_{l}\Gamma$ can be sorted
based on the same idea. To this end, we consider the following operations on
the residues if we can distinguish cases (1,2) and (3,4).

  * •

operation 1:

$\left\\{\begin{array}[]{lr}\hat{r}^{+}_{c,l}=\widetilde{r}^{+}_{c,l}\\\
\hat{r}^{-}_{c,l}=\widetilde{r}^{-}_{c,l}\end{array}\right.$ (28)

  * •

operation 2:

$\left\\{\begin{array}[]{lr}\hat{r}^{+}_{c,l}=\widetilde{r}^{+}_{c,l}\leavevmode\nobreak\
\text{when}\leavevmode\nobreak\
\widetilde{r}^{+}_{c,l}\in[0,\frac{\Gamma}{2}),\leavevmode\nobreak\
\text{or}\leavevmode\nobreak\
\hat{r}^{+}_{c,l}=\widetilde{r}^{+}_{c,l}-\Gamma\\\
\hat{r}^{-}_{c,l}=\widetilde{r}^{-}_{c,l}\leavevmode\nobreak\
\text{when}\leavevmode\nobreak\
\widetilde{r}^{-}_{c,l}\in[0,\frac{\Gamma}{2}),\leavevmode\nobreak\
\text{or}\leavevmode\nobreak\
\hat{r}^{-}_{c,l}=\widetilde{r}^{-}_{c,l}-\Gamma\end{array}\right.$ (29)

To arrange $\widetilde{r}^{\pm}_{c,l}$ in an ascending order in cases (1,2),
since $\delta=\frac{\Gamma}{4}$, we will apply operation 1, defined in (28),
on $\widetilde{r}^{\pm}_{c,l}$, where
$\hat{r}^{\pm}_{c,l}=\widetilde{r}^{\pm}_{c,l}={r}^{\pm}_{c}+\Delta^{\pm}_{l}$.
Consequently,
$(\langle\frac{\widetilde{r}^{+}_{l}-\hat{r}^{+}_{c,l}}{\Gamma}\rangle_{M_{l}},\langle\frac{\widetilde{r}^{-}_{l}-\hat{r}^{-}_{c,l}}{\Gamma}\rangle_{M_{l}})$
are residues of $(Y,-Y-1)$ modulo $M_{l}$, where $Y$ denotes the folding
number $\lfloor\frac{X}{\Gamma}\rfloor$.

Similarly, to sort $\widetilde{r}^{\pm}_{c,l}$ non-decreasingly in case 3),
since $r^{+}_{c,l}+\Delta^{+}_{l}\in(-\frac{\Gamma}{4},\frac{\Gamma}{2})$ and
$r^{-}_{c,l}+\Delta^{-}_{l}\in(\frac{\Gamma}{2},\frac{5\Gamma}{4})$, by
applying operation 2, shown in (29),
$(\langle\frac{\widetilde{r}^{+}_{l}-\hat{r}^{+}_{c,l}}{\Gamma}\rangle_{M_{l}},\langle\frac{\widetilde{r}^{-}_{l}-\hat{r}^{-}_{c,l}}{\Gamma}\rangle_{M_{l}})$
are residues of $(Y,-Y)$ modulo $M_{l}$.

In case 4), which is the dual circumstance of case 3),
$(\langle\frac{\widetilde{r}^{+}_{l}-\hat{r}^{+}_{c,l}}{\Gamma}\rangle_{M_{l}},\langle\frac{\widetilde{r}^{-}_{l}-\hat{r}^{-}_{c,l}}{\Gamma}\rangle_{M_{l}})$
are residues of $(Y+1,-Y-1)$ modulo $M_{l}$.

Provided sorted $\hat{r}^{\pm}_{c,l}$, since $\delta=\frac{\Gamma}{4}$, we
have

$|\hat{r}_{c,l_{1}}-\hat{r}_{c,l_{2}}|<2\delta=\frac{\Gamma}{2},$ (30)

where $\hat{r}_{c,l_{1}}$ and $\hat{r}_{c,l_{2}}$ denote the shifted common
residues derive from one integer. Similar to the complex waveform case, (30)
behaves as a criterion for a valid test. However, the main obstacle is the
misuse of operation 1 or 2, i.e., applying operation 1 on cases (3,4) or
operation 2 on cases (1,2). For example, if operation 1 is applied in case 3,
$\langle\frac{\widetilde{r}^{+}_{l}-\hat{r}^{+}_{c,l}}{\Gamma}\rangle_{M_{l}}$
can be either $\langle Y\rangle_{M_{l}}$ or $\langle Y-1\rangle_{M_{l}}$; and
$\langle\frac{\widetilde{r}^{-}_{l}-\hat{r}^{-}_{c,l}}{\Gamma}\rangle_{M_{l}}$
can be either $\langle-Y-1\rangle_{M_{l}}$ or $\langle-Y\rangle_{M_{l}}$,
which could produce an unpredictable reconstruction error. Indeed, it is
impossible to distinguish the four cases just from the locations of
$\widetilde{r}^{\pm}_{c,l}$. But anyway, applying the two operations on
$\widetilde{r}^{\pm}_{c,l}$, at least one of them can correctly recover $Y$.
The rest proof is developed by two parts. First, if one happens to apply the
proper operation, a robust estimation can be achieved. Second, if a wrong
operation is applied, one can distinguish that.

Now, we prove that with correct operation, $X\in[0,D)$ is a sufficient
condition that the folding number has a unique representation by residues
attained, i.e., $\\{Y,-Y-1\\}$, $\\{Y,-Y\\}$ and $\\{Y+1,-Y-1\\}$ can be
uniquely determined. Based on the conclusion from Section V-A, for
$\\{Y,-Y\\}$, the largest dynamic range of $Y$ is
$\min_{\mathbb{U}\subset\\{M_{1},M_{2},...,M_{L}\\}}\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}$.
Analogously, the residue representation of $\\{Y+1,-Y-1\\}$ is unique once
$Y<\min_{\mathbb{U}\subset\\{M_{1},M_{2},...,M_{L}\\}}\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1$.
Similarly, for $\\{Y,-Y-1\\}$, the largest dynamic range is
$\min_{\mathbb{U}\subset\\{M_{1},M_{2},...,M_{L}\\}}\
\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}-1}{2}$.
As a summary, the folding number $Y$ has a unique residue representation when

$Y<D_{q}=\min_{\mathbb{U}\subset\\{M_{1},M_{2},...,M_{L}\\}}\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1.$
(31)

Since
$X\in[0,\min_{\mathbb{U}\subset\\{M_{1},M_{2},\dots,M_{L}\\}}(\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1)\Gamma)$,
it’s clear that $Y<D_{q}$, satisfying (31).

In the following, we switch to the second part and take a closer look at the
scenario when a wrong operation is applied. To avoid tedious discussion, we
only elaborate on case 2) in the following proof. We list all the possible
situations below. Let $\mathbb{U}_{11}$ denote those
$l_{11}\in\mathbb{U}_{11}$ such that
$r^{+}_{c}+\Delta^{+}_{l_{11}}\geq\frac{\Gamma}{2}$. Similarly, when
$l_{12}\in\mathbb{U}_{12}$, $r^{+}_{c}+\Delta^{+}_{l_{12}}<\frac{\Gamma}{2}$;
for $l_{21}\in\mathbb{U}_{21}$,
$r^{-}_{c}+\Delta^{-}_{l_{21}}\geq\frac{\Gamma}{2}$; for
$l_{22}\in\mathbb{U}_{22}$, $r^{-}_{c}+\Delta^{-}_{l_{22}}<\frac{\Gamma}{2}$.
The wrong residue classifications violating (30) are elaborated below.

  * •

If $\hat{r}^{+}_{c,l_{11}}$ and $\hat{r}^{+}_{c,l_{12}}$ are clustered
together,
$|\hat{r}^{+}_{c,l_{11}}-\hat{r}^{+}_{c,l_{12}}|=|r^{+}_{c}+\Delta_{l_{11}}-\Gamma-(r^{+}_{c}+\Delta_{l_{12}})|=|\Gamma-(\Delta_{l_{11}}-\Delta_{l_{12}})|>\frac{\Gamma}{2}$.

  * •

If $\hat{r}^{-}_{c,l_{21}}$ and $\hat{r}^{-}_{c,l_{22}}$ are clustered
together,
$|\hat{r}^{-}_{c,l_{21}}-\hat{r}^{+}_{c,l_{22}}|=|r^{-}_{c}+\Delta_{l_{21}}-\Gamma-(r^{-}_{c}+\Delta_{l_{22}})|=|\Gamma-(\Delta_{l_{21}}-\Delta_{l_{22}})|>\frac{\Gamma}{2}$.

  * •

If $\hat{r}^{+}_{c,l_{12}}$ and $\hat{r}^{-}_{c,l_{21}}$ are clustered
together, since $\hat{r}^{+}_{c,l_{12}}\in(\frac{\Gamma}{4},\frac{\Gamma}{2})$
and $\hat{r}^{-}_{c,l_{21}}\in(-\frac{\Gamma}{2},-\frac{\Gamma}{4})$,
$|\hat{r}^{+}_{c,l_{12}}-\hat{r}^{-}_{c,l_{21}}|>\frac{\Gamma}{2}$.

Therefore, under the restriction of (30), the estimated folding number $q$
must be in one of the following forms in case 2) applied with operation 2.

  * 1

$q\equiv Y+1\mod lcm(M_{l}\in\mathbb{U}_{11})$ and $q\equiv-Y-1\mod
lcm(M_{l}\in\mathbb{U}_{22})$

  * 2

$q\equiv Y+1\mod lcm(M_{l}\in\mathbb{U}_{11})$ and $q\equiv-Y\mod
lcm(M_{l}\in\mathbb{U}_{21})$

  * 3

$q\equiv Y\mod lcm(M_{l}\in\mathbb{U}_{12})$ and $q\equiv-Y-1\mod
lcm(M_{l}\in\mathbb{U}_{22})$

With the above understanding, finally, we show $\hat{X}$ can be recovered
error-bounded by $\frac{3\Gamma}{4}$. Obviously, $q$ must be in the form under
operation 1:

  * •

$q\equiv Y\mod lcm(M_{l}\in\mathbb{U}_{1})$ and $q\equiv-Y-1\mod
lcm(M_{l}\in\mathbb{U}_{2})$.

The lower bound of
$\max\\{lcm(M_{l}\in\mathbb{U}_{1}),lcm(M_{l}\in\mathbb{U}_{2})\\}$ is
$\omega=\min_{\mathbb{U}\subset\\{M_{1},M_{2},\dots,M_{L}\\}}\max\\{\prod_{l\in\mathbb{U}}M_{l},\prod_{l\in\overline{\mathbb{U}}}M_{l}\\}$,
which is bigger than $D_{q}$. Without loss of generality, we assume that
$lcm(M_{l}\in\mathbb{U}_{1})\geq lcm(M_{l}\in\mathbb{U}_{2})$. Thus,

$q=k\cdot lcm(M_{l}\in\mathbb{U}_{1})+Y,k\geq 0$ (32)

If $Y\not\equiv-Y-1\mod lcm(M_{l}\in\mathbb{U}_{2})$ and
$\mathbb{U}_{2}\not=\emptyset$, $k\geq 1$ and $q$ clearly exceeds $D_{q}$, a
contradiction. Only when $Y\equiv-Y-1\mod lcm(M_{l}\in\mathbb{U}_{2})$ or
$\mathbb{U}_{2}=\emptyset$, can $q$ be within the dynamic range, i.e., $k=0$
and $q=Y$. Thus, the shifted common residues assigned to $q$ are
$\hat{r}^{+}_{c,l_{1}}$ and $\hat{r}^{-}_{c,l_{2}}$, where
$l_{1}\in\mathbb{U}_{1}$ and $l_{2}\in\mathbb{U}_{2}$. It is worth pointing
out that $\mathbb{U}_{2}=\emptyset$ may hold, i.e., the shifted common
residues assigned to $q$ are all in a form $\hat{r}^{+}_{c,l_{1}}$. For
$\hat{r}^{\pm}_{c,l}=\widetilde{r}^{\pm}_{c,l}$, we have
$|\widetilde{r}^{-}_{c,l_{2}}-\widetilde{r}^{+}_{c,l_{1}}|=|\widetilde{r}^{-}_{c,l_{2}}-(r^{+}_{c,l_{1}}+\Delta^{+}_{c,l_{1}})|<\frac{\Gamma}{2}$
based on (30). Since $|\Delta^{+}_{c,l_{1}}|<\frac{\Gamma}{4}$,
$|\widetilde{r}^{-}_{c,l_{2}}-r^{+}_{c}|<\frac{3\Gamma}{4}$ holds. Thus,

$\displaystyle|\hat{X}-X|$
$\displaystyle=|q\Gamma+\frac{\sum_{l_{1}\in\mathbb{U}_{1}}\widetilde{r}^{+}_{c,l_{1}}+\sum_{l_{2}\in\mathbb{U}_{2}}\widetilde{r}^{-}_{c,l_{2}}}{L}-X|$
(33)
$\displaystyle<|\lfloor\frac{X}{\Gamma}\rfloor\Gamma+r^{+}_{c}+\frac{3\Gamma}{4}-X|=\frac{3\Gamma}{4}.$

When we apply operation 2 to case 2), $q$ will be one of the three kinds.

  * 1

$q\equiv Y+1\mod lcm(M_{l}\in\mathbb{U}_{11})$ and $q\equiv-Y\mod
lcm(M_{l}\in\mathbb{U}_{21})$

  * 2

$q\equiv Y+1\mod lcm(M_{l}\in\mathbb{U}_{11})$ and $q\equiv-Y-1\mod
lcm(M_{l}\in\mathbb{U}_{22})$

  * 3

$q\equiv Y\mod lcm(M_{l}\in\mathbb{U}_{12})$ and $q\equiv-Y-1\mod
lcm(M_{l}\in\mathbb{U}_{22})$

We start with the first situation, where
$\Gamma>r^{+}_{c}+\Delta^{+}_{l_{11}}\geq\frac{\Gamma}{2}$ and
$\Gamma>r^{-}_{c}+\Delta^{-}_{l_{21}}\geq\frac{\Gamma}{2}$. Likewise, under
the assumption $lcm(M_{l}\in\mathbb{U}_{11})\geq
lcm(M_{l}\in\mathbb{U}_{21})$, $q$ would be in a form
$q=klcm(M_{l}\in\mathbb{U}_{11})+Y+1$. When $Y+1\not\equiv-Y\mod
lcm(M_{l}\in\mathbb{U}_{21})$ and $\mathbb{U}_{21}\not=\emptyset$, $q>\omega$
with $k>0$, i.e., $q=Y+1$ is the only solution. If $Y+1\equiv-Y\mod
lcm(M_{l}\in\mathbb{U}_{21})$ or $\mathbb{U}_{21}=\emptyset$, the shifted
common residues assigned to $q$ are $\hat{r}^{+}_{c,l_{11}}$ and
$\hat{r}^{-}_{c,l_{21}}$. For

$\begin{aligned}
\hat{X}&=(Y+1)\Gamma+\frac{\sum_{l_{11}\in\mathbb{U}_{11}}\hat{r}^{+}_{c,l_{11}}+\sum_{l_{21}\in\mathbb{U}_{21}}\hat{r}^{-}_{c,l_{21}}}{L}\\\
&=Y\Gamma+\frac{\sum_{l_{11}\in\mathbb{U}_{11}}\widetilde{r}^{+}_{c,l_{11}}+\sum_{l_{21}\in\mathbb{U}_{21}}\widetilde{r}^{-}_{c,l_{21}}}{L}\end{aligned},$

we have

$|\hat{X}-X|=|Y\Gamma+\frac{\sum_{l_{11}\in\mathbb{U}_{11}}\widetilde{r}^{+}_{c,l_{11}}+\sum_{l_{21}\in\mathbb{U}_{21}}\widetilde{r}^{-}_{c,l_{21}}}{L}-X|<\frac{3\Gamma}{4}.$

The proof is similar for the rest situations, which is omitted for simplicity.
∎

######  Corollary 3.

If $\min_{d=0,\pm 1}|r^{+}_{l}-r^{-}_{l}+dM_{l}\Gamma|>\Gamma$ holds for each
$l$, $X$ can be recovered robustly error bounded by $\frac{\Gamma}{4}$ for any
$X\in[0,D)$.

###### Proof.

Similarly, we only elaborate on case 2) for brevity. $\min_{d=0,\pm
1}|r^{+}_{l}-r^{-}_{l}+dM_{l}\Gamma|>\Gamma$ denotes the closest distance
between $r^{+}_{l}$ and $r^{-}_{l}$ on the circle of length $M_{l}\Gamma$.

$r^{+}_{l}=\langle X\rangle_{M_{l}\Gamma}=\langle
Y\Gamma\rangle_{M_{l}\Gamma}+r^{+}_{c}$ (34)
$r^{-}_{l}=\langle-X\rangle_{M_{l}\Gamma}=\langle(-Y-1)\Gamma\rangle_{M_{l}\Gamma}+r^{-}_{c}$
(35)

For $\min|r^{+}_{l}-r^{-}_{l}+dM_{l}\Gamma|>\Gamma$, replacing $r^{+}_{l}$ and
$r^{-}_{l}$ with the right hand of (34) and (35) respectively, we obtain

$\min|\langle
Y\Gamma\rangle_{M_{l}\Gamma}-\langle(-Y-1)\Gamma\rangle_{M_{l}\Gamma}+(r^{+}_{c}-r^{-}_{c})+dM_{l}\Gamma|>\Gamma.$

Since $0\leq r^{+}_{c}-r^{-}_{c}\leq\frac{\Gamma}{2}$,

$\min|\langle
Y\Gamma\rangle_{M_{l}\Gamma}-\langle(-Y-1)\Gamma\rangle_{M_{l}\Gamma}+dM_{l}\Gamma|\neq
0,$

i.e., for each $l$, $Y\not\equiv-Y-1\mod M_{l}$. Analogously,
$Y+1\not\equiv-Y\mod M_{l}$ and $Y+1\not\equiv-Y-1\mod M_{l}$. When $q$ is in
the form (32), since $Y\not\equiv-Y-1\mod M_{l}$, only if
$\mathbb{U}_{2}=\emptyset$, will $q$ be within the range, i.e., the shifted
common residues assigned to $q$ are in a form $\hat{r}^{+}_{c,l}$. Thus,

$|\hat{X}-X|=|Y\Gamma+\frac{\sum^{L}_{l=1}r^{+}_{c,l}}{L}-X|<\frac{\Gamma}{4}.$

One can similarly verify the claim in other situations, which is omitted for
simplicity. ∎

Hence, it is distinguishable if we classify residues of both $X$ and $-X$ into
one clustering when $\min_{d=0,\pm
1}|r^{+}_{l}-r^{-}_{l}+dM_{l}\Gamma|>\Gamma$ holds for each $l$. If
$r^{\pm}_{l}$ is uniformly distributed over the ring with a circumference of
$M_{l}\Gamma$,

$\Pr(\min_{d=0,\pm
1}|r^{+}_{l}-r^{-}_{l}+dM_{l}\Gamma|>\Gamma)=\prod^{L}_{l=1}\frac{M_{l}-2}{M_{l}},$

which can be close to $1$ given sufficiently large $M_{l}$.

##  Appendix E Proof of Theorem 6

###### Proof.

With a similar reasoning, the proof is organized into two parts. First, we
prove that when $X_{i}\in[0,D)$, their folding numbers have a unique residue
representation. Second, with proper residue classifications and sorting, one
can uniquely derive a robust estimation for each $X_{i}$.

The proof for the first part is straightforward, where the residue
representation of the folding number $\lfloor\frac{X_{i}}{\Gamma}\rfloor$ is
unique. In the following, we denote $\lfloor\frac{X_{i}}{\Gamma}\rfloor$ by
$Y_{i}$. Clearly, since $X_{i}<D$, we have

$Y_{i}<\min_{\mathbb{U}\subset\\{M_{1},M_{2},\dots,M_{\lceil\frac{L}{N}\rceil}\\}}\frac{\prod_{l\in\mathbb{U}}M_{l}+\prod_{l\in\overline{\mathbb{U}}}M_{l}}{2}-1=D_{q},$

where
$\mathbb{U}\cup\overline{\mathbb{U}}=\\{M_{1},M_{2},\dots,M_{\lceil\frac{L}{N}\rceil}\\}$,
satisfying (31).

Randomly selecting a residue from each
$\mathcal{R}_{l}=\\{\widetilde{r}^{+}_{i,l},\widetilde{r}^{-}_{i,l}|i=1,2,\dots,N\\}$,
we obtain a $L$-residue cluster
$S_{i}=\\{\widetilde{r}_{(i,l)}|i=1,2,\dots,N\\}$, where $(i,l)$ denotes the
index of the residues assigned to $S_{i}$. $S_{i}$ may contain residues from
different pairs $\\{X_{i},-X_{i}\\}$. Anyway, based on the pigeonhole
principle, at least $\lceil\frac{L}{N}\rceil$ residues are from some
$\\{X_{i_{0}},-X_{i_{0}}\\}$ in $S_{i}$. The estimated folding number obtained
from $S_{i}$ is denoted by $q_{i}$.

To avoid tedious discussion, we only consider the case where $r^{+}_{ic}\leq
r^{-}_{ic}$. Then, considering the tuple of residues
$(\langle\lfloor\frac{\widetilde{r}^{+}_{i,l}}{\Gamma}\rfloor\rangle_{M_{l}},\langle\lfloor\frac{\widetilde{r}^{-}_{i,l}}{\Gamma}\rfloor\rangle_{M_{l}})$,
it must be the residue of one of the four integer-pairs modulo $M_{l}$:

  1. 1.

$(Y_{i},-Y_{i}-1)$, when
$r^{+}_{ic,l}+\Delta^{+}_{i,l}\in[0,\Gamma),\leavevmode\nobreak\
r^{-}_{ic,l}+\Delta^{-}_{i,l}\in[0,\Gamma)$;

  2. 2.

$(Y_{i},-Y_{i})$, when
$r^{+}_{ic,l}+\Delta^{+}_{i,l}\in[0,\Gamma),\leavevmode\nobreak\
r^{-}_{ic,l}+\Delta^{-}_{i,l}\in[\Gamma,2\Gamma)$;

  3. 3.

$(Y_{i}-1,-Y_{i}-1)$, when
$r^{+}_{ic,l}+\Delta^{+}_{i,l}\in(-\Gamma,0),\leavevmode\nobreak\
r^{-}_{ic,l}+\Delta^{-}_{i,l}\in[0,\Gamma)$;

  4. 4.

$(Y_{i}-1,-Y_{i})$, when
$r^{+}_{ic,l}+\Delta^{+}_{i,l}\in(-\Gamma,0),\leavevmode\nobreak\
r^{-}_{ic,l}+\Delta^{-}_{i,l}\in[\Gamma,2\Gamma)$.

Perturbed by errors,
$\langle\lfloor\frac{\widetilde{r}^{+}_{i,l}}{\Gamma}\rfloor\rangle_{M_{l}}$
can be the residue of either
$\langle\lfloor\frac{X_{i}}{\Gamma}\rfloor\rangle_{M_{l}}$ or
$\langle\lfloor\frac{X_{i}}{\Gamma}\rfloor-1\rangle_{M_{l}}$. To test that
whether $\widetilde{r}_{(i,l)}$ in $S_{i}$ are all residues of one integer,
similarly, we introduce a binary variable $\tau_{(i,l)}\in\\{0,1\\}$ to each
$\widetilde{r}_{(i,l)}$. Then, the shifted common residues are
$\hat{r}_{(ic,l)}=\langle\widetilde{r}_{(i,l)}\rangle_{\Gamma}-\tau_{(i,l)}\Gamma=\widetilde{r}_{(ic,l)}-\tau_{(i,l)}\Gamma$.
Clearly, if all $\widetilde{r}_{(i,l)}$ are residues of $X_{i}$ assigned with
proper $\tau_{(i,l)}$, $\hat{r}_{(ic,l)}$ would satisfy

$|\hat{r}_{(ic,l_{1})}-\hat{r}_{(ic,l_{2})}|<\frac{\Gamma}{2}.$ (36)

Then,
$\langle\lfloor\frac{\widetilde{r}_{(i,l)}}{\Gamma}\rfloor+\tau_{(i,l)}\rangle_{M_{l}}$
must be the residue of one of the six integers modulo $M_{l}$: (a).$Y_{i}$,
for $l\in\mathbb{U}_{a}$; (b).$Y_{i}-1$, for $l\in\mathbb{U}_{b}$;
(c).$Y_{i}+1$, for $l\in\mathbb{U}_{c}$; (d).$-Y_{i}$, for
$l\in\mathbb{U}_{d}$; (e).$-Y_{i}-1$, for $l\in\mathbb{U}_{e}$;
(f).$-Y_{i}+1$, for $l\in\mathbb{U}_{f}$.

With a similar reasoning as previous proof, when
$r^{+}_{ic}\in[\frac{\Gamma}{4},\frac{\Gamma}{2}]$, only case 1) may happen in
(36), which results in that only cases a), c), d) and e) can happen
accordingly. It is noted that $|\mathbb{U}_{a}|\cdot|\mathbb{U}_{c}|=0$, i.e.,
one of $\mathbb{U}_{a}$ and $\mathbb{U}_{c}$ should be empty, since for any
$l_{1}\in\mathbb{U}_{a}$ and $l_{2}\in\mathbb{U}_{c}$,

$|\widetilde{r}_{ic,l_{1}}-(\widetilde{r}_{ic,l_{2}}-\Gamma)|\geq|(r^{+}_{ic}+\Delta^{+}_{i,l_{1}})-(r^{+}_{ic}+\Delta^{+}_{i,l_{2}}-\Gamma)|>\frac{\Gamma}{2},$
(37)

violating (36). Likewise, we have $|\mathbb{U}_{d}|\cdot|\mathbb{U}_{e}|=0$
and $|\mathbb{U}_{c}|\cdot|\mathbb{U}_{e}|=0$. Thus, the recovered folding
number $q_{i}$ will fall into one of the following forms in case 1):

  1. 1.

$q_{i}\equiv Y_{i}\mod lcm(M_{l},l\in\mathbb{U}_{a})$ and $q_{i}\equiv-
Y_{i}-1\mod lcm(M_{l},l\in\mathbb{U}_{e})$

  2. 2.

$q_{i}\equiv Y_{i}+1\mod lcm(M_{l},l\in\mathbb{U}_{c})$ and $q_{i}\equiv-
Y_{i}\mod lcm(M_{l},l\in\mathbb{U}_{d})$.

  3. 3.

$q_{i}\equiv Y_{i}\mod lcm(M_{l},l\in\mathbb{U}_{a})$ and $q_{i}\equiv-
Y_{i}\mod lcm(M_{l},l\in\mathbb{U}_{d})$.

Without loss of generality, we assume that at least $\lceil\frac{L}{N}\rceil$
residues in $S_{i}$ are from $\\{X_{1},-X_{1}\\}$, corresponding to the moduli
set $\mathbb{U}_{1}=\\{M_{1},M_{2},\dots,M_{l_{k}}\\}$, where
$l_{k}\geq\lceil\frac{L}{N}\rceil$. Based on the above understanding, $q_{i}$
must be in one of the forms, subject to the qualification of (36).

  * •

$q_{i}\equiv Y_{1}\mod lcm(M_{l},l\in\mathbb{U}_{1a})$ and $q_{i}\equiv-
Y_{1}-1\mod lcm(M_{l},l\in\mathbb{U}_{1e})$

  * •

$q_{i}\equiv Y_{1}+1\mod lcm(M_{l},l\in\mathbb{U}_{1c})$ and $q_{i}\equiv-
Y_{1}\mod lcm(M_{l},l\in\mathbb{U}_{1d})$

  * •

$q_{i}\equiv Y_{1}\mod lcm(M_{l},l\in\mathbb{U}_{1a})$ and $q_{i}\equiv-
Y_{1}\mod lcm(M_{l},l\in\mathbb{U}_{1d})$

We study the first situation first, where $q_{i}$ satisfies

$\left\\{\begin{array}[]{l}q_{i}\equiv Y_{1}\mod
lcm(M_{l}\in\mathbb{U}_{1a})\\\ q_{i}\equiv-Y_{1}-1\mod
lcm(M_{l}\in\mathbb{U}_{1e})\\\
q_{i}\equiv\lfloor\frac{\widetilde{r}_{(i,l^{\prime})}}{\Gamma}\rfloor+\tau_{(i,l^{\prime})}\mod
M_{l^{\prime}}\end{array}.\right.$ (38)

where $\mathbb{U}_{1a}\cup\mathbb{U}_{1e}=\mathbb{U}_{1}$ and
$l^{\prime}\in\\{l_{k}+1,l_{k}+2,\dots,L\\}$. Without loss of generality, we
assume that $lcm(M_{l}\in\mathbb{U}_{1a})>lcm(M_{l}\in\mathbb{U}_{1e})$. Then,
$q_{i}=klcm(M_{l}\in\mathbb{U}_{1a})+Y_{1}$. Clearly, when $Y_{1}\not\equiv-
Y_{1}-1\mod lcm(M_{l}\in\mathbb{U}_{1e})$, $q_{i}$ will contradict the dynamic
range assumed with $k\geq 1$. Only when $Y_{1}\equiv-Y_{1}-1\mod
lcm(M_{l}\in\mathbb{U}_{1e})$ or $lcm(M_{l}\in\mathbb{U}_{1e})=\emptyset$ is
satisfied will $q_{i}\in[0,D_{q})$, i.e., $q_{i}=Y_{1}$. So (38) changes to

$\left\\{\begin{array}[]{l}q_{i}\equiv Y_{1}\mod
lcm(M_{l}\in\mathbb{U}_{1})\\\
q_{i}\equiv\lfloor\frac{\widetilde{r}_{(i,l^{\prime})}}{\Gamma}\rfloor+\tau_{(i,l^{\prime})}\mod
M_{l^{\prime}}\end{array}\right.$ (39)

Thus, we have $q_{i}=klcm(M_{l}\in\mathbb{U}_{1})+Y_{1}$, where
$lcm(M_{l}\in\mathbb{U}_{1})>D_{q}$. Similarly, if at least one
$l^{\prime}\in[l_{k}+1,L]$ satisfies
$Y_{1}\not\equiv\lfloor\frac{\widetilde{r}_{(i,l^{\prime})}}{\Gamma}\rfloor+\tau_{(i,l^{\prime})}\mod
M_{l}^{\prime}$, $q_{i}$ will exceed the range $D_{q}$, leading to a
contradiction.

That is to say, $q_{i}=Y_{1}$ is the only efficient solution with at least
$\lceil\frac{L}{N}\rceil$ residues are from $\\{X_{1},-X_{1}\\}$ and
$lcm(M_{l}\in\mathbb{U}_{1a})>lcm(M_{l}\in\mathbb{U}_{1e})$. The shifted
common residues assigned to $q_{i}$ are $\hat{r}^{+}_{1c,l}$ and
$\hat{r}^{\pm}_{ic,l}$. Based on (36), we have
$|\hat{r}^{\pm}_{ic,l_{1}}-\hat{r}^{+}_{1c,l_{2}}|=|\hat{r}^{\pm}_{ic,l_{1}}-(r^{+}_{1c}+\Delta^{+}_{1c,l_{2}})|<\frac{\Gamma}{2}$,
i.e., $|\hat{r}^{\pm}_{ic,l_{1}}-r^{+}_{1c}|<\frac{3\Gamma}{4}$. Therefore,
$|\hat{X}_{1}-X_{1}|$ is equal to

$|Y_{1}\Gamma+\frac{\sum_{l\in\mathbb{U}_{1a}}\hat{r}^{+}_{1c,l}+\sum_{l\in\mathbb{U}_{1e}}\hat{r}^{-}_{1c,l}+\sum_{l\in\overline{\mathbb{U}}_{1}}\hat{r}^{\pm}_{ic,l}}{L}-X_{1}|<\frac{3\Gamma}{4}.$
(40)

In the following, we study the second situation, where $q_{i}\equiv
Y_{1}+1\mod lcm(M_{l},l\in\mathbb{U}_{1c})$ and $q_{i}\equiv-Y_{1}\mod
lcm(M_{l},l\in\mathbb{U}_{1d})$. Thus, (38) becomes

$\left\\{\begin{array}[]{l}q_{i}\equiv Y_{1}+1\mod
lcm(M_{l}\in\mathbb{U}_{1c})\\\ q_{i}\equiv-Y_{1}\mod
lcm(M_{l}\in\mathbb{U}_{1d})\\\
q_{i}\equiv\lfloor\frac{\widetilde{r}_{(i,l^{\prime})}}{\Gamma}\rfloor+\tau_{(i,l^{\prime})}\mod
M_{l^{\prime}}\end{array}\right.$ (41)

Similarly, under the assumption
$lcm(M_{l}\in\mathbb{U}_{1c})>lcm(M_{l}\in\mathbb{U}_{1d})$, $q_{i}=Y_{1}+1$
is the only efficient solution within $D_{q}$. So $|\hat{X}_{1}-X_{1}|$ is
equal to

$\displaystyle|(Y_{1}+1)\Gamma+\frac{\sum_{l\in\mathbb{U}_{1c}}\hat{r}^{+}_{1c,l}+\sum_{l\in\mathbb{U}_{1d}}\hat{r}^{-}_{1c,l}+\sum_{l\in\overline{\mathbb{U}}_{1}}\hat{r}^{\pm}_{ic,l}}{L}-X_{1}|$
(42)
$\displaystyle<|\lfloor\frac{X_{1}}{\Gamma}\rfloor\Gamma+\frac{\sum_{l\in\mathbb{U}}({r}^{+}_{1c}+\frac{3\Gamma}{4})}{L}-X_{1}|=\frac{3\Gamma}{4}.$

We can reach the same conclusion $|\hat{X_{i}}-X_{i}|<\frac{3\Gamma}{4}$ in
the last situation, which is omitted. Q.E.D. ∎

Algorithm 3 Robust Remaindering Decoding of Multiple Frequencies in Real
Waveform

Input: Moduli set: $\mathcal{M}=\\{m_{1}=M_{l}\Gamma|l=1,2,\dots,L\\}$, where
$M_{l}$ are sorted in ascending order;
Residue Sets: $\mathcal{R}_{l}=\\{\widetilde{r}^{\pm}_{i,l}|i=1,2,\dots,N\\}$,
$l=1,2,\dots,L$.

1: Repeat: Propose a clustering assignment

2: Following the proposed clustering by selecting one residue from each
$\mathcal{R}_{l}$ to obtain a $L$-residue clustering
$S_{i}=\\{\widetilde{r}_{(i,l)}|l=1,2,\dots,L\\}$, $i=1,2,\dots,N$, where
$(i,l)$ denotes the index of the residues assigned to $S_{i}$.

3: Assign a parameter $\tau_{(i,l)}\in\\{0,1\\}$ to each residue in $S_{i}$
randomly.

4: Calculate the shifted common residues
$\hat{r}_{(ic,l)}=\langle\widetilde{r}_{(i,l)}\rangle_{\Gamma}-\tau_{(i,l)}\Gamma$
corresponding to each $\widetilde{r}_{(i,l)}$ in $S_{i}$.

5: For each $S_{i}$, calculate
$q_{i}\equiv\frac{\widetilde{r}_{(i,l)}-\hat{r}_{(ic,l)}}{\Gamma}\mod M_{l}$
via CRT.

6: Until: There exist $N$ clusterings $S_{i}$ such that $q_{i}\in[0,D_{q})$
and each residue in $S_{i}$ satisfies (36) for $i=1,2,\dots,N$, where
$D_{q}=\min_{\mathbb{U}}\frac{lcm(M_{l}\in\mathbb{U})+lcm(M_{l}\in(\overline{\mathbb{U}})}{2}-1$
and
$\mathbb{U}\cup\overline{\mathbb{U}}=\\{M_{1},M_{2},\dots,M_{\lceil\frac{L}{N}\rceil}\\}$.

Output: $\hat{X_{i}}=q_{i}\Gamma+\frac{\sum^{L}_{l=1}\hat{r}_{(ic,l)}}{L}$.