<figure><img src="image_1.png" /></figure>

HAL
open science

# A new now complexity DOA estimation algorithm for massive MIMO systems

Yang Xiao, Liu Licheng, Yide Wang

► To cite this version:

Yang Xiao, Liu Licheng, Yide Wang. A new now complexity DOA estimation algorithm for massive MIMO systems. IEEE International Conference on Consumer Electronics China, Dec 2016, Guangzhou, China. 10.1109/ICCE-China.2016.7849735. [hal-01402284](https://hal.science/hal-01402284)

HAL Id: hal-01402284

https://hal.science/hal-01402284

Submitted on 11 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A New Low Complexity DOA Estimation Algorithm for Massive MIMO Systems

Xiao YANG

College of Information Engineering
Guangdong University of Technology
Guangzhou, China

Licheng LIU

College of Information Engineering
Guangdong University of Technology
Guangzhou, China

Yide WANG

IETR-UMR CNRS 6164
University of Nantes
Nantes, France

**Abstract**—Massive MIMO systems have attracted lots of attention, and its performance is closely related to the estimation of direction of arrivals (DOA) of incoming sources. The conventional DOA estimation algorithms like MUSIC and ESPRIT are not suitable for massive MIMO systems, due to the complexity brought by the large number of antennas. The DOA estimation algorithm named Propagator Method (PM) has been proved efficient in terms of computational complexity, because it does not require the eigenvalue decomposition (EVD) of the covariance matrix of the received signals. In this paper, a new version of PM is proposed. By partitioning the covariance matrix of the received signals, the data contaminated by the main contribution of noise can be avoided, and the reconstruction of the system model in not needed, so the computational complexity can be decreased significantly. Some numerical simulations are given to prove its efficiency.

**Keywords**—Massive MIMO; DOA Estimation; Propagator Method; Low Complexity

## I. INTRODUCTION

In the last years, significant attention has been attracted by massive multiple-input multiple-output (MIMO) systems. Compared with traditional MIMO systems, massive MIMO systems, equipping cellular base stations (BS) with a very large number of antennas, has a lot of advantages, such as the throughput, efficiency and reliability, and it has been considered as a promising key technology for the next generation mobile communication system. The performance of massive MIMO systems highly relies on the estimation of DOA of incoming signals, and DOA estimation is a crucial research domain of array signal processing. Conventional algorithms for DOA estimation like MUSIC [1] and ESPRIT [2] are not suitable for massive MIMO systems because they both need EVD of the covariance matrix of the received signals and/or (two-dimensional) 2-D peak searching, which have high computational complexity. The method for DOA finding named Propagator Method (PM) is proposed by Macros [3], [4], in which no EVD of the covariance matrix of the received signals is needed. Li extended PM into 2-D DOA estimation using a rectangular planar array [5], but the step of 2-D peak searching is still necessary. Wu proposed a new variant of PM for 2-D problem by using two parallel uniform linear arrays [6]. In this method, a set of linear operation is employed instead of the step of 2-D peak searching, so the computational complexity is greatly reduced. However, since

only the data from several rows of the propagator is used, the accuracy is influenced by the aperture loss. Li extended PM without peak searching into massive MIMO systems [7], in which all of the data from the propagator is employed and the aperture loss is avoided. In this paper, using the system model introduced in [7], a new version of PM for massive MIMO systems based on the partition of the covariance matrix of the received signals is proposed, aiming to further reduce the computational complexity. The main differences between the method in [7] and the proposed method are as follows: (1) by partitioning the covariance matrix of the received signals, the data contaminated by the main contribution of noise can be avoided; (2) the reconstruction of the system model is not needed. So the proposed method can reach a very low computational complexity and high responding speed, which are highly demanded in massive MIMO systems, with the least performance decrease.

The paper is organized as follows. The system model is introduced in Section 2. The proposed algorithm is presented in Section 3. In Section 4, some simulation results are shown and the conclusion is drawn in Section 5.

## II. SYSTEM MODEL

Consider a massive MIMO system with a rectangular planar array placed in X-Z plane. There are N antenna elements in X-axis and M antenna elements in Z-axis, with interelement spacing d and d=λ/2 (λ is the wavelength). The element located at the origin of the coordinate system, is referred to as the reference point, as shown in Fig.1. Suppose that there are K (K<<M, K<<N) uncorrelated, far-field and narrowband source signals, $s_i(t)$, $1 \le i \le K$, impinging on the array. The direction of arrival of the i-th signal is $(\theta_i, \varphi_i)$, which stands for the elevation angle and azimuth angle respectively.

The observed signals at the antenna array can be written as

$$X=AS+N \quad (1)$$

where S denotes the transmitted signal matrix with dimension of K×L, L is the number of snapshots and N is the Gaussian white noise matrix with dimension of MN×L.

A is the steering matrix of the antenna array, which can be defined as
<figure><img src="image_3.png" /><figcaption>Fig.1. System Model</figcaption></figure>

$$
\mathbf{A} = [\mathbf{a}(\theta_1, \varphi_1) \ \mathbf{a}(\theta_2, \varphi_2) \ \dots \ \mathbf{a}(\theta_K, \varphi_K)] \quad (2)
$$

where **a**(*θ*<sub>*i*</sub>, *φ*<sub>*i*</sub>) denotes the steering vector of the *i*<sup>th</sup> signal,
*i*=1,2,...,K.

Suppose that the system model is constructed following
the X-axis. **a**(*θ*<sub>*i*</sub>, *φ*<sub>*i*</sub>) can be defined as

$$
\begin{align*}
\mathbf{a}(\theta_i, \varphi_i) ={}& [1, e^{j\pi\nu_i}, \dots, e^{j\pi^{(N-1)}\nu_i}, e^{j\pi u_i}, \\
& e^{j\pi u_i + j\pi\nu_i}, \dots, e^{j\pi u_i + j\pi^{(N-1)}\nu_i}, \dots, e^{j\pi^{(M-1)}u_i}, \\
& e^{j(M-1)\pi u_i + j\pi\nu_i}, \dots, e^{j(M-1)\pi u_i + j\pi^{(N-1)}\pi\nu_i}]^T \tag{3}
\end{align*}
$$

with the following notations

$$
u_i = \cos \theta_i \tag{4}
$$

$$
v_i = \sin \theta_i \cos \varphi_i \quad (5)
$$

III. PROPOSED METHOD

With limited samples and noise, the covariance matrix of
the received signals can be estimated by

$$
\hat{\mathbf{R}} = \frac{1}{L} \mathbf{X} \mathbf{X}^{\mathrm{H}} \qquad (6)
$$

Introduce a partition of the covariance matrix $\hat{\mathbf{R}}$ like

$$
\hat{\mathbf{R}} = \begin{bmatrix} \mathbf{R}_{11} & \mathbf{R}_{12} \\ \mathbf{R}_{21} & \mathbf{R}_{22} \end{bmatrix} \quad (7)
$$

The dimensions of the four matrices **R**<sub>11</sub>, **R**<sub>12</sub>, **R**<sub>21</sub> and **R**<sub>22</sub> are (N×N), (N×(M-1)N), ((M-1)N×N) and ((M-1)N×(M-1)N) respectively. It is obvious that matrix **R**<sub>21</sub> avoids the diagonal elements of the covariance matrix, which contains the main contribution of noise that cannot be filtered by increasing the number of samples. So we focus on **R**<sub>21</sub> in the following.

Matrix $\mathbf{R}_{21}$ can be written as

$$
\mathbf{R}_{21} = \frac{1}{L} (\mathbf{A}_{\mathrm{P}} \mathbf{S}) (\mathbf{B}_{\mathrm{P}} \mathbf{S})^{\mathrm{H}} \quad (8)
$$

where $\mathbf{A}_P$ and $\mathbf{B}_P$ denotes the last $(M-1) \times N$ rows and the

first *N* rows of the steering matrix **A**.

The new steering matrix A<sub>P</sub> can be divided into two
sub-matrices A<sub>P1</sub> and A<sub>P2</sub> with dimension K×K and
[(M-1)N-K]×K respectively:

$$
\mathbf{A}_{\mathrm{P}} = \begin{bmatrix} \mathbf{A}_{\mathrm{P}1} \\ \mathbf{A}_{\mathrm{P}2} \end{bmatrix} \qquad (9)
$$

All the transmitted signals are supposed to be uncorrelated,
so the rank of A<sub>P</sub> is K and full, and the first K rows of A<sub>P</sub>
are linearly independent. Therefore there is a linear
transformation from A<sub>P1</sub> to A<sub>P2</sub> like

$$
\mathbf{A}_{\mathrm{P}2} = \mathbf{P}_{\mathrm{P}1}^{\mathrm{H}} \mathbf{A}_{\mathrm{P}1} \quad (10)
$$

Divide **B**<sub>P</sub> in the same way

$$
\mathbf{B}_{\mathrm{P}} = \begin{bmatrix} \mathbf{B}_{\mathrm{P}1} \\ \mathbf{B}_{\mathrm{P}2} \end{bmatrix} \qquad (11)
$$

with the dimensions of $\mathbf{B}_{\mathrm{P}1}$ and $\mathbf{B}_{\mathrm{P}2}$ are $K \times K$ and $(N-K) \times K$ respectively.

Similar, the relationship between them is

$$
\mathbf{B}_{\mathrm{P}2} = \mathbf{P}_{\mathrm{P}2}^{\mathrm{H}} \mathbf{B}_{\mathrm{P}1} \quad (12)
$$

where **P**<sub>P1</sub> and **P**<sub>P2</sub> are two propagator operators.

Divide **R**<sub>21</sub> in two ways

$$
\mathbf{R}_{21} = [\mathbf{R}_{21a} \quad \mathbf{R}_{21b}] \quad (13)
$$

and

$$
\mathbf{R}_{21} = \begin{bmatrix} \mathbf{R}_{21c} \\ \mathbf{R}_{21d} \end{bmatrix} \qquad (14)
$$

where

$$
\mathbf{R}_{21a} \in \mathbb{C}^{(M-1)N \times K}, \mathbf{R}_{21b} \in \mathbb{C}^{(M-1)N \times (N-K)}
$$

and

$$
\mathbf{R}_{21c} \in \mathbb{C}^{K \times N}, \quad \mathbf{R}_{21d} \in \mathbb{C}^{[(M-1)N-K] \times N}
$$

Then the two propagators **P**<sub>P1</sub> and **P**<sub>P2</sub> can be obtained
by

$$
\mathbf{P}_{\mathrm{P}1} = (\mathbf{R}_{21\mathrm{d}} \mathbf{R}_{21\mathrm{c}})^{\mathrm{H}} \quad (15)
$$

$$
\mathbf{P}_{\mathrm{P}2} = \mathbf{R}_{21a}^{+} \mathbf{R}_{21b} \quad (16)
$$

where $(\cdot)^{+}$ denotes the pseudoinverse, and then

$$
\mathbf{P}_{\mathrm{P}1} = (\mathbf{R}_{21\mathrm{d}}((\mathbf{R}_{21\mathrm{c}}^{\mathrm{H}}\mathbf{R}_{21\mathrm{c}})^{-1}\mathbf{R}_{21\mathrm{c}}))^{\mathrm{H}} \quad (17)
$$

$$
\mathbf{P}_{\mathrm{P}2} = (\mathbf{R}_{21a}^{\mathrm{H}} \mathbf{R}_{21a})^{-1} \mathbf{R}_{21a}^{\mathrm{H}} \mathbf{R}_{21b} \quad (18)
$$
Based on the two propagators, define two new matrices $P_{C1}$ and $P_{C2}$ with dimensions $(M-1)N \times K$ and $N \times K$ respectively, like

$$P_{PC1} = \begin{bmatrix} I_K \\ P_{P1}^H \end{bmatrix} \qquad (19)$$

$$P_{PC2} = \begin{bmatrix} I_K \\ P_{P2}^H \end{bmatrix} \qquad (20)$$

It is obvious that

$$P_{PC1}A_{P1} = \begin{bmatrix} I_K A_{P1} \\ P_{P1}^H A_{P1} \end{bmatrix} = \begin{bmatrix} A_{P1} \\ A_{P2} \end{bmatrix} = A_P \qquad (21)$$

$$P_{PC2}B_{P1} = \begin{bmatrix} I_K B_{P1} \\ P_{P2}^H B_{P1} \end{bmatrix} = \begin{bmatrix} B_{P1} \\ B_{P2} \end{bmatrix} = B_P \qquad (22)$$

Divide the matrix $P_{PC1}$ into two sub-matrices $P_{PC11}$ and $P_{PC12}$, composed by the first $N \times (M-2)$ rows and last $N \times (M-2)$ rows of $P_{PC1}$ respectively. Similarly, $A_{Pa}$ and $A_{Pb}$ are used to represent the $N \times (M-2)$ rows and last $N \times (M-2)$ rows of matrix $A_P$ respectively. By similar method, we use $P_{PC21}$ and $P_{PC22}$ to represent the first $N-1$ rows and last $N-1$ rows of $P_{PC1}$, and $B_{Pa}$ and $B_{Pb}$ to represent the first $N-1$ rows and last $N-1$ rows of $B_P$. It is easy to find that

$$\begin{bmatrix} P_{PC11} \\ P_{PC12} \end{bmatrix} A_{P1} = \begin{bmatrix} A_{Pa} \\ A_{Pb} \end{bmatrix} = \begin{bmatrix} A_{Pa} \\ A_{Pa} \emptyset_1 \end{bmatrix} \qquad (23)$$

$$\begin{bmatrix} P_{PC21} \\ P_{PC22} \end{bmatrix} B_{P1} = \begin{bmatrix} B_{Pa} \\ B_{Pb} \end{bmatrix} = \begin{bmatrix} B_{Pa} \\ B_{Pa} \emptyset_2 \end{bmatrix} \qquad (24)$$

where

$$\emptyset_1 = \text{diag}[e^{j\pi u_1}, e^{j\pi u_2}, \dots, e^{j\pi u_K}] \qquad (25)$$

$$\emptyset_2 = \text{diag}[e^{j\pi v_1}, e^{j\pi v_2}, \dots, e^{j\pi v_K}] \qquad (26)$$

and all the information of the direction of arrivals is contained in the above two diagonal matrices.

It can be deduced that

$$P_{PC11}^+ P_{PC12} = A_{P1} \emptyset_1 A_{P1}^{-1} \qquad (27)$$

$$P_{PC21}^+ P_{PC22} = B_{P1} \emptyset_2 B_{P1}^{-1} \qquad (28)$$

It is obvious that the diagonal elements of $\emptyset_1$ and $\emptyset_2$ can be obtained by performing the eigenvalue decomposition of the two matrices $P_{PC11}^+ P_{PC12}$ and $P_{PC21}^+ P_{PC22}$ respectively, and then we can obtain

$$u_i = \arg(\emptyset_1(i, i)) / \pi \qquad (29)$$

$$v_i = \arg(\emptyset_2(i, i)) / \pi \qquad (30)$$

The two sets of estimated results are obtained from two independent steps of eigenvalue decomposition, so a step of pair matching is necessary. Different methods of pair matching are proposed in [8] and [9]. After that, the DOA can be estimated by

$$\theta_i = \arccos(u_i) \qquad (31)$$

$$\varphi_i = \arccos\left(\frac{v_i}{\sqrt{1-(u_i)^2}}\right) \qquad (32)$$

## IV. SIMULATION AND ANALYSIS

The system model introduced in Section II is adopted in the simulation, with $M$ antenna elements in each column and $N$ antenna elements in each row of the antenna array. It is supposed that there are $K$ far-field uncorrelated signals impinging on the array and the number of snapshots is $L$. The major steps and complexity of the proposed method are shown as follows:

1) Construct the sub-covariance matrix of the received signals $R_{21}$ ...................................................$o(L(M-1)N^2)$;

