On Unlimited Sampling and Reconstruction

TL;DR

This paper introduces the Unlimited Sampling Framework, a novel approach that enables perfect reconstruction of bandlimited signals from modulo samples, overcoming ADC saturation issues by using a folding ADC technique.

Contribution

It develops a new sampling paradigm using modulo operations, providing conditions and algorithms for stable, perfect signal recovery independent of ADC voltage limits.

Findings

Recovery is possible under specific sampling conditions.

The method is robust to bounded noise and quantization effects.

Numerical experiments confirm the effectiveness of the approach.

Abstract

Shannon's sampling theorem is one of the cornerstone topics that is well understood and explored, both mathematically and algorithmically. That said, practical realization of this theorem still suffers from a severe bottleneck due to the fundamental assumption that the samples can span an arbitrary range of amplitudes. In practice, the theorem is realized using so-called analog-to-digital converters (ADCs) which clip or saturate whenever the signal amplitude exceeds the maximum recordable ADC voltage thus leading to a significant information loss. In this paper, we develop an alternative paradigm for sensing and recovery, called the Unlimited Sampling Framework. It is based on the observation that when a signal is mapped to an appropriate bounded interval via a modulo operation before entering the ADC, the saturation problem no longer exists, but one rather encounters a different type of information loss due to the modulo operation. Such an alternative setup can be implemented, for example, via so-called folding or self-reset ADCs, as they have been proposed in various contexts in the circuit design literature. The key task that we need to accomplish in order to cope with this new type of information loss is to recover a bandlimited signal from its modulo samples. In this paper we derive conditions when perfect recovery is possible and complement them with a stable recovery algorithm. The sampling density required to guarantee recovery is independent of the maximum recordable ADC voltage and depends on the signal bandwidth only. Our recovery guarantees extend to measurements affected by bounded noise, which includes the case of round-off quantization. Numerical experiments validate our approach. Applications of the unlimited sampling paradigm can be found in a number of fields such as signal processing, communication and imaging.

Figures

Equations

$${\mathcal{M}}_{\lambda }:f\mapsto 2\lambda \left(\left[\left[\frac{f}{2\lambda }+\frac{1}{2}\right]\right]-\frac{1}{2}\right),$$

$$\left(\Delta f\right)\left[k\right]=f\left[k+1\right]-f\left[k\right]$$

$$f\in {\mathcal{B}}_{\Omega }\iff \hat{f}\left(\omega \right)=\text{mathbbmtt}{1}_{\left[-\Omega ,\Omega \right]}\left(\omega \right)\hat{f}\left(\omega \right)$$

$$g\in {\mathsf{PW}}_{\Omega },g\left(t\right)$$

$$0<{\color[rgb]{0,0,0}\mathrm{T}}\leqslant\frac{\pi}{\Omega},\ \ g\left(t\right)=\sum\limits_{k=-\infty}^{k=+\infty}{g\left({k{\color[rgb]{0,0,0}\mathrm{T}}}\right)\operatorname{sinc}\left({\frac{t}{{\color[rgb]{0,0,0}\mathrm{T}}}-k}\right)},$$

$$z\left(t\right)\stackrel{\mathrm{def}}{=}{\mathcal{M}}_{\lambda }\left(g\left(t\right)\right).$$

$${y[k]}\stackrel{{\scriptstyle\rm{def}}}{{=}}z\left({k{\color[rgb]{0,0,0}\mathrm{T}}}\right)=\mathscr{M}_{\lambda}\left(g\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right)\right),\ \ k\in\mathbb{Z}$$

$${\mathcal{C}}_{B,\lambda }=\left\{\pm \frac{\left(2n+1\right)}{2^B}\lambda |n\in \left\{0,\dots ,2^{B-1}-1\right\}\right\}.$$

$$g\left(t\right)=z\left(t\right)+{\epsilon }_g\left(t\right)$$

$${\epsilon }_g\left(t\right)=2\lambda \sum _{m\in \mathbb{Z}}e[m]{\text{mathbbmtt}1}_{{\mathcal{D}}_m}\left(t\right),e[m]\in \mathbb{Z},$$

$${\epsilon }_g\left(t\right)\stackrel{(\text{ref}eg)}{=}g\left(t\right)-{\mathcal{M}}_{\lambda }\left(g\left(t\right)\right)$$

$$\stackrel{(\mathrm{b})}{=}2\lambda \lfloor h\left(t\right)\rfloor ,$$

