Parameter Estimation Based on Sigmoid Transform in Wideband Bistatic MIMO Radar System under Impulsive Noise Environment

Abstract

Since second-order statistics-based methods rely heavily on Gaussianity assumption and fractional lower-order statistics-based methods depend on a priori knowledge of non-Gaussian noise, there remains a void in wideband bistatic multiple-input/multiple-output (MIMO) radar systems under impulsive noise. In this paper, a novel method based on Sigmoid transform was used to estimate target parameters, which do not need a priori knowledge of the noise in an impulsive noise environment. Firstly, a novel wideband ambiguity function, termed Sigmoid wideband ambiguity function (Sigmoid-WBAF), is proposed to estimate the Doppler stretch and time delay by searching the peak of the Sigmoid-WBAF. A novel Sigmoid correlation function is proposed. Furthermore, a new MUSIC algorithm based on the Sigmoid correlation function (Sigmoid-MUSIC) is proposed to estimate the direction-of-departure (DOD) and direction-of-arrival (DOA). Then, the boundness of the Sigmoid-WBAF to the symmetric alpha stable ($S\alpha S$) noise, the feasibility analysis of the Sigmoid-WBAF, and complexity analysis of the Sigmoid-WBAF and Sigmoid-MUSIC are presented to evaluate the performance of the proposed method. In addition, the Cramér–Rao bound for parameter estimation was derived and computed in closed form, which shows that better performance was achieved. Simulation results and theoretical analyses are presented to verify the effectiveness of the proposed method.

1. Introduction

Recently, multiple-input/multiple-output (MIMO) systems attracted more and more attention [1,2,3]. There is rapidly growing literature concerning target parameter estimation in MIMO radar systems. Generally, target parameter estimation algorithms can be used in both narrowband and wideband MIMO radar systems. In narrowband MIMO radar systems, Yoe et al. proposed a computationally efficient method using properly designed projection filters to estimate the direction-of-departure (DOD) and direction-of-arrival (DOA) of targets [4]. Chen et al. proposed a method based on ESPRIT using the rotational factor produced by time-delay sampling to estimate the Doppler, DOD, and DOA [5]. Yao et.al. proposed a novel estimation method based on a non-uniform array configuration to estimate the DOD and DOA parameters, and analyzed their practical identifiability [6]. Jiang et al. proposed a joint estimation algorithm based on the canonical correlation decomposition (CCD) and exploits the shift-invariance properties in the Kronecker product structure of each column of the various steering matrices in unknown correlated noise [7]. Some improved algorithms based on MUSIC and ESPRIT were also proposed to estimate DOD and DOA parameters [8,9,10]. These methods were shown to yield good performance in terms of their parameters. In many applications, however, it is not appropriate to approximate the signal using a narrowband model.

In wideband radar systems, the echo from a wideband signal often contains a Doppler stretch (DS), not merely a Doppler shift, which results in parameter estimation difficulty [11,12]. To determine the range and relative velocity of a target, an accurate estimation of these parameters is crucial. Ma et al. proposed two novel methods for incoherent broadband chirp DOA estimation (BCD-I) and coherent broadband chirp DOA estimation (BCD-C) [13]. You et al. and Xu et al. presented DOA estimation methods of wideband signal based on spectral correlation signal subspace fitting (SC-SSF) [14,15]. Chen et al. proposed a method of DOA estimation for wideband coherent acoustic sources based on coherent signal subspace [16], and Yoon [17] presented a DOA estimation method based on the orthogonality of projected subspaces to estimate DOA. However, these methods did not estimate the Doppler shift and time delay, which are also very crucial for the determination of the range and velocity of the target in a wideband bistatic MIMO radar. At present, we seldom find the study of a joint estimation for the Doppler shift, time delay, DOD, and DOA in a wideband bistatic MIMO radar [18], which needs to be studied deeply because of its usage for target tracking and target localization.

In general, broadband chirp signals are frequently used in sonar and radar systems. The broadband chirp or linear frequency-modulated (LFM) signals are known to be insensitive to the Doppler of echoes and have better properties of low probability of interception. Therefore, in this paper, we utilized an LFM signal as a transmitted signal to study the wideband signal model parameter estimation in bistatic MIMO radars.

Until now, in most parameter estimation methods for array signal processing, additive noise is assumed to be Gaussian. Studies and experimental measurements showed that broad classes of noise such as underwater acoustic noise, atmospheric noise, multiuser interference, and radar clutters in real-world applications are non-Gaussian, primarily owing to impulsive phenomena [19,20,21,22,23,24,25,26]. Taking these scenarios into account, it is inappropriate to model the noise as Gaussian noise. Researchers studied this impulsive nature and showed that symmetric alpha stable ($S\alpha S$) processes are better models for impulsive noise than Gaussian processes. Conventional algorithms based on second-order statistics degenerate severely in an impulsive noise environment [12,21].

To reduce the alpha stable distribution noise interference, many parameter estimation algorithms based on the fractional lower-order statistics (FLOS) were proposed [18,19,20,21,22,23,24,25,26,27]. However, these algorithms have some limitations: (1) the noise characteristic exponent must be estimated to ensure $1\le p<\alpha$ or $0<p<\alpha /2$, where $p$ is the fractional lower order of moments and $\alpha$ is the characteristic exponent of the impulsive noise; (2) these methods cannot accurately estimate the parameters if there is no a priori knowledge of the characteristic exponent. Furthermore, the performance of these algorithms can degrade seriously and even become invalid when the fractional lower-order moment value is not appropriate. To handle this problem, a novel wideband ambiguity function based on Sigmoid transform, referred to as Sigmoid wideband ambiguity function (Sigmoid-WBAF), is proposed to estimate the Doppler stretch (DS) and time delay (TD) in this paper. Then, a novel correlation function based on the Sigmoid transform, referred to as Sigmoid correlation, is defined. Furthermore, a novel MUSIC method based on the Sigmoid correlation, referred as Sigmoid-MUSIC, is proposed to estimate the DOD and DOA.

This paper is organized as follows: Section 2 presents a signal model of a wideband bistatic MIMO radar system. Section 3 describes a novel Sigmoid wideband ambiguity function and a novel Sigmoid correlation. In Section 4, the Sigmoid-WBAF and Sigmoid-MUSIC methods are used to estimate the target parameters. In Section 5, the boundness of the Sigmoid-WBAF to the symmetric alpha stable ($S\alpha S$) noise, and the feasibility analysis of the Sigmoid-WBAF and the Cramér–Rao bound for parameter estimation are presented to evaluate the performance of the proposed method. In Section 6, the performance of the parameter estimation algorithm is studied through extensive numerical simulations. Finally, conclusions are drawn in Section 6.

2. Signal Model and Noise Model

2.1. Signal Model and Bandpass Matched Filter

2.1.1. Signal Model

Consider a bistatic MIMO radar system with a $Q$-transmitter and an $N$-receiver, both of which are half-wavelength spaced uniform linear arrays. Assume that there are $L$ uncorrelated targets in the far field, where the targets are located, and the signal wavefront is considered as a plane wave. The described bistatic MIMO radar system is illustrated as a two-dimensional (2D) configuration in Figure 1.

Each antenna transmits the linear frequency modulation signal $x_q(t)${$q=1,\dots ,Q$}. These transmitted signals can be reflected by L targets at positions $\left({\phi }_l,{\theta }_l\right)$ {$l=1,\dots ,L$}, where ${\phi }_l$ denotes the DOD and ${\theta }_l$ denotes the DOA. In wideband sonar and radar systems, the echo from a wideband signal often contains a Doppler stretch (DS), in addition to a Doppler shift, due to the moving of the target. Thus, the received signal $y_n\left(t\right)$ at the $n$-th receive antenna can be described as follows [18]:

$$y_n\left(t\right)=\begin{array}{c}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{q=1}^Q{\beta }_l{x}_q\left(\frac{t-{\tau }_l}{{\sigma }_l}\right)B_q\left({\phi }_l\right)A_n\left({\theta }_l\right)}}+w_n\left(t\right)\end{array},$$

