0 Introduction

Multiple-input multiple-output (MIMO) radar has drawn considerable attentions in radar signal processing field^[1-5]. According to the antennas configuration, MIMO radars can be categorized into two types. One is the case that the antennas are closely-spaced in both transmit array and receive array to achieve unambiguous angle estimation. The other is that the antennas are widely-spaced in both arrays to obtain coherent processing gain for solving target scintillation problem^[5]. Compared with conventional phase array radar, MIMO radar has more degrees of freedom and thus it has better interference rejection capability, parameter identifiability and high spatial resolution^[6-9]. The problem of joint estimation of direction of arrival (DOA) and Doppler frequency has become a research hotspot. Joint estimation algorithms are usually extended from one dimensional (1D) parameter estimation algorithms, such as multiple signal classification (MUSIC)^[10], estimation of signal parameters via rotational invariance technique (ESPRIT)^[11], propagator method (PM)^[12], trilinear decomposition^[13], etc. MUSIC algorithm estimates parameters via spectral peak searching and needs a large computational complexity. Ref.[14] proposed a root-MUSIC algorithm which achieves estimation with polynomial rooting instead of spectral search. ESPRIT algorithm^[11] utilized the rotational invariance to achieve automatically paired estimations of the angles. Refs.[15-16] discussed the application of ESPRIT in MIMO radar field. The PM^[12] has a significantly lower complexity since it needs no spectral peak searching nor eigen-value decomposition. The propagator is a linear operator which only depends on steering vectors and can be easily extracted from the received data. But it has a relatively worse performance, especially in low signal-to-noise ratio (SNR). Parallel factor (PARAFAC)^[17] and trilinear decomposition^[18] are iterative methods and can achieve better estimation accuracy than ESPRIT. Refs.[19-23] introduced several kinds of joint angle and Doppler frequency estimation algorithms. Ref.[19] was based on ESPRIT and can obtain automatically paired angle and Doppler frequency estimation. Ref.[20] used the rotational factor produced by time delay sampling. Ref.[21] presented the efficient joint angle and frequency estimation algorithm based on the PM for uniform linear arrays (ULA). Ref.[22] achieved the joint estimation via linking the estimation problem to the trilinear model and achieve better performance. Ref.[23] exploited the quadrilinear model and used the quadrilinear alternating least squares method to obtain joint estimation of angle and Doppler frequency.

Ref.[24] proposed a reduced-dimension MUSIC (RD-MUSIC) method for joint direction of departure (DOD) and DOA estimation. Though Ref.[24] discussed a different problem, it provides us an idea to overcome the shortcoming of conventional MUSIC algorithm for joint DOA and Doppler frequency estimation in a monostatic MIMO radar. Since conventional MUSIC algorithm requires an exhaustive two dimensional (2D) spectrum search, it needs high computational complexity. In this paper, we propose a computationally efficient MUSIC (CE-MUSIC) algorithm for joint DOA and Doppler frequency estimation in a monostatic MIMO radar. We firstly exploit a reduced-dimension transformation as a pre-processing step which can remove the redundant entries of steering matrix and further reduce the computational load. Then we extend the method in Ref.[24] to joint DOA and Doppler frequency problem. By means of the pre-processing step and 1D search, the proposed CE-MUSIC algorithm overcomes the shortcoming of heavy computational load in conventional 2D-MUSIC and provides very close estimation performance to it with much lower complexity. We also derive the mean square error (MSE) of the proposed CE-MUSIC algorithm and Cramer-Rao bound (CRB) for joint DOA and Doppler frequency estimation in monostatic MIMO radar.

1 Data Model

As shown in Fig. 1, we consider a monostatic MIMO radar system equipped with both ULAs for its transmit/receive array, which contains M/N elements respectively with half-wave length spacing between adjacent antennas. At the transmitter array, all elements emit different coded continues periodic signals simultaneously, which are orthogonal and have identical bandwidth and center frequency.

It is assumed that the Doppler frequencies have almost no effect on the orthogonality of the waveforms and the variance of the phase within repetition intervals. Each element of the transmitter emits orthogonal waveforms. There are K independent far-field targets with different Doppler frequencies. Since the transmit and receive array are close to each other, we assume that the DOA and DOD are of the same angle. For the kth source (k=1, …, K), θ_k is denoted to be DOA.

The output of the matched filters for the receive sensors can be expressed as ^[24]

$ x\left( t \right) = \mathit{\boldsymbol{As}}\left( t \right) + \mathit{\boldsymbol{n}}\left( t \right) $

(1)

where A=[a_r(θ₁)⊗a_t(θ₁), …, a_r(θ_K)⊗a_t(θ_K)] is a MN×K matrix composed of K MIMO steering vectors. a_t(θ_k)=[1, e^{－jπsinθ_k}, …, e^{－j(M－1)πsinθ_k}]^T] is the steering vector of the transmit array, a_r(θ_k)=[1, e^{－jπsinθ_k}, …, e^{－j(N－1)πsinθ_k}]^T is the steering vector of the receive array and ⊗ is the Kronecker product. s(t)=[s₁(t), s₂(t), …, s_K(t)]^T∈C^K×1, s_k(t)=β_ke^j2πf_kt/f_s with β_k, f_k and f_s being the radar cross section (RCS) fading coefficient, Doppler frequency of the kth target and the pulse repeat frequency, respectively. n(t)∈C^MN×1 is a noise vector assumed to be independent, zero-mean complex Gaussian distribution with covariance matrix σ²I.

For Doppler frequency estimation, we collect J snapshots, where the time interval between two adjacent snapshots is fixed. The received signals can be given by

$ \mathit{\boldsymbol{X}} = \left[ {\mathit{\boldsymbol{x}}\left( {{t_1}} \right),\mathit{\boldsymbol{x}}\left( {{t_2}} \right), \cdots ,\mathit{\boldsymbol{x}}\left( {{t_J}} \right)} \right] $

(2)

We divide the J sampled signals into P sections and each section contains L=J－P+1 snapshots. Then the set of all sections Y can be denoted as