$$\left\lfloor h\left(t\right)\right\rfloor=\sum\limits_{\ell\in\mathbb{Z}}{{e[\ell]}{{\mathbbmtt{1}}_{{\mathcal{D}_{\ell}}}}\left(t\right)},\quad e[\ell]\in\mathbb{Z},\ \ \mathop{\vphantom{\bigcup}\mathchoice{\leavevmode\vtop{\halign{\hfil$\m@th\displaystyle#$\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil$\m@th\textstyle#$\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil$\m@th\scriptstyle#$\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil$\m@th\scriptscriptstyle#$\hfil\cr\bigcup\cr\cdot\crcr}}}}\displaylimits\limits_{\ell\in\mathbb{Z}}\mathcal{D}_{\ell}=\mathbb{R}$$

$$\mathsf{S}:{\left\{{s\left[k\right]}\right\}_{k\in{\mathbb{Z}^{+}}}}\to\sum\limits_{m=1}^{k}{\color[rgb]{0,0,0}s\left[m\right]}.$$

$$\|\Delta^{N}\gamma\|_{\infty}\leqslant{\left({{\color[rgb]{0,0,0}\mathrm{T}}e}\right)^{N}}\|g^{(N)}\|_{\infty}.$$

$$\|\Delta^{N}\gamma\|_{\infty}\leqslant{\left({{\color[rgb]{0,0,0}\mathrm{T}}\Omega e}\right)^{N}}\|g\|_{\infty}.$$

$$g\left(t\right)=P_{N-1}\left(t\right)+R_{N-1}\left(t\right)$$

$$P_{N-1}\left(t\right)=\sum _{n=0}^{N-1}\frac{g^{\left(n-1\right)}\left(\tau \right)}{\left(n-1\right)!}{\left(t-\tau \right)}^{\left(n-1\right)}$$

$$R_{N-1}\left(t\right)=\frac{g^{\left(N\right)}\left(z\right)}{N!}{\left(t-\tau \right)}^N$$

$$t_{n}=\left(k+n\right){\color[rgb]{0,0,0}\mathrm{T}}\in\left[{k{\color[rgb]{0,0,0}\mathrm{T}},\left({n+k}\right){\color[rgb]{0,0,0}\mathrm{T}}}\right],\quad n=0,\ldots,N.$$

$$\gamma\left[k\right]={P_{N-1}}\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right)+{R_{N-1}}\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right)$$

$$\left({{\Delta^{N}}\gamma}\right)\left[k\right]=\left({{\Delta^{N}}{R_{N-1}}}\right)\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right)$$

$${\Vert {\Delta }^N\gamma \Vert }_{\infty }$$

$$\displaystyle\leqslant\frac{{{{\left({N{\color[rgb]{0,0,0}\mathrm{T}}}\right)}^{N}}}}{{N!}}\|{g^{\left(N\right)}}|{|_{\infty}}$$

$$\displaystyle\leqslant\frac{{{{\left({{\color[rgb]{0,0,0}\mathrm{T}}e}\right)}^{N}}}}{{\sqrt{2\pi N}}}\|{g^{\left(N\right)}}|{|_{\infty}}\hfill\leqslant\left({{\color[rgb]{0,0,0}\mathrm{T}}e}\right)^{N}\|{g^{\left(N\right)}}|{|_{\infty}}$$

$$\Vert g^{\left(N\right)}|{|}_p⩽{\Omega }^N\Vert g|{|}_p,1⩽p⩽\infty$$

$${\beta }_g\in 2\lambda \mathbb{Z}\text{mbox}and\Vert g{\Vert }_{\infty }⩽{\beta }_g$$

$$\quad N\geqslant\left\lceil{\frac{{\log\lambda-\log\beta_{g}}}{{\log\left({{\color[rgb]{0,0,0}\mathrm{T}}\Omega e}\right)}}}\right\rceil.$$

$${\Delta }^N\gamma \equiv {\mathcal{M}}_{\lambda }\left({\Delta }^N\gamma \right).$$

$${\mathcal{M}}_{\lambda }({\Delta }^{N}a)={\mathcal{M}}_{\lambda }({\Delta }^N({\mathcal{M}}_{\lambda }(a)).$$

$${\mathcal{M}}_{\lambda }\left(a_1+a_2\right)={\mathcal{M}}_{\lambda }\left({\mathcal{M}}_{\lambda }\left(a_1\right)+{\mathcal{M}}_{\lambda }\left(a_2\right)\right),$$

$${\mathcal{M}}_{\lambda }\left({\Delta }^N{\mathcal{M}}_{\lambda }\left(a\right)\right)$$

Full text

In this paper, we develop an alternative paradigm for sensing and recovery, called the Unlimited Sampling Framework. It is based on the observation that when a signal is mapped to an appropriate bounded interval via a modulo operation before entering the ADC, the saturation problem no longer exists, but one rather encounters a different type of information loss due to the modulo operation. Such an alternative setup can be implemented, for example, via so-called folding or self-reset ADCs, as they have been proposed in various contexts in the circuit design literature.

The key task that we need to accomplish in order to cope with this new type of information loss is to recover a bandlimited signal from its modulo samples. In this paper we derive conditions when perfect recovery is possible and complement them with a stable recovery algorithm. The sampling density required to guarantee recovery is independent of the maximum recordable ADC voltage and depends on the signal bandwidth only. Our recovery guarantees extend to measurements affected by bounded noise, which includes the case of round-off quantization.

Numerical experiments validate our approach. For example, it is possible to recover functions with amplitudes orders of magnitude higher than the ADC’s threshold from quantized modulo samples upto the unavoidable quantization error.

Applications of the unlimited sampling paradigm can be found in a number of fields such as signal processing, communication and imaging.

1. I Introduction

2. I-A The Sampling Theorem and a Practical Bottleneck

3. I-B A Solution via Modular Arithmetic

4. I-C Contributions and Overview of Results

5. I-D Related Problems in Other Fields and Recent Work

6. I-E Organization of this Paper

7. I-F Notation

8. II Sr-Adcs and Modulo Samples

9. II-A Background on Self-Reset ADCs

10. II-B Mathematical Model for Unlimited Sampling

11. II-C The Structure of Discontinuities in Modulo Representation

12. III On Identifiability in Unlimited Sampling

13. III-A Uniqueness of Modulo Representation

14. IV Recovery Conditions and a Reconstruction Algorithm

15. IV-A Overview of the Recovery Scheme

16. IV-B Towards a Recovery Guarantee for Unlimited Sampling

17. IV-B1 Recovering Higher Order Differences from Modulo Samples

18. IV-B2 Recovery of Samples from Higher Order Differences

19. IV-B3 Recovery Conditions for the Case of Bounded Noise

20. V Numerical Experiments

21. V-A Verifying the Unlimited Sampling Theorem

22. V-B Exploring Sharpness of the Unlimited Sampling Theorem

23. V-C Quantization and Bounded Noise

24. VI Conclusion

25. VI-A Summary of Results

26. VI-B Future Directions

I Introduction

The interplay between the continuous and the discrete realms is at the heart of all modern information processing systems. Thus, the sampling theory inevitably finds its way in different spheres of science and engineering.

Since its early days, the topic has grown far and wide [5, 6, 7]. Researchers and practitioners have proposed a number of extensions. When we think of sampling theory, in most cases, variation on the theme arises from diversity along the time-dimension. This could be due to priors (sparsity or smoothness) or sampling architecture (uniform or non-uniform sampling grid). On the other hand, a hypothesis on the amplitude dimension leads to a different class of sampling problems, namely, level cross sampling, amplitude sampling, quantization and phaseless reconstruction. Our work falls into the latter category. In some sense, it can be seen as an amplitude analogue of the time-warped sampling theory [8], as in, our work exploits warping along the amplitude axis.

I-A The Sampling Theorem and a Practical Bottleneck

The key result of sampling theory is Shannon’s sampling theorem [7], which states that functions of bounded frequency range can be perfectly recovered from equidistant samples. To realize this theorem in practice, one uses the so-called analog–to–digital converter, henceforth ADC. Unlike the input of the sampling theorem which are assumed to be perfect samples, the output of an ADC is clipped at some fixed maximum recordable voltage; ADCs are limited by their dynamic range. This means that only values in a particular range, say, $\left[{-\lambda,\lambda}\right]$, may be recorded. We will refer to $\lambda>0$ as the ADC threshold. Whenever the input signal amplitude exceeds the threshold, that is $|f_{\sf{In}}|>\lambda$, the ADC saturates and essentially outputs $f_{{\sf{Out}}}=\lambda$. This is shown in Fig. 1 (a) and Fig. 1 (b). Such cases cannot be handled by sampling theory.

Clipping or saturation poses a serious problem in a variety of applications. For example, a camera facing the sun leads to an all white photograph which is the basis of high dynamic range photography [9]. This effect is particularly important in the context of self-driving vehicles. Scenarios such as moving out of a tunnel in the day [10] or a headlight in the line-of-sight can cause the imaging sensors to saturate. Similarly, in scientific imaging systems such as ultrasound [11], radar [12] and seismic imaging [13], strong reflections or pulse echoes blind the sensor. In communication systems, clipping results in performance degradation [14]. In the context of music, clipped sound results in high frequency artifacts [15].

Due to the pervasiveness of the clipping effect, several papers have studied this problem in different contexts (cf. [19, 11, 15, 20]) and particular emphasis has been put on audio signals (cf. [15] and references therein)—specific examples of bandlimited signals. In general, clipping of a smooth signal results in non-smooth features which in turn corresponds to high-frequency distortions in the signal. Hence, clipped or saturated signals are prone to aliasing artifacts [21].

On one hand, purely computational approaches involve de-clipping or restoration algorithms [22, 23, 19, 20, 21] with the hope of recovering lost samples under certain assumptions. Such approaches suffer from two basic drawbacks:

• •

Recovery guarantees in terms of the sampling density are largely unexplored. Ideally, the required sampling rate should be completely governed by the signal bandwidth.

• •

The mapping is discrete to discrete rather than incorporating sampling theory and following the goal of minimizing the error of the reconstructed continuous signal.

On the other hand, hardware-only approaches tackle the dynamic range problem at the ADC level [24, 25, 26]. This results in sophisticated electronic architecture that is application specific and is agnostic to the benefits offered by computational and algorithmic methods.

I-B A Solution via Modular Arithmetic

Unlike the literature discussed above which relies either entirely on computational approaches or only optimizes the hardware, our work is based on a co-design approach; we aim to overcome the dynamic range barrier by repurposing hardware in conjunction with new recovery algorithms.

The key innovation of our approach is that instead of (potentially clipped) pointwise samples of the bandlimited function, we work with folded amplitudes with values in the range $\left[{-\lambda,\lambda}\right]$. Mathematically, this folding corresponds to injecting a non-linearity in the sensing process. This amounts to,

[TABLE]

where $\left[\kern-1.49994pt\left[f\right]\kern-1.49994pt\right]\stackrel{{\scriptstyle\rm{def}}}{{=}}f-\left\lfloor f\right\rfloor$ defines the fractional part of $f$ and $\lambda>0$ is the ADC threshold. Note that (1) is equivalent to a centered modulo operation since $\mathscr{M}_{\lambda}(f)\equiv f~{}\rm{mod}~{}2\lambda$. By implementing the mapping (1), it is clear that out-of-range amplitudes are folded back into the dynamic range $\left[{-\lambda,\lambda}\right]$.

To connect this mathematical conceptualization to real world applications, we capitalize on recent advances in the imaging and sensor design technology. Indeed, ADC architectures that are moving towards the goal of folded sampling are rapidly developing. We have dedicated Section II-A to the discussion of such ADCs. In essence, the so-called self-reset [16] or folding [27] ADCs implement folding of amplitudes via (1) using electronic circuitry. We compare the transfer function of the conventional ADC with the self-reset ADC (henceforth Sr-Adc) in Fig. 1 (a) and Fig. 1 (c), respectively. Folding the amplitudes ensures that the entire range of the ADC is utilized (cf. Fig. 1 (c) and Fig. 1 (d)). To give the reader an idea of the functionality of this new breed of ADCs, we show the raw output of the prototype111We thank Sasagawa [17], Yamaguchi [18] and co-workers for providing the data which allowed us to reproduce these images. Sr-Adc arising in the context of imaging in Fig. 1 (e) and Fig. 1 (f).

While the recent decade has seen remarkable progress on the hardware aspects of the new ADCs, theoretical and algorithmic aspects has not been a major research focus and did not make their way to other fields. Current hardware literature [16, 17, 18] employs a fairly elementary approach for recovery that uses both the modulo samples and the residuals or reset counts222As per current approaches, any function (signal or image) can be written as a sum of the folded function and the reset count map. By recording both the folded function as well as the reset count map simultaneously, the original function can be recovered by a straight-forward sum of the two entities. For example, the image in Fig. 1 (g) is obtained by adding Fig. 1 (e) and (f). (see Fig. 1 (g)). This requires,

• •

complex circuitry and ADC architectures as well as,

• •

additional power and storage.

Furthermore, the feasibility of the reset count for arbitrarily small modulo thresholds has not been investigated yet.

In view of this discussion, we aim for an approach that does not suffer from these shortcomings and recovers the signal without knowledge of residuals or reset-counts. In full generality, this is a very difficult problem, which may also explain the limited progress. Namely, without further assumption on the underlying image or signal, the problem is closely related to the phase unwrapping problem, which is known to be highly ill-posed [28] and typically relies on inversion of the first-order finite difference (cf. Section I-D and Itoh’s condition.). Indeed, in both cases, one seeks to “unwrap” a discrete function representation from modulo information. Consequently, algorithms proposed for phase unwrapping also yield solution strategies for the signal unfolding problem [29, 30].

As in the case of phase unwrapping, however, these approaches suffer from strong limitations of the dynamic range (cf. Fig. 6 for a numerical comparison). For more information on the phase unwrapping problem, we refer to the literature overview in Section I-D.

An important difference, however, is that for the scenario studied in this paper, one has a considerably larger degree of control of the data entering the sensing pipeline. In particular, it is possible to sense redundant information, which in many other setups has been shown to alleviate the ill-posedness and allow for guaranteed recovery. Examples include sigma-delta quantization of bandlimited functions [31] and compressed sensing [32, 33]. In analogy to these works, our goal is to explore redundant representations that allow to overcome the limitations of the conventional viewpoint of phase unwrapping.

I-C Contributions and Overview of Results

The goal of this paper is to pose and study the inverse problem of recovering a continuous-time bandlimited function from its noisy folded samples without requiring the knowledge of residuals or reset-counts. Our key contributions are as follows:

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

We take a first step towards a sampling theory for modulo samples. To this end,

• •

We establish that a finite energy bandlimited function is identifiable by its modulo samples even when the sampling rate is just slightly above Nyquist.

• •

We present a sampling theorem which describes sufficient conditions for sampling and reconstruction of bandlimited functions in the context of folded samples.

• •

We show that this sampling theorem includes stability with respect to bounded noise thus covering scenarios with round-off quantization. 2. $\textrm{C}_{2})$

On the algorithmic front, our sufficiency condition is complemented by a constructive recovery algorithm that is guaranteed to recover the underlying signal and is stable with respect to noise both in theory and in experiments.

Our main contribution leads to the Unlimited Sampling Theorem which is summarized below.

Theorem** (Unlimited Sampling Theorem).**

Let $g\left(t\right)$ be a finite energy, bandlimited signal with maximum frequency $\Omega$ and let ${y[k]},k\in\mathbb{Z}$ in (1) be the modulo samples of $g\left(t\right)$ with sampling rate {\color[rgb]{0,0,0}\mathrm{T}}. Then a sufficient condition for recovery of $g\left(t\right)$ from $\{y[k]\}_{k}$ is that {\color[rgb]{0,0,0}\mathrm{T}}\leqslant\frac{1}{2\Omega e} (up to additive multiples of $2\lambda$) where $e$ denotes Euler’s constant.

Remarkably, our theorem requires a sampling rate depending on the bandwidth only, independent of the ratio between ADC threshold, $\lambda$, and signal amplitude. This is why we refer to our method as unlimited sampling. Our numerical demonstrations in Section V clearly corroborate that it is possible to recover signals whose amplitude range significantly exceeds the dynamic range of the ADC under consideration, that is, $\max|f_{\sf{In}}|\gg\lambda$.

I-D Related Problems in Other Fields and Recent Work

To put our results into context, we briefly compare and contrast topics where the modulo operation arises naturally. A disparate survey of the literature is presented below.

1. $\bullet$

Phase Unwrapping. The phase unwrapping problem333We thank Prof. Laurent Jacques (UC Louvain) and Prof. Alan Oppenheim (MIT) for bringing the topic of phase unwrapping to our notice and clarifying the intricate differences with respect to sampling theory. [34, 35, 36], which finds widespread applications, specially in imaging, arises in a number of applications related to imaging (cf. [35] as well as Chapter 3 in [3]) where one only has access to a sinusoid of varying phase shift and the relevant information is encoded in these shifts, which are assumed to be varying smoothly. Due to the periodicity of the sinusoid, one can only extract this information modulo $2\pi$, so the goal is to recover a smooth phase function that corresponds to these observations. In this sense, the phase unwrapping problem is closely related to the problem studied in this paper, although it arises in quite a different way; the ambiguity is inherent in the sensing problem whereas we deliberately inject this non-linearity to enhance the reconstruction quality. This computational design aspect becomes important as it allows for further modifications such as the introduction of redundancy.

Given this similarity it should not come as a surprise that one of the cornerstone methods for phase unwrapping can be seen as a simplified version of the algorithm proposed in this paper. More precisely, when the max-norm of the first-order finite difference of the samples is bounded by $\lambda$ or, $\left|{\gamma\left[{k+1}\right]-\gamma\left[k\right]}\right|\leqslant\lambda,\gamma\left[k\right]=g\left({kT}\right)$, one can recover by inverting the first-order finite difference operator on the modulo sequence. This is an established result and is widely known as Itoh’s condition [34]. Please refer to Fig. 2 for a schematic explanation of Itoh’s condition and Fig. 6 for numerical comparison with our approach. Due to the redundancy via oversampling in our sensing setup, we are able to advance this method both in terms of admissible signal amplitudes and stability. 2. $\bullet$

Digital Communications. In the work dating back to early 1970s, Tomlinson [37] and Harashima [38] describe the use of modulo precoding444We thank Prof. Robert Gray (Boston University) who shared several historical facts regarding the role of modulo operations in the context of precoding and offering clarifications with reference to his work on modulo sigma-delta modulation that is mentioned in Section II-A. in the context of communication over channels with inter-symbol interference. In recent years, modulo operation has also been used in integer-forcing communication [39]. Although the role of such strategies in connection with sampling theory is unclear, there may be potentially interesting links for specific signal structures and reconstruction algorithms for the case of folded samples. 3. $\bullet$

Rate Distortion Analysis and Prediction Based Recovery. As early as in 1979, Ericson & Ramamoorthy [40] considered quantized modulo samples in the context of a source coding scheme they referred to as Modulo-PCM. Here, the authors, propose a heuristic decoder based on the Viterbi algorithm for the case of first-order Gauss-Markov processes. In early 90’s, Chou & Gray [41] studied the behaviour of quantization noise in the context of modulo sigma-delta modulation. More recently, Boufounos [42] developed a recovery scheme for modulo quantized measurements that is based on consistent reconstruction. However, all of these approaches are only feasible in limited dimensions and quite costly when the dimension is large. Our approach, in contrast, is stable and efficient and is easily scalable to a large number of samples in a mathematically guaranteed fashion. The approach in [42] can be made more practicable with the knowledge of additional side-information [43] via prediction based recovery. Following the initial ideas underlying this work as summarized in [1], a modulo sampling based simulated hardware implementation and quantization of bandlimited functions was studied by Ordentlich *et al. *in [44], again using additional side information. The authors considered oversampled representation for both single and multiple channel cases. Later, in [45], Romanov & Ordentlich studied recovery for non-quantized samples and showed that sampling faster than Nyquist rate is enough under some mild decay conditions. The side information used in [44] and [45] consists of a certain number of non-modulo samples based on which linear prediction using Chebyshev polynomials [46] can be employed to extrapolate the bandlimited function. In contrast to all these works, such side-information is not required for the method presented in this paper. This is of particular advantage when only a finite number of measurements is available, for instance, in the case of parametric signals [47, 48] as one can not expect that modulo non-linearity will leave any patch, unaffected. 4. $\bullet$

Modular Exponentiation. Consider the problem of computing and storing a number $a^{b}\bmod c$ where $a,b,c$ are all integers. This problem frequently arises in number theory, computer science and cryptography [49]. The literature in this area focuses on $a^{b}\bmod c$ without computing $a^{b}$ as the digital storage required by the former is much smaller than the latter. In this spirit, our work is similar to modular exponentiation as the ADC threshold $\lambda$ (in analogy to $c$) controls the amount of information/bits required for representing the dynamic range. 5. $\bullet$

Contemporary Literature (2018 onwards). During the completion of this manuscript, our first work [1] was followed up in several papers. Below, we summarize the key results.

1. —

Rudresh *et al. *[50] study a wavelet based algorithm for reconstruction of Lipschitz continuous functions in the context of unlimited sensing framework. 2. —

Cucuringu & Tyagi [51] investigated a more general setup based on Hölder continuous functions. They also presented a denoising approach based on a quadratically constrained quadratic program (QCQP) with non-convex constraints. 3. —

Compressed sensing of discrete-time sparse signals within the unlimited sampling architecture was investigated by Musa *et al. *in [52]. The authors developed a generalized approximate message passing approach to reconstruct discrete signals. Further bounds in context of Gaussian matrices were studied in [53]. 4. —

Unlimited sampling of continuous-time sparse signals in canonical and Fourier domain, was discussed in our works [47] and [48], respectively. 5. —

The idea of one-bit unlimited sampling was proposed by Graf *et al. *in [54]. This approach recovers high dynamic range signals from signed measurements thus adding to the functionality of the sigma-delta based ADCs. 6. —

Given modulo measurements of $n$ i.i.d samples drawn from a Gaussian distribution with unknown covariance matrix, Romanov & Ordentlich studied recovery methods in [55]. The case of known covariance matrix was covered in the work of Ordentlich *et al. *[44]. 7. —

In [56], the authors extend the idea of unlimited sensing for signals that live on a graph. To this end, the authors propose an integer programming based algorithm for recovery from modulo samples. 8. —

Modulo Radon Transform was presented in [57] and unlimited sampling based high dynamic range computed tomography was presented in [58]. 9. —

In recent works, [59, 60], the authors presented a multi-channel unlimited sampling approach in the context of sensor array signal processing.

I-E Organization of this Paper

The paper is organized as follows. In Section II, we discuss ADC architectures that lead to modulo samples. Section III is devoted to the study of uniqueness conditions. Recovery conditions and the reconstruction algorithm are presented in Section IV. Numerical examples that corroborate our theory are discussed in Section V. We conclude this work with future directions in Section VI.

I-F Notation

We use $\mathbb{N}$, $\mathbb{Z}$, $\mathbb{R}$ and $\mathbb{C}$ to denote the set of natural numbers, integers, reals and complex numbers, respectively. We use \mathop{\vphantom{\bigcup}\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\bigcup\cr\cdot\crcr}}}}\displaylimits to denote disjoint union of sets. Continuous functions are written as $f\left(t\right),t\in\mathbb{R}$ while their discrete counter-parts are represented as $f\left[k\right],k\in\mathbb{Z}$. The $L^{p}$ space equipped with the $p$-norm or $\|\cdot\|_{L^{p}}$ is the standard Lebesgue space. For instance, $L^{1}$ and $L^{2}$ denote the space of absolute and square-integrable functions, respectively. Spaces associated with sequences will be denoted by $\ell^{p}$. The max-norm ($p\to\infty$) of a function is defined as, ${\left\|f\right\|_{{\infty}}}=\inf\left\{{{c_{0}}\geqslant 0:\left|{f\left(t\right)}\right|\leqslant{c_{0}}}\right\}$ while for sequences, we have, ${\left\|f\right\|_{{\infty}}}={\max_{k}}\left|{f\left[k\right]}\right|$. The $N$–order derivative of a function is denoted by $f^{\left(N\right)}\left(t\right)$. The space of $N$–times differentiable, real-valued functions is denoted by $C^{N}\left(\mathbb{R}\right)$. For sequences,

