# Modulo Sampling in Shift-Invariant Spaces: Recovery and Stability
Enhancement

Yhonatan Kvich, Yonina C. Eldar This research was supported by the Tom and
Mary Beck Center for Renewable Energy as part of the Institute for
Environmental Sustainability (IES) at the Weizmann Institute of Science, by
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant No. 101000967) and by the Israel
Science Foundation (grant No. 536/22).

###### Abstract

Sampling shift-invariant (SI) signals with a high dynamic range poses a
notable challenge in the domain of analog-to-digital conversion (ADC). It is
essential for the ADC’s dynamic range to exceed that of the incoming analog
signal to ensure no vital information is lost during the conversion process.
Modulo sampling, an approach initially explored with bandlimited (BL) signals,
offers a promising solution to overcome the constraints of dynamic range. In
this paper, we expand on the recent advancements in modulo sampling to
encompass a broader range of SI signals. Our proposed strategy incorporates
analog preprocessing, including the use of a mixer and a low-pass filter
(LPF), to transform the signal into a bandlimited one. This BL signal can be
accurately reconstructed from its modulo samples if sampled at slightly above
its Nyquist frequency. The recovery of the original SI signal from this BL
representation is then achieved through suitable filtering. We also examine
the efficacy of this system across various noise conditions. Careful choice of
the mixer plays a pivotal role in enhancing the method’s reliability,
especially with generators prone to instability. Our approach thus broadens
the framework of modulo sampling’s utility in efficiently recovering SI
signals, pushing its boundaries beyond BL signals while sampling only slightly
above the rate needed for a SI signal.

###### Index Terms:

Modulo sampling, dynamic range, unlimited sampling, $B^{2}R^{2}$.

##  I Introduction

Analog-to-digital converters (ADCs) transform analog signals into a digital
format for processing in digital signal processing systems. The cost and power
requirements of ADCs escalate with an increase in the sampling rate, making it
preferable to operate at the minimum necessary rate for efficient sampling, as
highlighted in previous studies [1, 2]. The Shannon-Nyquist sampling theorem
is commonly employed in this context, asserting that bandlimited (BL) signals
can be accurately represented by their samples at a rate at least twice the
maximum frequency present in the signal. As the sampling rate increases, so do
the cost and power consumption associated with the analog-to-digital
conversion process. Therefore, sampling as close as possible to the signal’s
Nyquist rate is beneficial. Another important factor is the dynamic range of
ADCs. To avoid signal clipping and the consequent loss of information, as
illustrated in Fig. 1(a), the dynamic range of an ADC must surpass that of the
input analog signal.

Figure 1: (a) BL signal and its clipped version. (b) BL signal and output of
its modulo operation.

Numerous approaches have been developed to address clipping or to enhance the
dynamic range, with classification based on the presence or absence of
preprocessing prior to sampling. Techniques that forgo preprocessing exploit
the correlation inherent in BL signals’ samples taken above the Nyquist rate,
using this correlation to reconstruct any missing information due to clipping
by oversampling [3, 4]. Alternatively, leveraging the spectral vacancies
inherent in multiband communication systems can distinguish between original
and clipped signals [5, 6]. These methods, however, either necessitate
significant oversampling [3, 4] or presuppose knowledge of the spectral gaps
[5, 6], without offering theoretical assurances for their efficacy. Clipping
can also be circumvented through attenuation, though this approach risks
diminishing low-amplitude signals beneath the noise floor. Variable gain
attenuators, like automatic gain controls (AGCs) and companders, adjust to
preserve signal integrity without disproportionately affecting signals of
smaller amplitude. AGCs utilize a series of amplifiers regulated by feedback
to maintain a consistent output level [7, 8], while companding adjusts gain
inversely proportional to the signal amplitude. Companding, however, like
clipping, introduces nonlinear distortion and expands the signal’s bandwidth.
Beurling’s theorem and subsequent methodologies by Landau et al. have
facilitated the recovery of BL signals from their companded versions by
ensuring the compander’s output for finite energy inputs remains finite and
the compander function is monotone and differentiable across the input
signal’s dynamic range [9, 10]. Despite these advancements, the implementation
of companding-based solutions at minimal sampling rates faces challenges due
to requirements for monotonicity, differentiability, and finite energy output
[9, 10].

An alternate strategy involves applying a modulo operation to the input signal
before sampling to constrain its dynamic range, a technique that has found
application in high-dynamic-range ADCs, or self-reset ADCs, in imaging
contexts [11, 12, 13, 14]. This approach, alongside storing modulo signal
samples, often involves capturing additional data, such as the extent of
folding for each sample or the folding direction, complicating the sampling
circuitry while simplifying signal reconstruction from folded samples. See
Fig. 1(b) for a visual representation.

The concept of unlimited sampling, introduced by Bhandari et al., relies
solely on folded or modulo samples for signal recovery, demonstrating that
sampling above the Nyquist rate enables unique identification of BL signals
from their modulo samples [15]. This method, which extends Itoh’s unwrapping
algorithm, shows that a sufficient oversampling rate allows for the
computation of the original signal’s higher-order differences from modulo
samples, facilitating signal reconstruction through summation of these
differences. However, this technique’s effectiveness diminishes in noisy
environments, necessitating a significantly higher oversampling rate for
reliable recovery [15]. Subsequent improvements by Romanov and Ordentlich and
others have sought to reduce the required sampling rate, apply the technique
to various signal models, and explore hardware implementations for high-
dynamic-range ADCs using modulo operations [16, 17, 18, 19, 20, 21, 22, 23,
24, 25]. Despite these advancements, challenges remain in terms of missing
theoretical guarantees, stability concerns, the need for smooth and monotone
operators, and reliance on higher-than-Nyquist sampling rates.

Azar et al. [26, 27] introduced a recovery algorithm capable of reconstructing
BL signals from their modulo samples at a rate slightly above the Nyquist
rate, even under various noise conditions. This approach was further refined
in [28] through the incorporation of a sparsity assumption to enhance noise
robustness. Additionally, Mulleti et al. [29] explored the application of
modulo sampling to finite-rate-of-innovation (FRI) signals. The methodologies
developed for both BL and FRI signals have been successfully implemented in
hardware, as detailed in [30].

While the assumption of bandlimitation often serves as a practical
approximation, numerous signals exhibit more precise representations through
alternative bases or possess distinct structures within the Fourier domain
[31, 32, 20, 19, 22, 33, 29]. Shift-invariant (SI) spaces, in particular, hold
a crucial role in the theory of sampling. Signals within these spaces are
expressed as linear combinations of shifts from a collection of generating
functions [1, 32, 34, 35, 36, 37]. Any such signal can be accurately
reconstructed in an SI space, formed by shifting $k$ functions with period
$T$, using $k$ distinct sampling sequences. These samples are derived by
passing the signal through a filter bank consisting of $k$ filters and
sampling the output uniformly at intervals of $nT$. This technique allows for
a sampling rate of $k/T$, potentially much lower than the Nyquist rate for
signals with wideband generators, given certain conditions [32].

A unified theory for modulo sampling within SI spaces remains undeveloped.
Bhandari and Krahmer [24] considered 2D B-splines, theirs methodology
necessitates an oversampling rate determined by the modulo operation’s
requirements and the signal’s amplitude limits, and lacks noise robustness.

This study introduces an innovative recovery approach for SI signals utilizing
modulo sampling, which extends upon the methodology presented in [38]. To
enhance stability, we integrate a mixer into the signal processing framework.
This integration is designed to exploit the full spectral range of the signal
generator, significantly bolstering the stability of the recovery process. By
initially applying a low-pass filter (LPF) to transform the input signal into
a BL form suitable for modulo recovery, as described in prior work [26, 27],
and subsequently leveraging the characteristics of the SI space for signal
reconstruction, our approach requires only a slightly higher sampling rate
than that necessary for SI signals. The key improvement facilitated by the
mixer not only enables the efficient utilization of the entire spectrum of the
generator but also effectively addresses and mitigates potential stability
concerns, resulting in a more resilient and efficient signal recovery
mechanism.

We perform simulations with Lorentzian generators and achieve successful
signal recovery amidst noise interference. Our results reveal enhanced noise
robustness through the integration of the mixer. We also looked at scaled
B-spline generators, which typically present difficulties in SI recovery due
to inherent instability. Our observations confirm that the addition of a mixer
substantially improves both stability and the accuracy of recovery in noisy
environments, highlighting the effectiveness of our method in enhancing signal
reconstruction for various generator models. Following this, we apply modulo
sampling to Electrocardiogram (ECG) signals and successfully achieve recovery,
demonstrating the applicability of our method to biological signals.

