# Doppler Estimation for High–Velocity Targets
Using Subpulse Processing and the Classic
Chinese Remainder Theorem

Fernando Darío Almeida García, André Saito Guerreiro, Gustavo Rodrigues de
Lima Tejerina,
José Cândido S. Santos Filho, Gustavo Fraidenraich, and Michel Daoud Yacoub F.
D. A. García, A. S. Guerreiro, G. R. de Lima Tejerina, J. C. S. Santos Filho,
G. Fraidenraich, and M. D. Yacoub are with the Wireless Technology Laboratory,
Department of Communications, School of Electrical and Computer Engineering,
University of Campinas, 13083-852 Campinas, SP, Brazil, Tel:+55(19)3788-5106,
e-mail:$\left\\{\right.$ferdaral, andsaito, tejerina, candido, gf,
<EMAIL_ADDRESS>
###### Abstract

In pulsed Doppler radars, the classic Chinese remainder theorem (CCRT) is a
common method to resolve Doppler ambiguities caused by fast-moving targets.
Another issue concerning high-velocity targets is related to the loss in the
signal-to-noise ratio (SNR) after performing range compression. In particular,
this loss can be partially mitigated by the use of subpulse processing (SP).
Modern radars combine these techniques in order to reliably unfold the target
velocity. However, the presence of background noise may compromise the Doppler
estimates. Hence, a rigorous statistical analysis is imperative. In this work,
we provide a comprehensive analysis on Doppler estimation. In particular, we
derive novel closed-form expressions for the probability of detection (PD) and
probability of false alarm (PFA). To this end, we consider the newly introduce
SP along with the CCRT. A comparison analysis between SP and the classic pulse
processing (PP) technique is also carried out. Numerical results and Monte-
Carlo simulations corroborate the validity of our expressions and show that
the SP–plus–CCRT technique helps to greatly reduce the PFA compared to
previous studies, thereby improving radar detection.

###### Index Terms:

Classic Chinese remainder theorem, robust Chinese remainder theorem, Doppler
frequency estimation, subpulse processing, probability of detection.

##  I Introduction

One important concern in modern pulsed radars is related to the Doppler
frequency estimation of fast-velocity targets. Due to the high target’s radial
velocity, ambiguous estimates are more likely to occur. More specifically,
ambiguous estimates appear whenever the target’s Doppler shift is greater than
the pulse repetition frequency (PRF) [1]. It seems obvious to think that
increasing the PRF will overcome this problem. However, if we are interested
in detecting targets located at long distances, then the PRF will be
restricted to a maximum value. Therefore, the choice of PRF is a trade-off
between range and Doppler requirements [2]. Fortunately, there are some
techniques that can resolve ambiguities, although at the cost of extra
measurement time and processing load. These techniques make use of multiples
PRFs [3, 4, 5, 6, 7]. The most known and used technique is the classic Chinese
remainder theorem (CCRT). The CCRT is a fast and accurate method to resolve
the unambiguous Doppler frequency. This is accomplished by solving a set of
congruences, formed by the estimated measurements of each PRF [7, 8, 9].
Nevertheless, in this method, the number of PRFs will not be sufficient to
resolve a certain quantity of targets. In general, $L$ PRFs are required to
successfully disambiguate $L-1$ targets. If the number of targets exceeds
$L-1$, then ghosts can appear.111Ghosts are false targets resulting from false
coincidences of Doppler-ambiguous or range-ambiguous data [4]. Unless
additional data (e.g., tracking information) is available, the radar has no
way of recognizing possible false detections [4]. Care must be taken in the
analysis and design since the number of PRFs and the number of targets to be
detected have a direct relationship.

Another issue concerning high-velocity targets is related to the signal-to-
noise ratio (SNR) loss. This occurs because the Doppler shift of fast-moving
targets will provoke a mismatch between the received signal and its replica
[2]. Consequently, the SNR after range compression may be drastically
reduced.222Range compression refers to the convolution operation between the
received signal and the replica of the transmitted signal [10]. Some radar
systems estimate and remove the Doppler shift prior to applying range
compression. Nonetheless, some residual or uncompensated Doppler typically
remains. This concern was partially alleviated in [11, 12]. Specifically, in
[11], the authors proposed a subpulse processing (SP), which proved to have a
higher Doppler tolerance,333Doppler tolerance refers to the degree of
degradation in the compressed response due to uncompensated Doppler [13].
increasing the ability to detect fast-moving targets. The shortcomings of SP
are computation time (critical for most radars), processing load, and poor
velocity resolution.

As stated before, the CCRT and SP have hardware and physical limitations when
it comes to estimating high target velocities. In practice, modern pulsed
radars take advantage of these two techniques so as to improve the system’s
capability to accurately detect the target’s true Doppler frequency. Since SP
the CCRT are affected by the presence of background noise, then a thorough
statistical analysis involving these two estimation techniques must be carried
out. Recently in [14], the authors proposed a novel expression for the
probability to correctly estimate the unambiguous Doppler frequency
considering the CCRT and the common pulse processing (PP) technique [2].
However, to the best of our knowledge, there is no performance analysis
considering the SP–plus–CCRT technique.

The main objective of this research is to combine the statistical analysis
conducted in [14] along with the newly introduced SP and the CCRT. To do so,
we adopt a stochastic model that suits our Doppler estimation techniques.
Then, we derive novel and closed-form expressions for: i) the probability to
correctly estimate the Doppler frequency, also called probability of detection
(PD); and ii) the probability to erroneously estimate the Doppler frequency,
also called probability of false alarm (PFA).

The remainder of this paper is organized as follows. Section II introduces
some key concepts to understand how the velocity estimation is performed.
Section III describes the system model. Section IV analyzes the Doppler
estimation using multiple PRFs; Section V discusses the representative
numerical results. Finally, Section VI concludes this paper.

In what follows, $(a)\text{mod}(b)$ denotes the remainder of the euclidean
division of $a$ by $b$; $\left|\cdot\right|$, absolute value;
$\lfloor\cdot\rfloor$, floor operation; $\text{round}(\cdot)$, rounding
operation; $\text{Pr}\left[\cdot\right]$, probability; $\mathbb{E}(\cdot)$,
expectation; $\text{Var}(\cdot)$, variance; $(\cdot)^{*}$, complex conjugate;
$\bigcap$, intersection of events; $\bigcup$, union of events;
$\mathcal{N}(\mu,\sigma^{2})$ denotes a Gaussian distribution with mean $\mu$
and variance $\sigma^{2}$; $\mathcal{N}_{c}(\mu,\sigma^{2})$ denotes a complex
Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$, and
$j=\sqrt{-1}$ is the imaginary unit.

##  II Preliminaries

In this section, we present a brief introduction about the PP and SP
techniques. Latter, we describe the basis to understand the CCRT algorithm.
Finally, we show how the combined technique SP–plus–CCRT works in order to
improve Doppler estimation.

###  II-A Pulse Processing

PP is the common technique employed by the radar to estimate the target
velocity and improve the SNR. In this processing technique, the radar
transmits a sequence of $M$ pulses during a coherent processing interval (CPI)
[15]. Then, range compression is performed on each pulse to improve the
radar’s range resolution. Finally, the discrete Fourier transform (DFT) is
applied along the slow-time samples to increase the SNR and to estimate the
target Doppler frequency [2]. These samples are collected at a rate equal to
the PRF. The maximum Doppler frequency shift that the radar manages to detect
using PP is $\Psi_{max}=\pm\text{PRF}/2$. If the target Doppler frequency,
$\mathit{f}_{d}$, exceeds this value, then the radar will deliver ambiguous
Doppler measurements. The Doppler frequency shift will be positive for closing
targets and negative for receding targets. The target velocity,
$\mathit{v}_{t}$, and its corresponding Doppler shift are related by the
following equation [16]:

$\displaystyle\mathit{f}_{d}=\frac{2\mathit{v}_{t}\mathit{f}_{R}}{\mathit{c}}=\frac{2\mathit{v}_{t}}{\lambda},$
(1)

where $\mathit{f}_{R}$ is the radar’s operation frequency, $\mathit{c}$ is the
speed of light, and $\lambda$ is the radar frequency.

###  II-B Subpulse Processing

SP improves Doppler tolerance by mitigating the loss in SNR caused by the
uncompensated Doppler shift of fast-moving targets [2, 11]. Moreover, SP is
used to overcome the problem of ambiguous Doppler measurements. The SP
algorithm runs as follows:

  1. 1.

First, the replica of the transmitted signal is divided into $N$ subpulses –
unlike PP that used the entire replica.

  2. 2.

Latter, range compression is carried out between each subpulse and the
received signal (cf. [11, 12] for a detailed discussion on this). The use of
shorter replicas will enhance the system’s Doppler tolerance [10], increasing
the detection capability of fast-moving targets. Of course, this process leads
to a reduction in the peak amplitude of the sub-compression response (by a
factor of $1/N$). Here, the slow-time samples are collected at a rate of
$\Phi=N/\tau$, where $\tau$ is the pulse width. It is important to emphasize
that PP and SP are performed simultaneously, that is, for each of the $M$
compressions, the radar carried out $N$ sub-compressions [12].

  3. 3.