(1)

where $x_q(t)$ is the transmitted signal by the $q$-th transmit antenna, expressed as

$$x_q(t)=A_q\exp \left[j2\pi \left(f_{q0}t+{\mu }_{q0}{t}^2/2\right)\right],$$

(2)

where $A_q$ is the amplitude of the transmitted signal, $f_{q0}$ and ${\mu }_{q0}$ are the initial frequency and the chirp rate, respectively, ${\beta }_l$ is the radar cross-section (RCS), ${\sigma }_l$ is the Doppler stretch, ${\tau }_l$ is the time delay, $B_q\left({\phi }_l\right)=\exp \left(j\mathsf{\pi }\left(q-1\right)\sin {\phi }_l\right)$ is the q-th element of the transmitter steering vector, $A_n\left({\theta }_l\right)=\exp \left(j\mathsf{\pi }\left(n-1\right)\sin {\theta }_l\right)$ is the $n$-th element of the receiver steering vector, and $w_n\left(t\right)$ is a sequence of independent and identically distributed (i.i.d.) isotropic complex $S\alpha S$ random variables.

2.1.2. Bandpass Matched Filter

The fractional Fourier transform (FRFT) is a generalization of the FT, and can be interpreted as a rotation of the signal to any angles in the time–frequency plane [28]. The continuous FRFT of a signal $f(t)$ with angle $\rho$ is defined as

$$F\left(\rho ,m\right)=F^{\rho }\left[f\left(t\right)\right]\left(m\right)={\displaystyle {\int }_{-\infty }^{+\infty }f\left(t\right)}{K}_{\rho }\left(t,m\right)dt,$$

(3)

where $F^{\rho }$ denotes the FRFT operator and $K_{\rho }(t,m)$ is the kernel function of the fractional Fourier transform. $K_{\rho }(t,m)$ can be expressed as

$$K_{\rho }(t,m)=\{\begin{array}{l}\sqrt{(1-j\cot \rho )}\exp \left(j\pi (t^2\cot \rho -2mt\csc \rho +m^2\cot \rho )\right),\rho \ne n\pi \\ \delta (t-m),\rho =2n\pi \\ \delta (t+m),\rho =(2n+1)\pi \end{array},$$

(4)

where $\rho$ and $m$ are the rotation angle and the frequency in the FRFT domain, respectively.

According to Equations (2) and (3), the FRFT of transmitted signal $x_q(t)$ with an angle $\rho$ is defined as

$$\begin{array}{l}{X}_q\left(\rho ,m\right)=A\sqrt{\left(1-j\cot \rho \right)}\exp \left(j\pi m^2\cot \rho \right)\\ \cdot {\displaystyle {\int }_0^T\exp \left(j2\pi t\left(f_{q0}-m\csc \rho \right)\right)\exp \left(j\pi t^2\left(\cot \rho +{\mu }_{q0}\right)\right)}dt.\end{array}.$$

(5)

$X_q\left(\rho ,m\right)$ produces the peak value when ${\mu }_{q0}=-\cot {\rho }_{q0}$ and $f_{q0}=m_{q0}\csc {\rho }_{q0}$.

The FRFT $X_q\left({\rho }_{q0},m\right)$ of signal $x_q(t)$ with optimal angle ${\rho }_{q0}$ has an obvious peak value, and the energy of $X_q\left({\rho }_{q0},m\right)$ concentrates in a narrow band with the central frequency of $m_{q0}$.

A $q$-th bandpass matched filter with suitable bandwidth and central frequency $m_{q0}$ is designed. Let $R_{qn}\left(\alpha ,m\right)$ denote the output of the $q$-th matched filter at the $n$-th receive antenna. $y_{qn}\left(t\right)$ represents the output of the matched filter in the FRFT domain, which corresponds to the $q$-th transmit antenna signal in the n-th receive antenna. Through the inverse fractional Fourier transformation of $R_{qn}\left(\alpha ,m\right)$, $y_{qn}\left(t\right)$ can be expressed as

$$y_{qn}\left(t\right)={\displaystyle \sum _{l=1}^L{\beta }_l{x}_q\left(\frac{t-{\tau }_l}{{\sigma }_l}\right)A_n\left({\theta }_l\right)B_q\left({\phi }_l\right)+w_{qn}\left(t\right)},$$

(6)

where $w_{qn}\left(t\right)$ denotes the $S\alpha S$ noise in the received signal $y_{qn}\left(t\right)$.

2.2. $S\alpha S$ Distribution Noise Model

Symmetric $\alpha$ stable ($S\alpha S$) processes are the subject of increased attention as a suitable framework for efficient signal processing in impulsive environments [19,20,29,30]. A “zero-centered” symmetric $\alpha$ stable distribution is commonly described through its characteristic function as

$$\psi \left(\omega \right)=\exp \left(-\gamma {\left|\omega \right|}^{\alpha }\right),$$

(7)

where parameter $\alpha$ is usually called the characteristic exponent. It can be proven that, in order to define a characteristic function for Equation (7), the values of $\alpha$ must be restricted to the interval $\left(0,2\right]$. When $\alpha <2$, the distribution is algebraic-tailed with tail constant $\alpha$, implying infinite variance. When $\alpha =2$, the SαS distribution reduces to the Gaussian distribution implying lighter-than-algebraic tails. The parameter $\gamma$, usually called the dispersion, is a positive constant related to the scale of the distribution. For a fixed $\alpha$, larger values of $\gamma$ correspond to larger strengths of the process. It is easy to see that ${\gamma }^{1/\alpha }$ is, in fact, a scale parameter of the distribution.

Let $X$ be a logarithmic-order random variable. We define the geometric power of $X$ as follows [23]:

$$S_0=S_0\left(X\right)=\exp \left(E\left[\log \left|X\right|\right]\right).$$

(8)

The geometric power of a symmetric $\alpha$ stable random variable as defined in Equation (7) is given by

$$S_0={\left(C_g\gamma \right)}^{1/\alpha }/C_g,$$

(9)

where $C_g\approx 1.78$ is the exponential of the Euler constant.

Since the $\alpha$ stable distribution with $\alpha <2$ determines an infinite variance, we describe the signal-to-noise condition of $S\alpha S$ using the generalized signal-noise-ratio (GSNR) [19,29,30], which is defined as

$$\mathrm{GSNR}=10\lg \left({\sigma }_x^2/\gamma \right),$$

(10)

where ${\sigma }_x^2$ is the variance of the underlying signal.

3. Sigmoid Wideband Ambiguity Function and Sigmoid Correlation

3.1. Wideband Ambiguity Function

For a joint estimation of TD and DS, Swick [31,32] suggested the application of a wideband ambiguity function (WBAF), defined by

$$W_{s_{r}s}\left(\tau ,\sigma \right)=\frac{1}{\sqrt{\sigma }}{\displaystyle {\int }_{-\infty }^{\infty }{s}_r\left(t\right)}{s}^{*}\left(\frac{t-\tau }{\sigma }\right)\mathrm{d}t,$$

(11)

where $s_r\left(t\right)=s\left(\frac{t-{\tau }_0}{{\sigma }_0}\right)$, and ${\tau }_0$ and ${\sigma }_0$ are the time delay and Doppler stretch, respectively.

From the Schwarz inequality, we can see that

$$\begin{array}{ll}{\left|W_{s_{r}s}\left(\tau ,\sigma \right)\right|}^2& =\frac{1}{\sigma }{\left|{\displaystyle {\int }_{-\infty }^{\infty }{s}_r\left(t\right)}{s}^{*}\left(\frac{t-\tau }{\sigma }\right)\mathrm{d}t\right|}^2\\ & =\frac{1}{\sigma }{\left|{\displaystyle {\int }_{-\infty }^{\infty }{s}_r\left(\frac{t-{\tau }_0}{{\sigma }_0}\right)}{s}^{*}\left(\frac{t-\tau }{\sigma }\right)\mathrm{d}t\right|}^2\\ & \le \left[{{\displaystyle {\int }_{-\infty }^{\infty }\left|s_r\left(\frac{t-{\tau }_0}{{\sigma }_0}\right)\right|}}^2\mathrm{d}t\right]\left[\frac{1}{\sigma }{{\displaystyle {\int }_{-\infty }^{\infty }\left|s\left(\frac{t-\tau }{\sigma }\right)\right|}}^2\mathrm{d}t\right]\end{array}.$$

(12)

Equation (12) satisfies equality if and only if $s\left(\frac{t-\tau }{\sigma }\right)\propto s_r\left(\frac{t-{\tau }_0}{{\sigma }_0}\right)$ is equivalent to $\tau ={\tau }_0$ and $\sigma ={\sigma }_0$ when $\sigma >0$ $\sigma >0$ and ${\sigma }_0>0$. Since ${\left|W_{s_{r}s}\left(\tau ,\sigma \right)\right|}^2$ attains its maximum only at $\tau ={\tau }_0$ and $\sigma ={\sigma }_0$, the estimation of TD and DS becomes a problem in locating the maximum point of ${\left|W_{s_{r}s}\left(\tau ,\sigma \right)\right|}^2$, and the corresponding coordinate is the true TD and DS point $\left({\tau }_0,{\sigma }_0\right)$ and can be depicted as follows [12]:

$$\left(\tau ={\tau }_0,\sigma ={\sigma }_0\right)=\arg \max \left|W_{s_{r}s}\left(\tau ,\sigma \right)\right|.$$

(13)

According to Equation (11), we obtain the wideband ambiguity function $W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)$ of $x_q\left(t\right)$ and $x_{r,qnl}\left(t\right)$ as