The structure of this paper is outlined as follows: Section II introduces the
foundational concepts of modulo sampling, SI spaces, and the formulation of
the problem. Section III details the design of the recovery system, including
a theoretical framework for signal reconstruction. In Section IV, we
illustrate the system’s performance through simulations, particularly focusing
on its behavior in noisy conditions. Section V explores the application of SI
modulo sampling in the context of ECG signal analysis. Finally, Section VI
offers concluding remarks.

Throughout this paper, the following notations are used. The space of
sequences with finite norm is represented as $\ell^{2}$, and its corresponding
norm is denoted $\|\cdot\|_{2}$. For a given sequence $a[n]$, its Discrete-
Time Fourier Transform (DTFT) is
$A(e^{j\omega}):=\sum_{n\in\mathbb{Z}}a[n]e^{-j\omega n}$. Similarly, for a
function $x(t)$ with finite norm, its Continuous-Time Fourier Transform (CTFT)
is defined as $X(\omega):=\int_{t\in\mathbb{R}}x(t)e^{-j\omega t}$, we also
employ the notation $\mathcal{F}\\{x\\}(\omega)$. Denote $\ast$ as
convolution.

##  II Background and Problem Statement

###  II-A Preliminaries

The dynamic range of ADCs presents a significant challenge in accurately
capturing signals. An issue that often arises is clipping, which occurs when
the signal values exceed the ADC’s dynamic range, resulting in a loss of
information. One approach to circumvent this problem is to expand the dynamic
range, but this leads to increased quantization errors and necessitates the
use of higher-resolution ADCs, which in turn consume more power. An
alternative proposed solution is the application of a nonlinear operator to
the analog signal before sampling. This paper primarily focuses on the modulo
operator as the non-linear operator of interest, diverging from the
alternatives explored in [26, 27]. The modulo operator works by reducing the
dynamic range of a signal, “folding” it within a predetermined interval. For a
given $\lambda>0$, the operator maps real values to the interval
$[-\lambda,\lambda]$ as follows:

$\mathcal{M}_{\lambda}x:=((x+\lambda)\mod 2\lambda)-\lambda.$ (1)

Azar et al. [26, 27] have shown that by sampling slightly above the Nyquist
rate, it is feasible to “unfold” the samples of a BL signal. They also
developed an algorithm for signal recovery in such contexts, called
$B^{2}R^{2}$. Fig 2 illustrates the overall framework of the modulo recovery
process for BL signals. An input BL signal is first filtered through $s(t)$,
processed by the analog modulo operator. Finally, the signal is then sampled
at a rate $T_{s}<T$. Here, $T$ denotes the Nyquist rate, which is defined as
twice the highest frequency found within the bandwidth of $x(t)$. Various
other methods for recovering BL signals have been proposed, such as the one in
[28], which utilizes sparsity assumptions and fast-ISTA (FISTA) to improve the
method in [26, 27]. Techniques in [15, 16] and [25] have introduced
methodologies that use higher-order differences and iterative signal sieving,
respectively, for signal reconstruction. However these approaches require
higher-than-Nyquist sampling rates and are less robust to noise.

Figure 2: The modulo recovery system for BL signals, as introduced in [26, 27,
28].

This study aims to expand the findings of [26, 27, 28] to the domain of SI
signals [1]. A SI space with a single generator is defined as the span of a
generator signal $h(t)$ with time-shifts of $T$, meaning

$\mathcal{H}:=\Big{\\{}\sum_{n\in\mathbb{Z}}{a[n]h(t-nT)}|a\in\ell^{2}\Big{\\}}.$
(2)

Choosing $h(t)$ as the sinc function with a width of $\frac{1}{2T}$, the space
in (2) comprises $2\Omega$-BL signals, where $\Omega=\frac{1}{2T}$.

It is a well-established fact [1] that signals in this space can be
successfully recovered from samples taken at a rate of $\frac{1}{T}$,
independent of the Nyquist rate of the generator or whether it is BL. The
recovery algorithm involves frequency domain division using a correction
filter, the design of which is contingent upon the properties of the SI space,
the sampling rate, and the analog filter applied before sampling.

The domain of SI spaces offers increased flexibility in terms of signal
sampling, leading to a broader range of potential applications. The objective
here is to initiate a comprehensive study on modulo sampling of SI signals,
aiming for robust recovery while sampling at a rate close to $\frac{1}{T}$.

###  II-B Problem Formulation and Recovery Approach

Our goal is to reconstruct a SI signal

$x(t)=\sum_{n\in\mathbb{Z}}{a[n]h(t-nT)}\in\mathcal{H}$ (3)

from its modulo samples $\\{\mathcal{M}_{\lambda}x[nT_{s}]\\}$, where
$f_{s}=\frac{1}{T_{s}}$ is the sampling frequency, at least $\frac{1}{T}$,
nearing its value. Notably, $T$ can be significantly lower than the Nyquist
rate of both $h(t)$ and $x(t)$. To facilitate this, we first apply analog pre-
processing before implementing the modulo operator. Our objective is to
reconstruct the input signal from the modulo samples
$\\{\mathcal{M}_{\lambda}f(x)[nT_{s}]\\}$, with $f(\cdot)$ denoting the
necessary analog domain pre-processing. These pre-processing steps are
critical for ensuring unique recovery and augmenting the efficiency of the
modulo sampling method. We will demonstrate that modulo signal recovery is
achievable for $T_{s}<T$, for a broad class of generators $h(t)$.

##  III Recovery System and Theoretical Guarantees

###  III-A Recovery of SI Signal

In this section, we introduce the recovery system, its algorithm, and a
theorem that validates the uniqueness of recovery in SI space, as defined in
equation (3). We explain the system mechanics and mathematical principles,
followed by a discussion on method stability, highlighting its effectiveness
and areas for potential improvement.

Figure 3: Block diagram for modulo sampling for SI spaces as detailed in [38].
See theorem 1 for more details on the system and proof of recovery.

######  Theorem 1.

Let $\mathcal{H}$ be a SI space as in (2) and let $x(t)\in\mathcal{H}$ be the
input signal. LPF with a cutoff frequency of $\frac{\pi}{T}$ is applied to
yield $y(t)=\text{LPF}(x(t))$. Modulo operation is applied to $y(t)$,
resulting in $\mathcal{M}_{\lambda}y(t)$. The signal
$\mathcal{M}_{\lambda}y(t)$ is then sampled at rate $T_{s}<T$, producing the
samples $\mathcal{M}_{\lambda}y[nT_{s}]$.

The input signal $x(t)$ can be perfectly recovered from the samples
$\mathcal{M}_{\lambda}y[nT_{s}]$, provided that $H(\omega)$ is non-zero almost
everywhere within the frequency interval $[-\frac{\pi}{T},\frac{\pi}{T}]$.

###### Proof.

We can see that $y(t)$ is BL with Nyquist rate of $T$ and sampled at rate
$T_{s}<T$. Following the approach in [26, 27], this allows for the
reconstruction of the signal $y(t)$, meaning we obtain $Y(\omega)$ for
$\omega$ within $[-\frac{\pi}{T},\frac{\pi}{T}]$.

Since $y(t)=\text{LPF}(x(t))$, we know that $Y(\omega)=X(\omega)$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$, from the definition of LPF.
Following equation (3), we can apply a Fourier transform and using the
properties of DTFT and CTFT as shown in [1], obtaining $X(\omega)=A(e^{j\omega
T})H(\omega)$. So we conclude that $Y(\omega)=A(e^{j\omega T})H(\omega)$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$.

Given the theorem’s assumption, $A(e^{j\omega T})$ can be computed as
$\frac{Y(\omega)}{H(\omega)}$ for $\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$. A
discrete-time filter with a frequency response of $\frac{1}{H(\omega)}$ within
the frequency range $|\omega|\leq\frac{\pi}{T}$ is applied. This filter acts
as a correction filter, which ultimately leads to the extraction of $a[n]$ and
thereby the reconstruction of the original input signal $x(t)$ using (3). ∎

The crucial assumption in theorem 1 is that $H(\omega)\neq 0$ almost
everywhere within this frequency range, a condition met by most generators
used in practice. If $H(\omega)=0$ for any $\omega$ within the range
$[-\frac{\pi}{T},\frac{\pi}{T}]$, then the correction filter will enhance
noise at those frequencies.

To demonstrate stability within our recovery system, we explore the example of
a stable generator, specifically B-spline and Lorentzian. Define
$\beta^{(0)}(t)=\mathbf{1}_{[-0.5,0.5]}(t)$ as the $0$-th order B-spline [1].
The $n$-th order B-spline is recursively defined as
$\beta^{(n+1)}=\beta^{(n)}\ast\beta^{(0)}$. The CTFT of an $n$-th order spline
is $\text{sinc}^{n+1}(\omega)$. With no roots within
$[-\frac{\pi}{T},\frac{\pi}{T}]$, the correction filter becomes stable, with
frequency response of $\frac{1}{\text{sinc}^{n+1}(\omega)}$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$.

