# Covariance Recovery for One-Bit Sampled Data With Time-Varying Sampling
Thresholds—
Part ii@: Non-Stationary Signals

Arian Eamaz, , Farhang Yeganegi, and
Mojtaba Soltanalian This work was supported in part by the National Science
Foundation Grants CCF-1704401 and ECCS-1809225. Parts of this work were
presented at the International Conference on Acoustics, Speech and Signal
Processing (ICASSP) 2021, held in Toronto, Canada [1]. The first two authors
contributed equally to this work.A. Eamaz, F. Yeganegi and M. Soltanalian are
with the Dept. of Electrical and Computer Engineering, University of Illinois
Chicago, Chicago, USA (e-mails<EMAIL_ADDRESS>fyegan2@uic.edu).

###### Abstract

The recovery of the input signal covariance values from its one-bit sampled
counterpart has been deemed a challenging task in the literature. To deal with
its difficulties, some assumptions are typically made to find a relation
between the input covariance matrix and the autocorrelation values of the one-
bit sampled data. This includes the arcsine law and the modified arcsine law
that were discussed in Part i@ of this work [2]. We showed that by
facilitating the deployment of time-varying thresholds, the modified arcsine
law has a promising performance in covariance recovery. However, the modified
arcsine law also assumes input signals are stationary, which is typically a
simplifying assumption for real-world applications. In fact, in many signal
processing applications, the input signals are readily known to be non-
stationary with a non-Toeplitz covariance matrix. In this paper, we propose an
approach to extending the arcsine law to the case where one-bit ADCs apply
time-varying thresholds while dealing with input signals that originate from a
non-stationary process. In particular, the recovery methods are shown to
accurately recover the time-varying variance and autocorrelation values.
Furthermore, we extend the formulation of the Bussgang law to the case where
non-stationary input signals are considered.

###### Index Terms:

Arcsine law, Bussgang law, covariance matrix, one-bit quantization, modified
arcsine law, non-stationary signals, time-varying thresholds, time-varying
signal statistics.

##  I Introduction

Covariance matrix recovery plays an important role in statistical signal
processing applications such as directions of arrival (DOA) estimation, radar
waveform design, target parameter estimation, communication channel
estimation, and adaptive radar detection [3, 4, 5, 6, 7, 8, 9, 10]. When
digital signal processing is concerned, using one-bit quantization and
digitization, in which the input signals are compared with given threshold
levels, allows for sampling at a very high rate and with lower energy
consumption [11, 12, 13, 14]. As a result of employing one-bit sampling,
however, we can only use the sign data as partial available information to
recover the signal covariance, and second order statistics in general, making
it more challenging. In [15, 16, 17, 18], the authors have considered the
input signal as a stationary zero-mean Gaussian process, and with this
assumption, the input covariance is recovered by taking advantage of the
_arcsine law_ which connects the covariance of an unquantized signal with that
of its quantized counterpart [19, 13]. In [20], the relationship between the
cross-correlation matrix of the input signal and the one-bit output data is
characterized as the _Bussgang law_ for the stationary zero-mean Gaussian
signals. Note that the sampling threshold levels are considered to be zero in
these research efforts. The zero threshold values can give rise to some
difficulties in signal amplitude recovery and a considerable portion of signal
information may be lost. As a natural alternative, time-varying thresholds are
utilized in recent works which can lead to enhancements in recovery
performance [1, 21, 22, 23, 24, 25, 26].

Owing to the successful performance of time-varying sampling thresholds for
signals amplitude recovery, such time-varying thresholds were considered for
the covariance recovery problem [1], and more extensively in Part i@ of this
work [2], exhibiting a significantly improved performance in the estimation of
signal autocorrelation values via a _modified arcsine law_. Moreover, taking
time-varying thresholds into consideration for cross-correlation matrix
recovery, promising results were demonstrated with a _modified Bussgang law_
in Part i@ of this work [2].

A critical restriction of the arcsine law, as well as the modified arcsine
law, lies in the necessary assumption of a stationary input signal[15, 1, 2].
In real-world communication and digital signal processing applications,
however, input signals are non-stationary in general and have time-varying
variances [27, 28, 29]. In such scenarios, covarince recovery is an even more
prominent tool in the analysis of non-stationary processes and systems, and
can provide useful insights into their innate dynamics [30, 31, 32].
Nevertheless, in a non-stationary environment, the expected accuracy of
covariance recovery is typically diminished.

In this paper, we present an approach to extend our modified arcsine law for
time-varying sampling thresholds, discussed in Part i@ of this work [2], to
recover signal covariance matrices with an arbitrary non-Toeplitz structure.
Moreover, a Bussgang law with time-varying thresholds is established for the
non-stationary scenario.

###  I-A Contribution of the Paper

We study the covariance recovery for a non-stationary input signals in one-bit
quantization systems deploying time-varying thresholds. In particular, we
formulate an integral-based relation between the autocorrelation function of
the one-bit sampled data and the generic covariance matrix entries of the
input signal. Moreover, a closed-form formulation for the mean of the input
signal is obtained and the utilized to recover the time-varying signal
variances. It is demonstrated that to recover the autocorrelation values, we
should evaluate the obtained integral which appears to be intractable
analytically. To approach this problem, we first employ a one-point piece-wise
Padé approximation (PA) to recast the integrands as rational expressions which
are readily integrable. Next, we formulate an estimation criterion to recover
the desired input autocorrelation values. The accuracy of the PA is also
investigated. In the next step, two well-known numerical integration
techniques are employed to estimate the input autocorrelation values; namely,
the Gauss-Legendre quadrature and the Monte-Carlo integration techniques.
Interestingly, the proposed estimation criteria for these approaches take
convex form, facilitating an accelerated recovery. Lastly, a modified Bussgang
law for non-stationary input signals is presented. By using the modified
Bussgang law, the matrix elements associated with the cross-correlation
between the input signal and the one-bit sampled data can be recovered.
Several numerical results are presented to illustrate the effectiveness of the
proposed methodologies.

###  I-B Organization of the Paper

Section ii@ is dedicated to formulating the autocorrelation function of the
one-bit sampled data with time-varying thresholds in the case of non-
stationary inputs. In Section iii@, the time-varying variances are recovered
by using the proposed formula for the mean of the one-bit sampled data.
Sections iv@ presents our Padé Approximation (PA) to recover the input signal
autocorrelation sequence. Subsequently, Sections v@ and vi@ discuss two
widely-known numerical integration techniques, i.e. the Gauss-Legendre
quadrature and the Monte-Carlo integration methods, applied to our arbitrary
non-Toeplitz covariance matrix recovery problem. Section vii@ is where the
various methods proposed for covariance recovery are compared. A proper
thresholding for covariance recovery through the estimation of the threshold
mean is discussed in Section<EMAIL_ADDRESS>The modified Bussgang law for time-varying
thresholds in the case of non-stationary signals is presented in Section<EMAIL_ADDRESS>Finally, Section x@ concludes the paper.

_Notation:_ We use bold lowercase letters for vectors, and bold uppercase
letters for matrices and uppercase letters for matrix entries. For instance,
$\bm{R}_{\mathbf{x}}$ and $R_{\mathbf{x}}(i,j)$ denote the autocorrelation
matrix and the $ij$-th element of the autocorrelation matrix of the vector
$\mathbf{x}$, respectively. $(\cdot)^{\top}$ and $(\cdot)^{\mathrm{H}}$ denote
the vector/matrix transpose, and the Hermitian transpose, respectively.
$\mathbb{E}\left\\{.\right\\}$ denotes the expected value. The Frobenius norm
of a matrix $\bm{B}\in\mathbb{C}^{M\times N}$ is defined as
$\|\bm{B}\|_{\mathrm{F}}=\sqrt{\sum^{M}_{r=1}\sum^{N}_{s=1}\left|B(r,s)\right|^{2}}$
where $\\{B(r,s)\\}$ are entries of $\bm{B}$. For an event $\mathcal{E}$,
$\mathbb{I}_{(\mathcal{E})}$ is the indicator function for that event; i.e.
$\mathbb{I}_{(\mathcal{E})}$ is $1$ if $\mathcal{E}$ occurs and $0$ otherwise.
The $Q$-function is defined as

$\displaystyle Q(x)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}\exp\left(-\frac{z^{2}}{2}\right)\,dz.$
(1)

Further, $Q^{-1}$ is an inverse $Q$-function. The error function (erf) is
defined as
$\operatorname{erf}x=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-z^{2}}\,dz$. The
incomplete Gamma function is defined as

$\Gamma(s,x)=\int_{x}^{\infty}z^{s-1}e^{-z}\,dz.$ (2)

Finally, the cumulative distribution function (CDF) of a zero mean Gaussian
process $\bm{z}\sim\mathcal{N}(0,\zeta)$ is defined as

$\Psi(\bm{z})=\frac{1}{\sqrt{2\pi}}\int^{z}_{-\infty}e^{-\frac{t^{2}}{2\zeta^{2}}}\,dt.$
(3)

##  II Covariance Recovery In Non-Stationary Scenario

We assume that the input signal is a zero-mean non-stationary Gaussian process
$\mathbf{x}\sim\mathcal{N}\left(\mathbf{0},\bm{R}_{\mathbf{x}}\right)$, where
$\bm{R}_{\mathbf{x}}$ is the non-Toeplitz covariance matrix of $\mathbf{x}$
with the time-varying diagonal elements. Specifically, the input signal is
supposed to be non wide sense stationary (WSS) or weak stationary (for
abbreviation we use ”non-stationary” instead of ”non-WSS”). A signal is wide
sense stationary (weak stationary) if and only if:

  * a.

The signal mean is constant over time.

  * b.

$R_{\mathbf{x}}(i,j)=R_{\mathbf{x}}(l)$. This means that we have a constant
variance over time; i.e. $R_{\mathbf{x}}(i,i)=R_{\mathbf{x}}(0)$.

  * c.

$\mathbb{E}\left\\{x^{2}_{i}\right\\}<\infty$ or generally
$\mathbb{E}\left\\{x^{2n}_{i}\right\\}<\infty$ [33, 34, 35].

In our case, the input signal can satisfy (a) and (c), however, the
autocorrelation function does not rely on lag $R_{\mathbf{x}}(i,j)\neq
R_{\mathbf{x}}(l)$, therefore,
$\bm{R}_{\mathbf{x}}=\left[R_{\mathbf{x}}(i,j)\right]$ is a non-Toeplitz
covariance matrix with distinct diagonal elements (time-varying variance).
Consequently, our signal is non-stationary or marginal heteroskedastic [33,
35]. A simple but famous example for such signals are those originating in
_Wiener processes_ or _Brownian motion_
$\eta_{t}\sim\mathcal{N}\left(0,t\right)$. The input signal
$\mathbf{x}\in\mathbb{R}^{N}$ is considered to be an arbitrary temporal
sequence of the distribution ensembles $\left\\{\mathbf{x}(k)\right\\}$ whith
$k\in\left\\{1,\cdots,N_{\mathbf{x}}\right\\}$.