$$W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)=\frac{1}{\sqrt{\sigma }}{\displaystyle {\int }_{-\infty }^{\infty }{x}_{r,qnl}\left(t\right)}{x}_q^{*}\left(\frac{t-\tau }{\sigma }\right)\mathrm{d}t,$$

(14)

where $x_{r,qnl}\left(t\right)={\beta }_l{x}_q\left[\left(t-{\tau }_l\right)/{\sigma }_l\right]B_q\left({\phi }_l\right)A_n\left({\theta }_l\right)$, $l$ denotes the $l$-th target, and $l=1,\dots ,L$.

From the Schwarz inequality, we can also see that

$$\begin{array}{ll}{\left|W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)\right|}^2& =\frac{1}{\sigma }{\left|{\displaystyle {\int }_{-\infty }^{\infty }{x}_{r,qnl}\left(t\right)}{x}_q^{*}\left(\frac{t-\tau }{\sigma }\right)\mathrm{d}t\right|}^2\\ & \le {\left|{\beta }_l\right|}^2\left[{{\displaystyle {\int }_{-\infty }^{\infty }\left|x_q\left(\frac{t-{\tau }_l}{{\sigma }_l}\right)\right|}}^2\mathrm{d}t\right]\left[\frac{1}{\sigma }{{\displaystyle {\int }_{-\infty }^{\infty }\left|x_q\left(\frac{t-\tau }{\sigma }\right)\right|}}^2\mathrm{d}t\right]\end{array}.$$

(15)

Equation (15) holds with equality if and only if $x_q\left(\frac{t-\tau }{\sigma }\right)\propto x_q\left(\frac{t-{\tau }_l}{{\sigma }_l}\right)$ is equivalent to $\tau ={\tau }_l$ and $\sigma ={\sigma }_l$ when $\sigma >0$ and ${\sigma }_l>0$. Since ${\left|W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)\right|}^2$ attains its maximum only at $\tau ={\tau }_l$ and $\sigma ={\sigma }_l$, the estimation of TD and DS becomes a problem in locating the maximum point of ${\left|W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)\right|}^2$, and the corresponding coordinate is the true TD and DS point $\left({\tau }_l,{\sigma }_l\right)$.

Then, it follows directly from Equation (16) that the TD and DS are determined by

$$\begin{array}{l}\left({\tau }_l,{\sigma }_l\right)=\arg \underset{\tau ,\sigma }{\max }\left[\left|W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)\right|\right]\\ W_{x_{r,qnl}{x}_q}\left({\tau }_l,{\sigma }_l\right)=\max \left[\left|W_{x_{r,qnl}{x}_q}\left(\tau ,\sigma \right)\right|\right]\end{array}\}.$$

(16)

When the alpha stable distribution noise is added, this peak location algorithm may fail. The reason is that the alpha stable distribution does not have a finite $\alpha$-order moment and other higher-than-$\alpha$-order moments, and the wideband ambiguity function is based on a second-order moment. Accordingly, the WBAF algorithm becomes unbounded when the received noise contains a $\alpha$-stable distribution noise. Therefore, we present a nonlinear transform, the Sigmoid transform, to suppress the $\alpha$ stable distribution noise interference.

3.2. Sigmoid Transform

Sigmoid is a commonly used nonlinear transform [33,34,35]. Its definition is shown in Equation (17).

$$\mathrm{Sigmoid}\left[x\left(t\right)\right]=\frac{2}{1+\exp \left[-x\left(t\right)\right]}-1.$$

(17)

Two properties can be obtained from the analysis of the Sigmoid function transform.

Property 1.

If$x\left(t\right)$is an$S\alpha S$process with$\beta =0$and$a=0$, then$\mathrm{Sigmoid}\left[x\left(t\right)\right]$is a symmetric distribution with zero mean in its probability density function, and has the finite second-order moment with zero mean (referred to as a second-order moment process).

Property 2.

Set$X\left(t\right)=\mathrm{Sigmoid}\left[x\left(t\right)\right]$; then,$X\left(t\right)$has the same frequency shift as$x\left(t\right)$.

Since Properties 1 and 2 were proven in References [35] and [29,30], respectively, the relevant proof is skipped herein.

3.3. Sigmoid Wideband Ambiguity Function

A novel ambiguity function based on the Sigmoid transform, the Sigmoid wideband ambiguity function (referred to as Sigmoid-WBAF), $W_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$, is defined as

$$W_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)=\frac{1}{\sqrt{\sigma }}{\displaystyle {\int }_{-\infty }^{+\infty }\mathrm{Sigmoid}\left[s_r\left(t\right)\right]}{\mathrm{Sigmoid}}^{*}\left[s\left(\frac{t-\tau }{\sigma }\right)\right]\mathrm{d}t.$$

(18)

This is the form of Sigmoid-WBAF used throughout this paper. Note that $W_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$ is a 2D representation of the Sigmoid correlation between $s_r\left(t\right)$ and $s\left(\frac{t-\tau }{\sigma }\right)$ for various values of $\tau$ and $\sigma$.

We note that, in Equation (18), the limits of the integral are from $-\infty$ to $\infty$. However, in practice, we only use signals which are essentially time-limited to $[-T/2,T/2]$, meaning that the signal amplitude is negligible outside the essential duration. Thus, the practical implication of employing Equation (13) is that it will be applied only during the essential duration of $[-T/2,T/2]$. Also, in practice and for a time finite signal, $W_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$ can be estimated by Equation (19).

$${\widehat{W}}_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)=\frac{1}{\sqrt{\sigma }}{\displaystyle {\int }_{-T/2}^{T/2}\mathrm{Sigmoid}\left[s_r\left(t\right)\right]}{\mathrm{Sigmoid}}^{*}\left[s\left(\frac{t-\tau }{\sigma }\right)\right]\mathrm{d}t.$$

(19)

Similarly, the joint estimates of TD and DS can be obtained by

$$\begin{array}{r}\left({\tau }_0,{\sigma }_0\right)=\arg \underset{\tau ,\sigma }{\max }\left[\left|{\widehat{W}}_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)\right|\right]\\ {\widehat{W}}_{s_{r}s}^{\mathrm{Sigmoid}}\left({\tau }_0,{\sigma }_0\right)=\max \left[\left|{\widehat{W}}_{s_{r}s}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)\right|\right]\end{array}\}.$$

(20)

3.4. Sigmoid Correlation