A second example for stable generator is Lorentzian function, defined as

$h_{\gamma}(t)=\frac{1}{\pi\gamma\Big{(}1+(\frac{x}{\gamma})^{2}\Big{)}}$ (4)

where $\gamma>0$ is a scale parameter. The CTFT is

$H_{\gamma}(\omega)=\exp(-\gamma|\omega|).$ (5)

The correction filter have an impulse response given by
$\frac{1}{H_{\gamma}(\omega)}=\exp(\gamma|\omega|)$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$, thereby robust to noise. Fig. 9
displays the Fourier transform of Lorentzian generators with $\gamma=0.25$ and
$\gamma=0.5$. Fig. 6 displays our simulation results for Lorentzian
generators, which yield a stable correction filter. Notably, the final
application of the correction filter does not significantly alter the mean
squared error (MSE). We will delve into more details in Section IV.

However, some generators demonstrate instability. This is particularly evident
with scaled first-order splines, as discussed in section IV. Fig. 8
illustrates the instability of the correction filter, attributed to roots at
the central frequency in $H(\omega)$. This root presence amplifies the
recovery error in the final step, as illustrated in Fig. 7, where the
correction filter exacerbates the error by approximately 30 dB across the
varying noise levels. This highlights the necessity for a more robust approach
to address this instability, which we introduce in the subsequent section.

###  III-B Improving Robustness

In this section, we address the enhancement of robustness in the presence of
noise. During the recovery process discussed previously, we measure
$\mathcal{M}_{\lambda}y[nT_{s}]$. When noise is present, we assume the
measurement is $\mathcal{M}_{\lambda}y[nT_{s}]+e[n]$, where $e[n]$ represents
the noise. We begin with a BL modulo recovery technique as detailed in [26,
27], resulting in $\mathcal{M}_{\lambda}y[nT_{s}]+e[n]\to\hat{y}[nT_{s}]$.
Naturally, there is an error $e_{\text{BL}}[n]=\hat{y}[nT_{s}]-y[nT_{s}]$.
Note that the modulo recovery algorithm is nonlinear, and a comprehensive
analysis of its robustness is detailed in the cited works.

The subsequent phase, absent in the BL signal modulo recovery, involves
division in the frequency domain by $H(\omega)$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$. Let $\hat{a}[n]$ be the estimated
coefficients. From the recovery process, we know from Theorem 1 that
$A(e^{j\omega})=\frac{Y(e^{j\omega})}{H(\omega)}$ and from the proposed
recovery process $\hat{A}(e^{j\omega})=\frac{\hat{Y}(e^{j\omega})}{H(\omega)}$
for $\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$. The final reconstruction error
is $e_{\text{coef}}[n]=\hat{a}[n]-a[n]$. We can analyze the error’s DTFT as
follows:

$\begin{split}E_{\text{coef}}(e^{j\omega})=\hat{A}(e^{j\omega})-A(e^{j\omega})=\frac{\hat{Y}(e^{j\omega})}{H(\omega)}-\frac{Y(e^{j\omega})}{H(\omega)}=\\\
=\frac{1}{H(\omega)}E_{\text{BL}}(e^{j\omega})\end{split}$ (6)

for $\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$. According to Parseval’s
theorem, the MSE in discrete time is equivalent to the MSE in the frequency
domain, so:

$\|e_{\text{coef}}\|_{2}^{2}=\Big{\|}\frac{1}{H(\omega)}E_{\text{BL}}(e^{j\omega})\mathbf{1}_{[-\frac{\pi}{T},\frac{\pi}{T}]}\Big{\|}_{2}^{2}.$
(7)

Many common generators distribute their energy unevenly across the frequency
interval $[-\frac{\pi}{T},\frac{\pi}{T}]$, leading to high MSE in equation
(7). The main concern arises when $H(\omega)$ has roots within
$[-\frac{\pi}{T},\frac{\pi}{T}]$, causing the error from the BL modulo
recovery $e_{\text{BL}}$ to amplify the coefficients error $e_{\text{coef}}$,
as shown in equation (7).

To enhance the method’s stability for such generators, we incorporate a mixer
using a $T$-periodic function $p(t)$ prior to the LPF. This transforms the
input signal into a BL signal with a Nyquist rate of $T$ in the analog domain,
resulting in $y(t)=\text{LPF}\Big{(}p(t)x(t)\Big{)}$. The subsequent steps
mirror the recovery process described in Theorem 1, including the application
of an analog modulo operator, sampling above the Nyquist rate, and using a
proposed modulo recovery algorithm [26, 27, 28]. The final step applies a
correction filter different from $\frac{1}{H(\omega)}$ that will be specified.
The mixer enables the use of a broader frequency range of $H(\omega)$ to
create a favorable correction filter that minimizes noise amplification as
seen in equation (7). This avoids singularities and results in a more constant
filter across frequencies. The mixer selection is tailored to the specific
needs of the SI space, sampling rate, and noise levels.

The added mixer operates on the analog signal, allowing for the use of a wider
frequency range of $H(\omega)$ without necessitating a higher sampling rate.
The only difference in the recovery process is the use of a different digital
correction filter, which is known derives from system’s setting. Figure 8
provides a visual representation, showing the correction filter for a scaled
B-spline with two singularities, indicating unstable reconstruction. It also
illustrates the correction filter with an appropriately selected mixer, which
avoids the singularities and noise amplification as seen in equation (7),
leading to improved reconstruction.

Figure 4: Block diagram for modulo sampling for SI spaces with improved
stability. See theorem 2 for more details on the system and proof of recovery.

Let $\mathcal{H}$ represent a SI space as defined in (2). Consider a
$T$-periodic signal $p(t)$. It can be expressed using its Fourier series
representation as $p(t)=\sum_{l\in\mathbb{Z}}{c_{l}e^{-j\frac{2\pi}{T}lt}}$,
where the coefficients $c_{l}$ are given by
$c_{l}=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}{p(t)e^{j\frac{2\pi}{T}lt}}dt$.
Denote

$R(\omega)=\sum_{l\in\mathbb{Z}}{c_{l}H\Big{(}\omega-\frac{2\pi
l}{T}\Big{)}}.$ (8)

######  Theorem 2.

Let $\mathcal{H}$ be a SI space as in (2) and the input signal
$x(t)\in\mathcal{H}$. A mixer with $T$-periodic signal $p(t)$ and LPF with
cut-off frequency $\frac{\pi}{T}$ are applied, resulting in
$y(t)=\text{LPF}\Big{(}p(t)x(t)\Big{)}$. The modulo operator is applied,
giving $\mathcal{M}_{\lambda}y(t)$. The sampling rate is $T_{s}<T$, resulting
in the samples $\mathcal{M}_{\lambda}y[nT_{s}]$.

The input signal $x(t)$ can be perfectly recovered from the samples
$\mathcal{M}_{\lambda}y[nT_{s}]$, assuming that $R(\omega)$, as defined in
(8), is non-zero almost everywhere for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$.

###### Proof.

From the properties of CTFT and DTFT [1] we can tell that

$\mathcal{F}\Big{\\{}p(t)x(t)\Big{\\}}=\sum_{l\in\mathbb{Z}}{c_{l}X\Big{(}\omega-\frac{2\pi
l}{T}\Big{)}}$ (9)

additionally from (3) and [1], we can see that $X(\omega)=A(e^{j\omega
T})H(\omega)$. We conclude that

$\begin{split}\mathcal{F}\Big{\\{}p(t)x(t)\Big{\\}}=\sum_{l\in\mathbb{Z}}{c_{l}A(e^{j\omega
T})H\Big{(}\omega-\frac{2\pi l}{T}\Big{)}}=\\\ A(e^{j\omega
T})\Bigg{(}\sum_{l\in\mathbb{Z}}{c_{l}H\Big{(}\omega-\frac{2\pi
l}{T}\Big{)}}\Bigg{)}=A(e^{j\omega T})R(\omega).\end{split}$ (10)

After applying LPF on $p(t)x(t)$, we get a BL signal $y(t)$. From the
definition of LPF we know that
$\mathcal{F}\Big{\\{}p(t)x(t)\Big{\\}}(\omega)=Y(\omega)$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$. Thus, with the use of (10)

$Y(\omega)=A(e^{j\omega
T})R(\omega),\quad\omega\in\Big{[}-\frac{\pi}{T},\frac{\pi}{T}\Big{]}.$ (11)