###  II-A Modified Arcsine Law For Non-Stationary Input Signals

We consider a non-zero time-varying Gaussian threshold $\bm{\tau}$ that is
independent of the input signal with the distribution
$\bm{\tau}\sim\mathcal{N}\left(\mathbf{d}=\mathbf{1}d,\bm{\Sigma}\right)$, and
define a new random process $\mathbf{w}$ such that
$\mathbf{w}=\mathbf{x}-\bm{\tau}$. Clearly, $\mathbf{w}$ is a non-stationary
Gaussian process with
$\mathbf{w}\sim\mathcal{N}\left(-\mathbf{d},\bm{R}_{x}+\bm{\Sigma}=\bm{P}\right)$
where $\bm{P}$ is a non-Toeplitz matrix. The autocorrelation function of the
one-bit quantized output process is formulated in the following.

######  Theorem 1.

Suppose $p_{0i}=\bm{P}(i,i)$, $p_{0j}=\bm{P}(j,j)$ and $p_{ij}=\bm{P}(i,j)$,
where $\bm{P}$ is the covariance matrix of $\mathbf{w}$. Consider the one-bit
quantized random variable $\mathbf{y}=f(\mathbf{w})$ where $f(.)$ is the sign
function. Then, the autocorrelation function of $\mathbf{y}$ takes the form

$\displaystyle
R_{\mathbf{y}}(i,j)=\frac{e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}}{\pi\sqrt{\left(p_{0i}p_{0j}-p_{ij}^{2}\right)}}\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}+\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}\right.$
(4)
$\displaystyle\left.-\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{\beta_{n}}Q\left(\frac{\alpha_{n}}{\sqrt{2\beta_{n}}}\right)e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}d\theta\right\\}-1,$

where $\alpha_{n}$ and $\beta_{n}$ are evaluated as

$\displaystyle\alpha_{n}$
$\displaystyle=\frac{d\left(p_{0i}\sin\theta+p_{0j}\cos\theta-
p_{ij}(\cos\theta+\sin\theta)\right)}{(p_{0i}p_{0j}-p_{ij}^{2})},$ (5)
$\displaystyle\beta_{n}$
$\displaystyle=\frac{p_{0j}\cos^{2}\theta+p_{0i}\sin^{2}\theta-p_{ij}\sin
2\theta}{2(p_{0i}p_{0j}-p_{ij}^{2})}.$

###### Proof.

The autocorrelation value associated with lags $i$ and $j$ is given by

$R_{\mathbf{y}}(i,j)=\kappa\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(w_{i})f(w_{j})e^{\lambda(d)}dw_{i}dw_{j}$
(6)

where $\kappa$ and $\lambda(d)$ are defined as

$\kappa\triangleq\left(2\pi\sqrt{p_{0i}p_{0j}-p_{ij}^{2}}\right)^{-1},$ (7)
$\lambda(d)\triangleq\frac{(w_{i}+d)^{2}p_{0j}+(w_{j}+d)^{2}p_{0i}-2p_{ij}(w_{i}+d)(w_{j}+d)}{-2(p_{0i}p_{0j}-p_{ij}^{2})}.$
(8)

The autocorrelation function in (6) can be rewritten as

$\displaystyle R_{\mathbf{y}}(i,j)=$
$\displaystyle\kappa\left(\int_{0}^{\infty}\int_{0}^{\infty}e^{\lambda(d)}\,dw_{i}\,dw_{j}\right.$
(9)
$\displaystyle+\int_{-\infty}^{0}\int_{-\infty}^{0}e^{\lambda(d)}\,dw_{i}\,dw_{j}$
$\displaystyle-\int_{0}^{\infty}\int_{-\infty}^{0}e^{\lambda(d)}\,dw_{i}\,dw_{j}$
$\displaystyle\left.-\int_{-\infty}^{0}\int_{0}^{\infty}e^{\lambda(d)}\,dw_{i}\,dw_{j}\right).$

We can simplify (9) using the relation
$\kappa\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{\lambda(d)}\,dw_{i}\,dw_{j}=1$.
In fact, one can verify that

$R_{\mathbf{y}}(i,j)=2\kappa\int_{0}^{\infty}\int_{0}^{\infty}\left(e^{\lambda(d)}+e^{\lambda(-d)}\right)\,dw_{i}\,dw_{j}-1.$
(10)

By employing polar coordinates $w_{i}=\rho\cos\theta$, $w_{j}=\rho\sin\theta$,
we can recast the integral in (10) as

$\displaystyle
R_{\mathbf{y}}(i,j)=\chi\int_{0}^{\frac{\pi}{2}}\int_{0}^{\infty}e^{-\beta_{n}\rho^{2}}\left(e^{-\alpha_{n}\rho}+e^{\alpha_{n}\rho}\right)\rho\,d\rho\,d\theta-1,$
(11)

where $\alpha_{n}$ and $\beta_{n}$ are readily defined in (5), and

$\chi\triangleq 2\kappa
e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}.$ (12)

Integrating (11) with respect to $\rho$ leads to

$\displaystyle
R_{\mathbf{y}}(i,j)=\frac{e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}}{\pi\sqrt{\left(p_{0i}p_{0j}-p_{ij}^{2}\right)}}\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}+\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}\right.$
(13)
$\displaystyle\left.-\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{\beta_{n}}Q\left(\frac{\alpha_{n}}{\sqrt{2\beta_{n}}}\right)e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}d\theta\right\\}-1.$

∎

It remains to evaluate the integral in (4) in terms of $\\{p_{0i}\\}$,
$\\{p_{0j}\\}$ and $\\{p_{ij}\\}$, which have to be estimated—a task that is
central to our efforts in the rest of this paper. Finding $\\{p_{0i}\\}$,
$\\{p_{0j}\\}$ and $\\{p_{ij}\\}$ results in time-varying input variance and
autocorrelation recovery, which can be achieved by considering the relation:

$\bm{R}_{\mathbf{x}}(i,j)=\bm{P}(i,j)-\bm{\Sigma}(i,j).$ (14)

For $i=j$, the input variance is hence given by
$\bm{R}_{\mathbf{x}}(i,i)=r_{0i}=p_{0i}-\bm{\Sigma}(i,i)$ and
$\bm{R}_{\mathbf{x}}(j,j)=r_{0j}=p_{0j}-\bm{\Sigma}(j,j)$, while for $i\neq
j$, we have the input autocorrelation as
$\bm{R}_{\mathbf{x}}(i,j)=r_{ij}=p_{ij}-\bm{\Sigma}(i,j)$.

##  III Time-Varying Variance Recovery

To recover the time-varying variances $\\{r_{0i}\\}$, the following lemma
would be useful.

######  Lemma 1.

The first moment (mean) of the one-bit sampled data, depends on the threshold
distribution and the power of sampled data via the relation

$\mathbb{E}\left\\{y_{i}\right\\}=2Q\left(\frac{d}{\sqrt{p_{0i}}}\right)-1,\quad\forall
i\in\\{1,\cdots,N\\}.$ (15)

###### Proof.

We have

$\mathbb{E}\left\\{y_{i}\right\\}=\int_{-\infty}^{+\infty}f(w_{i})p(w_{i})\,dw_{i}$
(16)

for $i\in\\{1,\cdots,N\\}$, with $p(w_{i})=\left(\sqrt{2\pi
p_{0i}}\right)^{-1}e^{\frac{-\left(w_{i}+d\right)^{2}}{2p_{0i}}}$. We can
further simplify (16) as

$\displaystyle\mathbb{E}\left\\{y_{i}\right\\}$
$\displaystyle=-\int_{-\infty}^{0}p(w_{i})\,dw_{i}+\int_{0}^{\infty}p(w_{i})\,dw_{i}$
(17) $\displaystyle=2\int_{0}^{+\infty}p(w_{i})\,dw_{i}-1$
$\displaystyle=2Q\left(\frac{d}{\sqrt{p_{0i}}}\right)-1$

which completes the proof. ∎

In light of the above, a relation between the input variance and the mean of
one-bit sampled data is established, which provides an additional avenue to
estimate the variances $\left\\{p_{0i}\right\\}$. More precisely, according to
Lemma 1 and (14), the input time-varying variances $\left\\{r_{0i}\right\\}$
are given by

$\displaystyle r_{0i}^{\star}$
$\displaystyle=\left(\frac{d}{Q^{-1}\left(\frac{\mu_{i}+1}{2}\right)}\right)^{2}-\sigma^{2}_{\bm{\tau}},\quad
i\in\left\\{1,\cdots,N\right\\},$ (18)

where $\\{r_{0i}^{\star}\\}$ denote the optimal values of $\\{r_{0i}\\}$,
$\sigma^{2}_{\bm{\tau}}$ is the threshold variance, and $\\{\mu_{i}\\}$ denote
the entries of $\bm{\mu}$ which may be estimated via the sample mean
$\frac{1}{N_{\mathbf{x}}}\sum^{N_{\mathbf{x}}}_{k=1}\mathbf{y}(k)$ [34].

(a)

(b)

Figure 1: Average MSE for time-varying input variance recovery with (a) $i=2$
($r_{02}$), and (b) $j=8$ ($r_{08}$), for different one-bit sample sizes.

###  III-A Numerical Results

We will examine the effectiveness of (18) in estimating the time-varying input
variances. In all experiments, the input signals were generated as zero-mean
Wiener process sequences with time-varying variance ranging from $0.2$ to
$0.8$. Also, the number of states is $100$ (i.e.,
$\mathbf{x}\in\mathbb{R}^{100}$). Accordingly, we made use of the time-varying
thresholds with $d=0.5$ and diagonal $\bm{\Sigma}$ whose diagonal entries are
set to $0.2$. In the non-stationary input signal case, the true input variance
is not a constant number. Therefore, we define the experimental mean square
error (MSE) of the estimate $\hat{r}_{0i}$ of a variance $r_{0i}$ as

$\displaystyle\mathrm{MSE}=\frac{1}{E}\sum^{E}_{e=1}|r^{e}_{0i}-{\hat{r}^{e}_{0i}}|^{2},$
(19)

where $\left\\{r^{e}_{0i},\hat{r}^{e}_{0i}\right\\}$ are the time-varying
variances and their estimates in the $e$-th experiment. Also, the number of
experiments is assumed to be $E=15$. As can be seen in Fig. 1, we can
accurately estimate the time-varying variance elements of an input signal
based on (18) for $i=2$ and $j=8$ ($r_{02}\neq r_{08}$). The results are
obtained for the number of ensembles
$N_{\mathbf{x}}\in\left\\{1000,3000,6000,10000\right\\}$, with fixed $d$ and
$\bm{\Sigma}$ for each experiment. It is observed that the accuracy of time-
varying variance recovery will significantly improve as the number of one-bit
samples grows large.

