Joint Range-Doppler-Angle Estimation for OFDM-Based RadCom System via Tensor Decomposition

Abstract

Radar and communication (RadCom) systems have received increasing attention due to their high energy efficiency and spectral efficiency. They have been identified as green communications. This paper is concerned with a joint estimation of range-Doppler-angle parameters for an orthogonal frequency division multiplexing (OFDM) based RadCom system. The key idea of the proposed method is to derive different factor matrices by the tensor decomposition method and then extract parameters of the targets from these factor matrices. Different from the classical tensor decomposition method via alternating least squares or higher-order singular value decomposition, we adopt a greedy based method with each step constituted by a rank-1 approximation subproblem. To avoid local extremum, the rank-1 approximation is solved by using a multiple random initialized tensor power method with a comparison procedure followed. A parameterized rectification method is also proposed to incorporate the inherent structures of the factor matrices. The proposed algorithm can estimate all the parameters simultaneously without parameter pairing requirement. The numerical experiments demonstrate superior performance of the proposed algorithm compared with the existing methods.

This work is licensed under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

1. Introduction

The integrated radar and communication (RadCom) system has received much attention in recent years. By using a joint waveform, the occupied spectrum can be used efficiently and both radar and communication functions can be operated simultaneously. Due to the fact that the signal energy and frequency spectrum can be used efficiency and cognitively, the RadCom system is considered as a green communication system, which is a relatively new research discipline [1–4]. Such a RadCom system has been reported in many references [5–8]. Specifically, to embed the communication information into the radar waveforms efficiently, the performance of typical orthogonal frequency division multiplexing (OFDM) waveforms was analyzed [5]. Then, the single-input multiple-output (SIMO) and multiple-input multiple-output (MIMO) scenarios were extended in the RadCom system in order to estimate the directions of arrival (DOAs) of the targets [9, 10]. In general, the estimation methods can be classified into sequential methods and simultaneous methods. The existing radar systems usually employ the sequential methods for computation cost reduction. However, estimation in separate dimensions encounters the pair-matching problem for different parameters, as well as signal-to-noise ratio (SNR) loss. With the rapid development of computing power, the simultaneous methods receive more attention. Since the simultaneous method can recover multidimensional parameters at the same time, how to avoid the pair-matching procedure and improve the performance is of great importance.

The existing joint estimation methods mainly focus on the problems in Doppler and angle domains. The algorithms can be mainly divided into two groups: subspace-based algorithms and sparse representation (SR) based algorithms. In general, the subspace-based methods need a reasonably large number of snapshots and high enough SNR to implement the eigenvalue decomposition (EVD) with a desirable performance. In recent years, the SR-based techniques exploit the sparsity of the radar target scenarios. However, most of the SR-based techniques encounter the grid mismatch problem, which is caused by the solutions on discrete grids. The performance of such algorithms is directly affected by the grid resolution. On the other hand, these SR-based methods mainly focus on a one-dimensional problem and are usually extended to multidimensional cases by stacking operation [11, 12]. Nevertheless, the stacking operation ignores the inherent multidimensional structure of the received data.

Tensor based methods have been applied in radar applications [13]. With the benefits of multidimensional modeling and algorithms, the dimension of radar parameter estimation problem can be reduced and solved easily. The mainstream method is to convert the multidimensional problem into multiple one-dimensional problems with a low computational complexity. High-order singular value decomposition (HOSVD) algorithm and canonical polyadic (CP) decomposition algorithm, also known as CANDECOMP/PARAFAC decomposition, play an important role in processing multidimensional data. The alternating least squares (ALS) algorithm is still a workhorse for solving the CP decomposition problem [14, 15]. However, the ALS algorithm is faced with the troubles of local minimum and disappointing convergence properties. Another kind of algorithms is greedy algorithms, which is also known as rank-1 deflation. It is known that the greedy algorithms cannot generalize to tensor fields straightly [16]. In [17], the authors proposed a deflation method with a constraining procedure after each step. In [18], a similar method was proposed based on successive rank-1 approximations and an iterative process followed for eliminating the residue. The rank-1 approximation subproblem is usually computed by means of noniterative methods, including truncated high-order singular value decomposition (T-HOSVD) and sequential rank-one approximation and projection (SeROAP) [19]. However, they can only provide suboptimal solutions in spite of the low computational complexity.

In this paper, we introduce the tensor modeling for monostatic OFDM-SIMO based RadCom system. A data tensor is constructed from the demodulated OFDM symbols. Assuming a scenario with point scattering targets, the CP decomposition model is used to decompose the data tensor. Greedy CP decomposition (GCPD) algorithm combined with multiple random initialized tensor power method (TPM) is proposed for CP decomposition. Capitalizing on the inherent structure of the factor matrices, we present a parameterized rectification (PR) method to improve the target detection performance. The proposed algorithm deals with the received signals directly without multidimensional peak searching, covariance matrix estimation, or eigen-decomposition procedures which may bring error accumulation. Multidimensional parameter pairing is fulfilled automatically, avoiding the performance degradation caused by wrong pairing. The contributions of this paper can be summarized as follows:

(i)

A tensor model for OFDM-SIMO based RadCom system is proposed in order to jointly estimate target parameters in the range-Doppler-angle domain.

(ii)

A GCPD algorithm combined with multiple random initialized TPM is proposed for tensor decomposition. A globalization procedure is introduced to avoid the locally optimal solutions. This algorithm can achieve much better performance compared with the traditional algorithms.

(iii)

A PR method is proposed to take advantage of the inherent structures of the factor matrices. The PR method can significantly improve the target detection performance, even when there are coherent targets.

The rest of the paper is organized as follows. In Section 2, the system model for OFDM-SIMO based RadCom system and the problem formulation are introduced. A novel GCPD algorithm with multiple random initialized TPM and the PR method is presented in Section 3. In Section 4, the results of simulation in a typical multitargets scene are given to verify the performance of the proposed method. Finally, in Section 5, a conclusion is drawn.

Notation: We denote the scalars and vectors with lowercase letters $(a,b,\dots )$ and bold lowercase letters $(\mathbf{a},\mathbf{b},\dots )$. The matrices are written as bold uppercase letters $(\mathbf{A},\mathbf{B},\dots )$ and the symbol for tensors are calligraphic letters $(\mathcal{A},\mathcal{B},\dots )$. The symbols $\circ$, $\ast$ denote the outer and Hadamard (element-wise) products. The transpose, conjugate, and conjugate-transpose are denoted by ${•}^T$, ${•}^{\ast }$ and ${•}^H$, respectively. $‖•‖$ denotes the Euclidean ($l_{\mathrm{2}}$) norm of a vector. ${‖•‖}_F$ denotes the Frobenius norm of a tensor.

With respect to tensor $\mathcal{A}\in {\mathbb{C}}^{K\times N\times M}$ and vectors $\mathbf{u}\in {\mathbb{C}}^K$, $\mathbf{v}\in {\mathbb{C}}^N$ and $\mathbf{w}\in {\mathbb{C}}^M$, operator $A(\mathbf{u},\mathbf{v},\mathbf{w})$ is defined as $$\begin{array}{}\text{(1)}& A\left(\mathbf{u},\mathbf{v},\mathbf{w}\right)≔{\displaystyle \sum }_{i=\mathrm{1}}^K\hspace{0.17em}{\displaystyle \sum }_{j=\mathrm{1}}^N\hspace{0.17em}{\displaystyle \sum }_{l=\mathrm{1}}^M{u}_i^{\ast }{v}_j^{\ast }{w}_l^{\ast }\mathcal{A}\left(i,j,l\right)\in \mathbb{C}.\end{array}$$In particular, when one of these vectors is absent, we have $$\begin{array}{}\text{(2)}& A\left(I,\mathbf{v},\mathbf{w}\right)≔{\displaystyle \sum }_{j=\mathrm{1}}^N\hspace{0.17em}{\displaystyle \sum }_{l=\mathrm{1}}^M{v}_j^{\ast }{w}_l^{\ast }\mathcal{A}\left(\text{:},j,l\right)\in {\mathbb{C}}^K.\end{array}$$It is similar for the definitions of $A(\mathbf{u},I,\mathbf{w})$ and $A(\mathbf{u},\mathbf{v},I)$. Some preliminaries about tensor and its corresponding decomposition are given in the appendix.

