Computationally efficient MUSIC based DOA estimation algorithm for FMCW radar

Nowadays, radar sensors are receiving more and more attention from research communities as well as industry, and playing an essential role in various applications including defense and surveillance. Further, radar sensors are recently integrated with automotive vehicles due to the robustness of these sensors under diverse environmental conditions, such as fog, strong/no light, and humidity [1−4]. To meet various requirements, different modulation schemes are designed for automotive radar sensors [4,5]. Among these, the frequency-modulated continuous wave (FMCW) method is the most popular solution [6,7]. That is to say, FMCW radars are broadly employed among other radar systems because of their unique features, such as the small size, low power consumption, and low costs [8−11]. Although there have been great advancements in system design and hardware integration in the multiple input and multiple output (MIMO) automotive FMCW radar, high-resolution target detection along with low-cost signal processing remains as a challenging problem. In the case of massive MIMO radars, high-number channels are usually required, which results in high computational costs and eventually introduces the intolerable latency [12]. Thus, the real-time deployment of the (massive) MIMO automotive radar calls for high-accuracy yet low-cost signal processing techniques.

Various methodologies are presented in the literature to accurately estimate the directions of the intended sources. The minimum variance distortionless response (MVDR) [13] algorithm relies on the inversion of the measured covariance matrix (CM). However, the performance of this method degrades significantly, when CM is estimated over a limited number of snapshots and/or the incident signals are dependent. Additionally, the MVDR algorithm is computationally expensive due to the computation of inverse CM. Later, subspace-based algorithms were proposed to enhance the resolution of the MVDR technique. The multiple signal classification (MUSIC) [14] and estimation of signal parameters via rational invariance techniques (ESPRIT) [15] are the common representatives of supper resolution subspace techniques. By using the exact subspaces extracted from CM through eigenvalue decomposition (EVD), these algorithms achieve near-optimal performance [16]. Hence, they are widely employed in FMCW radar systems [17−20]. However, the high complexity and intolerable latency caused by subspace extraction limit their practical use [21]. The minimum norm (MinNorm) technique proposed by Reddi [22] and later enhanced by Kumaresan and Tufts [23] is another subspace-based algorithm, which tries to minimize the array output norm via the best array weights vector. However, the applicability of the MinNorm method is limited to the uniform linear array (ULA). Importantly, all the above-mentioned algorithms (i.e., MUSIC, ESPRIT, and MinNorm) share a common problem which is the high cost in the subspace acquisition stage. That is to say, these algorithms require $ \mathcal{O}\left({M}^{3}\right) $ time complexity, where $ M $ is the number of array elements, for either taking the inversion of CM (as in the case of MVDR) or decomposition of CM (in the cases of MUSIC and MinNorm). Thus, the high latency resulted from their computational burden makes these algorithms impractical [21], particularly when they are applied to massive MIMO automotive radar systems.

Simplifying the cubic time complexity $ \mathcal{O}\left({M}^{3}\right) $ has become the primary objective of several works. The fast Fourier transform (FFT) is widely employed as an alternative solution in (automotive) FMCW radar systems to detect the range, velocity, and directions of the targets, due to its promising low complexity, i.e., $ \mathcal{O}(M{\mathrm{log}}_{2}M) $ [20,24−26]. However, the resolution of the FFT method is not sufficient for detecting closely-spaced sources [27]. As one popular solution, the complexity of the subspace computation can also be reduced through a beamforming weight vector adapted to the beamspace method [28] which transforms the element space into the beamspace. Despite its favorable performance in the uniform circular array (UCA) [29], the beamspace method may suffer from severe performance degradation in the case of ULA. To avoid the eigen-decomposition of CM, the linear direction of arrivals (DOA) estimation algorithms, such as the propagator [30] and orthogonal propagator [31] schemes are introduced, whereby the noise subspace is estimated via the propagator method. The concept of the propagator method is built based on partitioning the steering vectors of the received signal into two submatrices to identify the subspaces. Comparing the eigenfeature based methods, this technique enjoys the complexity reduction from $ \mathcal{O}\left({M}^{3}\right) $ to $ \mathcal{O}\left({M}^{2}L\right) $ [32], where $ L $ is the number of targets to be detected. This simplification, however, comes at the cost of limited estimation capability, particularly when the targets are located within 70° to 90° with respect to the radar [33]. Furthermore, the complexity of the subspace computation can be simplified by processing a small set of columns that retain the overall trends of the original large matrix, known as columns subset selection (CSS) [34]. From this perspective, uniform sampling (US) [35] and non-uniform sampling (NUS) [36] methods are recently introduced to formulate the projection matrix (PM) based on a certain number of CM columns for localizing narrowband sources. The corresponding DOA algorithms are called the PM-US and PM-NUS methods, respectively. Both PM-US and PM-NUS algorithms guarantee better estimation accuracy than the propagator method, where they have the same complexity-simplification in the subspace acquisition step, i.e., $ \mathcal{O}\left({M}^{2}L\right) $.

In this paper, we aim to break down the complexity of the subspace computation to lower than that reached by the above techniques, i.e., $ \mathcal{O}\left({M}^{2}L\right) $, while maintaining the high-estimation accuracy. That is to say, the objective of this work is to realize real-time target detection yet with resolution comparable to the standard MUSIC algorithm. To this end, we apply a simple randomized low-rank approximation to CM $ \mathbf{R} $ and approximately compute its signal/noise subspace which is the major bottleneck in subspace-based algorithms. Specifically, we firstly approximate large $ \mathbf{R}\in {\mathbb{C}}^{M\times M} $ through three small sketches, in the form of $ \mathbf{R}\approx \mathbf{Q}\mathbf{B}{\mathbf{Q}}^{\mathrm{H}} $, where $ {\mathbb{C}}^{M\times M} $ denotes an $ M\times M $ complex matrix. Here, the matrix $\mathbf{Q}\in {\mathbb{C}}^{M\times z}$ contains the orthonormal basis for the range of the sketch matrix $ \mathbf{C}\in {\mathbb{C}}^{M\times z} $ which is extracted from $ \mathbf{R} $ by using the randomized uniform sampling technique (i.e., every column has a unique probability to be selected), and $ \mathbf{B}\in {\mathbb{C}}^{z\times z} $ is a weight-matrix reducing the approximation error. Relying on such randomized approximation of $ \mathbf{R} $, computations of subspace extraction can be performed more rapidly. That is, the proposed efficient MUSIC (E-MUSIC) algorithm accomplishes high-resolution target detection at the linear complexity $ \mathcal{O}\left({z}^{2}M\right) $, whereby $ z $ is a user-specified parameter which tradeoffs between the complexity and angular accuracy, which satisfies $ L\le z\ll M $. Moreover, we establish an error bound for the devised algorithm showing that the unitary matrix $ \mathbf{Q} $ accurately captures most of the actions of $ \mathbf{R} $. We comprehensively evaluate the performance of our E-MUSIC algorithm through numerical simulations and systematically compare it with the counterpart subspace-based algorithms. The achieved results validate the accuracy and efficiency of the proposed algorithm.