To further investigate the effectiveness of our proposed approach, we generate
a non-stationary Gaussian process according to
$\mathbf{x}\sim\mathcal{N}\left(0,\sigma^{2}_{t}\bm{I}\right)$, where
$\\{\sigma^{2}_{t}\\}$ are generated based on the generalized autoregressive
conditional heteroskedasticity (GARCH) model with order one, i.e.
$\operatorname{GARCH}\left(1,1\right)$, which may be written as [33],

$\sigma^{2}_{t}=\zeta_{0}+\zeta_{1}\sigma^{2}_{t-1}+\zeta_{2}\epsilon^{2}_{t-1},\quad
x_{t}=\epsilon_{t}|\psi_{t-1},$ (20)

where $\\{x_{t}\\}$ are elements of $\mathbf{x}$, $\epsilon_{t}|\psi_{t-1}$
denotes the conditional random variable $\epsilon_{t}$ given its previous
ensembles set $\psi_{t-1}$, and
$\left\\{\zeta_{0},\zeta_{1},\zeta_{2}\right\\}$ are our GARCH model
parameters. In Fig. 2, we present an example of time-varying variance sequence
recovery. The true input signal time-varying variance $\sigma^{2}_{t}$ and the
estimated values by our approach are presented when $t$ is a temporal sequence
of length $20$.

So far, we have obtained the time-varying variance elements of the input
covariance matrix. To recover the input autocorrelation values
($\\{r_{ij}\\},i\neq j$), the integral in (4) should be evaluated which
appears to be difficult to find in closed-form. Therefore, in the following,
we deploy various approximations to facilitate its evaluation, which enables
to the recovery of all elements of the covariance matrix
$\bm{R}_{\mathbf{x}}$.

Figure 2: Recovery of the input signal time-varying variance $\sigma^{2}_{t}$
generated by a $\operatorname{GARCH}\left(1,1\right)$ model based on (18),
when $t$ is a temporal sequence of length $20$, with the true values plotted
along the estimates.

##  IV Analytic Approach for Covariance Recovery

The first part of the integration in (4) can be analytically evaluated as

$\displaystyle\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}\,d\theta$
$\displaystyle=\sqrt{p_{0i}p_{0j}-p_{ij}^{2}}\left(\pi+2\sin^{-1}\left[\frac{p_{ij}}{\sqrt{p_{0i}p_{0j}}}\right]\right).$
(21)

If the threshold is considered as a zero-mean Gaussian process
($\bm{\tau}\sim\mathcal{N}\left(0,\bm{\Sigma}\right)$), one can resort to the
well-known _arcsine law_ relation for non-stationary signals, i.e.,
$R_{\mathbf{y}}(i,j)=\frac{2}{\pi}\sin^{-1}\left(\frac{p_{ij}}{\sqrt{p_{0i}p_{0j}}}\right)$.
However, in the general case, we evaluate all parts of the integration in (4).
Computing the integral in (4) with the integrands
$D_{1}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)=\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{\beta_{n}}Q\left(\frac{\alpha_{n}}{\sqrt{2\beta_{n}}}\right)e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}$
and
$D_{2}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)=\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}$
with respect to $\theta$ appears to be a difficult task. Thus, in the
following section, the Padé approximation (PA) [36, 37, 38] is utilized to
approximate said integrands $D_{1}$ and $D_{2}$. This facilitates the recovery
of $\\{p_{ij}\\}$ in Section IV-B.

###  IV-A Padé Approximation

As in Part i@ of this work [2], we use PA to approximate $D_{1}$ and $D_{2}$.
Note that the integration in (4) occurs in the interval
$\theta\in\left[0,\frac{\pi}{2}\right]$. To have a better fitness, we again
use the idea of piece-wise PA with three distinct intervals
$\left[0,\frac{\pi}{8}\right]$, $\left[\frac{\pi}{8},\frac{3\pi}{8}\right]$,
and $\left[\frac{3\pi}{8},\frac{\pi}{2}\right]$ with the expansion points
$\theta\in\left\\{0,\frac{\pi}{4},\frac{\pi}{2}\right\\}$. Consequently, the
function $D_{2}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)$ is approximated as

$\displaystyle\theta$
$\displaystyle\in\left[0,\frac{\pi}{8}\right]\cup\left[\frac{\pi}{8},\frac{3\pi}{8}\right]\cup\left[\frac{3\pi}{8},\frac{\pi}{2}\right]:$
(22) $\displaystyle D_{2}$
$\displaystyle=\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}\approx\frac{e+s\theta}{k+g\theta+h\theta^{2}},$

A similar approximation can be proposed for
$D_{1}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)$. It is straightforward to
verify that the two functions
$D_{1}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)$ and
$D_{2}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)$ are analytic at the
expansion points. Also, the $Q$-function in (4) is approximated by the
$\bar{Q}$-function as [39]:

$\bar{Q}\left(x\right)=\frac{1}{12}e^{\frac{-x^{2}}{2}}+\frac{1}{4}e^{\frac{-2x^{2}}{3}},\quad
x>0.$ (23)

Substituting $D_{2}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)$ with its
approximation and evaluating the integration in the associated parts of (4)
results in:

$\displaystyle\int_{0}^{\frac{\pi}{8}}$
$\displaystyle\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}\,d\theta\approx\frac{s}{2h}\ln{\left(\frac{\left|k+\frac{\pi
g}{8}+\frac{\pi^{2}h}{64}\right|}{\left|k\right|}\right)}+$ (24)
$\displaystyle\frac{2eh-sg}{h\sqrt{4hk-g^{2}}}\tan^{-1}\left(\frac{\pi
h\sqrt{4hk-g^{2}}}{16hk+\pi gh}\right),$
$\displaystyle\int_{\frac{\pi}{8}}^{\frac{3\pi}{8}}$
$\displaystyle\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}\,d\theta\approx\frac{s}{2h}\ln\left(\frac{\left|64k+9\pi^{2}h+24\pi
hg\right|}{\left|64k+\pi^{2}h+8\pi hg\right|}\right)+$ (25)
$\displaystyle\frac{2eh-sg}{h\sqrt{4kh-g^{2}}}\tan^{-1}\left(\frac{8\pi
h\sqrt{4hk-g^{2}}}{64kh+3\pi^{2}h^{2}+16\pi hg}\right),$
$\displaystyle\int_{\frac{3\pi}{8}}^{\frac{\pi}{2}}$
$\displaystyle\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}\,d\theta\approx\frac{s}{2h}\ln\left(\frac{\left|k+\frac{\pi
g}{2}+\frac{\pi^{2}h}{4}\right|}{\left|k+\frac{3\pi
g}{8}+\frac{9\pi^{2}h}{64}\right|}\right)+$ (26) $\displaystyle\frac{2eh-
sg}{h\sqrt{4kh-g^{2}}}\tan^{-1}\left(\frac{\pi
h\sqrt{4hk-g^{2}}}{16kh+3\pi^{2}h^{2}+7\pi hg}\right).$

Similar approximations can be obtained for terms associated with the function
$D_{1}\left(\theta;p_{0i},p_{0j},p_{ij},d\right)$.

###  IV-B Recovery Criterion

Based on our discussions in Section III, $p_{0i}^{\star}$ and $p_{0j}^{\star}$
may be immediately obtained by (18). Then, $\\{p_{ij}\\}$ are estimated by
formulating a minimization problem. For this purpose, one may consider the
following criterion:

$\displaystyle\bar{G}(p_{0i},p_{0j},p_{ij})\triangleq\log\left(\left|R_{\mathbf{y}}(i,j)-\chi\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}\right.\right.\right.$
(27)
$\displaystyle\left.\left.\left.+\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}-\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{\beta_{n}}Q\left(\frac{\alpha_{n}}{\sqrt{2\beta_{n}}}\right)e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}d\theta\right\\}+1\right|^{2}\right),$

where the autocorrelation of output signal ($R_{\mathbf{y}}$) can be estimated
via the sample covariance matrix (SCM) [35],

$\displaystyle\bm{R}_{\mathbf{y}}\approx\frac{1}{N_{\mathbf{x}}}\sum_{k=1}^{N_{\mathbf{x}}}\mathbf{y}(k)\mathbf{y}(k)^{\mathrm{H}},$
(28)

with $\\{\mathbf{y}(k)\\}$ being the observed sign vectors, and $\chi$ being
the same as defined in (12). Note that by now we have derived an approximated
version of (4) using PA. Let $H_{n}(p_{0i},p_{0j},p_{ij})$ denote this
approximation. Therefore, we can alternatively consider the criterion:

$\displaystyle G(p_{ij})$
$\displaystyle\triangleq\log\left(\left|R_{\mathbf{y}}(i,j)-H_{n}(p_{0i}^{\star},p_{0j}^{\star},p_{ij})\right|^{2}\right).$
(29)

A numerical investigation of (29) reveals that it is multi-modal, i.e. with
multiple local minima—see Fig. 3 for an example of the optimization landscape
of $G(p_{ij})$. Taking the feasible region of $\left\\{p_{ij}\right\\}$ into
account, we can formulate the recovery problem:

$\displaystyle\mathcal{P}_{i,j}$ $\displaystyle:$ $\displaystyle\min_{p_{ij}}$
$\displaystyle G(p_{ij}),$ s.t. $\displaystyle-p_{m}\leq p_{ij}\leq p_{m},$
(30)

where $p_{m}=\operatorname{min}\\{[p_{0i}^{\star},p_{0j}^{\star}]\\}$. The
non-convex problem in (30) may be solved via the gradient descent numerical
optimization approach by employing multiple random initial points. Once
$p_{ij}$ is estimated, we can estimate the autocorrelation values of
$\mathbf{x}$ via (14). The acquired optimum recovery results will be presented
in the following.

Figure 3: Example plot of the estimation criterion $G(p_{ij})$ with respect to
$p_{ij}$ showing its multi-modality, i.e. having multiple local optima.

###  IV-C Numerical Results

We will examine the effectiveness of the PA method by comparing its recovery
results with the true input signal autocorrelation values in the non-
stationary case. In all experiments, the input signals were generated as zero-
mean Wiener process sequences with time-varying variance ranging from $0.2$ to
$0.8$. The number of states is set to $100$ (i.e.,
$\mathbf{x}\in\mathbb{R}^{100}$). Accordingly, we make use of the time-varying
thresholds with $d=0.5$ and diagonal $\bm{\Sigma}$ whose diagonal entries are
set to $0.2$.

Figure 4: Recovery of the input signal autocorrelation $r_{ij}$ using the PA
approach for $i=2$ with $j$ being a temporal sequence of length $13$, with the
true values plotted along the estimates.