2. System Model and Problem Formulation

Consider a monostatic OFDM-SIMO based RadCom system equipped with a single antenna for transmitter and a uniform linear array for receiver, as shown in Figure 1. The receiver array consists of $N_R$ antennas uniformly spaced with the half wavelength separation, denoted by ${\Delta }_R=\lambda /\mathrm{2}$, where $\lambda$ is the wavelength of the transmitted signal.

[figure omitted; refer to PDF]

The steering vector of the array is represented as$$\begin{array}{}\text{(3)}& a_R\left(\theta \right)={\left[\begin{array}{cccc}{a}_R^{\mathrm{0}}\left(\theta \right)& a_R^{\mathrm{1}}\left(\theta \right)& \cdots & a_R^{N_R-\mathrm{1}}\left(\theta \right)\end{array}\right]}^T,\end{array}$$where $a_R^m\left(\theta \right)=e^{j\mathrm{2}\pi m{\Delta }_R\sin \left(\theta \right)/\lambda }$.

The transmitted waveform is modulated by Cyclic-Prefix OFDM (CP-OFDM) with the Quadrature Amplitude Modulation (QAM) or phase-shift keying (PSK) constellation mapping. The transmitted signal is given by$$\begin{array}{}\text{(4)}& x\left(t\right)={\displaystyle \sum }_{n=-\infty }^{\infty }\hspace{0.17em}{\displaystyle \sum }_{k\in K}{d}_k\left[n\right]e^{j{\omega }_k\left(t-nT\right)}{p}_T\left(t-nT\right),\end{array}$$where $d_k\left[n\right]$ denotes the transmitted data over the $k$th subcarrier of the $nth$ symbol with $k\in K$. ${\omega }_k$ is the angular frequency of the $k$th subcarrier. $T$ is the symbol period. $p_T(t-nT)$ is the shaping filter of the transmitter, which is usually a time-domain rectangular filter for the general OFDM realization.

Suppose that there are $H$ nonfluctuating (according to the Swerling-0 model) far-field point targets. The parameters of a target relative to the RadCom system are given by the quadruple $\{{\alpha }_h,r_h,v_h,{\theta }_h\}$, where ${\alpha }_h$, $r_h$, $v_h$, and ${\theta }_h$ represent the complex amplitude, the range, the radial velocity, and the angle of arrival of the $h$th target, respectively. It is more convenient to equivalently consider the quadruple $\{{\alpha }_h,{\tau }_h,f_h,\sin ({\theta }_h)\}$ of complex amplitude, time delay, Doppler shift, and normalized angle parameters, where ${\tau }_h=\mathrm{2}{r}_h/c$ and $f_h=\mathrm{2}{v}_h/\lambda$. The delay, Doppler shift, and angle information can always be transformed back into the physical coordinates. In this paper, we consider the gridless scenario. That is, the true targets are likely to locate at any position in the delay-Doppler-angle domain.

The received continuous signals $y_m(t)$ at antennas $m=\mathrm{0},\dots ,N_R-\mathrm{1}$ consist of the superposition of the reflections from the targets of the transmitted probing signals, as well as the additive noise (includes the thermal noise, jamming and clutter), which are given by$$\begin{array}{}\text{(5)}& y_m\left(t\right)={\displaystyle \sum }_{h=\mathrm{1}}^H{\alpha }_h{e}^{j\mathrm{2}\pi m{\Delta }_R\sin \left({\theta }_h\right)/\lambda }{e}^{j\mathrm{2}\pi f_{h}t}x\left(t-{\tau }_h\right)+w_m\left(t\right),\end{array}$$where $w_m(t)$ is additive white Gaussian noise (AWGN) with zero-mean and variance ${\sigma }^{\mathrm{2}}$ corresponding to the $m$th receiving antenna. ${\alpha }_h$ is the complex amplitude of the $h$th target affected by path loss, scattering, and processing gains.

Since the information is modulated in frequency domain with OFDM, the demodulated data corresponds to the $m$th receive antenna that is given by$$\begin{array}{}\text{(6)}& {\widehat{d}}_k^{\left(m\right)}\left[n\right]={\displaystyle \sum }_{h=\mathrm{1}}^H{d}_k\left[n\right]v+{\tilde{w}}_m\left[n\right],\end{array}$$where $v=e^{j\mathrm{2}\pi m{\Delta }_R\sin ({\theta }_h)/\lambda }{e}^{-j{\omega }_k{\tau }_h}{e}^{j\mathrm{2}\pi f_{h}nT}$.

By implementing the CP decomposition model in Definition 4, the demodulated data in (6) can be formatted as a third-order tensor:$$\begin{array}{}\text{(7)}& \widehat{\mathcal{D}}=\mathcal{D}\ast \mathcal{V}+\tilde{\mathcal{W}},\end{array}$$where $\widehat{\mathcal{D}}\in {\mathbb{C}}^{K\times N\times M}$. $\mathcal{D}$ is a tensor from the transmitted data that is duplicated in the third mode. $\mathcal{V}$ represents the response of the targets, expressed as$$\begin{array}{}\text{(8)}& \mathcal{V}={\displaystyle \sum }_{h=\mathrm{1}}^H{\mathbf{a}}_h\circ {\mathbf{b}}_h\circ {\mathbf{c}}_h,\end{array}$$where$$\begin{array}{c}\text{(9)}& {\mathbf{a}}_h={\left[\begin{array}{ccc}{e}^{-j{\omega }_{\mathrm{1}}{\tau }_h}& \cdots & e^{-j{\omega }_K{\tau }_h}\end{array}\right]}^T,{\mathbf{b}}_h={\left[\begin{array}{ccc}{e}^{j\mathrm{2}\pi f_h\mathrm{0}·T}& \cdots & e^{j\mathrm{2}\pi f_h\left(N-\mathrm{1}\right)T}\end{array}\right]}^T,\\ {\mathbf{c}}_h={\left[\begin{array}{ccc}{e}^{j\mathrm{2}\pi \mathrm{0}·{\Delta }_R\sin \left({\theta }_h\right)/\lambda }& \cdots & e^{j\mathrm{2}\pi \left(M-\mathrm{1}\right){\Delta }_R\sin \left({\theta }_h\right)/\lambda }\end{array}\right]}^T,\end{array}$$and $\tilde{\mathcal{W}}$ is the rearranged noise.

The objective is to estimate the target parameters $\{{\alpha }_h,{\tau }_h,f_h,\sin ({\theta }_h)\}$ from the demodulated data $\widehat{\mathcal{D}}$. Traditionally, the targets are firstly detected in the delay-Doppler plane with the matched filter, and then the angle estimation is performed. However, the sequential technique needs the pair-matching procedures to obtain the one-to-one relationship among the delay, Doppler shift, and angle.

Due to the sparse nature of the radar scene, the number of targets, $H$, is usually small relative to the dimensions of the tensor. Hence the response of the targets has an intrinsic low-rank structure. Different from the traditional sequential estimation methods, we hope to achieve a joint estimation procedure with the low-rank structure of the target response.