Finally, the slow-time samples are coherently integrated to estimate the
target Doppler frequency and to “restore” the peak amplitude of the sub-
compression response.

The number of subpulses can be chosen as high as needed, as long as it is
taken into consideration that each additional subpulse requires an extra range
compression operation, increasing the computational load and computation time.
The maximum Doppler frequency shift that the radar can now manage to detect is
$\Phi_{max}=\pm N/2\tau$ [11]. Since $\Phi_{max}>\Psi_{max}$, SP provides a
higher frequency range of detection for fast-moving targets.

Computation time is critical for most radars and depends strongly on the
radar’s operation mode (e.g. tracking, searching or imaging), thereby limiting
the number of subpulses. Commonly, the number of subpulses is set between 5
and 10. However, this small number yields to a poor discretization in the
frequency domain and, consequently, producing inaccurate estimates.

Figure 1: Block diagram for Doppler estimation.

###  II-C Classic Chinese Remainder Theorem

The use of multiples PRFs is a common approach to resolve range and Doppler
ambiguities [3, 8, 17, 4]. In this work, we only focus on solving Doppler
ambiguities. Consider for the moment a target with Doppler shift
$\mathit{f}_{d}>\Psi_{max}$. In this scenario, the radar will detect the
target with an apparent Doppler shift, $\mathit{f}_{d_{ap}}$, that satisfies

$\displaystyle\mathit{f}_{d}=\mathit{f}_{d_{ap}}+n\text{PRF},$ (2)

where $n$ is some integer. It is convenience to express the target’s Doppler
shift $\mathit{f}_{d}$ in terms of its corresponding Doppler bin, $b_{d}$.
Thus, (2) becomes

$b_{d}=b_{ap}+nM,$ (3)

in which $b_{{ap}}\in\left\\{0,1,2,\ldots,M-1\right\\}$ is the apparent
Doppler bin, defined as

$\displaystyle b_{ap}=$
$\displaystyle\left\lfloor\absolutevalue{\frac{\mathit{f}_{d_{ap}}}{\Delta
D}}\right\rfloor,\ \ \ \ \ \ \ \ \ \ \ \ \mathit{f}_{d_{ap}}\geq 0$ (4)
$\displaystyle b_{ap}=$ $\displaystyle
M-\left\lfloor\absolutevalue{\frac{\mathit{f}_{d_{ap}}}{\Delta
D}}\right\rfloor,\ \ \ \ \ \ \mathit{f}_{d_{ap}}<0$ (5)

with $\Delta D=\text{PRF}/M$ being the Doppler bin spacing. Under this
scenario, the radar is incapable to detect the target’s true Doppler
frequency.

Now, suppose that we have $L$ PRFs. Then, the unambiguous target’s Doppler bin
must satisfies the following congruences:

$\displaystyle b_{d}\equiv b_{{ap}_{i}}+n_{i}M_{i},\ \ \ \ 1\leq i\leq L$ (6)

The CCRT states that if all PRFs are pairwisely coprimes, then the set of
congruences in (6) will have a unique solution given by [17, 4, 18]

$\displaystyle
b_{d}=\left(\sum_{i=1}^{L}b_{{ap}_{i}}\beta_{i}\right)\text{mod}\left(\Theta\right),$
(7)

where $\Theta=\prod_{i=1}^{L}M_{i}$, $\beta_{i}=b_{i}\Theta/\text{PRF}_{i}$,
and $b_{i}$ is the smaller integer which can be computed by solving the
following expression:

$\displaystyle\left(\frac{b_{i}\Theta}{M_{i}}\right)\text{mod}\left(M_{i}\right)=1.$
(8)

###  II-D Doppler Estimation

Fig. 1 depicts the entire block diagram for Doppler estimation. First, the
received signal passes through two types of independent range compression
blocks, one for PP and one for SP. This process is performed in sequence for
each pulse repetition interval (PRI). The outputs of both blocks are combined
and stored in memory to form a datacube [2]. (The datacube’s data is organized
by range, number of pulses, and number of subpulses.) More datacubes are
needed when using more than one PRF, as shown in Fig. 1. Next, a 2D-DFT block
is applied to each datacube to perform coherent integration. (The 2D-DFT block
is referred to as a two-dimensional DFT applied along with pulses and
subpulses.) Latter, the output of the 2D-DFT block is a matrix with the same
size containing the estimated Doppler shifts. This new matrix is referred to
as Doppler datacube. Finally, the CCRT is applied over the Doppler datacubes.
This process will be clarified in Section V by means of simulation.

Noise, jammer, and clutter are major concerns in all radar systems. In this
work, we consider the presence complex white Gaussian noise (CWGN). Thus, the
Doppler spectrum of fast moving targets will be compromised due to the
intrinsic characteristics of noise. For example, a high noise power could mask
small target returns, degrading radar performance. Even if the target return
is entirely deterministic, the combined signal (target–plus–noise) is a random
process and must be treated as such. Therefore, we need to assess the
statistics underlying Doppler analysis, but first, we need to come up with a
specific stochastic model that suits the requirements and design of our
radar’s estimation scheme. This is discussed in the next section.

##  III System Model

In this section, we propose a stochastic model that fits our signal processing
schemes. In addition, we describe the premises (hypotheses) used for Doppler
estimation.

According to Sections II-A and II-B, the collected signals in the slow-time
domain corresponding to PP and SP can be expressed, respectively, as

$\displaystyle g_{1}\left[m\right]=$ $\displaystyle
s_{1}\left[m\right]+w_{1}\left[m\right]$ $\displaystyle=$ $\displaystyle
a_{1}\exp\left(j2\pi\mathit{f}_{d}m/\text{PRF}\right)+w_{1}\left[m\right],\ \
0\leq m\leq M-1$ (9) $\displaystyle g_{2}\left[n\right]=$ $\displaystyle
s_{2}\left[n\right]+w_{2}\left[n\right]$ $\displaystyle=$ $\displaystyle
a_{2}\exp\left(j2\pi\mathit{f}_{d}n/\Phi\right)+w_{2}\left[n\right],\ \ \ \ \
\ \ 0\leq n\leq N-1$ (10)

where $s_{1}\left[m\right]$ and $s_{2}\left[n\right]$ are discrete complex
sine signals444In most systems, the radio frequency (RF) signal is mixed to
baseband prior to compression, and a coherent detector is used in the
downconversion process to form in-phase (I) and quadrature (Q) receive
channels, thereby creating a complex baseband signal. originated by changes in
the target position; $w_{1}\left[m\right]$ and $w_{2}\left[n\right]$ are
discrete additive complex Gaussian noises; and finally, $a_{1}$ and $a_{2}$
are the amplitudes at the output of the matched filters. Depending on the
target velocity, the output amplitudes $a_{1}$ and $a_{2}$ maybe be greatly
attenuated. However, the attenuation in $a_{2}$ is partially mitigated by the
use of SP. In particular, it follows that $a_{2}>a_{1}$ for high-velocity
targets [11]. Additionally, we define $2\sigma_{t_{1}}^{2}$ and
$2\sigma_{t_{2}}^{2}$ as the total mean powers – in the time domain – for
$w_{1}\left[m\right]$ and $w_{2}\left[n\right]$, respectively. As seen in
practice, and due to the stationary characteristic of noise, we have that
$\sigma_{t_{1}}^{2}=\sigma_{t_{2}}^{2}$ [19]. However, we will remain using
separate notations for $\sigma_{t_{1}}^{2}$ and $\sigma_{t_{2}}^{2}$ so as to
distinguish the noise power from PP and SP. Of course, these separate
notations will not alter, in any form, our performance analysis.

The SNR measured in the time domain considering PP and SP, can be expressed,
respectively, as

$\displaystyle\text{SNR}_{t_{1}}=$
$\displaystyle\frac{\left|a_{1}\right|^{2}}{2\sigma_{t_{1}}^{2}}$ (11)
$\displaystyle\text{SNR}_{t_{2}}=$
$\displaystyle\frac{\left|a_{2}/N\right|^{2}}{2\sigma_{t_{2}}^{2}}.$ (12)

Observe in (12) that the fact of dividing the replica into $N$ subpulses
causes a reduction in the SNR by a factor of $1/N^{2}$, as mentioned in
Section II-B.

The DFT is the primary operation to implement coherent integration. More
precisely, the DFT provides a mechanism to test multiple candidate frequencies
to maximize the integration gain [2]. The corresponding DFTs for (III) and
(III) are given, respectively, by

$\displaystyle G_{1}\left[k^{\prime}\right]\triangleq$ $\displaystyle\
\mathscr{F}\left\\{g_{1}\left[m\right]\right\\}$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{M-1}g_{1}\left[m\right]\exp\left(-j2\pi
k^{\prime}m/M\right)$ $\displaystyle=$ $\displaystyle
S_{1}\left[k^{\prime}\right]+W_{1}\left[k^{\prime}\right],\ \ \ 0\leq
k^{\prime}\leq M-1$ (13) $\displaystyle
G_{2}\left[l^{\prime}\right]\triangleq$ $\displaystyle\
\mathscr{F}\left\\{g_{2}\left[n\right]\right\\}$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{N-1}g_{2}\left[n\right]\exp\left(-j2\pi
l^{\prime}n/N\right)$ $\displaystyle=$ $\displaystyle
S_{2}\left[l^{\prime}\right]+W_{2}\left[l^{\prime}\right],\ \ \ \ \ \ 0\leq
l^{\prime}\leq N-1$ (14)