[TABLE]

denotes the finite difference and recursively applying the same yields $N^{\textrm{th}}$ order difference, $\left({\Delta^{N}f}\right)\left[k\right]$. The Fourier transform of a function $f\in L^{1}$ is denoted by $\widehat{f}\left(\omega\right)=\int f\left(t\right)e^{-\jmath\omega t}dt$ with a natural extension to arbitrary distributions via duality. We say $f\in C^{N}\left(\mathbb{R}\right)$ is $\Omega$-bandlimited or,

[TABLE]

where $\mathbbmtt{1}_{\mathcal{D}}\left(t\right)$ is the indicator function on domain $\mathcal{D}$. This includes functions with infinite energy such as a pure sine wave. When dealing with finite-energy functions, it is convenient to work with the Paley-Wiener space, $\mathsf{PW}_{\Omega}=\mathcal{B}_{\Omega}\cap{L^{2}}$. Note that for $\epsilon>0$, one has $\mathsf{PW}_{\Omega}\subset\mathsf{PW}_{\Omega+\epsilon}$. In the context of numerical computations, $\left\lceil\cdot\right\rceil$ and $\left\lfloor\cdot\right\rfloor$ denote the floor and ceiling functions.

II Sr-Adcs and Modulo Samples

II-A Background on Self-Reset ADCs

