\section{Introduction and Background}

The concept of light detection and
ranging (LIDAR), which can be thought of as the optical analogue of radar,
is by no means new, and during the many years it has been in use,
it has found an extremely large variety of applications in a wide spectrum of
areas and disciplines, including: agriculture,
archaeology, biology, astronomy, geology and soil sciences, forestry,
meteorology, and military applications, just to name a few.
Most notably, there are several modern technologies that involve LIDAR, such as
autonomous vehicles, space flight devices, robots of many kinds,
systems with GRID based processing, and in the future, and face recognition in
biometric systems (e.g., at airports).

This background sets the stage and motivates renewed interest in optical signal detection and estimation.
The customary model of a direct--detection optical receiver (or detector) consists of a
photo-diode (PIN diode or avalanche photo-diode), that converts the intensity of the received
optical (laser) signal, modeled as the
rate function of a variable--rate Poisson process, into a train of current impulses
generated by the photo--electrons at random time instants, pertaining to the Poisson arrivals.
This current is then fed into some electronic circuitry, whose
first stage is normally a trans-impedance amplifier
(TIA), that amplifies the current signal and converts it into a relatively strong voltage signal, but this
amplification comes at the cost of some distortion as well as thermal noise associated with the
electronic circuitry.

The challenge in the development of a solid
theory of detection and estimation, for such a signal--plus--noise system, is that it is has a
rather complicated model due to the various types of noise involved.
The combination of shot noise, due to the photo-diode, and the thermal electronic noise is already not trivial.
In the case of an avalanche photo--diode (APD), which is the more relevant case,
there is an additional, third type of noise, namely, the {\it excess noise}, induced by the APD,
which is actually a multiplicative noise
process pertaining to fluctuations in the random gain associated with the
avalanche mechanism, as each primary
electron--hole pair (generated by a photon absorbed in the photo-diode),
may generate secondary electron--holes,
which in turn can generate additional electron--holes, and so on.
For more details, and additional aspects of the problem area, the interested reader is
referred to some earlier work, e.g.,
\cite{BL90}, \cite{Einarsson96}, \cite{EH81}, \cite{FGS75}, \cite{GJN09},
\cite{GP87}, \cite{MS76}, \cite{Olsson89},
\cite{Personick71a}, \cite{Personick71b}, \cite{Personick08}, \cite{Salz85},
\cite{SSH18}, \cite{TM94}, and \cite{YLC19},
which is by no means an exhaustive list of relevant articles (and a book).

To provide just a rough, preliminary view on the problem and to fix ideas,
we now give an informal presentation of the model and explain the
difficulties more concretely. The received signal is modeled as
\begin{equation}
\label{sigmod}
y(t)=\sum_{k=1}^K g_kh(t-t_k)+n(t),~~~~~~0\le t < T,
\end{equation}
where $K$ is the number of photo--electrons generated during the time interval $[0,T)$,
$h(\cdot)$ is the current pulse contributed by a single electron (which is nearly equal to the
charge of the electron multiplied by the Dirac delta
function), $\{t_k\}$ are the random Poisson arrival times induced by
the optical signal, $\lambda(t)$, sensed by the APD, $\{g_k\}$ are the random gains induced by
the APD photo-multiplication, and $\{n(t)\}$ is thermal noise, modeled here,
and in earlier works, to be
white Gaussian noise with spectral density $N_0/2$.\footnote{The flat spectrum
assumption is adopted here mainly for the sake of simplicity of the exposition. The extension to colored
noise is not difficult.}
The signal detection problem, in its basic form, is about
binary hypothesis testing. The null hypothesis is that $y(t)=n(t)$, whereas
the alternative is as in (\ref{sigmod}). Had $K$, $\{g_k\}$ and $\{t_k\}$ been known to the receiver,
the likelihood ratio (LR) would have been readily given by
(see, e.g., \cite{FGS75}):
\begin{eqnarray}
\label{lr}
L&=&\frac{\exp\left\{-\frac{1}{N_0}\int_0^T\left[y(t)-\sum_{k=1}^Kg_kh(t-t_k)\right]^2\mbox{d}t\right\}}
{\exp\left\{-\frac{1}{N_0}\int_0^Ty^2(t)\mbox{d}t\right\}}\nonumber\\
&=&\exp\left\{\frac{2}{N_0}\sum_{k=1}^K\int_0^Ty(t)h(t-t_k)\mbox{d}t-
\frac{1}{N_0}\sum_{k=1}^K\sum_{l=1}^Kg_kg_l
R(t_k-t_l)\right\},
\end{eqnarray}
where $R(\tau)=\int_0^Th(t)h(t-\tau)\mbox{d}t$.
Since these random parameters are unknown, the actual LR must be
obtained by taking the expectation of $L$ with
respect to (w.r.t.) their randomness. Deriving this expectation appears to be notoriously difficult,
mainly due to the second term at the exponent, i.e., the double sum over $k$ and $l$.

It is this difficulty that
triggered many researchers in the field to harness
their wisdom in the quest for satisfactory solutions, and
accordingly, there is rich literature on the subject,
dating back many years into the past. As far as general guidelines go,
a possible approach to alleviate this difficulty is
the {\it estimator--correlator approach} \cite{Kailath69}, which asserts that
the expected LR of detection of a random signal in Gaussian additive white noise
is given by the same expression as
if the desired signal, $\sum_kg_kh(t-t_k)$, was known (i.e., the same as if
$K$, $\{g_k\}$ and
$\{t_k\}$ were known), except that it is replaced by its causal,
minimum mean--square error (MMSE) estimator given $\{y(t)\}$.
The caveat, however, is clear:
deriving this MMSE causal estimator is an extremely difficult problem on its own.

To the best of the author's knowledge, the first article that is directly relevant to this kind of
study, for the above described specific signal model,
is the article by Foschini, Gilbert and Salz \cite{FGS75}. Their approach was to view the factor
associated with the double sum in the second line of (\ref{lr}), namely, the term,
$\exp\left\{-\frac{1}{N_0}\sum_{k,l}g_kg_lR(t_k-t_l)\right\}$,
as the characteristic function of the random variable $\sum_k g_kx(t_k)$ (for given
$K$, $\{g_k\}$ and $\{t_k\}$), where $\{x(t)\}$
is an auxiliary zero--mean, stationary Gaussian process with auto-correlation function $R(\tau)$.
At the next step, the expectation over the randomness of $\{x(t)\}$ was commuted with the
expectations over $K$, $\{t_k\}$ and $\{g_k\}$, which are easier to carry out for a given realization
of $\{x(t)\}$. The result is a more compact expression of the LR, but even after this simplification,
it is not explicit enough to be implementable in practice, or to analyze its
performance in full generality. At this point, the approach taken in \cite{FGS75} was
to carry out a series of approximations, yielding explicit
asymptotic forms of the optimal detector and its performance at least in the
limits of very low and very high signal--to--noise ratio (SNR). The resulting
approximate LRT for high SNR, however, was still rather complicated to implement.
Also, the behavior for moderate SNR was left open.

A year later, Mazo and Salz \cite{MS76} studied the performance of integrate--and--dump filters and
also obtained exact formulas for the random gain of the APD on the basis of the earlier study by
Personick \cite{Personick71a}, \cite{Personick71b}. See also \cite{TM94}.
Kadota \cite{Kadota88} has also derived an approximate LR test for a model like (\ref{sigmod}).
His approximation approach was different from that of \cite{FGS75}.
It was based on neglecting the effect of overlaps between
localized noise elements, which basically amounts to ignoring the cross terms of the double summation in
the exponent of (\ref{lr}) on the ground that $R(\cdot)$ is a very narrow function (see also \cite{Hero91}
who used the same approximation for the purpose of estimation).
More recently, Helstorm and Ho \cite{HH92} and
Ho \cite{Ho95} have applied saddle--point integration and thereby
studied the behavior of certain pulse shapes
at the optical receiver in terms of the performance of the decoder. Other studies are guided by the
approach of approximating the distribution of the shot noise of the
photo-diode by the Gaussian distribution, owing to considerations
in the spirit of the central limit theorem
(CLT), see e.g., \cite[Subsections 5.6.3, 5.8.4]{Einarsson96}, \cite{EH81}. In this context,
the well--known {\it optical matched filter}
(see, e.g., \cite{GP87}, \cite{Hero91}) is the main building block
of the optimal detector that simply maximizes the SNR at the sampling time, $t=T$.
This Gaussian approximation approach, however, raises some concerns since the CLT is not valid
for assessing the tails of the distribution and in particular, error exponents, which
are the relevant players when probabilities of large deviations events, like (the asymptotically
rare) FA and MD error
events, are studied.

In this paper, we take a different approach.
Motivated by considerations of the desired simplicity of
optical detectors for LIDAR
systems (especially when they need to be implemented on mobile devices), we consider
the class of optical signal detectors that are based on correlating the noisy received signal with
a given deterministic waveform, and we
characterize the waveform with the best trade--off between the
false--alarm (FA) probability and the missed--detection (MD) probability.
More precisely, our derivation addresses the trade-off between
the asymptotic error exponents of the FA and MD probabilities using
Chernoff bounds, without resorting to Gaussian approximations. We also provide numerical
results that compare the performance of the best correlator to that of
the optical matched filter (or, more precisely, the matched correlator),
which is coherent with the above--mentioned Gaussian approximation approach.
It is demonstrated that the proposed optimal
correlator outperforms the optical matched correlator, in terms of the
trade--off between the FA and the MD error exponents.
It should be pointed out that in addition to the random fluctuations of the
APD photo-multiplier, our model also incorporates the effect of
dark current that exists even under the null hypothesis.

In the same spirit and with a similar motivation,
we also study, for the same type of signal model, the
problem of estimating the delay of a received signal
on the basis of maximizing the correlation between the received signal and
a time--shifted waveform, as a function
of this shift. We characterize the optimal correlator waveform
that minimizes the mean square error (MSE)
in the regime of high SNR, as an extension of the analysis
provided by Bar-David \cite{BarDavid69}, who analyzed the high--SNR MSE of the
maximum likelihood (ML) estimator
for the pure Poissonian regime (i.e., without thermal noise).
Once again, the emphasis is on simplicity and
therefore, the performance of this estimator cannot be compared
to the much more complicated, approximate MAP estimator due to Hero \cite{Hero91}, which is based on
approximating the likelihood function, using the same approach as Kadota \cite{Kadota88}.

The optimal correlator waveforms for detection and for estimation turn out to be different,
but their limiting behavior is the same in
both detection and estimation problems:
when the thermal Gaussian noise is dominant,
the optimal correlator waveform becomes proportional to
the clean signal (like the classical matched filter for
additive white Gaussian noise), but when the thermal noise is
negligible compared to the other noises, then it becomes logarithmic
function of the clean signal, as expected in view of \cite{BarDavid69}.

The outline of the remaining part of the paper is as follows.
In Section \ref{thesignalmodel}, we present the model under discussion in full detail.
In Section \ref{detection}, we address the signal
detection problem, first and foremost, for the case of zero dark current. The case of positive
dark current, which follows the same general ideas
(but more complicated), is also outlined, but relatively briefly.
Finally, in Section \ref{estimation}, we address the problem of time delay estimation.