The rest of this paper is organized as follows. In Section 2, we explain the system model and signal generation of the FMCW radar. Section 3 briefly introduces the principle of subspace DOA estimation. The proposed E-MUSIC algorithm along with its complexity analysis and error bound is explained in Section 4. The obtained simulation results are presented and comprehensively discussed in Section 5. Finally, Section 6 concludes this paper.

In this paper, we use an automotive FMCW radar system which is composed of a single transmitting antenna ($ {T}_{x} $) and $ M $ receiving antennas ($ {R}_{x} $) as shown in the Fig. 1. The emitted FMCW radar signal from the $ {T}_{x} $ antenna is represented as follows:

Figure 1.  System model of the considered automative FMCW radar using single $ {T}_{x} $ and $M\;{R}_{x}$ antennas.

where $ {f}_{c} $ is the carrier frequency of the radar system, $ t $ is the time range, and $ \mu $ is a linear slope (i.e., the chirp slope) that rises according to the time $ \tau $ within the single chirp duration $ T $. Hence, $ \mu $ is defined as $ \mathrm{B}\mathrm{W}/T $, where $ \mathrm{B}\mathrm{W} $ represents the analog bandwidth of the radar.

The echo-signal reflected from the $ l $th target and received by the $ m $th $ {R}_{x} $ antenna is represented by

where $ {\rho }_{l} $ denotes the complex amplitude of the echo-signal reflected by the $ l $th target, $ {\tau }_{l} $ is a latency introduced to the signal due to the range between the $ l $th target and the radar, $ \lambda $ is the signal wavelength, $ d=\lambda /2 $ is the spatial distance between adjacent array elements, $ {\theta }_{l} $ represents DOA of the $ l $th target, and ${{z}}_{m}\left({t}\right)$ denotes the additive white Gaussian noise (AWGN) at the $ m $th $ {R}_{x} $ antenna.

One of the benefits of the FMCW radar system is that its signal can be de-modulated/de-chirped through a simple mixer, which significantly reduces the RF implementation costs [20]. Specifically, at the reception side, the $ {R}_{x} $ signals are multiplied with the conjugate of the signal of $ {T}_{x} $ to extract the so-called beat signal $ {y}_{m}\left(t\right) $ as follows:

After substituting (1) into (3) and using a low pass filter having an impulse response $ h\left(t\right) $, the target signal can be obtained and represented as follows [27]:

where ${\widetilde{{z}}}_{m}\left({t}\right)=h\left(t\right)\otimes \left[{{z}}_{m}\left({t}\right){s}_{{T}_{x}}\left(t\right)\right]$ denotes the noise after the de-chirping process and $ \otimes $ is the convolution operator. Then, an analog-to-digital convertor (ADC) with a sampling frequency $ {f}_{s}=1/{t}_{s} $ ($ {t}_{s} $ is the sampling-time interval) is applied to convert the analog signal $ {y}_{m}\left(t\right) $ to the digital form as

where $ 0\le m\le M-1 $ and $ 0\le n\le N-1 $ with $\left\lfloor {N\triangleq T/{t}_{s}} \right\rfloor$ being the total number of data samples, where $\left\lfloor {\cdot} \right\rfloor$ denotes the floor operator. Thus, the directions of the targets can be estimated from the digital beat signal $ {y}_{m}\left(n\right) $.

In this work, we do not focus on the channel environment, but the assumption of a line of side (LOS) channel is considered. The main focus will be on how to simplify subspace estimation. As the automotive radar ranges are limited to an order of hundred meters, the atmospheric effect is also marginal.

To achieve high resolution, subspace-based DOA algorithms are usually applied to the FMCW radar system [20,37]. MUSIC [14] and ESPRIT [15] algorithms are the representatives of the subspace methods [38]; all start from estimating CM from the signal matrix $ \mathbf{Y}\in {\mathbb{C}}^{M\times N} $, where

and CM $ \mathbf{R} $ is estimated by

Therefore, $ \mathbf{R}\in {\mathbb{C}}^{M\times M} $ is represented as follows:

It is noticeable that matrix-multiplication in (5) to estimate $ \mathbf{R} $ causes high computational complexity. Fortunately, this multiplication can be realized in parallel more efficiently [39]. We assume that there are $ L $ targets and the large CM $ \mathbf{R} $ is decomposed using singular value decomposition (SVD) [40] as follows:

Here $ \mathbf{\Sigma }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\left\{{\sigma }_{i}\right\}\right)$ for $i=1{\rm{,}}\;2{\rm{,}}\;\cdots{\rm{,}}\;M$, where $ {\sigma }_{i} $ is the $ i $th singular value of $ \mathbf{R} $ and ${\sigma }_{1} > {\sigma }_{2} > \cdots > \sigma _{L} > {\sigma }_{L+1}={\sigma }_{L+2}=\cdots ={\sigma }_{M}$. The columns of $ \mathbf{U} $ and $ \mathbf{V} $ represent the left and right singular vectors of $ \mathbf{R} $, respectively. By separating the signal and noise subspaces, (7) can be re-written as