When reaching the upper or lower saturation voltage, the Sr-Adcs reset to the respective threshold. This is what allows for capturing voltage variations beyond the conventional saturation limit. This mechanism aptly justifies the name self-reset or folding ADCs which is mathematically equivalent to using the modulo mapping in (1).

As is often the case, theoretical conceptualization of such ADCs predated its practical implementation. To bring the reader up to pace with the literature, we will quickly review the key references in the area. As early as in the late 1970’s, the concept of modulo limiters was used in the context of Tomlinson–Harashima decoders [37, 38] in communication theory. Chou and Gray [41] studied the resulting quantization noise but, neither the physical realization nor their recovery properties were investigated. It was only in the early 2000’s that the physical implementation of such ADCs started to develop. One of the earliest references is due to Rhee and Joo [16] where the authors proposed the Sr-Adc in context of CMOS imagers. Interested readers are also referred to a tutorial article on folding ADCs by Kester [27] which contains historical references. In a study on quantitative characterization of ADC architectures, Kavusi and Gamal [61] note that Sr-Adcs allow for simultaneous enhancement of both dynamic range as well as the signal-to-noise ratio. We remark that the revolution around the Sr-Adcs is mainly motivated by the goal of achieving improvement on the dynamic range. This is primarily because the dynamic range of natural scenes is typically much larger than what a conventional ADC or an imaging sensor may accommodate. While HDR imaging requires several exposures, Sr-Adc based imaging is inherently HDR due to folding of amplitudes. An alternative architecture for HDR imaging via modulo operations was recently used in [28]. Beyond consumer photography [9] and autonomous vehicle navigation [10], improvement on dynamic range is important for scientific and bio-imaging. To this end, Sasagawa [17] and Yamaguchi [18] recently developed an implantable Sr-Adc for functional brain imaging.

II-B Mathematical Model for Unlimited Sampling

Thanks to the Sr-Adcs [16, 17, 18], we can repurpose the resetting capability for obtaining modulo samples without recording the reset counts (cf. Section I-B and Fig. 1 (f)). The sampling process for obtaining modulo samples of a function is outlined in the block diagram in Fig. 3. The basic principles are similar to the conventional case except for the modulo mapping. A break down of the key steps is as follows:

1. $1)$

We start with a square-integrable function $f\in L^{2}$ (not-necessarily bandlimited) to be sampled. 2. $2)$

Pre-filtering of $f$ with $\phi\in\mathsf{PW}_{\Omega}$ results in a low-pass approximation which is given by,

[TABLE]

Also note that since $g\in\mathsf{PW}_{\Omega}$, it admits the standard sampling formula for bandlimited functions,

[TABLE]

where {\color[rgb]{0,0,0}\mathrm{T}} is the sampling rate and $\operatorname{sinc}\left(t\right)=\tfrac{{\sin\left({\pi t}\right)}}{{\left({\pi t}\right)}}$ is the ideal low-pass filter. 3. $3)$

The bandlimited function $g$ is folded in the range $\left[{-\lambda,\lambda}\right]$ via non-linear mapping (1) and results in,

[TABLE] 4. $4)$

The folded function $z\left(t\right)$ is sampled using impulse-train, {\otimes_{k{\color[rgb]{0,0,0}\mathrm{T}}}}\stackrel{{\scriptstyle\rm{def}}}{{=}}\sum\limits_{n\in\mathbb{Z}}{\delta\left({t-k{\color[rgb]{0,0,0}\mathrm{T}}}\right)}, with sampling rate {\color[rgb]{0,0,0}\mathrm{T}}>0 yielding uniform samples,

[TABLE]

as shown in Fig. 3. 5. $5)$

When considering quantization with a budget of $B$ bits per sample, each modulo sample $y[k]$ is rounded to the closest element in the set

[TABLE]

Note that this operation induces noise that is not following a random model, but is entry-wise bounded by $2^{-B}$, which is why our guarantee below considers bounded deterministic noise.

We plot $g\left(t\right)$, $z\left(t\right)$ and modulo samples ${\left\{{{y[k]}}\right\}_{k\in\mathbb{Z}}}$ in Fig. 1.

