A Computationally Efficient Estimation Algorithm for Direction of Arrival in Double Parallel Linear Array

Direction of arrival (DOA) estimation, a key issue in array signal processing, is common in various fields, namely, radar [1-3], sonar [4, 5], communications [6, 7], seismic survey [8], and radio astronomy [9]. Array antennas have been adopted in many methods to estimate the DOA, including autoregressive-moving average (ARMA) spectral analysis, maximum likelihood method, entropy spectral analysis, and eigenvalue decomposition.

Two subspace algorithms have been developed based on eigenvalue decomposition (EVD). In 1979, Schmidt proposed the multiple signal classification (MUSIC) algorithm [10-12], marking a breakthrough in DOA estimation theories. In 1986, Roy put forward the estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [13, 14]. Based on rotational invariance of subspaces, ESPRIT algorithm eliminates the need of global search and thus reduces the computing load. The key of subspace algorithms is to obtain signal subspace or noise subspace through EVD or singular value decomposition (SVD). Either EVD or SVD incurs a huge computing load, especially in massive multiple-input and multiple-output (MIMO) systems.

To simplify subspace algorithms, many excellent methods have been designed in recent years. For example, T.A.H. Bressner et al. [15] created the single snapshot algorithm for large arrays, which divides a large array into multiple subarrays for DOA estimation. Through low-rank matrix reconstruction and Vandermonde decomposition, Tian et al. [16] designed a gridless two-dimensional (2D) DOA estimation method that improves the running speed by avoiding EVD or SVD. But the applicable scope of this method is limited to uniform and sparse rectangular arrays. For 2D DOA estimation with L-shaped array, Nie and Li [17] proposed a computationally efficient subspace algorithm (CESA), which estimates 2D DOA efficiently with three vectors made from the elements in the first column, first row and diagonal of the covariance matrix, respectively. Since the two axes of the L-shaped array are composed of uniform linear arrays (ULAs), Yan et al. [18] established a CESA-based joint cross-covariance matrix (JCCM) algorithm for one-dimensional (1D) DOA estimation with ULA. Eliminating subspace decomposition, the JCCM algorithm sets up a JCCM based on the forward and backward cross-covariance matrices of input data, reconstructs the equivalent signal subspace, and estimates the DOA by solving polynomials

Drawing on ESPRIT algorithm [19], CESA algorithm, and JCCM algorithm, this paper proposes a 1D DOA estimation algorithm with double parallel linear array (DPLA). The DPLA [20] consists of two parallel ULAs with a known spacing. The DPLA can be viewed as an array of array elements pairs, each of which contains two sensors with identical response features. In the proposed algorithm, the noise subspace is constructed by processing the first column elements of the JCCM, and the DOA of signal is obtained by solving polynomials.

The remainder of this paper is organized as follows: Section II introduces the structure and model of the DPLA; Section III explains the DPLA algorithm and analyzes its performance; Section IV compares the DPLA algorithm with others through MATLAB simulation; Section V puts forward the conclusions.

3.1 DOA estimation

To mitigate the effect of additive noise, the JCCM between signal vectors ${{x}_{1}}\left( t \right)$ can be defined as:

$\begin{align}  & {{R}_{12}}=E\{{{x}_{1}}\left( t \right){{x}_{2}}^{H}\left( t \right)\} \\ & =A\left( \theta  \right)E\{s\left( t \right){{s}^{H}}\left( t \right)\}{{\psi }^{H}}{{A}^{H}}\left( \theta  \right)+E\{{{n}_{1}}\left( t \right){{n}_{2}}^{H}\left( t \right)\} \\ & =A\left( \theta  \right){{R}_{s}}{{\psi }^{H}}{{A}^{H}}\left( \theta  \right) \end{align}$     (6)

Similarly, the JCCM between signal vectors ${{x}_{1}}\left( t \right)$ and ${{x}_{2}}\left( t \right)$ can be defined as:

$\begin{align}  & {{R}_{21}}=E\{{{x}_{2}}\left( t \right){{x}_{1}}^{H}\left( t \right)\} \\ & =A\left( \theta  \right)\psi E\{s\left( t \right){{s}^{H}}\left( t \right)\}{{A}^{H}}\left( \theta  \right)+E\{{{n}_{2}}\left( t \right){{n}_{1}}^{H}\left( t \right)\} \\ & =A\left( \theta  \right)\psi {{R}_{s}}{{A}^{H}}\left( \theta  \right) \end{align}$     (7)

Since ${{n}_{1}}\left( t \right)$ and ${{n}_{2}}\left( t \right)$ are spatially independent and ${{s}_{k}}\left( t \right)k=1,2,\ldots ,K$ are uncorrelated, we have $E\{{{n}_{1}}\left( t \right){{n}_{2}}^{H}\left( t \right)\text{ }\!\!\}\!\!\text{ =}E\{{{n}_{2}}\left( t \right){{n}_{1}}^{H}\left( t \right)\text{ }\!\!\}\!\!\text{ =}0$,  ${{R}_{s}}=E\{s\left( t \right){{s}^{H}}\left( t \right)\}$ and ${{R}_{s}}=diag\left( \begin{matrix}   {{p}_{1}} & {{p}_{2}} & ... & {{p}_{K}}  \\\end{matrix} \right)$, with ${{p}_{i}}$$\left( i=1,2,...,K \right)$ being the power of the k-th source signal. The power of every source signal is a real number. Then, we have  $\forall {{p}_{i}}>0$ and ${{R}_{s}}={{R}_{s}}^{*}$.

According to the theory on matrix conjugation, ${{R}_{21}}$ is the conjugate of ${{R}_{12}}$. Then, matrix  $R$ can be defined as:

$\begin{align}  & R={{R}_{12}}+{{R}_{21}} \\ & =A\left( \theta  \right){{R}_{S}}\left( {{\psi }^{H}}+\psi  \right){{A}^{H}}\left( \theta  \right) \end{align}$     (8)

$\begin{align}  & {{R}_{S}}\left( {{\psi }^{H}}+\psi  \right)=\left( \begin{matrix}   {{p}_{1}} & {} & {}  \\   {} & \ddots  & {}  \\   {} & {} & {{p}_{k}}  \\\end{matrix} \right)\left( \begin{matrix}   {{e}^{-j{{\mu }_{1}}}}+{{e}^{j{{\mu }_{1}}}} & {} & {}  \\   {} & \ddots  & {}  \\   {} & {} & {{e}^{-j{{\mu }_{k}}}}+{{e}^{j{{\mu }_{K}}}}  \\\end{matrix} \right) \\ & =\left( \begin{matrix}   {{\alpha }_{1}} & {} & {}  \\   {} & \ddots  & {}  \\   {} & {} & {{\alpha }_{K}}  \\\end{matrix} \right)=\Lambda  \end{align}$     (9)

where, ${{\alpha }_{i}}={{p}_{i}}\left( {{e}^{-j{{\mu }_{i}}}}+{{e}^{j{{\mu }_{i}}}} \right),(i=1,2,...,K)$; $\forall {{\alpha }_{i}}>0$; $\Lambda \text{=}{{\Lambda }^{*}}$. Substituting formula (9) into formula (8), we have:

$R=A\left( \theta  \right)\Lambda {{A}^{H}}\left( \theta  \right)$     (10)

Since the first column elements of $A^{H}(\theta)$ are all ones, the $M \times 1$ vectors $r(:, 1)$ for the first column elements of $R$ can be expressed as $r(:, 1)=A(\theta) \Lambda 1,$ where 1 is a $K \times 1$ dimensional full 1 column vector. Based on $r(:, 1),$ a new column vector can be defined as $r^{n}(:, 1)=\operatorname{Jr}^{*}(:, 1),$ where

$J\text{=}\left( \begin{matrix}   {} & {} & 1  \\   {} &  & {}  \\   1 & {} & {}  \\\end{matrix} \right)$     (11)