\section{The Signal Model}
\label{thesignalmodel}

In this section, we provide a formal presentation of the signal model,
that was briefly described in the Introduction. As mentioned before,
we consider the model,
\begin{equation}
\label{model}
y(t)=\sum_{k=1}^K g_k h(t-t_k)+n(t),~~~~~~~~~~~~0\le t < T,
\end{equation}
whose various ingredients are described as follows. The variable
$K$ is a Poissonian random variable, distributed according to
\begin{equation}
\mbox{Pr}\{K=\kappa\}=e^{-\Lambda}\frac{\Lambda^\kappa}{\kappa!},
~~~~~~\kappa=0,1,2\ldots,~~~~\Lambda=\int_0^T\lambda(t)\mbox{d}t,
\end{equation}
where $\lambda(t)$ is the a rate function
that depends upon the intensity of the received optical signal. In particular,
\begin{equation}
\lambda(t)=\frac{\eta{\cal P}(t)}{\hbar\omega}+\lambda_{\mbox{\tiny d}},
\end{equation}
where $\eta$ is the quantum efficiency of the APD, ${\cal P}(t)$ is the
instantaneous power of the optical signal,
$\hbar$ is Planck's constant, $\omega$ is the angular frequency of the light wave, and
$\lambda_{\mbox{\tiny d}}$ is the dark current. The variables
$\{g_k\}$ are independently identically distributed (i.i.d.) positive integer
random variables that designate the avalanche gains. According to Personick \cite{Personick71a},
\cite{Personick71b}, the distribution of these random variables depends on
the physics of the APD, and its characteristic function obeys a certain implicit equation,
which is solvable in closed form when only the
electrons (and not holes) cause ionizing collisions. In this case,
the distribution of each $g_k$ is geometric:
\begin{equation}
\label{geom}
\mbox{Pr}\{g_k=g\}=(e^\zeta-1)\cdot e^{-\zeta g},~~~~~~g=1,2,\ldots,~~\zeta > 0.
\end{equation}
For the sake of concreteness, we will henceforth adopt the assumption of this geometric distribution.
It also includes the case of a deterministic gain ($g_k =1$ with probability
one), which corresponds to the case of the PIN diode, by taking the limit
$\zeta\to\infty$.
The function $h(t)$ is the current pulse contributed
by the passage of a single photo--electron and hence its integral must be equal to the
electric charge of the electron, $q_{\mbox{\tiny e}}$. Naturally,
this is a very narrow pulse, which for most practical purposes, can be approximated by $h(t)\approx
q_{\mbox{\tiny e}}\delta(t)$, where $\delta(t)$ is the Dirac delta function. However, $h(t)$ can also be
understood to include the convolution with some front--end
filter, which is part of the electronic circuitry (e.g., the TIA).
The times $\{t_k\}$ are the random Poissonian
photon arrival times, taking on values in $[0,T)$ and being induced by the optical waveform, $\lambda(t)$.
Finally, $\{n(t)\}$ is Gaussian white noise
with spectral density $N_0/2$, which is assumed to be independent of $K$, $\{g_k\}_{k=1}^K$ and
$\{t_k\}_{k=1}^K$. Also, given $K$, $\{g_k\}$ are statistically independent of $\{t_k\}$. Recall that for
the underlying Poissonian process defined, conditioned on the event $K=\kappa$, the
unordered random arrival times,
$t_1,\ldots,t_\kappa$, are i.i.d.\ and their common density function
is given by $f(t)=\lambda(t)/\Lambda$, for $0\le t\le T$, and $f(t)=0$ elsewhere.

Observe that if we present $g_k$ as $\bar{g}+\Delta g_k$, where
$\bar{g}=\bE\{g_k\}$ and then $\Delta g_k$ designates the fluctuation,
then the random input signal can be represented as
\begin{eqnarray}
\sum_k(\bar{g}+\Delta g_k)h(t-t_k)&=&\bar{g}\sum_kh(t-t_k)+\sum_k\Delta
g_kh(t-t_k)\nonumber\\
&=&\bar{g}[q_{\mbox{\tiny e}}\lambda(t)+n_{\mbox{\tiny s}}(t)]+\sum_k\Delta
g_kh(t-t_k)\nonumber\\
&=&\bar{g}q_{\mbox{\tiny e}}\lambda(t)+\bar{g}n_{\mbox{\tiny
s}}(t)+\sum_k\Delta g_kh(t-t_k).
\end{eqnarray}
The first term, $\bar{g}q_{\mbox{\tiny e}}\lambda(t)$, is the desired clean
signal (plus dark current), the second term, defined as
\begin{equation}
\bar{g}n_{\mbox{\tiny s}}(t)=\bar{g}\left[\sum_kh(t-t_k)-q_{\mbox{\tiny e}}\lambda(t)\right],
\end{equation}
is shot noise (amplified by $\bar{g}$), and the last
term is multiplicative noise. Thus, together with the Gaussian noise, $n(t)$,
of (\ref{model}), there are three types of noise in this model, as already mentioned in the Introduction.

\section{Signal Detection}
\label{detection}

In this section, we study the signal detection problem.
We begin with the case $\lambda_{\mbox{\tiny d}}=0$ (no dark current),
which is considerably simpler, and then
outline the extension to the more general case,
$\lambda_{\mbox{\tiny d}} > 0$. Before, we move into the technical details, a comment is in order,
and it applies even to the case of no dark current.
Consider then the signal detection problem of deciding between the two
hypotheses:
\begin{eqnarray}
& &{\cal H}_0:~y(t)=n(t)\\
& &{\cal H}_1:~y(t)=\sum_{k=1}^K g_k h(t-t_k)+n(t),
\end{eqnarray}
using a detector, that is based on a correlator, that is, calculating the quantity
$$\int_0^T w(t)y(t)\mbox{d}t$$
and comparing it to a threshold, $\theta T$. Here,
$\{w(t),~0\le t\le T\}$ is a deterministic waveform to be optimized,
and $\theta > 0$ is a threshold parameter that controls the trade-off between the
FA and MD probabilities.
Clearly, if
the noise was purely Gaussian white noise,
under both hypotheses, the optimal choice of $w(t)$ would
have been matched to the desired signal,
i.e., $w(t)=\lambda(t)$. Here, however, as explained in Section
\ref{thesignalmodel}, there are two additional types of non--Gaussian noise under ${\cal H}_1$.
Since the classical matched correlator, $w(t)=\lambda(t)$, is no longer
necessarily optimal under non--Gaussian noise, and since we would
still be interested in a detector that is relatively
easy to implement, the natural question is whether
there is a waveform $w(t)$ better than $w(t)=\lambda(t)$. If so, then what is the optimal waveform $w_*(t)$
for this detection problem?
In this context, we should also mention again the notion of the optical matched filter (see
e.g., \cite{GP87}, \cite{Hero91}), the optical analogue of the classical matched filter,
which maximizes the SNR at the sampling time, $t=T$, taking
into account that the intensity of the shot noise is proportional to the
desired signal, $\lambda(t)$ (unlike the case of pure thermal noise).
But the relevance of the SNR as the only
parameter that counts for detection performance is valid
only under the Gaussian regime, so the optical matched filter
is applied either under the Gaussian approximation, or when only second moments are important.
The Gaussian approximation, under this model,
is largely justified by CLT considerations. Here, however, we wish
to avoid CLT considerations, as we are interested in the FA and MD error exponents, which are,
in fact, given by large deviations rate functions.
As is well known, the CLT is not valid for the large-deviations regime and for tails
of distributions.
Indeed, the optimal $w_*(t)$ that we derive below will be different from the optical matched filter.

\subsection{The Case of No Dark Current}

Consider the hypothesis testing problem defined at the beginning of the
introductory part of this section, which
is for the case of no dark current.
Let $E\dfn\int_0^Tw^2(t)\mbox{d}t$ and $P\dfn E/T$. Then,
the FA probability is given by
\begin{eqnarray}
P_{\mbox{\tiny FA}}&=&\mbox{Pr}\left\{\int_0^T w(t)n(t)\mbox{d}t > \theta
T\right\}\nonumber\\
&=&Q\left(\frac{\theta T}{\sqrt{N_0E/2}}\right)\nonumber\\
&\exe&\exp\left\{-\frac{\theta^2T^2}{N_0E}\right\}\nonumber\\
&=&\exp\left\{-\frac{\theta^2T}{N_0P}\right\}\nonumber\\
&\dfn&e^{-E_{\mbox{\tiny FA}}(\theta)T}
\end{eqnarray}
where the notation $\exe$ designates asymptotic equivalence
in the exponential scale, i.e., for two positive functions of $T$, $a(T)$ and $b(T)$,
the assertion $a(T)\exe b(T)$,
stands for the assertion that
$\lim_{T\to\infty}\frac{1}{T}\log\frac{a(T)}{b(T)}=0$.
Since the FA error exponent, $E_{\mbox{\tiny FA}}(\theta)=\theta^2/N_0P$,
depends on the waveform $w(\cdot)$ only via its
power, $P$, it is obvious that the maximization of the MD error exponent for given
FA error exponent is equivalent to its maximization subject to a power
constraint imposed on $w(\cdot)$.
Denoting $Z\dfn\int_0^T w(t)n(t)\mbox{d}t$, we next assess the MD error exponent using the Chernoff bound,
assuming that $h(t)=q_{\mbox{\tiny e}}\delta(t)$.
\begin{eqnarray}
P_{\mbox{\tiny MD}}&=&\mbox{Pr}\left\{\int_0^T w(t)y(t)\mbox{d}t \le \theta
T\right\}\nonumber\\
&=&\mbox{Pr}\left\{q_{\mbox{\tiny e}}\sum_{k=1}^Kg_kw(t_k)+Z \le \theta
T\right\}\nonumber\\
&\le&\inf_{s\ge 0} e^{s\theta
T}\bE\left\{\exp\left[-sq_{\mbox{\tiny e}}\sum_{k=1}^Kg_kw(t_k)-sZ\right]\right\}\nonumber\\
&=&\inf_{s\ge 0} \exp\left\{s\theta
T+\frac{s^2N_0PT}{4}\right\}\bE\left\{\exp
\left[-sq_{\mbox{\tiny e}}\sum_{k=1}^Kg_kw(t_k)\right]\right\}\nonumber\\
&=&\inf_{s\ge 0} \exp\left\{s\theta
T+\frac{s^2N_0PT}{4}\right\}\bE\left\{\prod_{k=1}^Ke^{-sq_{\mbox{\tiny e}}g_kw(t_k)}\right\}.
\end{eqnarray}
To calculate the last expectation, we proceed in two steps. First, we average
each factor over $g_k$, assuming that it is geometrically distributed as in (\ref{geom}).
This gives
\begin{eqnarray}
\sum_{g=1}^\infty P(g)e^{-gsq_{\mbox{\tiny e}}w(t_k)}&=&(e^\zeta-1)
\sum_{g=1}^\infty e^{-g[sq_{\mbox{\tiny e}}w(t_k)+\zeta]}\nonumber\\
&=&\frac{e^{\zeta}-1}{e^{sq_{\mbox{\tiny e}}w(t_k)+\zeta}-1}.
\end{eqnarray}
As a second step, we average over $K$ and $\{t_k\}_{k=1}^K$, and get
\begin{eqnarray}
\bE\left\{\prod_{k=1}^Ke^{-sq_{\mbox{\tiny e}}g_kw(t_k)}\right\}&=&\bE\left\{\prod_{k=1}^K\frac{e^{\zeta}-1}
{e^{sq_{\mbox{\tiny e}}w(t_k)+\zeta}-1}\right\}\nonumber\\
&=&\exp\left\{\int_0^T\lambda(t)\left[\frac{e^{\zeta}-1}
{e^{sq_{\mbox{\tiny e}}w(t)+\zeta}-1}-1\right]\mbox{d}t\right\}\nonumber\\
&=&\exp\left\{-e^{\zeta}\int_0^T\lambda(t)\cdot\frac{e^{sq_{\mbox{\tiny e}}w(t)}-1}
{e^{sq_{\mbox{\tiny e}}w(t)+\zeta}-1}\mbox{d}t\right\},
\end{eqnarray}
where we have used the fact
\cite[eqs.\ (8)--(13)]{BarDavid69} that for an arbitrary positive function $f$,
\begin{equation}
\bE\left\{\prod_{k=1}^K f(t_k)\right\}=\exp\left\{\int_0^T\lambda(t)[f(t)-1]\mbox{d}t\right\}.
\end{equation}
Thus,
\begin{equation}
P_{\mbox{\tiny MD}}\le
\exp\left\{-T\sup_{s\ge
0}\left[\frac{e^{\zeta}}{T}\int_0^T\lambda(t)\cdot\frac{e^{sq_{\mbox{\tiny e}}w(t)}-1}
{e^{sq_{\mbox{\tiny e}}w(t)+\zeta}-1}\mbox{d}t-s\theta-s^2\frac{N_0P}{4}\right]\right\}.
\end{equation}
Consider first the special case where $g_k\equiv 1$ with probability one,
which is obtained in the limit
$\zeta\to\infty$. In this case, the above simplifies to
\begin{equation}
P_{\mbox{\tiny MD}}\le
\exp\left\{-T\sup_{s\ge 0}\left[\frac{1}{T}\int_0^T\lambda(t)[1-e^{-sq_{\mbox{\tiny e}}w(t)}]
\mbox{d}t-s\theta-s^2\frac{N_0P}{4}\right]\right\}.
\end{equation}

