**3. Proposed Algorithm**

The basic premise of the ESPRIT algorithm is that there are identical subarrays, the spacing between subarrays is known and the structure of subarrays is identical, which satisfies the rotational invariance in space [13]. Uniform linear arrays (ULAs) appear when it comes to one-dimensional DOA estimation using conventional ESPRIT [1,13]. Compared with ULAs, double parallel uniform linear arrays (DPULAs) can identify two-dimensional DOA because of the special construction, which consists of two parallel ULAs [36–38]; therefore, the algorithm discussed in this paper is based on DPULAs.

#### *3.1. Complex Representation Matrix and Related Calculations*

In view of the paucity of research on calculations with multi-vector, the Complex Representation Matrix (CRM) [20] is introduced because of the mature matrix theories. Consider a matrix **A** ∈ G*m*×*n* 3 , the CRM is defined by Ψ(**A**)

$$\Psi(\mathbf{A}) = \begin{pmatrix} \mathbf{A}\_0 + \mathbf{A}\_3 + (\mathbf{A}\_7 + \mathbf{A}\_4)\mathbf{e}\_{123} & -\mathbf{A}\_1 + \mathbf{A}\_6 + (\mathbf{A}\_2 - \mathbf{A}\_5)\mathbf{e}\_{123} \\ -\mathbf{A}\_1 - \mathbf{A}\_6 - (\mathbf{A}\_2 + \mathbf{A}\_5)\mathbf{e}\_{123} & \mathbf{A}\_0 - \mathbf{A}\_3 + (\mathbf{A}\_7 - \mathbf{A}\_4)\mathbf{e}\_{123} \end{pmatrix}. \tag{8}$$

Let *ν* = (−**e**<sup>1</sup> + **<sup>e</sup>**13)/2 ∈ G3, and its reversion is *ν*˜ = (−**e**<sup>1</sup> − **<sup>e</sup>**13)/2 ∈ G3. Then,

$$
\nu^2 = \vec{\nu}^2 = 0 \quad \text{and} \quad \nu \vec{\nu} + \vec{\nu}\nu = 1,\tag{9}
$$

which imply *ννν*˜ = *ν*, *νν*˜ *ν*˜ = *ν*˜, (*νν*˜)<sup>2</sup> = *νν*˜, (*νν*˜ )2 = *νν*˜ .

It immediately follows that, for every **A** ∈ G3, we have