II-C The Structure of Discontinuities in Modulo Representation

Clearly, the unlimited sampling architecture converts a smooth function into a discontinuous one. Recovering the unfolded function $g\left(t\right)$ from scrambled, low dynamic-range, samples then boils down to patching the discontinuities together. In fact, it turns out that the discontinuities admit a structure which is critical to our cause. Every bandlimited function, continuous or discrete, can be decomposed as a sum of a modulo function and a step-wise residual that we call simple function. This is shown in Fig. 4. We elaborate on this aspect in form of the following proposition.

Proposition 1** (Modular Decomposition Property).**

Let $g\in\mathcal{B}_{\Omega}$ and $\mathscr{M}_{\lambda}\left(\cdot\right)$ be defined in (1) where $\lambda$ is a fixed, positive constant. Then, the bandlimited function $g\left(t\right)$ admits a decomposition

[TABLE]

where $z\left(t\right)=\mathscr{M}_{\lambda}\left(g\left(t\right)\right)$ and ${\varepsilon_{g}}$ is a simple function,

[TABLE]

where \mathop{\vphantom{\bigcup}\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\bigcup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\bigcup\cr\cdot\crcr}}}}\displaylimits_{m\in\mathbb{Z}}\mathcal{D}_{m}=\mathbb{R} is a partition of the real line into intervals $\mathcal{D}_{m}$.

Proof.

Since $z\left(t\right)=\mathscr{M}_{\lambda}\left(g\left(t\right)\right)$, by definition, we write,

[TABLE]

where (a) is due to $h\stackrel{{\scriptstyle\rm{def}}}{{=}}\left(g/2\lambda\right)+1/2$ and in (b), we use $h=\left[\kern-1.49994pt\left[h\right]\kern-1.49994pt\right]+\left\lfloor h\right\rfloor$. Since, for an arbitrary function $h$, $\left\lfloor h\right\rfloor$ has the form,

[TABLE]

we obtain the desired result. $\blacksquare$

This proposition is the stepping stone towards the recovery algorithm. It allows us to tackle two problems at once:

1. $1)$

According to (5), if $\varepsilon_{g}$ is known, we can recover $g$ from $z$. In this paper we develop a method which allows for inferring $\varepsilon_{g}$ from $z$. 2. $2)$

According to (8), the fact that the amplitudes of $\varepsilon_{g}$ can only be integer multiples of $2\lambda$ allows us to enforce a consistency constraint555By the consistency constraint, we mean that in presence of noise or rounding errors, we may force the amplitudes of $\varepsilon_{g}$ to be on the grid of $2\lambda\mathbb{Z}$. This is implemented using (20) and step $4(b)$ in Algorithm 1. For numerical demonstration see Section V-C. which in turn leads to a robust recovery algorithm.

We emphasize that the decomposition property of Proposition 1 offers a general purpose utility whose potential benefits are typically not exploited. For instance, phase unwrapping [34] may directly benefit from this observation. In phase unwrapping literature, one recovers $g(nT)$ from its first-order finite difference using the anti-difference. Clearly, an arbitrary smooth function has far more degrees of freedom than (8) and hence, from robustness perspective, it is more effective to enforce the simple function structure.

III On Identifiability in Unlimited Sampling

In this section, we present identifiability conditions associated with the unlimited sampling strategy. Said differently, we answer the question of what sampling density leads to unique characterization of a bandlimited function in terms of its modulo samples? The key result of this section is that any sampling rate faster than critical sampling, that is, 0<{\color[rgb]{0,0,0}\mathrm{T}}<\pi/\Omega allows for a one-to-one mapping between a bandlimited function, $f\in\mathsf{PW}_{\Omega}$ and modulo samples {y[k]}=\mathscr{M}_{\lambda}\left(g\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right)\right). In what follows, without loss of generality, we will normalize the bandwidth such that $\Omega=\pi$ and the critical sampling rate is {\color[rgb]{0,0,0}\mathrm{T}}=1. In working with oversampled representations, we will denote the sampling rate by, {\color[rgb]{0,0,0}\mathrm{T}_{\epsilon}}=\frac{\pi}{\pi+\epsilon},\epsilon>0 implying that 0<{\color[rgb]{0,0,0}\mathrm{T}_{\epsilon}}<1.

III-A Uniqueness of Modulo Representation

The following Lemma is used for the proof of injectivity conditions.

Lemma 1**.**

Assume that $f\in\mathsf{PW}_{\pi}$, $\epsilon>0$ and $\mathbb{L}\subset\mathbb{Z}$ is some finite set. Then $f$ is uniquely characterized by its samples in {\color[rgb]{0,0,0}\mathrm{T}_{\epsilon}}\cdot(\mathbb{Z}\setminus\mathbb{L}) where {\color[rgb]{0,0,0}\mathrm{T}_{\epsilon}}=\frac{\pi}{\pi+\epsilon},\epsilon>0 is the sampling rate.

For proof of this Lemma, we refer the reader to [2]. Lemma 1 is used to prove that oversampling uniquely determines modulo samples. The formal statement is as follows.

Theorem 1** (Injectivity Theorem for Unlimited Sampling).**

Any $f\in\mathsf{PW}_{\pi}$ is uniquely determined by its modulo samples on the grid {\left\{{{t_{n}}=n{{\color[rgb]{0,0,0}\mathrm{T}_{\epsilon}}}}\right\}_{n\in\mathbb{Z}}} with $\epsilon>0$.

For proof of this Theorem, we refer the reader to [2]. A subsequent alternative form of proof can be found in [45].

IV Recovery Conditions and a Reconstruction Algorithm

The result of the previous section shows that any sampling rate faster than the Nyquist rate uniquely characterizes a bandlimited function in terms of its modulo samples. However, this does not yield a constructive algorithm. This section is concerned with the development of a constructive algorithm that recovers a bandlimited function from its modulo samples as defined in (4). Our algorithm is accompanied with a recovery guarantee. Akin to Shannon’s sampling theorem, our recovery guarantee purely depends on the signal bandwidth.

IV-A Overview of the Recovery Scheme

Our basic strategy for recovering functions from modulo samples is schematically explained in Fig. 5. The key observation at the heart of our approach is that finite differences commute with the modulo operation in a certain sense. With the first order difference given by, $\left({\Delta y}\right)\left[k\right]\stackrel{{\scriptstyle\rm{def}}}{{=}}y\left[{k+1}\right]-y\left[k\right]$, the $N^{\textrm{th}}$ order difference can be obtained by recursive application of the finite-difference operator, ${\Delta^{N}}y={\Delta^{N-1}}\left({\Delta y}\right)$. Starting with modulo samples $y[k]$ (cf. Fig. 5-(1)-b), its higher order differences are given by $\left({{\Delta^{N}}y}\right)\left[k\right]$ (cf. Fig. 5-(2)-b). Under appropriate conditions, as shown in Fig. 5-(3)-b, applying the modulo operation on $\left({{\Delta^{N}}y}\right)\left[k\right]$ reduces to $\left({{\Delta^{N}}g}\right)\left[k\right]$ and we obtain the equivalence $\mathscr{M}_{\lambda}\left({{\Delta^{N}}y}\right)\left[k\right]=\left({{\Delta^{N}}g}\right)\left[k\right]$. Equivalently, the residual function is annihilated that is, $\mathscr{M}_{\lambda}\left(\Delta^{N}\varepsilon_{g}\right)=0$. To eventually recover $g\left[k\right]$ from $\mathscr{M}_{\lambda}\left({{\Delta^{N}}y}\right)\left[k\right]=\left({{\Delta^{N}}g}\right)\left[k\right]$, we will work towards recovery of $\varepsilon_{g}$ which can be robustly estimated as it takes amplitudes on the grid of $2\lambda\mathbb{Z}$. To do so, we compute, $\left({{\Delta^{N}}y}\right)-\mathscr{M}_{\lambda}\left(\Delta^{N}y\right)=\Delta^{N}\varepsilon_{g}.$ For $N=1$, this is shown in Fig. 5-(2)-c. To recover $\varepsilon_{g}$ from \color[rgb]{0,0,0}{\Delta^{N}}\varepsilon_{g}, we will use the anti-difference operator defined as,

[TABLE]

followed by rounding of \mathsf{S}\left(\color[rgb]{0,0,0}{\Delta^{N}}\varepsilon_{g}\right) to the nearest multiple of $2\lambda\mathbb{Z}$ which adds to the stability. Although polynomials are in the kernel of the anti-difference operator, we will develop an approach that allows us to estimate $\varepsilon_{g}$ up to an unknown constant. Finally, having estimated the residual function $\widetilde{\varepsilon}_{g}$ up to an integer multiple of $2\lambda$, we obtain the bandlimited samples $g\left[k\right]$ by using (5). The continuous-time function is obtained by low-pass filtering the samples using (2).