where $ {\mathbf{U}}_{s} $ and $ {\mathbf{U}}_{n} $ are the signal and noise subspaces, respectively; $ {\mathbf{V}}_{s} $ and $ {\mathbf{V}}_{n} $ are the matrices containing the right singular vectors of $ \mathbf{R} $ associated with signal and noise subspaces, respectively. That is, ${\mathbf{U}}_{s}=[{\mathbf{u}}_{1}{\rm{,}}\;{\mathbf{u}}_{2}{\rm{,}}\;\boldsymbol{ }\cdots {\rm{,}}\;{\mathbf{u}}_{L}]$ and ${\mathbf{U}}_{n}=[{\mathbf{u}}_{L+1}{\rm{,}}\;{\mathbf{u}}_{L+2}{\rm{,}}\;\boldsymbol{ }\cdots{\rm{,}}\;{\mathbf{u}}_{M}]$, whereby $ {\mathbf{u}}_{i} $ represents the $ i $th singular vector of $ \mathbf{R} $ and $ {\sigma }_{n}^{2} $ denotes the noise power.

Next, we divide the elevation-angle range $\left[-\pi /2{\rm{,}}\;\pi /2\right]$ into a grid of $ L $ angles $\left\{{\theta }_{0}{\rm{,}}\;{\theta }_{1}{\rm{,}}\;\cdots {\rm{,}}\;{\theta }_{L-1}\right\}$. For each angle $ {\theta }_{l} $, we define an $ M\times 1 $ steering vector for the linear array as follows:

whereby

Now, we can use the computed subspaces to generate the pseudo-spectrum for the standard MUSIC algorithm as follows:

where $ {\mathbf{I}}_{M} $ is an $ M\times M $ identity matrix. The estimated targets’ directions are indicated by the peaks of pseudo- spectrum $ {\mathbf{F}}_{\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}\left(\theta \right) $. Thus, the unknown directions of several targets, ${\theta }_{l}\,(l=1{\rm{,}}\;2{\rm{,}}\;\cdots {\rm{,}}\;L)$, are estimated by

Due to the need to search through all the grids (i.e., computing $ {\boldsymbol{F}}_{\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}\left(\theta \right) $ for every $\theta \in \left\{{\theta }_{0}{\rm{,}}\;{\theta }_{1}{\rm{,}}\;\cdots {\rm{,}}\;{\theta }_{L-1}\right\}$), the time complexity of pseudo-spectrum generation is high, $ \mathcal{O}\left(M(M-L\right)) $ [35]. The evaluation of $ {\mathbf{F}}_{\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}\left(\theta \right) $ is not the focus of this paper, because it can be computed in an efficient manner [41]. Traditionally, acquiring $ {\mathbf{U}}_{s} $ needs $ \mathcal{O}\left({M}^{3}\right) $ time complexity due to involving the decomposition of large CM. It is worth mentioning that none of the existing parallel algorithms can reduce this complexity. Therefore, the emphasis of this work is on the fast computation of the signal subspace $ {\mathbf{U}}_{s} $, which is the primary bottleneck in subspace-based methods. To be specific, in the following section, we introduce a method to reduce the time complexity from $ \mathcal{O}\left({M}^{3}\right) $ to $ \mathcal{O}\left({z}^{2}M\right) $, where $ z\ll M $, by using an approximated technique.

In this section, the proposed E-MUSIC algorithm is presented and its computational complexity along with its error bound is analyzed.

The proposed method approximates large CM $ \mathbf{R} \in {\mathbb{C}}^{M\times M} $ with three small matrices in the form of $ \mathbf{R}\approx \mathbf{Q}\mathbf{B}{\mathbf{Q}}^{\mathrm{H}} $, where $ \mathbf{Q} $ is the orthonormal basis of the sketch $ \mathbf{C} $ shown in Fig. 2 and $ \mathbf{B}\in {\mathbb{C}}^{z\times z} $ is a weight matrix reducing the approximation error. Then, we compute approximated $ z $-SVD to the small sketches rather than exact SVD to original CM. Our method is inspired by the randomized low-rank matrix approximation studied in linear algebra and scientific computing [42−44]. To the best of authors’ knowledge, this is a unique attempt made to solve a long-standing complexity problem in the subspace-based DOA estimation through such approximation. We partition the algorithm procedure into three stages.

Figure 2.  Illustration of randomized low-rank matrix approximation of large $ \mathbf{R}\in {\mathbb{C}}^{M\times M} $ with the rank of $ \mathbf{R} $ being 3. Here $ z=3 $ and the indexing set $\mathcal{I}=\left\{1{\rm{,}}\;3{\rm{,}}\;6\right\}$. Thus, $ \mathbf{Q} $ has 10 rows with 3 columns, $ {\mathbf{Q}}^{\mathrm{H}} $ has 3 rows with 10 columns, and the square matrix $ \mathbf{B} $ has 3 columns and 3 rows.

To reduce the computational cost of subspaces acquisition, we leverage a simple random sampling method to obtain a matrix sketch $ \mathbf{C}\in {\mathbb{C}}^{M\times z} $ from R. Specifically, we select the sampling parameter $ z $ that satisfies $ L\le z\ll M $. Then, z items uniformly at random are selected within $\left\{0{\rm{,}}\;1{\rm{,}}\;\cdots {\rm{,}}\;M-1\right\}$ to form the indexing set $ \mathcal{I} $, i.e., $ \left|\mathcal{I}\right|=z $. Thus, the sketching matrix is formulated as $\mathbf{C}=\mathbf{R}\left(:{\rm{,}}\;\mathcal{I}\right)=\mathbf{R}{\Cambriabfont\text{Ω}}\in {\mathbb{C}}^{M\times z}$, where ${\Cambriabfont\text{Ω}}$ is a uniform sampling matrix having only one non-zero entry at each row and the location of this entry is selected randomly. Such that, $\mathbf{Y}={\mathbf{X}}_{N\times M}{\Cambriabfont\text{Ω}}_{M\times z}$ attains $ z $ uniformly sampled columns of $ \mathbf{X} $. Due to the randomness, the sampling matrix ${\Cambriabfont\text{Ω}}$ forms a linearly independent set and thus the sampled columns become linearly independent, which allows well-separation between noise and signal subspaces. As a result, the sampled columns span the range of original CM. Therefore, to compute the orthonormal basis for the range of $ \mathbf{R} $, we just need to orthonormalize the sampled matrix $ \mathbf{C} $. That is to say, the sampling matrix ${\Cambriabfont\text{Ω}}$ randomly selects $ z $ columns from $ \mathbf{R} $ (without replacement) to construct $ \mathbf{C} $. Hence, the sketching matrix $ \mathbf{C} $ has $ z $ columns of $ \mathbf{R} $ corresponding to the indices in $ \mathcal{I} $. Next, we compute the orthonormal basis of the obtained sketch $ \mathbf{C} $ as $ \mathbf{Q}\leftarrow \mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\left(\mathbf{C}\right)\in {\mathbb{C}}^{M\times z} $. Alternatively, one can obtain the orthonormal basis by computing the QR decomposition of $ \mathbf{C} $, so that $\left(\mathbf{Q}{\rm{,}}\; \sim \right)\leftarrow \mathrm{q}\mathrm{r}\left(\mathbf{C}\right)$. Now, the range of the orthonormal matrix Q captures most of the actions of $ \mathbf{R} $ (this claim will be proven in subsection 4.3). It can be noted that matrices C, ${\Cambriabfont\text{Ω}}$, and Q contain sufficient information to approximate $ \mathbf{R} $.