For an alpha stable distribution noise environment, the conditional MUSIC algorithm performance degrades seriously and even becomes invalid because the conditional MUSIC algorithm is based on a second-order moment. To suppress the alpha stable distribution noise interference, a novel correlation function is proposed in this section.

A novel correlation based on the Sigmoid transform $R_x^{\mathrm{Sigmoid}}\left(\tau \right)$, referred to as the Sigmoid correlation (SC), is defined in Equation (21).

$$\begin{array}{ll}{R}_x^{\mathrm{Sigmoid}}\left(\tau \right)& ={\displaystyle {\int }_{-\mathrm{T}/2}^{\mathrm{T}/2}\mathrm{Sigmoid}\left[x\left(t\right)\right]{\mathrm{Sigmoid}}^{*}\left[x\left(t-\tau \right)\right]}\mathrm{d}t\\ & ={\langle \mathrm{Sigmoid}\left[x\left(t\right)\right]{\mathrm{Sigmoid}}^{*}\left[x\left(t-\tau \right)\right]\rangle }_t\end{array},$$

(21)

where ${\langle \cdot \rangle }_t$ represents a time average.

According to the properties of the Sigmoid transform, we can deduce that the properties of correlation do not change if the Sigmoid transform is applied on both signals in advance. Therefore, the Sigmoid correlation function not only has the properties of a conventional autocorrelation function, but also has a suppression ability to impulsive noise. Therefore, the Sigmoid correlation function can be used to estimate DOD and DOA in an impulsive noise environment.

3.5. Sigmoid-MUSIC Algorithm

Let

$$y\left(t\right)=As\left(t\right)+n\left(t\right),$$

(22)

where $y\left(t\right)={\left[y_1\left(t\right)y_2\left(t\right)\cdots y_M\left(t\right)\right]}^{\mathrm{T}}$ is the $M\times 1$ vector of signals received by the array sensors, $s\left(t\right)={\left[s_1\left(t\right)s_2\left(t\right)\cdots s_P\left(t\right)\right]}^{\mathrm{T}}$ is the $P\times 1$ vector of the signals, $n\left(t\right)={\left[n_1\left(t\right)n_2\left(t\right)\cdots n_P\left(t\right)\right]}^{\mathrm{T}}$ is the $M\times 1$ vector of impulsive noise, and $A=\left[a\left({\theta }_1\right)\cdots a\left({\theta }_i\right)\cdots a\left({\theta }_P\right)\right]$ is the matrix of array steering vectors, in which $a\left({\theta }_i\right)=\left[1\cdots e^{j2\mathsf{\pi }\left(m-1\right)\mathrm{d}\sin {\theta }_i/\lambda }\cdots e^{j2\mathsf{\pi }\left(M-1\right)\mathrm{d}\sin {\theta }_i/\lambda }\right]$.

The Sigmoid-MUSIC algorithm can be implemented using the following procedure:

Step 1.

Compute Sigmoid correlation matrices $R_y^{\mathrm{Sigmoid}}\left(\tau \right)$ of the matrix $y\left(t\right)$, according to Equation (21).

Step 2.

Execute singular value decomposition (SVD) on $R_y^{\mathrm{Sigmoid}}\left(\tau \right)$, where the column vector $U_N$ describes the eigenvectors spanning the noise subspace.

Step 3.

Compute the corresponding Sigmoid-MUSIC spectrum $P_{\mathrm{Sigmoid}-\mathrm{MUSIC}}\left(\theta \right)$ as

$$P_{\mathrm{Sigmoid}-\mathrm{MUSIC}}\left(\theta \right)=\frac{1}{A^{\mathrm{H}}{U}_N{U}_N^{\mathrm{H}}A}.$$

(23)

Step 4.

The estimator of $\theta$ can be obtained by searching for peaks of the Sigmoid-MUSIC spectrum $P_{\mathrm{Sigmoid}-\mathrm{MUSIC}}\left(\theta \right)$.

4. Joint Estimation Parameter Based on Sigmoid-WBAF and Sigmoid-MUSIC

In this section, a study of parameter estimation is made by considering signal $y_{qnl}\left(t\right)$ as an example. Signal $y_{qnl}\left(t\right)$ denotes the received signals $y_{qn}\left(t\right)$ corresponding to the $l$-th target. $y_{qnl}\left(t\right)$ can be expressed as

$$y_{qnl}\left(t\right)={\beta }_{l}x\left(\frac{t-{\tau }_l}{{\sigma }_l}\right)A_n\left({\theta }_l\right)B_q\left({\phi }_l\right)+w_{qnl}\left(t\right),$$

(24)

where $w_{qnl}\left(t\right)$ denotes the $S\alpha S$ noise in the received signal $y_{qnl}\left(t\right)$.

4.1. Estimation of TD and DS Based on Sigmoid-WBAF

According to Equations (2), (19), and (24), we can obtain the Sigmoid-WBAF of $y_{qnl}\left(t\right)$ and $x_q\left(t\right)$ as follows:

$$\begin{array}{ll}{\widehat{W}}_{y_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)& =\frac{1}{\sqrt{\sigma }}{\displaystyle {\int }_{-T/2}^{T/2}\mathrm{Sigmoid}\left[y_{qnl}\left(t\right)\right]}{\mathrm{Sigmoid}}^{*}\left[x_q\left(\frac{t-\tau }{\sigma }\right)\right]\mathrm{d}t\\ & ={\widehat{W}}_{x_{r,qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)+{\widehat{W}}_{w_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)\end{array},$$

(25)

where ${\widehat{W}}_{w_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$ is the Sigmoid wideband ambiguity function of noise $w_{qnl}\left(t\right)$ and transmitted signal $x_q\left(t\right)$, and it is treated as a random interference.

Similarly, the estimators of TD and DS can be obtained from

$$\begin{array}{r}\left(\tau ={\widehat{\tau }}_l,\sigma ={\widehat{\sigma }}_l\right)=\arg \underset{\tau ,\sigma }{\max }\left[\left|{\widehat{W}}_{y_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)\right|\right]\\ {\widehat{W}}_{y_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left({\widehat{\tau }}_l,{\widehat{\sigma }}_l\right)=\max \left[\left|{\widehat{W}}_{y_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)\right|\right]\end{array}\}.$$

(26)

Accordingly, the estimation of the Doppler stretch and time delay in a wideband bistatic MIMO radar was achieved via the proposed Sigmoid wideband ambiguity function. The steps involved in this process were as follows:

Step 1.

Present the extracted signal $y_{qn}\left(t\right)$.

Step 2.

Compute the Sigmoid-WBAF function ${\widehat{W}}_{y_{qn}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$ from Equation (25).

Step 3.

Search for the peaks of ${\widehat{W}}_{y_{qn}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$ and obtain the locations of these peaks $\left({\widehat{\tau }}_l,{\widehat{\sigma }}_l\right)$, for $l=1,\dots ,L$.

Step 4.

Estimate the DS and TD according to Equation (26).

4.2. Estimation of DOD and DOA Based on Sigmoid-MUSIC

In this section, DOD and DOA are estimated by employing the proposed Sigmoid-MUSIC.

The vector form of the array output can be shown as

$$Y\left(t\right)=A\otimes BS+N\left(t\right),$$

(27)

where $A=\left[a_1,\dots ,a_L\right]$ and $B=\left[b_1,\dots ,b_L\right]$ are the receive and transmit steering matrices, where the receiver and transmitter steering vectors are given by $a_l={\left[1,e^{j\pi \sin {\theta }_l},\dots ,e^{j\pi \left(N-1\right)\sin {\theta }_l}\right]}^{\mathrm{T}}$ and $b_l={\left[1,e^{j\pi \sin {\phi }_l},\dots ,e^{j\pi \left(Q-1\right)\sin {\phi }_l}\right]}^{\mathrm{T}}$. $S=\left[s_1,\dots ,s_L\right]$, where $s_l={\left[{\beta }_l{x}_1\left(\frac{t-{\tau }_l}{{\sigma }_l}\right),\dots ,{\beta }_l{x}_Q\left(\frac{t-{\tau }_l}{{\sigma }_l}\right)\right]}^{\mathrm{T}}$ is a known vector since estimated values of ${\tau }_l$ and ${\sigma }_l$ are used, $\otimes$ denotes the Kronecker Matrix product, and $N\left(t\right)$ is a sequence of i.i.d isotropic complex $S\alpha S$ random variables.

According to Equation (27), two receive subarrays $Y_1$ and $Y_2$ can be constructed as follows:

$$Y_1=AS+N_1\hspace{1em}q=1;$$

(28)

$$Y_2=BS+N_2\hspace{1em}n=1.$$

(29)

According to Equation (21), we can obtain the Sigmoid correlation function $R_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)$ of the received signal vector $Y_1\left(t\right)$,

$$\begin{array}{ll}{R}_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)& ={\langle \mathrm{Sigmoid}\left[Y_1\left(t\right)\right]{\mathrm{Sigmoid}}^{*}\left[Y_1\left(t-\tau \right)\right]\rangle }_t\\ & =A{R}_S^{\mathrm{Sigmoid}}\left(\tau \right)A^{\mathrm{H}}+A{R}_{S{N}^{*}}^{\mathrm{Sigmoid}}\left(\tau \right)+R_{S^{*}N}^{\mathrm{Sigmoid}}\left(\tau \right)A^{\mathrm{H}}+R_{N{N}^{*}}^{\mathrm{Sigmoid}}\left(\tau \right)\end{array},$$

(30)

where $R_S^{\mathrm{Sigmoid}}\left(\tau \right)={\langle \mathrm{Sigmoid}\left[s\left(t\right)\right]{\mathrm{Sigmoid}}^{*}\left[s\left(t-\tau \right)\right]\rangle }_t$ denotes the Sigmoid autocorrelation of the matrix $S$. Note that the correlation matrix used in this algorithm is replaced by the Sigmoid correlation matrix.

Signal $s\left(t\right)$ is independent of noise $n\left(t\right)$; thus, Equation (30) can be rewritten as

$$R_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)=A{R}_S^{\mathrm{Sigmoid}}\left(\tau \right)A^{\mathrm{H}}.$$

(31)

Equation (31) can then be rewritten as

$$R_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)=\left[\begin{array}{cc}{U}_S& U_N\end{array}\right]\left[\begin{array}{cc}{\sum }_S& 0\\ 0& {\sum }_N\end{array}\right]{\left[\begin{array}{cc}{V}_S& V_N\end{array}\right]}^H=U_S{\sum }_S{V}_S^H,$$