We first present an example of autocorrelation sequence recovery. The true
input signal autocorrelation and the estimated autocorrelation values by our
approach are shown in Fig. 4, where $i=2$ and $j$ is a temporal sequence of
length $13$. Fig. 4 appears to confirm the possibility of recovering the
autocorrelation values from one-bit sampled data with time-varying thresholds
in the non-stationary case.

###  IV-D Fitness Analysis of the Proposed Approximations

In this section, we further examine the capability of the PA approach to
approximate the integrands in (4). In fact, it appears a precise approximation
of the integrands is connected to the one-bit comparison thresholds that are
used for sampling. To see how, note that the exponential term in (4), i.e.
$e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}$, should remain bounded to guarantee a
well-behaved Padé approximation. This is due to the fact that, in (4),
$\beta_{n}$ is non-zero and $e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}$ is the
only term that can grow very fast and create a round-off numerical error.
Consider a bounding of this term in the form

$e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}<\gamma_{1}.$ (31)

Since $\ln(\cdot)$ is a strictly increasing function, Eq. (31) can be
alternatively expressed as:

$\frac{\alpha_{n}^{2}}{4\beta_{n}}<\ln(\gamma_{1}).$ (32)

Note that $\alpha_{n}$ in (5) is _directly scaled by the sampling threshold
mean $d$_. By setting $\alpha_{n}=d\alpha_{e}$, (32) can be written as

$d^{2}\left(\frac{\alpha_{e}^{2}}{4\beta_{n}}\right)<\ln(\gamma_{1}),$ (33)

where $\alpha_{n}$ and $\beta_{n}$ are defined in (5). To guarantee the bound
in (33), we should have:

$d^{2}\max_{\theta}\left\\{\frac{\alpha_{e}^{2}}{4\beta_{n}}\mathbb{I}_{\left(0\leq\theta\leq\frac{\pi}{2}\right)}\right\\}<\ln(\gamma_{1}).$
(34)

The inner optimization problem
$\max_{\theta}\frac{\alpha_{e}^{2}}{4\beta_{n}}\mathbb{I}_{\left(0\leq\theta\leq\frac{\pi}{2}\right)}$
can be solved via golden section search and parabolic interpolation method
[40]. Note that $\gamma_{1}$ is greater than one. To prove this claim, it is
sufficient to show

$\exists\leavevmode\nobreak\
\theta\in\left[0,\frac{\pi}{2}\right]:\quad\beta_{n}(\theta)>0,$ (35)

since $\alpha_{e}^{2}$ is always non-negative. By plugging in $\theta=0$ in
$\beta_{n}$, we obtain
$\beta_{n}\bigg{|}_{\theta=0}=0.5p_{0j}\left(p_{0i}p_{0j}-p_{ij}^{2}\right)^{-1}$
which is always positive.

To assess the goodness of the considered approximations, consider the
integrand term:

$\displaystyle\Delta(\theta;p_{0i},p_{0j},p_{ij})$
$\displaystyle=\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{2\beta_{n}}e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}$
(36)
$\displaystyle-\sqrt{\frac{\pi}{\beta_{n}}}\frac{\alpha_{n}}{\beta_{n}}\bar{Q}\left(\frac{\alpha_{n}}{\sqrt{2\beta_{n}}}\right)e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}.$

We compare $\Delta(\cdot)$ with its approximated forms for
$\theta\in\left[0,\frac{\pi}{8}\right]\cup\left[\frac{\pi}{8},\frac{3\pi}{8}\right]\cup\left[\frac{3\pi}{8},\frac{\pi}{2}\right]$,
whose results are plotted in Fig. 5 with parameters $p_{0i}=0.8$,
$p_{0j}=0.7$, $p_{ij}=0.05$, and $d=0.7$. Note that, in this case, (34) is
satisfied by considering $\gamma_{1}=2$. As can be seen, PA appears to promise
good fitness, with a small mean square error (MSE) of $\sim 1.1699e-04$.

Figure 5: Comparing $\Delta$ with its approximated piece-wise function using
PA.

##  V Gaussian Quadrature Technique for Covariance Recovery

In this section, the Gauss-Legendre quadrature approach is adopted to evaluate
the integration in (4). This lays the ground for the recovery of
$\left\\{p_{ij}\right\\}$ since $p_{0i}$ and $p_{0j}$ are obtained by (18). At
first, an approximated version of (4) is obtained based on the Gauss-Legendre
quadrature technique in Section V-A. Then, a criterion will be presented to
recover $\left\\{p_{ij}\right\\}$, and subsequently, the input autocorrelation
values in Section V-B. Finally, the efficacy of this approach in estimating
the input autocorrelation values is numerically evaluated.

###  V-A Gauss-Legendre Quadrature Method for Integral Approximation

As discussed in Part i@ of this work, the central assumption to the use of the
Gauss-Legendre quadrature technique is that the integrand $f(x)$ should be
finite within the domain of integration. The integrands in (4) meet this
assumption; it is easy to verify that $\textbf{num}(\beta_{n})\neq 0$, where
$\textbf{num}(\cdot)$ denotes the numerator of the fractional argument.
Therefore, by employing the Gauss-Legendre quadrature technique, the relation
in (4) can be approximated as

$\displaystyle R_{\mathbf{y}}(i,j)$
$\displaystyle\approx\frac{e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}}{\pi\sqrt{\left(p_{0i}p_{0j}-p_{ij}^{2}\right)}}\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}d\theta\right.$
(37)
$\displaystyle\left.-\frac{\pi}{4}\sum_{\varepsilon=1}^{N_{q}}\omega_{\varepsilon}D_{1}\left(\frac{\pi}{4}(\theta_{\varepsilon}+1);p_{0i},p_{0j},p_{ij},d\right)\right.$
$\displaystyle\left.+\frac{\pi}{4}\sum_{\varepsilon=1}^{N_{q}}\omega_{\varepsilon}D_{2}\left(\frac{\pi}{4}(\theta_{\varepsilon}+1);p_{0i},p_{0j},p_{ij},d\right)\right\\}-1,$

where $\theta_{\varepsilon}$ denotes the $\varepsilon$-th Gauss node. Note
that the first part of the above integration was readily given in closed-form
in (21).

Figure 6: Example plot of the Gauss-Legendre quadrature approach-based
estimation criterion $\Psi(p_{ij})$ with respect to $p_{ij}$ showing its
convexity.

###  V-B Covariance Recovery via Convex Optimization

Based on our discussion in Section III, the values for $p_{0i}^{\star}$ and
$p_{0j}^{\star}$ are simply given by (18). The parameters of interest
$\left\\{p_{ij}\right\\}$ are then estimated by formulating a minimization
problem; namely, we consider the following criterion:

$\displaystyle\bar{\Psi}(p_{ij})\triangleq\log\left(\left|R_{\mathbf{y}}(i,j)-\chi\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}d\theta\right.\right.\right.$
(38)
$\displaystyle\left.\left.\left.-\frac{\pi}{4}\sum_{\varepsilon=1}^{N_{q}}\omega_{\varepsilon}D_{1}\left(\frac{\pi}{4}(\theta_{\varepsilon}+1);p_{0i}^{\star},p_{0j}^{\star},p_{ij},d\right)\right.\right.\right.$
$\displaystyle\left.\left.\left.+\frac{\pi}{4}\sum_{\varepsilon=1}^{N_{q}}\omega_{\varepsilon}D_{2}\left(\frac{\pi}{4}(\theta_{\varepsilon}+1);p_{0i}^{\star},p_{0j}^{\star},p_{ij},d\right)\right\\}+1\right|^{2}\right),$

for which the autocorrelation of output signal $R_{\mathbf{y}}$ can be
estimated using the SCM in (28), and $\chi$ is the same as that in (12).
Recall that we have obtained an approximated version of (4) using the Gauss-
Legendre quadrature in (37). Let $J_{n}(p_{ij})$ denote this approximation. As
a result, we can alternatively use the criterion:

$\Psi(p_{ij})\triangleq\log\left(\left|R_{\mathbf{y}}(i,j)-J_{n}(p_{ij})\right|^{2}\right).$
(39)

It is interesting to note that the criterion in (39) is a convex function with
respect to $p_{ij}$ (a proof is provided in Appendix A)—see Fig. 6 for an
example of the optimization landscape of $\Psi(p_{ij})$. By considering the
feasible region of $p_{ij}$, the following recovery problem is obtained:

$\displaystyle\mathcal{P}_{i,j}$ $\displaystyle:$ $\displaystyle\min_{p_{ij}}$
$\displaystyle\Psi(p_{ij}),$ s.t. $\displaystyle-p_{m}\leq p_{ij}\leq p_{m},$
(40)

where $p_{m}$ is defined in Section IV-B. The convex problem in (40) may be
solved efficiently using the golden section search and parabolic interpolation
approach. Once $\left\\{p_{ij}\right\\}$ is obtained, one can estimate the
autocorrelation values of $\mathbf{x}$ via (14). The recovery results will be
presented in the following.

###  V-C Numerical Results

We examine the usefulness of the Gauss-Legendre quadrature technique by
comparing its recovery results with the true input signal autocorrelation
values in the non-stationary case. In all experiments, the input signals were
generated as zero-mean Gaussian sequences with time-varying variance ranging
from $0.2$ to $0.8$. Accordingly, we made use of the time-varying thresholds
with $d=0.3$ and diagonal $\bm{\Sigma}$ whose diagonal entries are equal to
$0.1$.

We present an example of autocorrelation sequence recovery. The true input
signal autocorrelation and our estimated autocorrelation values are shown in
Fig. 7 with $i=2$ and $j$ being a temporal sequence of length $13$. Fig. 7
appears not only to confirm the possibility of recovering the autocorrelation
values from one-bit sampled data with time-varying thresholds in the non-
stationary case but also the effectiveness of the Gauss-Legendre technique.

Figure 7: Recovery of the input autocorrelation $r_{ij}$ using the Gauss-
Legendre technique, with $i=2$ and $j$ being a temporal sequence of length
$13$. The true values are plotted alongside the estimates.

##  VI Monte-Carlo Integration for Covariance Recovery

In this section, the Monte-Carlo integration approach is utilized to evaluate
the integral in (4). At first, we formulate an approximated version of (4)
based on the Monte-Carlo integration approach in Section VI-A. We then present
a new criterion to recover $\\{p_{ij}\\}$ based on this approximation. The
efficacy of this approach in estimating the input autocorrelation values is
numerically evaluated.

###  VI-A Monte-Carlo Method for Integral Evaluation

As mentioned in Part i@ of this work, the Monte-Carlo integration technique
can be utilized to recover the input signal covariance matrix. Here, the same
idea is adopted to evaluate the integration in (4) for the non-stationary
input signal scenario. More concretely, by employing the Monte-Carlo
integration technique, one can approximate (4) as follows:

$\displaystyle R_{\mathbf{y}}(i,j)$
$\displaystyle\approx\frac{e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}}{\pi\sqrt{\left(p_{0i}p_{0j}-p_{ij}^{2}\right)}}\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}d\theta\right.$
(41)
$\displaystyle\left.-\frac{\pi}{2N_{m}}\sum_{\varepsilon=1}^{N_{m}}D_{1}\left(\theta_{\varepsilon};p_{0i},p_{0j},p_{ij},d\right)\right.$
$\displaystyle\left.+\frac{\pi}{2N_{m}}\sum_{\varepsilon=1}^{N_{m}}D_{2}\left(\theta_{\varepsilon};p_{0i},p_{0j},p_{ij},d\right)\right\\}-1,$

where $\theta_{\varepsilon}$ denotes the $\varepsilon$-th random number
generated from the uniform distribution in the interval
$\left[0,\frac{\pi}{2}\right]$. Note that the first part of the above integral
was readily evaluated in closed-form in (21).

###  VI-B Convex Covariance Recovery

Similar the previous proposed approaches, we begin by estimating $p_{0i}$ and
$p_{0j}$ through (18). We then aim at estimating the unknown parameters
$\\{p_{ij}\\}$ by formulating a minimization problem. Namely, we consider the
following criterion:

$\displaystyle\bar{\Gamma}(p_{ij})$
$\displaystyle\triangleq\log\left(\left|R_{\mathbf{y}}(i,j)-\chi\left\\{\int_{0}^{\frac{\pi}{2}}\frac{1}{\beta_{n}}d\theta\right.\right.\right.$
(42)
$\displaystyle\left.\left.\left.-\frac{\pi}{2N_{m}}\sum_{\varepsilon=1}^{N_{m}}D_{1}\left(\theta_{\varepsilon};p_{0i}^{\star},p_{0j}^{\star},p_{ij},d\right)\right.\right.\right.$
$\displaystyle\left.\left.\left.+\frac{\pi}{2N_{m}}\sum_{\varepsilon=1}^{N_{m}}D_{2}\left(\theta_{\varepsilon};p_{0i}^{\star},p_{0j}^{\star},p_{ij},d\right)\right\\}+1\right|^{2}\right),$

where the autocorrelation of output signal $R_{\mathbf{y}}$ can be estimated
via (28). Let $F_{n}(p_{ij})$ denote the approximation of (4) using the Monte-
Carlo integration. Therefore, we can consider the following alternative
criterion:

$\Gamma(p_{ij})\triangleq\log\left(\left|R_{\mathbf{y}}(i,j)-F_{n}(p_{ij})\right|^{2}\right).$
(43)

Similar to the previous criterion in (39), $\Gamma(p_{ij})$ appears to be a
convex function with respect to $p_{ij}$, whose proof of convexity is similar
to that for $\Psi_{m}(.)$ in Appendix A—see Fig. 8 for an example of the
optimization landscape associated with $\Gamma(p_{ij})$. By considering the
feasible region of the parameter of interest $\\{p_{ij}\\}$, one can formulate
the following recovery problem:

$\displaystyle\mathcal{P}_{i,j}$ $\displaystyle:$ $\displaystyle\min_{p_{ij}}$
$\displaystyle\Gamma(p_{ij}),$ s.t. $\displaystyle-p_{m}\leq p_{ij}\leq
p_{m},$ (44)

where $p_{m}$ is defined in Section IV-B. The convex problem in (44) may be
tackled by the same tools as proposed in Section V-B. Recovery of
$\\{p_{ij}\\}$ leads to estimating the autocorrelation values of $\mathbf{x}$
via (14). The optimum recovery results will be presented in the following.

Figure 8: Example plot of the estimation criterion $\Gamma(p_{ij})$ with
respect to $p_{ij}$ showing its convexity.

###  VI-C Numerical Results

Herein, we examine the Monte-Carlo integration technique by comparing its
recovery results with the true input signal autocorrelation values in the non-
stationary case. In all experiments, the input signals were generated as zero-
mean Gaussian sequences with time-varying variances ranging from $0.2$ to
$0.8$. Accordingly, we made use of the time-varying thresholds with $d=0.3$
and diagonal $\bm{\Sigma}$ whose diagonal entries are set to $0.1$.

We present an example of autocorrelation sequence recovery. The true input
signal autocorrelation and the estimated autocorrelation values are shown in
Fig. 9, with $i=2$ and $j$ being a temporal sequence of length $13$. It can be
observed from Fig. 9 that the Monto-Carlo based approach presents satisfactory
recovery results in the non-stationary case as well.

Figure 9: Recovery of the input autocorrelation $r_{ij}$ using the Monte-Carlo
integration approach, with $i=2$ and $j$ being a temporal sequence of length
$13$. The true values are plotted alongside the estimates.

##  VII Comparing The Proposed Recovery Methods

It would be of interest to compare the discussed covariance recovery
approaches in the non-stationary setting: (i) employing the Padé approximation
of the integrands in (4), (ii) applying the Gauss-Legendre quadrature
technique, and (iii) applying the Monte-Carlo integration to evaluate the
integral in (4). To this end, we generate a non-stationary input signal
$\mathbf{x}\in\mathbb{R}^{5}$ with the ensemble length $N_{\mathbf{x}}=10000$
and the following non-Toeplitz covariance matrix:

$\displaystyle\bm{R}^{\star}_{\mathbf{x}}=\left(\begin{array}[]{ccccc}+0.5040&-0.0065&+0.0015&-0.0036&+0.0044\\\
-0.0065&+0.2565&-0.0034&+0.0086&+0.0031\\\
+0.0015&-0.0034&+0.3298&+0.0063&+0.0031\\\
-0.0036&+0.0086&+0.0063&+0.6376&-0.0062\\\
+0.0044&+0.0031&+0.0031&-0.0062&+0.4552\end{array}\right).$ (45)

The time-varying threshold is generated using the same settings as in Section
IV-C. Table I illustrates the squared Frobenius norm of the error, normalized
by the squared Frobenius norm of the desired covariance matrix
$\bm{R}^{\star}_{\mathbf{x}}$:

$\mathrm{NMSE}\triangleq\frac{\left\|\bm{R}^{\star}_{\mathbf{x}}-\hat{\bm{R}}_{\mathbf{x}}\right\|^{2}_{\mathrm{F}}}{\left\|\bm{R}^{\star}_{\mathbf{x}}\right\|^{2}_{\mathrm{F}}},$
(46)

where $\hat{\bm{R}}_{\mathbf{x}}$ is the recovered covariance matrix. The
presented results are averaged over $5$ experiments.

In the non-stationary input signal scenario, similar to the stationary case,
all three approaches show promising recovery results—see Table I. The Gauss-
Legendre method has a better performance in recovering the input signal
autocorrelation values in comparison with the PA technique and the Monte-Carlo
integration. It is also worth noting that the two proposed numerical
approaches other than the PA technique boil down to simplified convex
programs, hence ensuring convergence to the global optimum. However, a proper
selection of the number of nodes and quadrature points in the Gauss-Legendre
quadrature and the Monte-Carlo integration techniques is crucial and may
present itself as a bottleneck in an effective recovery. This is not an
obstacle in applying the PA technique. As a result, one may wish to run the
PA-based recovery to help with the proper deployment of the other two
techniques.

Table I: Average NMSE value for the Covariance Recovery. Covariance recovery approach | NMSE
---|---
PA technique | $7.813e-05$
Gauss-Legendre quadrature | $2.093e-05$
Monte-Carlo integration | $2.488e-05$

##  VIII Judicious Selection of Sampling Thresholds

While the use of time-varying thresholds for one-bit sampling has shown
promise in various signal recovery problems, tuning the applied thresholds
provides both an opportunity and a challenge. In this section, we will discuss
an approach to effectively set the sampling threshold mean $d$—whose
significance was already shown in our analysis in Section IV-D. The value of
the threshold mean $d$ is one of the parameters in our recovery cost
functions, and consequently, impacts the effectiveness of the autocorrelation
sequence recovery by various proposed approaches.

###  VIII-A Problem Formulation for Threshold Mean Optimization

Consider a set of thresholds distributed as
$\bm{\tau}\sim\mathcal{N}\left(\mathbf{d}=\mathbf{1}d,\bm{\Sigma}=\sigma^{2}_{\bm{\tau}}\bm{I}\right)$.
To design _optimal_ thresholds that are independent from the unknown zero-mean
Gaussian signal ($\mathbf{x}$), we use the CDF of the observed sign data
$\mathbf{y}$ and formulate a maximum likelihood estimation (MLE) problem. The
goal is to determine the threshold mean $d$ solely from the sign data
$\mathbf{y}$. The one-bit samples are generated as

$\displaystyle i\in\left\\{1,\cdots,N\right\\},\quad y_{i}$
$\displaystyle=\begin{cases}+1&x_{i}>\tau_{i},\\\
-1&x_{i}<\tau_{i}.\end{cases}$ (47)

The probability vector $\bm{p}$ for the one-bit measurement vector
$\mathbf{y}$ may be written as

$\displaystyle p_{\mathbf{y}}\left(y_{i}|\tau_{i}\right)=p_{i}$
$\displaystyle=\begin{cases}1-\Psi(\tau_{i})&\text{for}\quad\\{y_{i}=+1\\},\\\
\Psi(\tau_{i})&\text{for}\quad\\{y_{i}=-1\\},\end{cases}$ (48)

where $\Psi(.)$ is the CDF of $\mathbf{x}$, as defined in (3). The associated
log-likelihood function is hence given by

$\displaystyle\mathcal{L}_{\mathbf{y}}(\bm{r}_{0},\bm{\tau})$
$\displaystyle=\sum^{N}_{i=1}\left\\{\mathbb{I}_{(y_{i}=+1)}\log\left(1-\Psi(\tau_{i})\right)\right.$
(49)
$\displaystyle\left.+\mathbb{I}_{(y_{i}=-1)}\log\left(\Psi(\tau_{i})\right)\right\\},$

where $\bm{r}_{0}$ is a vector containing the diagonal entries of the
covariance matrix of the input signal $\mathbf{x}$. As mentioned earlier, in
the non-stationary scenario, the covariance matrix has an arbitrary non-
Toeplitz structure. The entries of $\bm{r}_{0}$, appearing in the CDF, are the
variances for the elements of the input signal. To immediately derive our
desired parameter $d$ from (49), we define the following statistical linear
model for our threshold $\bm{\tau}$:

$\bm{\tau}=\mathbf{1}d+\sigma^{2}_{\bm{\tau}}\leavevmode\nobreak\
\mathbf{z},\quad\mathbf{z}\sim\mathcal{N}\left(0,\bm{I}\right).$ (50)

Therefore, the MLE is formulated as

