\section{Introduction}\label{s_intro}
Consider the input/output linear system
\begin{equation}\label{system_td}\left\{\begin{array}{l} x_k=\mathrm q x_{k-1}+\mathrm w u_{k-1}~~k=1,\dots,K\\
y_k=c x_k+n_k\\
\end{array}\right.
\end{equation}
with $K\in\mathds{N}$ (possibly tending to infinity), $u_k\in\{0,1\}$ for $k=0,\dots, K-1$, $x_k\in\mathds{R}$ for $k=0,\dots,K$, $y_k,n_k\in\mathds{R}$ for $k=1,\dots,K$,  $\mathrm q,\mathrm w,c\in\mathds{R}$, and $\mathrm q\in (0,1)$ to preserve stability. Our aim is to recover the binary input $u_k$, in an online fashion, given the output $y_k$ corrupted by a noise $n_k$. To this purpose, we retrieve a low-complexity algorithm introduced in \cite{s09} and discussed in \cite{s10,s11}, and we propose a comprehensive theoretical analysis of its performance. As a result of the analysis, we will be able to evaluate the performance  as a function of the system's parameters.

The digital signal reconstruction problem is a paradigm in data transmissions, where signals arising from finite alphabets are sent over noisy continuous channels, and in hybrid frameworks, where digital and analog signals have to be merged in the same system. 
In \cite{s09}, a particular instance of model \eqref{system_td} was derived as time discretization of a convolution system and the input estimation described  as a deconvolution problem. The same can be achieved for model \eqref{system_td}:  if we consider the system
\begin{equation}\label{system}\left\{\begin{array}{l} x'(t)=ax(t)+b u(t)~~~~ t \in [0,T]\\
y(t)=cx(t)+n(t)~~~~x(0)=x_0\\
u(t), x(t), y(t),~ a, b, c \in \mathds{R}, a<0
\end{array}\right.
\end{equation} 
we have \begin{equation}\label{convolution} x(t)= e^{t a}x_0+ b\int_0^t e^{a(t-s)} u(s) \mathrm{d} s.\end{equation}
Given \begin{equation}\label{bin_input}u(t)=\sum_{k=0}^{K-1} u_k \mathds{1}_{[k\tau,(k+1)\tau[}(t),~~~u_k \in \mathcal{U}=\{0,1\},~\tau>0\end{equation} 
we can discretize in the following way: by defining
\begin{equation}\begin{split}
     x_k&:=x(k \tau )~~~ \text{ for } k=0,1,\dots, K=T/\tau\\
     \mathrm q&:=e^{\tau a}\in (0,1)\\
     \mathrm w&:=b\frac{1-e^{\tau a}}{-a}=-\frac{b}{a}(1-\mathrm q)
   \end{split}\end{equation}
we obtain
\begin{equation}\label{enc}\begin{split} x_k &=\mathrm q^k x_0+b\mathrm q^k \int_0^{k\tau}e^{-as}\sum_{h=0}^{K-1} u_h
\mathds{1}_{[h\tau,(h+1)\tau[}(s)\mathrm{d} s\\
&=\mathrm q^k x_0+b\mathrm q^k\sum_{h=0}^{k-1} u_h\int_{h\tau}^{(h+1)\tau}e^{-as}\mathrm{ d}s\\
&=\mathrm q^k x_0+\frac{b}{-a}\mathrm q^k\sum_{h=0}^{k-1} u_h e^{-a(h+1)\tau}(1-e^{a\tau})\\
&=\mathrm q^k x_0+\mathrm w\sum_{h=0}^{k-1} u_h \mathrm q^{k-1-h}\\
\end{split}\end{equation}
 from which we have the recursive formula  
\begin{equation}\label{rec_encoding}   x_k=\mathrm q x_{k-1}+\mathrm w u_{k-1}. \end{equation} 

In system \eqref{system}, recovering $u(t)$ basically consists  in the inversion of the convolution integral $ y(t)= c e^{t a}x_0+ cb\int_0^t e^{a(t-s)} u(s) \mathrm{d}+n(t)$ (where $n(t)$ represents an additive noise), which is a long-standing mathematical  ill-posed problem: small observation errors  may produce defective solutions. Several estimation approaches have been studied in the last fifty years and the literature on
deconvolution is widespread: we refer the reader to early papers \cite{tik63, tik77}  and to later  \cite{ary78,spa96,sta02,fag02},
which show also some possible applications in geophysics, astronomy, image processing and biomedical systems. For more references, see
\cite{SoThesis}.

Most of known deconvolution methods exploit the regularity of the input function to provide good estimations. This work  instead is a contribution for deconvolution in case of discontinuous input functions.

Considering a  binary alphabet, which has been chosen mainly  to keep the analysis straightforward, is  consistent with many applications: the output of several digital devices, such as computers and detection devices \cite{s10}, are binary. Nevertheless, the algorithm and the analysis presented in this paper could be generalized to larger alphabets with no much effort. 

In \cite{s09},  low-complexity decoding algorithms were introduced,  derived from the optimal BCJR \cite{BCJR} algorithm, and applied  to perform the 
deconvolution of the system \eqref{system} with $a=0$ and $b=c=1$. In this work, we apply the simplest of those algorithms, the so-called One State Algorithm (OSA for short) to the system \eqref{system_td}. We then describe the performance in terms of  Mean Square Error (MSE)
for long-time
transmissions, through a probabilistic analysis arising from the Markovian behavior of
the algorithm. The scheme of the analysis is the same proposed in \cite{s09}, but leads to completely 
different scenarios: while for  $a=0$, $b=c=1$  standard ergodic theorems for denumerable Markov Processes were sufficient to compute
the
MSE, in the present case the denumerable model does not proved the expected results, and more sophisticated
arguments are used, arising from Markov Processes, Iterated Random Functions (IRF for short) and sequences of probability measures.

The paper is organized as follows. In Section \ref{s_prob} we complete the description of the system, giving some observations and probabilistic assumptions; in Sections \ref{s_alg} and \ref{s_sim}, we present our algorithm and some simulations. The core of the paper is the performance analysis provided in Section \ref{s_analysis}. Finally we propose some concluding observations.  Notice that Sections \ref{s_prob} and \ref{s_alg} mainly retrieve the
model presented in \cite{s09} and \cite{s10}.

\subsection{Notation}
We use the following notation  throughout the paper: $\Prob$ indicates a discrete probability, while $P(\cdot,\cdot)$ is the
transition probability kernel of a Markov Process; $\mathbb{E}$ is the stochastic mean. 
$\mathcal{B}(\mathcal{S})$ indicates the Borel $\sigma$-field of a space $\mathcal{S}$. Given a bounded measurable function $v$ defined on a space $\mathcal{S}$, $Pv(x)=\int_{\mathcal{S}}v(y)P(x,\mathrm{d}y)$. For every measure $\mu$ on  $(\mathcal{S},\mathcal{B}(\mathcal{S}))$ and $F \in \mathcal{B}(\mathcal{S})$,  $\mu P(F)=\int_{\ms}P(x,F)\mu(\mathrm{d}x)$.
The complementary error function $\text{erfc}$ is defined as $\text{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_x^{+\infty} e^{-t^2}\mathrm{d}t$, $x\in\mathds{R}$; the indicator function $\mathds{1}_A(x)$ is equal to one if $x\in A$, and zero otherwise. Moreover we often use the following acronyms: OSA for One State Algorithm,  MSE for Mean Square Error, IRF for for Iterated Random Functions, MAP for Maximum a Posteriori.


\section{Problem Statement }\label{s_prob}
Let us develop a deeper understanding of the problem and  specify some assumptions.

First notice that, given $ x_k=\mathrm q x_{k-1}+\mathrm w u_{k-1}$, we have 
\begin{equation*} x_k=\mathrm q^k x_0+\mathrm w\sum_{h=0}^{k-1}u_{k-h-1}\mathrm q^h\end{equation*}
which shows that each $x_k$ is determined by the initial state $x_0$ and by a binary  polynomial  of degree $k-1$ in $\mathrm q$.
From now onwards, let $x_0=0$, so that, for any $k=0,1,\dots,K$,
\begin{equation} x_k\in\mathcal{X}=\mathrm w \left \{ \sum_{h=0}^{K} \mu_h \mathrm q^h,~\mu_h\in\{0,1\} \right\}.\end{equation}

Moreover, let us introduce some prior probabilistic information:

\emph{Assumption 1:} the additive noise $n_k$  is white Gaussian, that is,  $n_1,\dots, n_K$  are realizations of independent Gaussian
random
variables $N_1,\dots, N_K$, with null mean and variance $\sigma^2$.

\emph{Assumption 2:} the binary inputs $u_0,\dots, u_{K-1} $ are realizations of independent Bernoulli random variables $U_0, \dots
U_{k-1}$ with parameter $\frac{1}{2}$.

Input and noise are also supposed to be mutually independent. Under these assumptions the 
system  can be rewritten in probabilistic terms as follows  (capital letters are used instead of
small letters to indicate random quantities): for $ k=1,\dots, K$,
\begin{equation}
\left\{\begin{array}{l}
U_{k-1}\sim \text{Ber}\left(1/2\right)\\ 
N_k\sim \mathcal{N}(0,\sigma^2)\\ 
X_k=\mathrm q X_{k-1}+\mathrm w U_{k-1}~~(X_0=0)\\
Y_k=c X_k+N_k.\\
\end{array}\right.\end{equation}

While Assumption 1 is realistic in physical terms, Assumption 2 is less motivated by applications, where  source bits are often not independent (for example, they may be governed by a Markov Chain). We have however imposed it for simplicity of treatment, although extensions to more sophisticated prior distributions do not require much effort. Similarly, we have chosen to propose a one-dimensional problem to make the analysis more readable, while the structure would be almost the same also for multidimensional problems (see \cite{s10}).

Given this setting, we aim at providing a method to decode the bit $u_{k-1}$ at each time step $k=1,\dots, K$, based on the current measurement  $y_k$ and  the probabilistic properties of the system. In order to perform this online recovery, the algorithm is allowed to store a few information (just one real value) about the state $x_k$.

\section{Binary Input Reconstruction }\label{s_alg}
The One State Algorithm (OSA for short) introduced in \cite{s09} fits the requirements described in the previous section.

The OSA is a suboptimal version of the Bahl, Cocke, Jelinek, and Raviv algorithm (most known as BCJR \cite{BCJR}), a prominent decoding
algorithm  used for convolutional codes. The BCJR  performs a maximum a posteriori estimation (MAP) of the input bit sequence by
evaluating the probabilities of the states of the encoder, through a forward and a backward recursion; it is optimal in the sense that
it minimizes the Mean Square Error:
\begin{equation} \text{MSE}(\RU,\RUE):=\frac{1}{K}\sum_{k=0}^{K-1}\mathds{E}(U_k-\UE_{k})^2 \end{equation} where
$\RU=(U_0,\dots, U_{K-1})$ and $\RUE=(\UE_0,\dots, \UE_{K-1})$ is the estimated input sequence. In the binary case, this is equivalent
to
\begin{equation*} \text{MSE}(\RU,\RUE)=\frac{1}{K}\sum_{k=0}^{K-1}\Prob(U_k\neq\hat{U}_{k}).\end{equation*}
As shown in \cite{s09}, the BCJR can be adapted to the binary input deconvolution problem with optimal results, but with complexity drawbacks when
the transmission is long. This motivated the introduction of the OSA, a BCJR-based method that consists only in a  forward recursion and
that stores only one state at each iteration step. More precisely, the OSA pattern is as follows:
\begin{enumerate}
    \item Initialization of the state estimate $\xe_0$;
    \item For $k=1,\dots,K$: given $y_k \in \mathds{R}$ and $\xe_{k-1}$,
     \begin{equation}\label{OSApattern}\begin{split}&\ue_{k-1}=
\left\{\begin{array}{l}
0\text{ if }  |y_k-c\mathrm q\xe_{k-1}|\leq |y_k-(c\mathrm q\xe_{k-1}+c\mathrm w)|  \\
1\text{ otherwise.}
\end{array}
\right.\\
&\xe_{k}=\mathrm q\xe_{k-1}+\mathrm w \ue_{k-1}.\end{split}\end{equation}
\end{enumerate}
Typically, we assume to know $x_0$, so that we can initialize correctly by $\xe_0=x_0$. This point will be discussed later in Section \ref{s_analysis}. 
Given the probabilistic setting previously introduced, the OSA can also be written as:
\begin{equation}\left\{\begin{array}{l}
\UE_{k-1}=\mathds{1}_{(c\mathrm q \XE_{k-1}+\frac{c\mathrm w}{2},+\infty)}(Y_k)\\
\XE_k=\mathrm q \XE_{k-1}+\mathrm w \UE_{k-1}.\\
\end{array}\right.\end{equation}
While the BCJR estimates the probabilities of all the possible states at each step, the OSA  individuates the most likely state and
assumes it to be the correct one; on the basis of this state estimate, it decides on the current bit. The decoding is performed with
a MAP decision on the current bit, which in our probabilistic setting (Bernoulli input and Gaussian noise) reduces to the comparison between two Euclidean
distances.

The OSA is suboptimal, but presents two main good properties: (a) it is low-complexity,  both for number of computations and storage
locations; (b) it is causal, that is, it uses only the past and the present information to decode the current bit. Therefore, (a) it can
be applied to our case in which the number of states is (not countably) infinite and (b) it can be used online, making unnecessary the
complete transmission before starting deconvolution, this feature being fundamental to study long time transmissions.

In \cite{s09}, we introduced other causal algorithms: the Causal BCJR, which is a version of BCJR performing only the forward recursion,
and
the Two States Algorithm (TSA), which works basically the same as OSA, but estimates the two best states at each step with their probabilities of being the correct ones. The TSA is then oriented to soft decoding; in the differentiation
case $a=0$, $b=c=1$, it was proved to have similar performance to Causal BCJR and better than the OSA. Nevertheless,  neither the
Causal BCJR nor the TSA are efficient for system \eqref{system_td}. The first one, in fact, presents  complexity drawbacks due to the 
number of states. The second one, instead, has performance too much similar to the OSA, in spite of higher complexity: owing to the structure of the
state space $\mathcal{X}$,  the two best states turn out to be very close to each other,  which does not enhance the information
provided by the OSA.

\subsection{Similar algorithms}
We notice that our setting and decoding procedure \eqref{OSApattern} are very similar to the Decision Feedback Equalizer introduced to mitigate the effects of channels' intersymbol inference (ISI, see \cite{pulf_thesis} for a complete review). As in our case, the model considered in channel equalization is a linear system with digital input, and the goal is the input recovery for equalization purposes. Various methods have been proposed in literature and much effort has been addressed  towards complexity reduction (see, e.g., \cite{eyu88, due89, wil92, que07}). Typically, complexity is reduced by collecting information only from fixed time blocks, which is also our attempt; more precisely we consider only the current measurement and an estimation of the previous state, that is, ``minimal'' blocks of length one, but extensions to larger time blocks are possible to improve the performance.

The recovery techniques in the cited works present many analogies with ours. For example, the method introduced in \cite{que07}, if restricted to minimal blocks, differs from ours only for the introduction of a prior distribution on the state $x_k$.

Nevertheless, an outstanding difference lies in the model: channel equalization exploits the input estimate to provide a feedback equalizer to the system, while our final aim is just the input recovery.

Given the several connections, in our future work we will study possible implementations of our low-complexity algorithms, derived from BCJR, for channel equalization and propose more detailed comparisons.



\section{Simulations}\label{s_sim}

We now show  a few simulations' results, obtained by 2000 Monte Carlo Runs of 320 bit transmissions. 
\begin{figure}
   \begin{center}
   \includegraphics[width=8cm]{out_simu-crop.pdf}
   \caption{Mean Square Error in function of the Signal-To-Noise Ratio $\frac{c^2\mathrm w^2}{\sigma^2}$ (in dB); $b=c=1$, $a=-2,-1,-0.5,-0.3$}
\includegraphics[width=8cm]{out_simu_zoom-crop.pdf}\label{figF1}
   \caption{A zoom that highlights the gain obtained by decreasing $a$.}
   \label{figF1b}\end{center}
 \end{figure} We consider the system derived from \eqref{system}, with $\tau=1$, $\mathrm q=e^{ a}$, $\mathrm w=-\frac{b}{a}(1-\mathrm q)$.  We show the behavior of the
MSE with respect to $\frac{c^2\mathrm w^2}{\sigma^2}$, that can be interpreted as the Signal-To-Noise-Ratio (SNR for short) of the
transmission: since for each $k$, the transmitted signal is $cx_k\in\{c\mathrm q
x_{k-1},c\mathrm q x_{k-1}+c\mathrm w\}$ then $c^2\mathrm w^2$ is proportional to the signal power. As expected, the MSE tends to
zero when the SNR is large, while for small SNR tends to  $\frac{1}{2}$.

If we fix $b=c=1$ and let $a$ vary, we obtain a slight gain (that is, a lower MSE curve) by decreasing $a$, as shown in Figures
\ref{figF1}-\ref{figF1b}.
In other terms, more stable systems are preferable. This phenomenon will be retrieved in Section \ref{s_analysis}.


\section{Analysis of the Algorithm}\label{s_analysis}
For simplicity, in the next we will assume $\mathrm w>0$, the analysis in the case $\mathrm w<0$ being analogous.

The goal of this section is the analytic evaluation of the Mean Square Error for the One State Algorithm, in case of long-time
transmissions.

The analysis starts from the definition of the following Markov Process:
\begin{equation}\left\{\begin{array}{l}D_k=\XE_k-X_k=\mathrm q D_{k-1}+\mathrm w(\UE_{k-1}+U_{k-1})~~~~~k=1,2,\dots\\
D_0=\alpha\\
\end{array}\right.\end{equation}
For any $k=1,2,\dots$, if  $D_{k-1}=z$ then $D_{k}\in\{\mathrm q z,\mathrm q z+\mathrm w,\mathrm q z -\mathrm w\}$, and the transition probabilities are:

\begin{equation}\label{tpk}\begin{split}
P(z,\mathrm q z+\mathrm w)&=\Prob(\UE_k=1,U_k=0|D_k=z)\\
&=\frac{1}{4}\text{erfc}\left(\frac{c\mathrm q z+c\mathrm w/2}{\sigma\sqrt{2}}
\right)\\
P(z,\mathrm q z-\mathrm w)&=\Prob(\UE_k=0,U_k=1|D_k=z)\\
&=\frac{1}{4}\text{erfc}\left(\frac{-c\mathrm q z+c\mathrm w/2}{\sigma\sqrt{2}}
\right)\\
P(z,\mathrm q z)&=1-P(z,\mathrm q z+\mathrm w)-P(z,\mathrm q z-\mathrm w).\\
\end{split}\end{equation}


Since for any $k\in\mathds{N}$, $D_k\in\left \{\mathrm w  \sum_{h=0}^{k-1}
\mu_h \mathrm q^h,~\mu_h\in\{-1,0,1\} \right\}$, if we fix $D_0=0$, $(D_k)_{k\in\mathds{N}}$  is a Markov Process on the denumerable state space \begin{equation}\label{denumstsp}
\left \{ \mathrm w \sum_{h=0}^{\infty} \mu_h \mathrm q^h,~\mu_h\in\{-1,0,1\} \text{ such that } (\mu_h)_{h\in\mathds{N}} \text{ is definitely
null}\right\}.\end{equation}
By definitely null, we mean that for any $D_k$ the coefficients $\mu_h$ with $h\geq k$ are null. This set is denumerable since any $D_k$ can be seen as the ternary representation of  a non-negative integer. Notice that  fixing $D_0=0$ just means that $x_0$ is known.

The key point of the analysis is that, for large $k$, the MSE of the OSA can be computed using the ergodic properties of the Markov
Process $(D_k)_{k\in\mathds{N}}$; more precisely, we require the existence of a stationary distribution. In the next, we will propose two different ways to study the stationary distribution: the first one does not depend on the initial state $D_0\in \ms$ (where $\ms$ is compact set which will be defined shortly), but requires some contractive properties; the second one is valid even in the non-contractive case, but depends on the initial state. 

For both methods, the presented setting is not still
adequate to study the possible stationary distributions, since the states of $(D_k)_{k\in\mathds{N}}$ are transient: when the process visits a
state, then it leaves it definitely (except for a negligible set of states that have a periodic ternary representation, for example $0$, $\pm\mathrm w/(1-\mathrm q)$ ); moreover, the
process is not irreducible since there is no reciprocal communication between the states (see
\cite{s09} for more details). Thus we conclude that  no hypotheses are fulfilled to apply the
standard ergodic results for denumerable Markov
Processes  (see \cite{s09}). In other terms, if a stationary distribution exists, it does not concentrate on single states.

This suggests to consider  $(D_k)_{k\in\mathds{N}}$ on a non-denumerable state space. In particular, we can
 extend  \eqref{denumstsp} to
\begin{equation}\label{extstsp}\left \{ \mathrm w \sum_{h=0}^{\infty} \mu_h \mathrm q^h,~\mu_h\in\{-1,0,1\}\right\}.\end{equation}
It is interesting to notice that if $\mathrm q\geq\frac{1}{3}$, then the set \eqref{extstsp} coincides with the closed interval
$\left[-\frac{\mathrm w}{1-\mathrm q},+\frac{\mathrm w}{1-\mathrm q}\right]$, while if $\mathrm q<\frac{1}{3}$ it is a Cantor set included in
$\left[-\frac{\mathrm w}{1-\mathrm q},+\frac{\mathrm w}{1-\mathrm q}\right]$ (for a proof of this fact, see \cite{s10}).
Let us then consider as state space
\begin{equation*} \ms=\left[-\frac{\mathrm w}{1-\mathrm q},+\frac{\mathrm w}{1-\mathrm q}\right]\end{equation*}
and study the ergodic properties of the Markov Process $(D_k)_{k\in\mathds{N}}$ on $\ms$. 

Before continuing the analysis, let us introduce some rigorous notions that will be used in the next. 

Let $\mathcal{B}(\ms)$ be the Borel $\sigma$-field of $\ms$.  We call  \textit{transition probability kernel} (see, e.g.,
\cite[Section 3.4.1]{mey93})  an application $P:\ms \times \mathcal{B}(\ms)\to[0,1]$ such that\\ 
(i) for each $F \in\mathcal{B}(\ms)$, $P(\cdot, F)$ is a non-negative measurable function;\\
(ii) for each $x \in \ms$, $P(x, \cdot)$ is a probability measure  on $(\ms,\mathcal{B}(\ms))$.

Given a bounded measurable function $v$ on $\ms$, we denote by $Pv$ the bounded measurable function 
defined as
\begin{equation*}
Pv(x)=\int_{\ms}v(y)P(x,\mathrm{d}y).\end{equation*}
Furthermore, let $\mu$ be a measure on  $(\ms,\mathcal{B}(\ms))$: we define the measure $\mu P$ by 
\begin{equation*}\mu P(F)=\int_{\ms}P(x,F)\mu(\mathrm{d}x)\;\;\;\;\;F \in \mathcal{B}(\ms).\end{equation*}
We finally define the $n$-th power of the transition kernel $P$ by $P^1(x,F)=P(x,F)$ and 
$P^n(x,F)=\int_{\ms}P(x,\mathrm{d}y)P^{n-1}(y,F)$. It is easy to see that $P^n(x,F)$ are transition kernels, too.

At this point, we can make explicit the relationship between the MSE and $(D_k)_{k\in\mathds{N}}$. In the next, we will always consider $D_0=0$, if not differently specified.  Given the transition probability kernel $P$
of $(D_k)_{k\in\mathds{N}}$, defined by \eqref{tpk}, we have
\begin{equation}\begin{split}\label{acheserveilMP}\text{MSE}&(\RU,\RUE) =\frac{1}{K}\sum_{k=0}^{K-1}\Prob(\UE_k\neq U_k)=\\
&=\frac{1}{K}\sum_{k=0}^{K-1}\int_{\mathcal{D}}\Prob(\UE_k\neq
U_k|D_k=z)P^k(\alpha,\mathrm d z)=\frac{1}{K}\sum_{k=0}^{K-1}P^k g(\alpha)
\end{split}\end{equation}
where \begin{equation}\label{def_g}g(z)=\Prob(\UE_k\neq U_k|D_k=z)\end{equation}
and $D_0=\alpha$ is any initial state in $\ms$. Therefore $P^k g$ (and in particular its behavior for large $k$) will be the object of our further
analysis.

In the sequel, we will distinguish two main scenarios: when $(D_k)_{k\in\mathds{N}}$ has some
\emph{contractive} properties and when it has not. In the first scenario, we can exploit the theory of Iterated Random Functions to
prove that $P^k g$  converges, while in the second one we will use known results of convergence of probability measures.


\subsection{Contractive case}
Let $ l\text{-Lip}(\ms)$ be  the set of all the Lipschitz functions with Lipschitz constant equal to $l$ on $\ms$. We define the
Kantorovich (or Wasserstein) distance $d_W$  
between probability measures (see \cite[Section 2.1, Example 3.2.2]{rac91}) as
\begin{equation}\label{Kantorovich} d_W(\mu,\nu)=\sup_{f\in 1\text{-Lip}(\ms)}\bigg|\int_{\ms} f\mathrm d (\mu- \nu) \bigg|.\end{equation}

We can prove the following
\begin{theorem}\label{ERGT}
If $\frac{c^2\mathrm w^2}{\sigma^2}>4$ and $\mathrm q <\frac{1}{3+\sqrt{\frac{2}{e\pi}}}$  or if $\frac{c^2\mathrm w^2}{\sigma^2}\leq 4$ and
$\mathrm q\leq \frac{1}{1+\sqrt{\frac{2}{e\pi}}}$,
then,
\begin{equation}\label{ERG} \lim_{K\to\infty}\text{MSE}(\RU,\RUE)=\int_{\mathcal{D}}g \mathrm  d\mu~~~\text{ for any } D_0\in \ms\end{equation}
where $\mu$
is the unique probability measure such that  $\sup_{x\in\ms} d_W(P^k(x,\cdot),\mu(\cdot)) \stackrel{k\to\infty}{\longrightarrow}0$,
 $P$ being the kernel of $(D_k)_{k\in \mathds{N}}$.\end{theorem}

Notice that $g(z)$ is time-invariant (i.e., does not depend on $k$) and can be analytically computed. In fact, given $D_{k-1}=z$,
$D_k=\mathrm q z$ if and only if $\UE_k=U_k$,  $D_k=\mathrm q z+\mathrm w$ if and only if $\UE_k=1$ and $U_k=0$, $D_k=\mathrm q z-\mathrm w$ if and only if $\UE_k=0$ and
$U_k=1$ and
\begin{equation} g(z)=P(z,\mathrm q z+\mathrm w)+P(z,\mathrm q z -\mathrm w).\end{equation}
Furthermore, the probability measure $\mu$ can be numerically evaluated.

Recall  that, as already noticed, $\frac{c^2\mathrm w^2}{\sigma^2}$ can be interpreted as the Signal-To-Noise-Ratio. Moreover, the bounds on $\mathrm q$ can be interpreted as the necessity of a stronger stability for convergence.



In order to prove the theorem, let us introduce some elements from the Iterated Random Functions theory.
\subsubsection{Iterated Random Functions}\label{IRF}
Let $(\ms, d)$ be a complete metric space and $\mis$ be a measurable space. Consider a measurable function $w: \ms \times
\mis \to \ms$ and for each fixed $s \in \mis$, $w_s(x):=w(x,s)$, $x\in\ms$. Let $(I_k)_{k\in\mathds{N}}$ be a stochastic sequence in $\mis$ such
that $I_0, I_1, \dots$ are independent, identically distributed. Then, the set $\{w_{I_k}(x),~k\in\mathds{N}\}$ is a family of random
functions. 
The systems obtained by iterating such random functions, called Iterated Random Functions (IRF), are studied for diverse
purposes: for example, IRF with contractive properties are used to construct fractal sets, see \cite{hut81, dia99}. More interesting for
our study is the exploitation of IRF to study Markov Processes. Given an IRF and a starting state $x\in\ms$, we can  define the induced
Markov Process  $(Z_k(x))_{k\in\mathds{N}}$ as 
\begin{equation} Z_k(x):=w_{I_{k-1}}\circ w_{I_{k-2}}\circ w_{I_{k-3}}\circ\cdots\circ w_{I_{0}}(x)~~~~(k\geq 1)\end{equation}
and analyze its asymptotic behavior through the properties of $w_{I_k}(x), k\in\mathds{N}$. It has been proved that if the $w_{I_k}(x)$ have
\emph{some} contractive properties, the transition probability kernel of $Z_n(x)$  converges to a probability measure, unique for all
the initial states
$x\in \ms$. The required contractive properties may be slightly different: \cite{dia99} studied the case of Lipschitz functions
$w_{I_k}(x)$ ``contracting on average'', while similar results have been obtained by \cite{sten01} without the continuity requirement on
$w_{I_k}(x)$,  by \cite{stein99} for ``locally contractive'' functions, and  by \cite{jar01} for ``non-separating on average''
functions. A useful survey on the argument has been recently proposed by \cite{ios09}.

Let us show how to exploit the IRF theory in our framework.

The evolution of $(D_k)_{k\in\mathds{N}}$ can be modeled by IRF. We consider the complete metric space $\ms=\left[-\frac{\mathrm w}{1-\mathrm q},+\frac{\mathrm w}{1-\mathrm q}\right]$ naturally endowed with the
Euclidean metric $d$ of $\mathds{R}$, the measurable space $\mis=\{0,1\}\times \mathds{R}$ and the stochastic process $$I_k=(U_k,N_{k+1}),~~k\in\mathds{N}$$ on
$\mis$, and we define the random function
\begin{equation} w_{I_k}(x)=\mathrm q x+\mathrm w \mathds{1}_{(c\mathrm q x+c\mathrm w\left(\frac{1}{2}- U_k\right),+\infty)}(N_{k+1})-\mathrm w U_k,~~x\in\ms\end{equation} that describes
the
dynamics of $(D_k)_{k\in\mathds{N}}$. The key result for our purpose is the following theorem (here stated for compact spaces), which does not
require continuity:

\textbf{Stenflo's Theorem} \cite[Theorem 1]{sten01}
Suppose that there exists a constant $l<1$ such that \begin{equation}\label{av_conv} \mathbb{E}[d(w_{I_0}(x),w_{I_0}(y))]\leq l~d(x,y)\end{equation} for all
$x,y\in\ms$, $(\ms,d)$ being a compact metric space.
Then there exist a unique probability measure $\mu$ and a positive constant $\gamma_{\ms}$ such that
\begin{equation} \sup_{x\in\ms} d_W(P^n(x,\cdot),\mu(\cdot))\leq \frac{\gamma_{\ms}}{1-l}l^n~~~~n\geq 0\end{equation}
where $P^n(x,\cdot)$ is the $n$-step transition probability kernel of the Markov Process $Z_n(x)$.\\


Now, Theorem \ref{ERGT} can be proved by applying the Stenflo's Theorem.

\subsubsection{ Proof of Theorem \ref{ERGT}}
Let us analyze the condition (\ref{av_conv}). Consider $x,y\in\ms$ with $x>y$  (recall that $\mathrm q>0, \mathrm w>0$). Let $ H=H(x,y,I_0)$ and
$\mathcal{I}_u$ respectively be defined as
\begin{equation*}\begin{split}& H:= \mathds{1}_{\left(c\mathrm q y+c\mathrm w\left(\frac{1}{2}- U_0\right),c\mathrm q x+c\mathrm w\left(\frac{1}{2}-
U_0\right)\right)}(N_1)\\ 
&\mathcal{I}_u:=\frac{1}{\sqrt{2\pi}\sigma}\int_{c\mathrm q y+c\mathrm w\left(\frac{1}{2}-
u\right)}^{c\mathrm q x+c\mathrm w\left(\frac{1}{2}- u\right)}e^{-\frac{n^2}{2\sigma^2}}\mathrm d n\\  
&=\frac{1}{2}\text{erfc}\left(c\frac{\mathrm q y+\mathrm w\left(\frac{1}{2}-
u\right)}{\sigma\sqrt{2}}\right)-\frac{1}{2}\text{erfc}\left(c\frac{\mathrm q x+\mathrm w \left(\frac{1}{2}-
u\right)}{\sigma\sqrt{2}}\right).\\
\end{split}\end{equation*}
Hence,
\begin{equation}\label{my_contr}
\begin{split}
&\mathbb{E}\left[|w_{(U_0,N_1)}(x)-w_{(U_0,N_1)}(y)|\right]
=\mathbb{E}\left[ |\mathrm q(x-y)-\mathrm w H | \right]\\
&=\sum_{u\in\{0,1\}}\Prob(U_0=u) \frac{1}{\sqrt{2\pi}\sigma}\int_{\mathds{R}}e^{-\frac{n^2}{2\sigma^2}}|\mathrm q(x-y)-\mathrm w H | \mathrm{d} n\\
&=\frac{1}{2}\sum_{u\in\{0,1\}}\frac{1}{\sqrt{2\pi}\sigma}\int_{\mathds{R}}e^{-\frac{n^2}{2\sigma^2}}|\mathrm q(x-y)-\mathrm w H |\mathrm{d} n\\
&=\frac{1}{2}\sum_{u\in\{0,1\}}|\mathrm q(x-y)-\mathrm w|\mathcal{I}_u+\mathrm q(x-y)(1-\mathcal{I}_u).\\
\end{split}\end{equation} 
If $\mathrm q(x-y)>\mathrm w$, then $\mathbb{E}\left[|w_{(U_0,N_1)}(x)-w_{(U_0,N_1)}(y)|\right]<\mathrm q(x-y)$ and the contraction would be proved with $l=\mathrm q$.
This is never the case when $\mathrm q<\frac{1}{3}$ and  $|x-y|\leq 2\frac{\mathrm w}{1-\mathrm q}< \frac{\mathrm w}{\mathrm q}$.

Let us then consider  $\mathrm q(x-y)<\mathrm w$. We can write
\begin{equation}\begin{split}\mathbb{E}&\left[|w_{(U_0,N_1)}(x)-w_{(U_0,N_1)}(y)|\right]\\
&=\frac{1}{2}\sum_{u\in\{0,1\}}(\mathrm w-\mathrm q(x-y))\mathcal{I}_u+\mathrm q(x-y)(1-\mathcal{I}_u)\\
&\leq\frac{1}{2}\sum_{u\in\{0,1\}}\mathrm w\mathcal{I}_u+\mathrm q(x-y).\\
\end{split}\end{equation}
The last expression is obtained by neglecting $-\sum_{u\in\{0,1\}}~ \mathrm q$ $(x-y)\mathcal{I}_u$, which is the sum of two second degree
terms in $(x-y)$, since, by the integral mean value theorem,
 \begin{equation} \mathcal{I}_u= \frac{1}{\sqrt{2\pi}\sigma}c\mathrm q (x-y)e^{-\frac{n_0^2}{2\sigma^2}}\end{equation} for some $n_0\in\left(c\mathrm q
y+c\mathrm w\left(\frac{1}{2}-u\right),c\mathrm q x+c\mathrm w\left(\frac{1}{2}-u\right) \right)$, $(n_0 \neq 0).$ 
The remaining terms are of order one.
 Notice  that
\begin{equation}\frac{1}{2}\sum_{u\in\{0,1\}}\mathrm w\mathcal{I}_u+\mathrm q(x-y)=F(x)-F(y) \end{equation}
where \begin{equation*} F(x)=\mathrm q x-\frac{\mathrm w}{4}\text{erfc}\left(\frac{c\mathrm q x+c\frac{\mathrm w}{2}}{\sigma\sqrt{2}}\right)-\frac{\mathrm w}{4}\text{erfc}\left(\frac{c\mathrm q
x-c\frac{\mathrm w}{2}}{\sigma\sqrt{2}}\right).\end{equation*}
Therefore, the thesis is achieved if $F(x)$ is a contraction; since $F(x)$ is differentiable and monotone increasing, its Lipschitz
constant is the maximum of its first derivative:\begin{equation}\label{F1}\begin{split}&
F'(x)=\mathrm q+\\&\frac{c\mathrm w\mathrm q}{2\sigma\sqrt{2\pi}}\left[\exp\left(-\frac{\left(c\mathrm q x+c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)
+\exp\left(-\frac{\left(c\mathrm q x-c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)\right].\\\end{split}\end{equation}
In order to find the maximum of $F'(x)$, let us compute $F''(x)$:
\begin{equation*}\begin{split}&F''(x)=\\ &=-\frac{c\mathrm w\mathrm q}{2\sigma\sqrt{2\pi}}\frac{2\left(c\mathrm q x+c\frac{\mathrm w}{2}\right)c\mathrm q}{2\sigma^2}
\exp\left(-\frac{\left(c\mathrm q x+c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)+\\
&~~~-\frac{c\mathrm w\mathrm q}{2\sigma\sqrt{2\pi}}\frac{2\left(c\mathrm q x-c\frac{\mathrm w}{2}\right)c\mathrm q}{2\sigma^2} \exp\left(-\frac{\left(c\mathrm q
x-c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)\\
\end{split}\end{equation*}
which is null for  $x$ satisfying:
\begin{equation}\label{exp_eq} \left(\mathrm q x+\frac{\mathrm w}{2}\right)\exp\left(-\frac{c^2\mathrm q\mathrm w x}{\sigma^2}\right)+ \left(\mathrm q x-\frac{\mathrm w}{2}\right)=0 \end{equation}
a solution of which is $x=0$. Now, considering that $F'(x)$ is determined by a mixture of two Gaussians, two cases may occur: (a) $x=0$
is the maximum of $F'(x)$; (b) $x=0$ is a minimum for $F'(x)$ and there are two symmetric maxima ($F''(x)$ is an even function) at
$x_0\in (0,\frac{\mathrm w}{1-\mathrm q}]$ and $-x_0$, but $x_0$ cannot be analytically computed from the exponential equation (\ref{exp_eq}). Let us study the sign of $F''(x)$ for $x\to 0$ in order to determine the nature of the point $x=0$ for $F'(x)$. Notice that 
\begin{equation}\begin{split}&F''(x)>0 \Leftrightarrow 
-\left(\mathrm q x+\frac{\mathrm w}{2}\right)\exp\left(-\frac{c^2\mathrm q\mathrm w x}{\sigma^2}\right)- \left(\mathrm q x-\frac{\mathrm w}{2}\right)>0.\\
\end{split}\end{equation}
Moreover, if $x\to 0$, 
$\exp\left(-\frac{c^2\mathrm q\mathrm w x}{\sigma^2}\right)\sim 1 - \frac{c^2\mathrm q\mathrm w x}{\sigma^2} $ and 
\begin{equation}\begin{split}&-\left(\mathrm q x+\frac{\mathrm w}{2}\right)\exp\left(-\frac{c^2\mathrm q\mathrm w x}{\sigma^2}\right)- \left(\mathrm q x-\frac{\mathrm w}{2}\right)\sim-2\mathrm q x+\frac{c^2\mathrm q\mathrm w^2 x}{2\sigma^2}.\\
\end{split}\end{equation}
Finally, if $ \frac{c^2\mathrm w^2}{\sigma^2}>4$
\begin{equation} -2\mathrm q x+\frac{c^2\mathrm q\mathrm w^2 x}{2\sigma^2}\to 0^+ \text{ for } x\to 0^+\end{equation}
\begin{equation} -2\mathrm q x+\frac{c^2\mathrm q\mathrm w^2 x}{2\sigma^2} \to 0^- \text{ for } x\to 0^-.\end{equation}
In conclusion, $x=0$ is a maximum point for $F'(x)$ if and only if  $\frac{c^2\mathrm w^2}{\sigma^2}<4$, that is, only for large noise.


Let us now study $\frac{c^2\mathrm w^2}{\sigma^2}>4$ ($x=0$ is a minimum point) and let us state conditions that make $F(x)$  contractive.
In particular, consider
$x>0$ and $\sigma^2$ close to zero: by (\ref{exp_eq}),  $|x-\frac{\mathrm w}{2\mathrm q}|$ tends to zero more quickly than $\sigma^2$, hence
$\exp\left(-\frac{\left(c\mathrm q x-c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)$ tends to one and the maximum of $F'(x)$ (see \eqref{F1}) may
assume very large values.

More in general, we observe that the points $x=\pm\frac{\mathrm w}{2\mathrm q}$ are tricky as they are the unique points where the OSA fails: for
these values, the error probability given by (\ref{tpk}) is $\frac{1}{2}$, no matter which is the noise variance. This ``singular''
phenomenon is more evident when the noise is small; in terms of $F(x)$, it causes large variations (then the loss of the contractivity)
in a neighborhood of the point  $\pm\frac{\mathrm w}{2\mathrm q}$, the radius of the neighborhood being larger for smaller $\sigma^2$.

Let us set in the case  $\mathrm q<\frac{1}{3}$, so that $\pm\frac{\mathrm w}{2\mathrm q}\notin \ms$. Under this
assumption, for any $x\in\ms$,
\begin{equation}\begin{split}
&\exp\left(-\frac{\left(c\mathrm q
x+c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)<\exp\left(-\frac{\left(-c\mathrm q\frac{\mathrm w}{1-\mathrm q}+c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)\\
&\exp\left(-\frac{\left(c\mathrm q x-c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)<\exp\left(-\frac{\left(c\mathrm q
\frac{\mathrm w}{1-\mathrm q}-c\frac{\mathrm w}{2}\right)^2}{2\sigma^2}\right)\\
\end{split}\end{equation}
hence
\begin{equation} F'(x)\leq\mathrm q+\frac{c\mathrm w\mathrm q}{\sigma\sqrt{2\pi}}
\exp\left(-\frac{ c^2 \mathrm w^2\left(\frac{1-3\mathrm q}{2(1-\mathrm q)}\right)^2}{2\sigma^2}\right).\end{equation}
As for $t\geq 0$, $\max te^{-t^2}=\frac{1}{\sqrt{2e}}$,
\begin{equation*} \frac{\mathrm q}{\sqrt{\pi}}\frac{2}{1-3\mathrm q}\frac{1}{\sqrt{2e}}<1 \Longrightarrow F'(x)<1\end{equation*}


In conclusion a sufficient condition for average contractivity is

\begin{equation*}\mathrm q<\frac{1}{3+\sqrt{\frac{2}{e\pi}}}.\end{equation*} 


Let us now study the case $\frac{c^2\mathrm w^2}{\sigma^2}<4$ ($x=0$ is the maximum point):
\begin{equation}
\begin{split}
F'(x)&\leq F'(0)=\mathrm q+2\frac{c\mathrm w\mathrm q}{2\sigma \sqrt{2\pi}}\exp\left(-\frac{c^2\mathrm w^2}{8\sigma^2}\right) \\
       &\leq\mathrm q+2\frac{\mathrm q}{ \sqrt{2e\pi}}\\  
\end{split}
\end{equation}
then $F'(x)<1$ when $$\mathrm q<\frac{1}{1+\sqrt{\frac{2}{e\pi}}}.$$



In conclusion, we have stated that  if $\frac{c^2\mathrm w^2}{\sigma^2}>4$ and $\mathrm q<\frac{1}{3+\sqrt{\frac{2}{e\pi}}}$  or if 
$\frac{c^2\mathrm w^2}{\sigma^2}\leq 4$ and
$\mathrm q\leq \frac{1}{1+\sqrt{\frac{2}{e\pi}}}$, then the hypotheses of Stenflo's Theorem are fulfilled. 

Now, let us prove the
convergence of the Mean Square Error. 

Since $g\in L_g$-Lip$(\ms)$ where $L_g=\max_{z\in\ms}
|g'(z)|$ and 
\begin{equation} |g'(z)|=\frac{c\mathrm q}{2\sigma\sqrt{2\pi}}\left|-e^{-\frac{(c\mathrm q z+c\mathrm w)^2}{2\sigma^2}}
+e^{-\frac{(-c\mathrm q z+c\mathrm w)^2}{2\sigma^2}}\right|\leq\frac{c\mathrm q}{\sigma\sqrt{2\pi}} \end{equation}
then $L_g\leq \frac{c\mathrm q}{\sigma\sqrt{2\pi}}$ is finite. Since for any $L>0$ 
 \begin{equation*}
 \begin{split}
&\sup_{f\in L\text{-Lip}(\ms)}\bigg|\int_{\ms} f \mathrm d (\mu- \nu)\bigg| = \sup_{f\in L\text{-Lip}(\ms)}L\bigg|\int_{\ms}\frac{1}{L}
f \mathrm d
(\mu- \nu)\bigg|\\ &\leq  \sup_{f\in 1\text{-Lip}(\ms)}L\bigg|\int_{\ms} f \mathrm d (\mu- \nu)\bigg|= L~ d_W(\mu,\nu)   
 \end{split}
\end{equation*}
we have 
\begin{equation}
\begin{split}
  \sup_{x\in\ms}\bigg| P^k g(x)&-\int_{\ms} g\mathrm d \mu\bigg| =\sup_{x\in\ms}\bigg | \int_{\ms} g(z)P^k(x,\mathrm d z)-\int
_{\ms} g\mathrm d
\mu\bigg|
\\ &\leq \sup_{x\in\ms} \sup_{f\in L_g\text{-Lip}}\bigg|\int_{\ms} f(z)P^k(x,\mathrm d z)-\int_{\ms} f\mathrm d \mu\bigg|\\  
&= \sup_{x\in\ms} L_g d_W(P^k(x,\cdot), \mu(\cdot))\stackrel{k\to\infty}{\longrightarrow}0.\\  
\end{split}
\end{equation}
The convergence is then assured also for the Cesàro sum, for any initial state $D_0=\alpha\in\ms$:
\begin{equation} \frac{1}{K}\sum_{k=0}^{K-1}P^k g(\alpha)\stackrel{K\to\infty}{\longrightarrow}\int_{\ms} g\mathrm d \mu~~~~\forall \alpha\in\ms.\end{equation}
\qed

Notice that the initial state does not affect the convergence value if it is contained in $\ms$, but $\ms$ has been obtained by fixing $D_0=0$: this seems not coherent. However, even if we consider $D_0=\alpha\notin \ms$, given the dynamics of $D_k$ ($\alpha$ is multiplied by $\mathrm q$ at each step), $\ms$ turns out to be the ``limit'' state space, and with high probability $D_k$ enters $\ms$ for some finite $k$, so it makes sense to reduce to $\ms$. Further details are here omitted for brevity, but one must be convinced that considering the initial error lying in a compact set centered at $0$ is a suitable choice.

\subsection{Non contractive case}
If the hypotheses of Theorem \ref{ERGT} are not fulfilled, we can prove the following
\begin{theorem}\label{NonContrTheor}
For any initial state $D_0=\alpha$, there exists a unique probability measure $\phi(\alpha,\cdot)$ such that
\begin{equation}
    \label{NonContr}
     \lim_{K\to\infty}\text{MSE}(\RU,\RUE)=\phi g(\alpha)
  \end{equation}
where $\phi g(\alpha)=\int g(z)\phi(\alpha,\mathrm{d} z) $.
\end{theorem}
We recall that although this result holds also for the contractive case, the IRF argument is
preferable in that case since the convergence value is independent from $\alpha$.
\subsubsection{Proof of Theorem \ref{NonContrTheor}}
Let $A_{x,n}$ the set of the points that $D_k$ can reach in $n$ steps starting from $x$, i.e., $A_{x,n}=\{\mathrm q^n x +\mathrm w \sum_{i=0}^{n-1}\alpha_i\mathrm q^i,~\alpha_i\in\{-1,0,1\} \}$
\begin{lemma}\label{equilip}
For any $f\in l_f \text{-Lip}(\ms) $, there exists a positive constant $\mathrm{M}_f$ such that $P^n f \in  \frac{\mathrm{M}_f}{1-\mathrm q}\text{-Lip}(\ms)$ for any $n\in\mathds{N}$.
\end{lemma}
\proof
We have
\begin{equation}
\max_{x\in\ms}\max_{y\in A_{x,1}} \left|\frac{\mathrm{d}}{\mathrm{d}x}P(x,y)\right| \leq \frac{c\mathrm q}{\sigma\sqrt{2\pi}}
\end{equation}
then the three functions $P(x,\mathrm q x)$, $P(x, \mathrm q x+\mathrm w)$ and $P(x,\mathrm q x-\mathrm w)$ are Lipschitz with constant $l:=\frac{c\mathrm q}{\sigma\sqrt{2\pi}}$
in $\ms$. For any $x,x_0 \in \ms$ and any $(y,y_0)\in\{(\mathrm q x,\mathrm q x_0)$, $(\mathrm q x+\mathrm w, \mathrm q x_0+\mathrm w), (\mathrm q x-\mathrm w, \mathrm q x_0-\mathrm w) \}$,
\begin{equation}
  \begin{split}
&\left|P(x,y)f(y)-P(x_0,y_0)f(y_0) \right|=\\ 
&~~~~=\left|P(x,y)f(y)\pm  P(x_0,y_0)f(y) - P(x_0,y_0)f(y_0) \right| \\  
&~~~~\leq ||f||_{\infty}|P(x,y)-P(x_0,y_0)|+ P(x_0,y_0)|f(y) -f(y_0)| \\  
&~~~~\leq ||f||_{\infty} l |x-x_0|+ P(x_0,y_0)l_f|y-y_0| \\  
&~~~~\leq \left(||f||_{\infty} l+ P(x_0,y_0)l_f \mathrm q\right) |x-x_0|.  
  \end{split}
\end{equation}
Thus, 
\begin{equation}
  \begin{split}
&\left|Pf(x)-Pf(x_0) \right|=\\&=|P(x,\mathrm q x)f(\mathrm q x)-P(x_0,\mathrm q x_0)f(\mathrm q x_0)+ P(x,\mathrm q x+\mathrm w)f(\mathrm q x+\mathrm w)\\&~~~~-P(x_0,\mathrm q x_0+\mathrm w)f(\mathrm q
x_0+\mathrm w)+P(x,\mathrm q x-\mathrm w)f(\mathrm q x-\mathrm w)\\&~~~~-P(x_0,\mathrm q x_0-\mathrm w)f(\mathrm q x_0-\mathrm w)|  \\
&\leq \left(3 ||f||_{\infty} l |x-x_0|+l_f \mathrm q\right) |x-x_0|.    
  \end{split}
\end{equation}
In conclusion, $Pf\in L_1\text{-Lip}(\ms)$  where $L_1=3 l||f||_{\infty} +l_f \mathrm q$. 

Now, given any $n\in\mathds{N}$ and $x=x_0+\delta$,
 \begin{equation}\label{Pngeneral}
  \begin{split}
&\left|P^nf(x)-P^n f(x_0) \right|= \left|\sum_{z\in A_{x,1}}P(x,z)  P^{n-1}f(z)- \sum_{z_0\in A_{x_0,1}}P(x_0,z_0)  P^{n-1}f(z_0)
\right| \\
&=\bigg|\sum_{z\in A_{x,1}}P(x,z)  P^{n-1}f(z) \pm \sum_{z_0\in A_{x_0,1}}P(x_0,z_0)P^{n-1}f(z_0+\mathrm q \delta)\\&~~~~ - \sum_{z_0\in
A_{x_0,1}}P(x_0,z_0)  P^{n-1}f(z_0) \bigg| \\
&= \bigg|\sum_{z\in A_{x,1}} \left[P(x,z) - P(x-\delta,z-\mathrm q \delta)\right]  P^{n-1}f(z)\\&~~~~+ \sum_{z_0\in A_{x_0,1}}P(x_0,z_0)\left[
P^{n-1}f(z_0+\mathrm q \delta)  -  P^{n-1}f(z_0)\right] \bigg| \\
&\leq \sum_{z\in A_{x,1}} \left|P(x,z) - P(x-\delta,z-\mathrm q \delta)\right|~||f||_{\infty}\\&~~~~+ \sum_{z_0\in A_{x_0,1}}P(x_0,z_0)\left|
P^{n-1}f(z_0+\mathrm q \delta)  -  P^{n-1}f(z_0)\right| \\
&\leq 3 l\delta ||f||_{\infty}+  \sum_{z_0\in A_{x_0,1}}P(x_0,z_0)\left| P^{n-1}f(z_0+\mathrm q \delta)  -  P^{n-1}f(z_0)\right|.
 \end{split}
\end{equation}
If $n=2$,
\begin{equation*}
\left|P^2f(x)-P^2 f(x_0) \right|\leq 3 l \delta ||f||_{\infty}+  \sum_{z_0\in A_{x_0,1}}P(x_0,z_0)L_1\mathrm q\delta\leq  \left(3 l ||f||_{\infty}+L_1\mathrm q\right)\delta 
\end{equation*}
that is,  $P^2f\in L_2\text{-Lip}(\ms)$ where $L_2=3 l ||f||_{\infty}+L_1\mathrm q$.
At this point, by iterating \eqref{Pngeneral}, we obtain that for any $n\in\mathds{N}$, $P^nf\in L_n\text{-Lip}(\ms)$ where $L_n=3 l ||f||_{\infty}+L_{n-1}\mathrm q$.
Moreover, be recursion, 
\begin{equation*}
  L_n=3 l ||f||_{\infty}(1+\mathrm q+\dots \mathrm q^{n-1})+\mathrm q^n l_f \leq \frac{\mathrm{M}_f}{1-\mathrm q}~~~~~\mathrm{M}_f:=\max\{3 l  ||f||_{\infty}, l_f \}.
\end{equation*}
\qed

Let us recall that a sequence of measures $\{\mu_n\}_{n\in\mathds{N}}$ is said to be \emph{weakly convergent} to a measure $\mu$ if $\lim_{n\to\infty}\int f(x)\mathrm{d} \mu_n=\int f(x)\mathrm{d}\mu$ for every continuous and bounded function $f$. In the next, we will denote weakly convergence by $\mu_n\xrightarrow{~w~}\mu$.

\begin{lemma}\label{exist_subs}
Let $\overline{P}_N(x,\cdot)=\frac{1}{N}\sum_{n=0}^{N-1}P^n(x,\cdot)$, $N\in\mathds{N}$. For any $x\in\ms$, there exist  a subsequence $\overline{P}_{N_j}(x,\cdot)$, $j, N_j\in\mathds{N}$, and a probability measure $\phi(x,\cdot)$ such that $\overline{P}_{N_j}(x,\cdot)\xrightarrow{~w~} \phi(x,\cdot)$. 
\end{lemma}
\proof
This is a simple consequence of Prohorov's Theorem (see, e.g., \cite[Theorem 6.1]{bill68}; in our context tightness is trivial since the space $\ms$ is compact).
\qed


\begin{lemma}\label{all_subs}
 If all the convergent subsequences of $\overline{P}_{N}(x,\cdot)$ weakly converge to the same $\phi(x,\cdot)$, then also $\overline{P}_{N}(x,\cdot)$ weakly converges to $\phi(x,\cdot)$.
\end{lemma}
\proof
Again this is a consequence of Prohorov's Theorem (see, e.g.,\cite[Theorem 2.3]{bill68} )
\qed

Given  Lemmas \ref{exist_subs} and \ref{all_subs}, to prove Theorem \ref{NonContrTheor} it is sufficient to show that all the convergent subsequences of $\overline{P}_{N}(x,\cdot)$ converge to $\phi(x,\cdot)$. 

Let us suppose that there exist a subsequence $\{M_i\}_{i\in \mathds{N}} \neq \{N_j\}_{j\in\mathds{N}}$ and a probability measure $\psi(x,\cdot)\neq \phi(x,\cdot)$ on $\ms$ such that  $\overline{P}_{M_i}(x,\cdot )\xrightarrow{~w~} \psi(x,\cdot)$.

First notice that for any $m\in\mathds{N}$, by applying the  Dominated Convergence Theorem,
\begin{equation*}
  \begin{split}
P^m \phi f(x)&=\int_{y\in\ms}P^m(x,\mathrm{d}y)\phi f(y)\\&=\int_{y\in\ms}P^m(x,\mathrm{d}y)\lim_{j\to +\infty}\int_{z\in\ms}\frac{1}{N_j}\sum_{n=0}^{N_j-1}P^n(y,\mathrm{d}z)f(z)\\
&=\lim_{j\to +\infty}\frac{1}{N_j}\sum_{n=0}^{N_j-1}P^{n+m}f(x)\\&=\lim_{j\to +\infty}\frac{1}{N_j}\left(\sum_{n=0}^{N_j-1}P^{n}f(x)+\sum_{n=N_j}^{N_j-1+m}P^{n}f(x)-\sum_{n=0}^{m-1}P^{n}f(x)\right)\\&=\phi f(x).\\
\end{split}
\end{equation*}

Similarly, exploiting the continuity of  $P^mf$, we obtain
$$\phi P^m f(x)=\phi f (x).$$
The same can be clearly said for $\psi$.

Now, for any $M_i\in\mathds{N}$
$$\overline{P}_{M_i}\phi f(x)=\phi f(x)$$ and $$\lim_{i\to\infty}\overline{P}_{M_i}\phi f(x)=\phi f(x).$$
If $f$ is $l_f$-Lip$(\ms)$, then $\overline{P}_{N_j}f(x)$ are equicontinuous by Lemma \ref{equilip}, and clearly also equibounded by $||f||_{\infty}$. Therefore, by Ascoli-Arzelà Theorem, $\phi f(x)$ is continuous and $\lim_{i\to\infty}\overline{P}_{M_i}\phi f(x)=\psi \phi f(x)$
In conclusion,

$$\phi f (x)=\psi \phi f (x).$$
Analogously, one can prove that $\psi f(x)=\phi \psi f(x)$ and, since by Dominated Convergence $\psi \phi f(x)=\phi \psi f(x), $ we  obtain $\phi f(x)=\psi f(x)$. To summarize, we have proved that, for any $x\in \ms$, there exists a unique probability measure $\phi(x,\cdot)$ such that
$$\lim_{N\to\infty}\overline{P}_N f(x)=\phi f(x)$$ for any $f\in l_f$-Lip($\ms$).

The thesis of Theorem \ref{NonContrTheor} follows by considering $f=g$.\qed 

The arguments used in this proof partially arise from the proof of \cite[Theorem 12.4.1]{mey93}.
  
\subsection{Simulations vs Theoretical Results}
The convergence values $\int_{\mathcal{D}}g \mathrm  d\mu$ and $\phi g(\alpha)$ can be numerically evaluated. \begin{figure}
   \begin{center}
   \includegraphics[width=5.7cm]{out2-crop.pdf}
\includegraphics[width=5.7cm]{out1-crop.pdf}
\includegraphics[width=5.7cm]{out05-crop.pdf}
\includegraphics[width=5.7cm]{out03-crop.pdf}
   \caption{Simulations vs Theoretical Results: MSE for $b=c=1$,
$a=-2,-1,-0.5,-0.3$.}\label{F2}
   \end{center}
 \end{figure} In Figure \ref{F2}, we show  their consistency with the simulations previously presented. Notice that for simulations we have assumed to know the initial state $x_0$, so that $D_0=0$. Since analytical results are asymptotic,
while simulations' results are obtained by averaging  transmissions of 320 bits, we intuitively conclude that the rate of
convergence is fast. 
\section{Conclusion}\label{s_concl}
In this paper, we have proposed using the One State decoding algorithm to recover the binary input of a linear system and  we have
analyzed its behavior. When the system has particular contractive properties, the analysis is based on Iterated Random Functions, while
in the non-contractive case known results from convergence of probability measures can be exploited. The theoretical results allow to predict the Mean Square
Error of the One State Algorithm for long-time transmissions, given the parameters of the system and some prior probabilistic information. Simulations and theoretical results are consistent.

The One State Algorithm could be extended to multi-dimensional problems and to the recovery of digital inputs arising from larger source alphabets and with different probabilistic distributions. Moreover, its use for problems with feedback, such as channel equalization, should be further studied.  

\section{Acknowledgements}
The author wishes to thank Prof. Fabio Fagnani who suggested the problem and the possible solutions.








\bibliographystyle{model1-num-names}
