Mechanism of tensor decomposition

The target information is implicit in the received sample, and the received sample or sample covariance needs to be decomposed. Tensor decomposition is a means to achieve low-rank approximation of tensors. Compared with matrix factorization, the difference is: because tensor factorization mines the structure information of the data more, it has better identification or achieves more accurate decomposition. These advantages make it be introduced into the field of array signal processing to improve the performance of DOA estimation. The canonical polyadic decomposition (CPD) used in this paper is introduced below.

Consider a tensor ℓ with rank K and dimension R.

(1)$$\ell = \sum\limits_{k = 1}^K {{\bf l}_1^k \circ {\bf l}_2^k \circ \ldots \circ {\bf l}_R^k } , \;$$

where ${\bf l}_r^k$ is the kth component k = 1, 2, …, K corresponding to the rth dimension.

The purpose of CPD is to represent a tensor as a sum of rank 1 tensors. The number of rank 1 tensors is the rank of the entire tensor. The schematic diagram of CPD based on the third-order tensor is shown in Fig. 1.

Fig. 1.

Schematic diagram of the third-order tensor CPD.

CPD directly decomposes the tensor into the sum of each rank 1 tensor, and uses the structural information of the tensor to decompose it into the various factors themselves. Taking a third-order tensor as an example, if it can be expressed by (1), then the r-module matrix expansion along the three dimensions can be expressed as

(2a)$$\ell _{( 1) } = {\bf L}_1( {\bf L}_2\odot {\bf L}_3) ^T$$

(2b)$$\ell _{( 2) } = {\bf L}_2( {\bf L}_3\odot {\bf L}_1) ^T$$

(2c)$$\ell _{( 3) } = {\bf L}_3( {\bf L}_2\odot {\bf L}_1) ^T$$

where ${\bf L}_1 = [ {{\bf l}_1^1 , \;{\bf l}_1^2 , \;\ldots , \;{\bf l}_1^K } ] , \;{\bf L}_2 = [ {{\bf l}_1^1 , \;{\bf l}_2^2 , \;\ldots , \;{\bf l}_2^K } ] , \;{\bf L}_3 = [ {{\bf l}_3^1 , \;{\bf l}_3^2 , \ldots , {\bf l}_3^K } ]$ are the matrix factors that constitute the tensor. Equation (2) shows the relationship between the tensor and the matrix factors obtained after the modulus r is expanded. For the received sample of the multi-dimensional array, assuming that the array size is I ₁ × I ₂ × ⋅ ⋅ ⋅ × I _R, the sample received by an element in the array at time t can be expressed as [Reference Cheng, Wu and Poor32]

(3)$$x_{i_1, \ldots , i_R}( t) = \sum\limits_{k = 1}^K {l_{i1} \ldots l_{iR}( \mu _k) {\bf s}_k( t) + n_{i1, \ldots , iR}} , \;$$

where l _ir denotes the i _rth element in the rth dimension in the array, s_k(t) denotes the sampling at the time t of the kth signal, n(t) denotes the noise item received at the time t, and μ _k denotes the parameter vector of the kth signal. The specific form is decided. The steering vector of the kth signal corresponding to the rth dimension of the array is l_r(μ _r,k). Then the sample vector ${\boldsymbol x}( t) \in {\mathbb C}^{\Pi _{i = 1}^R I_i \times 1}$ received by the entire array at time t can be expressed as

(4)$${\boldsymbol x}( {\boldsymbol t}) {\rm} = \sum\limits_{k = 1}^K {{\bf l}( \mu _k) {\bf s}_k( t) + {\bf n}( t) } .$$

The steering vector of the entire array can be expressed by the Khatri–Rao product of the steering vectors of each dimension as

(5)$${\bf l}( \mu _k) = \mathop \odot \limits_{r = 1}^R {\bf l}_r( \mu _{r, k}) , \;$$

where l_r is the steering vector ${\bf l}_r\in {\mathbb C}^{I_r \times 1}$ corresponding to the kth signal in the rth dimension of the multi-dimensional array, I _r is the size of the rth dimensional array, and μ _r,k denotes the parameter of the kth signal in the rth dimension of the multi-dimensional array.

The received data of the entire array can be expressed as a tensor form ${\cal X} ( t) \in {\mathbb C}^{I_1 \times I_2 \times \cdots \times I_R}$. For a total of T time samples, the tensor mode of the array received data can be expressed by the cascade of sub-tensors at each time as

(6)$${\cal X} = {\cal X} ( 1) \cup _{R + 1}{\cal X} ( 2) \cup _{R + 1} \ldots \cup _{R + 1}{\cal X} ( T) , \;$$

where ${\cup} _{R + 1}$ denotes that two tensors are concatenated along the rth dimension. Thus, the received data can be constructed as

(7)$${\cal X}{\rm} = {\cal L} \times _{R + 1}\mathbf{S} + {\cal N}$$