IV-B Towards a Recovery Guarantee for Unlimited Sampling

Our scheme relies on conditions under which the higher order differences of the bandlimited samples are small enough in amplitude so that the modulo non-linearity has no effect on the same. To analyze how much can one shrink the amplitudes, we will begin our analysis with a bound that relates max-norm of higher order differences with its continuous-time counterpart, the derivative. Below, we summarize some well-known consequences of Taylor’s theorem and Bernšteĭn’s inequality in form of the following lemma.

Lemma 2**.**

For any $g\in C^{m}(\mathbb{R})\cap L^{\infty}(\mathbb{R})$ and $N\in\mathbb{N}$, its samples \gamma[k]\stackrel{{\scriptstyle\rm{def}}}{{=}}g(k{\color[rgb]{0,0,0}\mathrm{T}}) satisfy

[TABLE]

Furthermore, whenever $g\left(t\right)$ is an $\Omega$–bandlimited function, its samples \gamma[k]\stackrel{{\scriptstyle\rm{def}}}{{=}}g\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right),{\color[rgb]{0,0,0}\mathrm{T}}\leqslant\pi/\Omega satisfy,

[TABLE]

Proof.

Using Taylor’s theorem, we can write,

[TABLE]

where ${P_{N-1}}\left(t\right)$ is the Taylor polynomial of degree $N-1$ around $\tau$,

[TABLE]

and the remainder of the Taylor series in Lagrange form is given by,

[TABLE]

where $z$ is some number between $\tau$ and $t$. Given a sequence of samples \gamma[k]=g\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right), its higher order differences $\left({{\Delta^{N}}\gamma}\right)\left[k\right]$ use $N$ contiguous values of $g\left(t_{n}\right)$ where,

[TABLE]

By letting $\tau$ to be the center of the above and hence \tau=\left(k+\frac{N}{2}\right){\color[rgb]{0,0,0}\mathrm{T}} and sampling (12) on both sides, we can now write

[TABLE]

and consequently,

[TABLE]

because $\Delta^{N}$ annihilates polynomials of degree $N-1$. Noting that ${\left\|{{\Delta^{N}}{R_{N}}}\right\|_{\infty}}\leqslant{2^{N}}{\left\|{{R_{N}}}\right\|_{\infty}}$ and using (13), we have,

[TABLE]

where the last but one inequality is obtained by applying the Stirling’s approximation, that is, $N!\sim\sqrt{2\pi N}{\left({\frac{N}{e}}\right)^{N}}$ and this proves (10).

Next, we use the well known Bernšteĭn’s inequality to bound the right hand side of (10). For $\Omega$–bandlimited functions $g\in\mathsf{PW}_{\Omega}$, the Bernšteĭn’s inequality (cf. pg. 116, [62]) asserts that,

[TABLE]

and by plugging the same in (10), we obtain (11). $\blacksquare$

IV-B1 Recovering Higher Order Differences from Modulo Samples

The result of Lemma 2 and in particular, the inequality in (11) will be the key to our proposed recovery method. By letting {\color[rgb]{0,0,0}\mathrm{T}}<1/\Omega e and choosing $N$ logarithmically in $\|g\|_{\infty}$, the right hand side of (11) can be made to shrink arbitrarily, ultimately ensuring that $\|{\Delta^{N}}\gamma|{|_{\infty}}<\lambda$. More precisely, assuming that we know some constant $\beta_{g}$ such that,

[TABLE]

a suitable choice for $N$ that satisfies {\left({{\color[rgb]{0,0,0}\mathrm{T}}\Omega e}\right)^{N}}\beta_{g}\leqslant\lambda is,

[TABLE]

For this choice of $N$ and {\color[rgb]{0,0,0}\mathrm{T}}\leqslant{\left({2\Omega e}\right)^{-1}}, (11) entails that $\left({{\Delta^{N}}\gamma}\right)\left[k\right]$ is unaffected by the modulo operation, that is,

[TABLE]

The following proposition relates the right hand side of (15) to the modulo samples $y[k]$ in (4).

Proposition 2**.**

For any sequence $a[k]$ it holds that

[TABLE]

Proof.

In view of the Modular Decomposition Property, we write, $a=\mathscr{M}_{\lambda}\left(a\right)+\varepsilon_{a}$ to obtain $\Delta^{N}\mathscr{M}_{\lambda}\left(a\right)=\Delta^{N}{a}-\Delta^{N}\varepsilon_{a}$ where $\varepsilon_{a}$ is a simple function. With the observation that

[TABLE]

we obtain,

[TABLE]

Now since $\varepsilon_{a}$ in (6) take values in $2\lambda\mathbb{Z}$ and the finite difference filter has integer valued coefficients, it follows that $\Delta^{N}\varepsilon_{a}\in 2\lambda\mathbb{Z}$. Consequently, $\Delta^{N}\varepsilon_{a}$ is in the kernel of $\mathscr{M}_{\lambda}\left(\cdot\right)$ or $\mathscr{M}_{\lambda}\left(\Delta^{N}\varepsilon_{a}\right)=0$. Hence,

[TABLE]

which proves the result in (16). $\blacksquare$

By letting $a=\gamma$ in (16), we obtain, $\mathscr{M}_{\lambda}\left(\Delta^{N}\gamma\right)=\mathscr{M}_{\lambda}\left(\Delta^{N}y\right)$. Combining this with (15) yields,

[TABLE]

Hence, starting with modulo samples $y$ in (4), we are able to relate the same with higher order differences of bandlimited samples $\gamma$.

IV-B2 Recovery of Samples from Higher Order Differences

The result of Proposition 1 yields that,

[TABLE]

where the sequence $e$ satisfies $e[m]\in\mathbb{Z}$ for all $m\in\mathbb{Z}$. Hence given modulo samples $y[k]$, recovering \gamma[k]=g\left(k{\color[rgb]{0,0,0}\mathrm{T}}\right) is equivalent to recovering $\varepsilon_{\gamma}[k]$. Thanks to Proposition 2, with $N$ in (14), we can use the relation in (15) to estimate \color[rgb]{0,0,0}{\Delta^{N}{\varepsilon_{\gamma}}} directly from the modulo samples. This is because,

[TABLE]

It now remains to recover $\varepsilon_{\gamma}$ from \color[rgb]{0,0,0}\Delta^{N}{\varepsilon}_{\gamma}. We remind the reader that \color[rgb]{0,0,0}\Delta^{N}{\varepsilon}_{\gamma}\in 2\lambda\mathbb{Z} and this massive restriction on the range of values that \color[rgb]{0,0,0}\Delta^{N}{\varepsilon}_{\gamma} can take makes the recovery procedure considerably less ill-posed than recovering $\gamma$ from \color[rgb]{0,0,0}{\Delta^{N}{\gamma}}. Algorithmically, recursively applying the anti-difference operator $\mathsf{S}$ in (9) and then rounding the same to the nearest multiple of $2\lambda$, that is,

[TABLE]

recovers \color[rgb]{0,0,0}{\Delta^{N-1}}\varepsilon_{\gamma} up to an additive constant. The remaining ambiguity is due to the constant sequence being in the kernel of the first order finite difference operator. The unknown constant can again only take values in $2\lambda\mathbb{Z}$. More precisely,

[TABLE]

where $l[k]=2\lambda\mathbb{Z}$ is the constant sequence and $\kappa_{n}\in\mathbb{Z}$. Clearly, when $n=1$, the ambiguity cannot be resolved as $\gamma[k]+2\lambda\mathbb{Z}$ leads to the same modulo samples. However, for ${1<n\leqslant N}$, we can resolve this ambiguity. By applying $\mathsf{S}$ to (21), we obtain,

[TABLE]

where $\left(\mathsf{S}l\right)[k]={{{\left\{{2\lambda k}\right\}}_{k\in{\mathbb{Z}}}}}$ is a linear sequence with slope $2\lambda$. Next, observe that,

[TABLE]

where $(\mathrm{a})$ uses the modulo decomposition property (5). Furthermore, $\|y\|_{\infty}\leqslant\lambda$ implies that $\left({{\Delta^{n-2}}y}\right)\leqslant{2^{n-2}}\lambda$ and from Lemma 2, (11), we have \|{\Delta^{n-2}}\gamma{\|_{\infty}}\leqslant{\left({{\color[rgb]{0,0,0}\mathrm{T}}\Omega e}\right)^{n-2}}\beta_{g}. Hence,

[TABLE]