Note that $y(t)$ is BL with Nyquist rate of $T$, modulo operator is applied
and than sampled at rate $T_{s}<T$, resulting in the samples
$\mathcal{M}_{\lambda}y[nT_{s}]$. Following [26, 27], we can recover the
signal $y(y)$, meaning we get $Y(\omega)$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$.

Using equation (11), we can compute $A(e^{j\omega
T})=\frac{Y(\omega)}{R(\omega)}$. Note that we required in the theorem that
$R(\omega)\neq 0$ almost everywhere for $\omega$ in the interval
$[-\frac{\pi}{T},\frac{\pi}{T}]$. We apply a discrete-time filter with a
frequency response of $\frac{1}{R(\omega)}$ over the interval
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$. This operator serves as the
correction filter. Thereby obtaining $a[n]$ and consequently recovering the
input signal $x(t)$ using (3).

∎

Note that if we select $p(t)=1$ than $R(\omega)=H(\omega)$. From the
assumption in the theorem, we can see that the correction filter
$\frac{1}{R(\omega)}=\frac{1}{H(\omega)}$ in properly defined. Fig. 5
illustrates a noise-free recovery scenario, showcasing the complete sequence
of system operations. The input signal resides within the SI space formed by a
Lorentzian function, as defined in (4), having a scale parameter $\gamma=0.5$,
and the space’s shift is defined as $T=1$. The SI space is
$\mathcal{H}=\Big{\\{}\sum_{n\in\mathbb{Z}}{a[n]h_{0.5}(t-n)}|a\in\ell^{2}\Big{\\}}$.
Note that here we do not use a mixer, we set $\lambda=0.1$ and the
oversampling rate is set to $5$. In Fig. 5(a), the progression of the signal,
its state after undergoing LPF, and the outcome after applying the modulo
operation is depicted. Note that the original signal is not BL. The
application of a LPF, however, enables subsequent modulo recovery. Fig. 5(b)
displays both the true parameters $a[n]$ and the parameters recovered through
the process.

Figure 5: Signal recovery for a Lorentzian-based signal within the SI space.
(a) the initial input signal (blue)
$x(t)=\sum_{n\in\mathbb{Z}}{a[n]h_{0.5}(t-n)}\in\mathcal{H}$, followed by LPF,
$y(t)=\text{LPF}(x(t))$ (orange), and modulo operation
$\mathcal{M}_{\lambda}(y(t))$ (green). Subfigure (b) the parameter recovery
phase, meaning $a[n]$ and the reconstruction $\hat{a}[n]$.

##  IV Simulations and Noise Resilience

This section illustrates the modulo sampling process within SI spaces,
particularly focusing on the effects of noise. A key issue is the influence of
the LPF on system energy loss. We utilize a Lorentzian function generator
$h_{\gamma}(t)$, as defined in (4) with scales $\gamma=0.5$ and $\gamma=0.25$,
the Fourier transform is shown in (5). The SI space, as in (3) will use $T=1$.
The study compares two scenarios: one without a mixer and another with a non-
trivial mixer, which exploits the generator’s higher frequencies to stabilize
the correction filter.

Our simulations involved 500 distinct signals, each generated with 50
coefficients $a[n]$ drawn from a uniform distribution over $[-1,1]$. For the
scenario incorporating a mixer, the signals $x(t)$ were combined with $p(t)$,
followed by LPF application with a cutoff frequency $\frac{\pi}{T}$, resulting
in $y(t)=\text{LPF}(p(t)x(t))$. The mixer design was tailored to the SI space,
while in the non-mixer scenario (effectively a constant mixer), LPF was
directly applied. Post-normalization ensured the maximum absolute amplitude of
the samples was one, before applying the modulo operator with $\lambda=0.2$ to
obtain $\mathcal{M}_{\lambda}y(t)$. Sampling was then conducted at an
oversampling rate of 5, denoted as $T_{s}=\frac{T}{5}$, considering the
Nyquist rate of the post-LPF signal $y(t)$ is $T$. This setup provided samples
$\mathcal{M}_{\lambda}y[nT_{s}]$. We introduced noise at various Signal-to-
Noise Ratio (SNR) levels and utilized the $B^{2}R^{2}$ algorithm [26, 27] for
the modulo recovery of the BL signal $y[nT_{s}]$. To deduce the coefficients
$a[n]$, the correction filter from theorem 2 was applied, converting
$Y(\omega)\to\frac{Y(\omega)}{R(\omega)}=A(e^{j\omega T})$ for
$\omega\in[-\frac{\pi}{T},\frac{\pi}{T}]$, with $R(\omega)$ adjusted based on
the mixer selection.

The approximated coefficients were denoted as $\hat{a}[n]$. We evaluated MSE
metrics, calculating $\frac{\|a-\hat{a}\|_{2}^{2}}{\|a\|_{2}^{2}}$ for
coefficient recovery and
$\frac{\|y[nT_{s}]-\hat{y}[nT_{s}]\|_{2}^{2}}{\|y[nT_{s}]\|_{2}^{2}}$ for
modulo recovery, with $\hat{y}[nT_{s}]$ indicating the reconstructed modulo
samples. This method facilitated an independent assessment of modulo recovery
stability.

Fig. 6 graphically presents the results at various SNR levels, setting
$p(t)=1+200\cos(\frac{2\pi}{T}t)+200\cos(\frac{4\pi}{T}t)$. Analysis shows
similar levels of MSE in both the modulo recovery process and the final
recovery of coefficients, regardless of the mixer’s use. The consistency in
results is remarkable, considering the different approaches taken in
generating the signals and the use of distinct correction filters in each
instance, where the filters are expected to be stable. At lower SNR levels,
the results are closely matched, with the mixer scenario showing a marginal
improvement. However, at higher SNR levels, the performance without a mixer is
notably better.

The LPF application incurs energy loss, impacting the system’s performance
while not affecting the modulo recovery step directly. For $\gamma=0.5$, the
energy loss was about 2.2%, in contrast to $\gamma=0.25$, which experienced a
higher energy loss of 11%. Fig. 9 displays the CTFT of the generators and
their lower energy bounds.

Figure 6: MSE in dB for varying SNR values using a Lorentzian function
generator with scales $\gamma=0.5$ (top) and $\gamma=0.25$ (bottom). “BL
error” refers to the recovery discrepancy of $\hat{y}[nT_{s}]$, and
“coefficient error” relates to the precision in reconstructing $\hat{a}[n]$,
showcased for scenarios with and without mixer usage.

Additional simulations on an alternative generator aimed to analyze scenarios
with an unstable correction filter. Denote the first order B-spline as

$\beta^{1}(t)=\begin{cases}1-|t|,\quad-1\leq|t|\leq 1\\\
0,\quad\text{o.w}\end{cases}.$ (12)

Employing a scaled first-order B-spline by $2.5$ for a generator
$h(t)=\beta^{1}(\frac{t}{2.5})$ created a scenario where
$H(\omega)=\text{sinc}^{2}(2.5\omega)$, with singularities within
$[-\frac{\pi}{T},\frac{\pi}{T}]$. $\lambda$ remains he same but the mixer is
set to $p(t)=1+1000\cos(\frac{2\pi}{T}t)+1000\cos(\frac{4\pi}{T}t)$ with Fig.
7 showcasing the outcomes. Using this particular generator, a significant
decline in MSE is observed transitioning from the modulo recovery stage to the
final coefficient error, a consequence of the correction filter’s instability.
The introduction of a mixer addresses this problem by producing a considerably
more stable correction filter, thereby enhancing overall performance. Notably,
while the BL recovery shows improved results in the absence of a mixer at
higher SNR levels, the instability of the filter leads to a substantial
worsening of the final MSE compared to that achieved with the mixer’s
application.

Figure 7: MSE in decibels across SNR values for a scaled first-order B-spline
generator. Error metrics follow the same definition as in Fig. 6, with results
indicating the mixer’s contribution to improved recovery outcomes at lower SNR
levels.

Fig. 8 displays the $R(\omega)$ function for both scenarios in the
simulations, following $\ell_{2}$ normalization. In the no-mixer situation,
the presence of roots is evident, whereas they are absent in the mixer case.
Fig. 8 illustrates the correction filter, represented as
$\frac{1}{R(\omega)}$. In the no-mixer scenario, significant peaks are
observable, indicating that noise at those frequencies will be significantly
amplified. Conversely, the correction filter in the mixer-utilized scenario
exhibits much greater stability.

Figure 8: On the top, the $R(\omega)$ function as determined in the
simulations for both the mixer and non-mixer scenarios, post-$\ell_{2}$
normalization. On the bottom, the associated correction filter,
$\frac{1}{R(\omega)}$, note the peaks observed in the non-mixer case.