where ${\cal X}$ is the sample tensor form ${\cal X}\in {\mathbb C}^{I_1 \times I_2 \times \cdots \times I_R \times T}$ of the multi-dimensional array, ${\cal L}$ is the steering tensor of the multi-dimensional array, whose dimension is ${\cal L}\in {\mathbb C}^{I_1 \times I_2 \times \cdots \times I_R \times K}$, S denotes the signal matrix, and ${\cal N}$ denotes the noise contained in the T samples received by the multi-dimensional array. Then the steering tensor of the kth signal can be expressed by the outer product of the steering vectors of each dimension as

(8)$${\cal L} _k( \mu _k) = {\bf l}_1( \mu _{1, k}) \circ {\bf l}_2( \mu _{2, k}) \circ \ldots \circ {\bf l}_R( \mu _{R, k}) .$$

Thus, the estimation of the array steering vector can be obtained by CPD.

Signal model of FDA-MIMO radar

Consider an FDA-MIMO scheme in which both transmit and receive arrays are uniform linear arrays, as shown in Fig. 2. The transmitted carrier frequency is f ₀, and the transmitted frequency of the nth array element can be expressed as

(9)$$f_n = f_0 + n\Delta f.$$

Fig. 2.

Schematic diagram of FDA-MIMO radar.

The transmitted signal of the nth element can be expressed as

(10)$$s_n( t) = rect\left({\displaystyle{t \over {T_p}}} \right)\psi _n( t) e^{\,j2\pi ( f_0 + n\Delta f) t}, \;$$

where T _p is the pulse width, ψ _n(t) is the baseband envelope of the nth transmitting array element signal, which satisfies the orthogonality condition as

(11)$$\int {\psi _{n1}^\ast ( t) \cdot \psi _{n2}( t-\tau ) } \,dt = 0, \;\quad n1\ne n2, \;\;\forall \tau .$$

Assuming that there are K targets in the far field, the signal returned from the position (r, θ) and received by the mth received element can be expressed as

(12)$$s_{n, m}( t) = rect\left({\displaystyle{{t-\tau_{n, m}} \over {T_p}}} \right)\psi _n( t-\tau _{n, m}) e^{\,j2\pi ( f_0 + n\Delta f) ( t-\tau _{n, m}) }, \;$$

where τ _n,m is the propagation delay from the nth transmitted element to the mth received element, which can be expressed as

(13)$$\tau _{n, m} = {{2r_k} / c}-{{nd_T\sin \theta _k} / c}-{{md_R\sin \theta _k} / c}\;( 1 \le k \le K) .$$

Thus, after signal processing at the receiving end, the output signal is obtained as

(14)$$ \eqalign{& x_{n, m}( t) \approx rect\left({\displaystyle{{t-\tau_{n, m}} \over {T_p}}} \right) \cr& \xi _se^{\,j2\pi ( {-}f_0( 2r_k/c) -n\Delta f( 2r_k/c) + f_0( ( nd_T\sin \theta _k) /c) + f_0( ( md_R\sin \theta _k) /c) ) }, \;} $$

where ξ _s is the signal complex coefficient after matched filtering. When only considering the equivalent transmit beampattern, the above formula can be simplified as

(15)$$x_{n, m}( t) \approx rect\left({\displaystyle{{t-\tau_{n, m}} \over {T_p}}} \right)\xi _se^{\,j2\pi ( {-}f_0( 2r_k/c) -n\Delta f( 2r_k/c) + f_0( ( nd_T\sin \theta _k) /c) ) }.$$

Then the transmit steering vector can be expressed as

(16)$$\eqalign{{\bf S}_t( r_k, \;\theta _k) =\; & {\bf a}( r_k, \;\theta _k) \cr =\; & [ 1, \;e^{\,j2\pi [ -\Delta f 2r_k/c + f_0 d_T\sin \theta _k /c ] }, \;\ldots , \; \cr& e^{\,j2\pi [ -( N-1) \Delta f 2r_k/c + ( N-1) f_0 d_T\sin \theta _k /c ] }] ^{\rm T}.} $$

When the frequency increment Δf = 0 is considered, the transmit steering vector is simplified to

(17)$${\bf S}_{t0}( \theta _k) = [ 1, \;e^{\,j2\pi f_0 d_T\sin \theta _k /c }, \;\ldots , \;e^{\,j2\pi ( N-1) f_0 d_T\sin \theta _k /c }] ^{\rm T}.$$

Double-pulse scheme for the localization of targets

A method for target localization using a double pulse frequency diverse array is proposed in [Reference Wang and Shao12]. In this method, two steps are required to separately estimate the angle and range parameters of the target.

The first step is to use the FDA pulse with Δf = 0 to estimate the target angle parameter. In this case, the data model of the received signal can be written as

(18)$${\boldsymbol x} = \alpha _0{\bf S}_{t0}( \theta _k) + {\bf n}{\rm , \;}$$

where α ₀ represents the complex amplitude of the signal, and ${\bf \boldsymbol n}$ is the additive white noise vector. The weight vector using non-adaptive beamforming is expressed as

(19)$${\bf w} = {\bf S}_{t0}( \theta _k) .$$

Thus, the formula used to estimate the angle parameter of the target is

(20)$$\hat{\theta }_k = \arg \left\{{\mathop {\max }\limits_\theta {\vert {{\bf w}^{\rm H}\cdot {\bf S}( \theta ) } \vert }^2} \right\}.$$