From the above inequality, it is clear that oversampling or {\color[rgb]{0,0,0}\mathrm{T}}<1/\Omega e keeps the max-norm of the higher order differences bounded. Using this with (23), we conclude,

[TABLE]

where,

[TABLE]

Here, the last inclusion follows as {\color[rgb]{0,0,0}\mathrm{T}}\leqslant\left(2\Omega e\right)^{-1}. Using $n\leqslant N$ and {\color[rgb]{0,0,0}\mathrm{T}}\leqslant\left(2\Omega e\right)^{-1}, again, we have

[TABLE]

With (25), the above inequality yields ${\rho_{n}}\left({\lambda,\beta_{g}}\right)=\frac{3\beta_{g}}{2\lambda J}$. Let us define the sequence,

[TABLE]

Setting ${\rho_{n}}\left({\lambda,\beta_{g}}\right)=\frac{3\beta_{g}}{2\lambda J}$ in (24) entails,

[TABLE]

For $J={{6\beta_{g}}}/{\lambda}$, the right hand side in the above is an interval of length ${1}/{2}$ and hence only contains one integer, namely,

[TABLE]

Together with the rounding operation defined in (20), the estimate of ${\kappa_{n}}$ allows us to recursively compute,

[TABLE]

where the sequence $s_{\left(n\right)}[k]$ is initialized with initial condition s_{\left({0}\right)}[k]\stackrel{{\scriptstyle(\ref{eq:egy})}}{{=}}\color[rgb]{0,0,0}\Delta^{N}{\varepsilon}_{\gamma}[k]. This completes the recovery procedure which is summarized in Algorithm 1. The following theorem provides a recovery guarantee for this algorithm.

Theorem 2** (Unlimited Sampling Theorem).**

Let $g\left(t\right)\in\mathsf{PW}_{\Omega}$ and {y[k]}={\left.{\mathscr{M}_{\lambda}\left({g\left(t\right)}\right)}\right|_{t=k{\color[rgb]{0,0,0}\mathrm{T}}}},k\in\mathbb{Z} in (4) be the modulo samples of $g\left(t\right)$ with sampling rate {\color[rgb]{0,0,0}\mathrm{T}}. Then a sufficient condition for recovery of $g\left(t\right)$ from the $\{y[k]\}_{k}$ up to additive multiples of $2\lambda$ is that

[TABLE]

Provided that this condition is met and assuming that $\beta_{g}\in 2\lambda\mathbb{Z}$ is known with $\|g\|_{\infty}\leqslant\beta_{g}$, then choosing

[TABLE]

yields that {\left({{\color[rgb]{0,0,0}\mathrm{T}}\Omega e}\right)^{N}}{\left\|{g}\right\|_{\infty}}<\lambda and Algorithm 1 recovers $g$ from $y$ again up to the ambiguity of adding multiples of $2\lambda$.

The proof of this Theorem is deferred to the next section (Section IV-B3), which considers the case of bounded noise. The same proof applies in the absence of noise.

IV-B3 Recovery Conditions for the Case of Bounded Noise

In this section, we consider the effect of bounded noise on the recovery procedure. More precisely, we assume that the modulo samples $y[k]$ are affected by noise $\eta$ of amplitude bounded by a constant $b_{0}>0$. That is,

[TABLE]

Note that due the presence of noise, it may happen that {\color[rgb]{0,0,0}y_{\eta}}[k]\not\in[-\lambda,\lambda]. Nonetheless, for $b_{0}$ below some fixed threshold, our recovery method provably recovers noisy bandlimited samples $\gamma[k]$ from the associated noisy modulo samples {\color[rgb]{0,0,0}y_{\eta}}[k] up to an unknown additive constant, where the noise appearing in the recovered samples is in entry-wise agreement with the one affecting the modulo samples. That is, $\widetilde{\gamma}\left[k\right]=\gamma\left[k\right]+\eta\left[k\right]+2m\lambda,m\in\mathbb{Z}$ .

Theorem 3** (Unlimited Sampling Theorem with Noise).**

Let $g\left(t\right)\in\mathsf{PW}_{\Omega}$ and assume that $\beta_{g}\in 2\lambda\mathbb{Z}$ is known with $\|g\|_{\infty}\leqslant\beta_{g}$. For the dynamic range we work with the normalization $\mathsf{DR}={\beta_{g}}/{\lambda}$.

Let the noisy modulo samples be of the form (30) with a noise bound given in terms of the dynamic range as

[TABLE]

Then a sufficient condition for Algorithm 1 with

[TABLE]

to approximately recover the bandlimited samples $\gamma[k]$ in the sense that it returns $\widetilde{\gamma}\left[k\right]=\gamma\left[k\right]+\eta\left[k\right]+2m\lambda,m\in\mathbb{Z}$ is that

[TABLE]

Proof.

We start by observing that the modified choice of $N$ as given in (32) yields that ${(T\Omega e)^{N}}\leqslant{\left({{\lambda}/{{{(2\beta_{g})}}}}\right)}$. Together with (11), this yields that $\|\Delta^{N}\gamma\|_{\infty}\leqslant\frac{\lambda}{2}$. Similarly, the modified choice of $T$ yields that

[TABLE]

which yields via Young’s inequality and (31) that

[TABLE]

These observations now entail that even the noisy modulo measurements {\color[rgb]{0,0,0}y_{\eta}} allow for exact computation of $\varepsilon_{\gamma}[k]$. Indeed, we observe that

[TABLE]

Here in (a), we used Proposition 2. Equality (b) follows from (17) and (34). For (c) we again use (34) combined with the bound $\|\Delta^{N}\gamma\|_{\infty}\leqslant\frac{\lambda}{2}$ derived above.

We will show now by induction in $m\in[0,N-1]$ that ${s_{\left(m\right)}}={\Delta^{N-m}\varepsilon_{\gamma}}$ for ${s_{\left(m\right)}}$ as in (27) and $\varepsilon_{\gamma}$ as in (18). The choice $m=N-1$ then ensures that ${\mathsf{S}s_{\left(N-1\right)}}=\varepsilon_{\gamma}+2k\lambda$ for some $k\in\mathbb{Z}$. As per modulo decomposition in Proposition 1, this ensures that step 5 of Algorithm 1 returns noisy bandlimited samples $\widetilde{\gamma}[k]=\gamma[k]+2m\lambda+\eta[k]$, that is, one recovers up to the measurement noise and an unavoidable constant ambiguity of $2\lambda\mathbb{Z}$. The result then follows from Shannon’s sampling theorem which is implemented in step 6 of Algorithm 1.

To prove that ${s_{\left(m\right)}}={\Delta^{N-m}\varepsilon_{\gamma}},\quad m\in[0,N-1]$, we first note that the induction seed $m=0$ reduces to showing that s_{\left({0}\right)}=\color[rgb]{0,0,0}{\Delta^{N}{\varepsilon_{\gamma}}}, which holds by definition.

For the induction step, assume that for some $m\in[0,N-2]$, one has ${s_{\left(m\right)}}={\Delta^{N-m}\varepsilon_{\gamma}}$. Recall that we have derived above that for $J=6\beta_{g}/\lambda$ as chosen in Algorithm 1, $\kappa_{m}$ given by (26) is the unique multiple of $2\lambda$ for which $s_{\left(m+1\right)}[k]$ defined in (27) is a bounded sequence. Hence, combining (21) with the induction hypothesis and the observation that ceiling and floor in (27) only have an effect in case of round-off errors, we conclude that ${s_{\left({m+1}\right)}}={\Delta^{N-\left({m+1}\right)}}\varepsilon_{\gamma}$ as desired. $\blacksquare$

V Numerical Experiments

In this section, we numerically verify several results linked with the topic of unlimited sampling and reconstruction. In particular, (a) we verify our sampling theorem and explore its limitations, (b) compare Algorithm 1 with Itoh’s method for phase unwrapping and (c) confirm that our approach is stable in the presence of quantization noise. These results are plotted in Fig. 6 (a)–(d) and Fig. 7 (a)–(c).

V-A Verifying the Unlimited Sampling Theorem