The primary concern in modulo recovery $\mathcal{M}_{\lambda}f(x)[nT_{s}]\to
f(x)[nT_{s}]$ varies with mixer use, as seen in the scaled-spline scenario.
The correction filter’s step $Y(\omega)\to\frac{Y(\omega)}{R(\omega)}$
emphasizes the importance of $h(t)$’s energy distribution across frequencies,
as regions with diminished energy magnify noise through the filter.

Figure 9: Fourier transforms of three generators within the interval
$[-\frac{\pi}{T},\frac{\pi}{T}]$. Scaled first-order spline (blue), Lorentzian
with $\gamma=0.5$ (orange), and with $\gamma=0.25$ (green).

The emphasis is on the impact of the LPF and mixer on the system’s performance
and the detailed examination of energy loss’s role in the overall efficiency
of the signal recovery process.

##  V Exploring Modulo Sampling in ECG Signal Analysis

ECG signals are vital for monitoring the heart’s electrical activity, offering
insights into cardiac health and identifying various heart conditions. These
signals capture the sequential heartbeats, each comprising distinct components
like the P wave, QRS complex, and T wave, which reflect different phases of
the heart cycle [39, 40].

ECG signals exhibit a wide dynamic range due to the amplitude variations in
heartbeat waveforms, stemming from physiological differences and the heart’s
activity level. The broad dynamic range poses challenges in signal acquisition
and processing, where capturing the full scale of signal variations without
distortion or loss of detail is critical. Modulo sampling presents a promising
approach to address these challenges by enabling the effective compression of
ECG signals’ dynamic range, facilitating high-fidelity capture and
reconstruction of the original waveform without the need for high-resolution
ADCs.

For this study, we utilized data from [41], a public dataset comprising ECG
recordings from 30 healthy adults under various measurement conditions. This
dataset, primarily focused on exploring radar-based vital sign detection,
includes high-quality ECG signals sampled at 1000 Hz. Given that the spectrum
of ECG signals is mostly contained up to 50 Hz, they can effectively be
sampled at lower rates, such as 250 Hz or 500 Hz, without significant loss of
information. However, the availability of ECG data sampled at 1000 Hz is
advantageous, providing a high sampling rate that ensures comprehensive
capture of the signal’s details, making it particularly suitable for detailed
analysis and processing tasks.

From this dataset, we extracted a single heartbeat from the first ECG channel
of the first participant to serve as our base pulse. This pulse represents a
typical cardiac cycle, embodying the characteristic waveform components
essential for accurate ECG analysis.

Mathematically, the SI space will be modeled using the base pulse represented
by $h(t)$, with a defined time shift $T>0$. This involves summing instances of
the base pulse, each shifted by integer multiples of $T$. We will specify
$T=0.05$ seconds, corresponding to a frequency $f=\frac{1}{T}=20$ Hz. The
coefficients will comprise a sparse array where selected indices are set to
one—representing the occurrences of heartbeats—and the remaining values are
zero, thus indicating the precise moments of heartbeats.

Applying the modulo operator to this ECG signal model allows us to compress
its dynamic range, facilitating the recording and analysis of ECG signals
under limited resolution conditions. Subsequently, we employ a recovery
algorithm designed to reconstruct the original ECG waveform from its modulo-
sampled version. Through this process, we aim to accurately recover precise
timing of each heartbeat, from the modulo samples. This methodology
underscores the potential of modulo sampling in enhancing ECG signal analysis,
promising improvements in both diagnostic tools and wearable health monitoring
devices.

We apply a LPF to the SI signal, setting the cutoff frequency at
$\frac{5\pi}{T}$. Following this, we implement the modulo operator, selecting
$\lambda$ to be $0.1$ times the maximum amplitude of the signal. The signal is
then sampled at an oversampling rate of 5. For the recovery of the modulo
operation, we employ the methodology developed by Bahandari et al. [15]. The
final step involves applying a correction filter, which is formulated based on
the base pulse and the time shift $T$. Refer to Fig. 10 for a visual
representation of the process.

Figure 10: The figure demonstrates the application of modulo sampling to ECG
signals in a shift-invariant context. (a) The base pulse used, shown at the
original sampling rate of $1000$ Hz. (b) The function $|H(\omega)|$, given
these parameters, demonstrates that the correction filter is equivalent to
$\frac{1}{H(\omega)}$. It is important to note that $|H(\omega)|$ is bounded
away from zero, ensuring some stability in the correction process. (c) A
segment of the input signal (in blue) alongside its modulo-folded version (in
orange), illustrating the folding process. (d) Compares the original signal
(in blue) with their reconstruction (in orange), underscoring the precise
recovery of coefficients and thereby the accurate reconstruction of the ECG
signal.

##  VI Conclusion

In this study, we introduced an innovative approach for the recovery of SI
signals through modulo sampling, highlighting the integration of a mixer and a
LPF to produce a BL signal. This BL signal is then accurately reconstructed
from its modulo samples if sampled above the Nyquist rate, building upon
findings from prior research. The original SI signal is effectively retrieved
from this BL counterpart. Our methodology enables recovery at any sampling
rate below $T$, which is the theoretical minimum rate for SI signals. Through
extensive simulations, we evaluated the system’s robustness against varying
noise levels, showcasing the significant advantage of incorporating a mixer.
This addition proves particularly beneficial for enhancing stability in cases
involving unstable generators, thereby expanding the practical applications
and effectiveness of modulo sampling in SI signal recovery.