The second step is to estimate the range parameter based on the angle parameter information of the first step. The data model of the received signal can be written as

(21)$${\bf x} = \alpha _0{\bf S}_t( r_k, \;\hat{\theta }_k) + {\bf n} .$$

The weight vector of beamforming can be expressed as

(22)$${\bf w} = {\bf S}_t( r_k, \;\hat{\theta }_k) .$$

Therefore, the formula that can be used to solve the target parameter is

(23)$$\hat{r}_k = \arg \left\{{\mathop {\max }\limits_r {\vert {{\bf w}^{\rm H}\cdot {\bf x}} \vert }^2} \right\}.$$

Double-pulse localization scheme based on CPD

The data model of the double pulses signal of the lth snapshot can be written as

(24)$${\boldsymbol x}( l) = \left[\matrix{{\bf S}_{t0}( \theta_1) \otimes {\bf S}_t( r_1, \;\theta_1) , \;{\bf S}_{t0}( \theta_2) \otimes {\bf S}_t( r_2, \;\theta_2) , \;\hfill \cr \ldots , \;{\bf S}_{t0}( \theta_K) \otimes {\bf S}_t( r_K, \;\theta_K) \hfill} \right]{\boldsymbol s}( l) + {\boldsymbol n}( l) , \;$$

where ${\boldsymbol s}( l) = [ {s_1( l) , \;s_2( l) , \;\ldots , \;s_K( l) } ] ^T$ is composed of the amplitude and phase generated by the K targets in the lth pulse. The noise vector ${\bf n}( l)$ obeys the Gaussian distribution whose mean value is zero and the variance is σ ²I_NN. The entire array flow matrix can be further expressed as

(25)$${\bf S} = [ {{\bf S}_{t0}( \theta_1) \otimes {\bf S}_t( r_1, \;\theta_1) , \;\ldots , \;{\bf S}_{t0}( \theta_K) \otimes {\bf S}_t( r_K, \;\theta_K) } ] .$$

Considering that the total sampling times is L, the data matrix form outputted by matched filtering can be expressed as

(26)$$X = ( {{\bf S}_{t0}\odot {\bf S}_t} ) \boldsymbol{S}^T + N{\rm , \;}$$

where X = [x(1), x(2), …, x(L)] ∈ ℂ^NN×L, S = [s(1), s(2), …, s(L)] ∈ ℂ^L×K and noise matrix N ∈ ℂ^NN×L.

Since (26) satisfies the modulo n expansion form of the tensor third-order PARAFAC model, canonical polyadic decomposition (CPD) can be introduced to achieve target localization. The third-order tensor constructed by the array received data represented by (26) can be expressed as

(27)$$\gamma ( g, \;h, \;l) = \sum\limits_{k = 1}^K {S_{t0}( g, \;k) S_t( h, \;k) S( l, \;k) + \eta ( g, \;h, \;l) } , \;$$

where g = 1, 2, …, G, h = 1, 2, …, H, l = 1, 2, …, L. S _t0(g, k) represents the element at position (g, k) of the factor matrix S_t0, S _t(h, k) represents the element at position (h, k) of the factor matrix S_t, and S(l, k) represents the element at position (l, k) of the factor matrix $\boldsymbol{S}$. η is the third-order tensor representation of additive white noise.

In the first step, since there is a coupling of range and angle in the FDA's equivalent transmit steering vector, it is necessary to estimate the target angle parameter first. The CPD directly decomposes the tensor into the sum of each rank-1 tensor, which can use the structural information of the tensor to decompose it into the various factors. Considering that the kth column vector of the factor matrix S_t0 is $\varsigma _{tk}$, the phase obtained after normalization to eliminate scale scaling can be expressed as

(28)$$\mu _{tk} = angle( \varsigma _{tk}) = [ {0, \;\pi \sin \theta_k, \;\ldots , \;( G-1) \pi \sin \theta_k} ] ^T.$$

The weight vector of beamforming can be expressed as

(29)$${\bf w}_{t0} = \exp ( j\mu _{tk}) $$

Thus, the formula that can be used to estimate the angle parameter of the target can be written as

(30)$$\hat{\theta }_k = \arg \left\{{\mathop {\max }\limits_\theta {\vert {{\bf w}_{t0}^{\rm H} \cdot {\bf S}( \theta ) } \vert }^2} \right\}.$$

The second step is to combine the factor matrix S_t and μ _tk in the first step to estimate the range parameter. The kth column vector of the factor matrix S_t is ξ _tk. The phase obtained after normalization to eliminate scale scaling can be expressed as

(31a)$$\eqalignno{\nu _{tk} =\; & angle( \xi _{tk}) \cr =\; & [ {0, \;{{-4\Delta f\pi r_k} / c} + \pi \sin \theta_k, \;\ldots , \;{{-4( H-1) \Delta f\pi r_k} / c}} \cr & + ( G-1) \pi \sin \theta_k ] ^T,\;}$$