\noindent
{\it Remark.} Note that had the channel been purely Gaussian, that is, without the shot noise, and the
desired signal was $q_{\mbox{\tiny e}}\lambda(t)$, we would have obtained
\begin{equation}
P_{\mbox{\tiny MD}}\le
\exp\left\{-T\sup_{s\ge 0}\left[\frac{1}{T}\int_0^Tq_{\mbox{\tiny e}}\lambda(t)\cdot sw(t)
\mbox{d}t-s\theta-s^2\frac{N_0P}{4}\right]\right\}.
\end{equation}
This means that the difference
$$\frac{1}{T}\int_0^T\lambda(t)\cdot[sq_{\mbox{\tiny e}}w(t)-(1-e^{-sq_{\mbox{\tiny e}}w(t)})]\mbox{d}t$$
designates the loss due to the additional shot noise.\footnote{The integrand
is, of course, non--negative due to the inequality $x\ge 1-e^{-x}$.}

We therefore need to maximize the exponent over both $s$ and $\{w(t),~0\le t\le T\}$.
For a given $s$, the optimal $\{w(t)\}$ minimizes
$\int_0^T\lambda(t)e^{-sq_{\mbox{\tiny e}}
w(t)}\mbox{d}t$ subject to the power constraint
$\int_0^Tw^2(t)\mbox{d}t\le PT$, which is equivalent to minimizing
\begin{equation}
\int_0^T\lambda(t)e^{-sq_{\mbox{\tiny e}}w(t)}\mbox{d}t+
\frac{s^2q_{\mbox{\tiny e}}^2}{2c}\left[\int_0^Tw^2(t)\mbox{d}t-PT\right],
\end{equation}
where $s^2q_{\mbox{\tiny e}}^2/2c > 0$ is a Lagrange
multiplier.\footnote{This form of the Lagrange multiplier is adopted as it allows a convenient
representation of the solution. It is legitimate since there is complete freedom to control
its value by the choice of the constant $c > 0$.}
Finding the optimal function $w(\cdot)$ is a standard problem in
calculus of variations, whose solution is characterized as follows.
Let the function $p[\cdot]$ denote the inverse of the monotonically increasing function
$b[x]\dfn xe^x$, $x\ge 0$, i.e., $p[y]$ is the solution $x$ to the equation $xe^x=y$,
$y\ge 0$. The optimal $w(t)$ is given by
\begin{equation}
w_*(t)=\frac{1}{sq_{\mbox{\tiny e}}}p[c\cdot\lambda(t)],
\end{equation}
where $c > 0$ is chosen such that
$\int_0^Tp^2[c\cdot\lambda(t)]\mbox{d}t=s^2q_{\mbox{\tiny e}}^2PT$.
The MD exponent is then given by
\begin{equation}
E_{\mbox{\tiny MD}}(\theta)=\sup_{s\ge
0}\left[\frac{1}{T}\int_0^T\lambda(t)(1-\exp\{-p[c\cdot\lambda(t)]\})
\mbox{d}t-s\theta-s^2\frac{N_0P}{4}\right],
\end{equation}
where it should be kept in mind that $c$ depends on $s$.

