LPI Radar Signal Detection Based on Autocorrelation Function and Wigner-Ville Distribution

In this study, PSK-modulated CW signals and FMCW signals are analyzed. In the analysis, signals are used in the analytic form. By forming the analytic signal, the negative frequency components of a real signal may be eliminated from the signal representation without losing information. This elimination provides sampling below the Nyquist rate and avoids interference terms by the interaction of positive and negative components in quadratic TFDs [35], the analytic form of the transmitted radar waveform expressed as

$s(t)=A e^{j\left(2 \pi f_c t+\emptyset_k\right)}$    (1)

The in phase (I) and quadrature (Q) component of s(t) is defined as

$I=A \cos \left(2 \pi f_c t+\emptyset_k\right)$    (2)

$Q=A \sin \left(2 \pi f_c t+\emptyset_k\right)$    (3)

where, A is the amplitude, $f_c$ is the carrier frequency of the transmit signal and $\emptyset_k$ is the phase of the signal. For FMCW signal, phase is constant and carrier frequency changes in a symmetric triangular waveform such as:

for the up-ramp of the signal: $f_c-\frac{\Delta F}{2}+\frac{\Delta F}{t_m} t$      (4)

for the down-ramp of the signal: $f_c+\frac{\Delta F}{2}+\frac{\Delta F}{t_m} t$         (5)

where, $\Delta F$ is the transmit modulation bandwidth and $t_m$ is the modulation period.

For PSK signals, the carrier frequency is constant. Instead of frequency, phase values are changed by specific codes such as binary phase codes, polyphase codes, or polytime codes. PSK modulated signal's phase changes after every subcode period $t_b$. Consequently, the code period T and code rate $C_R$ are defined as:

$T=N_c t_b$    (6)

$C_R=\frac{1}{N_c t_b}$   (7)

where, $N_c$ represents the number of phase codes and represents the code length as well. The duration of the transmitted signal is$P \times T$, where P  is the number of periods.

If cpp  is the number of cycles of the carrier frequency per subcode, the bandwidth of the transmitted signal is

$B=\frac{f_c}{c p p}=\frac{1}{t_b} \mathrm{~Hz}$       (8)

For polyphase-encoded waveforms, using $c p p=1$ results in the maximum bandwidth achievable with any carrier frequency. The shorter the duration of the subcode period, the greater the bandwidth. On the other hand, for polytime codes $t_b$ is not constant. Polytime code sequences use fixed phase states with varying periods at each phase state. Polytime coded waveforms bandwidth is measured due to the shortest phase change duration. Polytime coding has the advantage that arbitrary time bandwidth waveforms can be generated with only a few phase states. An increase in code lengths of PSK waveform does not affect bandwidth but increases the processing gain. The processing gain is equal to the compression ratio of the signal. For both the frequency modulation and phase modulation LPI radar, the transmitted CW signal is coded with a reference signal to spread the transmitted energy in frequency [3]. The codes with zero sidelobes in periodic ACF (PACF) are called perfect codes. The low sidelobe levels of the ACF and PACF help quantify the LPI waveform's ability to detect targets without interfering with sidelobe targets. High sidelobes might be caused nearby targets to hide in a sidelobe [36]. Frank, P1, P3, and P4 codes have perfect PACF [3, 37]. Polyphase is the most frequently used code because of easy digital implementation, versatility, high range resolution, and Doppler tolerance [13]. The specific codes used in LPI radar waveform PSK modulation is given below.

2.1 Binary barker codes

A binary Barker sequence is a finite length, discrete time sequence $A=\left\{a_0, a_1, \ldots, a_n \right\}$ of $+1^{\prime}$ s and $-1^{\prime}$ s of length $n \geq 2$, such that the aperiodic autocorrelation coefficients (or sidelobes):

$r_k=\sum_{j=1}^{n-k} a_j a_{j+k}$     (9)

Satisfies $\left|r_k\right| \leq 1$ for $k \neq 0$ and similarly $r_{-k}=r_k$. The Binary Barker sequences correspond the phase values of $0^{\circ}$ or $180^{\circ}$  used for BPSK modulation.