To see the intuition behind the proposed method, let us define the (currently unknown) matrix $ \mathbf{B}={\mathbf{Q}}^{\mathrm{H}}\mathbf{R}\mathbf{Q} $. Post multiplying B by ${\mathbf{Q}}^{\mathrm{H}}{\Cambriabfont\text{Ω}}$ results the following identity:

As $ \bf{Q} $ has orthonormal columns (i.e., $ \mathbf{R}\mathbf{Q}{\mathbf{Q}}^{\mathrm{H}}\approx \mathbf{R} $) and $\mathbf{R}{\Cambriabfont\text{Ω}}=\mathbf{C}$, the unknown matrix $\bf{B}$ must satisfy

All the matrices ${\Cambriabfont\text{Ω}}$, $ \mathbf{Q} $, and $ \mathbf{C} $ exist, therefore, (12) can be easily solved to obtain the unknown matrix $ \mathbf{B} $. Then, we simply achieve the following low-rank factorization of $ \mathbf{R} $:

Referring to the fact that the spectrum and Frobenius norms are unitary invariant [45] (i.e., ${\left\|\mathbf{X}\mathbf{W}{\mathbf{X}}^{\mathrm{H}}\right\|}_{{\epsilon}}={\left\|\mathbf{W}\right\|}_{{\epsilon}}$, for ${\epsilon}\in \{2{\rm{,}}\;F\}$, where $ \mathbf{X} $ is a unitary matrix and $ \mathbf{W} $ is a matrix with a suitable dimension) and relying on the low-rank factorization in (13), one can approximately identify the subspaces of $ \mathbf{R} $ by computing SVD of the sketch matrix $ \mathbf{B} $. Suppose that $\mathbf{B}={\mathbf{U}}_{B}{\Cambriabfont\text{Σ}}_{B}{\mathbf{V}}_{B}^{\mathrm{H}}$, and we achieve the rank-restricted approximate SVD of $ \mathbf{R} $ as follows:

where $ {\mathbf{U}}_{B}\in {\mathbb{C}}^{z\times z} $ and ${\bf{V}}_{B}\in {\mathbb{C}}^{z\times z}$ contains left and right singular vectors of $ \mathbf{B} $ and ${\Cambriabfont\text{Σ}}_{B}\in {\mathbb{R}}^{z\times z}$ is a diagonal matrix containing nonnegative eigenvalues of $ \mathbf{B} $ sorted in the descending order.

One can further simplify and re-write rank-restricted SVD of $ \mathbf{R} $ in (14) as follows:

where $\tilde {\mathbf{U}}\triangleq \mathbf{Q}{\mathbf{U}}_{B}{\in \mathbb{C}}^{M\times z}$ and ${\tilde {\mathbf{V}}}^{\rm{H}}\triangleq {\bf{V}}_{{B}}^{\rm{H}}{\bf{Q}}^{\rm H}{\in \mathbb{C}}^{z\times M}$. Relying on the above procedure, we are now ready to estimate the approximated signal subspace by ${\tilde {\mathbf{U}}}_{s}=\tilde {\mathbf{U}}\left(:{\rm{,}}\;1:L\right)$. Finally, the approximated pseudo-spectrum can be computed as

Based on the above subspace approximation, the directions of multiple targets are estimated as

The main steps of the proposed method are summarized in Table 1.

Relying on the above analysis and the algorithms’ procedure presented in Table 1, the overall complexity of the proposed method is $ \mathcal{O}\left({z}^{2}M\right) $ [42], which comes from the matrix multiplication in Stage A. We will analyze it in detail as follows.

Uniform sampling constructs a matrix sketch $ \mathbf{C} $ without visiting the whole matrix [46]. Therefore, randomized uniform sampling performed in Step 1 has constant complexity, $\mathcal{O}(1) $. In addition to its low complexity, it is shown that such uniform random sampling is more effective in terms of increasing the degree of freedom (DOF) with the same array size, removing the dependency between selected columns, and eventually making the singular values more distinguishable [35]. The superiority of this sampling method (i.e., randomized uniform sampling) has also been proven over other sampling techniques in Ref. [47].

Computing the orthonormal matrix $\bf{Q}$ in Step 2, which involves matrix multiplication between ${\mathbf{C}}^{\mathrm{H}} \in {\mathbb{C}}^{z\times M}$and $ \mathbf{C}\in {\mathbb{C}}^{M\times z} $, needs $ \mathcal{O}\left({z}^{2}M\right) $ time complexity. The complexity of computing the matrix $ \mathbf{B} $ in Step 3 is identical to that of Step 2.

Decomposition of the matrix $ \mathbf{B}\in {\mathbb{C}}^{z\times z} $ in Step 1 requires $ \mathcal{O}\left({z}^{3}\right) $ times of flops. Combined all together, the overall complexity of the proposed technique to estimate the approximated signal subspace ${\tilde{\mathbf{U}}}_{s}$ is $ \mathcal{O}\left({z}^{2}M\right) $. This is significantly lower than applying exact SVD to large CM $ \mathbf{R} $, which needs $ \mathcal{O}\left({M}^{3}\right) $ arithmetic operations. Note that, dealing with the computation in Stage C is not the focus of this paper as it can be computed more efficiently [41].