2) Partition the matrix $R_{21}$ to estimate the first Propagator $P_{P1}$ and define the new matrix $P_{PC1}$ ...................................................$o(N^3+2KN^2+KN[(M-1)N-K])$;

3) Partition the matrix $R_{21}$ to estimate the other Propagator $P_{P2}$ and define the other new matrix $P_{PC2}$ ...................................................$o(K^3+2(M-1)NK^2+NK(M-1)(N-K))$;

4) Obtain the two diagonal matrices $\emptyset_1$ and $\emptyset_2$ via (27), (28) and get the estimation results...................................................$o(3K^2N(M-2)+3K^2(N-1)+4K^2)$;

It is can be seen that due to the partition of the covariance matrix, the reconstruction of the system model is not needed and the dimensions of the matrices used in computation are decreased, so the proposed method has a much lower computational complexity. Fig.2 shows the complexity comparison of the method in [7] and the proposed method with $K=3$.

<figure><img src="image_4.png" /><figcaption>Fig.2. Complexity Comparison with Different Number of Antenna Elements (L=200)</figcaption></figure>
For the performance measure, Root Mean Square Error (RMSE) is introduced, which is defined as

$$RMSE = \frac{1}{K} \sum_{k=1}^{K} \sqrt{\frac{1}{Q} \sum_{q=1}^{Q} (\hat{\beta}_{k,q} - \beta_k)^2} \quad (33)$$