(31b)$$\eqalignno{\varepsilon _{tk} & = \nu _{tk}-\mu _{tk} \cr & = [ {0, \;{{-4\Delta f\pi r_k} / c}, \;\ldots , \;{{-4( H-1) \Delta f\pi r_k} / c}} ] ^T.}$$

Thus, the weight vector of beamforming can be obtained as

(32)$${\bf w}_t = \exp \{ {\,j\varepsilon_{tk}( r_k) } \} .$$

The formula for solving the target parameter can be written as

(33)$$\hat{r}_k = \arg \left\{{\mathop {\max }\limits_r {\vert {{\bf w}_t^{\rm H} \cdot {\bf S}( r, \;\theta ) } \vert }^2} \right\}.$$

Stepped frequency pulses localization scheme based on CPD

The time frequency schematic diagram of the stepped frequency pulses localization scheme is shown in Fig. 3. The total number of pulses is P. The frequency of the pth pulse of the nth element can be expressed as

(34)$$f_{n, k} = f_0 + n\Delta f + p\Delta f_t, \;\quad p = 0, \;1, \;\ldots , \;P-1, \;$$

where Δf _t is the unit increment of frequency between pulses.

Fig. 3.

Schematic diagram of stepped frequency pulse.

Thus, the transmitted signal of the nth element can be expressed as

(35)$$s_n( t) = \sum\limits_{\,p = 0}^{P-1} {rect\left({\displaystyle{{t-pT_p} \over {T_p}}} \right)\psi _n( t-pT_p) e^{\,j2\pi ( f_0 + n\Delta f + p\Delta f_t) ( t-pT_p) }} .$$

When observing a target point at (r, θ), the received signal related to the nth transmit array element and the mth receive array element can be expressed as

(36)$$ \eqalign{& s_{n, m}( t) = \sum\limits_{\,p = 0}^{P-1} {rect \bigg({\displaystyle{{t-pT_p-\tau_{n, m}} \over {T_p}}} \bigg) }\cr& \psi _n( t-pT_p-\tau _{n, m}) e^{\,j2\pi ( f_0 + n\Delta f + p\Delta f_t) ( t-pT_p-\tau _{n, m}) }} .$$

Thus, after signal processing and matched filtering at the receiver, the approximate signal can be expressed as

(37)$$ \eqalign{& x_{n, m}( t) = \cr& \sum\limits_{\,p = 0}^{P-1} {\xi _se^{\,j2\pi ( {-}f_0( 2r/c) -n\Delta f( 2r/c) -p\Delta f_t( 2r/c) \;+\; f_0( ( nd_T\sin \theta ) /c) \;+\; f_0( ( md_R\sin \theta ) /c) ) }}} .$$

The array factor can be expressed as

(38)$$\eqalign{AF = \;& \sum\limits_{m = 0}^{M-1} {\sum\limits_{n = 0}^{N-1} {x_{n, m}( t) } } \cr = \;&\xi _se^{{-}j2\pi f_0( 2r/c) } \sum\limits_{m = 0}^{M-1} {\sum\limits_{n = 0}^{N-1} {\sum\limits_{\,p = 0}^{P-1} }}\cr& {\xi _se^{\,j2\pi ( {-}n\Delta f( 2r/c) -p\Delta f_t( 2r/c) + f_0( ( nd_T\sin \theta ) /c) + f_0( ( md_R\sin \theta ) /c) ) .}} } $$

Thus, the steering vector of the entire array can be written as

(39)$${\bf S} = {\bf a}( r, \;\theta ) \otimes {\bf p}( r) \otimes {\bf b}( \theta ) , \;$$

where ${\bf a}( r, \;\theta )$, p(r), and ${\bf b}( \theta )$ can be denoted as

(40a)$$ \eqalign{ {\bf a}( r, \;\theta ) = & [ 1, \;\quad e^{\,j2\pi [ -\Delta f 2r/c + f_0 d_T\sin \theta /c ] }, \;\ldots , \; \cr& e^{\,j2\pi [ -( N-1) \Delta f 2r/c + ( N-1) f_0 d_T\sin \theta /c] }] ^{\rm T}, \;} $$

(40b)$${\bf p}( r) = [ 1, \;\quad e^{{-}j2\pi \Delta f_t( 2r/c) }, \;\ldots , \;e^{{-}j2\pi ( P-1) \Delta f_t( 2r/c) }] ^{\rm T}, \;$$

(40c)$${\bf b}( \theta ) = [ 1, \;\quad e^{\,j2\pi f_0( ( d_R\sin \theta ) /c) }, \;\ldots , \;e^{\,j2\pi ( M-1) f_0( ( d_R\sin \theta ) /c) }] ^{\rm T}.$$

The transmit steering vector can be expressed as

(41)$${\bf S}_t = {\bf a}( r, \;\theta ) \otimes {\bf p}( r) .$$

To ensure the fairness of the comparison, according to the analysis in [Reference Wang and Shao12], the equivalent transmit beampattern is still analyzed here. The data model of the received signal can be expressed as

(42)$${\bf x}_{r, \theta } = \alpha _0{\bf S}_t( r, \;\theta ) + {\bf n}.$$

After matched filtering, the data model of the lth snapshot equivalent transmitter signal can be rewritten as