Note that if the optimal $s$ is very small (which is the case when $N_0$ and/or $\theta$
are large, then the solution to the equation
$sq_{\mbox{\tiny e}}w(t)e^{sq_{\mbox{\tiny e}}w(t)}=c\cdot\lambda(t)$
is found near the origin, where $p[x]\approx x$, which means that $w_*(t)$ is
nearly proportional to $\lambda(t)$, namely, the classical matched correlator. If, on
the other hand, the optimal $s$ is very large (which is the case when $N_0$ and $\theta$ are both
small), then the solution is found away
from the origin, where $p[x]\approx \ln x$. In this case,
$w_*(t)\approx\frac{1}{sq_{\mbox{\tiny e}}}\ln[c\cdot\lambda(t)]$, in agreement of optimal
photo--counting detector (see, e.g., \cite{BarDavid69}),
which is obtained in the absence of Gaussian noise.\\

\noindent
{\it Example 1.}
Consider the frequently--encountered case where $\lambda(t)$ is a two--level signal,
where half of the time $\lambda(t)=\lambda_1$ and in the other half,
$\lambda(t)=\lambda_2$. Then, $w_*(t)$ must also be a two--level signal. Owing to the power constraint,
we may denote these two levels by $w$ and $\sqrt{2P-w^2}$, respectively, and it remains to maximize the
exponent over $s$ and $w$ alone. Specifically, we have
\begin{equation}
E_{\mbox{\tiny MD}}(\theta)=\sup_{s\ge
0}\max_{0\le w\le\sqrt{2P}}\left[\frac{\lambda_1}{2}(1-e^{-sq_{\mbox{\tiny e}}w})+
\frac{\lambda_2}{2}(1-e^{-sq_{\mbox{\tiny e}}\sqrt{2P-w^2}})
-s\theta-s^2\frac{N_0P}{4}\right].
\end{equation}
Fig.\ \ref{graph1} compares the exponent of the optimal waveform, $w_*$, to that of
the optical matched filter \cite[eq.\ (8)]{GP87},
\begin{equation}
w_{\mbox{\tiny omf}}(t)=\frac{\lambda(t)}{\lambda(t)+N_0/(2q_{\mbox{\tiny e}}^2\overline{g^2})},
\end{equation}
for the following values of the parameters of the problem:
$P=10$, $N_0/q_{\mbox{\tiny e}}^2=0.0001$, $\lambda_1=1$ and $\lambda_2=10$.
The numerical value of the spectral density
of the thermal noise was deliberately chosen extremely small, in order
to demonstrate a situation where the noise is far from being Gaussian, and
thereby examine sharply the validity of the Gaussian approximation. In
this case, $w_{\mbox{\tiny omf}}(t)$ is nearly equal to unity for all $t$,
which means a pure, unweighted integrator \cite[p.\ 1291, Remark 2]{GP87}.
As can be seen, the optimal correlator $w_*$ improves upon the optical matched
correlator fairly significantly, especially for large values of the threshold
parameter, $\theta$. This concludes Example 1. $\Box$

\begin{figure}[h!t!b!]
\centering
\includegraphics[width=8.5cm, height=8.5cm]{md-exponent-dgain-opticalmf.eps}
\caption{MD error exponents of the optical matched correlator (dashed blue
curve) and the
optimal correlator $w_*$ (solid red curve) as functions of
$\theta$, for a deterministic gain
($\zeta\to\infty$), $P=10$,
$N_0/q_{\mbox{\tiny e}}^2=0.0001$, $\lambda_1=1$ and $\lambda_2=10$.}
\label{graph1}
\end{figure}

Returning to the case of a general, finite $\zeta$,
and carrying out a similar optimization, we find that
the optimal $w_*(t)$ is now given by
\begin{equation}
w_*(t)=\frac{1}{sq_{\mbox{\tiny e}}}\cdot p_\zeta[c\cdot\lambda(t)],
\end{equation}
where $p_\zeta$ is the inverse of the function
\begin{equation}
b_\zeta[x]=\frac{x(e^{x+\zeta}-1)^2}{e^{x+2\zeta}-e^{x+\zeta}},
\end{equation}
and where, once again, $c$ is chosen such that
\begin{equation}
\int_0^Tp_\zeta^2[c\cdot\lambda(t)]\mbox{d}t=s^2q_{\mbox{\tiny e}}^2PT.
\end{equation}
Here too, if $N_0$ and/or $\theta$ are large, then $s$ must be small, and then
due to the power constraint, $c$ must be small too, which means that the
functions $p_\zeta$ and $b_\zeta$ operate near the origin,
where they are roughly linear, as
\begin{equation}
b_\zeta[x]\approx x\cdot\frac{(e^\zeta-1)^2}{e^{2\zeta}-e^\zeta}.
\end{equation}
At the other extreme, on the other hand, $p_\zeta$ and $b_\zeta$ operate away from the origin,
where
\begin{equation}
b_\zeta[x]\approx
\frac{xe^{2x+2\zeta}}{e^{x+2\zeta}-e^{x+\zeta}}=\frac{xe^{x+\zeta}}{e^\zeta-1},
\end{equation}
which is again, nearly exponential, and so, $p_\zeta$ is approximately logarithmic,
as before. The MD exponent is therefore given by
\begin{equation}
E_{\mbox{\tiny MD}}(\theta)=
\sup_{s\ge
0}\left[\frac{e^{\zeta}}{T}\int_0^T\lambda(t)\cdot\frac{\exp\{p_\zeta[c\cdot\lambda(t)]\}-1}
{\exp\{p_\zeta[c\cdot\lambda(t)]+\zeta\}-1}\mbox{d}t-s\theta-s^2\frac{N_0P}{4}\right].
\end{equation}

\noindent
{\it Example 2.} Consider again the setting of Example 1 above, except that here $\zeta$ is finite.
Specifically, Fig.\ 2 displays a comparison analogous to that of Fig.\ 1
for the case of a random gain with
parameter $\zeta=0.1$, where for the red solid curve, $s$ and $w$ are chosen to maximize
\begin{equation}
E_{\mbox{\tiny MD}}(\theta)=\sup_{s\ge
0}\max_{0\le
w\le\sqrt{2P}}\left[\frac{\lambda_1e^\zeta}{2}\cdot\frac{e^{sq_{\mbox{\tiny e}}w}-1}
{e^{sq_{\mbox{\tiny e}}w+\zeta}-1}+
\frac{\lambda_2e^\zeta}{2}\cdot\frac{e^{sq_{\mbox{\tiny e}}\sqrt{2P-w^2}}-1}
{e^{sq_{\mbox{\tiny e}}\sqrt{2P-w^2}+\zeta}-1}
-s\theta-s^2\frac{N_0P}{4}\right],
\end{equation}
and for the blue, dashed
curve, $w$ is chosen to be the optical matched correlator as before.
As can be seen, here too, the optimal $w$ improves upon the optical matched
correlator and the gap is rather considerable, especially as $\theta$ grows.

\begin{figure}[h!t!b!]
\centering
\includegraphics[width=8.5cm, height=8.5cm]{md-exponent-rgain-opticalmf.eps}
\caption{MD error exponents of the optical matched correlator (blue dashed
curve) and
the optimal correlator (red solid curve) as functions of $\theta$, for a random gain
($\zeta=0.1$), $P=10$,
$N_0/q_{\mbox{\tiny e}}^2=0.0001$, $\lambda_1=1$ and $\lambda_2=10$.}
\label{graph2}
\end{figure}

\subsection{The Case of Positive Dark Current}

The case of $\lambda_{\mbox{\tiny d}} > 0$ is analyzed on the basis of similar ideas, and we
therefore cover it relatively briefly, highlighting mostly the points where there is a substantial
difference relative to the zero dark--current case.

Extending the analysis to the case of positive dark current, a similar
derivation yields the following FA exponent for a given correlator waveform, $w(t)$:
\begin{equation}
E_{\mbox{\tiny FA}}(\theta,w)=\sup_{0\le s<
\zeta/w_{\max}}\left[s\theta-\frac{s^2N_0P}{4}-\frac{\lambda_{\mbox{\tiny d}}
e^\zeta}{T}\int_0^T\frac{e^{sq_{\mbox{\tiny e}}w(t)}-1}{e^\zeta-e^{sq_{\mbox{\tiny e}}w(t)}}\mbox{d}t\right],
\end{equation}
where $w_{\max}=\sup_{0\le t\le T}w(t)$.
Here, the trade--off between the FA and the MD error exponents is somewhat more
involved than in the zero dark--current case,
since the FA error exponent depends on $w(\cdot)$ in a more complicated manner
than just via its power $P$.
In particular, we would now like
to find the pair $(\theta,w)$ that maximizes $E_{\mbox{\tiny MD}}(\theta,w)$ over all pairs
such that $E_{\mbox{\tiny FA}}(\theta,w)\ge S$, for a given $S >
0$ that designates the target FA exponent.\footnote{Previously,
$S$ was given by $\theta^2/N_0P$, so for a given $P$, $S$ was
proportional to $\theta^2$.} Equivalently, we would like to solve the problem
\begin{equation}
\sup_{\theta,w}\left\{E_{\mbox{\tiny
MD}}(\theta,w)+\mu E_{\mbox{\tiny FA}}(\theta,w)\right\},
\end{equation}
where $\mu$ is a Lagrange multiplier chosen to meet the constraint, $E_{\mbox{\tiny FA}}(\theta,w)\ge S$.
More specifically, using Chernoff bounds for both error exponents, this
amounts to solving the problem,
\begin{eqnarray}
& &\sup_{s\ge 0}\sup_{\sigma\ge 0}\sup_{\theta\ge 0}\sup_{w:~w_{\max} < \zeta/s}\left[
\frac{e^{\zeta}}{T}\int_0^T\lambda(t)\cdot\frac{e^{\sigma  q_{\mbox{\tiny e}} w(t)}-1}
{e^{\sigma q_{\mbox{\tiny e}}
w(t)+\zeta}-1}\mbox{d}t-\sigma\theta-\sigma^2\frac{N_0P}{4}+\right.\nonumber\\
& &\left.\mu s\theta-\frac{\mu s^2N_0P}{4}-\frac{\mu\lambda_0
e^\zeta}{T}\int_0^T\frac{e^{sq_{\mbox{\tiny e}}w(t)}-1}{e^\zeta-
e^{sq_{\mbox{\tiny e}}w(t)}}\mbox{d}t\right]\nonumber\\
&=&\sup_{s\ge 0}\sup_{\sigma\ge 0}\sup_{\theta\ge 0}\sup_{w:~w_{\max} <
\zeta/s}\left[
\frac{e^{\zeta}}{T}\int_0^T\left(\lambda(t)\cdot\frac{e^{\sigma q_{\mbox{\tiny e}}w(t)}-1}
{e^{\sigma q_{\mbox{\tiny e}}
w(t)+\zeta}-1}-\mu\lambda_0
\frac{e^{sq_{\mbox{\tiny e}}w(t)}-1}{e^\zeta-e^{sq_{\mbox{\tiny e}}w(t)}}\right)\mbox{d}t+\right.\nonumber\\
& &\left.(\mu s-\sigma)\theta-(\sigma^2+\mu s^2)\frac{N_0P}{4}\right].
\end{eqnarray}
We will not continue any further to the full, detailed solution of this
problem, beyond the
following comment which applies to the case of a deterministic gain.
In the limit of $\zeta\to\infty$, the above trade-off yields
\begin{equation}
w_*(t)=\ln\left[\frac{\sigma}{\lambda_{\mbox{\tiny d}}\mu s}\cdot\lambda(t)\right],
\end{equation}
in other words, in the presence of dark-current, the relation is always
logarithmic.

\section{Time Delay Estimation}
\label{estimation}

In this section, we consider the problem of time
delay estimation. The underlying model is the same as
before, except that the optical signal
is time--shifted, i.e., $\lambda(t-\theta)$, where $\theta$ is the
delay. It is assumed that the support of
$\lambda(t-\theta)$ in included in the interval $[0,T]$, for every
$\theta$ in the range of uncertainty.
As mentioned in the Introduction, this is a relevant problem in LIDAR systems, where distances to
certain objects have to be estimated, similarly as in classical radar systems.
Here too, our basic building block is a correlator. Consider an estimator of the
form,
\begin{equation}
\hat{\theta}=\arg\max_\theta Q(\theta),
\end{equation}
where
\begin{equation}
Q(\theta)=\int_0^Ty(t)w(t-\theta)\mbox{d}t,
\end{equation}
and where $w(\cdot)$ is a twice--differentiable waveform to be optimized, with the property
that the temporal cross--correlation function,
\begin{equation}
R_{w\lambda}(\tau)=\int_0^T\lambda(t)w(t-\tau)\mbox{d}t
\end{equation}
achieves its maximum at $\tau=0$. This implies
\begin{equation}
\dot{R}_{w\lambda}(0)=\int_0^T\lambda(t)\dot{w}(t)\mbox{d}t=0,
\end{equation}
where dotted functions designate derivatives. Similarly as the assumption concerning $\lambda(\cdot)$,
we also assume that the support of the waveform
$w(t-\theta)$ is fully included in $[0,T]$ for the entire range of search of the estimated delay.

Our analysis begins similarly as in \cite{BarDavid69},
which assumes high SNR and small estimation errors. Accordingly, consider the
Taylor series expansion of first derivative, $\dot{Q}(\theta)$, around the
true parameter value, $\theta_0$:
\begin{equation}
0=\dot{Q}(\hat{\theta})\approx\dot{Q}(\theta_0)+(\hat{\theta}-\theta_0)\ddot{Q}(\theta_0),
\end{equation}
where $\ddot{Q}(\theta)$ is the second derivative of $Q(\theta)$. This yields
\begin{eqnarray}
\label{approxerror}
\hat{\theta}-\theta_0&\approx&-\frac{\dot{Q}(\theta_0)}{\ddot{Q}(\theta_0)}\nonumber\\
&=&-\frac{\bar{g}q_{\mbox{\tiny
e}}\int_0^T\lambda(t-\theta_0)\dot{w}(t-\theta_0)\mbox{d}t+\int_0^Tn_T(t)\dot{w}(t-\theta_0)\mbox{d}t}
{\bar{g}q_{\mbox{\tiny e}}\int_0^T\lambda(t-\theta_0)\ddot{w}(t-\theta_0)\mbox{d}t
+\int_0^Tn_T(t)\ddot{w}(t-\theta_0)\mbox{d}t}\nonumber\\
&=&-\frac{\bar{g}q_{\mbox{\tiny
e}}\int_0^T\lambda(t)\dot{w}(t)\mbox{d}t+\int_0^Tn_T(t)\dot{w}(t-\theta_0)\mbox{d}t}
{\bar{g}q_{\mbox{\tiny e}}\int_0^T\lambda(t)\ddot{w}(t)\mbox{d}t
+\int_0^Tn_T(t)\ddot{w}(t-\theta_0)\mbox{d}t}\nonumber\\
&=&-\frac{
\int_0^Tn_T(t)\dot{w}(t-\theta_0)\mbox{d}t}
{\bar{g}q_{\mbox{\tiny e}}\cdot\int_0^T\lambda(t)\ddot{w}(t)\mbox{d}t
+\int_0^Tn_T(t)\ddot{w}(t-\theta_0)\mbox{d}t}\nonumber\\
&\dfn&\frac{U_n}{A+U_d},
\end{eqnarray}
where $n_T(t)$ is the total noise,
composed of the shot noise, the multiplicative noise and the thermal noise,
i.e.,
\begin{equation}
n_{\mbox{\tiny T}}(t)=n_{\mbox{\tiny s}}(t)+n_{\mbox{\tiny m}}(t)+n(t),
\end{equation}
where
\begin{equation}
n_{\mbox{\tiny s}}(t)=\bar{g}\cdot\left[\sum_kh(t-t_k)-q_{\mbox{\tiny
e}}\lambda(t-\theta)\right],
\end{equation}
and
\begin{equation}
n_{\mbox{\tiny m}}(t)=\sum_k\Delta g_kh(t-t_k).
\end{equation}
Thus,
\begin{equation}
\bE\{(\hat{\theta}-\theta)^2\}\approx\frac{\bE\{U_n^2\}}{A^2},
\end{equation}
where we have made a further approximation by neglecting the contribution of
the random variable $U_d$ relative to the deterministic constant $A$
(at the denominator of the last line of (\ref{approxerror})) since $A$ is proportional to $T$, whereas
the standard deviation of $U_d$ is proportional
to $\sqrt{T}$ (see also \cite[eqs.\ (38)--(40)]{BarDavid69} for a more
detailed justification of a similar approximation).
It is easy to see that all three noise components are uncorrelated with each other.
The auto-correlation function of the thermal noise is
$R_n(t,s)=\frac{N_0}{2}\delta(t-s)$. The auto-correlation function of the (amplified) shot noise is
\begin{eqnarray}
R_{\mbox{\tiny
s}}(t,s)&=&\bar{g}^2\cdot\bE\left\{\left[\sum_kh(t-t_k)-q_{\mbox{\tiny
e}}\lambda(t)\right]\cdot\left[\sum_kh(s-t_k)-q_{\mbox{\tiny
e}}\lambda(s)\right]\right\}\nonumber\\
&=&\bar{g}^2\cdot\bE\left\{\sum_{k,l}h(t-t_k)h(s-t_l)\right\}-\bar{g}^2q_{\mbox{\tiny
e}}^2\lambda(t)\lambda(s)\nonumber\\
&=&\bar{g}^2\cdot\bE\left\{\sum_{k}h(t-t_k)h(s-t_k)\right\}+
\bar{g}^2\cdot\bE\left\{\sum_{k\ne l}h(t-t_k)h(s-t_l)\right\}
-\bar{g}^2q_{\mbox{\tiny
e}}^2\lambda(t)\lambda(s)\nonumber\\
&=&\bar{g}^2\cdot\bE\{K\}\frac{1}{\Lambda}
\int_0^T\lambda(\tau-\theta)h(t-\tau)h(s-\tau)\mbox{d}\tau
+\nonumber\\
& &\bar{g}^2\cdot\bE\{K(K-1)\}\frac{1}{\Lambda^2}
\left[\int_0^T\lambda(\tau-\theta)h(t-\tau)\mbox{d}\tau\right]
\left[\int_0^T\lambda(\tau-\theta)h(s-\tau)\mbox{d}\tau\right]-\bar{g}^2q_{\mbox{\tiny
e}}^2\lambda(t)\lambda(s)\nonumber\\
&=&\bar{g}^2\cdot\int_0^T\lambda(\tau-\theta)h(t-\tau)h(s-\tau)\mbox{d}\tau.
\end{eqnarray}
Similarly,
\begin{eqnarray}
R_{\mbox{\tiny
m}}(t,s)&=&\bE\{n_{\mbox{\tiny m}}(t)n_{\mbox{\tiny m}}(s)\}\nonumber\\
&=&\bE\left\{\sum_{k,l}\Delta g_k\Delta g_lh(t-t_k)h(s-t_l)\right\}\nonumber\\
&=&\mbox{Var}\{g\}\cdot\bE\left\{\sum_{k}h(t-t_k)h(s-t_k)\right\}+
\bE\left\{\sum_{k\ne l}\Delta g_k\Delta g_lh(t-t_k)h(s-t_l)\right\}\nonumber\\
&=&\mbox{Var}\{g\}\cdot\bE\left\{\sum_{k}h(t-t_k)h(s-t_k)\right\}\nonumber\\
&=&\mbox{Var}\{g\}\cdot\int_0^T\lambda(\tau-\theta)h(t-\tau)h(s-\tau)\mbox{d}\tau.
\end{eqnarray}
It follows that
\begin{equation}
R_{\mbox{\tiny T}}(t,s)\dfn R_{\mbox{\tiny
s}}(t,s)+R_{\mbox{\tiny
m}}(t,s)+R_n(t,s)=\overline{g^2}\cdot\int_0^T\lambda(\tau-
\theta)h(t-\tau)h(s-\tau)\mbox{d}\tau+\frac{N_0}{2}\delta(t-s).
\end{equation}
Thus,
\begin{eqnarray}
\bE\{U_n^2\}&=&\int_0^T\int_0^T\mbox{d}t\mbox{d}sR_{\mbox{\tiny
T}}(t,s)\dot{w}(t-\theta)\dot{w}(s-\theta)\nonumber\\
&=&\frac{N_0}{2}\int_0^T\dot{w}^2(t-\theta)\mbox{d}t+\overline{g^2}
\int_0^T\int_0^T\mbox{d}t\mbox{d}s\dot{w}(t-\theta)\dot{w}(s-\theta)\int_0^T
\lambda(\tau-\theta)h(t-\tau)h(s-\tau)\mbox{d}\tau\nonumber\\
&=&\frac{N_0}{2}\int_0^T\dot{w}^2(t-\theta)\mbox{d}t+\overline{g^2}
\int_0^T\mbox{d}\tau\lambda(\tau-\theta)\left[\int_0^T\mbox{d}t\dot{w}(t-\theta)h(t-\tau)\right]
\left[\int_0^T\mbox{d}s\dot{w}(s-\theta)h(s-\tau)\right]\nonumber\\
&\approx&\frac{N_0}{2}\int_0^T\dot{w}^2(t-\theta)\mbox{d}t+\overline{g^2}q_{\mbox{\tiny
e}}^2\int_0^T\lambda(t)\dot{w}^2(t)\mbox{d}t\nonumber\\
&=&\int_0^T\left[\frac{N_0}{2}+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)\right]\dot{w}^2(t)\mbox{d}t,
\end{eqnarray}
which yields
\begin{equation}
\bE\{(\hat{\theta}-\theta)^2\}\approx
\frac{\int_0^T\left[\frac{N_0}{2}+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)\right]\dot{w}^2(t)\mbox{d}t}{\overline{g}^2q_{\mbox{\tiny
e}}^2\left[\int_0^T\lambda(t)\ddot{w}(t)\mbox{d}t\right]^2}.
\end{equation}
It is desired to find the optimal waveform $\{w(t)\}$.
First, observe that the MSE is invariant to scaling of $\{w(t)\}$, so the
problem is equivalent to maximizing the absolute value of
\begin{equation}
\int_0^T\lambda(t)\ddot{w}(t)\mbox{d}t
\end{equation}
subject to a given value of
\begin{equation}
\label{cons}
\int_0^T\left[\frac{N_0}{2}+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)\right]\dot{w}^2(t)\mbox{d}t.
\end{equation}
Let $\omega(t)=\dot{w}(t)$ and
\begin{equation}
v(t)=\omega(t)\cdot\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}.
\end{equation}
Then, the problem is equivalent to that of
maximizing the absolute value of
\begin{equation}
\int_0^T\lambda(t)\cdot\left[\frac{v(t)}{\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}}\right]'\mbox{d}t,
\end{equation}
or equivalently,
\begin{equation}
\lambda(t)\cdot\frac{v(t)}{\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}}\Bigg|_0^T-
\int_0^T\dot{\lambda}(t)\cdot\frac{v(t)}{\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}}\mbox{d}t,
\end{equation}
subject to a given value of (\ref{cons}).
Now, assuming that
$\lambda(0)=\lambda(T)=0$, the contribution of the
first term in the last expression vanishes, and so,
we wish to maximize the absolute value of
\begin{equation}
\int_0^T\dot{\lambda}(t)\cdot\frac{v(t)}{\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}}\mbox{d}t
\end{equation}
for a given energy of $\{v(t)\}$. The maximum is achieved for
\begin{equation}
v(t)=\frac{\dot{\lambda}(t)}{\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}},
\end{equation}
which yields
\begin{equation}
\omega(t)=\frac{v(t)}{\sqrt{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)}}=\frac{\dot{\lambda}(t)}{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(t)},
\end{equation}
and so, finally,
\begin{equation}
w_*(t)=\int_0^t\frac{\dot{\lambda}(s)\mbox{d}s}{N_0/2+\overline{g^2}q_{\mbox{\tiny
e}}^2\lambda(s)}=\frac{1}{\overline{g^2}q_{\mbox{\tiny e}}^2}
\ln\left[\frac{N_0}{2}+\overline{g^2}q_{\mbox{\tiny e}}^2\lambda(t)\right]+c,
\end{equation}
where $c$ is an integration constant (not to be confused with the earlier
defined constant, $c$).
Equivalently, after a simple manipulation
of the integration constant, $w_*(t)$ can be presented as
\begin{equation}
w_*(t)\propto\ln\left[1+\frac{2\overline{g^2}q_{\mbox{\tiny
e}}^2}{N_0}\cdot\lambda(t)\right].
\end{equation}
As can be seen, for very large $N_0$, $w(t)$ is approximately proportional to $\lambda(t)$, as expected.
For very small $N_0$, $w(t)$ is proportional to
$\ln\lambda(t)$, which is also expected. Note that if $\lambda(t)$ includes a
dark--current component, $\lambda_{\mbox{\tiny d}}$, then it simply adds to the spectral
density of the thermal noise, $N_0/2$. In other words, there is no distinction
here between the dark--current shot noise and the Gaussian thermal noise.
This should not be surprising, because when calculating the MSE, only second
moments are relevant, and not finer details of the distributions of the
various kinds of noise.