## References

  * [1] Y. C. Eldar, _Sampling theory: Beyond bandlimited systems_. Cambridge University Press, 2015.
  * [2] M. Mishali and Y. C. Eldar, “Sub-Nyquist sampling,” _IEEE Signal Processing Magazine_ , vol. 28, no. 6, pp. 98–124, 2011.
  * [3] R. Marks, “Restoring lost samples from an oversampled band-limited signal,” _IEEE transactions on acoustics, speech, and signal processing_ , vol. 31, no. 3, pp. 752–755, 1983.
  * [4] R. Marks and D. Radbel, “Error of linear estimation of lost samples in an oversampled band-limited signal,” _IEEE transactions on acoustics, speech, and signal processing_ , vol. 32, no. 3, pp. 648–654, 1984.
  * [5] J. S. Abel, “Restoring a clipped signal,” in _Acoustics, Speech, and Signal Processing, IEEE International Conference On_. IEEE Computer Society, 1991, pp. 1745–1748.
  * [6] R. Rietman, J.-P. Linnartz, and E. P. de Vries, “Clip correction in wireless lan receivers,” in _2008 European Conference on Wireless Technology_. IEEE, 2008, pp. 174–177.
  * [7] J. P. A. Pérez, S. C. Pueyo, and B. C. López, _Automatic gain control_. Springer, 2011.
  * [8] D. Mercy, “A review of automatic gain control theory,” _Radio and Electronic Engineer_ , vol. 11, no. 51, pp. 579–590, 1981.
  * [9] H. Landau, “On the recovery of a band-limited signal, after instantaneous companding and subsequent band limiting,” _Bell System Technical Journal_ , vol. 39, no. 2, pp. 351–364, 1960.
  * [10] H. Landau and W. Miranker, “The recovery of distorted band-limited signals,” _Journal of Mathematical Analysis and Applications_ , vol. 2, no. 1, pp. 97–104, 1961.
  * [11] D. Park, J. Rhee, and Y. Joo, “A wide dynamic-range cmos image sensor using self-reset technique,” _IEEE Electron Device Letters_ , vol. 28, no. 10, pp. 890–892, 2007.
  * [12] K. Sasagawa, T. Yamaguchi, M. Haruta, Y. Sunaga, H. Takehara, H. Takehara, T. Noda, T. Tokuda, and J. Ohta, “An implantable cmos image sensor with self-reset pixels for functional brain imaging,” _IEEE Transactions on Electron Devices_ , vol. 63, no. 1, pp. 215–222, 2015.
  * [13] J. Yuan, H. Y. Chan, S. W. Fung, and B. Liu, “An activity-triggered 95.3 db dr $-$75.6 db thd cmos imaging sensor with digital calibration,” _IEEE journal of solid-state circuits_ , vol. 44, no. 10, pp. 2834–2843, 2009.
  * [14] A. Krishna, S. Rudresh, V. Shaw, H. R. Sabbella, C. S. Seelamantula, and C. S. Thakur, “Unlimited dynamic range analog-to-digital conversion,” _arXiv preprint arXiv:1911.09371_ , 2019.
  * [15] A. Bhandari, F. Krahmer, and R. Raskar, “On unlimited sampling and reconstruction,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 3827–3839, 2020.
  * [16] E. Romanov and O. Ordentlich, “Above the nyquist rate, modulo folding does not hurt,” _IEEE Signal Processing Letters_ , vol. 26, no. 8, pp. 1167–1171, 2019.
  * [17] G. Lu and H. Liu, “High dynamic range sensing using multi-channel modulo samplers,” in _2020 IEEE 11th Sensor Array and Multichannel Signal Processing Workshop (SAM)_. IEEE, 2020, pp. 1–5.
  * [18] A. Bhandari, F. Krahmer, and T. Poskitt, “Unlimited sampling from theory to practice: Fourier-prony recovery and prototype ADC,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 1131–1141, 2021.
  * [19] S. Rudresh, A. Adiga, B. A. Shenoy, and C. S. Seelamantula, “Wavelet-based reconstruction for unlimited sampling,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 4584–4588.
  * [20] A. Bhandari, F. Krahmer, and R. Raskar, “Unlimited sampling of sparse signals,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 4569–4573.
  * [21] O. Musa, P. Jung, and N. Goertz, “Generalized approximate message passing for unlimited sampling of sparse signals,” in _2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_. IEEE, 2018, pp. 336–340.
  * [22] D. Prasanna, C. Sriram, and C. R. Murthy, “On the identifiability of sparse vectors from modulo compressed sensing measurements,” _IEEE Signal Processing Letters_ , vol. 28, pp. 131–134, 2020.
  * [23] F. Ji, P. Pratibha, and W. P. Tay, “On folded graph signals,” in _2019 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_. IEEE, 2019, pp. 1–5.
  * [24] A. Bhandari and F. Krahmer, “HDR imaging from quantization noise,” in _2020 IEEE International Conference on Image Processing (ICIP)_ , Oct 2020, p. 101–105.
  * [25] R. Guo and A. Bhandari, “ITER-SIS: Robust unlimited sampling via iterative signal sieving,” in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Jun 2023, p. 1–5.
  * [26] E. Azar, S. Mulleti, and Y. C. Eldar, “Robust unlimited sampling beyond modulo,” _arXiv preprint arXiv:2206.14656_ , 2022.
  * [27] ——, “Residual recovery algorithm for modulo sampling,” in _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2022, p. 5722–5726.
  * [28] S. B. Shah, S. Mulleti, and Y. C. Eldar, “Lasso-based fast residual recovery for modulo sampling,” in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Jun 2023, p. 1–5.
  * [29] S. Mulleti and Y. C. Eldar, “Modulo sampling of FRI signals,” _arXiv preprint arXiv:2207.08774_ , 2022.
  * [30] S. Mulleti, E. Reznitskiy, S. Savariego, M. Namer, N. Glazer, and Y. C. Eldar, “A hardware prototype of wideband high-dynamic range analog-to-digital converter,” _IET Circuits, Devices & Systems_, vol. 17, no. 4, pp. 181–192, 2023.
  * [31] M. Unser, “Sampling-50 years after shannon,” _Proceedings of the IEEE_ , vol. 88, no. 4, pp. 569–587, 2000.
  * [32] Y. C. Eldar, “Compressed sensing of analog signals in shift-invariant spaces,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 8, pp. 2986–2997, 2009.
  * [33] F. Ji, W. P. Tay _et al._ , “Unlimited dynamic range signal recovery for folded graph signals,” _Signal Processing_ , vol. 198, p. 108574, 2022.
  * [34] C. Deboor, R. A. DeVore, and A. Ron, “The structure of finitely generated shift-invariant spaces in l2 (rd),” _Journal of Functional Analysis_ , vol. 119, no. 1, pp. 37–78, 1994.
  * [35] O. Christensen and Y. C. Eldar, “Oblique dual frames and shift-invariant spaces,” _Applied and Computational Harmonic Analysis_ , vol. 17, no. 1, pp. 48–68, 2004.
  * [36] A. Aldroubi and K. Gröchenig, “Nonuniform sampling and reconstruction in shift-invariant spaces,” _SIAM review_ , vol. 43, no. 4, pp. 585–620, 2001\.
  * [37] A. Bhandari and A. I. Zayed, “Shift-invariant and sampling spaces associated with the fractional fourier transform domain,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 4, pp. 1627–1637, 2011.
  * [38] Y. Kvich and Y. C. Eldar, “Modulo sampling and recovery in shift-invariant spaces,” in _ICASSP 2024 - 2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2024, pp. 11–15.
  * [39] E. Ajdaraga and M. Gusev, “Analysis of sampling frequency and resolution in ecg signals,” in _2017 25th Telecommunication Forum (TELFOR)_. Belgrade: IEEE, Nov. 2017, p. 1–4. [Online]. Available: http://ieeexplore.ieee.org/document/8249438/
  * [40] O. Kwon, J. Jeong, H. B. Kim, I. H. Kwon, S. Y. Park, J. E. Kim, and Y. Choi, “Electrocardiogram sampling frequency range acceptable for heart rate variability analysis,” _Healthcare Informatics Research_ , vol. 24, no. 3, p. 198–206, Jul. 2018.
  * [41] S. Schellenberger, K. Shi, T. Steigleder, A. Malessa, F. Michler, L. Hameyer, N. Neumann, F. Lurz, R. Weigel, C. Ostgathe, and A. Koelpin, “A dataset of clinically recorded radar vital signs with synchronised reference sensor signals,” _Scientific Data_ , vol. 7, no. 11, p. 291, Sep. 2020\.


# DIFAI: Diverse Facial Inpainting using StyleGAN Inversion

###### Abstract

Image inpainting is an old problem in computer vision that restores occluded
regions and completes damaged images. In the case of facial image inpainting,
most of the methods generate only one result for each masked image, even
though there are other reasonable possibilities. To prevent any potential
biases and unnatural constraints stemming from generating only one image, we
propose a novel framework for diverse facial inpainting exploiting the
embedding space of StyleGAN. Our framework employs pSp encoder and SeFa
algorithm to identify semantic components of the StyleGAN embeddings and feed
them into our proposed SPARN decoder that adopts region normalization for
plausible inpainting. We demonstrate that our proposed method outperforms
several state-of-the-art methods.

Index Terms—  Facial Image Inpainting, Pluralistic Image Inpainting, StyleGAN
Inversion

##  1 Introduction

Image inpainting is one of the tasks in computer vision that removes undesired
objects or restores occluded regions. Since it is a well-known problem in
computer vision, numerous approaches have been proposed in the past. Among
them, traditional approaches [1, 2] propagate small patches from the
background area to the missing regions using similarity. However, unlike
natural or landscape image inpainting, facial images have unique parts such as
nose or mouth, so these methods could not be used.

With the recent development of generative adversarial learning [3], GAN-based
image inpainting is capable of synthesizing plausible results even when novel
objects are present. GLCIC [4] utilized dilated convolutional layers and two
auxiliary discriminators in order to complete the missing parts. CA [5] was
the first attention-based method that uses coarse-to-fine networks for
context-aware inpainting. More lately, LBAM [6] proposed learnable
bidirectional attention maps that enable more realistic inpainting for
irregular masks.

For inpainting an image with backgrounds and objects, there are a number of
plausible ways to fill in the missing part. Thus, the methods capable of only
suggesting one of these many possibilities would have some built-in biases and
constraints that may pose limitations on the utilities of such methods. For
instance, it would be highly desirable in GAN-based data augmentation to have
a method capable of suggesting a variety of possibilities in a missing region
of an image.

To overcome this limitation, PIC [7] used a short+long term attention layer to
synthesize pluralistic inpainting results. Furthermore, this method proposed a
probabilistic principled framework comprised of two parallel paths called
generative path and reconstructive path. PD-GAN [8] provides diverse
inpainting results using the proposed SPDNorm Resblock. PD-GAN also adds
perceptual diversity loss for various inpainting possibilities in the training
process. Contrary to the well-known diversity loss [9] perceptual diversity
loss is calculated on the perceptual space for keeping the context unchanged
and training stable.

On a related front, numerous researches in GAN have been conducted to
synthesize diverse high-resolution images. StyleGAN-based methods [10, 11, 12]
have been highly successful in generating some astonishing images by
controlling the latent space. Pixel2style2pixel (pSp) [13] proposed an encoder
network that directly generates a sequence of style vectors which are input
into a StyleGAN decoder, forming the extended latent space. SeFa [14] or
GANSpace [15] uncover relevant directions in the latent space of pre-trained
StyleGAN that affect the semantic properties of the decoded image in an
unsupervised manner. The aforementioned methods are capable of generating a
variety of images in a range of details from coarse features such as shapes or
poses to fine properties such as lighting, background attributes, or feature
variations.

Fig. 1: Summary of the proposed framework. $I_{style^{\prime}}$ and $M$ are
the input of SPARN decoder. In the SPARN decoder, each region normalization
layer uses the $I_{style^{\prime}}$ and $M$ to modulate the layer activations.