$ \mathit{\boldsymbol{Y}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}\\ {{\mathit{\boldsymbol{X}}_2}}\\ \vdots \\ {{\mathit{\boldsymbol{X}}_P}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {x\left( {{t_1}} \right),x\left( {{t_2}} \right), \cdots ,x\left( {{t_L}} \right)}\\ {x\left( {{t_2}} \right),x\left( {{t_3}} \right), \cdots ,x\left( {{t_{L + 1}}} \right)}\\ \vdots \\ {x\left( {{t_P}} \right),x\left( {{t_{P + 1}}} \right), \cdots ,x\left( {{t_J}} \right)} \end{array}} \right] \in {{\bf{C}}^{MNP \times L}} $

(3)

Define a P×K matrix

$ \mathit{\boldsymbol{\psi }} = \left[ {\begin{array}{*{20}{c}} 1&1& \cdots &1\\ {{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\tau {f_1}}}}&{{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\tau {f_2}}}}& \cdots &{{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\tau {f_K}}}}\\ \vdots&\vdots&\vdots&\vdots \\ {{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\left( {P - 1} \right)\tau {f_1}}}}&{{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\left( {P - 1} \right)\tau {f_2}}}}& \cdots &{{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\left( {P - 1} \right)\tau {f_K}}}} \end{array}} \right] $

where τ denotes the time delay between two adjacent snapshots. Then for the pth section (p=1, 2, …, P), we have

$ {\mathit{\boldsymbol{X}}_p} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{D}}_p}\left( \mathit{\boldsymbol{\psi }} \right)\mathit{\boldsymbol{S}} + {\mathit{\boldsymbol{N}}_p} $

(4)

where S=[s(t₁), s(t₂), …, s(t_L)]∈C^K×L, N_p∈C^MN×L is the noise matrix corresponding to the pth section, and D_p(ψ) denotes a diagonal matrix whose diagonal elements are the p-th row of ψ. Y can be compactly written as

$ \mathit{\boldsymbol{Y}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}\\ {{\mathit{\boldsymbol{X}}_2}}\\ \vdots \\ {{\mathit{\boldsymbol{X}}_P}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{D}}_1}\left( \mathit{\boldsymbol{\psi }} \right)}\\ {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{D}}_2}\left( \mathit{\boldsymbol{\psi }} \right)}\\ \vdots \\ {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{D}}_P}\left( \mathit{\boldsymbol{\psi }} \right)} \end{array}} \right]\mathit{\boldsymbol{S}} + \left[ {\begin{array}{*{20}{c}} {{N_1}}\\ {{N_2}}\\ \vdots \\ {{N_P}} \end{array}} \right] = \left[ {\mathit{\boldsymbol{\psi }} \circ \mathit{\boldsymbol{A}}} \right]\mathit{\boldsymbol{S}} + \mathit{\boldsymbol{N'}} $

(5)

where N′=[N₁]^T, N₂^T, …, N_P^TT denotes the noise part of Y and "°" is the Khatri_Rao product.

The covariance matrix R=E[YY^H] can be decomposed as^[11]

$ \mathit{\boldsymbol{R}} = {\mathit{\boldsymbol{E}}_{\rm{s}}}{\mathit{\boldsymbol{D}}_{\rm{s}}}\mathit{\boldsymbol{E}}_{\rm{s}}^{\rm{H}} + {\mathit{\boldsymbol{E}}_{\rm{n}}}{\mathit{\boldsymbol{D}}_{\rm{n}}}\mathit{\boldsymbol{E}}_{\rm{n}}^{\rm{H}} $

(6)

where D_s denotes a K×K diagonal matrix formed by K largest eigenvalues, and D_n a diagonal matrix formed by remaining (MNP-K) smallest eigenvalues. E_s is the matrix composed of the eigenvectors corresponding to the K largest eigenvalues. E_n includes the remaining eigenvectors. E_s and E_n represent the signal subspace and noise subspace, respectively.

For the signal model in Eq.(1), R can be estimated by $\mathit{\boldsymbol{\hat R}} = \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{Y}}^{\text{H}}}/L$.

2 Computationally Efficient MUSIC (CE-MUSIC) Algorithm for Joint DOA and Doppler Frequency Estimation 2.1 Conventional 2D-MUSIC algorithm

Since the targets are assumed to be incoherent, it is straightforward to know that the columns of [ψ°A] are orthogonal to the noise subspace E_n^[10]. Then the 2D-MUSIC spatial spectrum function for joint DOA and Doppler frequency estimation can be obtained as ^[24]

$ \begin{array}{*{20}{c}} {{F_{{\rm{2D - MUSIC}}}}\left( {\theta ,f} \right) = }\\ {\frac{1}{{{{\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) \otimes a\left( \theta \right)} \right]}^{\rm{H}}}{\mathit{\boldsymbol{E}}_{\rm{n}}}\mathit{\boldsymbol{E}}_{\rm{n}}^{\rm{H}}\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) \otimes \mathit{\boldsymbol{a}}\left( \theta \right)} \right]}}} \end{array} $

(7)

where Ω(f)=[1, e^j2πτf/f_s, …, e^{j2π(P－1)τf/f_s}]^T∈C^P×1, a(θ)=a_r(θ)⊗a_t(θ)∈C^MN×1, a_r(θ)=[1, e^－jπsinθ, …, e^{－j(N－1)πsinθ}]^T, a_t(θ)=[1, e^－jπsinθ, …, e^{－j(M－1)πsinθ}]^T.

After 2D spectral peak searching of f, θ, we take the K largest peaks of F_2D-MUSIC(θ, f) to obtain the DOA and Doppler frequency estimation. However, the 2D search results in large computational complexity. Aiming at this insufficiency, we propose the CE-MUSIC algorithm, which uses a reduced-dimension transformation and only needs 1D search, with much lower complexity.

2.2 CE-MUSIC algorithm for joint DOA and Doppler frequency estimation 2.2.1 Reduced-dimension transformation

Define

$ {\mathit{\boldsymbol{a}}_{\rm{r}}}\left( {{\theta _k}} \right) \otimes {\mathit{\boldsymbol{a}}_t}\left( {{\theta _k}} \right) = \mathit{\boldsymbol{Gb}}\left( {{\theta _k}} \right) $