(43)$$\eqalign{{\boldsymbol x}( l) =\; & [ {{\bf S}_t( r_1, \;\theta_1) , \;{\bf S}_t( r_2, \;\theta_2) , \;\ldots , \;{\bf S}_t( r_K, \;\theta_K) } ] {\boldsymbol s}( l) + {\boldsymbol n}( l) \cr =\; & \left[\matrix{{\bf a}( r_1, \;\theta_1) \otimes {\bf p}( r_1) , \;{\bf a}( r_2, \;\theta_2) \otimes {\bf p}( r_2) , \;\hfill \cr \ldots , \;{\bf a}( r_K, \;\theta_K) \otimes {\bf p}( r_K) \hfill} \right]{\boldsymbol s}( l) + {\boldsymbol n}( l) , \;} $$

where ${\bf s}( l) = [ {s_1( l) , \;s_2( l) , \;\ldots , \;s_K( l) } ] ^T$ is composed of the amplitude and phase generated by the K targets in the lth pulse. The noise vector ${\boldsymbol n}( l)$ obeys the Gaussian distribution whose mean value is zero and the variance is $\sigma ^2{\bf I}_{NN}$. Thus the array flow matrix can be further expressed as

(44)$${\bf S} = [ {{\bf a}( r_1, \;\theta_1) \otimes {\bf p}( r_1) , \;\ldots , \;{\bf a}( r_K, \;\theta_K) \otimes {\bf p}( r_K) } ] .$$

To introduce tensor operations, the above formula can be rewritten in the form of the Khatri–Rao product as

(45)$${\bf S} = {\bf a}\odot {\bf p} .$$

Considering that the total number of samples is L, the data matrix form outputted by matched filtering can be expressed as

(46)$${\boldsymbol X} = ( {{\bf a}\odot {\bf p}} ) {\boldsymbol S}^T + {\boldsymbol N}, \;$$

where ${\boldsymbol X} = [ {x( 1) , \;x( 2) , \;\ldots , \;x( L) } ] \in {\mathbb C}^{NN \times L}, \;{\boldsymbol S} = [ {s( 1) , \;s( 2) , \;\ldots , s( L) } ] \in {\mathbb C}^{L \times K}$ and noise matrix ${\boldsymbol N}\in {\mathbb C}^{NN \times L}$.

Using CPD to obtain localization of the targets, the third-order tensor generated from the data received by the array of (46) can be constructed as

(47)$$\varpi ( i, \;j, \;l) = \sum\limits_{k = 1}^K {a( i, \;k) p( j, \;k) S( l, \;k) + \eta ( i, \;j, \;l) } , \;$$

where i = 1, 2,.…, I, j = 1, 2, …, J, l = 1, 2, …, L. a(i, k) represents the element at position (i, k) of the factor matrix ${\bf a}$, p(j, k) represents the element at position (j, k) of the factor matrix ${\bf p}$, and S(l, k) represents the element at position (l, k) of the factor matrix ${\boldsymbol S}$. η is the third-order tensor representation of additive white noise.

In the first step, due to the coupling of range and angle equivalent transmit steering vector in FDA, it is necessary to estimate the target range parameter first. Considering that the kth column vector of the factor matrix ${\bf p}$ is ϕ _k, the phase obtained after normalization to eliminate scale scaling can be expressed as

(48)$$\eqalign{\varphi _k & = angle( \phi _k) \cr & = [ {0, \;{{-4\pi \Delta f_tr} / c}, \;\ldots , \;{{-4\pi ( J-1) \Delta f_tr} / c}} ] ^T.}$$

The weight vector of beamforming can be expressed as

(49)$${\bf w} = \exp ( j\varphi _k) .$$

Then the formula that can be used to estimate the range parameter of the target can be written as

(50)$$\hat{r}_k = \arg \left\{{\mathop {\max }\limits_r {\vert {{\bf w}^{\rm H}\cdot {\bf x}} \vert }^2} \right\}.$$

The second step is to estimate the angle parameter based on the factor matrix a and φ_k in the first step. The kth column vector of the factor matrix a is φ_k. Thus, the phase obtained after normalization to eliminate scale scaling can be expressed as

(51a)$$\eqalign{\zeta _k = \;& angle( \varphi _k) \cr = \;& [ {0, \;{{-4\Delta f\pi r_k} / c} + \pi \sin \theta_k, \;\ldots , \; }\cr& { {{-4( J-1) \Delta f\pi r_k} / c} + ( I-1) \pi \sin \theta_k} ] ^T, \;} $$

(51b)$$\delta _k = \zeta _k-\varphi _k = [ {0, \;\pi \sin \theta_k, \;\ldots , \;( I-1) \pi \sin \theta_k} ] ^T.$$

The weight vector of beamforming can be obtained as

(52)$${\bf w} = \exp \{ {\,j\delta_k( \theta_k) } \} .$$

The formula for solving the target parameter can be written as

(53)$$\hat{\theta }_k = \arg \left\{{\mathop {\max }\limits_\theta {\vert {{\bf w}^{\rm H}\cdot {\bf x}} \vert }^2} \right\}.$$