The SNR measured in the frequency domain considering PP and SP, are given,
respectively, by [2, Eq. (17.37)]

$\displaystyle\text{SNR}_{1}=$
$\displaystyle\frac{|Ma_{1}|^{2}}{2\sigma_{1}^{2}}$ (15)
$\displaystyle\text{SNR}_{2}=$
$\displaystyle\frac{|a_{2}|^{2}}{2\sigma_{2}^{2}},$ (16)

in which $\sigma_{1}^{2}=M\sigma_{t_{1}}^{2}$ and
$\sigma_{2}^{2}=N\sigma_{t_{1}}^{2}$ are half of the noise powers – in the
frequency domain – for $W_{1}\left[k^{\prime}\right]$ and
$W_{2}\left[l^{\prime}\right]$, respectively.

The Doppler estimates are based on the absolute values of
$G_{1}\left[k^{\prime}\right]$ and $G_{2}\left[l^{\prime}\right]$. That is,
(III) and (III) will provide estimates for $\mathit{f}_{d}$, say
$\hat{\mathit{f}}_{1}$ and $\hat{\mathit{f}}_{2}$, by searching $k^{\prime}$
and $l^{\prime}$, in which the absolute values of
$G_{1}\left[k^{\prime}\right]$ and $G_{2}\left[l^{\prime}\right]$ are maximum.
It is worth mentioning that if $\Psi_{max}<\mathit{f}_{d}$ and
$\Phi_{max}<\mathit{f}_{d}$, then $\hat{\mathit{f}}_{1}$ and
$\hat{\mathit{f}}_{2}$ will display ambiguous Doppler estimates.

Now, considering $L$ PRFs (say, $\text{PRF}_{1},\ldots,\text{PRF}_{L}$), we
can define the absolute values for $G_{1}\left[k^{\prime}\right]$ and
$G_{2}\left[l^{\prime}\right]$ at the $i$-th PRF, respectively, as

$\displaystyle H_{1,i}\left[k^{\prime}\right]\triangleq$
$\displaystyle|G_{1,i}\left[k^{\prime}\right]|\ \ \ \ \ 0\leq k^{\prime}\leq
M_{i}-1$ (17) $\displaystyle H_{2,i}\left[l^{\prime}\right]\triangleq$
$\displaystyle|G_{2,i}\left[l^{\prime}\right]|\ \ \ \ \ \ 0\leq l^{\prime}\leq
N_{i}-1$ (18)

where the subscript $i\in\\{1,\ldots,L\\}$ denotes the association to the
$i$-th PRF.

Herein, we assume that $G_{1,i}\left[k^{\prime}\right]$ is composed of
$M_{i}-1$ independent and identically distributed noise samples and one
target–plus–noise sample, denoted as $\mathcal{G}_{1,i}$. On the other hand,
$G_{2,i}\left[l^{\prime}\right]$ is composed of $N_{i}-1$ independent and
identically distributed noise samples and one combined sample, denoted as
$\mathcal{G}_{2,i}$. The target–plus–noise samples $\mathcal{G}_{1,i}$ and
$\mathcal{G}_{2,i}$ can be modeled, respectively, by [20, Eq. (1)]

$\displaystyle\mathcal{G}_{1,i}=$
$\displaystyle\sigma_{1,i}\left(\sqrt{1-\lambda_{1,i}^{2}}A_{1,i}+\lambda_{1,i}A_{0,i}\right)$
$\displaystyle+j\sigma_{1,i}\left(\sqrt{1-\lambda_{1,i}^{2}}B_{1,i}+\lambda_{1,i}B_{0,i}\right)$
(19) $\displaystyle\mathcal{G}_{2,i}=$
$\displaystyle\sigma_{2,i}\left(\sqrt{1-\lambda_{2,i}^{2}}A_{2,i}+\lambda_{2,i}A_{0,i}\right)$
$\displaystyle+j\sigma_{2,i}\left(\sqrt{1-\lambda_{2,i}^{2}}B_{2,i}+\lambda_{2,i}B_{0,i}\right),$
(20)

where $A_{p,i}$ and $B_{p,i}$ ($p=1,2$) are mutually independent random
variables (RVs) distributed as $\mathcal{N}(0,\frac{1}{2})$, and
$\lambda_{p,i},\in(0,1]$. Then, for any $p$ and $q$ ($q=1,2$), it follows that
$\mathbb{E}(A_{p,i}B_{q,i})=0$ and
$\mathbb{E}(A_{p,i}A_{q,i})=\mathbb{E}(B_{p,i}B_{q,i})=\frac{1}{2}\delta_{pq}$.
($\delta_{pq}=1$ if $p=q$, and $\delta_{pq}=0$ otherwise.) In addition,
$A_{0,i}$ and $B_{0,i}$ are mutually independent RVs distributed as
$\mathcal{N}(m_{\textbf{Re},i},\frac{1}{2})$ and
$\mathcal{N}(m_{\textbf{Im},i},\frac{1}{2})$, respectively. Thus,
$\mathcal{G}_{1,i}$ and $\mathcal{G}_{2,i}$ are non-zero mean complex Gaussian
RVs with probability density functions (PDFs) given, respectively, by
$\mathcal{N}_{c}(\lambda_{1,i}(m_{\textbf{Re},i}+jm_{\textbf{Im},i}),\sigma_{1,i}^{2})$
and
$\mathcal{N}_{c}(\lambda_{2,i}(m_{\textbf{Re},i}+jm_{\textbf{Im},i}),\sigma_{2,i}^{2})$.
The correlation coefficient between any pair of ($\mathcal{G}_{1,i}$,
$\mathcal{G}_{2,i}$), can be calculated as [20, Eq. (2)]

$\displaystyle\rho_{kl,i}\triangleq$
$\displaystyle\frac{\mathbb{E}(\mathcal{G}_{1,i}\mathcal{G}_{2,i}^{*})-\mathbb{E}(\mathcal{G}_{1,i})\mathbb{E}(\mathcal{G}_{2,i}^{*})}{\sqrt{\text{Var}(\mathcal{G}_{1,i})\text{Var}(\mathcal{G}_{2,i})}}$
$\displaystyle=$ $\displaystyle\lambda_{1,i}\lambda_{2,i}.$ (21)

This correlation exists because both PP and SP use the same received signal
when performing range compression [2]. Observe that the parameters
$\lambda_{1,i}^{2}$, $\lambda_{2,i}^{2}$, $m_{\textbf{Re},i}$ and
$m_{\textbf{Im},i}$ can be used to model the compressed responses
$\left|M_{i}a_{1,i}\right|^{2}$ and $\left|a_{2,i}\right|^{2}$. This can be
done by making the following substitutions:
$|M_{i}a_{1,i}|^{2}=\lambda_{1,i}^{2}(m_{\textbf{Re},i}^{2}+m_{\textbf{Im}}^{2})$
and
$|a_{2,i}|^{2}=\lambda_{2,i}^{2}(m_{\textbf{Re},i}^{2}+m_{\textbf{Im}}^{2})$.
On the other hand, $\lambda_{1,i}$ and $\lambda_{2,i}$ can be chosen to meet a
desire correlation coefficient.

By the above, it follows that $H_{1,i}\left[k^{\prime}\right]$ is composed of
$M_{i}-1$ Rayleigh distributed samples, denoted as $X_{k,i}$
$\left(k\in\left\\{1,2,\ldots,M_{i}-1\right\\}\right)$, and one Rice
distributed sample, denoted as $R_{1,i}$. Similarly,
$H_{2,i}\left[l^{\prime}\right]$ is composed of $N_{i}-1$ Rayleigh distributed
samples, denoted as $Y_{l,i}$
$\left(l\in\left\\{1,2,\ldots,N_{i}-1\right\\}\right)$, and one Rice
distributed sample, denoted as $R_{2,i}$. The PDFs of $X_{k,i}$ and $Y_{l,i}$
are given, respectively, by

$\displaystyle f_{X_{k,i}}(x_{k,i})=$
$\displaystyle\frac{x_{k,i}\exp\left(-\frac{x_{k,i}^{2}}{2\sigma_{k,i}^{2}}\right)}{\sigma_{k,i}}$
(22) $\displaystyle f_{Y_{l,i}}(y_{l,i})=$
$\displaystyle\frac{y_{l,i}\exp\left(-\frac{y_{l,i}^{2}}{2\sigma_{l,i}^{2}}\right)}{\sigma_{l,i}}.$
(23)

Moreover, since $R_{2,i}$ and $R_{2,i}$ bear a certain degree of correlation,
they are governed by a bivariate Rician distribution, given by [20, 21]