Binary Barker codes are only known for lengths $N_c=2,3,4,5,7,11$ and 13  given in Table 1 [38].

Table 1. Binary codes

2.2 Polyphase barker codes

Polyphase code sequences correspond to more than two phase values. Their alphabet size $N_c>2$. The idea of increasing the size of the alphabet for constructing Barker sequences was completed in 1965 by Golomb and Scholz. Golomb and Scholz introduced the concept of Generalized Barker codes, defined as a multiphase sequence [39]. Allowing any phase values (nonbinary) can lead to lower sidelobes. However, the outermost sidelobe is always 1 [37]. Consider the generalized Barker sequences $\left\{a_j\right\}$ of finite length n where the terms $a_j$ are allowed complex numbers of absolute value of 1 where the correlation is now Hermetian dot product.

$r_k=\sum_{j=1}^{n-k} a_j a_{j+k}^*$    (10)

where, $a^*$ represents the complex conjugate of a . If all the sidelobes of the ACF of any polyphase sequence are bounded by

$|r(k)| \leq 1$ for $k \neq 0$   (11)

Then the sequence is called generalized Barker sequence or polyphase Barker sequence. The binary Barker sequences can be regarded as a special case of polyphase Barker sequences. If the sequence elements are taken from an alphabet of size M , consisting of the $M^{t h}$ roots of unity,

$a_m=\exp \left\{2 \pi i \cdot \frac{m}{M}\right\}=: \exp \left(i \emptyset_m\right) \ 0 \leq m \leq M-1$     (12)

The sequence is alternatively named an M-phase Barker sequence [40]. Polyphase Barker signal is difficult to intercept and correctly classify using the ES system [2]. In many applications, the phases are restricted to values that are the kth roots of unity (e.g, $k^{}= 2$ gives the original Barker codes or $k^{}= 6$  for six phase/sextic Barker codes). As for polyphase Barker codes, systematic methods to construct polyphase Barker sequences have not been found yet. Searches for generalized Barker sequences with no restriction on the values of the sequence phases have been carried out using numerical optimization [37].

2.3 Frank code

Frank code derived from the phase history of a linearly frequency stepped pulse, using M  frequency steps and M  samples per frequency. The number of code elements is equal to $M^2$. The elements of the original Frank code $a_{(n-1) M+k}=\exp \left(j \emptyset_{n, k}\right)$, for $1 \leq n \leq M$ and $1 \leq k \leq M$. The phase values are expressed as:

$\emptyset_{n, k}=\frac{2 \pi}{M} (n-1)(k-1), \quad n=1,2, \ldots, M$ and $k=1,2, \ldots, M$   (13)

where, n  is the number of the sample in a given frequency and k  is the number of the frequency.

2.4 P1, P2, P3 and P4 code

P1 and P2 signals are derived from step approximation to a linear frequency modulated waveform. In these codes using M frequency steps and M samples per frequency are produced. P3 and P4 signals are derived from linear frequency modulated waveform. Phase equations of P-codes given in Table 2.

Table 2.  P1-P2-P3-P4 Codes and phase equation

2.5 Polytime codes

$T 1(n)$ and $T 2(n)$ types of polytime code can be generated from the stepped frequency model, $T 3(n)$ and $T 4(n)$ types are approximations to a linear frequency modulation.

Phase versus time relation of $T 1(n)$ and $T 3(n)$ code $\emptyset_{T 1}(t), \emptyset_{T 3}(t)$, are expressed as follows, respectively,

$\emptyset_{T 1}(t)=mod \left\{\frac{2 \pi}{n} I N T\left[(k t-j T) \frac{j n}{T}\right], 2 \pi\right\}$    (14)

$\emptyset_{T 3}(t)=mod \left\{\frac{2 \pi}{n} I N T\left[\frac{n \Delta F t^2}{2 t_m}\right], 2 \pi\right\}$   (15)