(8)

where a_r(θ_k)⊗a_t(θ_k) is the steering vector of the kth source and b(θ_k)=[1, e^{－jπsinθ_k}, …, e^{－jπ(M+N－2)sinθ_k}]^T∈C^(M+N－1)×1. G∈C^{MN×(M+N－1)} can be expressed as^[25]

$ \mathit{\boldsymbol{G}} = \left[ {\left. {\begin{array}{*{20}{c}} {\left. {\begin{array}{*{20}{c}} 1&0& \cdots &0&0&0& \cdots &0\\ 0&1& \cdots &0&0&0& \cdots &0\\ \vdots&\vdots&\vdots&\vdots&\vdots &{}& \vdots &{}\\ 0&0& \cdots &1&0&0& \cdots &0 \end{array}} \right\}M}\\ {\left. {\begin{array}{*{20}{c}} 0&1&0& \cdots &0&0& \cdots &0\\ 0&0&1& \cdots &0&0& \cdots &0\\ \vdots&\vdots &{}& \ddots&\vdots&\vdots&\vdots &{}\\ 0&0&0& \cdots &0&1& \cdots &0 \end{array}} \right\}M}\\ \vdots \\ {\left. {\begin{array}{*{20}{c}} 0& \cdots &0&1&0&0& \cdots &0\\ 0& \cdots &0&0&1&0& \cdots &0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots &{}\\ 0& \cdots &0&0&0&0& \cdots &1 \end{array}} \right\}M} \end{array}} \right\}M \times N} \right] $

(9)

Define W=G^HG, which is given as

$ \mathit{\boldsymbol{W}} = {\rm{diag}}\left( {1, \cdots ,\underbrace {\min \left( {M,N} \right), \cdots ,\min \left( {M,N} \right)}_{\left| {M - N} \right| + 1}, \cdots ,1} \right) $

(10)

We use the reduced-dimension transformation W^－1/2G^H for the received signal and can obtain

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{x'}}\left( t \right) = {\mathit{\boldsymbol{W}}^{ - 1/2}}{\mathit{\boldsymbol{G}}^{\rm{H}}}\mathit{\boldsymbol{x}}\left( t \right) = {\mathit{\boldsymbol{W}}^{ - 1/2}}}\\ {{\mathit{\boldsymbol{G}}^{\rm{H}}}\left( {\mathit{\boldsymbol{G}}\left[ {\mathit{\boldsymbol{b}}\left( {{\theta _1}} \right), \cdots ,\mathit{\boldsymbol{b}}\left( {{\theta _K}} \right)} \right]\mathit{\boldsymbol{s}}\left( t \right) + \mathit{\boldsymbol{n}}\left( t \right)} \right) = }\\ {{\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{Bs}}\left( t \right) + {\mathit{\boldsymbol{W}}^{ - 1/2}}{\mathit{\boldsymbol{G}}^{\rm{H}}}\mathit{\boldsymbol{n}}\left( t \right)} \end{array} $

(11)

where B=[b(θ₁), …, b(θ_K)]∈C^(M+N－1)×K is a Vandermonde matrix.

This transformation can remove the redundant entries of the steering matrix and reserve the distinct entries meanwhile, which means it reserves all the effective information and avoids performance loss. Since the reduced-dimension transformation matrix is sparse, it only adds a minor computational load. From Ref. [25], we find that the reduced-dimension transformation does not cause additional coloring and can also obtain SNR gain.

2.2.2 Joint DOA and Doppler frequency estimation using only one-dimensional search

After using the reduced-dimension transformation, the received signal can be transformed as

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_{{\rm{rd}}}} = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{X'}}}_1}}\\ {{{\mathit{\boldsymbol{X'}}}_2}}\\ \vdots \\ {{{\mathit{\boldsymbol{X'}}}_P}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x'}}\left( {{t_1}} \right),\mathit{\boldsymbol{x'}}\left( {{t_2}} \right), \cdots ,\mathit{\boldsymbol{x'}}\left( {{t_L}} \right)}\\ {\mathit{\boldsymbol{x'}}\left( {{t_2}} \right),\mathit{\boldsymbol{x'}}\left( {{t_3}} \right), \cdots ,\mathit{\boldsymbol{x'}}\left( {{t_{L + 1}}} \right)}\\ \vdots \\ {\mathit{\boldsymbol{x'}}\left( {{t_P}} \right),\mathit{\boldsymbol{x'}}\left( {{t_{P + 1}}} \right), \cdots ,\mathit{\boldsymbol{x'}}\left( {{t_J}} \right)} \end{array}} \right] \in }\\ {{{\bf{C}}^{\left( {M + N - 1} \right)P \times L}}} \end{array} $

(12)

For the pth section (p=1, 2, …, P), it can be expressed as

$ {{\mathit{\boldsymbol{X'}}}_P} = {\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_p}\left( \mathit{\boldsymbol{\psi }} \right)\mathit{\boldsymbol{S}} + {{\mathit{\boldsymbol{N'}}}_p} $

(13)

where N′_p=W^1/2B[n(t_p), n(t_p+1), …, n(t_L－1+p)]^T is the noise of the corresponding section after reduced-dimension transformation.