$$\mathbf{A} = \mathbf{E}\_{2m} \Psi(\mathbf{A}) \mathbf{E}\_{2n'}^H \tag{10}$$

$$\Psi(\mathbf{A}) = \mathbf{Q}\_{2m} \begin{bmatrix} \mathbf{A} & \mathbf{0} \\ \mathbf{0} & \mathbf{A} \end{bmatrix} \mathbf{Q}\_{2n\prime} \tag{11}$$

where in (10) and (11) we have

$$\mathbf{E}\_{2k} = \begin{bmatrix} \nu \boldsymbol{\upnu} \mathbf{I}\_k & \nu \mathbf{I}\_k \end{bmatrix} \in \mathbb{G}\_3^{k \times 2k},\tag{12}$$

$$\mathbf{Q}\_{2k} = \begin{bmatrix} \ \nu \vec{\upsilon} \mathbf{I}\_k & \vec{\upsilon} \mathbf{I}\_k \\ \ \nu \mathbf{I}\_k & \vec{\upsilon} \boldsymbol{\upsilon} \mathbf{I}\_k \end{bmatrix} \in \mathbb{G}\_3^{2k \times 2k}. \tag{13}$$

**I***k* denotes the *k* × *k* identity matrix. It is not difficult to prove that

$$\mathbf{Q}\_{2k} = \mathbf{Q}\_{2k}^H = \mathbf{Q}\_{2k}{}^{-1},\tag{14a}$$

$$\Psi(\mathbf{A}^H) = \left(\Psi(\mathbf{A})\right)^H,\tag{14b}$$

$$\Psi(\mathbf{A}^+) = (\Psi(\mathbf{A}))^+,\tag{14c}$$

where {+} denotes the pseudo-inverse. Referring to (10) and (14c), the pseudo-inverse of any **A** ∈ G3 is

$$\mathbf{A}^{+} = \mathbf{E}\_{2u} (\mathbf{\varPsi}(\mathbf{A}))^{+} \mathbf{E}\_{2u}^{H}. \tag{15}$$

Since **e**2123 = −1 and **e**123 commutes with all elements in G3, one can identify it with the complex imaginary unit *j* [20], and so we can view Ψ(**A**) given in (8) as a complex matrix.

#### *3.2. Model for DPULAs*

Consider a DPULA with 2*M* + 2 sensors, as shown in Figure 2, in which *d* and *M* refer to the spacing between two adjacent sensors and the number of sensors in per subarray, respectively. The array is divided into three subarrays. The 1st to *M*th sensors on the *x*-axis compose the first subarray, the 2nd to (*M* + 1)th sensors form the second subarray and the (*M* + 2)th to (2*M* + 1)th that located on a straight line parallel to the *x*-axis make up the third subarray. The reason for the division can be found in Figure 3, that is, there are two unknown DOA parameters in the model, which need two rotational invariance relations.

**Figure 2.** Double parallel uniform linear array.

**Figure 3.** Schematic diagram of GA-ESPRIT.

Since the three subarrays have the same structure and the same number of sensor, each output of them has only one phase difference for the same signal. Signals received by subarray one, two and three are defined as **<sup>Y</sup>**1*EH*, **<sup>Y</sup>**2*EH* and **<sup>Y</sup>**3*EH*, respectively. According to the above array model, the outputs of the three subarrays at time *t* are as follows

$$\begin{aligned} \mathbf{Y}\_{EH}^1(t) &= \mathbf{AS}(t) + \mathbf{N}^1(t), \\ \mathbf{Y}\_{EH}^2(t) &= \mathbf{AFS}(t) + \mathbf{N}^2(t), \\ \mathbf{Y}\_{EH}^3(t) &= \mathbf{AGS}(t) + \mathbf{N}^3(t), \end{aligned} \tag{16}$$

where

$$\begin{aligned} \mathbf{Y}\_{EH}^{1}(t) &= \left[ \mathbf{Y}\_{EH}^{(1)}(t), \dots, \mathbf{Y}\_{EH}^{(M)}(t) \right]^T, \\ \mathbf{Y}\_{EH}^{2}(t) &= \left[ \mathbf{Y}\_{EH}^{(2)}(t), \dots, \mathbf{Y}\_{EH}^{(M+1)}(t) \right]^T, \\ \mathbf{Y}\_{EH}^{3}(t) &= \left[ \mathbf{Y}\_{EH}^{(M+2)}(t), \dots, \mathbf{Y}\_{EH}^{(2M+1)}(t) \right]^T, \end{aligned} \tag{17}$$

and

$$\begin{aligned} \mathbf{A} &= [a(\Gamma\_1), \dots, a(\Gamma\_K)], \\ a(\Gamma\_k) &= \left[1, \mathbf{x}(\theta\_k, \phi\_k), \dots, \mathbf{x}^{M-1}(\theta\_k, \phi\_k)\right]^T V\_k P\_{\mathbf{k}}, \\ \mathbf{x}(\theta\_{k'} \phi\_k) &= \epsilon^{\mathbf{r} \times 23} \frac{2\pi}{\lambda} \text{diag}\,\theta\_k \sin\phi\_{k'} \\ \mathbf{F} &= \text{diag}(f\_1, \dots, f\_K), \quad \mathbf{G} = \text{diag}(g\_1, \dots, g\_K). \end{aligned} \tag{18}$$

According to (18), we find that the DOA information is contained in matrix **A**, **F** and **G**. Because **F** and **G** are diagonal matrices that only contain direction information of incident signals, the focus is the two matrices, i.e.,

$$\begin{cases} f\_k = e^{\mathfrak{e}\_{123} \frac{2\mathfrak{g}}{\Lambda} \text{d} \cos \theta\_k \sin \phi\_{k\_f}}\\ g\_k = e^{\mathfrak{e}\_{123} \frac{2\mathfrak{g}}{\Lambda} \text{d} \sin \theta\_k \sin \phi\_k} \end{cases} \tag{19}$$

Clearly, it is easy to figure out the DOA in the light of (19) if we obtain the two ideal matrices **F** and **G**. From the rules of subarray division, we can see that the latter (*M* − 1) sensors of subarray one and the former (*M* + 1) sensors of subarray two are overlapped. Thus, in order to reduce the computational complexity, subarray one and subarray two can be merged to form a new matrix **Y***EH*, that is,

$$\mathbf{Y}\_{EH}(t) = \begin{bmatrix} y\_1(t), y\_2(t), \dots, y\_{M+1}(t) \end{bmatrix}^T. \tag{20}$$

After merging, the (2*M* + 2)*th* redundant sensor is added to subarray three to form a new subarray **P***EH*, so that the third subarray has the same dimension as **Y***EH*

$$\mathbf{P}\_{EH}(t) = \begin{bmatrix} y\_{M+2}(t), y\_{M+3}(t), \dots, y\_{2M+2}(t) \end{bmatrix}^T. \tag{21}$$

Let **A** be the array flow pattern of **Y***EH*, then

$$\begin{aligned} \mathbf{A} &= \left[ \overline{\mathfrak{a}}(\Gamma\_1), \overline{\mathfrak{a}}(\Gamma\_2), \dots, \overline{\mathfrak{a}}(\Gamma\_K) \right]\_\prime \\ \overline{\mathfrak{a}}(\Gamma\_k) &= \left[ 1, \mathbf{x}(\theta\_k, \phi\_k), \dots, \mathbf{x}^M(\theta\_k, \phi\_k) \right]^T V\_k P\_k. \end{aligned} \tag{22}$$

**Y***EH* and **P***EH* can be written as

$$\begin{aligned} \mathbf{Y}\_{EH}(t) &= \overline{\mathbf{A}} \mathbf{S}(t) + \mathbf{N}\_d(t), \\ \mathbf{P}\_{EH}(t) &= \overline{\mathbf{A}} \mathbf{G} \mathbf{S}(t) + \mathbf{N}\_b(t), \end{aligned} \tag{23}$$

where

$$\mathbf{N}\_a(t) = \begin{bmatrix} \mathbf{N}^1(t) \\ \boldsymbol{\eta}\_{M+1}(t) \end{bmatrix}, \quad \mathbf{N}\_b(t) = \begin{bmatrix} \mathbf{N}^3(t) \\ \boldsymbol{\eta}\_{2M+2}(t) \end{bmatrix}.$$

Then, **B**(*t*) is defined as

$$\mathbf{B}(t) = \begin{bmatrix} \mathbf{Y}\_{EH}(t) \\ \mathbf{P}\_{EH}(t) \end{bmatrix} = \mathbf{C}\mathbf{S}(t) + \mathbf{N}(t), \tag{24}$$

where

$$\mathbf{C} = \left[ \begin{array}{cc} \overline{\mathbf{A}} \\ \overline{\mathbf{A}} \mathbf{G} \end{array} \right], \quad \mathbf{N}(t) = \left[ \begin{array}{cc} \mathbf{N}\_a(t) \\ \mathbf{N}\_b(t) \end{array} \right].$$

Finally, the output of the whole array is denoted by

$$\mathbf{B}(t) = \mathbf{C}\mathbf{S}(t) + \mathbf{N}(t). \tag{25}$$