$\displaystyle\mathit{f}_{R_{1,i},R_{2,i}}$
$\displaystyle\left(r_{1,i},r_{2,i}|\mathcal{H}_{1}\right)=\int_{0}^{\infty}\exp\left(-t\xi_{i}\right)$
$\displaystyle\times\exp\left(-\textbf{m}_{i}\right)I_{0}\left(2\sqrt{\textbf{m}_{i}t}\right)\prod_{p=1}^{2}\frac{r_{p,i}}{\Omega_{p,i}^{2}}$
$\displaystyle\times\exp\left(-\frac{r_{p,i}^{2}}{2\Omega_{p,i}^{2}}\right)I_{0}\left(\frac{r_{p,i}\sqrt{t\sigma_{p,i}^{2}\lambda_{p,i}^{2}}}{\Omega_{p,i}^{2}}\right)\text{d}t,$
(24)

where $I_{0}(\cdot)$ is the modified Bessel function of the first kind and
order zero [22, Eq. (9.6.16)],
$\textbf{m}_{i}=m_{\textbf{Re},i}^{2}+m_{\textbf{Im},i}^{2}$, and

$\displaystyle\Omega_{p,i}^{2}$
$\displaystyle=\sigma_{p,i}^{2}\left(\frac{1-\lambda_{p,i}^{2}}{2}\right)$
(25a) $\displaystyle\xi_{i}$
$\displaystyle=1+\sum_{p=1}^{2}\frac{\sigma_{p,i}^{2}\lambda_{p,i}^{2}}{2\Omega_{p,i}^{2}}.$
(25b)

##  IV Doppler Analysis

In this section, we provide a comprehensive statistical analysis on Doppler
estimation. To do so, we derive the performance metrics for both SP and
SP–plus–CCRT.

###  IV-A SP Analysis

First, let us define the following events:

$\displaystyle\mathcal{A}_{k,i}=$
$\displaystyle\left\\{R_{1,i}>X_{k,i}\right\\}$ (26)
$\displaystyle\mathcal{B}_{l,i}=$
$\displaystyle\left\\{R_{2,i}>Y_{l,i}\right\\}$ (27)
$\displaystyle\mathcal{C}_{k,i}=$
$\displaystyle\left\\{X_{k,i}>R_{1,i}\right\\}$ (28)
$\displaystyle\mathcal{D}_{l,i}=$
$\displaystyle\left\\{Y_{l,i}>R_{2,i}\right\\}.$ (29)

######  Proposition I.

Let $\text{PD}_{i}$ be the probability of detection at the $i$-th PRF.
Specifically, $\text{PD}_{i}$ is defined as the probability that $R_{1,i}$ is
greater than $X_{k,i}$ and, simultaneously, that $R_{2,i}$ is greater than
$Y_{l,i}$, i.e.,

$\displaystyle\text{PD}_{i}\triangleq\text{Pr}$
$\displaystyle\left[\left(\bigcap_{k=1}^{M_{i}-1}\mathcal{A}_{k,i}\right)\bigcap\left(\bigcap_{l=1}^{N_{i}-1}\mathcal{B}_{l,i}\right)\right].$
(30)

Then, from (22)–(III), (30) can be expressed in closed-form as

$\displaystyle\text{PD}_{i}=$
$\displaystyle\sum_{k=0}^{M_{i}-1}\sum_{l=0}^{N_{i}-1}\left(\begin{array}[]{c}M_{i}-1\\\
k\\\ \end{array}\right)\left(\begin{array}[]{c}N_{i}-1\\\ l\\\
\end{array}\right)$ (35)
$\displaystyle\times\frac{(-1)^{-k-l+M_{i}+N_{i}}\mathcal{V}_{i}(k,l)}{\mathcal{U}_{i}(k,l)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{U}_{i}(k,l)}\right),$
(36)

wherein $\mathcal{U}_{i}(k,l)$ and $\mathcal{V}_{i}(k,l)$ are auxiliary
functions defined, respectively, as

$\displaystyle\mathcal{U}_{i}(k,l)=$ $\displaystyle\
\xi_{i}-\frac{\xi_{i}\lambda_{1,i}^{2}\sigma_{1,i}^{4}}{2\Omega_{1,i}^{2}\left(\Omega_{1,i}^{2}(k-M_{i}+1)-\sigma_{1,i}^{2}\right)}$
$\displaystyle-\frac{\xi_{i}\lambda_{2,i}^{2}\sigma_{2,i}^{4}}{2\Omega_{2,i}^{2}\left(\Omega_{2,i}^{2}(l-N_{i}+1)-\sigma_{2,i}^{2}\right)}$
(37a) $\displaystyle\mathcal{V}_{i}(k,l)=$ $\displaystyle\
\frac{\sigma_{1,i}^{2}}{\left(\Omega_{1,i}^{2}(-k+M_{i}-1)+\sigma_{1,i}^{2}\right)}$
$\displaystyle\times\frac{\sigma_{2,i}^{2}}{\left(\Omega_{2,i}^{2}(-l+N_{i}-1)+\sigma_{2,i}^{2}\right)}.$
(37b)

###### Proof.

See Appendix A. ∎

######  Corollary I.

Let $\text{PFA}_{i}$ be the probability of false alarm at the $i$-th PRF. More
precisely, $\text{PFA}_{i}$ is defined as the probability that at least one of
$X_{k,i}$ is greater than $R_{1,i}$ and, simultaneously, that at least one of
$Y_{l,i}$ is greater than $R_{2,i}$, i.e.,

$\displaystyle\text{PFA}_{i}\triangleq\text{Pr}$
$\displaystyle\left[\underset{k=1}{\overset{M_{i}-1}{\bigcup}}\underset{l=1}{\overset{N_{i}-1}{\bigcup}}\left(\mathcal{C}_{k,i}\bigcap\mathcal{D}_{l,i}\right)\right].$
(38)

Then, from (22)–(III), (38) can be written in closed-form as in (Corollary I),
shown at the top of the next page, where $\mathcal{P}_{i}\left(k,l\right)$ and
$\mathcal{Q}_{i}\left(k,l\right)$ are auxiliary functions defined,
respectively, by

$\displaystyle\mathcal{P}_{i}\left(k,l\right)=$ $\displaystyle\
\xi_{i}-\frac{\lambda_{1,i}^{2}\sigma^{4}_{1,i}}{2\Omega_{1,i}^{2}\left(k\
\Omega_{1,i}^{2}+\sigma^{2}_{1,i}\right)}$
$\displaystyle-\frac{\lambda_{2,i}^{2}\sigma^{4}_{2,i}}{2\Omega_{2,i}^{2}\left(l\
\Omega_{2,i}^{2}+\sigma^{2}_{2,i}\right)}$ (39a)
$\displaystyle\mathcal{Q}_{i}\left(k,l\right)=$ $\displaystyle\
\frac{\sigma^{2}_{1,i}\sigma^{2}_{2,i}}{\left(k\
\Omega_{1,i}^{2}+\sigma^{2}_{1,i}\right)\left(l\
\Omega_{2,i}^{2}+\sigma^{2}_{2,i}\right)}.$ (39b)

$\displaystyle\textit{PFA}_{i}=$
$\displaystyle\frac{\left(M_{i}-1\right)\left(N_{i}-1\right)\mathcal{Q}_{i}\left(1,1\right)}{\mathcal{P}_{i}\left(1,1\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(1,1\right)}\right)-\binom{M_{i}-1}{2}\binom{N_{i}-1}{2}\frac{\mathcal{Q}_{i}\left(2,2\right)}{\mathcal{P}_{i}\left(2,2\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(2,2\right)}\right)+\ldots$
$\displaystyle+(-1)^{M_{i}-N_{i}-1}\frac{\mathcal{Q}_{i}\left(M_{i}-1,N_{i}-1\right)}{\mathcal{P}_{i}\left(M_{i}-1,N_{i}-1\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(M_{i}-1,N_{i}-1\right)}\right)$
(40)

###### Proof.

See Appendix B. ∎

It is worth mentioning that (35) and (Corollary I) are novel and original
contributions of this work, derived in closed-form even though (III) is given
in integral form.

###  IV-B SP–Plus–CCRT Analysis

Similar to [14], we assume that each individual pulse on each sweep results in
an independent random value for the target returns.

Now, using (35) and taking into account the $\mathcal{M}$–of–$L$ detection
criterion,555Instead of detecting a target on the basis of at least one
detection in $L$ tries, system designers often require that some number
$\mathcal{M}$ or more detections be required in $L$ tries before a target
detection is accepted [2]. the probability of detection for the combined
technique SP–plus–CCR can be calculated as follows [23]

$\displaystyle\text{PD}_{\text{CCRT}}\triangleq$
$\displaystyle\sum_{l=\mathcal{M}}^{L}\sum_{\mathcal{E}\in\mathcal{F}_{l}}\left\\{\left(\prod_{i\in\mathcal{E}}\text{PD}_{i}\right)\left(\prod_{j\in\mathcal{E}^{c}}\left(1-\text{PD}_{j}\right)\right)\right\\},$
(41)

where $\mathcal{F}_{l}$ is the set of all subsets of $l$ integers that can be
selected from $\left\\{1,2,\ldots,L\right\\}$, and $\mathcal{E}^{c}$ is the
complement of $\mathcal{E}$. For example, if $l=2$ and $L=3$, then
$\mathcal{F}_{2}=\left\\{\left\\{1,2\right\\},\left\\{1,3\right\\},\left\\{2,3\right\\}\right\\}$,
and $\mathcal{E}^{c}=\left\\{1,2,\ldots,L\right\\}\backslash\mathcal{E}$.