In Eq. (14), $j=0,1,2, \ldots, k-1$ is the segment number in the stepped frequency waveform, k is the number of the segments in the T1 code sequence, t is the time, T is the overall code duration, n is the number of phase states in the code sequence. In Eq. (15), $t_m$ is the modulation period, and $\Delta F$ is modulation bandwidth. As more phase states are added to the polytime sequence, the agreement in time sidelobe behavior improves [3].

An example of LPI radar signal which is PSK modulated by Frank code where $N_c=64, \mathrm{f}_{\mathrm{c}}=1 \mathrm{kHz}$, cpp=1, $\mathrm{SNR}=-6 \mathrm{~dB}$ shown in Figure 1.

12.png

Figure 1. Frank coded waveform, SNR = -6 dB

4.1 Autocorrelation based signal detection

This section presents spread spectrum DSSS and radar signals detection methods. DSSS signal detection by autocorrelation is subcategorized as single-channel autocorrelation, dual-channel autocorrelation, and compounded autocorrelation [18]. Detection results from the exploitation of cyclostationary of DSSS signals and similar communication transmissions by correlating received signals. Detection is performed by analyzing correlation peaks. The fluctuation of correlation estimators method is presented to detect DSSS signals [43]. This method is the same as compounded autocorrelation mentioned in the study of [18]. The idea behind this technique is that a DSSS signal's autocorrelation is like a noise's autocorrelation. On the other hand, the fluctuations of estimators are different. The detection is based on the peaks of fluctuations.

The autocorrelation sequence of a periodic signal is periodic too. This property is used for the detection and parameter estimation of LPI radar signals [19]. The autocorrelation function is also used to filter out the noise. In [44], for LPI radar signal detection, autocorrelation processing of the signal is used to filter out the noise, and radar signal classification is performed based on features obtained from the autocorrelation function and directed graphical model.

Periodic heritage and noise reduction property of the ACF explained below:

When $x(t)$ is the observable signal at the receiver output is given in Eq. (16), the autocorrelation function of $x(t)$ is

$\begin{aligned} & R_{x x}(\tau)=E[x(t) x(t-\tau)] \\ & =E\{[s(t)+n(t)][s(t-\tau)+n(t-\tau)]\} \\ & =E[s(t) s(t-\tau)]+E[n(t) n(t-\tau)] \\ & +E[s(t) n(t-\tau)] \\ & +E[n(t) s(t-\tau)] \\ & =R_{s s}(\tau)+R_{n n}(\tau)+R_{s n}(\tau)+R_{n s}(\tau)\end{aligned}$     (17)

$R_{s n}(\tau)=R_{n s}(\tau)=0$ since the noise is not related to the s(t). Hence

$R_{x x}(\tau)=R_{s s}(\tau)+R_{n n}(\tau)$    （18）

For zero mean AWGN with wide bandwidth, the autocorrelation function $R_{n n}(\tau)$ mainly affects the nearby position $\tau=0$.

$R_{n n}(\tau)=\sigma^2 e^{-\alpha|\tau|}$ for $|\tau| \gg 1$   （19）

Hence $R_{x x}(\tau) \approx R_{s s}(\tau)$ [19, 34, 44].

4.2 Time frequency analysis based detection

Time-frequency $(t, f)$ methods can be grouped into amplitude (linear) "atomic" decompositions and "energy" quadratic (or bilinear) distributions. The linear class represents a signal as a linear combination of a family of elementary signals (called atoms) that are well localized in time and frequency. The linear class includes the short-time Fourier transform (STFT) and the wavelet transform (WT). The bilinear class includes quadratic TFDs (QTFDs) that distribute the energy of the signal between two variables, namely time and frequency (or time and scale). Wigner Ville Distribution (WVD) is a bilinear distribution. All other QTFDs are filtered versions of the WVD. Windowed versions of the WVD are the pseudo WVD (PWVD), and smoothed pseudo WVD used to improve TF representations [35]. The WVD of a signal $s(t)$ denoted by $W_z(t, f)$, is defined as:

$W_z(t, f)=\underset{\tau \rightarrow f}{\mathcal{F}}\left\{z\left(t+\frac{\tau}{2}\right) z^*\left(t-\frac{\tau}{2}\right)\right\}$   （20)