CPD-based localization methods have an additional step of tensor decomposition compared to the traditional method. Its complexity is mainly reflected in iteratively solving the three factor matrices in equations (27) and (47). Assuming that the number of signals is K, the total number of iterations is Q and the number of stepped frequency pulses is L. Then the total computational complexity of the three factor matrices after one iteration in two CPD-based methods are O(3K ³ + 3KN ²Q + K(N ² + 2NQ) + 2K ²(N ² + 2NQ)) and O(3K ³ + 3KNLQ + K(NL + NQ + LQ) + 2K ²(NL + NQ + LQ)), respectively. Finally, all three of the methods require a one-dimensional search. Assuming that the step size of the search is γ, the computational complexity of the search step is O(γ[N ²(N ² − K) + N ² − K]). It indicates that the effect of the variation in the number of array elements on the computational complexity of the CPD algorithm is extremely obvious, and the performance advantage of the method is derived from the higher computational complexity. Thus, it is important to explore the low computational complexity of the decomposition algorithm for the real-time implementation of the method in this paper.

CRLB

Equation (42) can be rewritten as

(54)$$\bf{x}_u = \sqrt {{\rm SNR}} \cdot {\bf u}( r, \;\theta ) + {\bf n}, \;$$

where SNR is the signal-to-noise power ratio. The mean and variance of the noise vector ${\bf n}$ denote all-zero vector and the identity matrix, respectively. When the double pulses FDA scheme in [Reference Wang and Shao12] is adopted, ${\bf u}( r, \;\theta )$ can be expressed as

(55)$${\bf u}_D( {r, \;\theta } ) = [ {{ {{\bf u}( {r, \;\theta } ) } \vert }_{\omega_i}, \;{ {{\bf u}( {r, \;\theta } ) } \vert }_{\omega_0}} ] ^T.$$

When the stepped frequency pulses FDA scheme is adopted, ${\bf u}( r, \;\theta )$ can be denoted as

(56)$${\bf u}_S( {r, \;\theta } ) = [ {{ {{\bf a}( {r, \;\theta } ) } \vert }_{\omega_0}, \;\ldots , \;{ {{\bf a}( {r, \;\theta } ) } \vert }_{\omega_{P-1}}} ] ^T, \;$$

where

(57a)$$ {{\bf u}( {r, \;\theta } ) } \vert _{\omega _i} = [ {1, \;e^{{-}j 2\pi f_0\Delta d\sin \theta /c }, \;\ldots , \;e^{{-}j 2\pi ( N-1) f_0\Delta d\sin \theta /c}} ] , \;$$

(57b)$$ {{\bf u}( {r, \;\theta } ) } \vert _{\omega _0} = \left[\matrix{1, \;e^{{-}j[ { 2\pi ( f_0 + \Delta f) \Delta d\sin \theta /c - 2\pi \Delta fr/c } ] }, \;\hfill \cr \ldots , \;e^{{-}j[ {( N-1) 2\pi ( f_0 + \Delta f) \Delta d\sin \theta /c -( N-1) 2\pi \Delta fr/c } ] } \hfill} \right], \;$$

(57c)$$ \eqalign{& {{\bf u}( {r, \;\theta } ) } \vert _{\omega _{P-1}} = \cr & \hskip-2pt\left[\matrix{1, \;e^{{-}j[ { { {2\pi [ {\,f_0 + \Delta f + ( P-1) \Delta f_t} ] \Delta d\sin \theta } /c} - 2\pi \Delta fr/c -( {P-1} ) 2\pi \Delta f_tr/c } ] }, \;\hfill \cr \ldots , \;e^{{-}j[ {( N-1) { {2\pi [ {\,f_0 + \Delta f + ( P-1) \Delta f_t} ] \Delta d\sin \theta } /c} -( N-1) 2\pi \Delta fr/c -( {P-1} ) 2\pi \Delta f_tr/c } ] } \hfill} \right].}$$

The mean ${\bf \mu }$ and variance ${\bf \Gamma }$ of the data model x_u can be expressed as

(58)$${\bf \mu } = {\bf u}( r, \;\theta ) \sqrt {{\rm SNR}} , \;$$

(59)$${\bf \Gamma } = {\bf I}.$$

Thus, the Fisher information matrix (FIM) can be obtained as

(60)$${\bf F} = 2{\mathop{\rm Re}\nolimits} \left[{\displaystyle{{d{\it \mu }^\ast } \over {d{\it \beta }}}{\bf \Gamma }^{{-}1}\displaystyle{{d{\it \mu }} \over {d{\it \beta }^T}}} \right], \;$$

where ${\mathop{\rm Re}\nolimits} [{\cdot} ]$ stands for real value operation, ${\it \beta } = [ {r, \;\theta } ] ^T$. Thus, the FIM of stepped frequency pulses FDA scheme can be written as