Let us compare the MSE of the optimal correlator $w_*$ to that of the ordinary
matched correlator, $w(t)=\lambda(t)$. Denote
$\lambda_0=2\overline{g^2}q_{\mbox{\tiny e}}^2/N_0$. Then, the MSE of the latter
is given by
\begin{equation}
\mbox{MSE}
=\frac{1}{\lambda_0}\cdot\frac{\int_0^T[1+\lambda(t)/\lambda_0]\dot{\lambda}^2(t)\mbox{d}t}
{\left[\int_0^T\dot{\lambda}^2(t)\mbox{d}t\right]^2},
\end{equation}
whereas the former is given by
\begin{equation}
\mbox{MSE}_{\mbox{\tiny opt}}
=\frac{1}{\lambda_0}\cdot\frac{1}{\int_0^T\frac{\dot{\lambda}^2(t)\mbox{d}t}{1+\lambda(t)/\lambda_0}}.
\end{equation}
It is easy to see that $\mbox{MSE}_{\mbox{\tiny opt}}$ is indeed smaller than
$\mbox{MSE}$, due to the Schwartz--Cauchy inequality, as
\begin{eqnarray}
\left[\int_0^T\dot{\lambda}^2(t)\mbox{d}t\right]^2&=&
\left[\int_0^T\frac{\dot{\lambda}(t)}{\sqrt{1+\lambda(t)/\lambda_0}}\cdot
\dot{\lambda}(t)\sqrt{1+\frac{\lambda(t)}{\lambda_0}}\mbox{d}t\right]^2\nonumber\\
&\le&\left[\int_0^T\frac{\dot{\lambda}^2(t)}{1+\lambda(t)/\lambda_0}\mbox{d}t\right]\cdot
\left[\int_0^T\dot{\lambda}^2(t)\left(1+\frac{\lambda(t)}{\lambda_0}\right)\mbox{d}t\right].
\end{eqnarray}
The difference becomes smaller when
$\dot{\lambda}(t)\sqrt{1+\lambda(t)/\lambda_0}$ becomes
closer to be proportional to
$\dot{\lambda}(t)/\sqrt{1+\lambda(t)/\lambda_0}$, namely, when
$\lambda(t)\ll \lambda_0$, as expected.