$\displaystyle\min_{d,\sigma^{2}_{\bm{\tau}}}$
$\displaystyle-\mathcal{L}_{\mathbf{y}}(\bm{r}_{0},\bm{\tau})$ (51) s.t.
$\displaystyle
r_{0i}=\left(\frac{d}{Q^{-1}\left(\frac{\mu_{i}+1}{2}\right)}\right)^{2}-\sigma^{2}_{\bm{\tau}},$
$\displaystyle i\in\left\\{1,\cdots,N\right\\},$

where the equality constraint is obtained from (15).

###  VIII-B Numerical Illustrations for Threshold Mean Design

To numerically scrutinize our approach, the input signal is generated using
the same settings as described in Section IV-C. The desired time-varying
threshold is generated as a Gaussian process with $d^{\star}=0.3$ and
$\Sigma=0.1\bm{I}$. The sign data $\mathbf{y}$ was generated accordingly to be
utilized in order to estimate the desired threshold mean by the MLE problem in
(51). The results are presented in Fig. 10 based on the NMSE between the
desired threshold mean $d^{\star}$ and the recovered mean $\hat{d}$, defined
as:

$\displaystyle\mathrm{NMSE}\triangleq\frac{|d^{\star}-{\hat{d}}|^{2}}{|d^{\star}|^{2}}.$
(52)

Each presented data point is averaged over $5$ experiments. As can be seen in
Fig. 10, the proposed method can accurately estimate the mean of a time-
varying threshold. The results are obtained with sequence lengths
$N_{\mathbf{x}}\in\left\\{1000,3000,6000,10000\right\\}$. Moreover,
$\sigma^{2}_{\bm{\tau}}$ is estimated using (51) with the average NMSE of
$\sim 1e-03$ by considering $N_{\mathbf{x}}=10000$.

Figure 10: Average NMSE for the estimated sampling threshold mean $\hat{d}$,
based on MLE problem presented in (51).

##  IX Modified Bussgang Law For Time-Varying Sampling Thresholds

The modified Bussgang law for stationary input signals was derived in Part i@
of this work. This modified Bussgang law presents a useful relation in
stochastic analysis of stationary input signals when they are sampled with
time-varying thresholds. In this section, the modified Bussgang law is
extended to the case when the non-stationary input signals are considered in
such settings.

(a) Padé approximation technique

(b) Gauss-Legendre quadrature technique

(c) Monte-Carlo integration technique

Figure 11: The recovery of the cross-correlation between the input signal and
the one-bit sampled data by the modified Bussgang law applied in conjunction
with various one-bit autocorrelation recovery approaches for a sequence of
length $13$, with the true values plotted alongside the estimates.

###  IX-A Modified Bussgang Law for Non-Stationary Input Signals

By considering time-varying thresholds, the cross-correlation matrix between
the one-bit sampled data and the non-stationary input signal can be formulated
as follows.

######  Theorem 2.

Suppose
$\bm{\tau}\sim\mathcal{N}\left(\mathbf{d}=\mathbf{1}d,\bm{\Sigma}\right)$ is a
time-varying threshold, and $\mathbf{x}$ is a non-stationary input signal. Let
$\mathbf{y}=g\mathbf{(w)}$ denote the one-bit sampled data, where
$\mathbf{w}=\mathbf{x}-\bm{\tau}$ is distributed as
$\mathbf{w}\sim\mathcal{N}\left(-\mathbf{d},\bm{R}_{\mathbf{x}}+\bm{\Sigma}=\bm{P}\right)$,
with $p_{0j}=\bm{P}(j,j)$ and $p_{ij}=\bm{P}(i,j)$. Then, the cross-
correlation matrix between $\mathbf{y}$ and $\mathbf{x}$ satisfies the
relation,

$\displaystyle
R_{\mathbf{yx}}(i,j)=R_{\mathbf{y}\bm{\tau}}(i,j)+\left[\varepsilon_{1}p_{ij}-\varepsilon_{2}d(p_{0j}-p_{ij})\right],$
(53)

where $\varepsilon_{1}$ and $\varepsilon_{2}$ are given by

$\displaystyle\varepsilon_{1}$ $\displaystyle=\sqrt{\frac{2}{\pi
p_{0j}}}\Gamma\left(1,\dfrac{d^{2}}{2p_{0j}}\right)-\frac{d}{\sqrt{\pi
p^{2}_{0j}}}\left(\Gamma\left(\dfrac{1}{2},\dfrac{d^{2}}{2p_{0j}}\right)-\sqrt{{\pi}}\right),$
(54) $\displaystyle\varepsilon_{2}$
$\displaystyle=-\frac{1}{p_{0j}}\operatorname{erf}\left(\frac{d}{\sqrt{2p_{0j}}}\right).$

In particular, for $i=j$, the relation in (53) yields

$\displaystyle R_{\mathbf{yx}}(i,i)$
$\displaystyle=R_{\mathbf{y}\bm{\tau}}(i,i)+\sqrt{\frac{2p_{0i}}{\pi}}\leavevmode\nobreak\
\Gamma\left(1,\frac{d^{2}}{2p_{0i}}\right)$ (55)
$\displaystyle-\frac{d}{\sqrt{\pi}}\Gamma\left(\frac{1}{2},\frac{d^{2}}{2p_{0i}}\right)+d.$

###### Proof.

Suppose $w_{i}$ and $w_{j}$ are the $i$-th and the $j$-th entries of
$\mathbf{w}$ ($i\neq j$) with
$\mathbb{E}\\{w_{i}\\}=\mathbb{E}\\{w_{j}\\}=-d$, and that
$p_{0i}=\bm{P}(i,i)$, $p_{0j}=\bm{P}(j,j)$, and $p_{ij}=\bm{P}(i,j)$, where
$\bm{P}$ denotes the covariance matrix of $\mathbf{w}$. Consider the quantized
random variables $y_{i}=g(w_{i})$ and $y_{j}=g(w_{j})$, where $g(.)$ denotes a
non-linear transformation function. Therefore, the cross-correlation function
between $w_{i}$ and $y_{j}$ can be obtained as

$\displaystyle R_{\mathbf{yw}}(i,j)$
$\displaystyle=\kappa\int^{\infty}_{-\infty}\int^{\infty}_{-\infty}w_{i}g(w_{j})e^{\lambda(d)}dw_{i}dw_{j},$
(56)

where $\kappa$ and $\lambda(d)$ are defined in (7) and (8), respectively. We
first calculate the integral in (56) with respect to $w_{i}$, i.e.,

$\displaystyle R_{\mathbf{yw}}(i,j)$
$\displaystyle=\frac{e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}}{2\pi\sqrt{p_{0i}p_{0j}-p_{ij}^{2}}}\int^{\infty}_{-\infty}g(w_{j})e^{\frac{2d(p_{0i}-p_{ij})w_{j}+w^{2}_{j}p_{0i}}{-2(p_{0i}p_{0j}-p_{ij}^{2})}}$
(57)
$\displaystyle\int^{\infty}_{-\infty}w_{i}e^{\frac{2d(p_{0j}-p_{ij})w_{i}+w^{2}_{i}p_{0j}-2p_{ij}w_{i}w_{j}}{-2(p_{0i}p_{0j}-p_{ij}^{2})}}dw_{i}dw_{j}$
$\displaystyle=\left[\varepsilon_{1}p_{ij}-\varepsilon_{2}d(p_{0j}-p_{ij})\right],$

where $\varepsilon_{1}$ and $\varepsilon_{2}$ are given as

$\displaystyle\varepsilon_{1}(j)$ $\displaystyle=\frac{1}{\sqrt{2\pi
p^{3}_{0j}}}\int^{\infty}_{-\infty}w_{j}g(w_{j})e^{\frac{-(w_{j}+d)^{2}}{2p_{0j}}}dw_{j},$
(58) $\displaystyle\varepsilon_{2}(j)$ $\displaystyle=\frac{1}{\sqrt{2\pi
p^{3}_{0j}}}\int^{\infty}_{-\infty}g(w_{j})e^{\frac{-(w_{j}+d)^{2}}{2p_{0j}}}dw_{j}.$

A detailed proof of the results in (57) and (58) is presented in Appendix B.
Based on (58), it can be seen that the values of $\varepsilon_{1}$ and
$\varepsilon_{2}$ are dependent on the entry index number $j$. As a result,
the modified Bussgang law for the non-stationary input signal can be presented
as:

$\displaystyle R_{\mathbf{yx}}(i,j)-R_{\mathbf{y}\bm{\tau}}(i,j)=$
$\displaystyle\left[(\varepsilon_{1}(j)+d\varepsilon_{2}(j))\right.$ (59)
$\displaystyle\left.\left(R_{\mathbf{x}}(i,j)+\bm{\Sigma}(i,j)\right)-d\varepsilon_{2}(j)p_{0j}\right].$

If the non-linear function $g(.)$ is the sign function, which is the case in
one-bit quantization, $\varepsilon_{1}$ and $\varepsilon_{2}$ can be obtained
using similar steps as presented in Part i@ of our work. However, in this
case, $p_{0j}$ should be utilized in lieu of $p_{0}$. For $i=j$, Eq. (59)
boils down to (55), a proof of which is presented in Appendix C. Based on
(59), the cross-correlation matrix between the input and the output one-bit
data is computed, where $\\{p_{0j}\\}$ is obtained by (18) and $\\{p_{ij}\\}$
can be either recovered using (30), (40) or (44). Note that the cross-
correlation matrix between the threshold vector $\bm{\tau}$ and the output
vector $\mathbf{y}$ can be estimated via a _sample cross-correlation matrix_ :

$\bm{R}_{\mathbf{y}\bm{\tau}}\approx\frac{1}{N_{\mathbf{x}}}\sum_{k=1}^{N_{\mathbf{x}}}\mathbf{y}(k)\bm{\tau}(k)^{\mathrm{H}}.$
(60)

∎

###  IX-B A Numerical Investigation of the Modified Bussgang Law

We now examine the modified Bussgang law for the non-stationary input signals
by comparing its recovery results with the true cross-correlation values
between the input signal and one-bit quantized data. In all experiments, the
input signal settings are the same as Section IV-C. The time-varying threshold
settings are as follows: (a) PA: $d=0.5$ and $\bm{\Sigma}=0.2\bm{I}$, (b)
Gauss-Legendre: $d=0.3$ and $\bm{\Sigma}=0.1\bm{I}$, (c) Monte-Carlo: $d=0.3$
and $\bm{\Sigma}=0.1\bm{I}$ , where $\bm{I}$ denotes the identity matrix.

The true cross-correlation between the input signal and the one-bit sampled
data and the estimated cross-correlation values obtained using our approach
are shown in Fig. 11, for $i=2$ and $j$ as a random sequence of length $13$.
Our results appear to confirm the possibility of recovering the cross-
correlation values from one-bit sampled data with time-varying thresholds by
employing any of the three recovery methods (PA, Gauss-Legendre method and
Monte-Carlo integration).