In order to verify our sampling theorem, we generate $1000$ realizations of a bandlimited function $g\in\mathcal{B}_{\pi}$ with piecewise constant Fourier spectrum taking values chosen uniformly at random, $\widehat{g}\left(\omega\right)\in\mathrm{unif}\left(0,1\right)$ where $\mathrm{unif}\left(a,b\right)$ denotes the uniform random distribution. The bandlimited function is normalized such that $||g||_{\infty}=1$. Furthermore, for each realization, we use $\lambda\in\mathrm{unif}\left(\frac{1}{100},\frac{1}{10}\right)$ and $\beta_{g}=\left\lfloor{{{\left\|g\right\|}_{\infty}}/2\lambda}\right\rfloor$. We obtain modulo samples $y[k]$ in (4) using sampling rate {\color[rgb]{0,0,0}\mathrm{T}}=\frac{11}{200}<\frac{1}{2\pi e}. There on, we use Algorithm 1 to recover bandlimited samples. Although Algorithm 1 recovers bandlimited samples up to an unknown additive constant of $2\lambda\mathbb{Z}$, we use the knowledge of ground truth to estimate this unknown constant so that we can compute the mean squared error. For each of the $1000$ realizations, the samples $\widetilde{\gamma}$ are recovered up to machine precision ($10^{{-33}}$).

A specific realization of the experiment is shown in Fig. 6 (a). With $\lambda=809/3125,\beta_{g}=20\lambda$ and {\color[rgb]{0,0,0}\mathrm{T}}={11}/{200}, we estimate $N=5$. The recovered samples $\widetilde{\gamma}$ are also shown in Fig. 6 (a). The mean squared error between the bandlimited samples and recovered samples in this case is $6.2455\times 10^{-33}$. The residue function corresponding to this realization, that is $\varepsilon_{\gamma}$, is computed using the ground truth and its estimated version which is obtained using Algorithm 1. This is shown in Fig. 6 (b). We also compare our reconstruction with the phase unwrapping method based on Itoh’s condition [34]. As explained in Fig. 2, Itoh’s condition succeeds when $|\Delta y|<\lambda$, which is not the case here. As shown in Fig. 6 (c), phase unwrapping based reconstruction fails.

V-B Exploring Sharpness of the Unlimited Sampling Theorem

On one hand, the injectivity condition in Theorem 1 proves that any sampling rate above the Nyquist rate guarantees unique representation of a bandlimited function in terms of its modulo samples. On the other hand, for Algorithm 1 to succeed, the sampling rate must be a factor of $2\pi e$ faster than Nyquist rate. Here we set up a demonstration that shows that Algorithm 1 succeeds even when the sampling rate is much slower than what is prescribed by Theorem 2, however, the recovery is not guaranteed. For any $g\in\mathcal{B}_{\pi}$, we denote the critical sampling rate by {\color[rgb]{0,0,0}\mathrm{T}_{\sf{Shannon}}}=1 and the corresponding sampling rate prescribed by the unlimited sampling theorem is {\color[rgb]{0,0,0}\mathrm{T}_{\sf{US}}}=1/2\pi e. We consider modulo samples acquired with $\lambda=2/10$, so the theorem is valid for all orders $N\geqslant 1$. To assess the sharpness of Theorem 2, we use $1000$ realizations of randomly generated bandlimited functions with $||g||_{\infty}=1$. Each realization is then sampled with sampling rate {\color[rgb]{0,0,0}\mathrm{T}}=0.01,\ldots,{\color[rgb]{0,0,0}\mathrm{T}_{\sf{Shannon}}} in steps of $1/100$ and we set $\beta_{g}=1$. For sampling rates beyond the applicability of the theorem, that is, {\color[rgb]{0,0,0}\mathrm{T}}>{\color[rgb]{0,0,0}\mathrm{T}_{\sf{US}}}, we run Algorithm 1 for each of $1000$ realizations by varying from $N=1$ to $N=5$. Then for each combination of {\color[rgb]{0,0,0}\mathrm{T}} and $N$, we compute the success rate which is defined as the fraction of trials in which the algorithm reconstructs up to machine precision (in the sense of mean-squared error). As shown in Fig. 6 (d), indeed the recovery is possible even when {\color[rgb]{0,0,0}\mathrm{T}}\in\left({\color[rgb]{0,0,0}\mathrm{T}_{\sf{US}}},{\color[rgb]{0,0,0}\mathrm{T}_{\sf{Shannon}}}\right). Furthermore, higher order differences extend recover to smaller over-sampling factors, but there seems to be a minimum oversampling rate around {\color[rgb]{0,0,0}\mathrm{T}}=0.4 above which the algorithm always fails for the chosen value of $\lambda$.

V-C Quantization and Bounded Noise

As discussed in Section IV-B3, Theorem 3 provides recovery guarantees in the presence of bounded noise. A main motivation for this noise model is that it includes the error caused by round-off quantization as in Section II-B, Step 5. In Fig. 7 we illustrate the recovery guarantees for $3$–bit quantized modulo samples. That is, each modulo measurement is rounded to the closest element of

[TABLE]

In other words, the quantized measurements $y_{\eta}$ are incurring quantization noise $\|\eta\|_{\infty}\leqslant\tfrac{\lambda}{8}$. In the experiment, we use a bandlimited function with $||g||_{\infty}=12.5$ and $\lambda=1$. As before, {\color[rgb]{0,0,0}\mathrm{T}}={11}/{200} with $\beta_{g}=14\lambda$ and we set $N=2$. The experiment shows that our method is even more robust than predicted by Theorem 3: despite the quantization noise $\left\|\eta\right\|_{\infty}\approx\lambda/8$, Algorithm 1 recovers the underlying (noisy) bandlimited samples which are shown in Fig. 7 (b). The reconstruction mean squared error (MSE) between $\gamma$ and its estimate $\widetilde{\gamma}$ (using Algorithm 1) is is observed to be $5.1\times 10^{-3}$ which is the same as the MSE between $y$ and {\color[rgb]{0,0,0}y_{\eta}} and much smaller than the MSE of $3\times 10^{-1}$ incurred with $3$-bit quantization of non-modulo samples.

VI Conclusion

VI-A Summary of Results

In this work, we have described a new sensing and recovery scheme that overcomes the dynamic range barrier fundamental to digital sensing and sampling theory. Our approach harnesses a co-design between hardware and algorithms. On the hardware front, a modulo non-linearity injected in the sensing process scrambles high-dynamic-range information into low-dynamic range samples. Theorem 1 addresses the question of uniqueness; a bandlimited function is uniquely characterized by its modulo samples provided that the sampling rate is faster than critical sampling density. On the algorithmic front, given modulo samples of an $\Omega$-bandlimited function, the sampling rate {\color[rgb]{0,0,0}\mathrm{T}}\leqslant 1/\left(2\Omega e\right) guarantees recovery. This Unlimited Sampling Theorem is formally stated in Theorem 2; stability with respect to quantization noise and other bounded noise is established in Theorem 3. The proposed reconstruction approach hinges on the observation that the modulo and the (higher-order) difference operators commute in a certain sense.

VI-B Future Directions

1. $\bullet$

Implementation Issues Implementing our approach in practice remains an integral aspect of future research. In particular, we are working to design computational techniques to address inaccuracies and instabilities arising in the implementation of our techniques in hardware. We are currently exploring such instabilities in the context of a prototype ADC. Important challenges include timing jitter and that our sensing pipeline requires ADCs without anti-aliasing filter, which makes them more sensitive to noise. For this reason, noise control needs to be implemented in alternative ways. 2. $\bullet$

Wider Classes of Signals and Function Spaces A natural question to ask is how our results can be extended to a wider class of signals such as smooth functions (shift-invariant spaces), sparse signals, parametric classes of signals among others. Some partial answers are provided in [50, 51, 52, 53, 47, 48] but this is an open area with several unanswered questions. Unlimited sampling in spline-spaces was recently discussed in [63]. Another interesting variation is to consider multi-dimensional signal models defined on arbitrary lattices, such as [64]. 3. $\bullet$

Wider Classes of Inverse Problems The unlimited sensing framework leads to a new class of inverse problems due to the modulo non-linearity. It can be further combined with several interesting inverse problems where dynamic range poses a natural limitation. More generally, one may consider $y_{n}=\left(\mathscr{M}_{\lambda}\circ\mathscr{T}\right)\left(x_{n}\right)$ where $\mathscr{T}$ is some operator. For instance, using $\mathscr{T}$ as the Radon Transform, recently we have introduced the Modulo Radon Transform in the context of high-dynamic-range computed tomography [57]. 4. $\bullet$

New Sampling Architectures The unlimited sensing framework can be combined with different sampling architectures. For instance, in [54], one-bit recovery is proposed based on the sigma-delta modulation scheme. Similarly, this work can be combined with multi-channel acquisition approaches which naturally find applications in imaging, for example, time-of-flight imaging [3] and our upcoming work on direction-of-arrival estimation from modulo samples [59]. 5. $\bullet$

Robustness to Outliers While our results yield recovery for bounded noise, our approach is not yet stable with respect to outliers. Given the redundancy of the data, we expect that it will be possible to address this issue and, subsequently, study unbounded noise models such as Gaussian and Poissonian noise.