Then X_rd can be written clearly as

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_{{\rm{rd}}}} = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{X'}}}_1}}\\ {{{\mathit{\boldsymbol{X'}}}_2}}\\ \vdots \\ {{{\mathit{\boldsymbol{X'}}}_P}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_1}\left( \mathit{\boldsymbol{\psi }} \right)}\\ {{\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_2}\left( \mathit{\boldsymbol{\psi }} \right)}\\ \vdots \\ {{\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_P}\left( \mathit{\boldsymbol{\psi }} \right)} \end{array}} \right]\mathit{\boldsymbol{S}} + \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{N'}}}_1}}\\ {{{\mathit{\boldsymbol{N'}}}_2}}\\ \vdots \\ {{{\mathit{\boldsymbol{N'}}}_P}} \end{array}} \right] = }\\ {\left[ {\mathit{\boldsymbol{\psi }} \circ \left( {{\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{B}}} \right)} \right]\mathit{\boldsymbol{S}} + \mathit{\boldsymbol{N'}} = \mathit{\boldsymbol{B'S}} + \mathit{\boldsymbol{N'}}} \end{array} $

(14)

where B′=ψ°(W^1/2B)∈C^{(M+N－1)P×K}, and N′ denotes the corresponding noise. The noise subspace of X_rd is obtained in the same way as Eq.(6) and E′_n is used to denote it. Similarly from Eq.^[10], the columns of ψ°(W^1/2B) are orthogonal to E′_n. Thus, the spatial spectrum function can be denoted as

$ \begin{array}{l} {F_{{\rm{RD - MUSIC}}}}\left( {\theta ,f} \right) = \\ \frac{1}{{{{\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) \otimes \mathit{\boldsymbol{b'}}\left( \theta \right)} \right]}^{\rm{H}}}{{\mathit{\boldsymbol{E'}}}_{\rm{n}}}\mathit{\boldsymbol{E'}}_{\rm{n}}^{\rm{H}}\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) \otimes \mathit{\boldsymbol{b'}}\left( \theta \right)} \right]}} \end{array} $

(15)

where b′(θ)=W^1/2b(θ).

We define

$ \begin{array}{l} \mathit{\boldsymbol{W}}\left( {\theta ,f} \right) = \\ {\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) \otimes \mathit{\boldsymbol{b'}}\left( \theta \right)} \right]^{\rm{H}}}{{\mathit{\boldsymbol{E'}}}_{n}}\mathit{\boldsymbol{E'}}_{\rm{n}}^{\rm{H}}\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) \otimes \mathit{\boldsymbol{b'}}\left( \theta \right)} \right] \end{array} $

(16)

which can be also denoted as

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}\left( {\theta ,f} \right) = \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{{\left( f \right)}^{\rm{H}}}{{\left[ {{\mathit{\boldsymbol{I}}_P} \otimes \mathit{\boldsymbol{b'}}\left( \theta \right)} \right]}^{\rm{H}}}}\\ {{{\mathit{\boldsymbol{E'}}}_{\rm{n}}}\mathit{\boldsymbol{E'}}_{\rm{n}}^{\rm{H}}\left[ {{\mathit{\boldsymbol{I}}_P} \otimes \mathit{\boldsymbol{b'}}\left( \theta \right)} \right]\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) = \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{{\left( f \right)}^{\rm{H}}}\mathit{\boldsymbol{Q}}\left( \theta \right)\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right)} \end{array} $

(17)

where Q(θ)=[I_P⊗b′(θ)]^HE′_nE′_n^H[I_P⊗b′(θ)].

Eq.(17) can be regarded as an optimization problem with the constraint of e₁^HΩ(f)=1, where e₁=[1, 0, …, 0]^T∈C^P×1. According to Ref. [24], the optimization problem of Eq.(18) can be reconstructed as

$ \mathop {\min }\limits_\theta \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\left( f \right)^{\rm{H}}}\mathit{\boldsymbol{Q}}\left( \theta \right)\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right),{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{e}}_1^{\rm{H}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) = 1 $

(18)

From Eq.(7), the costing function is constructed as L(θ, f)=Ω(f)^HQ(θ)Ω(f)^H－λ(e₁^HΩ(f)－1), where λ is a constant. We have

$ \frac{{\partial L\left( {\theta ,f} \right)}}{{\partial \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right)}} = 2\mathit{\boldsymbol{Q}}\left( \theta \right)\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) + \lambda {\mathit{\boldsymbol{e}}_1} = 0 $

(19)

From Eq.(19), we get Ω(f)=μQ^－1(θ)e₁, where μ is a constant. From e₁^HΩ(f)=1, $\mu = \frac{1}{{\mathit{\boldsymbol{e}}_1^{\text{H}}{Q^{-1}}\left( \theta \right){\mathit{\boldsymbol{e}}_1}}}$. Then Ω(f) can be obtained via

$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right) = \frac{{{\mathit{\boldsymbol{Q}}^{ - 1}}\left( \theta \right){\mathit{\boldsymbol{e}}_1}}}{{\mathit{\boldsymbol{e}}_1^{\rm{H}}{\mathit{\boldsymbol{Q}}^{ - 1}}\left( \theta \right){\mathit{\boldsymbol{e}}_1}}} $

(20)

Substituting Eq.(20) into $\mathop {\min }\limits_\theta \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\left( f \right)^{\text{H}}}\mathit{\boldsymbol{Q}}\left( \theta \right)\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( f \right)$, the DOA can be estimated via

$ \hat \theta = \arg \mathop {\min }\limits_\theta \frac{1}{{\mathit{\boldsymbol{e}}_1^{\rm{H}}{\mathit{\boldsymbol{Q}}^{ - 1}}\left( \theta \right){\mathit{\boldsymbol{e}}_1}}}\arg \mathop {\max }\limits_\theta \mathit{\boldsymbol{e}}_1^{\rm{H}}{\mathit{\boldsymbol{Q}}^{ - 1}}\left( \theta \right){\mathit{\boldsymbol{e}}_1} $

(21)

Considering e₁=[1, 0, …, 0]^T∈C^P×1, we can search θ∈[－90°, 90°] and estimate DOAs by finding K largest peaks of the (1, 1) element of Q^－1(θ). We use (${\hat \theta _1}, \cdots, {\hat \theta _K}$) and form K vectors $\mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_1}} \right), \cdots, \mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_K}} \right)$ according to Eq.(20). For Ω(f)=[1, e^j2πf/f_s, …, e^{j2π(P－1)f/f_s}]^T, let g_k denotes the phase angle of Ω(f_k), k=1, 2, …, K

$ {\mathit{\boldsymbol{g}}_k} = {\rm{angle}}\left( {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( {{f_k}} \right)} \right) $

(22)

Then we get g_k=[0, 2πf_k/f_s, …, 2fπf(P－1)f_k/f_s]^T and adopt the least square (LS) principle to estimate f_k. After the normalization for $\mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_1}} \right), \cdots, \mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_K}} \right)$, we obtain the estimated ${\mathit{\boldsymbol{\hat g}}_k}$, k=1, 2, …, K. Then the LS principle is used to estimate the Doppler frequency f_k. LS fitting is