(32)

where $\left[\begin{array}{cc}{U}_S& U_N\end{array}\right]$ and $\left[\begin{array}{cc}{V}_S& V_N\end{array}\right]$ are unitary, and the diagonal elements of the diagonal matrix ${\sum }_S$ are positive. The column vectors of $U_S$ and $U_N$ are the eigenvectors spanning the signal subspace and noise subspace of $R_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)$, respectively, with the associated eigenvalues as the diagonals of ${\sum }_S$ and ${\sum }_N$. Since signal $S$ is independent of the noise $N$,

$$R_S^{\mathrm{Sigmoid}}\left(\tau \right)U_N=0,$$

(33)

and the signal and noise subspaces are orthogonal, i.e.,

$$A^{\mathrm{H}}{U}_N=0.$$

(34)

Therefore, the spatial spectrum of Sigmoid-MUSIC can be obtained based on the classical MUSIC algorithm, which can be expressed as

$$P_{Y_1}\left(\theta \right)=\frac{1}{A^{\mathrm{H}}{U}_N{U}_N^{\mathrm{H}}A}.$$

(35)

Searching for the spectral peak of $P_{Y_1}\left(\theta \right)$, we can get the DOA estimator ${\theta }_l$.

Similarly, we can obtain the Sigmoid correlation function $R_{Y_2}^{\mathrm{Sigmoid}}\left(\tau \right)$ of the received signal vector $Y_2\left(t\right)$ according to Equations (29)–(32) as

$$R_{Y_2}^{\mathrm{Sigmoid}}\left(\tau \right)=\left[\begin{array}{cc}{U}_{S_2}& U_{N_2}\end{array}\right]\left[\begin{array}{cc}{\sum }_{S_2}& 0\\ 0& {\sum }_{N_2}\end{array}\right]{\left[\begin{array}{cc}{V}_{S_2}& V_{N_2}\end{array}\right]}^H=U_{S_2}{\sum }_{S_2}{V}_{S_2}^H,$$

(36)

where $\left[\begin{array}{cc}{U}_{S_2}& U_{N_2}\end{array}\right]$ and $\left[\begin{array}{cc}{V}_{S_2}& V_{N_2}\end{array}\right]$ are unitary, and the diagonal elements of the diagonal matrix ${\sum }_{S_2}$ are positive. The column vectors of $U_{S_2}$ and $U_{N_2}$ are the eigenvectors spanning the signal subspace and noise subspace of $R_{Y_2}^{\mathrm{Sigmoid}}\left(\tau \right)$, respectively, with the associated eigenvalues as the diagonals of ${\sum }_{S_2}$ and ${\sum }_{N_2}$.

Therefore, the Sigmoid-MUSIC algorithm was employed on matrix $Y_2$; we can also obtain the spatial spectrum of Sigmoid-MUSIC for the subarray $Y_2$ as

$$P_{Y_2}\left(\phi \right)=\frac{1}{B^{\mathrm{H}}{U}_{N_2}{U}_{N_2}^{\mathrm{H}}B}.$$

(37)

Searching for the spectral peak of $P_{Y_2}\left(\phi \right)$, we can get the DOD estimator ${\phi }_l$.

Accordingly, the estimation of DOD and DOA in a wideband bistatic MIMO radar was achieved via the proposed Sigmoid-MUSIC algorithm. The steps involved in this process were as follows:

Step 1.

Construct two matrices $Y_1$ and $Y_2$.

Step 2.

Substitute the time average with the statistic average, two Sigmoid correlation matrices $R_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)$ and $R_{Y_2}^{\mathrm{Sigmoid}}\left(\tau \right)$ are then constructed according to Equation (21).

Step 3.

Apply the singular value decomposition (SVD) to $R_{Y_1}^{\mathrm{Sigmoid}}\left(\tau \right)$ and $R_{Y_2}^{\mathrm{Sigmoid}}\left(\tau \right)$, where the column vectors $U_N$ and $U_{N_2}$ are formed from the eigenvectors spanning the noise subspace.

Step 4.

Compute the corresponding Sigmoid-MUSIC spectra $P_{Y_1}\left(\theta \right)$ and $P_{Y_2}\left(\phi \right)$ from Equations (36) and (37).

Step 5.

The DOA and DOD estimates can be obtained by identifying the peaks of the spatial spectra $P_{Y_1}\left(\theta \right)$ and $P_{Y_2}\left(\phi \right)$.

5. Analysis of Sigmoid-WBAF and Sigmoid-MUSIC

5.1. Boundness of Sigmoid-WBAF

According to the properties of the Sigmoid transform, the SαS process with $a=0$ can be transformed to a second-order moment process by the Sigmoid transform. Therefore, ${\widehat{W}}_{y_{qnl}{x}_q}^{\mathrm{Sigmoid}}\left(\tau ,\sigma \right)$ is bounded for the $S\alpha S$ process because it is only involved with $\mathrm{Sigmoid}\left[x\left(t\right)\right]$ under the $S\alpha S$ noise. Furthermore, the transformation does not change the estimation results of the time delay and Doppler frequency shift. Therefore, the Sigmoid-WBAF method can be used to estimate the parameters of wideband echoes $y\left(t\right)$ under $S\alpha S$ stable distribution noise.

5.2. Feasibility Analysis of Sigmoid-WBAF

Simulation results and theoretical analyses illustrated that the Sigmoid transform does not change the modulation characteristics of the signal [29]. Therefore, the characteristics of WBAF do not change if the Sigmoid transform is applied to the signal in advance. This simulation result is illustrated in Figure 2 to verify this property. Figure 2 shows the spectra of WBAF and Sigmoid-WBAF.