where K is the number of signals and Q is the number of Monte Carlo trials. K=3 and Q=500 are used in the simulations with the number of antenna elements M=N=30. $\hat{\beta}_{k,q}$ is the estimation of elevation angle $\theta_k$ or azimuth angle $\varphi_k$ of the qth Monte Carlo trial.

The comparison of performance in different SNRs between the two algorithms is shown in Fig.3 and Fig.4. It can be seen that due to the partition of the covariance matrix of the received signals, some information needed in elevation angles estimation is lost, but the data contaminated by the main contribution of noise is avoided simultaneously, so the proposed method can reach a RMSE which is close to the other method's. However, more information needed in azimuth angles estimation is lost due to the partition of the covariance matrix, so the performance of the azimuth angles estimation will decrease.

<figure><img src="image_5.png" /><figcaption>Fig.3. Elevation Angles Estimation Performance</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig.4. Azimuth Angles Estimation Performance</figcaption></figure>

## V. CONCLUSION

In this paper, a new low complexity DOA estimation method for massive MIMO systems is proposed based on the Propagator Method. Taking the advantage of the partition of the covariance matrix of the received signals, the proposed approach does not require the reconstruction of the system model, and the data contaminated by the main contribution of noise is avoided to be used. The simulation results show that the computational complexity is significantly decreased, and a good elevation angles estimation can be reached, at the price of slight degradation in azimuth angles estimation. So this method can be widely used in the system which demands low computational complexity and fast responding, or is not sensitive to the estimation of azimuth angles.