$ \mathop {\min }\limits_{{c_k}} \left\| {\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{c}}_k} - {{\mathit{\boldsymbol{\hat g}}}_k}} \right\|_F^2 $

(23)

where c_k=[c_k0, c_k1]^T∈C^2×1 is an unknown parameter vector, and c_k1 the estimated value of f_k/f_s. P is given by

$ \mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} 1&0\\ 1&{2{\rm{ \mathsf{ π} }}\tau {f_k}/{f_{\rm{s}}}}\\ \vdots&\vdots \\ 1&{2{\rm{ \mathsf{ π} }}\left( {P - 1} \right)\tau {f_k}/{f_{\rm{s}}}} \end{array}} \right] $

(24)

The LS solution for c_k is ${\left[{{{\hat c}_{k0}}, {{\hat c}_{k1}}} \right]^{\text{T}}} = {\left( {{\mathit{\boldsymbol{P}}^{\text{T}}}\mathit{\boldsymbol{P}}} \right)^{ -1}}{\mathit{\boldsymbol{P}}^{\text{T}}}{\mathit{\boldsymbol{\hat g}}_k}$, and the Doppler frequency f_k can be estimated via

$ {{\hat f}_k} = {{\hat c}_{k1}} \cdot {f_{\rm{s}}} $

(25)

Till now, we have proposed the CE-MUSIC algorithm for joint DOA and Doppler frequency estimation in monostatic MIMO radar. The major steps and corresponding computational complexity of our algorithms are listed as follows:

(1) Form the received data matrix Y as Eq.(3), performing reduced-dimension transformation and obtaining noise subspace E′_n.O{((M+N－1)P)²(L+(M+N－1)P)};

(2) Searching θ∈[－90°, 90°] and getting the estimate of DOA by finding K largest peaks of the (1, 1) element of Q^－1(θ) as Eq.(21).O{n[P²((M+N－1)P－K)(M+N)+(M+N－1)²]}, where n is the total search times within the search range;

(3) According to Eq.(20), form K vectors $\mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_1}} \right), \cdots, \mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_K}} \right)$, getting the estimate of Doppler frequency via LS principle. O {K(P³+P²+10P)}.

2.2.3 Complexity analysis

The complexity of CE-MUSIC, 2D-MUSIC and ESPRIT^[19] algorithm have been shown in Table 1, where H=M+N－1, n is the total search times for CE-MUSIC within the searching range and n_f, n_DOAare the search times that 2D-MUSIC needs for Doppler frequency and DOA estimation respectively. Our proposed algorithm has much lower computational complexity than 2D-MUSIC and higher complexity than ESPRIT and this property has been proven in Fig. 2. (Herein we assume n=n_f=n_DOA=3 000 and the other parameters have been shown in each title). Notably, with the number of search times increasing, the improvement will become greater.

2.2.4 Advantages of the proposed CE-MUSIC algorithm

The proposed CE-MUSIC algorithm has advantages as follows:

(1) It only needs 1D search and has much lower complexity than conventional 2D-MUSIC algorithm for joint DOA and Doppler frequency estimation.

(2) It has very close joint DOA and Doppler frequency estimation performance to the conventional 2D-MUSIC algorithm, which will be shown in detail in Section 4.

(3) It has much better joint estimation performance than the ESPRIT algorithm, which will be proved in detail in Section 4.

(4) It can obtain automatically paired DOA and Doppler frequency estimation.

Remark A  In this paper, the number of source signals K is assumed to be known. If not, it can be pre-estimated via the method in Ref.[26].

Remark B   From Eq.(20), after searching θ∈[－90°, 90°] and finding K peaks of the (1, 1) element of Q^－1(θ), the corresponding $\mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_1}} \right), \cdots, \mathit{\boldsymbol{ \boldsymbol{\hat \varOmega} }}\left( {{f_K}} \right)$ are determined. The one-to-one correspondence ensures that the DOA and Doppler frequency estimation can be automatically paired.

Remark C   There are some differences between the proposed algorithm and the RD-MUSIC algorithm^[24]. The RD-MUSIC algorithm works for two-dimensional angle estimation in a bistatic MIMO radar, whereas the proposed algorithm is used for joint DOA and Doppler frequency estimation in a monostatic MIMO radar. Since the proposed algorithm is studied for a different problem, it uses a different spectrum function from that in Ref. [24]. We derive the new function and reduce the complexity significantly, while not debasing the performance of joint DOA and Doppler frequency estimation.

3 Performance Analysis

The aim of this section is to analyze estimation performance of the proposed CE-MUSIC algorithm. We derive the mean square error of joint DOA and Doppler frequency estimation, and propose the Cramer-Rao bound (CRB) of joint estimation in a monostatic MIMO radar system.

3.1 Error analysis

Define^[27]