##  X Conclusion

We studied a generalization of the modified arcsine law discussed in Part i@
of our work through Padé approximations, Gauss-Legendre quadrature approach,
and Monte-Carlo integration, to cases where the input signal is assumed to be
non-stationary. The numerical results present the efficacy of all three
approaches in the covariance matrix recovery. Moreover, a modified Bussgang
law was established for the one-bit sampling of non-stationary input signals
with time-varying thresholds.

##  Appendix A Proof of The Convexity of $\Psi(p_{ij})$ in (39)

Since $\log(\cdot)$ is a strictly increasing function, we can analyze the
criterion
$\Psi_{m}(p_{ij})=\left(R_{\mathbf{y}}(i,j)-J_{n}(p_{ij})\right)^{2}$ to show
the convexity of $\Psi(p_{ij})$. The derivative of $\Psi_{m}(p_{ij})$ with
respect to $p_{ij}$ is computed as

$\Psi_{m}^{\prime}(p_{ij})=-2\left(R_{\mathbf{y}}(i,j)-J_{n}(p_{ij})\right)J_{n}^{\prime}(p_{ij}),$
(61)

where $J_{n}$ is the approximated version of (4) using the Gauss-Legendre
quadrature with the following close-form formula:

$\displaystyle
J_{n}(p_{ij})=\chi(p_{ij})\Bigg{(}\sqrt{p_{0i}p_{0j}-p_{ij}^{2}}$
$\displaystyle\left(\pi+2\sin^{-1}\left(\frac{p_{ij}}{\sqrt{p_{0i}p_{0j}}}\right)\right)$
(62) $\displaystyle\frac{\pi}{4}I\Bigg{)}-1=\chi(p_{ij})T(p_{ij})-1,$

where $\chi$ is defined in (12), and $I$ is given by

$I=\sum_{\varepsilon=1}^{N_{q}}\omega_{\varepsilon}\sqrt{\frac{\pi}{\beta_{n}}}\left(\frac{\alpha_{n}}{\beta_{n}}\right)\left(\frac{1}{2}-Q\left(\frac{\alpha_{n}}{\sqrt{2\beta_{n}}}\right)\right)e^{\frac{\alpha_{n}^{2}}{4\beta_{n}}}.$
(63)

Based on (62) and (63), $J_{n}^{\prime}(p_{ij})$ can be formulated as

$\displaystyle J_{n}^{\prime}(p_{ij})$
$\displaystyle=\chi(p_{ij})\left(2-\frac{p_{ij}\left(\pi+2\sin^{-1}\left(\frac{p_{ij}}{\sqrt{p_{0i}p_{0j}}}\right)\right)}{\sqrt{p_{0i}p_{0j}-p_{ij}^{2}}}\right)$
(64) $\displaystyle+\chi(p_{ij})\left(\frac{\partial
I}{\partial\alpha_{n}}\frac{\partial\alpha_{n}}{\partial
p_{ij}}+\frac{\partial I}{\partial\beta_{n}}\frac{\partial\beta_{n}}{\partial
p_{ij}}\right)$ $\displaystyle+\frac{\partial\chi}{\partial p_{ij}}T(p_{ij}),$

where $\frac{\partial\alpha_{n}}{\partial p_{ij}}$ and
$\frac{\partial\beta_{n}}{\partial p_{ij}}$ are given according to (5). As can
be seen in (61), (62), and (64), analyzing the convexity of $\Psi_{m}(p_{ij})$
depends on the parameters $d$, $p_{0i}$, $p_{0j}$, $N_{q}$, and
$\\{\theta_{\varepsilon}\\}$, which indicates that the analysis is restricted
to the case where the mentioned parameters are known; i.e. the parameters must
be specified for the covariance matrix recovery. Generally speaking, based on
(61), (62) and (64), $\Psi_{m}(p_{ij})$ is convex when
$J_{n}^{\prime}(p_{ij})>0$, or equivalently when $J_{n}(p_{ij})$ is a strictly
increasing function in the feasible region of $p_{ij}$; i.e. $-p_{m}\leq
p_{ij}\leq p_{m}$ where $p_{m}=\operatorname{min}\\{[p_{0i},p_{0j}]\\}$. As a
result, $\Psi_{m}^{\prime}(p_{ij})=0$ has only one solution which is the value
of $p_{ij}$ that satisfies $R_{\mathbf{y}}(i,j)=J_{n}(p_{ij})$. Therefore, the
convexity of $\Psi_{m}(p_{ij})$ can be easily concluded in light of (61). For
instance, one may easily verify that the selected parameters for the recovery
of the input covariance matrix in Section V-C makes $J_{n}(p_{ij})$ a strictly
increasing function, and thus, $\Psi_{m}(p_{ij})$ a convex function.

##  Appendix B Proof of The Modified Bussgang Law Formula

Note that