From the above analysis, it seems that the complexity of our method becomes identical to that of the standard MUSIC algorithm for $ z=M $. However, it is worth mentioning that the actual complexity of our algorithm is noticeably lower than $ \mathcal{O}\left({z}^{2}M\right) $. Because the polynomial complexity on $ z $ comes from the matrix multiplication. Fortunately, in many cases, the matrix multiplication can be implemented in parallel. Therefore, its complexity is dramatically lower than SVD, matrix inversion, or QR decomposition even when $ z=M $. For example, the Strassen-like algorithms [48] reduce the matrix multiplication cost between the matrices $ \mathbf{A}\in {\mathbb{F}}^{z\times z} $ and $ \mathbf{B}\in {\mathbb{F}}^{z\times z} $ from $ \mathcal{O}\left({z}^{3}\right) $ to $ \mathcal{O}\left({z}^{2.376}\right) $, $\mathbb{F}\in \left\{\mathbb{C}{\rm{,}}\;\mathbb{R}\right\}$. From this perspective, our method is more efficient than the presented algorithms in Ref. [19] and Ref. [49], as their costs are dominated by the matrix decomposition rather than matrix multiplication. The subspace computational costs in several algorithms along with the proposed method are listed in Table 2.

As the signal processing cost of subspace acquisition is dominated by the array size $ M $, we conduct a numerical example to calculate the number of operations required for subspace acquisition based on various array sizes for different algorithms as shown in Fig. 3.

Figure 3.  Number of arithmetic operations needed for subspace acquisition vs. array elements (M), where $ L=5 $ and $ z=L+2 $.

Fig. 3 confirms that the subspace extraction based on the proposed method requires the minimum number of operations (i.e., its complexity is linearly scalable with $ M $) compared with other well-known algorithms. Particularly, a significant gap is observed between the cost of applying exact SVD to CM adopted in standard MUSIC and the approximated SVD used in the proposed method.

In this subsection, we establish the error bound for the proposed algorithm presented in Table 1.

Notations. Let $\mathbf{U}{\Cambriabfont\text{Σ}}{\mathbf{V}}^{\mathrm{H}}$ be SVD for observed CM $ \mathbf{R}\in {\mathbb{C}}^{M\times M} $. To conduct the analysis, we partition SVD as follows:

We then factorize the random sampling matrix, defined in subsection 4.1 into two parts as

The above decomposition allows us to re-formulate the sketching matrix C as

Intuitively, ${\Cambriabfont\text{Σ}}_{2}{{\Cambriabfont\text{Ω}}}_{2}$ is a perturbation of $ \mathbf{R} $. As shown in Algorithm 1, the proposed method forms the orthonormal basis $ \mathbf{Q} $ for the range of the sketching matrix $ \mathbf{C} $. Our objective here is to quantify how well the range of the matrix $ \mathbf{Q} $ captures the action of $ \mathbf{R} $. Thus, the problem at hand is to achieve the bounds on the following approximation error:

where $ {\mathbf{P}}_{\mathbf{C}} $ is an orthogonal projector (OP) that satisfies $ {\mathbf{P}}_{\mathbf{C}}={\mathbf{P}}_{\mathbf{C}}^{2} $. For a given full column-rank matrix ${\Cambriabfont\text{Ξ}}$, we can express $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left({\mathbf{P}}_{{\Cambriabfont\text{Ξ}}}\right)=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left({\Cambriabfont\text{Ξ}}\right)$ where $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(\cdot)$ is the range of the matrix, and ${\mathbf{P}}_{{\Cambriabfont\text{Ξ}}}$ is explicitly formulated by

Particular OP is determined by its range. Since $ \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathbf{Q}\right)=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathbf{C}\right) $, we are allowed to express $ {\mathbf{P}}_{\mathbf{C}} $ as

Theorem 1. Assuming the notations defined above and that ${\Cambriabfont\text{Ω}}_{1}$ attains the full row-rank property, the following approximation error is guaranteed:

where $ {\sigma }_{L+1} $ denotes the $ (L+1) $th singular value of $ \mathbf{R} $.

The accuracy of subspaces obtained from Algorithm 1 is dominated by the matrix $ \mathbf{Q} $, therefore, this bound is a reasonable approximation error for the generated pseudo-spectrum of the proposed method. In the high signal noise ratio (SNR) situation, $ {\sigma }_{L+1} $ approaches zero. Consequently, the approximation error becomes marginal. This directly implies that ${\tilde{\mathbf{F}}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{d}}\left(\theta \right)$ well approximates and closely tracks $ {\mathbf{F}}_{\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}\left(\theta \right) $. The proof of Theorem 1 can be found in Appendix A.

It is intuitive to compare (22) with that of the traditional MUSIC algorithm. As the MUSIC algorithm uses exact CM (i.e., all the columns of CM) and thus computes exact (rather than approximated) SVD to extract subspaces, one can model its approximation error as follows:

where $\widehat{\mathbf{R}}=\mathbf{U}{\Cambriabfont\text{Σ}}{\mathbf{V}}^{\mathrm{H}}$. Therefore, it is observable that the approximation error of the standard MUSIC algorithm seems to be lower than that for the proposed method due to adopting the best rank-$ L $ approximation. i.e.,

However, the optimal rank-$ L $ approximation of $ \mathbf{R} $ used in the MUSIC algorithm requires the cubic running time of $ \mathcal{O}\left({M}^{3}\right) $ as previously shown in subsection 4.2.

In this section, the performance of the proposed E-MUSIC algorithm is evaluated and systematically compared with the well-known subspace-based algorithms in the literature. Additionally, we include the MVDR algorithm in our comparison because it is widely used in radar signal processing [50−52]. As a measure of accuracy, we compute the average root mean square error (ARMSE) and probability of successful detection (PSD) for all the algorithms as follows:

where $ K $ represents the number of Monte Carlo iterations; $ {\theta }_{k} $ and $ {\widehat{\theta }}_{k} $ are true and estimated DOA, respectively; $ {\mathcal{S}}_{i} $ denotes the number of successful detection at the $ i $th iteration.

To meet the requirements of the long-range radar (LRR), we consider the following parameters during FMCW signal generation as shown in Table 3. The maximum range of the radar is selected as 250 m and the chirp duration (T) is 5 times the corresponding time of the radar range. The carrier frequency $ {f}_{c} $ is set to be $77\;\mathrm{G}\mathrm{H}\mathrm{z}$ which is sufficient for millimeter-wave MIMO radar systems. We configure the signal BW to be $ 0.15\;\mathrm{G}\mathrm{H}\mathrm{z} $ which is satisfactory for the 1-m range resolution. We also take the target speed into account, such that the maximum speed of the target is 200 km/h and the maximum possible Doppler frequency is calculated accordingly. The signal is sampled at the rate equaling to the signal BW, $ {f}_{s}=0.15\;\mathrm{G}\mathrm{H}\mathrm{z} $, and the maximum received intermediate frequency (IF) is $ 30\;\mathrm{M}\mathrm{H}\mathrm{z} $. We use the runtime of the central processing unit (CPU) (CPU $ 2.20\;\mathrm{G}\mathrm{H}\mathrm{z} $, 4-GB Installed random access memory (RAM)), which corresponds to the needed number of multiplication flops, to calculate the time complexity of the investigated algorithms.

The objective of this section is twofold: First, to evaluate the effect of the number of selected columns $ z $ on the estimation performance of the proposed algorithm; second, to determine how many computations are saved due to applying the proposed methodology compared with the conventional method. Toward that end, we first compare the (approximated) pseudo-spectrum of the proposed algorithm based on different values of $ z $ with that of the MUSIC algorithm. The simulation parameters are set as follows: $ M=16 $, $ L=4 $, $\mathrm{S}\mathrm{N}\mathrm{R}=0\;\mathrm{d}\mathrm{B}$, $ N=100 $, and we investigate a various number of sampled columns as $z=4{\rm{,}}\;6{\rm{,}}\;\mathrm{a}\mathrm{n}\mathrm{d}\;8$ as shown in Fig. 4.

Figure 4.  Pseudo-spectra of standard MUSIC and proposed algorithms based on different numbers of sampled columns $z=4{\rm{,}}\;6{\rm{,}}\;\mathrm{a}\mathrm{n}\mathrm{d}\;8$ ($ M=16 $, $ L=4 $, $\mathrm{S}\mathrm{N}\mathrm{R}=0\;\mathrm{ }\mathrm{d}\mathrm{B}$, and $ N=100) $.

The figure shows that, under different values of $ z $ (i.e., z = L, z = L + 2, and z = L + 4), the proposed algorithm detects DOAs of all the available targets. The approximation error relative to the pseudo-spectrum of the standard MUSIC algorithm can be defined as $\sum _{l=1}^{L}{\left[{\widetilde {\mathbf{F}}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{d}}\left({\theta }_{l}\right)-{\mathbf{F}}_{\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}\left({\theta }_{l}\right)\right]}^{2}$. Expectably, the approximation error drops by increasing $ z $, whereby it converges to zero when z = 8.

Using the same simulation parameters, we then compute ARMSE and PSD of the MUSIC algorithm and the proposed method over a wide range of $ \mathrm{S}\mathrm{N}\mathrm{R} $s for different numbers of sampled columns $ z $. $ \mathrm{S}\mathrm{N}\mathrm{R} $ is varied from $ -10\mathrm{ }\mathrm{d}\mathrm{B} $ to $ 10\mathrm{ }\mathrm{d}\mathrm{B} $ with the step of $ 5\mathrm{ }\mathrm{d}\mathrm{B} $. To guarantee fair results, we apply one-thousand Monte Carlo simulations ($ K=1000 $), and the calculated estimation errors and successful detection for both algorithms are shown in Fig. 5 and Fig. 6.

Figure 5.  Impact of the number of sampled columns on estimation errors for $ M=16 $.

Figure 6.  Impact of the number of sampled columns on successful detection for $ M=16 $.

Unsurprisingly, the figures show that the estimation performance of the proposed technique is improved by including a larger number of sampled columns $ z $. As demonstrated, the suggested algorithm achieves near-optimal performance, when the number of the processed columns is half (i.e., $ z=8 $), compared with the standard MUSIC algorithm. In other words, the standard MUSIC algorithm needs to process $ 16\times 16 $ CM to achieve such a result, whereas the proposed method obtained the same accuracy by only processing an $ 16\times 8 $ matrix sketch, $ \mathbf{C} $. This reduces the overall complexity of the algorithm. To demonstrate that, we secondly compute the complexity of the proposed algorithm with different values of z as shown in Fig. 7.

Figure 7.  Complexity vs. the number of sampled columns, $ M=16 $.

It is worth mentioning that when $ z=16=M $, the complexity of our method becomes identical to that of the standard MUSIC algorithm. However, the suggested algorithm obtains the near-optimal performance with $ z=8 $ as confirmed by the previous results (i.e., Fig. 4, Fig. 5, and Fig. 6). Hence, Fig. 7 demonstrates that the proposed method saves 75% computation compared with the standard MUSIC algorithm without performance degradation. This is accomplished by using the randomized low-rank approximation of CM, which accelerates the subspace computation by the orders of the magnitude. To conclude, the presented results in this subsection confirm that the devised technique breaks the contradictory relation between complexity and estimation accuracy, so that the overall complexity is reduced without sacrificing angular accuracy.

In this section, we compare the runtime of the proposed method with popular subspace algorithms. Here, we only compare the runtime for computing signal subspace which dominates the overall cost of subspace-based DOA estimation techniques. In this simulation, we assume L = 2, z = L + 2, N = 100, and $ \mathrm{S}\mathrm{N}\mathrm{R}=0\mathrm{ }\mathrm{d}\mathrm{B} $. The reported runtime is averaged over 500 realizations. The computation latency against different numbers of array elements M is shown in Fig. 8. Although standard MUSIC obtains almost optimal performance in the target detection, it has $ \mathcal{O}\left({M}^{3}\right) $ time complexity because of applying SVD to large CM. Thus, with increasing $ M $, it produces an intolerable latency which makes it impractical for real-time target detection. The same problem (i.e., the high latency) is observed with MVDR and MinNorm algorithms due to their needs for taking the inversion of CM and decomposition of large CM, respectively. The processing time of PM-US, PM-NUS, and propagator methods are almost comparable and noticeably lower than previous algorithms due to their $ \mathcal{O}\left({M}^{2}L\right) $ time complexity.