3. Parameter Estimation via Low-Rank Tensor Approximation

In this section, we present the joint parameter estimation algorithm using tensor decomposition. The estimation procedure includes two basic stages: (1) target separation, and (2) parameter estimation. The former is achieved by CP decomposition and the latter by a correlation-based estimation with the decomposed factor matrices.

In order to eliminate the influences of the transmitted user information in the RadCom system, we normalize the received data tensor with the transmitted data tensor by element-wise product as follows:$$\begin{array}{}\text{(10)}& \check{\mathcal{D}}=\frac{\widehat{\mathcal{D}}\ast {\mathcal{D}}^{\ast }}{\mathcal{D}\ast {\mathcal{D}}^{\ast }}=\mathcal{V}+\check{\mathcal{W}},\end{array}$$where $\check{W}=\left(\tilde{\mathcal{W}}\ast {\mathcal{D}}^{\ast }\right)/\left(\mathcal{D}\ast {\mathcal{D}}^{\ast }\right)$ is an i.i.d. white Gaussian noise with the same distribution with $w_m(t)$.

Since the number of targets is usually small relative to the dimensions of the data tensor, $\check{\mathcal{D}}$ is low-rank but contaminated by additive noise. Therefore the target separation can be obtained by performing a CP decomposition of the normalized data tensor $\check{\mathcal{D}}$.

Here we assume that the number of targets $H$ is known or estimated in advance. Then, the CP decomposition can be accomplished by solving$$\begin{array}{}\text{(11)}& \underset{\widehat{\alpha },\widehat{\mathbf{A}},\widehat{\mathbf{B}},\widehat{\mathbf{C}}}{\min }\hspace{1em}{‖\check{\mathcal{D}}-\sum _{h=\mathrm{1}}^H{\widehat{\alpha }}_h{\widehat{\mathbf{a}}}_h\circ {\widehat{\mathbf{b}}}_h\circ {\widehat{\mathbf{c}}}_h‖}_F^{\mathrm{2}},\text{(12)}& \text{s.t.}\hspace{1em}‖{\widehat{\mathbf{a}}}_h‖=‖{\widehat{\mathbf{b}}}_h‖=‖{\widehat{\mathbf{c}}}_h‖=\mathrm{1},\hspace{1em}\forall h=\mathrm{1},\dots ,H,\end{array}$$where $\widehat{\alpha }=\left[{\widehat{\alpha }}_{\mathrm{1}},{\widehat{\alpha }}_{\mathrm{2}},\dots ,{\widehat{\alpha }}_H\right]$, $\widehat{\mathbf{A}}=\left[{\widehat{\mathbf{a}}}_{\mathrm{1}},{\widehat{\mathbf{a}}}_{\mathrm{2}},\dots ,{\widehat{\mathbf{a}}}_H\right]$, $\widehat{\mathbf{B}}=\left[{\widehat{\mathbf{b}}}_{\mathrm{1}},{\widehat{\mathbf{b}}}_{\mathrm{2}},\dots ,{\widehat{\mathbf{b}}}_H\right]$, $\widehat{\mathbf{C}}=\left[{\widehat{\mathbf{c}}}_{\mathrm{1}},{\widehat{\mathbf{c}}}_{\mathrm{2}},\dots ,{\widehat{\mathbf{c}}}_H\right]$.

Remark 1.

The uniqueness under mild conditions is a key feature of CP decomposition. The CP decomposition of the tensor $\stackrel{-}{\mathcal{D}}$ is essentially unique when the following sufficient condition is satisfied ([21], Theorem 10)$$\begin{array}{}\text{(13)}& k\left(\mathbf{A}\right)+k\left(\mathbf{B}\right)+k\left(\mathbf{C}\right)\ge \mathrm{2}H+\mathrm{2},\end{array}$$where $k(\mathbf{M})$ is the Kruskal rank of the matrix $\mathbf{M}$ and is defined as the largest integer $k$ such that any $k$ columns of $\mathbf{M}$ are linearly independent. When multiple noncoherent targets exist, i.e., different ranges, velocities, and angles, those three factor matrices are always full column rank so that (13) always satisfies.

Remark 2.

It is clear that when there exists coherent targets, e.g., targets with the same range but different velocities and angles, (13) will not satisfy any more. In fact, the CP decomposition cannot ensure uniqueness and correctness in this case with third-order tensors even under a much milder condition ([21], Theorem 9). To improve the percentage of successful decomposition, we utilize inherent structures of the factor matrices. Refer to Section 3.3 for detail.

The most commonly used algorithm for solving the CP decomposition is ALS [22], which is quite simple and can be executed by updating each factor matrix alternately in each iteration. The ALS algorithm is extremely fast but not stable, and the global optimal solution is hard to obtain. It has been shown that tensors of order 3 or higher fail to have best rank-$r$ approximation for $r\ge \mathrm{2}$ [23]. Fortunately, the rank-1 tensor approximation always exists and can be calculated by the high-order power method (HOPM), also known as tensor power method (TPM) [24].

Our proposed tensor decomposition method includes two main parts: (1) rank-1 tensor decomposition based on multiple random initialized TPM and (2) greedy iterations for removing the estimated components as well as residual error.

3.1. Tensor Rank-1 Approximation

There are several rank-1 tensor approximation algorithms, such as T-HOSVD and SeROAP. However, they are all suboptimal algorithms in spite of the low computational complexity. In this paper, we employ the TPM algorithm for the quasi-optimal rank-1 tensor approximation.

The problem of best rank-1 approximation of tensor $\mathcal{R}$ can be formulated as$$\begin{array}{}\text{(14)}& \underset{\widehat{\alpha },\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}}}{\min }\hspace{1em}{‖\check{\mathcal{R}}-\widehat{\alpha }\widehat{\mathbf{a}}\circ \widehat{\mathbf{b}}\circ \widehat{\mathbf{c}}‖}_F^{\mathrm{2}},\text{(15)}& \text{s.t.}\hspace{1em}‖\widehat{\mathbf{a}}‖=‖\widehat{\mathbf{b}}‖=‖\widehat{\mathbf{c}}‖=\mathrm{1}.\end{array}$$

This can be efficiently solved by TPM iterations in the complex domain, given as$$\begin{array}{}\text{(16)}& {\widehat{\mathbf{a}}}^{\left(t+\mathrm{1}\right)}=\frac{R\left(I,{\widehat{\mathbf{b}}}^{\left(t\right)},{\widehat{\mathbf{c}}}^{\left(t\right)}\right)}{‖R\left(I,{\widehat{\mathbf{b}}}^{\left(t\right)},{\widehat{\mathbf{c}}}^{\left(t\right)}\right)‖},\text{(17)}& {\widehat{\mathbf{b}}}^{\left(t+\mathrm{1}\right)}=\frac{R\left({\widehat{\mathbf{a}}}^{\left(t+\mathrm{1}\right)},I,{\widehat{\mathbf{c}}}^{\left(t\right)}\right)}{‖R\left({\widehat{\mathbf{a}}}^{\left(t+\mathrm{1}\right)},I,{\widehat{\mathbf{c}}}^{\left(t\right)}\right)‖},\text{(18)}& {\widehat{\mathbf{c}}}^{\left(t+\mathrm{1}\right)}=\frac{R\left({\widehat{\mathbf{a}}}^{\left(t+\mathrm{1}\right)},{\widehat{\mathbf{b}}}^{\left(t+\mathrm{1}\right)},I\right)}{‖R\left({\widehat{\mathbf{a}}}^{\left(t+\mathrm{1}\right)},{\widehat{\mathbf{b}}}^{\left(t+\mathrm{1}\right)},I\right)‖}.\end{array}$$