$\displaystyle R_{\mathbf{yw}}(i,j)$
$\displaystyle=\frac{e^{\frac{-d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}}{2\pi\sqrt{p_{0i}p_{0j}-p_{ij}^{2}}}\int^{\infty}_{-\infty}g(w_{j})e^{\frac{2d(p_{0i}-p_{ij})w_{j}+w^{2}_{j}p_{0i}}{-2(p_{0i}p_{0j}-p_{ij}^{2})}}$
(65)
$\displaystyle\int^{\infty}_{-\infty}w_{i}e^{\frac{2d(p_{0j}-p_{ij})w_{i}+w^{2}_{i}p_{0j}-2p_{ij}w_{i}w_{j}}{-2(p_{0i}p_{0j}-p_{ij}^{2})}}dw_{i}dw_{j},$

where the inner integral and the outer integral are called $\mathcal{L}_{1}$
and $\mathcal{L}_{2}$, respectively. The inner integral may be evaluated as

$\displaystyle\mathcal{L}_{1}$
$\displaystyle=\int^{\infty}_{-\infty}w_{i}e^{\frac{2d(p_{0j}-p_{ij})w_{i}+w^{2}_{i}p_{0j}-2p_{ij}w_{i}w_{j}}{-2(p_{0i}p_{0j}-p^{2}_{ij})}}dw_{i}$
(66)
$\displaystyle=e^{\frac{\left(p_{0j}d-p_{ij}\left(w_{j}+d\right)\right)^{2}}{2p_{0j}(p_{0i}p_{0j}-p^{2}_{ij})}}\int^{\infty}_{-\infty}w_{i}e^{-\frac{\left(w_{i}+\left(d-\frac{p_{ij}}{p_{0j}}\left(w_{j}+d\right)\right)\right)^{2}}{2(p_{0j}-\frac{p^{2}_{ij}}{p_{0j}})}}dw_{i}$
$\displaystyle=e^{\frac{\left(p_{0j}d-p_{ij}\left(w_{j}+d\right)\right)^{2}}{2p_{0j}(p_{0i}p_{0j}-p^{2}_{ij})}}\sqrt{2\pi\left(p_{0j}-\frac{p^{2}_{ij}}{p_{0j}}\right)}\times\cdots$
$\displaystyle\left(\frac{p_{ij}}{p_{0j}}\left(w_{j}+d\right)-d\right).$

Next, the outer integral is evaluated as

$\displaystyle\mathcal{L}_{2}$
$\displaystyle=\sqrt{2\pi\left(p_{0j}-\frac{p^{2}_{ij}}{p_{0j}}\right)}e^{\frac{p^{2}_{0j}d^{2}+p^{2}_{ij}d^{2}-2dp_{0j}p_{ij}}{2p_{0j}(p_{0i}p_{0j}-p^{2}_{ij})}}\times\cdots$
(67)
$\displaystyle\int^{\infty}_{-\infty}g(w_{j})\left(\frac{p_{ij}}{p_{0j}}\left(w_{j}+d\right)-d\right)e^{\frac{\alpha_{ij}\left(w^{2}_{j}+2dw_{j}\right)}{\beta_{ij}}}dw_{j},$

where $\alpha_{ij}=p_{0i}p_{0j}-p^{2}_{ij}$ and
$\beta_{ij}=-2p_{0j}\alpha_{ij}$. The above integral can be simplified as
below:

$\displaystyle\mathcal{L}_{2}$
$\displaystyle=\sqrt{2\pi\left((p_{0j}-\frac{p^{2}_{ij}}{p_{0j}})\right)}e^{\frac{d^{2}(p_{0i}+p_{0j}-2p_{ij})}{2(p_{0i}p_{0j}-p_{ij}^{2})}}\times\cdots$
(68)
$\displaystyle\int^{\infty}_{-\infty}g(w_{j})\left(\frac{p_{ij}}{p_{0j}}\left(w_{j}+d\right)-d\right)e^{-\frac{\left(w_{j}+d\right)^{2}}{2p_{0j}}}dw_{j}.$

When the values from (66) and (68) are inserted in Eq. (65), we have the final
form of the modified Bussgang law, i.e.

$\displaystyle
R_{\mathbf{yw}}(i,j)=\varepsilon_{1}p_{ij}-\varepsilon_{2}d(p_{0j}-p_{ij}),$
(69)

where $\varepsilon_{1}$ and $\varepsilon_{2}$ are obtained as below:

$\displaystyle\varepsilon_{1}$ $\displaystyle=\frac{1}{\sqrt{2\pi
p^{3}_{0j}}}\int^{\infty}_{-\infty}w_{j}g(w_{j})e^{-\frac{(w_{j}+d)^{2}}{2p_{0j}}}dw_{j},$
(70) $\displaystyle\varepsilon_{2}$ $\displaystyle=\frac{1}{\sqrt{2\pi
p^{3}_{0j}}}\int^{\infty}_{-\infty}g(w_{j})e^{-\frac{(w_{j}+d)^{2}}{2p_{0j}}}dw_{j}.$

##  Appendix C Proof of Equation (55)

To prove (55), we may formulate $\mathbb{E}\\{w_{i}y_{i}\\}$ as

$\displaystyle\mathbb{E}\\{w_{i}y_{i}\\}$
$\displaystyle=\int_{-\infty}^{\infty}w_{i}f(w_{i})p(w_{i})\,dw_{i}$ (71)
$\displaystyle=\int_{0}^{\infty}w_{i}p(w_{i})\,dw_{i}-\int_{-\infty}^{0}w_{i}p(w_{i})\,dw_{i}$
$\displaystyle=\int_{0}^{\infty}w_{i}\left(p(w_{i})+p(-w_{i})\right)\,dw_{i},$

where $p(w_{i})=\left(\sqrt{2\pi
p_{0i}}\right)^{-1}e^{\frac{-\left(w_{i}+d\right)^{2}}{2p_{0i}}}$. Let
$I_{1}=\int_{0}^{\infty}w_{i}e^{\frac{-\left(w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}$
and
$I_{2}=\int_{0}^{\infty}w_{i}e^{\frac{-\left(-w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}$.
Then, we have

$\displaystyle I_{1}$
$\displaystyle=\int_{0}^{\infty}w_{i}e^{\frac{-\left(w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}=\int_{0}^{\infty}\left(w_{i}+d\right)e^{\frac{-\left(w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}$
(72)
$\displaystyle-d\int_{0}^{\infty}e^{\frac{-\left(w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}=d\sqrt{\frac{\pi
p_{0i}}{2}}\left(-1+\operatorname{erf}\left(\frac{d}{\sqrt{2p_{0i}}}\right)\right)$
$\displaystyle+p_{0i}e^{-\frac{d^{2}}{2p_{0i}}},$ $\displaystyle I_{2}$
$\displaystyle=\int_{0}^{\infty}w_{i}e^{\frac{-\left(-w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}=\int_{0}^{\infty}\left(w_{i}-d\right)e^{\frac{-\left(-w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}$
$\displaystyle+d\int_{0}^{\infty}e^{\frac{-\left(-w_{i}+d\right)^{2}}{2p_{0i}}}\,dw_{i}=d\sqrt{\frac{\pi
p_{0i}}{2}}\left(1+\operatorname{erf}\left(\frac{d}{\sqrt{2p_{0i}}}\right)\right)$
$\displaystyle+p_{0i}e^{-\frac{d^{2}}{2p_{0i}}}.$

Therefore, based on (72), a simplified form of the integration in (71) can be
proposed as

$\mathbb{E}\\{w_{i}y_{i}\\}=d\operatorname{erf}\left(\frac{d}{\sqrt{2p_{0i}}}\right)+\sqrt{\frac{2p_{0i}}{\pi}}e^{-\frac{d^{2}}{2p_{0i}}}.$
(73)

Since $\Gamma\left(1,x\right)=e^{-x}$ and
$\Gamma\left(\frac{1}{2},x\right)=\sqrt{\pi}\left(1-\operatorname{erf}\left(\sqrt{x}\right)\right)$,
(73) is identical to (55).

## References

  * [1] Arian Eamaz, Farhang Yeganegi, and Mojtaba Soltanalian,  “Modified arcsine law for one-bit sampled stationary signals with time-varying thresholds,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 5459–5463.
  * [2] Arian Eamaz, Farhang Yeganegi, and Mojtaba Soltanalian,  “Covariance recovery for one-bit sampled data with time-varying sampling thresholds—Part I: Stationary signals,”  IEEE Transactions on Signal Processing, Submitted, 2022.
  * [3] Ricardo Antonio Maronna,  “Robust m-estimators of multivariate location and scatter,”  The Annals of Statistics, pp. 51–67, 1976.
  * [4] Amir Weiss,  “Blind direction-of-arrival estimation in acoustic vector-sensor arrays via tensor decomposition and kullback-leibler divergence covariance fitting,”  IEEE Transactions on Signal Processing, vol. 69, pp. 531–545, 2020\.
  * [5] Arindam Bose, Shahin Khobahi, and Mojtaba Soltanalian,  “Efficient waveform covariance matrix design and antenna selection for mimo radar,”  Signal Processing, vol. 183, pp. 107985, 2021.
  * [6] Ziyang Cheng, Zishu He, Bin Liao, and Min Fang,  “Mimo radar waveform design with papr and similarity constraints,”  IEEE Transactions on Signal Processing, vol. 66, no. 4, pp. 968–981, 2017.
  * [7] Ilya Bekkerman and Joseph Tabrikian,  “Target detection and localization using mimo radars and sonars,”  IEEE Transactions on Signal Processing, vol. 54, no. 10, pp. 3873–3883, 2006.
  * [8] Luzhou Xu, Jian Li, and Petre Stoica,  “Target detection and parameter estimation for mimo radar systems,”  IEEE Transactions on Aerospace and Electronic Systems, vol. 44, no. 3, pp. 927–939, 2008.
  * [9] Hamza Djelouat, Markus Leinonen, Lucas Ribeiro, and Markku Juntti,  “Joint user identification and channel estimation via exploiting spatial channel covariance in mmtc,”  IEEE Wireless Communications Letters, vol. 10, no. 4, pp. 887–891, 2021.
  * [10] Mélanie Mahot, Frédéric Pascal, Philippe Forster, and Jean-Philippe Ovarlez,  “Asymptotic properties of robust complex covariance matrix estimates,”  IEEE Transactions on Signal Processing, vol. 61, no. 13, pp. 3348–3356, 2013.
  * [11] Texas Instruments,  “Analog embedded processing,”  Active Filter Design Techniques, Literature Number SLOA088 Texas Instruments, excerpted from OP-Amps for Everyone, Literature Number SLOD006A.
  * [12] Amine Mezghani and A Lee Swindlehurst,  “Blind estimation of sparse broadband massive MIMO channels with ideal and one-bit ADCs,”  IEEE Transactions on Signal Processing, vol. 66, no. 11, pp. 2972–2983, 2018.
  * [13] Aria Ameri, Jian Li, and Mojtaba Soltanalian,  “One-bit radar processing and estimation with time-varying sampling thresholds,”  in 2018 IEEE 10th Sensor Array and Multichannel Signal Processing Workshop (SAM). IEEE, 2018, pp. 208–212.
  * [14] Saeid Sedighi, Bhavani Shankar, Mojtaba Soltanalian, and Björn Ottersten,  “One-bit DoA estimation via sparse linear arrays,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 9135–9139.
  * [15] J. H. Van Vleck and David Middleton,  “The spectrum of clipped noise,”  Proceedings of the IEEE, vol. 54, no. 1, pp. 2–19, 1966.
  * [16] Giovanni Jacovitti, Alessandro Neri, and Roberto Cusani,  “Methods for estimating the autocorrelation function of complex Gaussian stationary processes,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 35, no. 8, pp. 1126–1138, 1987.
  * [17] Giovanni Jacovitti and Alessandro Neri,  “Estimation of the autocorrelation function of complex Gaussian stationary processes by amplitude clipped signals,”  IEEE Transactions on Information Theory, vol. 40, no. 1, pp. 239–245, 1994.
  * [18] Ofer Bar-Shalom and Anthony J Weiss,  “DoA estimation using one-bit quantized measurements,”  IEEE Transactions on Aerospace and Electronic Systems, vol. 38, no. 3, pp. 868–884, 2002.
  * [19] Chun-Lin Liu and PP Vaidyanathan,  “One-bit sparse array DoA estimation,”  in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2017, pp. 3126–3130.
  * [20] Julian Jakob Bussgang,  “Crosscorrelation functions of amplitude-distorted gaussian signals,”  1952\.
  * [21] Cheng Qian and Jian Li,  “ADMM for harmonic retrieval from one-bit sampling with time-varying thresholds,”  in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2017, pp. 3699–3703.
  * [22] Christopher Gianelli, Luzhou Xu, Jian Li, and Petre Stoica,  “One-bit compressive sampling with time-varying thresholds for sparse parameter estimation,”  in 2016 IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM). IEEE, 2016, pp. 1–5.
  * [23] Shahin Khobahi and Mojtaba Soltanalian,  “Model-based deep learning for one-bit compressive sensing,”  IEEE Transactions on Signal Processing, vol. 68, pp. 5292–5307, 2020.
  * [24] Shahin Khobahi and Mojtaba Soltanalian,  “Signal recovery from 1-bit quantized noisy samples via adaptive thresholding,”  in 2018 52nd Asilomar Conference on Signals, Systems, and Computers. IEEE, 2018, pp. 1757–1761.
  * [25] Pu Wang, Jian Li, Milutin Pajovic, Petros T Boufounos, and Philip V Orlik,  “On angular-domain channel estimation for one-bit massive MIMO systems with fixed and time-varying thresholds,”  in 51st Asilomar Conference on Signals, Systems, and Computers. IEEE, 2017, pp. 1056–1060.
  * [26] Feng Xi, Yijian Xiang, Shengyao Chen, and Arye Nehorai,  “Gridless parameter estimation for one-bit MIMO radar with time-varying thresholds,”  IEEE Transactions on Signal Processing, vol. 68, pp. 1048–1063, 2020.
  * [27] Simon Haykin and Liang Li,  “Nonlinear adaptive prediction of nonstationary signals,”  IEEE Transactions on signal processing, vol. 43, no. 2, pp. 526–535, 1995.
  * [28] Yue Hu, Fucai Li, Hongguang Li, and Chengliang Liu,  “An enhanced empirical wavelet transform for noisy and non-stationary signal processing,”  Digital Signal Processing, vol. 60, pp. 220–229, 2017.
  * [29] Yves Grenier,  “Time-dependent ARMA modeling of nonstationary signals,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 31, no. 4, pp. 899–911, 1983.
  * [30] Katherine Tsai,  A non-convex framework for structured non-stationary covariance recovery theory and application,  Ph.D. thesis, 2020.
  * [31] Robert F Engle, Olivier Ledoit, and Michael Wolf,  “Large dynamic covariance matrices,”  Journal of Business & Economic Statistics, vol. 37, no. 2, pp. 363–375, 2019.
  * [32] Emily B Fox and David B Dunson,  “Bayesian nonparametric covariance regression,”  The Journal of Machine Learning Research, vol. 16, no. 1, pp. 2501–2542, 2015.
  * [33] Tim Bollerslev, Robert F Engle, and Daniel B Nelson,  “Arch models,”  Handbook of econometrics, vol. 4, pp. 2959–3038, 1994.
  * [34] Maurice George Kendall, Alan Stuart, and J Keith Ord,  Kendall’s advanced theory of statistics,  Oxford University Press, Inc., 1987.
  * [35] Monson H Hayes,  Statistical digital signal processing and modeling,  John Wiley & Sons, 2009.
  * [36] JL Basdevant,  “The Padé approximation and its physical applications,”  Fortschritte der Physik, vol. 20, no. 5, pp. 283–331, 1972.
  * [37] Claude Brezinski and J Van Iseghem,  “A taste of Padé approximation,”  Acta numerica, vol. 4, pp. 53–103, 1995.
  * [38] Pedro Gonnet, Stefan Guttel, and Lloyd N Trefethen,  “Robust Padé approximation via SVD,”  SIAM review, vol. 55, no. 1, pp. 101–117, 2013.
  * [39] Marco Chiani, Davide Dardari, and Marvin K Simon,  “New exponential bounds and approximations for the computation of error probability in fading channels,”  IEEE Transactions on Wireless Communications, vol. 2, no. 4, pp. 840–845, 2003.
  * [40] Wenyu Sun and Ya-Xiang Yuan,  Optimization theory and methods: nonlinear programming, vol. 1,  Springer Science & Business Media, 2006.