Figure 8.  Execution time against the number of array elements $ M $.

As shown by Fig. 8, our E-MUSIC algorithm is the fastest among the compared algorithms and its time complexity is linearly scalable with the number of array elements, $ \mathcal{O}\left({z}^{2}M\right) $. In particular, the computation latency of the suggested method is less than that of the standard MUSIC algorithm by the orders of the magnitude, meanwhile, a significant gap of processing time exists between the proposed and the other compared algorithms, when the array size increases. So the proposed method can be seen as a promising candidate for the real-time target detection in FMCW radar systems.

Here we conduct a numerical example to compare the proposed method with the well-known subspace algorithms in terms of the resolution of DOA estimation. We configure the simulation parameters as follows: M = 32, L = 5, z = L + 2, N=100, and $ \mathrm{S}\mathrm{N}\mathrm{R}=0\mathrm{ }\mathrm{d}\mathrm{B} $. In this scenario, the two sets of two closely-spaced targets located in $\left[\pm 49{\circ}{\rm{,}}\;\pm 55{\circ}\right]$ represent a severe test for estimation resolution. As shown in Fig. 9, the standard MUSIC algorithm obtains the highest resolution due to the use of the exact noise subspace achieved by applying full SVD to the factorization of large CM. The proposed algorithm computes an approximated spectrum that closely follows the pseudo-spectrum of the standard MUSIC algorithm. As such, the suggested algorithm detects the positions of all the targets with the resolution close enough to that of the standard MUSIC algorithm. Although the MinNorm algorithm generates five peaks, its detection capability is less accurate and seems to deviate from the actual target’s positions. Besides, the generated side-fluctuations in its spectrum and its low-noise immunity could result in high estimation errors, particularly in low $ \mathrm{S}\mathrm{N}\mathrm{R} $ regions. Furthermore, other compared algorithms (i.e., MVDR, propagator, PM-US, and PM-NUS) severely suffer from separating closely-located targets. That is to say, they detect the two closely-located targets as a single target.

Figure 9.  Comparison between the estimated pseudo-spectrums of the proposed algorithm and other popular subspace methods.

To further confirm that the resolution of the proposed algorithm closely tracks standard MUSIC, we assume that the two neighbor targets illustrated in Fig. 9 become even closer to each other and separated by only 1°. That is to say, the two signals impinge upon ULA from directions ${\theta }_{1}$ = 0° and ${\theta }_{2}$ = 1°. Using the same simulation parameters mentioned above, except L = 2, we apply both standard MUSIC and E-MUSIC algorithms to detect these two adjacent targets over multiple trials. As verified by Fig. 10, the performance of the proposed algorithm is almost the same as standard MUSIC in this challenging test. Therefore, in addition to its low complexity, the proposed algorithm is highly accurate in detecting closely located target positions. In contrast to its counterparts, our method does not sacrifice accuracy for complexity, hence, it resolves the contraction problem between complexity and accuracy in FMCW radar signal processing.

Figure 10.  Comparing the resolution capabilities of standard MUSIC and proposed algorithms for detecting two targets separated by 1° over multiple trials.

In this section, we examine the accuracy of the proposed algorithm along with several subspace methods based on different values of $ \mathrm{S}\mathrm{N}\mathrm{R}\mathrm{s} $. In this scenario, we assume that M = 32, L = 4, z = L + 2, and N = 100. To investigate the impact of $ \mathrm{S}\mathrm{N}\mathrm{R} $ on the estimation accuracy, $ \mathrm{S}\mathrm{N}\mathrm{R} $ is postulated to vary from −5 dB to 20 dB with 5 dB increments. To guarantee fair results, one-thousand trials of Monte Carlo simulations (K = 1000) are conducted independently at each step. The 1000 sets of four DOAs are randomly generated within the angular range (−90°, 90°) and these angles are equally applied to all the algorithms to ensure the equitable comparison. Then, ARMSE and corresponding PSD at each $ \mathrm{S}\mathrm{N}\mathrm{R} $ level are calculated and the obtained results are plotted as shown in Fig. 11 and Fig. 12, respectively.

Figure 11.  ARMSE against $ \mathrm{S}\mathrm{N}\mathrm{R} $ variations.

Figure 12.  PSD against $ \mathrm{S}\mathrm{N}\mathrm{R} $ variations.

The figures illustrate that for all the examined methods, the accuracy of the DOA estimation is proportional to the $ \mathrm{S}\mathrm{N}\mathrm{R} $ values. It is noteworthy that the accuracy of the proposed method closely tracks standard MUSIC and they become almost identical in high $ \mathrm{S}\mathrm{N}\mathrm{R} $ regions. Furthermore, compared with other methods, the devised algorithm attains the highest accuracy through all the considered $ \mathrm{S}\mathrm{N}\mathrm{R} $ values. Note that, our method achieves such high accuracy with substantially reducing the time complexity.

As illustrated in Fig. 3, a significant gap of the signal processing cost exists between the existing subspace-based algorithms and the suggested method, when the array size is large. So, it is interesting to observe the estimation accuracy with various values of $ M $. Thus, in this scenario, we compute ARMSE and PSD as the function of the array elements, that is, $M=10{\rm{,}}\;15{\rm{,}}\;20{\rm{,}}\;25{\rm{,}}\;\mathrm{a}\mathrm{n}\mathrm{d}\;30$, and other simulation parameters are set to be the same as those shown in subsection 5.5. Applying one-thousand Monte Carlo trials, ARMSE and PSD are computed for each value of $ M $. As illustrated in Fig. 13 and Fig. 14, increasing the array size means increasing the array aperture and this ends up with the higher accuracy (i.e., lower ARMSE and higher PSD) for all the considered algorithms. Importantly, the accuracy of the proposed algorithm is approximately as high as that of standard MUSIC and outclasses those of other subspace-based methods over the entire range of $ M $ with significantly reduced time complexity. This favorable property of our algorithm makes it an appropriate candidate for (massive) MIMO automotive radar systems.