$ \left[ {{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}\left| {{{\mathit{\boldsymbol{E'}}}_{\rm{n}}}} \right.} \right] = \left[ {{\mathit{\boldsymbol{s}}_1}, \cdots ,{\mathit{\boldsymbol{s}}_K}\left| {{\mathit{\boldsymbol{g}}_1}, \cdots ,{\mathit{\boldsymbol{g}}_{\left( {M + N - 1} \right)P - K}}} \right.} \right] $

(26)

$ \left[ {{{\mathit{\boldsymbol{\hat E'}}}_{\rm{s}}}\left| {{{\mathit{\boldsymbol{\hat E'}}}_{\rm{n}}}} \right.} \right] = \left[ {{{\mathit{\boldsymbol{\hat s}}}_1}, \cdots ,{{\mathit{\boldsymbol{\hat s}}}_K}\left| {{{\hat g}_1}, \cdots ,{{\mathit{\boldsymbol{\hat g}}}_{\left( {M + N - 1} \right)P - K}}} \right.} \right] $

(27)

as the matrices composed of eigenvectors associated with the eigenvalues {λ_i}_i=1^(M+N－1)P of reduced-dimension transformed signal covariance matrix R_rd and {λ_i}_i=1^(M+N－1)P of ${\mathit{\boldsymbol{\hat R}}_{{\text{rd}}}}$ in ascending order, respectively. We have

$ \begin{array}{*{20}{c}} {E\left[ {\left( {{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}\mathit{\boldsymbol{E'}}_{\rm{s}}^{\rm{H}}{\mathit{\boldsymbol{g}}_i}} \right){{\left( {{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}{\mathit{\boldsymbol{g}}_j}} \right)}^{\rm{H}}}} \right] = }\\ {\frac{{{\sigma ^2}}}{L}\left[ {\sum\limits_{k = 1}^K {\frac{{{\lambda _k}}}{{{{\left( {{\sigma ^2} - {\lambda _k}} \right)}^2}}}{\mathit{\boldsymbol{s}}_k}\mathit{\boldsymbol{s}}_k^{\rm{H}}} } \right]{\delta _{i,j}} = \frac{1}{L}U{\delta _{i,j}}{\sigma ^2}} \end{array} $

(28)

$ E\left[ {\left( {{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}\mathit{\boldsymbol{E'}}_{\rm{s}}^{\rm{H}}{{\mathit{\boldsymbol{\hat g}}}_i}} \right){{\left( {{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}{{\mathit{\boldsymbol{\hat g}}}_j}} \right)}^{\rm{T}}}} \right] = 0 $

(29)

where ${\delta _{i, j}} = \left\{ \begin{gathered} 1\;\;\;\;i = j \hfill \\ 0\;\;\;\;{\text{other}} \hfill \\ \end{gathered} \right.$, $\mathit{\boldsymbol{U}} = \sum\limits_{k = 1}^K {\frac{{{\lambda _k}}}{{{{\left( {{\sigma ^2}-{\lambda _k}} \right)}^2}}}} {\mathit{\boldsymbol{s}}_k}\mathit{\boldsymbol{s}}_k^{\text{H}}$, σ² is the power of the noise.

Define $\mathit{\boldsymbol{r}} \triangleq {\left[{u, v} \right]^{\text{T}}}$ and ${\mathit{\boldsymbol{v}}_k} \triangleq \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\left( {{u_k}} \right) \otimes \left( {{\mathit{\boldsymbol{W}}^{1/2}}\mathit{\boldsymbol{b}}\left( {{v_k}} \right)} \right)$. The cost function becomes

$ w\left( \mathit{\boldsymbol{r}} \right) = \mathit{\boldsymbol{v}}_k^{\rm{H}}{{\mathit{\boldsymbol{E'}}}_{\rm{n}}}\mathit{\boldsymbol{E'}}_{\rm{n}}^{\rm{H}}{\mathit{\boldsymbol{v}}_k} $

(30)

Define ${\nabla _r} \triangleq {\left[{\partial /\partial u, \partial /\partial v} \right]^{\text{T}}}$, $\mathit{\boldsymbol{d}}\left( {{u_k}} \right) = \partial {\mathit{\boldsymbol{v}}_k}/\partial {u_k}$ and $\mathit{\boldsymbol{d}}\left( {{v_k}} \right) = \partial {\mathit{\boldsymbol{v}}_k}/\partial {v_k}$. For $\left\{ {{{\mathit{\boldsymbol{\hat r}}}_i}} \right\}$ is the minimum point of w(r), we have $w'\left( {{{\mathit{\boldsymbol{\hat r}}}_i}} \right) = 0$, Using a first order Taylor series expansion, we have ^[28]

$ \mathit{\boldsymbol{0}} \cong {\nabla _r}w\left( \mathit{\boldsymbol{r}} \right)\left| {_{{\mathit{\boldsymbol{r}}_i}}} \right. + {\nabla _r}{\left( {{\nabla _r}w\left( \mathit{\boldsymbol{r}} \right)} \right)^{\rm{T}}}\left| {_{{\mathit{\boldsymbol{r}}_{i\xi }}}} \right.\left( {{{\mathit{\boldsymbol{\hat r}}}_i} - {\mathit{\boldsymbol{r}}_i}} \right) $

(31)

where $\cong $ is a symbol used to denote items that are approximately equal. r_iξ is a vector of some value on the line segment joining r_i and ${\mathit{\boldsymbol{\hat r}}_i}$. According to Eq.(31), we get the estimation error of the vector r_i as

$ {{\mathit{\boldsymbol{\hat r}}}_i} - {\mathit{\boldsymbol{r}}_i} = \frac{{ - {\nabla _r}w\left( \mathit{\boldsymbol{r}} \right)\left| {_{{\mathit{\boldsymbol{r}}_i}}} \right.}}{{{\nabla _r}{{\left( {{\nabla _r}w\left( \mathit{\boldsymbol{r}} \right)} \right)}^{\rm{T}}}\left| {_{{\mathit{\boldsymbol{r}}_{i\xi }}}} \right.}} $

(32)

We define ${\mathit{\boldsymbol{H}}_i} = \mathop {\lim }\limits_{N \to \infty } {\nabla _r}{\left( {{\nabla _r}w\left( \mathit{\boldsymbol{r}} \right)} \right)^{\text{T}}}{|_{{\mathit{\boldsymbol{r}}_i}}}$, and the asymptotic covariance matrices is

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{ik}} = \mathop {\lim }\limits_{N \to \infty } E\left[ {\left( {{{\mathit{\boldsymbol{\hat r}}}_i} - {\mathit{\boldsymbol{r}}_i}} \right){{\left( {{{\mathit{\boldsymbol{\hat r}}}_k} - {\mathit{\boldsymbol{r}}_k}} \right)}^{\rm{T}}}} \right] = }\\ {{{\left( {{\mathit{\boldsymbol{H}}_i}} \right)}^{ - 1}}\mathop {\lim }\limits_{N \to \infty } E\left[ {{\nabla _{{\mathit{\boldsymbol{r}}_i}}}w\left( \mathit{\boldsymbol{r}} \right){{\left( {{\nabla _{{\mathit{\boldsymbol{r}}_k}}}w\left( \mathit{\boldsymbol{r}} \right)} \right)}^{\rm{T}}}} \right]\left| {_{{\mathit{\boldsymbol{r}}_i},{\mathit{\boldsymbol{r}}_k}}} \right.{{\left( {{\mathit{\boldsymbol{H}}_k}} \right)}^{ - 1}}} \end{array} $

(33)

We have

$ \begin{array}{l} \mathit{\boldsymbol{a}}_k^{\rm{H}}{{\mathit{\boldsymbol{{\hat E}'}}}_{\rm{n}}}\mathit{\boldsymbol{{\hat E}'}}_{\rm{n}}^{\rm{H}}\mathit{\boldsymbol{d}}\left( {{u_k}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\sum\limits_{p = 1}^{\left( {M + N - 1} \right)P - K} {\left[ {\mathit{\boldsymbol{g}}_p^{\rm{H}}\mathit{\boldsymbol{d}}\left( {{u_k}} \right))} \right]\left[ {\mathit{\boldsymbol{a}}_k^{\rm{H}}{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}\mathit{\boldsymbol{E'}}_{\rm{s}}^{\rm{H}}{{\mathit{\boldsymbol{\hat g}}}_p}} \right]} \end{array} $

(34)

$ \begin{array}{l} \mathit{\boldsymbol{a}}_k^{\rm{H}}{{\mathit{\boldsymbol{{\hat E}'}}}_{\rm{n}}}\mathit{\boldsymbol{{\hat E}'}}_{\rm{n}}^{\rm{H}}\mathit{\boldsymbol{d}}\left( {{v_k}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\sum\limits_{p = 1}^{\left( {M + N - 1} \right)P - K} {\left[ {\mathit{\boldsymbol{g}}_p^{\rm{H}}\mathit{\boldsymbol{d}}\left( {{v_k}} \right)} \right]\left[ {\mathit{\boldsymbol{a}}_k^{\rm{H}}{{\mathit{\boldsymbol{E'}}}_{\rm{s}}}\mathit{\boldsymbol{E'}}_{\rm{s}}^{\rm{H}}{{\mathit{\boldsymbol{\hat g}}}_p}} \right]} \end{array} $

(35)

So the covariance matrix of error estimation can be obtained as

$ \begin{array}{l} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{ik}} = \frac{{{\sigma ^2}}}{{2Lw\left( i \right)w\left( k \right)}}\\ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{d}}{{\left( {{v_i}} \right)}^{\rm{H}}}\left( i \right)\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_i}} \right)}&{\rho \left( i \right)}\\ {\rho \left( i \right)}&{\mathit{\boldsymbol{d}}{{\left( {{u_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_i}} \right)} \end{array}} \right] \times \\ {\mathop{\rm Re}\nolimits} \left\{ {\left( {\mathit{\boldsymbol{v}}_i^{\rm{H}}U{v_k}} \right)\left[ \begin{array}{l} \mathit{\boldsymbol{d}}{\left( {{u_k}} \right)^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_i}} \right)\mathit{\boldsymbol{d}}{\left( {{v_k}} \right)^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_i}} \right)\\ \mathit{\boldsymbol{d}}{\left( {{u_k}} \right)^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_i}} \right)\mathit{\boldsymbol{d}}{\left( {{v_k}} \right)^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_i}} \right) \end{array} \right]} \right\} \times \\ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{d}}{{\left( {{v_k}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_k}} \right)}&{\rho \left( k \right)}\\ {\rho \left( k \right)}&{\mathit{\boldsymbol{d}}{{\left( {{u_k}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_k}} \right)} \end{array}} \right] \end{array} $

(36)

where

$ \begin{array}{*{20}{c}} {w\left( i \right) = \mathit{\boldsymbol{d}}{{\left( {{u_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_i}} \right)\mathit{\boldsymbol{d}}{{\left( {{v_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_i}} \right) - }\\ {\frac{1}{4}{{\left[ {\mathit{\boldsymbol{d}}{{\left( {{u_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_i}} \right)\mathit{\boldsymbol{ + d}}{{\left( {{v_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_i}} \right)} \right]}^2},}\\ {\rho \left( i \right) = \frac{{ - \left[ {\mathit{\boldsymbol{d}}{{\left( {{u_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{v_i}} \right)\mathit{\boldsymbol{ + d}}{{\left( {{v_i}} \right)}^{\rm{H}}}\mathit{\Pi }_{B'}^ \bot \mathit{\boldsymbol{d}}\left( {{u_i}} \right)} \right]}}{2} \cdot }\\ {\mathit{\Pi }_{B'}^ \bot = {\mathit{\boldsymbol{I}}_{\left( {M + N - 1} \right)P}} - \mathit{\boldsymbol{B'}}{{\left( {{{\mathit{\boldsymbol{B'}}}^{\rm{H}}}\mathit{\boldsymbol{B'}}} \right)}^{ - 1}}{{B'}^{\rm{H}}}.} \end{array} $

According to Eq.(36), the variance of u_k estimate error is $ E\left[{\partial u_k^2} \right] = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{kk}}\left( {1, 1} \right)$, and that of v_k estimate error $E\left[{\partial v_k^2} \right] = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{kk}}\left( {2, 2} \right)$.

According to Eqs.(21) and (25), we get the MSE of the DOA and Doppler frequency

$ E\left[ {\partial \theta _k^2} \right] = \frac{1}{{{{\cos }^2}{\theta _k}}}E\left[ {\partial v_k^2} \right] = \frac{1}{{{{\cos }^2}{\theta _k}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{kk}}\left( {2,2} \right) $

(37)

$ E\left[ {\partial f_k^2} \right] = \frac{{f_{\rm{s}}^2}}{{4{{\rm{ \mathsf{ π} }}^2}}}E\left[ {\partial {u^2}} \right] = \frac{{f_{\rm{s}}^2}}{{4{{\rm{ \mathsf{ π} }}^2}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{kk}}\left( {1,1} \right) $

(38)

3.2 Cramer-Rao bound (CRB)

According to Ref. [29], the CRB for joint DOA and Doppler frequency estimation in a monostatic MIMO radar is given by

$ {\rm{CRB}} = \frac{{{\sigma ^2}}}{2}{\left[ {{\mathop{\rm Re}\nolimits} \left( {{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}^{\rm{H}}}\mathit{\Pi }_G^ \bot \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}} \right)} \right]^{ - 1}} $

(39)

where Δ=[Δ₁, Δ₂], ${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_1} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{d}}_1}{\mathit{\boldsymbol{s}}_1}\left( {{t_1}} \right)}& \cdots &{{\mathit{\boldsymbol{d}}_K}{\mathit{\boldsymbol{s}}_K}\left( {{t_1}} \right)} \\ \vdots&\ddots&\vdots \\ {{\mathit{\boldsymbol{d}}_1}{\mathit{\boldsymbol{s}}_1}\left( {{t_J}} \right)}& \cdots &{{\mathit{\boldsymbol{d}}_K}{\mathit{\boldsymbol{s}}_K}\left( {{t_J}} \right)} \end{array}} \right]$, ${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_2} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{a}}_1}{{\mathit{\boldsymbol{s'}}}_1}\left( {{t_1}} \right)}& \cdots &{{\mathit{\boldsymbol{a}}_K}{{\mathit{\boldsymbol{s'}}}_K}\left( {{t_1}} \right)} \\ \vdots&\ddots&\vdots \\ {{\mathit{\boldsymbol{a}}_1}{{\mathit{\boldsymbol{s'}}}_1}\left( {{t_J}} \right)}& \cdots &{{\mathit{\boldsymbol{a}}_K}{{\mathit{\boldsymbol{s'}}}_K}\left( {{t_J}} \right)} \end{array}} \right]$, $\mathit{\boldsymbol{G}} = \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}&{}&{\bf{0}} \\ \vdots&\ddots&\vdots \\ {\bf{0}}&{}&\mathit{\boldsymbol{A}} \end{array}} \right]{d_k} = \partial {\mathit{\boldsymbol{a}}_k}/\partial {\theta _k}, {\mathit{\boldsymbol{s}}_k}\left( t \right) = \partial {\mathit{\boldsymbol{s}}_k}\left( t \right)/\partial {f_k}$, $\varPi _G^ \bot = \mathit{\boldsymbol{I}}-\mathit{\boldsymbol{G}}{\left( {{\mathit{\boldsymbol{G}}^{\text{H}}}\mathit{\boldsymbol{G}}} \right)^{-1}}{\mathit{\boldsymbol{G}}^{\text{H}}}$.

4 Simulation Results

In this section, Monte Carlo simulations are conducted to measure the angle and frequency estimation performance of the proposed algorithm. In the following simulations except for Fig. 3, the number of Monte Carlo trials is 500. Define root mean square error (RMSE) of DOA/Doppler frequency as

$ {\rm{RMSE}} = \frac{1}{K}\sum\limits_{k = 1}^K {\sqrt {\frac{1}{{500}}\sum\limits_{l = 1}^{500} {{{\left( {{{\hat a}_{k,l}} - {a_k}} \right)}^2}} } } $

(40)

where ${\hat a_{k, l}}$ is the estimate of DOA/Doppler frequency a_k at the lth Monte Carlo trial (l=1, …, 500).

It is assumed that there are K = 3 incoherent targets located at θ₁=10°, θ₂=20°, θ₃=30°, and the Doppler frequencies are f₁=1000Hz, f₂=2000Hz, f₃=3000Hz, separately. f_s=10kHz is considered in the simulations. M, N, J, P and K are the number of transmit antennas, the receive antennas, the snapshots, the number of sections and the number of targets, respectively.

Fig. 3 depicts the joint angle and frequency estimation of our CE-MUSIC algorithm in SNR=0dB with M=8, N=8, J=100 and P=2. The number of Monte Carlo simulations is 50. It is illustrated that CE-MUSIC algorithm is effective even in low SNR.

Fig. 4 shows the joint DOA and Doppler frequency estimation performance comparison of the proposed CE-MUSIC with 2D-MUSIC, ESPRIT and CRB (assuming M=8, N=8, P=6 and J=100). It is illustrated that the CE-MUSIC algorithm has very close performance to the conventional 2D-MUSIC and both of them have better performance than ESPRIT. Furthermore, CE-MUSIC has much lower complexity than 2D-MUSIC.

Fig. 5 shows DOA and frequency estimation performance of our algorithm with M=6, N=6, J=50 and different values of P. It is clearly shown in Fig. 5 that when P increases, the DOA estimation performance changes very little. The frequency estimation performance gets improved obviously when P is increased from 3 to 9 and the improvement becomes very minor when P is bigger than 9. The increase of P results in better frequency estimation performance and the limit of this improvement exists. Meanwhile it has minor improving effects on DOA estimation performance.

Fig. 6 shows angle and frequency estimation performance of CE-MUSIC algorithm with M=6, N=6, P=3 and different snapshots. From Fig. 6, the DOA and frequency performance of CE-MUSIC algorithm is improved significantly when the number of snapshots increases. More snapshots means more samples, which can make the algorithm perform better.

Figs. 7—8 present angle and frequency estimation performance of CE-MUSIC algorithm with J=100, P=3 and different M/N. It is shown that the DOA and frequency estimation performance of the proposed algorithm get improved with the number of transmit/receive antennas increasing. Multiple transmit/receive antennas improve angle estimation performance because of the diversity gain.

5 Conclusions

In this paper, we have proposed a computationally efficient MUSIC-based algorithm for joint DOA and Doppler frequency estimation in a monostatic MIMO radar. The proposed algorithm has much lower computationally complexity than conventional 2D-MUSIC algorithm. The simulation results prove that our CE-MUSIC algorithm outperforms ESPRIT algorithm and has very close performance to conventional MUSIC algorithm. Furthermore the proposed algorithm requires no pair matching.

Acknowledgements

This work was supported in part by the Funding for Outstanding Doctoral Dissertation in NUAA (No.BCXJ15-03), the Funding of Jiangsu Innovation Program for Graduate Education (No.KYLX15_0281), and the Fundamental Research Funds for the Central Universities.