Given this, we propose a novel framework for diverse facial inpainting based
on controlling StyleGAN’s latent space for generating a set of plausible
inpainted regions while maintaining the remaining regions. Our approach
requires only an image with a masked region as input. Our framework first
coarsely completes the input (masked) image from the pre-trained inpainting
network so that pSp encoder will extract style vectors in the latent space.
Afterward, we manipulate the latent space in meaningful directions to
transform the semantic attributes of the decoded images. By feeding the
manipulated latent space into the StyleGAN decoder, we could generate images
with transformed facial shapes or attributes. Additionally, we feed decoded
images as a condition into our proposed spatially adaptive region
normalization (SPARN) decoder.

Our proposed decoder adopts region normalization [16] in each layer to allow
synthesizing realistic inpainting results. Thus, the proposed generator can be
trained to perform more diverse image inpainting using StyleGAN without any
prior condition. Based on our experiments of the publicly released dataset
CelebA-HQ, we demonstrate that the proposed approach delivers superior
performance compared to various state-of-the-art approaches specialized in
inpainting tasks.

##  2 Proposed Method

###  2.1 Proposed Framework

This section introduces our proposed facial image inpainting framework. As
shown in Fig. 1, our framework consists of four parts: a pre-trained
inpainting network, a pSp encoder, a StyleGAN decoder, and the proposed
generator. We first apply a customized MLGN [17] model for coarse inpainting.
Our customized MLGN differs from the original model by an adjustment made to
the lambda parameter for synthesizing blurry results. The blurry results would
promote more diverse embeddings by pSp encoder, thereby allowing StyleGAN to
generate more diverse image inpainting. We generate ground truth and masked
image pairs as

$I_{masked}=I_{gt}\odot M,$ (1)

where $I_{gt}$ is the ground truth image, $M$ is the mask applied to erase
portions of the ground truth image, and $I_{masked}$ is the masked image. We
input $I_{masked}$ into a pre-trained customized MLGN as

$I_{coarse}=\text{cMLGN}(I_{masked}).$ (2)

|  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---
Input |  $\text{PIC}_{1}$[7] | $\text{PIC}_{2}$ | $\text{PIC}_{3}$ | $\text{PIC}_{4}$ | $\text{Ours}_{1}$ | $\text{Ours}_{2}$ | $\text{Ours}_{3}$ | $\text{Ours}_{4}$

Fig. 2: Pluralistic qualitative comparison with PIC [7] and Ours. |  |  |  |  |
---|---|---|---|---|---
Input | LBAM[6] | EC[18] | MLGN[17] | Ours |

Fig. 3: Qualitative comparison with LBAM [6], EC [18], MLGN [17], and Ours.
Fig. 4: Illustration of continuous inpainting results using SeFa that adjusts
delta gradually.

Most image inpainting methods generate only one result for each masked image,
even though there are many other possibilities. As such, there are always
possibilities of unrealistic biases and constraints due to the network being
forced to produce only one of many plausible results. To prevent such
artificial biases, we use StyleGAN-based [11] image augmentation that is
capable of synthesizing a variety of images that has a similar structure as
the ground truth but with changed facial attributes. $I_{coarse}$ is applied
to pSp [13] encoder which maps the embedding vector $w$ in a latent space
$W+$. The extracted $w$ is then decoded to produce an initial set of diverse
images using StyleGAN. Afterward, we use the SeFa [14] algorithm that performs
eigen-decomposition of StyleGAN’s weight matrix to discover principal
components that span dominant changes in the decoded images. We feed the
proposed generator with the embeddings $w_{\delta_{i}}$ perturbed by
$\delta_{i}$ in a number of principal directions during the training process
for synthesizing multiple images.

$\displaystyle\begin{split}I_{style}&=\text{StyleGAN}(\text{pSp}(I_{corase})),\\\
I_{style+}&=\\{\text{StyleGAN}(w_{\delta_{1}}),...,\text{StyleGAN}(w_{\delta_{\alpha}})\\},\end{split}$
(3)

where $\alpha$ is the number of $I_{style+}$ for training. Our generator
$G(\cdot)$ is comprised of a SPADE [19] encoder and the proposed SPARN
decoder. The proposed SPARN decoder may look similar to the existing SPADE
decoder, but it can maintain consistency in the masked and unmasked regions by
using region normalization [16] for image inpainting. As an input condition of
SPARN decoder, we define $I_{style^{\prime}}$ as,

$I_{style^{\prime}}=I_{masked}+M_{r}\odot\\{I_{style},I_{style+}\\},$ (4)

where, $M_{r}$ is the reversed mask. Additionally, we input $I_{masked}$ into
the SPADE encoder to ensure that features present in the masked image are
maintained in the output images. Our SPARN decoder consists of SPARN residual
block and following upsampling layers. Since each residual block runs at a
different scale, we downsample the input $M$ and $I_{style^{\prime}}$ to match
the spatial resolution. Thereby, we could perform more diverse facial image
inpainting using as conditions various images that transformed several facial
attribute detail.

$I_{out}=G(I_{masked},I_{style^{\prime}},M)$ (5)

Furthermore, we encode $I_{out}$ and $I_{gt}$ into the discriminator for
calculating adversarial loss. We use Spectral Normalization [20, 21] in the
discriminator since it is faster and more stable than other normalizations by
a simple formulation.

###  2.2 Loss Function

In order to synthesize plausible and realistic image inpainting, we define our
loss function in two parts: inpainting loss and adversarial loss. Our proposed
inpainting loss is composed of four components: reconstruction loss, VGG style
loss, perceptual loss, and MS-SSIM loss. Reconstruction loss completes
occluded regions using $l1$-norm error. By comparing the generated image to
the ground truth, we calculate the hole region loss and valid region loss,
respectively. Additionally, we define perceptual loss and VGG style loss with
the VGG-19 network [22] pre-trained on ImageNet. As the name indicates,
perceptual loss [23] measures the feature map distance between the generated
image and the ground truth image. Our perceptual loss looks similar to the
aforementioned perceptual loss, but we further measure the distance between
$I_{out}$ with $I_{style}$. We adopt $M_{r}$ to reflect $I_{style}$ more
plausible in erased regions. We defined our perceptual loss as,

CC L_per=∑_i∑_j=1^α——F_i(I_out_j)⋅M_r-F_i(I_style_j)⋅M_r——_1
+∑_i——F_i(I_out)⋅M-F_i(I_gt)⋅M——_1, where $F_{i}$ denotes the feature maps of
the $i\,^{\prime}th$ layer of a VGG-19 network. We use VGG style loss, as
defined by [24], which alleviates “checkerboard” artifacts caused by
upsampling convolution layers. Our VGG style loss also compares $I_{out}$ with
$I_{style}$, using $M_{r}$.

CC L_style=∑_k∑_j=1^α——G^F_k(I_out_j)⋅M_r-G^F_k(I_style_j)⋅M_r——_1
+∑_k—G^F_k(I_out)⋅M-G^F_k(I_gt)⋅M——_1, where $G^{F}_{k}$ is a gram matrix
consisting of feature maps $F_{k}$. Additionally, we customize another loss
function by utilizing MS-SSIM [25, 17], which is one of the image quality
comparison approaches.

$L_{\text{MS-SSIM}}=1-\frac{1}{N}\sum_{i=1}^{N}\text{MS-SSIM}_{n}$ (6)

We calculate adversarial loss using WGAN-GP which optimizes the Wasserstein
distance. We define adversarial loss $L_{G}$ and $L_{D}$ as, CCL
L_G=E_I_masked[D(G(I_masked, I_style^′, M))],
L_D=E_I_gt[D(I_gt)]-E_I_out[D(I_out)]
-λ_gpE_^I[(—— ∇_^ID(^I)——_2-1)^2].

Our overall loss denoted as, CC L_all=λ_advL_adv+λ_ssimL_MS-SSIM+λ_styL_style
+L_per+λ_holeL_hole+λ_validL_valid, where, $\lambda$ are hyper-parameters that
control the terms’ relative importance.

| Mask | Ours | LBAM | EC | MLGN
---|---|---|---|---|---
SSIM | Quickdraw | 0.833 | 0.821 | 0.817 | 0.832
10-20% | 0.837 | 0.814 | 0.827 | 0.839
20-30% | 0.777 | 0.744 | 0.761 | 0.777
30-40% | 0.709 | 0.667 | 0.681 | 0.706
40-50% | 0.633 | 0.583 | 0.595 | 0.624
LPIPS | Quickdraw | 0.049 | 0.047 | 0.047 | 0.063
10-20% | 0.052 | 0.057 | 0.051 | 0.065
20-30% | 0.084 | 0.087 | 0.082 | 0.103
30-40% | 0.124 | 0.126 | 0.128 | 0.148
40-50% | 0.170 | 0.172 | 0.183 | 0.201
FID | Quickdraw | 25.95 | 25.79 | 27.49 | 28.45
10-20% | 24.74 | 27.63 | 25.65 | 26.73
20-30% | 34.55 | 36.51 | 34.80 | 38.34
30-40% | 46.87 | 48.47 | 47.14 | 52.54
40-50% | 64.86 | 64.40 | 63.75 | 73.07
Table 1: Quantitative comparison on CelebA-HQ. In each row, the best results
are shown in bold text.