Then, a $2M\times 1$ vector $r$ can be constructed from $r\left( :,1 \right)$ and ${{r}^{n}}\left( :,1 \right)$ :

$r=\left[ \begin{matrix}   {{r}^{n}}\left( :,1 \right)  \\   r\left( :,1 \right)  \\\end{matrix} \right]=\left[ \begin{matrix}   J{{A}^{*}}\left( \theta  \right)  \\   A\left( \theta  \right)  \\\end{matrix} \right]\Lambda 1$     (12)

Since the first row elements of $A\left( \theta  \right)$ are all ones, the last row elements of $J{{A}^{*}}\left( \theta  \right)$ must be all ones. Removing the first row of $A\left( \theta  \right)$ or the last row of $J{{A}^{*}}\left( \theta  \right)$, a $\left( 2M-1 \right)\times 1$ vector ${{r}_{new}}$ can be obtained:

$\begin{align}  & {{r}_{new}}=\left( \begin{matrix}   {{e}^{{j2\pi \left( M-1 \right)d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{j2\pi \left( M-1 \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\   \begin{matrix}   \vdots   \\   \begin{matrix}   1  \\   \vdots   \\\end{matrix}  \\\end{matrix} & \begin{matrix}   \cdots   \\   \begin{matrix}   \cdots   \\   \cdots   \\\end{matrix}  \\\end{matrix} & \begin{matrix}   \vdots   \\   \begin{matrix}   1  \\   \vdots   \\\end{matrix}  \\\end{matrix}  \\   {{e}^{{-j2\pi \left( M-1 \right)d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{-j2\pi \left( M-1 \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\\end{matrix} \right)\left[ \begin{matrix}   {{\alpha }_{1}}  \\   {{\alpha }_{2}}  \\   \vdots   \\   {{\alpha }_{K}}  \\\end{matrix} \right] \\ & \text{=}\overline{A\left( \theta  \right)}\alpha  \end{align}$     (13)

where

$\overline{A\left( \theta  \right)}=\left( \begin{matrix}   {{e}^{{j2\pi \left( M-1 \right)d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{j2\pi \left( M-1 \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\   \begin{matrix}   \vdots   \\   \begin{matrix}   1  \\   \vdots   \\\end{matrix}  \\\end{matrix} & \begin{matrix}   \cdots   \\   \begin{matrix}   \cdots   \\   \cdots   \\\end{matrix}  \\\end{matrix} & \begin{matrix}   \vdots   \\   \begin{matrix}   1  \\   \vdots   \\\end{matrix}  \\\end{matrix}  \\   {{e}^{{-j2\pi \left( M-1 \right)d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{-j2\pi \left( M-1 \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\\end{matrix} \right)$$\alpha \text{=}{{\left[ \begin{matrix}   {{\alpha }_{1}} & {{\alpha }_{2}} & \cdots  & {{\alpha }_{K}}  \\\end{matrix} \right]}^{T}}$

Obviously, ${{r}_{new}}$ is the data received by $2M-1$ elements through a single snapshot. On this basis, a $\left( 2M-K \right)\times K$ matrix $S$ can be constructed:

$S=\left[ \begin{matrix}   {{r}_{1}} & {{r}_{2}} & \cdots  & {{r}_{k}}  \\\end{matrix} \right]$     (14)

where, ${{r}_{i}}={{r}_{new}}\left( i:2M-K-1+i \right),(i=1,2,...,K)$ , ${{r}_{new}}\left( u:v \right)$is the column vector of the $u\text{th}$-$v\text{th}$ row of ${{r}_{new}}$. According to formulas (13) and (14), matrix $S$ can be rewritten as

$S=G\Lambda Q$     (15)

where

$G=\left( \begin{matrix}   {{e}^{{j2\pi \left( M-1 \right)d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{j2\pi \left( M-1 \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\   \begin{matrix}   \vdots   \\   \begin{matrix}   1  \\   \vdots   \\\end{matrix}  \\\end{matrix} & \begin{matrix}   \cdots   \\   \begin{matrix}   \cdots   \\   \cdots   \\\end{matrix}  \\\end{matrix} & \begin{matrix}   \vdots   \\   \begin{matrix}   1  \\   \vdots   \\\end{matrix}  \\\end{matrix}  \\   {{e}^{{-j2\pi \left( M-K \right)d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{-j2\pi \left( M-K \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\\end{matrix} \right)$

$Q=\left[ \begin{matrix}   1 & {{e}^{{-j2\pi d\sin {{\theta }_{1}}}/{\lambda }\;}} & \cdots  & {{e}^{{-j2\pi \left( K-1 \right)d\sin {{\theta }_{1}}}/{\lambda }\;}}  \\   \vdots  & \cdots  & \cdots  & \vdots   \\   1 & {{e}^{{-j2\pi d\sin {{\theta }_{K}}}/{\lambda }\;}} & \cdots  & {{e}^{{-j2\pi \left( K-1 \right)d\sin {{\theta }_{K}}}/{\lambda }\;}}  \\\end{matrix} \right]$

The relevant definitions are provided as follows [22]:

• If $A=\left[ {{\alpha }_{1}},{{\alpha }_{2}},\cdots ,{{\alpha }_{n}} \right]\in {{C}^{m\times n}}$ is a complex matrix, then the subspace $Col(A)$of all linear combinations of column vectors is called the column space or column span of matrix $A$:

$\begin{align}  & Col(A)=Span\{{{\alpha }_{1}},{{\alpha }_{2}},\cdots ,{{\alpha }_{n}}\} \\ & =\{y\in {{C}^{m}}:y=\sum\limits_{j=1}^{n}{{{x}_{j}}{{\alpha }_{j}}}:{{x}_{j}}\in C\} \end{align}$

The column space matrix $A$ is often represented by $Span\{A\}$:

$Col(A)=Span(A)=Span\{{{a}_{1}},{{a}_{2}},\cdots ,{{a}_{n}}\}$

• If $A$ is a complex matrix, the range of $A$ can be defined as$Range(A)=\{y\in {{C}^{m}}:Ax=y,x\in {{C}^{n}}\}$.

If $A=\left[ {{\alpha }_{1}},{{\alpha }_{2}},\cdots ,{{\alpha }_{n}} \right]$ is a column of $A$, and if $x={{\left[ {{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}} \right]}^{T}}$, then $Ax=\sum\limits_{j=1}^{n}{{{x}_{j}}{{\alpha }_{j}}}$. Thus, we have:

$\begin{align}  & Range(A)=\{y\in {{C}^{m}}:y=\sum\limits_{j=1}^{n}{{{x}_{j}}{{\alpha }_{j}}}:{{x}_{j}}\in C\} \\ & =Span\{{{a}_{1}},{{a}_{2}},\cdots ,{{a}_{n}}\} \end{align}$

This means the range of matrix $A$ is the column space of $A$:

$Range(A)=Col(A)=Span\{{{a}_{1}},{{a}_{2}},\cdots ,{{a}_{n}}\}$

• For matrix ${{A}_{m\times n}}$, the rank of $A$ has the following correlation with its column space and range:

$rank(A)=\dim[Col(A)]=\dim[Range(A)]$

Obviously, $G$ and $Q$ are Vandermonde matrices, and the source signals are uncorrelated. Therefore, the ranks of $G$, $Q$ and $\Lambda $ are all equal to the source number $K$. Then, we have $rank(G)=rank(S)=K$. Thus, $G$ and $S$ have $K$ linear independent columns, span the same K-dimensional space, and cover the same range: $Range(G)=Range(S)$. According to the relationship between column space and range, we have  $Span\{G\}=Span\{S\}$.

As a result, matrix $S$ contains all DOA information of source signals. Then, $\overline{S}$ can be obtained through Schmidt orthogonalization of the column vector of $S$.

Inspired by the root-MUSIC algorithm [23], it is assumed that $z={{e}^{jw}}$. Then, a polynomial can be constructed as:

$f\left( z \right)\text{=}{{z}^{2N-K}}{{p}^{H}}\left( z \right){{U}_{N}}p\left( z \right)$     (16)

where, $p\left( z \right)={{\left[ \begin{matrix}   1 & z & \cdots  & {{z}^{2N-K-1}}  \\\end{matrix} \right]}^{T}}$; ${{U}_{N}}={{I}_{2N-K}}-\overline{S}{{\overline{S}}^{H}}$.

Solving the polynomial, the $K$ roots closest to the unit circle can be obtained. Then, the DOA of the signal can be calculated by:

${{\theta }_{i}}=\arcsin \left( \frac{\lambda }{2\pi d}\arg \left\{ {{z}_{i}} \right\} \right)i=1,2,...,K$     (17)

In actual applications, the theoretical covariance matrices can be replaced with sample covariance matrices:

${{\overset{\scriptscriptstyle\frown}{R}}_{12}}=\frac{1}{N}\sum\limits_{t=1}^{N}{{{x}_{1}}}\left( t \right){{x}_{2}}^{H}\left( t \right)={{{X}_{1}}{{X}_{2}}^{H}}/{N}\;$     (18)

${{\overset{\scriptscriptstyle\frown}{R}}_{21}}=\frac{1}{N}\sum\limits_{t=1}^{N}{{{x}_{2}}}\left( t \right){{x}_{1}}^{H}\left( t \right)\text{=}{{{X}_{2}}{{X}_{1}}^{H}}/{N}\;$     (19)

$\overset{\scriptscriptstyle\frown}{R}\text{=}{{\overset{\scriptscriptstyle\frown}{R}}_{12}}+{{\overset{\scriptscriptstyle\frown}{R}}_{21}}$     (20)

where, ${{X}_{1}}=\left[ \begin{matrix}   {{x}_{1}}\left( 1 \right) & {{x}_{1}}\left( 2 \right) & \cdots  & {{x}_{1}}\left( N \right)  \\\end{matrix} \right]$ and ${{X}_{2}}=\left[ \begin{matrix}   {{x}_{2}}\left( 1 \right) & {{x}_{2}}\left( 2 \right) & \cdots  & {{x}_{2}}\left( N \right)  \\\end{matrix} \right]$ are $N$ snapshots of the signals outputted from subarrays 1 and 2, respectively.

Since only the first column of matrix $R$ is involved, the DPLA algorithm can be further simplified to reduce computing load. According to formula (20), the first column of $\overset{\scriptscriptstyle\frown}{R}$ equals the sum of the first column ${{\overset{\scriptscriptstyle\frown}{R}}_{12}}\left( :,1 \right)$ of ${{\overset{\scriptscriptstyle\frown}{R}}_{12}}$ and the first column ${{\overset{\scriptscriptstyle\frown}{R}}_{21}}\left( :,1 \right)$ of ${{\overset{\scriptscriptstyle\frown}{R}}_{21}}$:

${{\overset{\scriptscriptstyle\frown}{R}}_{12}}\left( :,1 \right)={\left( {{X}_{1}}{{X}_{2}}^{H}\left( 1,: \right) \right)}/{N}\;$     (21)

${{\overset{\scriptscriptstyle\frown}{R}}_{21}}\left( :,1 \right)={\left( {{X}_{2}}{{X}_{1}}^{H}\left( 1,: \right) \right)}/{N}\;$     (22)

Similarly, ${{X}_{2}}^{H}\left( 1,: \right)$ and ${{X}_{1}}^{H}\left( 1,: \right)$are the conjugate transposes of the first rows of ${{X}_{2}}$ and ${{X}_{1}}$, respectively. Then, $\overset{\scriptscriptstyle\frown}{R}\left( :,1 \right)$ can be written as:

$\overset{\scriptscriptstyle\frown}{R}\left( :,1 \right)\text{=}{{\overset{\scriptscriptstyle\frown}{R}}_{12}}\left( :,1 \right)+{{\overset{\scriptscriptstyle\frown}{R}}_{21}}\left( :,1 \right)$     (23)

where, $\overset{\scriptscriptstyle\frown}{R}\left( :,1 \right)$ is the first column of $R$. Thus, $\overset{\scriptscriptstyle\frown}{r}\left( :,1 \right)$ and ${{\overset{\scriptscriptstyle\frown}{r}}^{n}}\left( :,1 \right)$ can be respectively expressed as:

$\overset{\scriptscriptstyle\frown}{r}\left( :,1 \right)\text{=}\overset{\scriptscriptstyle\frown}{R}\left( :,1 \right)\text{=}{{\overset{\scriptscriptstyle\frown}{R}}_{12}}\left( :,1 \right)+{{\overset{\scriptscriptstyle\frown}{R}}_{21}}\left( :,1 \right)$     (24)

${{\overset{\scriptscriptstyle\frown}{r}}^{n}}\left( :,1 \right)=J{{\overset{\scriptscriptstyle\frown}{r}}^{*}}\left( :,1 \right)$     (25)

To sum up, the DPLA algorithm can be implemented in the following steps:

Step 1. Derive ${{\overset{\scriptscriptstyle\frown}{R}}_{12}}\left( :,1 \right)$ and ${{\overset{\scriptscriptstyle\frown}{R}}_{21}}\left( :,1 \right)$ from formulas (21) and (22).

Step 2. Obtain $\overset{\scriptscriptstyle\frown}{r}\left( :,1 \right)$ from formula (24).

Step 3. Obtain new vector from formula (25).

Step 4. Remove the same row in ${{\overset{\scriptscriptstyle\frown}{r}}^{n}}\left( :,1 \right)$ and $\overset{\scriptscriptstyle\frown}{r}\left( :,1 \right)$ to construct the column vector ${{r}_{\text{new}}}$.

Step 5. Obtain matrix $S$ by formula (14) and obtain matrix $\overline{S}$through Schmidt orthogonalization.

Step 6. Construct polynomial $f\left( z \right)$ by formula (16), and solve the polynomial to find the $K$ roots closest to the unit circle; calculate the DOA of the signals by formula (17).

3.2 Computing complexity analysis

The DPLA algorithm assumes that the signals are rotational invariant, due to the shift invariance of the DPLA, and derives the noise subspace from column vectors. Compared with traditional subspace algorithms, the DPLA algorithm eliminates the complex EVD or SVD, without sacrificing estimation accuracy. Hence, the algorithm could effectively simplify the computation, and shorten the estimation time.

The DPLA algorithm has 2M array elements, N snapshots, and K signal sources. The computing complexity of the algorithm mainly arises from  $\overset{\scriptscriptstyle\frown}{r}(:,1)$ acquisition, Schmidt orthogonalization, and polynomial rooting. The computing complexities of these three operations are $O\left( 4MN \right)$, $O\left( K\left( K+1 \right)\left( 2M-K \right)/2 \right)$, and $O\left( {{\left( 2M-K \right)}^{3}} \right)$, respectively.

For comparison, the computing compelxity of the MUSIC algorithm with ULA comes from EVD, the estaimation of covariance matrix, and the search for spectral peak. When the spectral peak is searched for at 500 equally spaced intervals as in our research, the computing complexities of the three operations in the MUSIC alogrithm are $O\left( {{\left( \text{2}M \right)}^{3}} \right)$, $O\left( \text{4}{{M}^{2}}N \right)$, and $O\left( {{\left( \text{2}M-K \right)}^{\text{1}5\text{0}0}} \right)$, respectively.

In actual applications, $N\gg M>K$. Therefore, the DPLA algorithm enjoys a much lower computing complexity than the MUSIC algorithm.

This paper puts forward a novel DPLA algorithm for DOA estimation, and verifies its superiority over traditional algorithms through MATLAB simulation. There are two major innovations in our algorithm: First, the rotational invariance of the two subarray signals was guaranteed by the shift invariance of the two subarrays, which eliminates the effect of array aperture and improves estimation accuracy. Secondly, our algorithm rearranges the first column elements of the JCCM into the noise subspace, and thus greatly reducing the computing load. The research results have great application potential in DOA estimation tasks.