(61)$$\eqalign{{\bf F}_S = & 2\cdot SNR\cdot \cr & \left[{\matrix{ {\sum\limits_{\,p = 0}^{P-1} {\sum\limits_{n = 0}^{N-1} {[ {\vartheta_p^2 ( n ) } ] } } } \hfill & {\sum\limits_{\,p = 0}^{P-1} {\sum\limits_{n = 0}^{N-1} {[ {\vartheta_p( n ) \cdot \rho_p( n ) } ] } } } \hfill \cr {\sum\limits_{\,p = 0}^{P-1} {\sum\limits_{n = 0}^{N-1} {[ {\vartheta_p( n ) \cdot \rho_p( n ) } ] } } } \hfill & {\sum\limits_{\,p = 0}^{P-1} {\sum\limits_{n = 0}^{N-1} {[ {\rho_p^2 ( n ) } ] } } } \hfill \cr } } \right], \;}$$

where

(62a)$$\vartheta _p( n ) = n\displaystyle{{2\pi [ {\,f_0 + \Delta f + p\Delta f_t} ] \Delta d\cos \theta } \over c}, \;$$

(62b)$$\rho _p( n ) = {-}n\displaystyle{{2\pi \Delta fr} \over c}-p\displaystyle{{2\pi \Delta f_tr} \over c}.$$

Then the inverse matrix of the FIM can be expressed as

(63)$${\bf F}_S^{{-}1} = \displaystyle{1 \over {F_{1, 1}F_{2, 2}-F_{1, 2}F_{2, 1}}}\left[{\matrix{ {F_{2, 2}} & {-F_{2, 1}} \cr {-F_{1, 2}} & {F_{1, 1}} \cr } } \right].$$

Thus, the CRLB of the angle dimension and the range dimension can be solved, which are expressed as

(64a)$${\rm CRLB}_\theta = [ {\bf F}_S^{{-}1} ] _{1, 1} = \displaystyle{{F_{2, 2}} \over {F_{1, 1}F_{2, 2}-F_{1, 2}F_{2, 1}}}, \;$$

(64b)$${\rm CRLB}_r = [ {\bf F}_S^{{-}1} ] _{2, 2} = \displaystyle{{F_{1, 1}} \over {F_{1, 1}F_{2, 2}-F_{1, 2}F_{2, 1}}}, \;$$

where [ ⋅ ]_i,j represents the element in the ith row and jth column of the matrix.

Simulation results

In this section, the performance of proposed target localization scheme is verified by numerical simulation. Both single-target and dual-target situations are considered, and a comparison with the traditional double pulse FDA localization effects is provided.

Simulation parameters are considered as follows: the number of radar elements is N = 16, the reference frequency of the FDA transmit signal is f ₀ = 6 GHz, the frequency offset increment between elements is Δf = 2 kHz, and the frequency increment between pulses is Δf _t = 2 kHz. In the simulation, Gaussian white noise with the same variance and zero mean is used.

Detecting a single target is considered first. Assuming that the target's position information is $( r_1, \;\theta _1) = ( 80\,{\rm km}, \;0^\circ )$ under the condition of a signal-to-noise ratio (SNR) of −10 dB. The equivalent transmit beampattern of the traditional double-pulse FDA (DP-FDA) is shown in Fig. 4. It can be seen that the first step of DP-FDA is to use the phased array with Δf = 0 to estimate the angle parameters of the target. Since the transmit beampattern during Δf ≠ 0 is coupled with range and angle parameters, the angle parameter from the first step is used as the prior information to estimate the range parameters of the target.

Fig. 4.

Traditional double-pulse FDA equivalent transmit beampattern. (a) The first step. (b) The second step.

Figure 5 shows the equivalent transmit beampattern of CPD-based double-pulse FDA (CPD-DP-FDA) and CPD-based stepped frequency pulses FDA (CPD-SP-FDA). CPD-DP-FDA can decouple the range and angle of the equivalent transmit beampattern according to the data received by the radar. Thus, the maximum peak position of the beam can be used to locate the target. Since the radar's multi-carrier frequency time-division transmission provides a higher degree of freedom, the CPD-SP-FDA can form an equivalent beampattern of decoupling range and angle. Thus, the range parameter and angle parameter of the target can be estimated directly according to the maximum peak position of the beam.

Fig. 5.

FDA equivalent transmit beampattern based on CPD. (a) CPD-DP-FDA. (b) CPD-SP-FDA.

Figure 6 shows the range and angle dimension beampatterns for parameter estimation with the three schemes. It can be observed that in the range dimension, the mainlobe of the CPD-SP-FDA is narrower and the sidelobes are lower, which is more conducive to the estimation of range parameters. In the angle dimension, the tensor decomposition effectively extracts the target angle information in the double pulse data, which achieve the selectivity in the angle dimension. Thus, the mainlobe of range-angle decoupling is formed in the beampattern, which shows that CPD-DP-FDA has certain advantages in angle parameter estimation. By deriving CRLB, it can be found that the number of pulses P, the number of array elements N, the frequency offset increment between array elements Δf, and the frequency increment between pulses Δf _t will all affect the estimation accuracy of the range parameter.

Fig. 6.

Comparison of beam performance. (a) Range dimension. (b) Angle dimension.