On the other hand, the probability of false alarm for the combined technique
SP–plus-CCRT can be calculated as [23]

$\displaystyle\text{PFA}_{\text{CCRT}}\triangleq$
$\displaystyle\sum_{l=\mathcal{M}}^{L}\sum_{\mathcal{E}\in\mathcal{F}_{l}}\left\\{\left(\prod_{i\in\mathcal{E}}\text{PFA}_{i}\right)\left(\prod_{j\in\mathcal{E}^{c}}\left(1-\text{PFA}_{j}\right)\right)\right\\}.$
(42)

For the case where $\mathcal{M}=L$, (41) and (42) reduce, respectively, to

$\displaystyle\text{PD}_{\text{CCRT}}=$
$\displaystyle\prod_{i=1}^{L}\text{PD}_{i}$ (43)
$\displaystyle\text{PFA}_{\text{CCRT}}=$
$\displaystyle\prod_{i=1}^{L}\text{PFA}_{i}.$ (44)

##  V Numerical Results

Figure 2: $\text{PD}_{i}$ vs $\text{SNR}_{1}$ using $N_{i}=8$,
$\lambda_{1,i}=0.5$, $\lambda_{2,i}=0.99$, and different values of $M_{i}$
($i\in\left\\{1,2,3,4\right\\}$). Figure 3: $\text{PMD}_{i}$ vs
$\text{SNR}_{1}$ using $N_{i}=8$, $\lambda_{1,i}=0.5$, $\lambda_{2,i}=0.99$,
and different values of $M_{i}$ ($i\in\left\\{1,2,3,4\right\\}$).

---
(a–1) $\text{PRF}_{1}=1700$ [Hz], $M_{1}=11$

---
(a–2) $\text{PRF}_{2}=1900$ [Hz], $M_{2}=13$

---
(a–3) $\text{PRF}_{3}=2100$ [Hz], $M_{3}=17$

---
(a–4) $\text{PRF}_{4}=2300$ [Hz], $M_{4}=19$

---
(b–1) $\text{PRF}_{1}=1700$ [Hz], $N_{1}=8$

---
(b–2) $\text{PRF}_{2}=1900$ [Hz], $N_{2}=8$

---
(b–3) $\text{PRF}_{3}=2100$ [Hz], $N_{3}=8$

---
(b–4) $\text{PRF}_{4}=2300$ [Hz], $N_{4}=8$

---
(c–1) $\text{PRF}_{1}=1700$ [Hz]

---
(c–2) $\text{PRF}_{2}=1900$ [Hz]

---
(c–3) $\text{PRF}_{3}=2100$ [Hz]

---
(c–4) $\text{PRF}_{4}=2300$ [Hz]

Figure 4: Doppler estimation.

In this section, we illustrate through Fig. 4 how the Doppler estimation
process is carried out. Latter, we validate our derived expressions by means
of Monte-Carlo simulations666The number of realizations in Monte-Carlo
simulations was set to $10^{6}$.. To do so, we make use of the following radar
setup: $\text{PRF}_{1}=700$ [Hz], $\text{PRF}_{2}=1100$ [Hz],
$\text{PRF}_{3}=1300$ [Hz], $\text{PRF}_{4}=1700$ [Hz], $L=\mathcal{M}=4$,
$\mathit{f}_{R}=6\ [\text{GHz}]$, $\tau=25\ [\mu s]$, $\lambda=0.05\
[\text{m}]$, $M_{1}=11$, $M_{2}=13$, $M_{3}=17$, $M_{4}=19$, and $N_{i}=8\
\forall i\in\left\\{1,2,3,4\right\\}$. In addition, we consider a linear
frequency-modulated pulse with bandwidth $B=2\ [\text{MHz}]$.

Fig. 4 illustrates the output data after the 2D-DFT blocks. In this simulation
example, we placed a target at an initial range of 10 [Km], traveling with a
constant velocity of $\mathit{v}_{t}=900$ [m/s] in the opposite direction of
the radar (i.e., the target is receding). Fig. 4(a) shows the normalized
output data – Velocity vs Range – using PP. Observe that in all 4 scenarios,
the target at 10 [Km] is unlikely to be detected due to the high loss in SNR.
On the other hand, Fig. 4(b) shows the normalized output data – Velocity vs
Range – using SP. Observe that the loss in SNR is partially mitigated by means
of SP. Therefore, the target located at 10 [Km] can now be easily be detect
without further processing. At last, Fig. 4(c) shows the combined pulse and
subpulse information. Note in Fig. 4(c) that SP provides a better intuition
about the target location, but due to its poor discretization, it is not
sufficient to determine the exact velocity. Conversely, PP provides a better
discretization but, unfortunately, its velocity estimation is more likely to
be ambiguous. Thus, by combining SP and the CCRT, we provide the system a high
capability to unfold the target’s true velocity.

Fig. 2 shows $\text{PD}_{i}$ versus $\text{SNR}_{1}$ using different values of
$M_{i}$. Note how radar performance improves as $M_{i}$ increases, requiring a
lower SNR for a given PD. This is because when increasing $M_{i}$, we are, in
fact, increasing the compressed response of PP by means of coherent
integration. In particular, for a fixed $\text{SNR}_{1}=10$ [dB], we obtain
the following probabilities of detection: $\text{PD}_{1}=0.66$ for $M_{1}=7$;
$\text{PD}_{2}=0.78$ for $M_{2}=11$; $\text{PD}_{3}=0.85$ for $M_{3}=13$; and
$\text{PD}_{4}=0.93$ for $M_{4}=17$. Also, observe that for the high and
medium SNR regime, our derived expression matches perfectly the PD of [14, Eq.
(28)]. Nevertheless, there is a small difference in the PD for the low SNR
regime. This occurs because if the compressed response of PP is less than the
background noise, then the intersection probability in (30) will be less than
the probability of $\bigcap_{k=1}^{M_{i}-1}\mathcal{A}_{k,i}$ . For example,
given $\text{SNR}_{1}=4$ [dB] and $M_{1}=7$, we obtain $\text{PD}_{1}=0.15$
with our proposed SP–plus–CCRT technique, and $\text{PD}_{1}=0.18$ with [14,
Eq. (28)]. However, this small reduction in the PD is compensated by a greater
reduction in the PFA, as shall be seen next.

Figure 5: $\text{PD}_{\text{CCRT}}$ vs $\text{SNR}_{1}$ using $N_{i}=8$,
$\lambda_{1,i}=0.5$, $\lambda_{2,i}=0.99$, $\mathcal{M}=4$, and different
values of $M_{i}$ ($i\in\left\\{1,2,3,4\right\\}$).  Figure 6:
$\text{PMD}_{\text{CCRT}}$ vs $\text{SNR}_{1}$ using $N_{i}=8$,
$\lambda_{1,i}=0.5$, $\lambda_{2,i}=0.99$, $\mathcal{M}=4$, and different
values of $M_{i}$ ($i\in\left\\{1,2,3,4\right\\}$).

Fig. 3 shows $\text{PFA}_{i}$ versus $\text{SNR}_{1}$ using different values
for $M_{i}$. Observe how $\text{PFA}_{i}$ decreases as $M_{i}$ increases. This
occurs because as we increase $M_{i}$, the received target echo becomes
stronger compared to the noise background. For example, for a fixed
$\text{SNR}_{1}=5$ [dB], we obtain the following probabilities of false alarm:
$\text{PFA}_{1}=0.83$ for $M_{1}=7$; $\text{PFA}_{2}=0.77$ for $M_{2}=11$;
$\text{PFA}_{3}=0.73$ for $M_{3}=13$; and $\text{PFA}_{4}=0.60$ for
$M_{4}=17$. More interesting, observe how $\text{PFA}_{i}$ decays rapidly
compared to [14]. This difference in $\text{PFA}_{i}$ is because intuitively
SP acts as a backup detection process. That is, since the compressed response
of SP is greater of the PP response (for high-velocity targets), then the
probability in (38) is lower than the probability of
$\underset{k=1}{\overset{M_{i}-1}{\bigcup}}\mathcal{C}_{k,i}$ . For example,
using the classic PP technique [14], we obtain the following probabilities of
false alarm: $\text{PFA}_{1}=0.96$ for $M_{1}=7$; $\text{PFA}_{2}=0.97$ for
$M_{2}=11$; $\text{PFA}_{3}=0.98$ for $M_{3}=13$; and $\text{PFA}_{4}=0.99$
for $M_{3}=17$.

Finally, Figs. 5 and 6 show $\text{PD}_{\text{CCRT}}$ and
$\text{PFA}_{\text{CCRT}}$ versus $\text{SNR}_{1}$, respectively. Observe in
Fig. 5, the perfect agreement between (41) and [14, Eq. (29)]. Hence, in this
case, we have no advantage when using SP–plus–CCRT. On the other hand, observe
in Fig. 6, the high difference in the PFA between of (42) and that in [14]. In
this case, the use of SP–plus–CCRT improves radar performance by considerably
reducing the false alarms. For instance, for given $\text{SNR}_{1}=2$ [dB], we
obtain probabilities of $\text{PFA}_{\text{CCRT}}=0.94$ using PP–plus–CCRT,
and $\text{PFA}_{\text{CCRT}}=0.54$ using SP–plus–CCRT.