Figure 13.  ARMSE vs. the number of array elements.

Figure 14.  PSD vs. the number of array elements.

In this paper, a low-cost yet high-resolution target detection method for the automotive FMCW radar is investigated. The fundamental problem with existing subspace algorithms is the computational complexity resulted from the decomposition of large CM. Whereas, other simplified techniques, such as MinNorm, propagator, and PMs, suffer from limited accuracy. Therefore, accurate target detection with low computational costs remains a long-established dilemma in radar signal processing. To address this, we have proposed an efficient MUSIC based algorithm through adopting the randomized matrix approximation, whereby observed CM is approximated by three low-dimensional sketches. Relying on such matrix approximation, the approximated (rather than exact) signal subspace is computed and the subspace computation is accelerated by the orders of the magnitude without sacrificing much accuracy. The theoretical analysis shows that, in high SNR scenarios, the approximated pseudo-spectrum generated by our algorithm is almost as accurate as standard MUSIC. The conducted numerical study also justifies that although the accuracy of the proposed algorithm in acquiring random DOAs just closely tracks standard MUSIC, it is tremendously faster than standard MUSIC and other popular simplified techniques. That is to say, our algorithm enjoys not only linear complexity but also high accuracy, resolving the contradictory relation between accuracy and complexity in radar signal processing. This will guarantee high-resolution real-time target detection in emerging automotive radar systems.

Here, the bound is established for the spectrum norm error. Let us begin with a definition of an auxiliary $\small {\boldsymbol{\mathfrak{R} }}$, and a corresponding sketching matrix $\small {\bf{Ͼ}}$

Proposition 1. For a given matrix M and a unitary matrix U, the following equality holds

where $\small \mathbf{U}=\left[{\mathbf{U}}_{1}\;\; {\mathbf{U}}_{2}\right]$.

Proof: Let $\small {\mathbf{P}=\mathbf{U}}^{\mathrm{H}}{\mathbf{P}}_{\mathbf{M}}\mathbf{U}$. Needless to say, P is a Hermitian matrix and OP that satisfies $\small {\mathbf{P}}^{2}=\mathbf{P}$. Consequently,

As OP is strictly identified by its range, it can be concluded that $\small \mathbf{P}={\mathbf{P}}_{{\mathbf{U}}^{\mathrm{H}}\mathbf{M}}$.■

Recalling the unitary invariance property spectrum norm [45] and according to Proposition 1, the following equality is obtained:

Since $\small {\bf{Ͼ}}$ $\small ={\boldsymbol{\mathfrak{R}}}{\Cambriabfont\text{Ω}}={\mathbf{U}}^{\mathrm{H}}\mathbf{R}{\Cambriabfont\text{Ω}}={\mathbf{U}}^{\mathrm{H}}\mathbf{C}$, then

Thus, our goal now is to prove

Assuming that the value of $\small L$ is selected such that all the elements of $\small {\Cambriabfont\text{Σ}}_{1}$ are strictly non-negative. This implies that the entries of $\small {\Cambriabfont\text{Σ}}_{2}$ become zero. Owing to the full row-rank property of $\small {\mathbf{V}}_{1}^{\mathrm{H}}$ and $\small {\Cambriabfont\text{Ω}}_{1}$, we obtain the following equality:

Our argument is primarily built based on the principles of the perturbation theory [53]. To understand these principles, we start with a matrix $\small \mathbf{W}$ related to $\small {\bf{Ͼ}}$:

Due to the full row-rank of $\small {{\Cambriabfont\text{Σ}}_{1}{\Cambriabfont\text{Ω}}}_{1}$,

where $\small {\mathbf{I}}_{L}$ is an $\small L\times L$ identity matrix and 0 is an $\small (M-L)\times L$ zero matrix. The right-hand side of (29) is full column rank, therefore, we are allowed to apply (21) for OP, which directly yields

More specifically, $\small \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathbf{W}\right)$ is within the first $\small L$-coordinates. This follows that $\small \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathbf{W}\right)$ spans the same subspace as the $\small L$-dominant singular vectors of $\small {\boldsymbol{\mathfrak{R}}}$. Consequently, the action of $\small {\boldsymbol{\mathfrak{R}}}$ is captured by the range of $\small \mathbf{W}$, this is what we want from range $\small \left(\bf{Ͼ}\right)$.

Re-visiting $ \mathbf{Ͼ} $ in (25), we obtain the following matrix:

whereby

Equation (31) guarantees that $\small \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\mathbf{Z}\right)\subset \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\left(\bf{Ͼ}\right)$. Following the guidance of Proposition 8.5 in Ref. [42], the approximation error satisfies

We then square the last inequality to obtain

To understand the last identity, one may need to relook (25). Thus, we can complete the proof by constructing a bound for the last right-hand side identified in (32). To proceed, (31) shows that Z has full column ranks. Therefore,

After some manipulations, the complementary projector verifies

According to Theorem 7.2.6 in Ref. [54], the following relation holds:

Let us re-write (33) as

where

To approach the identity in (32), we modify the above expression as

From the rule of conjugation, the left-hand side of (34) is positive semi-definite. This directly implies that the right-hand side is also positive semi-definite. It is shown by Theorem 7.7.7 in Ref. [54] that

then

From the above relation

Subsuming $\small \mathbf{D}$ with $\small {{\Cambriabfont\text{Σ}}_{2}{\Cambriabfont\text{Ω}}}_{2}{{\Cambriabfont\text{Ω}}}_{1}^{\dagger }{\Cambriabfont\text{Σ}}_{1}^{-1}$ yields

Referring to the fact that $\small \left\|{\Cambriabfont\text{Σ}}_{2}\right\|={\sigma }_{L+1}$, one can re-write the last inequality as

Finally, substituting (35) into (32) completes the proof. ■

The authors declare no conflicts of interest.