The initialization vectors $\left({\widehat{\mathbf{a}}}^{(\mathrm{0})},{\widehat{\mathbf{b}}}^{(\mathrm{0})},{\widehat{\mathbf{c}}}^{(\mathrm{0})}\right)$ are usually randomly selected, where the values of ${\widehat{\mathbf{a}}}^{(\mathrm{0})}$ and ${\widehat{\mathbf{b}}}^{(\mathrm{0})}$ are uniformly chosen in the complex domain and vector ${\widehat{\mathbf{c}}}^{(\mathrm{0})}$ is calculated by (18). Also, they can be given by any other algorithms, such as HOSVD and SeROAP. The iterations with (16)-(18) will repeat until $t$ is larger than the maximum number of iterations $t_{max}$. In order to reduce the computational complexity, an additional option of the algorithm is used to stop the iterations when the following criteria satisfy $$\begin{array}{}\text{(19)}& ‖{\widehat{\mathbf{a}}}^{\left(t\right)}-{\widehat{\mathbf{a}}}^{\left(t-\mathrm{1}\right)}‖•‖{\widehat{\mathbf{b}}}^{\left(t\right)}-{\widehat{\mathbf{b}}}^{\left(t-\mathrm{1}\right)}‖•‖{\widehat{\mathbf{c}}}^{\left(t\right)}-{\widehat{\mathbf{c}}}^{\left(t-\mathrm{1}\right)}‖\le t_S,\end{array}$$where $t_S$ is the stopping threshold.

On account that the rank-1 approximation is a nonconvex problem and many local optima exist, careful initialization is required for TPM iterations to ensure the convergence to the true rank-1 tensor components. Here we consider an approximate globalization procedure. Multiple randomly generated initializations are used for the TPM iterations. In order to identify the best one among these initializations, we need a projection procedure to obtain the final estimates of the vectors. This procedure is performed with the estimated vectors that is projected to the original data tensor. The vectors corresponding to the largest absolute value of these projections are selected as the final results. It can be formulated as $$\begin{array}{}\text{(20)}& \left(\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}}\right)=\arg \underset{{\widehat{\mathbf{a}}}_{\tau },{\widehat{\mathbf{b}}}_{\tau },{\widehat{\mathbf{c}}}_{\tau }}{\min }\left|\check{D}\left({\widehat{\mathbf{a}}}_{\tau },{\widehat{\mathbf{b}}}_{\tau },{\widehat{\mathbf{c}}}_{\tau }\right)\right|,\hspace{1em}\tau =\mathrm{1},\dots ,L,\end{array}$$where $L$ is the number of initializations.

Because the vectors $(\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}})$ are unit norm, the amplitude is estimated as $$\begin{array}{}\text{(21)}& \widehat{\alpha }=R\left(\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}}\right).\end{array}$$

Note that, in order to obtain the approximate global optimal solution, the operator in (20) is different from that in (21). Tensor rank-1 approximation via TPM algorithm is summarised in Algorithm 1.

Algorithm 1: Tensor rank-1 approximation via TPM.

Input: Tensor $\check{\mathcal{D}}\in {\mathbb{C}}^{K\times N\times M}$, $\mathcal{R}\in {\mathbb{C}}^{K\times N\times M}$, maximum

number of iterations $t_{max}$, number of initializations $L$

1:  for  $\tau =\mathrm{1}$  to  $L$  do

2:  Initialize unit vectors  ${\widehat{\mathbf{a}}}_{\tau }^{(\mathrm{0})}\in {\mathbb{C}}^K$, ${\widehat{\mathbf{b}}}_{\tau }^{(\mathrm{0})}\in {\mathbb{C}}^N$, and ${\widehat{\mathbf{c}}}_{\tau }^{(\mathrm{0})}\in {\mathbb{C}}^M$ as

  $\ast$ Option 1: random initialization.

  $\ast$ Option 2: Preset values or the last estimated results.

3:  for  $t=\mathrm{1}$  to  $t_{max}$  do

4:     $${\widehat{\mathbf{a}}}_{\tau }^{(t+\mathrm{1})}=\frac{R(I,{\widehat{\mathbf{b}}}_{\tau }^{(t)},{\widehat{\mathbf{c}}}_{\tau }^{(t)})}{‖R(I,{\widehat{\mathbf{b}}}_{\tau }^{(t)},{\widehat{\mathbf{c}}}_{\tau }^{(t)})‖}$$,

  $${\widehat{\mathbf{b}}}_{\tau }^{(t+\mathrm{1})}=\frac{R({\widehat{\mathbf{a}}}_{\tau }^{(t+\mathrm{1})},I,{\widehat{\mathbf{c}}}_{\tau }^{(t)})}{‖R({\widehat{\mathbf{a}}}_{\tau }^{(t+\mathrm{1})},I,{\widehat{\mathbf{c}}}_{\tau }^{(t)})‖}$$,

  $${\widehat{\mathbf{c}}}_{\tau }^{(t+\mathrm{1})}=\frac{R({\widehat{\mathbf{a}}}_{\tau }^{(t+\mathrm{1})},{\widehat{\mathbf{b}}}_{\tau }^{(t+\mathrm{1})},I)}{‖R({\widehat{\mathbf{a}}}_{\tau }^{(t+\mathrm{1})},{\widehat{\mathbf{b}}}_{\tau }^{(t+\mathrm{1})},I)‖}$$.

5:    end for

6:  end for

7:  Choose $(\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}})$ in $\left\{({\widehat{\mathbf{a}}}_{\tau }^{(t+\mathrm{1})},{\widehat{\mathbf{b}}}_{\tau }^{(t+\mathrm{1})},{\widehat{\mathbf{c}}}_{\tau }^{(t+\mathrm{1})}),\tau =\mathrm{1},...,L\right\}$

that correspond to the largest $\left|\check{D}(\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}})\right|$.

8:  Amplitude estimation: $\widehat{\alpha }=R(\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}})$

9:  return  $(\widehat{\alpha },\widehat{\mathbf{a}},\widehat{\mathbf{b}},\widehat{\mathbf{c}})$

3.2. Greedy CPD

In this section, we present the GCPD algorithm, which solves the problem of tensor decomposition in a greedy manner. The GCPD algorithm calculates the best rank-1 approximation and then removes the extracted component at each step. Since the best rank-1 approximation may not be the actual component of the tensor decomposition [16], additional iterations for refinement are employed. The idea of refinement is common for the greedy-like algorithms in the compressive sensing community.

The decomposition of tensor $\check{\mathcal{D}}$ can be formulated as $$\begin{array}{}\text{(22)}& \check{\mathcal{D}}={\displaystyle \sum }_{h=\mathrm{1}}^H{\alpha }_h{\mathbf{a}}_h\circ {\mathbf{b}}_h\circ {\mathbf{c}}_h+\mathcal{E},\end{array}$$where ${\alpha }_h{\mathbf{a}}_h\circ {\mathbf{b}}_h\circ {\mathbf{c}}_h$ corresponds to the $h\mathrm{t}\mathrm{h}$ rank-1 component with $\mathrm{1}\le h\le H$, and $\mathcal{E}$ is the residual error tensor, e.g., the additive noise. In the first round iterations, we compute the $H$ rank-1 components one by one and remove the extracted components after each computation. Let ${\alpha }_h^{\prime }{\mathbf{a}}_h^{\prime }\circ {\mathbf{b}}_h^{\prime }\circ {\mathbf{c}}_h^{\prime }$ be the extracted $h\mathrm{t}\mathrm{h}$ rank-1 component and the iterations can be formulated as $$\begin{array}{}\text{(23)}& \left({\alpha }_h^{\prime },{\mathbf{a}}_h^{\prime },{\mathbf{b}}_h^{\prime },{\mathbf{c}}_h^{\prime }\right)=\arg \underset{\alpha ,\mathbf{a},\mathbf{b},\mathbf{c}}{\min }{‖{\mathcal{E}}_{h-\mathrm{1}}-\alpha \mathbf{a}\circ \mathbf{b}\circ \mathbf{c}‖}_F^{\mathrm{2}},\hspace{1em}\mathrm{1}\le h\le H,\end{array}$$where ${\mathcal{E}}_h=\check{\mathcal{D}}-{\sum }_{p=\mathrm{1}}^h{\alpha }_p^{\prime }{\mathbf{a}}_p^{\prime }\circ {\mathbf{b}}_p^{\prime }\circ {\mathbf{c}}_p^{\prime }$.