##  VI Conclusion

In this work, we provided a thorough statistical analysis on Doppler
estimation when both SP and the CCRT were employed. To do so, we derived novel
and closed-form expressions for the PD and PFA. Moreover, a comparison
analysis between our proposed SP–plus–CCRT technique and the classic
PP–plus–CCRT was carried out. Numerical results and Monte-Carlo simulations
corroborated the validity of our expressions and showed that the PFA when
using SP–plus–CCRT technique was greatly reduced compared to [14], thereby
enhancing radar detection.

##  Appendix A Proof of Proposition I

Applying [24, Eq. (5.48)] and using the fact that $X_{k,i}$ and $Y_{l,i}$ are
independent RVs, (30) can be rewritten as follows:

$\displaystyle\text{PD}_{i}=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\left(\prod_{k=1}^{M_{i}-1}\text{Pr}\left[X_{k,i}<r_{1,i}|R_{1,i}=r_{1,i}\right]\right)$
$\displaystyle\times\left(\prod_{l=1}^{N_{i}-1}\text{Pr}\left[Y_{l,i}<r_{2,i}|R_{2,i}=r_{2,i}\right]\right)$
$\displaystyle\times\mathit{f}_{R_{1,i},R_{2,i}}(r_{1,i},r_{2,i})\
\text{d}r_{1,i}\ \text{d}r_{2,i}.$ (45)

Now, with the aid of [24, Eq. (4.11)] and taking into account that $X_{k,i}$
and $Y_{l,i}$ are identically distributed RVs, yields

$\displaystyle\text{PD}_{i}=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\left(\int_{0}^{r_{1,i}}f_{X_{1,i}}(x_{1,i})\
\text{d}x_{1,i}\right)^{M_{i}-1}$
$\displaystyle\times\left(\int_{0}^{r_{2,i}}f_{Y_{1,i}}(y_{1,i})\
\text{d}y_{1,i}\right)^{N_{i}-1}$
$\displaystyle\times\mathit{f}_{R_{1,i},R_{2,i}}(r_{1,i},r_{2,i})\
\text{d}r_{1,i}\ \text{d}r_{2,i}.$ (46)

Replacing (22)–(III) in (A), we obtain

$\displaystyle\text{PD}_{i}=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\underbrace{\left(\int_{0}^{r_{1,i}}\frac{x_{1,i}\exp\left(-\frac{x_{1,i}^{2}}{2\sigma_{1,i}^{2}}\right)}{\sigma_{1,i}}\text{d}x_{1,i}\right)^{M_{i}-1}}_{\triangleq\
\mathcal{I}_{1}}$
$\displaystyle\times\underbrace{\left(\int_{0}^{r_{2,i}}\frac{y_{1,i}\exp\left(-\frac{y_{1,i}^{2}}{2\sigma_{2,i}^{2}}\right)}{\sigma_{2,i}}\text{d}y_{1,i}\right)^{N_{i}-1}}_{\triangleq\
\mathcal{I}_{2}}$
$\displaystyle\times\int_{0}^{\infty}\exp(-\xi_{i}t)\exp\left(-\textbf{m}_{i}\right)I_{0}\left(2\sqrt{\textbf{m}_{i}t}\right)$
$\displaystyle\times\prod_{p=1}^{2}\frac{r_{p,i}}{\Omega_{p,i}^{2}}\exp\left(-\frac{r_{p,i}^{2}}{2\Omega_{p,i}^{2}}\right)$
$\displaystyle\times
I_{0}\left(\frac{r_{p,i}\sqrt{t\sigma_{p,i}^{2}\lambda_{p,i}^{2}}}{\Omega_{p,i}^{2}}\right)\text{d}t\
\text{d}r_{1,i}\ \text{d}r_{2,i}.$ (47)

In order to solve (A), we must first evaluate $\mathcal{I}_{1}$ and
$\mathcal{I}_{2}$. In particular, $\mathcal{I}_{1}$ can be calculated as
follows:

$\displaystyle\mathcal{I}_{1}$
$\displaystyle\overset{(a)}{=}\left(1-\exp\left(-\frac{r_{1,i}^{2}}{2\sigma_{1,i}^{2}}\right)\right)^{M_{i}-1}$
$\displaystyle\overset{(b)}{=}\sum_{k=0}^{M_{i}-1}\left(\begin{array}[]{c}M_{i}-1\\\
k\\\
\end{array}\right)\left(-\exp\left(-\frac{r_{1,i}^{2}}{2\sigma_{1,i}^{2}}\right)\right)^{M_{i}-1-k},$
(50)

where in step (a), we have developed the inner integral; and in step (b), we
have used the binomial Theorem [24].

Using a similar approach to that used in (A), $\mathcal{I}_{2}$ can be
calculated as

$\displaystyle\mathcal{I}_{2}=\sum_{l=0}^{N_{i}-1}\left(\begin{array}[]{c}N_{i}-1\\\
l\\\
\end{array}\right)\left(-\exp\left(-\frac{r_{2,i}^{2}}{2\sigma_{2,i}^{2}}\right)\right)^{N_{i}-1-l}.$
(53)

Inserting (A) and (53) in (A), followed by changing the order of
integration777The change in the order of integration was performed without
loss of generality since (22), (23) and (III) are non-negative real functions
[25]. and along with minor manipulations, we obtain (58), displayed at the top
of the next page.

$\displaystyle\text{PD}_{i}=$
$\displaystyle\sum_{k=0}^{M_{i}-1}\sum_{l=0}^{N_{i}-1}1^{k+l}\left(\begin{array}[]{c}M_{i}-1\\\
k\\\ \end{array}\right)\left(\begin{array}[]{c}N_{i}-1\\\ l\\\
\end{array}\right)\int_{0}^{\infty}\exp(-\xi_{i}t)\exp\left(-\textbf{m}_{i}\right)I_{0}\left(2\sqrt{\textbf{m}_{i}t}\right)$
(58)
$\displaystyle\times\underbrace{\int_{0}^{\infty}\left(-\exp\left(-\frac{r_{1,i}^{2}}{2\sigma_{1,i}^{2}}\right)\right)^{M_{i}-1-k}\frac{r_{1,i}}{\Omega_{1,i}^{2}}\exp\left(-\frac{r_{1,i}^{2}}{2\Omega_{1,i}^{2}}\right)I_{0}\left(\frac{r_{1,i}\sqrt{t\sigma_{1,i}^{2}\lambda_{1,i}^{2}}}{\Omega_{1,i}^{2}}\right)\text{d}r_{1,i}}_{\triangleq\
\mathcal{I}_{3}}$
$\displaystyle\times\underbrace{\int_{0}^{\infty}\left(-\exp\left(-\frac{r_{2,i}^{2}}{2\sigma_{2,i}^{2}}\right)\right)^{N_{i}-1-l}\frac{r_{2,i}}{\Omega_{2,i}^{2}}\exp\left(-\frac{r_{2,i}^{2}}{2\Omega_{2,i}^{2}}\right)I_{0}\left(\frac{r_{2,i}\sqrt{t\sigma_{2,i}^{2}\lambda_{2,i}^{2}}}{\Omega_{2,i}^{2}}\right)\text{d}r_{2,i}}_{\triangleq\
\mathcal{I}_{4}}\text{d}t.$ (59)

Now, it remains to find $\mathcal{I}_{3}$ and $\mathcal{I}_{4}$. More
precisely, $\mathcal{I}_{3}$ can be computed as

$\displaystyle\mathcal{I}_{3}\overset{(a)}{=}$
$\displaystyle\int_{0}^{\infty}\left(-\exp\left(-\frac{r_{1,i}^{2}}{2\sigma_{1,i}^{2}}\right)\right)^{M_{i}-1-k}\frac{r_{1,i}}{\Omega_{1,i}^{2}}$
$\displaystyle\times\exp\left(-\frac{r_{1,i}^{2}}{2\Omega_{1,i}^{2}}\right)\sum_{q=0}^{\infty}\frac{\left(\frac{r_{1,i}\sqrt{t\lambda_{1,i}^{2}\sigma_{1,i}^{2}}}{2\Omega_{1,i}^{2}}\right)^{2q}}{q!\
\Gamma(q+1)}\text{d}r_{1,i}$ $\displaystyle\overset{(b)}{=}$
$\displaystyle\frac{(-1)^{-k+M_{i}+1}}{\Omega_{1,i}^{2}\left(\frac{-k+M_{i}-1}{\sigma_{1,i}^{2}}+\frac{1}{\Omega_{1,i}^{2}}\right)}$
$\displaystyle\times\sum_{q=0}^{\infty}\frac{\left(\frac{t\lambda_{1,i}^{2}\sigma_{1,i}^{4}}{2\Omega_{1,i}^{2}\left(\Omega_{1,i}^{2}(-k+M_{i}-1)+\sigma_{1,i}^{2}\right)}\right)^{q}}{q!}$
$\displaystyle\overset{(c)}{=}$
$\displaystyle\frac{(-1)^{-k+M_{i}+1}}{\Omega_{1,i}^{2}\left(\frac{-k+M_{i}-1}{\sigma_{1,i}^{2}}+\frac{1}{\Omega_{1,i}^{2}}\right)}$
$\displaystyle\times\exp\left(\frac{t\lambda_{1,i}^{2}\sigma_{1,i}^{4}}{2\Omega_{1,i}^{2}\left(\Omega_{1,i}^{2}(-k+M_{i}-1)+\sigma_{1,i}^{2}\right)}\right),$
(60)