This high--SNR analysis can be extended to handle non--white noise and
even non--stationary thermal noise of auto-correlation function
$R_n(s,t)=\bE\{n(s)n(t)\}$. The only difference is that now the numerator of
(\ref{approxerror}) will be replaced by
\begin{equation}
\int_0^T\dot{w}(s)\dot{w}(t)[R_n(s,t)+\overline{g^2}q_{\mbox{\tiny
e}}^2\sqrt{\lambda(s)\lambda(t)}\delta(t-s)]\mbox{d}s\mbox{d}t
\dfn\int_0^T\dot{w}(s)\dot{w}(t)R_0(s,t)\mbox{d}s\mbox{d}t.
\end{equation}
The resulting optimum correlator would then be given by
\begin{equation}
w_*(t)=\int_0^t\mbox{d}\tau\int_0^T\mbox{d}s\dot{\lambda}(s)R_0^{-1}(s,\tau),
\end{equation}
where $R_0^{-1}(s,\tau)$ is the formal inverse
of the kernel $R_0(s,t)$ (provided that it exists), namely, the
kernel that satisfies
\begin{equation}
\int_0^TR_0(s,t_1)R_0^{-1}(s,t_2)\mbox{d}s=\delta(t_1-t_2).
\end{equation}

The other side of the coin of high--SNR estimation performance is the
probability of anomaly (see, e.g., \cite[Chap. 8]{WJ65}).
The common practice in assessing this probability is
to divide the interval $[0,T]$ into small bins whose sizes are about the pulse
width and to apply a union bound on the
(essentially identical) probabilities of the pairwise events that the maximum
$Q(\theta)$ in each bin exceeds the one that contains the correct delay.
Since the number of bins grows only linearly in $T$, the exponent of this
probability is the same as that of each individual bin. The analysis of each
pairwise event can be carried out very similarly to the earlier analysis of
the FA probability (just with a few minor twists), whose dependence on the mismatched correlator $w$ is
different that of the above derived approximated MSE.
Thus, the signal design should seek a compromise that trades off
the weak--noise MSE with the probability of anomaly. In the high--SNR regime, a
logarithmic correlator is good both for the MSE and for keeping the anomaly
probability small.

\section*{Acknowledgment}

Interesting discussions with Moshe Nazarathy, at the early stages of this work, are
acknowledged with thanks.


\section{Introduction}
Rice (Oryza sativa L.)  is a primary staple food in the world \cite{tubiello2013faostat}. China and India are the leading rice-producing nation in the world. These two countries with other Asian nations contribute 90\% of world rice production \cite{muthayya2014overview}. In current climate changing scenarios, food security is the prime concern in these highly populated developing nations, especially India and China \cite{tubiello2013faostat}, \cite{muthayya2014overview}.  To cater to the need of the rising population in these countries, constant efforts are being made to achieve high yield through genetic improvements programs. To mitigate the issues arising due to change in the climate in recent decades (IPPC, 2007, \cite{cruz2007asia}) demands climate-resilient genotype for these crops. So there is a need to identify the donor rice genotypes with higher precision in breeding programs. Conventionally, in traditional phenotyping, the often-prevalent recurrent selection breeding techniques are applied for the initial screening and evaluation of genotypes to develop a new variety  \cite{dhondt2013cell}. Moreover, traditional phenotyping involves mostly human senses, which limits high efficiency, and are prone to human error. These approaches are destructive, invasive, time consuming, less effective and restricts to achieve high precision. Furthermore, for high precision observations and combating these bottlenecks of traditional phenotyping, the High Throughput Plant Phenotyping evolved and has many advantages \cite{cobb2013next}, \cite{seelig2008assessment}, \cite{araus2014field}. The HTPP applies various sensors, imaging techniques, and Remote Sensing (RS) principles for recording plant observations and weather parameters (environmental variability) for phenotyping in an automatic and semi-automatic manner \cite{araus2014field}, \cite{bai2016multi}, \cite{ghanem2015physiological}. Many morphological traits like plant height, plant canopy spreading type, convex hull, plant leaf area, number of leaves in a plant, color, etc. are used for phenotyping \cite{bardenas1965morphology}, \cite{upov1987guidelines}. These traits are also used for grouping of the genotypes and understanding the diversity among genotypes as well as species.

In this study, we have considered many image-based traits like the number of leaves per plant, the ratio of the number of leaves per plant to the per unit area of its convex hull, area of the convex hull for each plant, the plant spread (canopy spread) based on the top view of the plant images, to classify the diverse rice (Oryza sativa L.) genotypes based on the top view of the plant. For leaf counting, we used deep learning techniques, You Only Look Once (YOLO) \cite{redmon2013darknet}. YOLO is convolutional neural networks based state-of-the-art object detection algorithm \cite{redmon2017yolo9000}, \cite{redmon2018yolov3}, \cite{ren2015faster}. We employed YOLO for leaf detection owing to its speed for real-time object detection. Apart from this, the convex hull of all genotype, vertical spread, and horizontal spread of the canopy is calculated for diversity analysis.
\begin{figure}[ht]
\centering
  \subfigure[Rice genotypes in a controlled greenhouse]
    \includegraphics[width=6cm]{a.jpg}
  }
  \quad
  \subfigure[Rice plants on the car, moving for imaging]
    \includegraphics[width=6cm] {b.jpg}
  }
  \caption{Experimental set up  at the Nanaji Deshmukh Plant Phenomics Centre (NDPPC)}
  \label{fig:setup}
\end{figure}
\section{Dataset and Methods}
\subsection{Dataset}
The experimental setup is shown in Figure \ref{fig:setup}. In this study, datasets are generated for diverse rice genotypes grown in pots in the climate-controlled greenhouses at Nanaji Deshmukh Plant Phenomics Centre (NDPPC), Indian Council of Agricultural Research-Indian Agricultural Research Institute (ICAR-IARI) located in New Delhi, India. The data is a subset of 150 rice genotypes grown for Genome-Wide-Association-Study (GWAS) to identify the major Quantitative Trait Loci (QTL) responsible for drought tolerance. Pot culture experiment was conducted with one set of plants (3 replications for each genotype) with moisture-deficit stress (~12.5 \% soil moisture content) during maximum tillering stage (35 to 55 days after transplanting), while the other set grown in well watered (25\%) condition. The RGB visual images of pot-grown rice plants of each selected genotypes were captured using a visual camera (Prosilica GT6600 series for monochrome and color image,RGB of spectral range 400-700 nm)installed in scanalyzer3D, LemnaTec, imaging platform (LemnaTec GmbH, Aachen, Germany) at NDPPC. We capture the top view and side view images of plants having a size $6576 \times 4384$ recorded over an interval of five days. The rice genotypes were selected for high polymorphic variation in the leaf number.
\begin{figure}[ht]
\centering
  \subfigure[RGB image of a rice plant, top view]
    \includegraphics[width=4cm]{cfinal3.png}
  }
  \quad
  \subfigure[Annotations of leaf tip for a rice plant]
    \includegraphics[width=4cm]{d.jpg}
  }
  \caption{Data annotation procedure}
  \label{fig:annotation}
\end{figure}
\subsection{Methods}
The original images were pre-processed and resized for the study.  After pre-processing, around 300 images were annotated, as shown in Figure \ref{fig:annotation}. These images are used for the training deep learning model, YOLO \cite{redmon2013darknet}, \cite{redmon2017yolo9000}, \cite{redmon2018yolov3}. YOLO has been implemented for leaves tip detection and estimation of the number of leaves. In the case of a rice plant, a leaf has a single tip. Therefore, a hypothesis was made that the number of leaves is equal to the number of tips for a rice plant. Few leaf counting works were done using convolution networks and de-convolution networks \cite{dobrescu2017leveraging}, \cite{paul2017count}, \cite{segui2015learning}, \cite{aich2017leaf}. Rice plant has a diverse plant architecture. Variation from sparse to dense leaves (canopy spread), and overlapping leaves create challenges for detection and segmentation of leaves in a rice plant. However, tips are rarely overlapped and are easily detectable, except leaves and their tips from few unproductive tillers that are crowded and hidden near culms and leaf sheath area in a rice plant just above-ground. The convex hull of each plant is calculated from predicted bounding boxes and their coordinates of leaves tips. This methodology has the advantage of drawing a convex hull simultaneously with predicted bounding boxes and their coordinates for train as well as test data. Few samples of leaves tips predicted by YOLO and convex hull formed by corresponding leaves tips are shown in Figure \ref{fig:results_1_samples}. The ratio of the number of leaves in a plant to the area of convex hull formed by leaf tips is considered to analyze plant architecture and understand canopy spreading architecture of rice plants under drought stress. When a plant goes under stress, then leaf rolling occurs and thus reduces the convex hull and compactness of the plant for spread type plants. Based on the number of leaves for each genotype and the number of leaves per unit area of the convex hull, we clustered the genotypes. This will help the breeder to infer about the phenotypic behavior of the genotype and assist in selecting the donors for further breeding programme. The clustering of genotype based on the number of leaves per unit area of the convex hull also has the potential for drought and agronomic study as we observe the vigour and canopy spread of the plant.
\begin{figure}[t]
\centering
  \subfigure[]
    \includegraphics[width=6.5cm]{results_1.png}
  }
  \quad
  \subfigure[]
    \includegraphics[width=6.5cm]{Graph_ICLR_After_Review.jpg}
  }
  \caption{(a) Variations in the number of leaves, and the number of leaves per unit area of the convex hull of a plant for a genotype in control (top) and drought (bottom) condition. (b) Figure shows variation in the number of leaves for Anjali genotype in control and drought conditions during stress period (35 to 55 DAT)}
  \label{fig:results_1_samples}
\end{figure}
\section{Results and Discussion}
\begin{figure}[ht]
\centering
  \subfigure[Considering all parameters]
    \includegraphics[width=6cm]{dendo1.png}
    \label{dendo1}
  }
  \quad
  \subfigure[Considering the number of leaves per plant]
    \includegraphics[width=6cm]{dendo2.png}
    \label{dendo2}
  }
  \caption{Dendogram constructed with Ward clustering method for ten (10) genotypes of rice to discriminate drought and non-drought genotypes (0-ANJALI, 1-BLACKGORA, 2-ABHAYA X DAGADESI, 3-HEERA, 4-NAGINA 22, 5-RASI, 6-DULAR, 7-PMK-2, 8-SERATOES, 9-KALINGA-1)}
  \label{fig:dendo}
\end{figure}
A subset of 10 genotypes out of 150 genotypes was taken for current work, which has a total of approximately 400 images of these genotypes. After pre-processing the acquired images, we trained the YOLO for leaves tips detection. Leaf tip annotation, as shown in Figure 2. (b) for a rice plant as tip (object) detection, annotated tips are inside bounding boxes. These same predicted tips of each plants are used for calculation of convex hull for each plant for each replication of  genotypes i.e. plant in the experiment. We propose the number of leaves per unit area of the convex hull as a 2D trait in plant phenotyping. Figure \ref{fig:results_1_samples} depicts the variation of the number of leaves and convex hull for a genotype, Anjali, in controlled and drought conditions. We observed that the number of leaves setting in stress condition plants i.e. drought susceptible genotypes is less as compared to controlled plants as well as the drought tolerant plants. Also, for spread type plants, the convex hull was reduced in stress conditions. We observe that the number of leaves per unit area of the convex hull will be a potential trait to understand plant architecture and understanding different stages of plants, as result confirms the difference between drought and control plants based on architecture and number of leaves. Apart from this, we extracted features such as the  number of leaves detected and estimated by YOLO, area of convex hull formed by tips of leaves, horizontal spread of the plant, vertical spread of the plant and the ratio of the  the numbers of leaves to the area of convex hull to discriminate these genotypes. Figure \ref{fig:results_1_samples} shows detected leaf tips by YOLO and convex hull for control and drought conditions. Our model detects upto  90 percent of tips and leaves accurately.