From Figure 2, the WBAF of the LFM signal and Sigmoid-WBAF of the LFM signal have the same peak location in the WBAF domain. Therefore, the parameters of the Doppler stretch and time delay can be estimated by searching for the peak of Sigmoid-WBAF.

5.3. The Cramer–Rao Bound

In this section, we derive a novel explicit expression for the exact Cramer–Rao Bound (CRB) on the accuracy of estimating the signal model parameters [29].

The received signal can be expressed as

$$y\left(t\right)=\Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta +N\left(t\right),$$

(38)

where

$$\Re \left(\sigma ,\tau ,\phi ,\theta \right)=\left[\Re \left({\sigma }_1,{\tau }_1,{\phi }_1,{\theta }_1\right),\dots ,\Re \left({\sigma }_L,{\tau }_L,{\phi }_L,{\theta }_L\right)\right],$$

(39)

$$\Re \left({\sigma }_l,{\tau }_l,{\phi }_l,{\theta }_l\right)=\left[x\left({\sigma }_l,{\tau }_l,t\right)\circ B\left({\phi }_l\right)\right]\odot A\left({\theta }_l\right),$$

(40)

and $\beta \triangleq {\left[{\beta }_1,\dots ,{\beta }_L\right]}^T$, $x\left(\sigma ,\tau ,t\right)=\left[x\left(\left(t-{\tau }_1\right)/{\sigma }_1\right),x\left(\left(t-{\tau }_2\right)/{\sigma }_2\right),\dots ,x\left(\left(t-{\tau }_L\right)/{\sigma }_L\right)\right]$.

The four parameters to be estimated are the time delay $\tau \triangleq \left[{\tau }_1,{\tau }_2,\dots ,{\tau }_L\right]$, the Doppler stretch $\sigma \triangleq \left[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_L\right]$, the DOD $\phi =\left[{\phi }_1,{\phi }_2,\dots ,{\phi }_L\right]$, and the DOA $\theta =\left[{\theta }_1,{\theta }_2,\dots ,{\theta }_L\right]$, which form the parameter vector $\xi$ as $\xi ={\left[a,\tau ,\phi ,\theta \right]}^T$.

The element at $i,j$ in the Fisher information matrix (FIM) for estimating the vector $\xi$ can be depicted as follows [21,32,36,37,38]:

$${\mathsf{\Gamma }}_{ij}\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\left(\frac{\partial \Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta }{\partial {\xi }_i}\right)}^H{Q}_n^{-1}\left(\frac{\partial \Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta }{\partial {\xi }_j}\right)\right\}},$$

(41)

where $N$ denotes the number of snapshots, and $Q_n=S_0{I}_N$.

For Equation (41), we may calculate the partial derivations.

$$\begin{array}{ll}\frac{\partial \left\{{\beta }^{\mathrm{H}}{\Re }^{\mathrm{H}}\left(\sigma ,\tau ,\phi ,\theta \right)\right\}}{\partial {\sigma }_l}& ={\beta }^{\mathrm{H}}{A}^{\mathrm{H}}\left(\theta \right)e_l{e}_l^{\mathrm{T}}\odot {\left[{x^{\prime }}_{\left(\sigma \right)}\left(\sigma ,\tau ,t\right)\circ B\left(\phi \right)\right]}^{\mathrm{H}}\\ & ={\beta }^{\mathrm{H}}{e}_l{e}_l^{\mathrm{T}}{A}^{\mathrm{H}}\left(\theta \right)\odot {\left[{x^{\prime }}_{\left(\sigma \right)}\left(a,\tau ,t\right)\circ B\left(\phi \right)\right]}^{\mathrm{H}}\\ & ={\beta }^{\mathrm{H}}{e}_l{e}_l^{\mathrm{T}}{{\Re }^{\prime }}_a^{\mathrm{H}}\left(\sigma ,\tau ,\phi ,\theta \right)\end{array},$$

(42)

where $e_l$ is the $l$-th column of the identity matrix (i.e., vector containing one in the $l$-th position and zeroes elsewhere). The following notations are introduced:

$${{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)=\frac{\partial \Re \left(\sigma ,\tau ,\phi ,\theta \right)}{\partial \sigma }=\left[{x^{\prime }}_{\left(\sigma \right)}\left(\sigma ,\tau ,t\right)\circ B\left(\phi \right)\right]\odot A\left(\theta \right),$$

(43)

where $\circ$ is the Schur–Hadamard matrix product and $\odot$ is Khatri–Rao matrix product.

$${x^{\prime }}_{\left(\sigma \right)}\left(\sigma ,\tau ,t\right)=\left[\frac{\partial x\left(\left(t-{\tau }_1\right)/{\sigma }_1\right)}{\partial {\sigma }_1},\frac{\partial x\left(\left(t-{\tau }_2\right)/{\sigma }_2\right)}{\partial {\sigma }_2},\dots ,\frac{\partial x\left(\left(t-{\tau }_L\right)/{\sigma }_L\right)}{\partial {\sigma }_L}\right].$$

(44)

From Equations (42)–(44), we obtain that

$$\begin{array}{ll}\frac{\partial {\left\{\Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta \right\}}^{\mathrm{H}}}{\partial \sigma }& =\frac{\partial \left\{{\beta }^{\mathrm{H}}{\Re }^{\mathrm{H}}\left(\sigma ,\tau ,\phi ,\theta \right)\right\}}{\partial \sigma }\\ & ={\left[{{\Re }^{\prime }}_a\left(\sigma ,\tau ,\phi ,\theta \right)e_1{e}_1^{\mathrm{T}}\beta ,\dots ,{{\Re }^{\prime }}_a\left(\sigma ,\tau ,\phi ,\theta \right)e_L{e}_L^{\mathrm{T}}\beta \right]}^{\mathrm{H}}\\ & ={\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\alpha }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}\end{array},$$

(45)

where $\Delta =diag\left\{{\beta }_1,{\beta }_2,\dots ,{\beta }_L\right\}$.

$$\begin{array}{ll}{\mathsf{\Gamma }}_{\sigma \sigma }\left(\xi \right)& =2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\left(\frac{\partial \Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta }{\partial \sigma }\right)}^{\mathrm{H}}{Q}_n^{-1}\left(\frac{\partial \Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta }{\partial \sigma }\right)\right\}}\\ & =2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{\frac{\partial \left\{{\beta }^{\mathrm{H}}{\Re }^{\mathrm{H}}\left(\sigma ,\tau ,\phi ,\theta \right)\right\}}{\partial \sigma }{Q}_n^{-1}\frac{\partial \left\{\Re \left(\sigma ,\tau ,\phi ,\theta \right)\beta \right\}}{\partial \sigma }\right\}}\\ & =2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}}\end{array}.$$

(46)

Using Equation (46), the following explicit expressions for the blocks of the FIM were derived for the proposed signal model:

$${\mathsf{\Gamma }}_{\tau \tau }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(47)

$${\mathsf{\Gamma }}_{\phi \phi }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\phi }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(48)

$${\mathsf{\Gamma }}_{\theta \theta }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\theta }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\theta }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(49)

$${\mathsf{\Gamma }}_{\sigma \tau }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(50)

$${\mathsf{\Gamma }}_{a\phi }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\phi }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(51)

$${\mathsf{\Gamma }}_{a\theta }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\theta }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(52)

$${\mathsf{\Gamma }}_{\tau \phi }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\phi }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(53)

$${\mathsf{\Gamma }}_{\tau \theta }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\theta }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(54)

$${\mathsf{\Gamma }}_{\phi \theta }\left(\xi \right)=2\mathrm{Re}{\displaystyle \sum _{t=1}^N\left\{{\Delta }^{\mathrm{H}}{\left({{\Re }^{\prime }}_{\phi }\left(\sigma ,\tau ,\phi ,\theta \right)\right)}^{\mathrm{H}}{Q}_n^{-1}\left({{\Re }^{\prime }}_{\theta }\left(\sigma ,\tau ,\phi ,\theta \right)\right)\Delta \right\}};$$