On account that the extracted rank-1 components may not be the actual component of the tensor decomposition, the residual error ${\mathcal{E}}_H$ usually contains the decomposition error and is not identical to the original residual error, i.e., $\mathcal{E}$ in (22).

As a result, the refinement iterations are formulated as $$\begin{array}{}\text{(24)}& \left({\alpha }_h^{\prime },{\mathbf{a}}_h^{\prime },{\mathbf{b}}_h^{\prime },{\mathbf{c}}_h^{\prime }\right)=\arg \underset{\left(\alpha ,\mathbf{a},\mathbf{b},\mathbf{c}\right)}{\min }{‖{\mathcal{E}}_h-\alpha \mathbf{a}\circ \mathbf{b}\circ \mathbf{c}‖}_F^{\mathrm{2}},\hspace{1em}\mathrm{1}\le h\le H,\end{array}$$where ${\mathcal{E}}_h$ is the extracted rank-1 component in the last round iterations as well as the residual error given by ${\mathcal{E}}_h={\mathcal{E}}_H+{\alpha }_h^{\prime }{\mathbf{a}}_h^{\prime }\circ {\mathbf{b}}_h^{\prime }\circ {\mathbf{c}}_h^{\prime }$. The refinement iterations will be implemented multiple rounds. Clearly, the refinement iterations play a role on correction of the extracted components. Note that the estimated results in the last round are used as initialization for the next round iterations.

The refinement iterations will repeat until a stopping criterion is met. In this paper, we use the following criteria: $$\begin{array}{}\text{(25)}& {‖{\mathcal{E}}_H-{\mathcal{E}}_H^{\prime }‖}_F<t_G,\end{array}$$where ${\mathcal{E}}_H$ and ${\mathcal{E}}_H^{\prime }$ are the residual errors corresponding to the current and the last round of refinement iterations, respectively. $t_G$ is the stopping threshold.

The GCPD algorithm is summarised in Algorithm 2.

Algorithm 2: Greedy CPD.

Input: Tensor $\check{\mathcal{D}}\in {\mathbb{C}}^{K\times N\times M}$, number of targets $H$

1: $\mathcal{R}\leftarrow \check{\mathcal{D}}$

2:  for  $h=\mathrm{1}$  to  $H$  do

3:  Calculate $({\widehat{\alpha }}_h,{\widehat{\mathbf{a}}}_h,{\widehat{\mathbf{b}}}_h,{\widehat{\mathbf{c}}}_h)$ via Algorithm 1 (option 1).

4:    $\mathcal{R}\leftarrow \mathcal{R}-{\widehat{\alpha }}_h{\widehat{\mathbf{a}}}_h\circ {\widehat{\mathbf{b}}}_h\circ {\widehat{\mathbf{c}}}_h$

5:  end for

6:  repeat

7:    for  $h=\mathrm{1}$  to  $H$  do

8:    $\mathcal{R}\leftarrow \mathcal{R}+{\widehat{\alpha }}_h{\widehat{\mathbf{a}}}_h\circ {\widehat{\mathbf{b}}}_h\circ {\widehat{\mathbf{c}}}_h$

9:    Update $({\widehat{\alpha }}_h,{\widehat{\mathbf{a}}}_h,{\widehat{\mathbf{b}}}_h,{\widehat{\mathbf{c}}}_h)$ via Algorithm 1 (option 2).

10:  $\mathcal{R}\leftarrow \mathcal{R}-{\widehat{\alpha }}_h{\widehat{\mathbf{a}}}_h\circ {\widehat{\mathbf{b}}}_h\circ {\widehat{\mathbf{c}}}_h$

11:   end for

12:  until  a stopping criterion is met

13:  return  $\left({\widehat{\alpha }}_h,{\widehat{\mathbf{a}}}_h,{\widehat{\mathbf{b}}}_h,{\widehat{\mathbf{c}}}_h\right),\hspace{0.17em}\hspace{0.17em}h=\mathrm{1},\dots ,H.$

3.3. Target Parameter Estimation and Parameterized Rectification

We now discuss how to estimate the target parameters based on the estimated vectors from Algorithm 2. According to the definitions of these vectors in (8), each vector is characterized by the associated delay, Doppler shift, or angle of one target.

Hence, the delay of the $h\mathrm{t}\mathrm{h}$ target ${\tau }_h$ can be estimated via a correlation-based method given by$$\begin{array}{}\text{(26)}& {\widehat{\tau }}_h=\arg \underset{{\tau }_h}{\max }\frac{\left|{\widehat{\mathbf{a}}}_h^H{\mathbf{a}}_h\left({\tau }_h\right)\right|}{‖{\widehat{\mathbf{a}}}_h‖‖{\mathbf{a}}_h\left({\tau }_h\right)‖}.\end{array}$$

With the additive white Gaussian noise, the correlation-based method is indeed a maximum likelihood (ML) estimator and provides the optimal solution.

The Doppler shift and angle of each target can be obtained similarly as$$\begin{array}{}\text{(27)}& {\widehat{f}}_h=\arg \underset{f_h}{\max }\frac{\left|{\widehat{\mathbf{b}}}_h^H{\mathbf{b}}_h\left(f_h\right)\right|}{‖{\widehat{\mathbf{b}}}_h‖‖{\mathbf{b}}_h\left(f_h\right)‖},\text{(28)}& {\widehat{\theta }}_h=\arg \underset{{\theta }_h}{\max }\frac{\left|{\widehat{\mathbf{c}}}_h^H{\mathbf{c}}_h\left({\theta }_h\right)\right|}{‖{\widehat{\mathbf{c}}}_h‖‖{\mathbf{c}}_h\left({\theta }_h\right)‖}.\end{array}$$

The maximization problems in (26)-(28) involve one-dimensional search and can be performed by zero-padded FFT efficiently combining the inherent structures of these vectors.

Vectors ${\mathbf{a}}_h$, ${\mathbf{b}}_h$ and ${\mathbf{c}}_h$ are inherently determined by a few parameters. However, Algorithms 1 and 2 do not take this into account. In order to make use of this structural information, we propose a method named parameterized rectification (PR). The PR method is based on parameter estimation and structure reconstruction. That is, an estimation process performed by (26)-(28) is inserted at the end of Algorithm 1. The returned vectors are regenerated with the desired structures and the estimated parameters. The PR method is summarised in Algorithm 3. The GCPD algorithm combined with PR is abbreviated as PR-GCPD.

Algorithm 3: Parameterized rectification.

Procedure  $\mathrm{P}\mathrm{R}({\widehat{\mathbf{a}}}_h,{\widehat{\mathbf{b}}}_h,{\widehat{\mathbf{c}}}_h)$

Estimate parameters ${\widehat{\tau }}_h$, ${\widehat{f}}_h$ and ${\widehat{\theta }}_h$ via Eq. (26)-(28).