where, $z(t)$ is complex and defined as the analytic associate of $s(t)$. So, the WVD of signal $s(t)$ is given as:

$W_z(t, f)=\int_{-\infty}^{+\infty} z\left(t+\frac{\tau}{2}\right) z^*\left(t-\frac{\tau}{2}\right) e^{-j 2 \pi f \tau} d \tau$     (21)

WVD is also defined as the Fourier transform of the instantaneous autocorrelation function of a signal $z(t, f)$ [45]. The marginal integration of TFD of WVD over frequency and time gives instantaneous power and energy spectrum of the signal respectively that are given Eqns. (22) and (23).

$\int_{-\infty}^{+\infty} W_z(t, f) d f=|z(t)|^2$    (22)

$\int_{-\infty}^{+\infty} W_z(t, f) d f=|z(f)|^2$   (23)

The integration of the WVD over the entire $(t, f)$ plane yields the signal energy $E_z$.

$\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} W_z(t, f) d t d f=E_z$   (24)

The discrete WVD expressed as:

$W(l, w)=2 \sum_{-\infty}^{+\infty} z(l+n) z^*(l-n) e^{-j 2 w n}$   (25)

Windowing the data results in the pseudo-WVD and is defined by:

$\begin{array}{r}W(l, w)=2 \sum_{n=-N+1}^{N+1} z(l+n) z^*(l-n) \ldots \\ \ldots w (n) w(-n) e^{-j 2 w n}\end{array}$    (26)

where, $w(n)$ is a length $2 N-1$  real window function with $w(0)=1$. The two drawbacks of WVD are negative values and cross terms. Because of its bilinear structure WVD has cross terms appear in the middle of every pair of signal components. WVD can have large negative values. Most TFDS used in practical applications can take negative values as they are not necessarily positive definite, so they do not represent and instantaneous energy spectral density at any arbitrary $(t, f)$ location [35]. One of the features that can be extracted from the Time-frequency distributions is joint time frequency moments. The two-dimensional moment of order $(p+q)$ of a time frequency distribution $T F D(t, f)$ is defined as:

$M_{p q}=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} t^p f^q T F D(t, f) d t d f$    (27)

The central moments of $T F D(t, f)$ is defined as:

$\mu_{p q}=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty}(t-\bar{t})^p(f-\bar{f})^q T F D(t, f) d t d f$   (28)

where, $\bar{t}$ and $\bar{f}$ are defined as:

$\bar{t}=\frac{M_{10}}{M_{00}}, \quad \bar{f}=\frac{M_{01}}{M_{00}}$   (29)

4.3 Energy detector

The energy detector computes the energy in the received data and compares it to a threshold. The energy detector (ED) for an unknown deterministic signal $x[n]$ is expresses as

$T(x)=\sum_{n=0}^{N-1} x^2[n]>\gamma$     (30)

where, $x[n]$ is signal pluse noise under hypothesis H1, only noise signal under hypothesis H0. $T(x)$ is called the test statistics. ED decides H1 if $T(x)$ is greater than the threshold. Under H0, the received signal is chi-squared distributed which probability density function given as [28, 46].

$p(x)= \begin{cases}\frac{1}{2^{\frac{v}{2}-} \Gamma\left(\frac{v}{2}\right)} x^{\frac{v}{2}-1} \exp \left(-\frac{1}{2} x\right) & x>0 \\ 0 & x<0\end{cases}$      (31)

where, v is degree of freedom.

The threshold determined by a given probability of false alarm rate $\left(P_{F A}\right)$.

It is given that if the number of sampled signals (N) is large enough and the channel is nonfading [47], the chi-square distribution with N degrees of freedom can be considered as a normal distribution by using the central limit theorem, and $P_{F A}$ written as

$P_{F A}=Q\left(\frac{\gamma-M N \sigma^2}{\sigma^2 \sqrt{2 M N}}\right)$    (32)

where, $\sigma^2$ is the noise variance and M is antenna number.