(55)

where

$$\Delta =diag\left\{{\beta }_1,{\beta }_2,\dots ,{\beta }_L\right\};$$

(56)

$${{\Re }^{\prime }}_{\sigma }\left(\sigma ,\tau ,\phi ,\theta \right)=\frac{\partial {\Re }_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)}{\partial \sigma }=\left[{x^{\prime }}_{\left(\sigma \right)}\left(\sigma ,\tau ,t\right)\circ B\left(\phi \right)\right]\odot A\left(\theta \right);$$

(57)

$${x^{\prime }}_{\left(\sigma \right)}\left(\sigma ,\tau ,t\right)\triangleq \left[\frac{\partial x\left(\left(t-{\tau }_1\right)/{\sigma }_1\right)}{\partial {\sigma }_1},\frac{\partial x\left(\left(t-{\tau }_2\right)/{\sigma }_2\right)}{\partial {\sigma }_2},\dots ,\frac{\partial x\left(\left(t-{\tau }_L\right)/{\sigma }_L\right)}{\partial {\sigma }_L}\right];$$

(58)

$${{\Re }^{\prime }}_{\tau }\left(\sigma ,\tau ,\phi ,\theta \right)=\left[{x^{\prime }}_{\left(\tau \right)}\left(\sigma ,\tau ,t\right)\circ B\left(\phi \right)\right]\odot A\left(\theta \right);$$

(59)

$${x^{\prime }}_{\left(\tau \right)}\left(\sigma ,\tau ,t\right)\triangleq \left[\frac{\partial x\left(\left(t-{\tau }_1\right)/{\sigma }_1\right)}{\partial {\tau }_1},\frac{\partial x\left(\left(t-{\tau }_2\right)/{\sigma }_2\right)}{\partial {\tau }_2},\dots ,\frac{\partial x\left(\left(t-{\tau }_L\right)/{\sigma }_L\right)}{\partial {\tau }_L}\right];$$

(60)

$${{\Re }^{\prime }}_{\phi }\left(\sigma ,\tau ,\phi ,\theta \right)=\left[x\left(\sigma ,\tau ,t\right)\circ B^{\prime }\left(\phi \right)\right]\odot A\left(\theta \right);$$

(61)

$$B^{\prime }\left(\phi \right)=\left[\frac{\partial B\left({\phi }_1\right)}{\partial {\phi }_1},\frac{\partial B\left({\phi }_2\right)}{\partial {\phi }_2},\dots ,\frac{\partial B\left({\phi }_L\right)}{\partial {\phi }_L}\right];$$

(62)

$${{\Re }^{\prime }}_{\theta }\left(\sigma ,\tau ,\phi ,\theta \right)=\left[x\left(\sigma ,\tau ,t\right)\circ B\left(\phi \right)\right]\odot A^{\prime }\left(\theta \right);$$

(63)

$$A^{\prime }\left(\theta \right)=\left[\frac{\partial A\left({\theta }_1\right)}{\partial {\theta }_1},\frac{\partial A\left({\theta }_2\right)}{\partial {\theta }_2},\dots ,\frac{\partial A\left({\theta }_L\right)}{\partial {\theta }_L}\right].$$

(64)

The expression for the CRB, shown in Equation (65), is obtained by substituting Equations (46)–(64) into Equation (39).

$$\mathrm{CRB}\left(\xi \right)={\mathsf{\Gamma }}^{-1}.$$

(65)

5.4. Complexity Analysis

In this section, we evaluate the computation complexity of the proposed method.

5.4.1. Doppler Stretch and Time Delay

The method based on the WBAF can estimate Doppler stretch and time delay by searching the peak of the WBAF. Denoting the number of snapshots, time delay, and Doppler stretch to be searched as $N_s$, $N_{\tau }$, and $N_{\sigma }$, respectively, the computational complexity of the WBAF is then around $O\left(N_{\tau }{N}_{\sigma }{N}_s\right)$ [39,40]. The computation complexity of the FLOS-WBAF algorithm is $O\left(N_{\tau }{N}_{\sigma }{N}_s\right)$. The parameter estimation based on Sigmoid-WBAF method is a two-step process. The first step carries the Sigmoid transformation and the second step computes the WBAF. Accordingly, the computational complexity of the Sigmoid-WBAF algorithm is $O\left(N_{\tau }{N}_{\sigma }{N}_s+N_s\right)\approx O\left(N_{\tau }{N}_{\sigma }{N}_s\right)$. Through the computational complexity analysis, the Sigmoid-WBAF method not only has the same computation complexity as the WBAF method and FLOS-WBAF method, but also can suppress impulsive noise interference and does not need any a priori knowledge of the noise.

5.4.2. DOD and DOA

The computational complexities of the proposed Sigmoid-MUSIC, MUSIC, and FLOM-MUSIC algorithms were compared. All these methods include the eigen decomposition step which is represented by the term $O\left(N^3\right)$. The computation of $J$ samples of the MUSIC spectrum function requires $O\left(JNL\right)$. Thus, the computational complexity of spectral MUSIC is $O\left(N^3+JNL\right)$, where $N$ denotes the number of the received antennas, and $L$ denotes the number of the targets [41]. The computational complexity of FLOM-MUSIC is also $O\left(N^3+JNL\right)$. The Sigmoid-MUSIC method needs to carry the Sigmoid transformation in advance. Thus, the computational complexity of Sigmoid-MUSIC is $O\left(N^3+JNL+N_s\right)\approx O\left(N^3+JNL\right)$.

Through the computational complexity analysis, we can deduce that the Sigmoid-MUSIC method not only has the same computation complexity as other methods, but also can suppress impulsive noise interference and does not need any a priori knowledge of the noise.

6. Simulation Results

We performed simulation experiments to assess the relative performance of the different methods, including the WBAF [12], FLOS-WBAF [18], and Sigmoid-WBAF to estimate TD and DS; and the conditional MUSIC [21], FLOM-MUSIC [22], FLOM-SC-SSF [14], $l_p$-MUSIC [24], and Sigmoid-MUSIC to estimate DOD and DOA, under $\alpha$ stable noise condition.

The wideband bistatic MIMO radar was composed of $Q=6$ transmit antennas and $N=8$ receive antennas. Suppose that the target is located at the positions $\left({\mathsf{\phi }}_1,{\mathsf{\theta }}_1\right)=\left({10}^{\mathrm{o}},{50}^{\mathrm{o}}\right)$, $\left({\phi }_2,{\theta }_2\right)=\left({40}^{\circ },{30}^{\circ }\right)$; DS and TD are ${\sigma }_1=1.2$, ${\tau }_1=30/f_s$, ${\sigma }_2=0.9$, ${\tau }_2=10/f_s$, respectively, where the sampling frequency $f_s$ is $1\mathrm{KHz}$. The root-mean-square error (RMSE) is defined as

$$\mathrm{RMSE}=\frac{1}{2}\left(\sqrt{\frac{1}{K}{{\displaystyle \sum _{k=1}^K\left[{\widehat{x}}_1\left(k\right)-x_1\right]}}^2}+\sqrt{\frac{1}{K}{{\displaystyle \sum _{k=1}^K\left[{\widehat{x}}_2\left(k\right)-x_2\right]}}^2}\right),$$

(66)

where ${\widehat{x}}_1$ and ${\widehat{x}}_2$ are the estimations of $x_1$ and $x_2$, and $K$ is the total number of successful runs. For each simulation, the numbers of Monte Carlo runs and snapshots were 300 and 1000, respectively.

6.1. Simulation 1: Spectra of WBAF, FLOS-WBAF, and Sigmoid-WBAF for a Single Estimation for Two Targets

In this simulation, we discuss the estimation results obtained from the WBAF, FLOS-WBAF, and Sigmoid-WBAF for two targets in the Gaussian noise and impulsive noise environment. The corresponding results are illustrated in Figure 3, Figure 4 and Figure 5.