Regenerate  ${\widehat{\mathbf{a}}}_h$, ${\widehat{\mathbf{b}}}_h$ and  ${\widehat{\mathbf{c}}}_h$ with the structures in Eq. (8) and

estimated parameters.

return  $({\widehat{\mathbf{a}}}_h,{\widehat{\mathbf{b}}}_h,{\widehat{\mathbf{c}}}_h)$

Note that the PR method generally increases the rank, especially in the noisy case, and several initializations and iterations are necessary to obtain the global optimal solutions.

3.4. Computational Complexity Analysis

We use the number of complex multiplications (operations) as the complexity metric. Since GCPD is an iterative algorithm, the total complexity is unbounded. The complexity is mainly dominated by the rank-1 approximation, which is repeatedly computed several times. The major computing task of the rank-1 approximation is the TPM iteration. The computation of vectors ${\widehat{\mathbf{a}}}^{(t+\mathrm{1})}$, ${\widehat{\mathbf{b}}}^{(t+\mathrm{1})}$, and ${\widehat{\mathbf{b}}}^{(t+\mathrm{1})}$ in (16)-(18) needs $\mathrm{3}KNM+KN+KM+NM$ operations. The TPM iterations are repeated with the maximum $t_{max}$ times. Thus, the complexity of the rank-1 approximation is given by $\mathcal{O}\{\mathrm{3}{t}_{max}KNM\}$, where $\mathcal{O}\{·\}$ is the big-O notation. Since the iterations with multiple initializations can be performed in parallel, the number of initializations is not considered in the computational complexity. Assume that the maximum number of iterations for refinement in the GCPD agorithm is $t_{max}^{\prime }$, the total complexity of GCPD is given by $\mathcal{O}\{\mathrm{3}(t_{max}^{\prime }+\mathrm{1})t_{max}HKNM\}$. Note that the stopping criteria in (19) and (25) are usually applied, the actual number of operations is much smaller.

In Table 1, we summarize the computational complexity of the algorithms presented in the next section, where $k$ is a user-defined parameter. Since the PR method is performed on vectors, the number of operations per iteration needed by the PR method is negligible compared with the TPM iterations. However, on account that the PR method generally increases the rank, it usually slows down the convergence and may lead to more run time. The computational complexity is evaluated via simulation in the next section.

Table 1

Computational complexity of different algorithms.

Algorithms	Complexity
FB-RootMUSIC	$$\mathcal{O}\left\{KNM\left(K+H\left(K+N\right)\left(K+M\right)\right)\right\}$$   [20]
ALS	$$\mathcal{O}\left\{\mathrm{3}{t}_{max}HKNM\right\}$$
GCPD (SeROAP init)	$$\mathcal{O}\left\{\mathrm{2}\left(k+\mathrm{1}\right)HK{N}^{\mathrm{2}}{M}^{\mathrm{2}}+3\left(t_{max}^{\prime }+\mathrm{1}\right)t_{max}HKNM\right\}$$   [19]
GCPD (random init)	$$\mathcal{O}\left\{\mathrm{3}\left(t_{max}^{\prime }+\mathrm{1}\right)t_{max}HKNM\right\}$$

4. Numerical Results

In this section, some numerical results are used to illustrate the performance of the proposed method. Different from the related literature [18], the dimensions of the data are much larger and the simulated performance characteristics may be obviously distinct.

4.1. RadCom Parameters and Performance Metrics

The transmitted signal is modulated with CP-OFDM. The parameters for OFDM waveform generation are listed in Table 2. The receiver of the RadCom system is equipped with a uniform linear array with $M_R=\mathrm{15}$ receive antennas spaced at a distance ${\Delta }_R=\mathrm{0.5}{\lambda }_{\mathrm{0}}$. The surveillance field is in the far field of the RadCom system. In each experiment, targets are placed randomly on the predefined unambiguous range-Doppler-angle parameter space. The target amplitudes are chosen with constant absolute value and random phase. White Gaussian noise is added to the data tensor with variance ${\sigma }^{\mathrm{2}}$ corresponding to the specified output SNR. The SNR is changed from -40 dB to -20 dB in step of 1 dB.

Table 2

OFDM waveform parameters.

Parameters	Values
Carrier frequency	5.9 GHz
Subcarrier spacing	90.909 kHz
OFDM symbol length	11 us
Cyclic prefix length	1.375 us
Bandwidth	93.1 MHz
Constellation Mapping	4 QAM

In all simulations, the maximum number of iterations $t_{max}$ for the TPM method is set to be 1000. Also, the stopping threshold $t_S$ in (19) is assigned to be $\mathrm{1}{\mathrm{0}}^{-\mathrm{10}}$. For the GCPD algorithm, the stopping threshold $t_G$ in (25) is set to be $\mathrm{1}{\mathrm{0}}^{-\mathrm{8}}$. In order to avoid the lengthy refinement iterations in the GCPD algorithm, we restrict the number of refinement rounds no greater than 10. In each experiment, 64 OFDM symbols are collected for signal processing. A total number of 1000 experiments are conducted.

Here we use the following performance metrics:

(1)