Figure \ref{fig:dendo} shows dendogram obtained using all the features i.e. the number of leaves in plant, area of convex hull of the plant, the number of leaves per plant per unit area convex hull of the plants, vertical and horizontal spread (spacing) of the plant Figure \ref{dendo1} and using only the number of leaves in a plants Figure, \ref{dendo2} based on controlled and drought conditions. Dendogram in Figure \ref{dendo1} is based on all parameters like the number of leaves,vertical and horizontal plant spread (considering top view of RGB image), the convex hull i.e. area under convex hull for each plants, and the number of leaves per unit area of convex hull of each plant. Results shows that genotypes  6 (DULAR), 7  (PMK-2) and 8 (SERATOES) fall in one group are very similar genotypes, same observation is confirmed in dendogram constructed based on the number of leaves, where 6, 7 and 8 are very in one group or similar genotypes. Genotype 9 (KALINGA-1) is distinct based on all parameters wheres considering the number of leaves genotype 4 (NAGINA 22) is distinct to others. These genotypes can be a basis for further genetic improvement programs.

\section{Conclusions and Future Work}
Image-based phenotyping has an advantage over traditional phenotyping. Object detection techniques have potential in leaves counting tasks as we trained and demonstrated YOLO considering tip as an object in the current study. The number of leaves per convex hull of a plant is an inherent trait for image-based phenotyping. This can be the basis for monitoring the rate of leaf development in field crops like rice, wheat, maize in the early vegetative stages, and simulate it for different growth stages for different genotypes in a controlled environment. This can be used for further abiotic stress studies and phenotyping.
\subsubsection*{Acknowledgment}
Authors acknowledge scientists of Nanaji Deshmukh Plant Phenomics Center (NDPPC), IARI, New Delhi, India, and the ICAR-IARI for providing the high throughput image data to undertake this research work. The NDPPC is created and supported by National Agriculture Science Fund (NASF), ICAR, Grant No.NASF/Phen-6005/2016–17.


\section{Introduction}
Speech emotion recognition is currently gaining interest from both academia
and commercial industries. Researchers in the affective computing field progressively
proposed new methods to improve the accuracy of automatic emotion recognition.
Commercial industries are trying to make this technology
available to the market since its potential applications. Previously,
researchers has attempted to implement speech-based emotion recognition for
wellbeing detection \cite{Gao2018}, call center application \cite{Petrushin1998}, and automotive safety
\cite{Nass2005}.

One of the common requirements in computing speech emotion recognition is the
availability of high-performance computing since the dataset usually is very large
in size, and the classifying methods are complicated.
Graphical processing units (GPU)-based computers are often used
over CPU-based computers due to its processing speed to process the data, particularly,
image-like data.

This paper proposes the use of deep multilayer perceptrons (MLP) to overcome
the requirement of high computing power required by modern deep learning
architectures. The inputs are high-level statistical functions (HSF), which are
used to reduce the dimension of input features. The outputs are emotion
dimensions, i.e.,
degree of valence, arousal, and dominance.

According to research in psychology, dimensional emotion is another view
in emotion theories apart from categorical emotion.
Russel \cite{Russell1979} argued that emotion categories could be derived from this
dimensional emotion, particularly in valence-arousal space.
Given the benefit of the ability to convert dimensional
emotion to categorical emotion, but not vice versa,
predicting the emotion dimension is more beneficial than
predicting the emotion category. We added dominance, since it is suggested in
\cite{Bakker2014}, and the availability of those labels in the datasets.
Dimensional emotion recognition are evaluated with deep MLP
from speech data since the target applications are speech-based technology like call center
and voice assistant applications.

The contribution of this paper is the evaluation of the classical MLP technique
with deep layers compared with modern deep learning techniques, i.e., LSTM and CNN,
in terms of concordance correlation coefficient (CCC).
The deep neural network is an extension of the neural network with deeper layers,
commonly five or more layers \cite{Atmaja2019d}.
Our results show that on both speaker-dependent and speaker-independent,
in three datasets,
deep MLP obtained higher performances than LSTM and CNN.
The proposed method worked effectively on the small size
of feature, in which, this may be a limitation of our proposed deep MLP method.

\section{Data and Feature Sets}
In this section, we describe data and feature sets used in this
research.
\paragraph{IEMOCAP} The interactive emotional dyadic motion capture database
is used in this research \cite{Busso2008}. Although the database consists of the measurement of
speech, facial expression, head, and movements of affective dyadic sessions,
only speech data are processed. The database contains approximately 12 h of data
with 10039 utterances. All of those data are used.
The dimensional emotion labels are given in
valence, arousal, and dominance, in range [1-5] score and normalized
into [-1, 1] as in \cite{parthasarathy2017jointly} when those labels are fed into the neural network system.
For speech data, two versions are available in the dataset, stereo data per dialog
and mono data per sentence (utterance). We used the mono version since it is easy to
process with the labels. The sampling rate of the data was 16 kHz and 16-bit PCM.

We arranged the IEMOCAP dataset into two scenarios, speaker-dependent
(SD) and speaker-independent. On speaker-independent, we split the dataset with
ratio 80/20 for training/test set, while in speaker-independent, the last
session, i.e., session five, is left for the test set (leave one session out, LOSO).
The ratio of dataset partition in the speaker-independent scenario is similar
to speaker-dependent split which is shown in Table \ref{tab:data-partition}.

\paragraph{MSP-IMPROV} We used the MSP-IMPROV corpus to generalize the
impact of the evaluated methods. MSP-IMPROV is an acted corpus of dyadic
interactions to study emotion perception. This dataset consists of speech and
visual recording of 18 hours of affective dyadic sessions. Same as IEMOCAP dataset,
we only used the speech data, with 8438 utterances. The same split ratio
is used in speaker-dependent scenario while the last
session six is used for test set in in speaker-independent scenario,
with the same labels scale and normalization.
While IEMOCAP labels are annotated by
at least two subjects, these MSP-IMPROV labels are annotated by at least
five annotators. The speech data were mono, 44 kHz,
and 16-bit PCM.

Table \ref{tab:data-partition} shows the number of utterances allocated
for each set partition for both speaker-dependent and speaker-independent,
including MSP-IMPROV dataset.

\paragraph{Mixed-corpus} In addition to the two datasets above, we mixed
those two datasets to create a new category of dataset namely mixed-corpus.
In mixed-corpus, we concatenated speaker-dependent from IEMOCAP with
speaker-dependent from MSP-IMPROV for each, training, development
and test sets. The same rules also applied for the speaker-independent scenario.

\begin{table}[htbp]
\caption{Number of utterances used in the datasets on each partition}
\begin{center}
\begin{tabular}{l c c c}
\hline
Scenarios   &    Training   & Development   & Test  \\
\hline
IEMOCAP-SD  &   6431   &    1608    &   2000\\
IEMOCAP-LOSO    &   6295    &   1574    &   2170\\
IMPROV-SD   &   5256    &   1314    &   1868\\
IMPROV-LOSO     &   5452    &   1364    &   1622\\
\hline
\end{tabular}
\label{tab:data-partition}
\end{center}
\end{table}

\paragraph{Acoustic Feature Set}
We used high statistical functions of the low-level descriptor (LLD) from Geneva Minimalistic Acoustic Parameter Set (GeMAPS), which is developed by Eyben et al. \cite{eyben}. The HSF features are extracted per utterance depend on the given labels,
while the LLDs are processed on a frame-based level with 25 ms window size and 10 ms of hop size. The use of HSF feature reduce computation complexity since the feature size decrease from (3409 $\times$ 23) to  (1 $\times$ 23 features), that is, for
IEMOCAP dataset. To obtain the HSF feature, however, LLDs must obtained first. Then, HSF can be calculated as statistics of those LLDs for
a fixed time, in this case, per utterance.
\begin{table}
\centering
\caption{GeMAPS feature \cite{eyben} and its functionals; only functionals
        (HSFs) used in this dimensional SER.}
\begin{tabular}{l p{6.5cm}}
\hline
LLDs & loudness, alpha ratio, hammarberg index, spectral slope 0-500 Hz, spectral slope
500-1500 Hz, spectral flux, 4 MFCCs, F0, jitter, shimmer, Harmonics-to-Noise Ratio (HNR),
Harmonic difference H1-H2, Harmonic difference H1-A3, F1, F1 bandwidth, F1
amplitude, F2, F2 amplitude, F3, and F3 amplitude. \\
\hline
HSFs & mean (of LLDs), standard deviation (of LLDs), silence \\
\hline
\end{tabular}
\label{tab:feature}
\end{table}

To add those functionals, we used a silence feature, which is also extracted per utterance.
Silence feature, in this paper, is defined as the portion of the silence frames compared
to the total frames in an utterance. In human-human communication, this portion of
silence in speaking depends on the speaker's emotion. For example,
high arousal emotion category like happy may have fewer silences (or pauses)
than a sad emotion category. The ratio of silence per utterance is calculated as

\begin{equation} \label{eq:silence}
    S = \frac{N_{s}}{N_{t}},
\end{equation}
where $N_s$ is the number of frames to be categorized as silence (silence frames),
and $N_t$ is the number of total frames within an utterance. To be categorized as
silence, a frame is checked whether it is less than a threshold, which is a
multiplication of a factor with a root mean
square (RMS) energy ($X_{rms}$). This RMS energy is formulated as

\begin{equation}
    th = 0.3 \times \overline{X_{rms}}
\end{equation}
and $X_{rms}$ is defined as

\begin{equation}
    X_{rms} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}x[i]^2}
\end{equation}
where silence factor of $0.3$ is obtained from experiments.
These equations are similar to what is proposed in \cite{Sahu}
and \cite{moore2014word}. In \cite{Sahu}, the
author of that paper used a fixed threshold, while we evaluated some factors of
to find the best factor for silence feature in speech emotion
recognition. In \cite{moore2014word}, the author divided the total duration of
disfluency over the total utterance length on $n$ words and
counted it as a disfluency feature.

\section{Benchmarked and Proposed Method}
We evaluated three different methods in this paper, LSTM, CNN and MLP.
LSTM and CNN are used as baselines, while MLP with deep layers is the proposed
method. All evaluated methods used the same numbers of layers, units
and value of parameters.

\subsection{Benchmarked Methods: LSTM and CNN}
LSTM and CNN are two common deep learning architectures widely used in
speech emotion recognition \cite{Xie2019, Schmitt2019, Atmaja2019}.
We used those two architectures as the baselines due
to its reported effectiveness on predicting valence, arousal, and dominance.
Both LSTM and CNN evaluated here have the same five layers with
the same number of units. For the size of the kernel in CNN architecture, we
determined it in order that the number of trainable parameters is similar.
The other parameters, like batch size, feature and label standardization,
loss function, number of iterations, and callback criteria, are same
for both architectures.