In the Gaussian noise, all three methods revealed two clear peaks, as illustrated in Figure 3. However, from Figure 4 and Figure 5, it can clearly be seen that the WBAF algorithm failed for the impulsive noise case. The reason is that the WBAF method does not have the ability to suppress impulsive noise. On the other hand, the FLOS-WBAF algorithm, combining the fractional lower-order statistics theory with the wideband ambiguity function, can effectively suppress the alpha stable noise interference, where clear peaks were obtained for the impulsive noise case with $\mathrm{GSNR}=2 \mathrm{dB}$, $\alpha =1.2$, and $p=1.1$. However, FLOS-WBAF failed to obtain the correct spectrum peaks for impulsive noise with $\alpha =1.2$ and $p=1.4$, mainly due to the fact that the fractional lower-order moment $p$ value was not appropriate, as illustrated in Figure 4. Based on the fractional lower-order statistics theory, the characteristic exponent of the noise must be estimated to ensure $1\le p<\alpha$ or $0<p<\alpha /2$. The methods employing the FLOS theory cannot accurately estimate the parameters if there is no a priori knowledge of the characteristic exponent. Furthermore, the algorithm’s performance degrades seriously and even becomes invalid when the fractional lower-order moment value is outside the specified range. On the contrary, Sigmoid-WBAF had clear peaks regardless of the value of the fractional lower-order moment $p$. As illustrated in Figure 5, FLOS-WBAF failed to obtain the correct peaks for impulsive noise with $\mathrm{GSNR}=-2\mathrm{dB}$, $\alpha =1.2$, and $p=1.1$. However, Sigmoid-WBAF had clear peaks. The reason is that the Sigmoid transform can suppress impulsive noise better than any of the FLOS-type methods, including FLOS-WBAF.

6.2. Simulation 2: Spectrum Performances of the Four Algorithms

In this simulation, the generalized signal-to-noise ratio was set to GSNR = 0 dB and the characteristic exponent $\alpha$ was set to 1.2.

Figure 6 shows the spatial spectra of the MUSIC, FLOM-MUSIC, $l_p$-MUSIC, FLOM-SC-SSF, and Sigmoid-MUSIC algorithms. The performance of the MUSIC algorithm degraded seriously under the α stable distribution noise. Although the spatial spectra of the FLOM-MUSIC, FLOM-SC-SSF, and $l_p$-MUSIC algorithms had two spectrum peaks, the spectrum peak of FLOM-MUSIC showed a deviation from the true location, and the spectrum peaks of FLOM-SC-SSF and $l_p$-MUSIC did not have two spiked peaks. The Sigmoid-MUSIC algorithm based on Sigmoid transform can suppress the impulse noise interference, especially for low GSNR. Therefore, the Sigmoid-MUSIC algorithm had a more accurate spectrum peak and better estimation performance.

6.3. Simulation 3: Generalized Signal-to-Noise Ratio (GSNR)

In this simulation and to evaluate the performances of TD and DS, the characteristic exponent $\alpha$ was set to $\alpha =1.2$, and the fractional lower-order moment $p$ was set to $p=1.0$ and $p=1.6$ for FLOS-WBAF. The parameter p was set to $p=1.4$ for FLOM-MUSIC while estimating the performances of DOD and DOA. The root CRBs and RMSEs of parameter estimation under different GSNRs are shown in Figure 7.

From Figure 7a,b, we can find that the WBAF method had poor estimation performance. The estimation performance of FLOS-WBAF was affected by the fractional-lower order moment $p$. The estimation performance of FLOS-WBAF with $p>\alpha$ had poor performance. The RMSE of Sigmoid-WBAF was obviously lower than that of FLOS-WBAF and WBAF when $\mathrm{GSNR}<5$. When $\mathrm{GSNR}\ge 5$, the performance of Sigmoid-WBAF was similar to the other methods. From Figure 7c,d, the conventional MUSIC method was inferior to the other methods. For low GSNRs, both Sigmoid-MUSIC and $l_p$-MUSIC could get lower RMSEs than other methods; however, Sigmoid-MUSIC yielded more accurate DOA estimation results. Therefore, the estimation performance of the proposed method was superior to other methods. As shown in Figure 7, the root CRB decreased with increasing GSNR. Moreover, the root CRB could be apparently affected by the GSNR for a given characteristic exponent $\alpha$.

6.4. Simulation 4: Characteristic Exponent $\alpha$

In this simulation and to measure the estimation performances of TD and DS, the fractional lower-order moment $p$ was set to $p=\alpha -0.2$ and $p=1.4$ for FLOS-WBAF. The GSNR was set to $\mathrm{GSNR}=5\mathrm{dB}$. While evaluating the performances of DOD and DOA, the fractional lower-order moment $p$ was set to $p=1.4$ for FLOM-MUSIC. The root CRBs and RMSEs on parameter estimation under different characteristic exponents $\alpha$ are shown in Figure 8.

From Figure 8a,b, we can find that WBAF and FLOS-WBAF with $p=1.6$ had poor estimation performance. When $0.5<\alpha <1.5$, the estimation performance of Sigmoid-WBAF was obviously better than that of the FLOS-WBAF method with $p<\alpha$. Although FLOS and Sigmoid transform can both suppress impulsive noise, the suppression capacity of FLOS was not enough because ${\left|x_2\left(t\right)\right|}^p>{\left|x_1\left(t\right)\right|}^p>1$ for any $\left|x_2\left(t\right)\right|>\left|x_1\left(t\right)\right|>1$. For this reason, large outliers of impulsive noise cannot be restrained sufficiently when the impulsiveness is extremely intensive. In contrast, the Sigmoid function can be assumed to be approximately linear. $\left|x\left(t\right)\right|$ can reach significant suppression when $\left|x\left(t\right)\right|$ is far away from zero. Because the signals are often assumed to have a zero-mean value, the Sigmoid function can be used to suppress the outliers [42]. For any $\left|x\left(t\right)\right|>1$, ${\left|x_2\left(t\right)\right|}^p>{\left|x_1\left(t\right)\right|}^p>1>\left|\mathrm{Sigmoid}\left[x\left(t\right)\right]\right|$. Therefore, the Sigmoid function suppresses the outliers much better than FLOS, and the estimation performance of the Sigmoid-WBAF algorithm was better than that of the FLOS-WBAF algorithm. From Figure 8c,d, we can find that the conventional MUSIC method was inferior to other methods. For highly impulsive noise, Sigmoid-MUSIC could get more accurate DOA estimation results. When $\alpha \ge 1$, Sigmoid-MUSIC and $l_p$-MUSIC could obtain smaller RMSEs than other methods. When $\alpha =2$, the $\alpha$ stable distribution became a Gaussian distribution. As shown in Figure 8, in this case, the three algorithms had similar performance. As illustrated in Figure 8, root CRB decreased with increasing characteristic exponent $\alpha$; however, root CRB was less affected by $\alpha$ for a given GSNR.

7. Conclusions

In this paper, we considered the problem of parameter estimation for a wideband bistatic MIMO radar in impulsive noise environments. Based on the WBAF and MUSIC methods, we developed a novel method without any a priori knowledge of the noise. A robust estimator, termed Sigmoid-WBAF, was employed to estimate the time delay and Doppler stretch in the presence of impulsive noise. Then, the properties of the Sigmoid transform and Sigmoid correlation were presented. A novel MUSIC-based Sigmoid correlation (Sigmoid-MUSIC) was developed to estimate DOD and DOA. Furthermore, the boundness of Sigmoid-WBAF to the symmetric alpha stable ($S\alpha S$) noise, the feasibility analysis of Sigmoid-WBAF, and complexity analysis of Sigmoid-WBAF and Sigmoid-MUSIC algorithms were presented to evaluate the performance of the proposed method. In addition, the Cramér–Rao bound for parameter estimation was derived and computed in closed form, which showed that better performance was achieved. Lastly, comprehensive simulations were carried out to evaluate the performance of different methods. Simulation results and theoretical analyses showed that the proposed Sigmoid-WBAF and Sigmoid-MUSIC had better estimation performance, especially in highly impulsive noise environments.