Root mean-square errors (RMSE) in the delay, Doppler, and angle estimation are given by $$\begin{array}{}\text{(29)}& {\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\tau }={\left(\frac{\mathrm{1}}{HL}\sum _{l=\mathrm{1}}^L\hspace{0.17em}\sum _{h=\mathrm{1}}^H{\left({\widehat{\tau }}_{h,l}-{\tau }_{h,l}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\text{(30)}& {\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_f={\left(\frac{\mathrm{1}}{HL}\sum _{l=\mathrm{1}}^L\hspace{0.17em}\sum _{h=\mathrm{1}}^H{\left({\widehat{f}}_{h,l}-f_{h,l}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\text{(31)}& {\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\theta }={\left(\frac{\mathrm{1}}{HL}\sum _{l=\mathrm{1}}^L\hspace{0.17em}\sum _{h=\mathrm{1}}^H{\left({\widehat{\theta }}_{h,l}-{\theta }_{h,l}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

 

where $H$ is the number of targets and $L$ is the total number of experiments.

(2)

Detection probability ($P_d$): the fraction of the total number of targets that are correctly detected in an ellipsoid area. A target is correctly detected when the estimated location of the target in the delay-Doppler-angle space falls within the ellipsoid with axes equivalent to ±3 times the classical delay, Doppler, and angle resolution bins, ${\Delta }_{\tau }$, ${\Delta }_f$, and ${\Delta }_{\theta }$, respectively. $$\begin{array}{}\text{(32)}& P_d=\frac{\mathrm{1}}{HL}\left|\left\{{\widehat{\tau }}_{h,l},{\widehat{f}}_{h,l},{\widehat{\theta }}_{h,l}\mid {\Phi }_{h,l}\le {ϵ}_{th},\right\}\right|,\end{array}$$

 

where$$\begin{array}{}\text{(33)}& {\Phi }_{h,l}={\left({\left(\frac{{\widehat{\tau }}_{h,l}-{\tau }_{h,l}}{{\Delta }_{\tau }}\right)}^{\mathrm{2}}+{\left(\frac{{\widehat{f}}_{h,l}-f_{h,l}}{{\Delta }_f}\right)}^{\mathrm{2}}+{\left(\frac{{\widehat{\theta }}_{h,l}-{\theta }_{h,l}}{{\Delta }_{\theta }}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

 

$H$ and $L$ are the same as those in (29)-(31). ${ϵ}_{th}$ is the normalized tolerance factor, determined by application. We chose ${ϵ}_{th}=\mathrm{3}$, which is three times the resolution bins.

When calculating (29)-(32) in the multiple target scenarios, ${\widehat{\tau }}_{h,l}$, ${\widehat{f}}_{h,l}$, and ${\widehat{\theta }}_{h,l}$ are selected as the ones closest to ${\tau }_{h,l}$, $f_{h,l}$, and ${\theta }_{h,l}$, respectively, among all the estimated parameters of $H$ targets, i.e., $\{{\widehat{\tau }}_{h,l},{\widehat{f}}_{h,l},{\widehat{\theta }}_{h,l}\mid h=\mathrm{1,2},\dots ,H\}$.

4.2. Performance with Single Target

In this case, the problem corresponds to the rank-1 situation. Figure 2 compares different rank-1 approximation methods as well as different number random initializations for TPM method. Because the dimensions of the tensor are high in this paper, we do not compare the results with the best rank-1 approximation described in [25].

[figure omitted; refer to PDF]

In Figure 2, the number of random initializations for TPM method is selected in the collection $\mathcal{I}=\left\{\mathrm{1,30,100}\right\}$. The methods for comparison include T-HOSVD and SeROAP. The former is a classical approach for tensor decomposition and the latter was recently proposed by Alex P. da Silva et al. in [19]. TPM initialized by the results of SeROAP is also considered. In addition, the subspace-based forward-backward root-MUSIC (FB-RootMUSIC) method [20] incorporating the inherent signal structures is selected for comparison.

The target detection probability versus SNRs for different methods is shown in Figure 2. This simulation shows that the TPM method with 100 random initializations performs much better than the others. TPM with 30 random initializations performs very close to the case with 100 random initializations. The TPM method with multiple random initializations is better than that with initialization generated by the SeROAP method.

In order to further analyze the impact of the number of random initializations, another experiment is performed. The number of random initializations is selected in the collection $\mathcal{I}=\left\{\mathrm{1,5},\mathrm{10,15,20,25,30,35,40,45,50,75,100}\right\}$. Figure 3 shows the detection probability versus the number of random initializations for SNR = -30, -31, and -32 dB, respectively. The detection probability tends to a fixed value with the increasing number of initializations. When the initialization number is sufficiently large, the output is close to the global optimal solution. However, in view of the computational burden, we choose the number of initializations to be 30 as an acceptable suboptimal solution. It is worth noting that the TPM solving processes with multiple initializations can be easily realized by parallel computing, so that this is not a limiting factor for practical applications.

[figure omitted; refer to PDF]

To simplify the analysis, the number of random initializations of TPM is fixed to be 30 in the subsequent analyses.

4.3. Performance with Multiple Noncoherent Targets

Consider five targets selected uniformly at random in the RadCom’s unambiguous region. The delay, Doppler shift, and angle of each target do not overlap with others. The spacings between different targets in each dimension are larger than three times of the corresponding resolution bin. The amplitude of the targets are chosen such that they are with the same magnitude and random phases.

Figure 4 depicts the RMSE performance of the aforementioned methods in the range, Doppler and angle domains. When the SNR is low, the proposed GCPD and PR-GCPD methods with 30 random initializations have much better performance than the other methods. As the SNR increases, GCPD can achieve similar performance compared to PR-GCPD. The FB-RootMUSIC method cannot achieve similar performance compared to the GCPD and PR-GCPD methods, even at high SNRs. The ALS method performs the worst due to its instability. Figure 6 presents the detection probability of multiple targets. It can be seen that our proposed GCPD and PR-GCPD methods have better detection probability than the classical ALS and FB-RootMUSIC methods. Also, the PR-GCPD method initialized by SeROAP has an acceptable detection probability but is slightly lower than the multiple random initialized ones. The facts show that the globalization procedure with multiple random initializations, as well as the PR method, can improve performance significantly.

[figures omitted; refer to PDF]

4.4. Influence of Coherent Targets

Consider scenarios with five targets and two of them are coherent targets in the range dimension. The coherence indicates that the parameters of different targets in one dimension are the same, e.g., the same delay, Doppler shift, or angle. The coherence destroys the uniqueness conditions of CP decomposition and significantly influences the performance of parameter estimation and target detection.

From Figure 5, we can see that, in this condition, the proposed PR-GCPD method combined with multiple random initialized TPM provides the best performance. Methods except for the PR-GCPD have much worse performance on account that they cannot resolve the coherent targets robustly. Figure 7 shows the detection probability of all the methods. It is observed that the classical ALS and FB-RootMUSIC methods all have miss detection even when the SNR is high. This is primarily because these methods cannot distinguish the coherent targets. The GCPD with 30 random initializations performs better, though it cannot achieve the best performance at high SNRs. The proposed PR-GCPD method performs much better than all the other methods, and the detection probability is 1 at high SNRs.

[figures omitted; refer to PDF]

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

Figure 8 is given to evaluate the run times of different algorithms relative to the scenario existing five targets and two of them are coherent. The scenarios with noncoherent targets have similar results and will not be shown here. The run times are obtained by using a computer with Intel(R) Xeon(R) CPU E5-2682 v4 CPU, 16 GB RAM. As has been mentioned before, the proposed GCPD algorithm with multiple random initializations has the ability to run in parallel. Here we calculate the equivalent run time of GCPD by selecting the largest run time among all these initializations.

[figure omitted; refer to PDF]

From Figure 8 we can see that the required run time of FB-RootMUSIC is constant as it is a finite algorithm. The required run times of GCPD (SeROAP init) and GCPD (1 rand-init) are both smaller than that of the ALS algorithm, especially when the SNR is low. This indicates that the complexities of GCPD (SeROAP init) and GCPD (1 rand-init) are lower than that of ALS. The required run time of the GCPD (30 rand-inits) is slightly larger than that of the ALS algorithm. The reason is that the run times with different initializations are not all equal and the largest one determines the required time.

When the PR method is applied, the run times of all three PR-GCPD algorithms become larger in the low SNR region. This is mainly because that the PR method generally increases the rank and slows down the convergence. However, when the SNR is high, the PR-GCPD algorithm performs much faster. As has been mentioned, the PR method increases the detection and estimation performance of the GCPD algorithm both in low and in high SNR regions, although it may result in a higher computational cost.

5. Conclusion

In this article, we investigated joint range-Doppler-angle estimation in an OFDM-SIMO based RadCom system using CP decomposition. The signal model with tensor algebra was developed and a novel algorithm for CP decomposition was presented. Different from the classical ALS algorithm, the proposed one adopts a greedy strategy with each step solved by TPM with multiple random generated initializations and a globalization procedure. This globalization procedure alleviates the local optimal problem to some extent. A PR method was proposed to make use of the inherent structures of the factor matrices. We demonstrated that our methods can estimate parameters for multiple targets, both in noncoherent and in coherent cases, and require no pair matching. The multiple random initialized TPM can be easily realized by parallel computing and it is beneficial for realistic applications.

Conflicts of Interest

The authors declare that they have no conflicts of interest and the received funding did not lead to any conflicts of interest regarding the publication of this paper.

Acknowledgments

The research reported in this article was supported in part by the National Natural Science Foundation of China (61661028 and 61561034).

Appendix

Appendix

Definition 3 (rank-1 tensor).

A third-order tensor $\mathcal{T}$ has rank 1 if it can be expressed by the following form:$$\begin{array}{}\text{(A.1)}& \mathcal{T}=\mathbf{a}\circ \mathbf{b}\circ \mathbf{c},\end{array}$$where $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$ are three vectors, with $\mathbf{a}\in {\mathbb{C}}^{I\times \mathrm{1}}$, $\mathbf{b}\in {\mathbb{C}}^{J\times \mathrm{1}}$, $\mathbf{c}\in {\mathbb{C}}^{K\times \mathrm{1}}$, and $\mathcal{T}\in {\mathbb{C}}^{I\times J\times K}$.

Definition 4 (CP decomposition).

The CP decomposition of the third-order tensor $\mathcal{T}\in {\mathbb{C}}^{I\times J\times K}$ is a decomposition of $\mathcal{T}$ with a summation of minimal number of rank-1 tensors:$$\begin{array}{}\text{(A.2)}& \mathcal{T}={\displaystyle \sum }_{r=\mathrm{1}}^R{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r,\end{array}$$where ${\mathbf{a}}_r$, ${\mathbf{b}}_r$, ${\mathbf{c}}_r$ are the $r$th columns of factor matrices $\mathbf{A}\in {\mathbb{C}}^{I\times R}$, $\mathbf{B}\in {\mathbb{C}}^{J\times R}$, and $\mathbf{C}\in {\mathbb{C}}^{K\times R}$.

References

[1] P. Gandotra, R. K. Jha, S. Jain, "Green Communication in Next Generation Cellular Networks: A Survey," IEEE Access, vol. 5, pp. 11727-11758, DOI: 10.1109/ACCESS.2017.2711784, 2017.

[2] F. Zhou, Y. Wang, D. Qin, Y. Wang, Y. Wu, "Secure EE maximisation in green CR: guaranteed SC," IET Communications, vol. 11 no. 16, pp. 2507-2513, DOI: 10.1049/iet-com.2016.1349, 2017.

[3] F. Zhou, N. C. Beaulieu, Z. Li, J. Si, P. Qi, "Energy-efficient optimal power allocation for fading cognitive radio channels: ergodic capacity, outage capacity, and minimum-rate capacity," IEEE Transactions on Wireless Communications, vol. 15 no. 4, pp. 2741-2755, DOI: 10.1109/TWC.2015.2509069, 2016.

[4] F. Zhou, Z. Chu, H. Sun, R. Q. Hu, L. Hanzo, "Artificial Noise Aided Secure Cognitive Beamforming for Cooperative MISO-NOMA Using SWIPT," IEEE Journal on Selected Areas in Communications,DOI: 10.1109/JSAC.2018.2824622, 2018.

[5] C. Sturm, W. Wiesbeck, "Waveform design and signal processing aspects for fusion of wireless communications and radar sensing," Proceedings of the IEEE, vol. 99 no. 7, pp. 1236-1259, DOI: 10.1109/jproc.2011.2131110, 2011.

[6] J. Moghaddasi, K. Wu, "Multifunctional transceiver for future radar sensing and radio communicating data-fusion platform," IEEE Access, vol. 4, pp. 818-838, DOI: 10.1109/access.2016.2530979, 2016.

[7] X. Tian, T. Zhang, Q. Zhang, Z. Song, "High Accuracy Doppler Processing With Low Complexity in OFDM-Based RadCom Systems," IEEE Communications Letters, vol. 21 no. 12, pp. 2618-2621, DOI: 10.1109/LCOMM.2017.2746563, 2017.

[8] F. Liu, L. Zhou, C. Masouros, A. Li, W. Luo, A. Petropulu, "Towards Dual-functional Radar-Communication Systems: Optimal Waveform Design," IEEE Transactions on Signal Processing,DOI: 10.1109/TSP.2018.2847648, .

[9] W. Wiesbeck, L. Sit, M. Younis, T. Rommel, G. Krieger, A. Moreira, "Radar 2020: the future of radar systems," Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS '15), pp. 188-191, DOI: 10.1109/igarss.2015.7325731, .

[10] J. Zhang, I. Podkurkov, M. Haardt, A. Nadeev, "Efficient multidimensional parameter estimation for joint wideband radar and communication systems based on OFDM," Proceedings of the 2017 IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP 2017, pp. 3096-3100, .

[11] X. Wang, M. Huang, X. Wu, G. Bi, "Direction of Arrival Estimation for MIMO Radar via Unitary Nuclear Norm Minimization," Sensors, vol. 17 no. 4,DOI: 10.3390/s17040939, 2017.

[12] R. Heckel, "Super-resolution MIMO radar," Proceedings of the 2016 IEEE International Symposium on Information Theory, ISIT 2016, pp. 1416-1420, .

[13] D. Nion, N. D. Sidiropoulos, "Tensor algebra and multidimensional harmonic retrieval in signal processing for MIMO radar," IEEE Transactions on Signal Processing, vol. 58 no. 11, pp. 5693-5705, DOI: 10.1109/TSP.2010.2058802, 2010.

[14] T. G. Kolda, B. W. Bader, "Tensor decompositions and applications," SIAM Review, vol. 51 no. 3, pp. 455-500, DOI: 10.1137/07070111X, 2009.

[15] A. Uschmajew, "Local convergence of the alternating least squares algorithm for canonical tensor approximation," SIAM Journal on Matrix Analysis and Applications, vol. 33 no. 2, pp. 639-652, DOI: 10.1137/110843587, 2012.

[16] A. Stegeman, P. Comon, "Subtracting a best rank-1 approximation may increase tensor rank," Linear Algebra and its Applications, vol. 433 no. 7, pp. 1276-1300, DOI: 10.1016/j.laa.2010.06.027, 2010.

[17] A.-H. Phan, P. Tichavský, A. Cichocki, "Tensor deflation for CANDECOMP/PARAFAC-part I: Alternating subspace update algorithm," IEEE Transactions on Signal Processing, vol. 63 no. 22, pp. 5924-5938, DOI: 10.1109/TSP.2015.2458785, 2015.

[18] A. P. Da Silva, P. Comon, A. L. F. De Almeida, "An iterative deflation algorithm for exact CP tensor decomposition," Proceedings of the 40th IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP 2015, pp. 3961-3965, .

[19] A. P. Da Silva, P. Comon, A. L. F. De Almeida, "A Finite Algorithm to Compute Rank-1 Tensor Approximations," IEEE Signal Processing Letters, vol. 23 no. 7, pp. 959-963, DOI: 10.1109/LSP.2016.2570862, 2016.

[20] L. Huang, Y. Wu, H. C. So, Y. Zhang, L. Huang, "Multidimensional sinusoidal frequency estimation using subspace and projection separation approaches," IEEE Transactions on Signal Processing, vol. 60 no. 10, pp. 5536-5543, DOI: 10.1109/TSP.2012.2206590, 2012.

[21] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, C. Faloutsos, "Tensor decomposition for signal processing and machine learning," IEEE Transactions on Signal Processing, vol. 65 no. 13, pp. 3551-3582, DOI: 10.1109/TSP.2017.2690524, 2017.

[22] P. Comon, X. Luciani, A. L. F. de Almeida, "Tensor decompositions, alternating least squares and other tales," Journal of Chemometrics, vol. 23 no. 7-8, pp. 393-405, DOI: 10.1002/cem.1236, 2009.

[23] V. de Silva, L.-H. Lim, "Tensor rank and the ill-posedness of the best low-rank approximation problem," SIAM Journal on Matrix Analysis and Applications, vol. 30 no. 3, pp. 1084-1127, DOI: 10.1137/06066518X, 2008.

[24] L. Wang, M. T. Chu, "On the global convergence of the alternating least squares method for rank-one approximation to generic tensors," SIAM Journal on Matrix Analysis and Applications, vol. 35 no. 3, pp. 1058-1072, DOI: 10.1137/130938207, 2014.

[25] J. B. Lasserre, "Global optimization with polynomials and the problem of moments," SIAM Journal on Optimization, vol. 11 no. 3, pp. 796-817, DOI: 10.1137/S1052623400366802, 2000/01.