Fig. \ref{fig:lstm-cnn} shows both structures of LSTM and CNN. On the
first layer, 256 neurons are used and decreased half for the next layers since
the models are supposed to learn better along with those layers. Five LSTM layers
used tanh as activation function, while five CNN layers used ReLU activation
function. We kept all output from the last LTM layer and flatten it before splitting
into three dense layers for obtaining the prediction of valence, arousal,
and dominance. For the CNN architecture, the same flatten layer
is used before entering three one-unit dense layers.

Both architectures used the same standardization for input features and labels.
The z-score normalization is used to standardize feature, while minmax scaler
is used to scale the labels into [0, 1] range. We used CCC
\cite{lawrence1989concordance} loss as the cost function
with multitask (MTL) approach in which the prediction of valence, arousal,
and dominance are done simultaneously. While CCC loss (CCCL) is used as the cost
function, the following CCC is used to evaluate the performance of recognition.
\begin{align}
CCC &= \dfrac{2\rho \sigma_x \sigma_y}{\sigma_x^2 + \sigma_y^2 + (\mu_x - \mu_y)^2} \\
CCCL &= 1 - CCC
\end{align}
where $\rho$ is Pearson's correlation between gold standard $y$ and and
predicted score $x$, $\sigma$ is standard deviation, and $\mu$ is the mean score.
The total loss ($L_T$) for three variables is then defined as the sum of CCCL for
those three with corresponding weighting factors,
\begin{align}
L_T  = \alpha CCCL_V + \beta CCCL_A + (1-\alpha-\beta) CCCL_D
\end{align}
where $\alpha$ and $\beta$ are weighting factors for loss function of valence
($CCCL_V$) and arousal ($CCCL_A$), respectively. The weighting factors of
loss function for dominance ($CCCL_D$) is obtained by subtracting 1 by the sum
of those weighting factors. The CCC is used to evaluate all architectures
including the proposed deep MLP method.

All architectures used a mini-batch size of 200 utterances by shuffling
its orders, maximum number iteration of 180, and 10 iterations patience
of early stopping criteria (callback).
An adam optimizer \cite{Kingma2014} is used to adjust the learning
rate during the training process. Both architectures run on GPU-based
machine using CuDNN implementation
\cite{Chetlur2014} within Keras toolkit \cite{chollet2015keras}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.45\textwidth]{./fig/eusipco_cnnlstm.pdf}}
\caption{Diagram of LSTM and CNN used for benchmarking of proposed deep MLP; the number inside the
bracket denoted the number of units (nodes); the second number on CNN denotes kernel size,
\# denotes number of trainable parameters.}
\label{fig:lstm-cnn}
\end{figure}

\subsection{Proposed Method: Deep MLP}
Fig. \ref{fig:deepmlp} shows our proposed deep MLP structure. The
MLP used here similar to the definition of connectionist learning proposed
by Hinton \cite{Hinton1989}. As the benchmarked methods,
deep MLP also has five layers with the same number of units as previous.
The only difference of the layer structure from the previous is the absent of
flatten layer since the output of the last MLP layers can be coupled directly to
three one-unit dense layers. While the previous LSTM and CNN used batch
normalization layer in the beginning (input) to speed-up computation process,
this deep MLP structure did not use that layer since the implementation
only need a minute to run on a CPU-based machine.

We used the same batch size, tolerance for early stopping criteria, optimizer,
and maximum number of iteration as the benchmarked methods. While the
benchmarked methods used CCC the as loss function, the proposed deep MLP method used
a mean square error (MSE) as the cost function,
\begin{equation}
MSE = \dfrac{1}{n} \sum_{i=1}^n (y_i-x_i)^2.
\end{equation}
The total loss function is given as an average of MSE scores from valence, arousal,
and dominance,
\begin{equation}
MSE_T = \dfrac{1}{3} (MSE_V + MSE_A + MSE_D).
\end{equation}
There are no weighting factors used here since we do not find a way to
implement it via scikit-learn toolkit \cite{scikit-learn}, in which
the proposed deep MLP is implemented. The same reason applied for the selection
of MSE over CCC for the loss function. The Python implementation codes for
both proposed and benchmarked methods are available at
\url{https://github.com/bagustris/deep_mlp_ser}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.4\textwidth]{./fig/deep_mlp.pdf}}
\caption{Diagram of proposed deep MLP with five layers; the number inside the
bracket denoted the number of units.}
\label{fig:deepmlp}
\end{figure}

\section{Experiment Results and Discussion}
CCC is the standard metric used in affective computing to measure the
performance of dimensional emotion recognition. We presented our results
in that metric in two different groups; within-corpus and mixed-corpus
evaluation. The results are shown in Table \ref{tab:result-within} and \ref{tab:result-mixed}.

Table \ref{tab:result-within} shows CCC scores of valence (V), arousal (A), dominance (D)
and its average from different datasets, scenarios, and methods. The proposed
deep MLP method outperforms benchmarked methods by remarkable margins.
On every emotion dimensions and averaged score, the proposed deep MLP gained
the highest CCC score for both speaker-dependent and speaker-independent
scenarios (typed in bold). On IEMOCAP dataset, the score of speaker-dependent
is only slightly
higher than speaker-independent due to the nature of dataset structure.
The utterances in IEMOCAP dataset is already in order by its session when it is
sorted by file names. The change from speaker-dependent to speaker-independent is done by
changing the number of train/test partitions. In contrast, the file naming
of utterances in MSP-IMPROV made the arrangement of the sessions not in order
when utterances are sorted by its file names. We did the sorting process
to assure the pair of features and labels. The gap between
speaker-dependent and speaker-independent in MSP-IMPROV is larger than in
IEMOCAP which may be caused by those different files organization. A case where our deep MLP method gained a lower score
is on dominance part of MSP-IMPROV speaker-dependent scenario,
however, the averaged CCC score, in that case, is still the highest.

\begin{table}[htbp]
\caption{Results of CCC scores on within-corpus evaluation}
\begin{center}
\begin{tabular}{c l c c c c}
\hline
Dataset &   Method   &    V   & A   &   D   & Mean\\
\hline
\parbox[t]{2mm}{\multirow{8}{*}{\rotatebox[origin=c]{90}{IEMOCAP}}}
& \multicolumn{5}{c}{\textit{speaker-dependent}} \\
& LSTM    &       0.222    &       0.508    &   0.432    & 0.387 \\
& CNN        &   0.086        &   0.433    &   0.401    & 0.307 \\
& MLP        &   \textbf{0.335}        &   \textbf{0.599}    &   \textbf{0.473}    & \textbf{0.469} \\
& \multicolumn{5}{c}{\textit{speaker-independent (LOSO)}} \\
& LSTM    & 0.210    & 0.474    & 0.394    & 0.359 \\
& CNN    & 0.113    & 0.460    & 0.410    & 0.328 \\
& MLP    & \textbf{0.316}    & \textbf{0.488}    & \textbf{0.454}    & \textbf{0.453} \\
\hline
\parbox[t]{2mm}{\multirow{8}{*}{\rotatebox[origin=c]{90}{MSP-IMPROV}}}
& \multicolumn{5}{c}{\textit{speaker-dependent}} \\
& LSTM    & 0.392    & 0.629    & \textbf{0.524}    & 0.515 \\
& CNN    & 0.346    & 0.623    & 0.522    & 0.497 \\
& MLP    & \textbf{0.438}    & \textbf{0.650}    & 0.519    & \textbf{0.536} \\
& \multicolumn{5}{c}{\textit{speaker-independent (LOSO)}} \\
& LSTM    & 0.210    & 0.483    & 0.313    & 0.335 \\
& CNN    & 0.216    & 0.524    & 0.387    & 0.375 \\
& MLP    & \textbf{0.269}    & \textbf{0.551}    & \textbf{0.401}    & \textbf{0.407} \\
\hline
\end{tabular}
\label{tab:result-within}
\end{center}
\end{table}

Table \ref{tab:result-mixed} shows the results from the mixed-corpus dataset. This
corpus is concatenation of IEMOCAP with MSP-IMPROV as listed in Table
\ref{tab:data-partition}, for both speaker-dependent and speaker-independent
scenarios. In this mixed-corpus, the proposed deep MLP method also outperformed
LSTM and CNN in all emotion dimensions and averaged CCC scores. The score
on speaker-dependent in that mixed-corpus is in between the score of
speaker-dependent in IEMOCAP and MSP-IMPROV within-corpus. For
speaker-independent, the score is lower than in within-corpus. This low score
suggested that speaker variability (in different sessions) affected the
result, even with the z-normalization process. Instead of predicting
one different session (LOSO), the test set in the mixed-corpus consists
of two different sessions, each from IEMOCAP and MSP-IMPROV, which made
regression task more difficult.

We showed that our proposed deep MLP functioned to overcome
the requirement of modern neural network architectures since it surpassed
the results obtained by those architectures. Using a small dimension of feature size,
i.e., 47-dimensional data, our deep MLP with five layers, excluding input
and output layers, achieved the highest performance. Modern deep learning
architectures require high computation hardware, e.g., GPU card, which
costs expensive. We showed that using a small deep MLP architecture, which
does not require high computation load, better performance can be achieved.
Our proposed deep MLP method gained a higher performance than benchmarked methods
not only on both within-corpus and mixed-corpus but also on both
speaker-dependent and speaker-independent scenarios. Although the proposed
method used the different loss function from the benchmarked methods,
i.e., MSE versus CCC loss, we presumed that our proposed deep MLP will
achieve higher performance if it used the CCC loss since the evaluation
metric is CCC.

\begin{table}[htbp]
\caption{Results of CCC scores on mixed-corpus evaluation}
\begin{center}
\begin{tabular}{l c c c c}
\hline
Method   &    V   & A   &   D   & Mean\\
\hline
\multicolumn{5}{c}{\textit{speaker-dependent}} \\
LSTM    &   0.262    &   0.518    &   0.424    &   0.401 \\
CNN    &   0.198    &   0.494    &   0.424    &   0.372 \\
MLP    &   \textbf{0.395}    &   \textbf{0.640}    &   \textbf{0.461}    &   \textbf{0.499} \\
\hline
\multicolumn{5}{c}{\textit{speaker-independent (LOSO)}} \\
LSTM    &   0.118    &   0.270    &   0.242    &   0.210 \\
CNN    &   0.073    &   0.265    &   0.249    &   0.196 \\
MLP    &   \textbf{0.212}    &   \textbf{0.402}    &   \textbf{0.269}    &   \textbf{0.294} \\
\hline
\end{tabular}
\label{tab:result-mixed}
\end{center}
\end{table}

\section{Conclusions}
This paper demonstrated that the use of deep MLP with proper parameter choices
outperformed the more modern neural network architectures with the same number
of layers. For both speaker-dependent and speaker-independent, the proposed
deep MLP gained the consistent highest performance among evaluated methods.
The proposed deep MLP also gained the highest score on both within-corpus and mixed-corpus
scenarios. Based on the results of these investigations, there is no high requirements
on computing power to obtain outstanding results on dimensional
speech emotion recognition. The proper choice of
feature (i.e., using small size feature) and the classifier can leverage the
performance of conventional neural networks.
Future research can be directed to
investigate the performance of the proposed method on cross-corpus evaluations,
which is not evaluated in this paper.

\bibliographystyle{IEEEtran}