where in step (a), we have used the series representation of the modified
Bessel function of the first kind and order zero [26, Eq. (03.02.02.0001.01)];
in step (b), we have solved the integral by first changing the order of
integration; finally, in step (c), we have used [26, Eq. (01.03.06.0002.01)]
and performed some algebraic manipulations.

In like manner as in (A), $\mathcal{I}_{4}$ can be computed as

$\displaystyle\mathcal{I}_{4}=$
$\displaystyle\frac{(-1)^{-l+N_{i}+1}}{\Omega_{2,i}^{2}\left(\frac{-l+N_{i}-1}{\sigma_{2,i}^{2}}+\frac{1}{\Omega_{2,i}^{2}}\right)}$
$\displaystyle\times\exp\left(\frac{t\lambda_{2,i}^{2}\sigma_{2,i}^{4}}{2\Omega_{2,i}^{2}\left(\Omega_{2,i}^{2}(-l+N_{i}-1)+\sigma_{2,i}^{2}\right)}\right).$
(61)

Now, replacing (A) and (A) in (58), we obtain

$\displaystyle\text{PD}_{i}=$
$\displaystyle\sum_{k=0}^{M_{i}-1}\sum_{l=0}^{N_{i}-1}\left(\begin{array}[]{c}M_{i}-1\\\
k\\\ \end{array}\right)\left(\begin{array}[]{c}N_{i}-1\\\ l\\\
\end{array}\right)$ (66)
$\displaystyle\times\exp\left(-\textbf{m}_{i}\right)\left(\frac{\sigma_{1,i}^{2}(-1)^{-k+M_{i}+1}}{\Omega_{1,i}^{2}(-k+M_{i}-1)+\sigma_{1,i}^{2}}\right)$
$\displaystyle\times\left(\frac{\sigma_{2,i}^{2}(-1)^{-l+N_{i}+1}}{\Omega_{2,i}^{2}(-l+N_{i}-1)+\sigma_{2,i}^{2}}\right)$
$\displaystyle\times\int_{0}^{\infty}\exp(-\xi_{i}t)I_{0}\left(2\sqrt{\textbf{m}_{i}t}\right)$
$\displaystyle\times\exp\left(\frac{t\lambda_{1,i}^{2}\sigma_{1,i}^{4}}{2\Omega_{1,i}^{2}\left(\Omega_{1,i}^{2}(-k+M_{i}-1)+\sigma_{1,i}^{2}\right)}\right)$
$\displaystyle\times\exp\left(\frac{t\lambda_{2,i}^{2}\sigma_{2,i}^{4}}{2\Omega_{2,i}^{2}\left(\Omega_{2,i}^{2}(-l+N_{i}-1)+\sigma_{2,i}^{2}\right)}\right)\text{d}t.$
(67)

Finally, using the following identity [27, Eq. (1.11.2.4)]

$\int_{0}^{\infty}\exp(tb)I_{0}(\sqrt{t}a)\
\text{d}t=-\frac{\exp\left(-\frac{a^{2}}{4b}\right)}{b},$ (68)

and after performing some minor simplifications, we can express (66) in
closed-form as in (35), which completes the proof.

##  Appendix B Proof of Corollary I

By making use of [24, Coroll. 6], we can express (38) as

$\displaystyle\text{PFA}_{i}=\sum_{k=1}^{M_{i}-1}\sum_{l=1}^{N_{i}-1}\text{Pr}\left[\mathcal{C}_{k,i}\bigcap\mathcal{D}_{l,i}\right]$
$\displaystyle\
-\underset{k<p,l<q}{\sum_{k=1}^{M_{i}-1}\sum_{l=1}^{N_{i}-1}\sum_{p=2}^{M_{i}-1}\sum_{q=2}^{N_{i}-1}}\text{Pr}\left[\mathcal{C}_{k,i}\bigcap\mathcal{D}_{l,i}\bigcap\mathcal{C}_{p,i}\bigcap\mathcal{D}_{q,i}\right]+\ldots$
$\displaystyle\
+(-1)^{M_{i}-N_{i}-1}\text{Pr}\left[\mathcal{C}_{1,i}\bigcap\mathcal{D}_{1,i}\bigcap\ldots\bigcap\mathcal{C}_{M_{i}-1,i}\bigcap\mathcal{D}_{N_{i}-1,i}\right].$
(69)

Now, we need to find the event probabilities. First, let us derive the last
event probability of (B), that is,

Pr
$\displaystyle\left[\mathcal{C}_{1,i}\bigcap\mathcal{D}_{1,i}\bigcap\ldots\bigcap\mathcal{C}_{M_{i}-1,i}\bigcap\mathcal{D}_{N_{i}-1,i}\right]$
$\displaystyle\overset{a}{=}$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\left(\prod_{k=1}^{M_{i}-1}\text{Pr}\left[X_{k,i}>r_{1,i}|R_{1,i}=r_{1,i}\right]\right)$
$\displaystyle\times\left(\prod_{l=1}^{N_{i}-1}\text{Pr}\left[Y_{l,i}>r_{2,i}|R_{2,i}=r_{2,i}\right]\right)$
$\displaystyle\times\mathit{f}_{R_{1,i},R_{2,i}}(r_{1,i},r_{2,i})\
\text{d}r_{1,i}\ \text{d}r_{2,i}$ $\displaystyle\overset{b}{=}$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\left(\int_{r_{1,i}}^{\infty}f_{X_{1,i}}(x_{1,i})\
\text{d}x_{1,i}\right)^{M_{i}-1}$
$\displaystyle\times\left(\int_{r_{2,i}}^{\infty}f_{Y_{1,i}}(y_{1,i})\
\text{d}y_{1,i}\right)^{N_{i}-1}$
$\displaystyle\times\mathit{f}_{R_{1,i},R_{2,i}}(r_{1,i},r_{2,i})\
\text{d}r_{1,i}\ \text{d}r_{2,i},$ (70)

where in step (a) we have used [24, Eq. (5.48)]; and in step (b) we have used
[24, Eq. (4.11)] along with the fact that $X_{k,i}$ and $Y_{l,i}$ are
identically distributed RVs.

Replacing (22)–(III) in (B), yields