##  3 Experiments

###  3.1 Implement Details

For the implementation, we used the Pytorch library. Our hyper-parameters
$\lambda_{adv}$, $\lambda_{ssim}$, $\lambda_{sty}$, $\lambda_{hole}$, and
$\lambda_{valid}$ are set to 0.5, 120, 3, and 0.5 respectively. In this paper,
we evaluate all the models using CelebA-HQ dataset and split them into two
groups: 28,000 selected for training and 2000 for testing. We used 256×256
images with irregular holes to train and evaluate the proposed methods. In
addition, we combine Quickdraw irregular mask dataset [26] with 85×85 square
holes in random positions to create more irregular holes. By combining square
holes with the Quickdraw dataset, the model becomes more robust to irregular
holes.

| Mask | PIC | Ours
---|---|---|---
Diversity (LPIPS) $\uparrow$ | 20-30% | 0.0714 | 0.0849
30-40% | 0.1079 | 0.1313
Table 2: Quantitative comparison of diversity on CelebA-HQ. For the diversity
comparison, higher LPIPS is better.

###  3.2 Qualitative Comparisons

First, we compare the image inpainting quality of our baseline against three
state-of-the-art methods. Fig. 3 describes the images generated by the
proposed method and those generated by the other methods. Our model is
superior to all the others in the aspect of image quality and plausibility.
Fig. 2 compares diverse images generated by PIC [7] and ours. Compares to the
PIC, our method accomplishes more diverse and pluralistic instances.

|  | SPADE
---
w/o RN | Ours
SSIM | 0.816 | 0.817 | 0.833
LPIPS | 0.056 | 0.057 | 0.049
FID | 27.96 | 28.45 | 25.95
Table 3: Quantitative comparison of ablation study on CelebA-HQ. We choose the
Quickdraw irregular mask here.

###  3.3 Quantitative Comparisons

We implement a quantitative comparison of image inpainting to three existing
methods and our own, using different types and sizes of masks. As shown in
Table. 1, our method outperforms three metrics SSIM, LPIPS [27], and FID [28]
to existing methods that specialize in only image inpainting tasks. Table. 2
shows that our method achieves a relatively higher diversity score than
another method. The diversity score is calculated between 4K pairs synthesized
from a sampling of 1K images. Overall, we feed the $w$ into StyleGAN to
calculate Table. 1 and feed the $w_{\delta}$ to calculate the diversity score.

To justify the effectiveness of the proposed SPARN decoder, we conduct the
ablation study as follows: 1) Using SPADE [19] decoder; 2) replacing all the
region normalization [16] with batch normalization (w/o RN). As shown in
Table. 3, each of the proposed sub-modules performs a very important role in
the overall architecture.

##  4 Conclusion

We propose a novel method of generating diverse facial inpainted images based
on manipulating StyleGAN embedding space. To properly discover meaningful
direction and the associated variations for diverse facial inpaintings, we
utilize pSp encoder and SeFa algorithm. These embedded vectors and the
variations are then fed into our proposed SPARN decoder as conditions for
diverse inpainting. From the proposed framework, we demonstrated that our
method synthesizes plausible diverse images from a single masked input while
maintaining high inpainting quality.

## References

  * [1] Alexei A Efros and Thomas K Leung,  “Texture synthesis by non-parametric sampling,”  in ICCV, 1999.
  * [2] Connelly Barnes, Eli Shechtman, Adam Finkelstein, and Dan B Goldman,  “Patchmatch: A randomized correspondence algorithm for structural image editing,”  ToG, 2009.
  * [3] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio,  “Generative adversarial nets,”  in NeurIPS, 2014.
  * [4] Satoshi Iizuka, Edgar Simo-Serra, and Hiroshi Ishikawa,  “Globally and locally consistent image completion,”  ToG, 2017.
  * [5] Jiahui Yu, Zhe Lin, Jimei Yang, Xiaohui Shen, Xin Lu, and Thomas S Huang,  “Generative image inpainting with contextual attention,”  in CVPR, 2018.
  * [6] Chaohao Xie, Shaohui Liu, Chao Li, Ming-Ming Cheng, Wangmeng Zuo, Xiao Liu, Shilei Wen, and Errui Ding,  “Image inpainting with learnable bidirectional attention maps,”  in ICCV, 2019.
  * [7] Chuanxia Zheng, Tat-Jen Cham, and Jianfei Cai,  “Pluralistic image completion,”  in CVPR, 2019.
  * [8] Hongyu Liu, Ziyu Wan, Wei Huang, Yibing Song, Xintong Han, and Jing Liao,  “Pd-gan: Probabilistic diverse gan for image inpainting,”  in CVPR, 2021.
  * [9] Qi Mao, Hsin-Ying Lee, Hung-Yu Tseng, Siwei Ma, and Ming-Hsuan Yang,  “Mode seeking generative adversarial networks for diverse image synthesis,”  in CVPR, 2019.
  * [10] Tero Karras, Samuli Laine, and Timo Aila,  “A style-based generator architecture for generative adversarial networks,”  in CVPR, 2019.
  * [11] Tero Karras, Samuli Laine, Miika Aittala, Janne Hellsten, Jaakko Lehtinen, and Timo Aila,  “Analyzing and improving the image quality of stylegan,”  in CVPR, 2020.
  * [12] Tero Karras, Miika Aittala, Samuli Laine, Erik Härkönen, Janne Hellsten, Jaakko Lehtinen, and Timo Aila,  “Alias-free generative adversarial networks,”  in NeurIPS, 2021.
  * [13] Elad Richardson, Yuval Alaluf, Or Patashnik, Yotam Nitzan, Yaniv Azar, Stav Shapiro, and Daniel Cohen-Or,  “Encoding in style: a stylegan encoder for image-to-image translation,”  in CVPR, 2021.
  * [14] Yujun Shen and Bolei Zhou,  “Closed-form factorization of latent semantics in gans,”  in CVPR, 2021.
  * [15] Erik Härkönen, Aaron Hertzmann, Jaakko Lehtinen, and Sylvain Paris,  “Ganspace: Discovering interpretable gan controls,”  Arxiv, 2020.
  * [16] Tao Yu, Zongyu Guo, Xin Jin, Shilin Wu, Zhibo Chen, Weiping Li, Zhizheng Zhang, and Sen Liu,  “Region normalization for image inpainting,”  in AAAI, 2020.
  * [17] Jie Liu and Cheolkon Jung,  “Facial image inpainting using multi-level generative network,”  in ICME, 2019.
  * [18] Kamyar Nazeri, Eric Ng, Tony Joseph, Faisal Qureshi, and Mehran Ebrahimi,  “Edgeconnect: Structure guided image inpainting using edge prediction,”  in ICCVW, 2019.
  * [19] Taesung Park, Ming-Yu Liu, Ting-Chun Wang, and Jun-Yan Zhu,  “Semantic image synthesis with spatially-adaptive normalization,”  in CVPR, 2019.
  * [20] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida,  “Spectral normalization for generative adversarial networks,”  Arxiv, 2018.
  * [21] Jiahui Yu, Zhe Lin, Jimei Yang, Xiaohui Shen, Xin Lu, and Thomas S Huang,  “Free-form image inpainting with gated convolution,”  in ICCV, 2019.
  * [22] Karen Simonyan and Andrew Zisserman,  “Very deep convolutional networks for large-scale image recognition,”  Arxiv, 2014.
  * [23] Justin Johnson, Alexandre Alahi, and Li Fei-Fei,  “Perceptual losses for real-time style transfer and super-resolution,”  in ECCV, 2016.
  * [24] Mehdi SM Sajjadi, Bernhard Scholkopf, and Michael Hirsch,  “Enhancenet: Single image super-resolution through automated texture synthesis,”  in ICCV, 2017.
  * [25] Zhou Wang, Eero P Simoncelli, and Alan C Bovik,  “Multiscale structural similarity for image quality assessment,”  in ACSSC, 2003.
  * [26] Karim Iskakov,  “Semi-parametric image inpainting,”  Arxiv, 2018.
  * [27] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang,  “The unreasonable effectiveness of deep features as a perceptual metric,”  in CVPR, 2018.
  * [28] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter,  “Gans trained by a two time-scale update rule converge to a local nash equilibrium,”  in NeurIPS, 2017.