## ACKNOWLEDGMENT

The authors are gratefully appreciative of the reviewers for their careful reading, helpful comments and valuable suggestions.

## REFERENCES

[1] Schmidt, R. O. (1986). Multiple emitter location and signal parameter estimation. IEEE Transactions on Antennas & Propagation, 34(3), 276-280.

[2] Roy, R., & Kailath, T. (1989). Esprit-estimation of signal parameters via rotational invariance techniques. IEEE Transactions on Acoustics Speech & Signal Processing, 37(7), 984-995.

[3] Marcos, S., Marsal, A., & Benidir, M. (1995). The propagator method for source bearing estimation. Signal Processing, 42(2), 121-138.

[4] Marcos, S., & Benidir, M. (1990). On a high resolution array processing method non-based on the eigenanalysis approach. International Conference on Acoustics, 5, 2955-2958.

[5] Li, P., Yu, B., & Sun, J. (1995). A new method for two-dimensional array signal processing in unknown noise environments. Signal Processing, 47(3), 319-327.

[6] Wu, Y., Liao, G., & So, H. C. (2003). A fast algorithm for 2-d direction-of-arrival estimation. Signal Processing, 83(8), 1827-1831.

[7] Li, X., Wu, W., & Zhang, X. (2015). A Low-Complexity DOA Estimation Algorithm for 2D Massive MIMO System. International Conference on Information Technology and Management Innovation.

[8] Dong, Y. Y., Dong, C. X., Xu, J., & Zhao, G. Q. (2016). Computationally efficient 2-d doa estimation for 1-shaped array with automatic pairing. IEEE Antennas & Wireless Propagation Letters, 15, 1669-1672.

[9] Kikuchi, S., Tsuji, H., & Sano, A. (2007). Pair-matching method for estimating 2-d angle of arrival with a cross-correlation matrix. IEEE Antennas & Wireless Propagation Letters, 5(1), 35-40.