Pr
$\displaystyle\left[\mathcal{C}_{1,i}\bigcap\mathcal{D}_{1,i}\bigcap\ldots\bigcap\mathcal{C}_{M_{i}-1,i}\bigcap\mathcal{D}_{N_{i}-1,i}\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\underbrace{\left(\int_{r_{1,i}}^{\infty}\frac{x_{1,i}\exp\left(-\frac{x_{1,i}^{2}}{2\sigma_{1,i}^{2}}\right)}{\sigma_{1,i}}\text{d}x_{1,i}\right)^{M_{i}-1}}_{\triangleq\
\mathcal{I}_{5}}$
$\displaystyle\times\underbrace{\left(\int_{r_{2,i}}^{\infty}\frac{y_{1,i}\exp\left(-\frac{y_{1,i}^{2}}{2\sigma_{2,i}^{2}}\right)}{\sigma_{2,i}}\text{d}y_{1,i}\right)^{N_{i}-1}}_{\triangleq\
\mathcal{I}_{6}}$
$\displaystyle\times\int_{0}^{\infty}\exp(-\xi_{i}t)\exp\left(-\textbf{m}_{i}\right)I_{0}\left(2\sqrt{\textbf{m}_{i}t}\right)$
$\displaystyle\times\prod_{p=1}^{2}\frac{r_{p,i}}{\Omega_{p,i}^{2}}\exp\left(-\frac{r_{p,i}^{2}}{2\Omega_{p,i}^{2}}\right)$
$\displaystyle\times
I_{0}\left(\frac{r_{p,i}\sqrt{t\sigma_{p,i}^{2}\lambda_{p,i}^{2}}}{\Omega_{p,i}^{2}}\right)\text{d}t\
\text{d}r_{1,i}\ \text{d}r_{2,i}.$ (71)

After some mathematical manipulations, $\mathcal{I}_{5}$ and $\mathcal{I}_{6}$
can be calculated, respectively, as

$\displaystyle\mathcal{I}_{5}=$
$\displaystyle\exp\left(-\frac{r_{1,i}^{2}(M_{i}-1)}{2\sigma_{1,i}^{2}}\right)$
(72) $\displaystyle\mathcal{I}_{6}=$
$\displaystyle\exp\left(-\frac{r_{2,i}^{2}(N_{i}-1)}{2\sigma_{2,i}^{2}}\right).$
(73)

$\displaystyle\text{PFA}_{i}=$
$\displaystyle\binom{M_{i}-1}{1}\binom{N_{i}-1}{1}\frac{\mathcal{Q}_{i}\left(1,1\right)}{\mathcal{P}_{i}\left(1,1\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(1,1\right)}\right)-\binom{M_{i}-1}{2}\binom{N_{i}-1}{2}\frac{\mathcal{Q}_{i}\left(2,2\right)}{\mathcal{P}_{i}\left(2,2\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(2,2\right)}\right)+\ldots$
$\displaystyle+(-1)^{M_{i}-N_{i}-1}\binom{M_{i}-1}{M_{i}-1}\binom{N_{i}-1}{N_{i}-1}\frac{\mathcal{Q}_{i}\left(M_{i}-1,N_{i}-1\right)}{\mathcal{P}_{i}\left(M_{i}-1,N_{i}-1\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(M_{i}-1,N_{i}-1\right)}\right)$
(74)

Now, replacing (72) and (73) in (B), and after solving remaining three
integrals by applying the same procedure as in (66), we obtain

Pr
$\displaystyle\left[\mathcal{C}_{1,i}\bigcap\mathcal{D}_{1,i}\bigcap\ldots\bigcap\mathcal{C}_{M_{i}-1,i}\bigcap\mathcal{D}_{N_{i}-1,i}\right]$
$\displaystyle=\frac{\mathcal{Q}_{i}\left(M_{i}-1,N_{i}-1\right)}{\mathcal{P}_{i}\left(M_{i}-1,N_{i}-1\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(M_{i}-1,N_{i}-1\right)}\right),$
(75)

where $\mathcal{P}_{i}\left(k,l\right)$ and $\mathcal{Q}_{i}\left(k,l\right)$
are auxiliary functions defined in (39), and the parameters
$k\in\left\\{1,2,\ldots,M_{i}-1\right\\}$ and
$l\in\left\\{1,2,\ldots,N_{i}-1\right\\}$ denote the number of events for
$\mathcal{C}_{k,i}$ and $\mathcal{D}_{l,i}$, respectively. Thus, the remaining
event probabilities in (B) can be easily obtained by a proper choice of the
parameters $k$ and $l$. For example, for $k=1$ and $l=3$, we obtain

Pr
$\displaystyle\left[\mathcal{C}_{1,i}\bigcap\mathcal{D}_{1,i}\bigcap\mathcal{D}_{2,i}\bigcap\mathcal{D}_{3,i}\right]$
$\displaystyle=\frac{\mathcal{Q}_{i}\left(1,3\right)}{\mathcal{P}_{i}\left(1,3\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(1,3\right)}\right),$
(76)

whereas for $k=3$ and $l=2$, we have

Pr
$\displaystyle\left[\mathcal{C}_{1,i}\bigcap\mathcal{D}_{1,i}\bigcap\mathcal{C}_{2,i}\bigcap\mathcal{D}_{2,i}\bigcap\mathcal{C}_{3,i}\right]$
$\displaystyle=\frac{\mathcal{Q}_{i}\left(3,2\right)}{\mathcal{P}_{i}\left(3,2\right)}\exp\left(-\textbf{m}_{i}+\frac{\textbf{m}_{i}}{\mathcal{P}_{i}\left(3,2\right)}\right).$
(77)

Later, with the aid of (B) and after some algebraic manipulations, we can
rewrite (B) as in (B), displayed at the top of the next page. Finally, and
after minor simplifications, (B) reduces to (Corollary I), which completes the
proof.

## References

  * [1] G. Morris and L. Harkness, _Airborne Pulsed Doppler Radar_ , 2nd ed. Norwood, MA, USA: Artech House, 1996.
  * [2] M. A. Richards, J. Scheer, W. A. Holm, and W. L. Melvin, _Principles of Modern Radar: Basic Principles_ , 1st ed. West Perth, WA, Australia: SciTech, 2010.
  * [3] G. V. Trunk, “Range resolution of targets using automatic detectors,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. AES-14, no. 5, pp. 750–755, Sept. 1978.
  * [4] S. A. Hovanessian, “An algorithm for calculation of range in a multiple PRF radar,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. AES-12, no. 2, pp. 287–290, Mar. 1976.
  * [5] X.-G. Xia and G. Wang, “Phase unwrapping and a robust chinese remainder theorem,” _IEEE Signal Process. Lett._ , vol. 14, no. 4, pp. 247–250, Apr. 2007.
  * [6] X. Li, H. Liang, and X. Xia, “A robust chinese remainder theorem with its applications in frequency estimation from undersampled waveforms,” _IEEE Trans. Signal Process._ , vol. 57, no. 11, pp. 4314–4322, Nov. 2009\.
  * [7] W. Wang and X. Xia, “A closed-form robust chinese remainder theorem and its performance analysis,” _IEEE Trans. Signal Process._ , vol. 58, no. 11, pp. 5655–5666, Nov. 2010.
  * [8] G. V. Trunk and W. M. Kim, “Ambiguity resolution of multiple targets using pulse-Doppler waveforms,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 30, no. 4, pp. 1130–1137, Oct. 1994.
  * [9] F. D. A. García, A. S. Guerreiro, G. R. L. Tejerina, J. C. S. Santos Filho, G. Fraidenraich, M. D. Yacoub, M. A. M. Miranda, and H. Cioqueta, “Probability of detection for unambiguous doppler frequencies in pulsed radars using the chinese remainder theorem and subpulse processing,” in _Proc. 53rd Asilomar Conference on Signals, Systems, and Computers_ , Pacific Grove, CA, USA, Nov. 2019, pp. 138–142.
  * [10] M. I. Skolnik, _Introduction to Radar Systems_ , 3rd ed. Ney York, NY, USA: McGraw-Hill, 2001.
  * [11] G. Beltrao, L. Pralon, M. Menezes, P. Vyplavin, B. Pompeo, and M. Pralon, “Subpulse processing for long range surveillance noise radars,” in _Proc. International Conference on Radar Systems (Radar 2017)_ , Belfast, UK, Oct. 2017, pp. 1–4.
  * [12] A. Barreto, L. Pralon, B. Pompeo, G. Beltrao, and M. Pralon, “FPGA design and implementation of a real-time subpulse processing architecture for noise radars,” in _Proc. 2019 International Radar Conference (RADAR)_ , Toulon, France, Sept. 2019, pp. 1–6.
  * [13] D. S. Doviak and R. J. Zrnic, _Doppler Radar and Weather Observations_ , 2nd ed. San Diego, CA, USA: Academic Press, 2001.
  * [14] B. Silva and G. Fraidenraich, “Performance analysis of the classic and robust chinese remainder theorems in pulsed doppler radars,” _IEEE Trans. Signal Process._ , vol. 66, no. 18, pp. 4898–4903, Sept. 2018.
  * [15] M. A. Richards, _Fundamentals of Radar Signal Processing_ , 2nd ed. Ney York, NY, USA: McGraw-Hill, 2014.
  * [16] D. K. Barton, _Radar Equations for Modern Radar_ , 1st ed. Massachusetts, MA, USA: Artech House, 2013.
  * [17] G. Trunk and S. Brockett, “Range and velocity ambiguity resolution,” in _Proc. Record IEEE Nat. Radar Conf._ , Lynnfield, MA, USA, Apr. 1993, pp. 146–149.
  * [18] A. Ferrari, C. Berenguer, and G. Alengrin, “Doppler ambiguity resolution using multiple PRF,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 33, no. 3, pp. 738–751, Jul. 1997.
  * [19] A. Papoulis, _Probability, Random Variables, and Stochastic Processes_ , 4th ed. Ney York, NY, USA: McGraw-Hill, 2002.
  * [20] N. C. Beaulieu and K. T. Hemachandra, “Novel representations for the bivariate rician distribution,” _IEEE Trans. Commun._ , vol. 59, no. 11, pp. 2951–2954, Nov. 2011.
  * [21] A. Behnad, N. C. Beaulieu, and K. T. Hemachandra, “Correction to “Novel representations for the bivariate rician distribution”,” _IEEE Trans. Commun._ , vol. 60, no. 6, pp. 1486–1486, Jun. 2012.
  * [22] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_. Washington, DC: US Dept. of Commerce: National Bureau of Standards, 1972.
  * [23] Y. H. Wang, “On the number of successes in independent trials,” _Statistica Sinica_ , vol. 3, no. 2, pp. 295–312, 1993.
  * [24] A. Leon-Garcia, _Probability and Random Processes for Electrical Engineering_ , 3rd ed. New Jersey, NJ, USA: Pearson Prentice Hall, 1994.
  * [25] H. Friedman, “A consistent Fubini-Tonelli theorem for nonmeasurable functions,” _Illinois J. Math._ , vol. 24, no. 3, pp. 390–395, 1980.
  * [26] Wolfram Research, Inc. (2018), _Wolfram Research_ , Accessed: Sept. 19, 2018\. [Online]. Available: http://functions.wolfram.com
  * [27] A. P. Prudnikov, Y. A. Bryčkov, and O. I. Maričev, _Integral and Series: Vol. 2_ , 2nd ed., Fizmatlit, Ed. Moscow, Russia: Fizmatlit, 1992.