Next, parameter estimation for multiple targets is considered. Assuming that the coordinates of the two targets are $( r_1, \;\theta _1) = ( 70\,{\rm km}, \;20^\circ )$ and $( r_2, \;\theta _2) = ( 80\,{\rm km}, \;0^\circ )$ respectively. Figures 7 and 8 show the equivalent transmit beampatterns of the three schemes when two targets are considered. Since the two targets have different angle parameters, two angle parameters will be obtained in the first step of the traditional double pulse processing, as shown in Fig. 7(a). In the second step of estimating the range parameter, it cannot be used as a priori information to estimate the range parameter directly because the angle dimension peak information is not unique, as shown in Fig. 7(b). In Fig. 8(a), since the factor matrix containing the target information after CPD is corresponding, CPD-DP-FDA realizes the automatic matching of the two target parameter estimates. Thus, there is no problem of ambiguity in parameter estimates. CPD-SP-FDA can also correspond to the positions of the two targets to the peak points, as shown in Fig. 8(b), which is conducive to parameter estimation of multiple targets.

Fig. 7.

Traditional double-pulse FDA equivalent transmit beampattern. (a) The first step. (b) The second step.

Fig. 8.

FDA equivalent transmit beampattern based on CPD. (a) CPD-DP-FDA. (b) CPD-SP-FDA.

DP-FDA will obtain two angle parameters $\hat{\theta }_1 = 0^\circ$ and $\hat{\theta }_2 = 20^\circ$ in the first step, and obtain four target points $( 45\,{\rm km}, \;0^\circ ) , \;( 70\,{\rm km}, \;20^\circ ) , \;( 80\,{\rm km}, \;0^\circ ) , \;( 105\,{\rm km}, \;20^\circ )$ in the subsequent parameter estimation in the second step. Thus, neither the number of real targets nor the position information under this scheme can be effectively distinguished. When CPD-DP-FDA scheme is used, automatic matching of target parameters is achieved during the solution process so that the parameters of the two targets can be estimated. CPD-SP-FDA realizes the decoupling of the beampattern by time-division transmission. Thus, it can also estimate the parameters of the two targets. Figure 9 is a schematic diagram of the estimation effect of the first target. It can be found that the two CPD-based schemes have a unique mainlobe in both the range dimension and the angle dimension, which can directly estimate the range and angle parameters of the first target.

Fig. 9.

Performance comparison of the first target beampattern. (a) Range dimension $( \theta = 20^\circ )$. (b) Angle dimension (r = 70 km).

Figure 10 shows a schematic diagram of the estimation effect of the second target. Similarly, the two CPD-based schemes have a unique mainlobe in the range dimension and the angle dimension, which have lower sidelobe levels. Thus, CPD-DP-FDA can effectively solve the problem of unclear target indication in DP-FDA. Then CPD-SP-FDA can also effectively estimate the range and angle parameters of multiple targets. Moreover, the half-power (−3 dB) contours of the equivalent target indication patterns corresponding to DP-FDA, CPD-DP-FDA, and CPD-SP-FDA are compared in Fig. 11. It is clearly shown that the scheme based on tensor decomposition can effectively solve the problem of uncertain multi-target parameter estimation, and the improvement of target parameter estimation performance brought by the stepped frequency pulse is also indicated. Table 1 shows the performance of beampattern in range and angle dimensions with different FDA methods. It can be seen from Table 1 that CPD-DP-FDA has the same BW with DP-FDA due to the same parameter while it can accomplish the discrimination of multiple targets. The BW _r of CPD-SP-FDA is 3.1 km, which is much better than that of other methods. In addition, due to the complex solution of factor matrices, the solution time of the CPD-based methods also grow.

Fig. 10.

Performance comparison of the second target beampattern. (a) Range dimension $( \theta = 0^\circ )$. (b) Angle dimension (r = 80 km).

Fig. 11.

Comparison of the −3 dB profile indicated by the equivalent beampattern.

Table 1.

Comparison of different FDA methods

Finally, Fig. 12 shows the CRLB performance comparison between the range dimension and the angle dimension. Because the traditional DP-FDA has ambiguity in multi-target parameter estimation, the two CPD-based schemes are compared in the figure. It can be seen that the parameter estimation effect of CPD-SP-FDA in the range dimension and the angle dimension is better than that of CPD-DP-FDA. Considering that SP-FDA uses two pulses, the performance of CRLB is still better than DP-FDA due to the higher degree of freedom provided by multiple carrier frequencies. As the number of pulses increases to 8, the selectivity of SP-FDA in the range dimension is further improved, and CRLB performance is further optimized. However, it should be noted that the increase in the number of pulses increases the computational complexity of tensor decomposition, which means that the number of pulses needs to be reasonably selected according to the calculation ability in the actual application of the solution.

Fig. 12.

Performance comparison of CRLB. (a) Range dimension. (b) Angle dimension.

Although the proposed method is superior to other methods, the paper only considers scenarios with static targets. It means that extending the model to moving target scenarios requires the incorporation of new techniques. The introduction of cognitive radar ideas into FDA radar has recently attracted scholarly attention [Reference Wang33–Reference Basit, Wang and Nusenu36]. The cognitive-based processing model enables the detection and tracking of dynamic targets. Thus, it provides a reference for extending the model to dynamic target scenarios, which is the next research focus of system applicability enhancement.