\section{Introduction}

Massive MIMO (multiple input multiple output) is a key technology for current generation cellular networks \cite{Hoydis2013,Bjornson2018,Bjornson2019a}. Base stations (BSs) are deploying hundreds of antennas, enabling spatial multiplexing of several user equipments (UEs) per cell, both in uplink and downlink, with simple linear signal processing methods \cite{Rusek2013}. The increasing number of UEs, and as a consequence the sharing of the available pilots whenever the coherence blocks are limited in size, causes interference between UEs in neighboring cells. To make it possible for a BS to separate the UEs that it serves from the interference, second-order statistical information present in the channel covariance matrices has to be exploited during the channel estimation.

The channel covariance matrices are commonly assumed to be perfectly known in massive MIMO literature, which is a strong assumption since the matrix dimension grows with the number of antennas and furthermore the statistics are changing over time. Existing covariance matrix estimators either use additional pilot overhead, sacrifice samples for data transmission \cite{Bjornson2017a}, use strong separability conditions that might not always hold \cite{Haghighatshoar2017}  or require knowledge of the used pilots of all the UEs in the network \cite{Neumann2018, Xie2016}. Therefore, a new channel covariance matrix estimator for low-rank channel covariance matrices is presented in this paper, using a generalized eigenvalue decomposition (GEVD) of two covariance matrices that can be estimated from the available uplink data. The requirements for the system are minimal and, except for synchronization, there is no need for communication between the different cells and no prior knowledge on the background noise is required. Approximate MMSE estimators are also derived based on the newly proposed channel covariance matrix estimator. The effectiveness of the proposed methods is demonstrated in numerical simulations.


\section{Channel and Signal Model}
\label{cemMIMO:sec:channel_signal}

\subsection{Channel Model}
A cellular massive MIMO system with $L$ cells is considered where each cell contains a BS with $N$ antennas and $K$ single-antenna UEs. The channel from UE $k$ in cell $l$ to the BS in cell $j$ is denoted by $\mathbf{h}_{jlk} \in \mathbb{C}^N$  and is assumed to remain constant during a coherence block of $\tau_c$ samples. Superscript $.^t$ will be used to denote the quantities in coherence block $t$, e.g.  $\mathbf{h}^t_{jkl}$. A common channel model for $\mathbf{h}^t_{jkl}$ is that it is drawn from a correlated Rayleigh fading realization $\mathcal{NC}(\mathbf{0},\overline{\mathbf{R}}_{jlk})$, where $\overline{\mathbf{R}}_{jlk} \in \mathbb{C}^{N \times N}$ is the positive semi-definite channel covariance matrix describing the large-scale fading, including geometric pathloss, shadowing, antenna gains, and spatial channel correlation \cite{Bjornson2017}. The complex Gaussian distribution around $\overline{\mathbf{R}}_{jlk}$ models the small-scale fading. Due to the different positions of the UEs, it is assumed that the channels for different UEs are uncorrelated, i.e. $E\{\mathbf{h}^t_{jlk}  \mathbf{h}^{t,H}_{jmi}\} = \mathbf{0}$ if $\{l,k\} \ne \{m,i\}$, where $E\{.\} $ denotes the expected value operator, with respect to different channel realizations and $.^H$ is the Hermitian transpose operator. Measurements in \cite{Viering2002} show that the matrices $\{\overline{\mathbf{R}}_{jlk}\}$ remain constant over several thousands of coherence blocks, where this large number is denoted by $T$.


The scattering is mostly localized around the UE, since the BSs are elevated and hence have limited scattering in their near-field \cite{Xie2016}. The multipath components defining the total channel $\mathbf{h}_{jlk}$ are therefor arriving from a specific localized region in the network, so that the matrices $\{ \overline{\mathbf{R}}_{jlk} \}$ can be approximated by low-rank matrices, where the rank is denoted by $R$.  $\overline{\mathbf{R}}_{jlk}$ can thus be approximated as a sum of $R$ rank-1 terms:
\begin{equation}
\label{cemMIMO:eq:R_rank}
\overline{\mathbf{R}}_{jlk} = \sum_{r=1}^{R} \overline{\mathbf{q}}_{jlk,r} \overline{\mathbf{q}}_{jlk,r}^H
\end{equation}
where $\overline{\mathbf{q}}_{jlk,r}$ will be defined later. This implies that the channel $\mathbf{h}^t_{jlk}$ in coherence block $t$ can be modeled  as
\begin{equation}
\label{cemMIMO:eq:model_h}
\mathbf{h}^t_{jlk} = \sum_{r=1}^{R} \overline{\mathbf{q}}_{jlk,r} \overline{z}^t_{jlk,r} = \overline{\bm{\mathds{Q}}}_{jlk} \overline{\mathbf{z}}_{jlk}^t
\end{equation}
where the normal complex variables $\overline{\mathbf{z}}_{jlk}^t = [\overline{z}^t_{jlk,1} \ ... \ \overline{z}^t_{jlk,R}]^T$ with $.^T$ the transpose operator, are independently drawn from $\mathcal{NC}(0,1)$ and $\overline{\bm{\mathds{Q}}}_{jlk} = [\overline{\mathbf{q}}_{jlk,1} \ ... \ \overline{\mathbf{q}}_{jlk,R}]$.


\subsection{Signal Model}
The samples in one coherence block of $\tau_c$ samples are divided in $\tau_p$ samples for uplink channel estimation and $\tau_u$ samples for uplink data transmission with $\tau_c = \tau_p + \tau_u$. With  $\tau_p$ samples available, the network can predefine $\tau_p$ different orthogonal and unitary pilot sequences $\{\mathbf{s}_b\}$ where
\begin{equation}
\mathbf{s}_b = [s_b(1) \ ... \ s_b(\tau_p)]^T, \ b=1...\tau_p
\end{equation}
with $\mathbf{s}^H_b \mathbf{s}_b = \tau_p$ and $\mathbf{s}^H_b \mathbf{s}_c = 0$ if $b \ne c$. In the uplink channel estimation phase, each UE $k$ in each cell $l$ selects randomly one of the $\tau_p$ pilot sequences in each coherence block $t$ , the index of this pilot sequence is denoted by $b^t_{lk}$. An example is provided in Table \ref{cemMIMO:tab:training_sequences}. The UEs transmits the chosen pilot sequence with powers $\{p_{lk}\}$. The signal received at BS $j$ is then given as
\begin{equation}
\label{cemMIMO:eq:y_l-training}
\mathbf{y}^t_{j}(p) = \sum_{l=1}^L \sum_{k=1}^K \sqrt{p_{lk}} \mathbf{h}^t_{jlk} s_{b_{lk}^t}(p) +  \mathbf{n}^t_j(p), \ p=1...\tau_p
\end{equation}
where $p$ is the sample index and $\mathbf{n}^t_j(p)$ is the background noise at BS $j$. 

\begin{table}
    \hspace*{-0.16in}
    \resizebox{1.04\columnwidth}{!}{
    \begin{tabular}{l|c|c:c:c:c|c:c:c:c|c:c:c:c|cccc|c:c:c:c|}
        \multicolumn{1}{c}{} & \multicolumn{1}{c}{} & \multicolumn{20}{c}{{\large Coherence block}} \\ 
        \cline{3-22}
        \multicolumn{1}{c}{\begin{sideways}\end{sideways}} &  & \multicolumn{4}{c|}{} & \multicolumn{4}{c|}{} & \multicolumn{4}{c|}{} & \multicolumn{4}{c|}{} & \multicolumn{4}{c|}{} \\ [-1em]
        \multicolumn{1}{c}{\begin{sideways}\end{sideways}} &  & \multicolumn{4}{c|}{{\large 1}} & \multicolumn{4}{c|}{{\large 2}} & \multicolumn{4}{c|}{{\large 3}} & \multicolumn{4}{c|}{...} & \multicolumn{4}{c|}{{\large $T$}} \\ 
        \cline{3-22}
        \multicolumn{1}{c}{\begin{sideways}\end{sideways}} &  & $\mathbf{s}_1$ & $\mathbf{s}_2$ & $\mathbf{s}_3$ & $\mathbf{s}_4$ & $\mathbf{s}_1$ & $\mathbf{s}_2$ & $\mathbf{s}_3$ & $\mathbf{s}_4$ & $\mathbf{s}_1$ & $\mathbf{s}_2$ & $\mathbf{s}_3$ & $\mathbf{s}_4$ & \multicolumn{4}{c|}{...} & $\mathbf{s}_1$ & $\mathbf{s}_2$ & $\mathbf{s}_3$ & $\mathbf{s}_4$ \\ 
        \cline{2-22}
        \multirow{5}{*}{\rotcell{{\large UE}}} & 1 &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  &  & \multicolumn{2}{c}{...} & &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} \\ 
        \cline{2-22}
        & 2 &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  \multicolumn{4}{c|}{...}  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  \\ 
        \cline{2-22}
        & 3 & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  \multicolumn{4}{c|}{...}  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  \\ 
        \cline{2-22}
        & 4 &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  & & {\cellcolor[rgb]{0.753,0.753,0.753}}   &  &  &  &  \multicolumn{4}{c|}{...}  & {\cellcolor[rgb]{0.753,0.753,0.753}}  &  &  &  \\ 
        \cline{2-22}
        & ... & \multicolumn{1}{c}{} & \multicolumn{2}{c}{...} &  & \multicolumn{1}{c}{} & \multicolumn{2}{c}{...} &  & \multicolumn{1}{c}{} & \multicolumn{2}{c}{...} &  &  \multicolumn{4}{c|}{...}  & \multicolumn{1}{c}{} & \multicolumn{2}{c}{...} &  \\ 
        \cline{2-22}
        & $K$ &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} & {\cellcolor[rgb]{0.753,0.753,0.753}} &  &  &  &  \multicolumn{4}{c|}{...}  &  &  & {\cellcolor[rgb]{0.753,0.753,0.753}} &  \\
        \cline{2-22}
    \end{tabular}
    }
    \caption{Example of random allocation of 4 pilot sequences over $K$ UEs in a certain cell.}
    \label{cemMIMO:tab:training_sequences}
\end{table}

In the uplink data transmission phase, each UE $k$ in each cell $l$ transmits the unitary signal $s_{lk}(u)$ with zero mean and power $p_{lk}$. The received signal at BS $j$ is then given as
\begin{equation}
\label{cemMIMO:eq:y_l-data}
\mathbf{y}^t_{j}(u) = \sum_{l=1}^{L} \sum_{k=1}^K \sqrt{p_{lk}} \mathbf{h}^t_{jlk} s_{lk}(u) + \mathbf{n}^t_j(u), \ u=1...\tau_u.
\end{equation}
Note that the indices $p$ and $u$ are used to distinguish between the different phases. BS $j$ has to detect the uplink data from the UEs in its cell, i.e. $s_{jk}(u)$ for $k=1...K$. To perform receive combining strategies like maximum-ratio (MR), (regularized) zero-forcing (RZF) or minimum mean squared error (MMSE) combining \cite{Bjornson2017,Ozdogan2019}, good estimates of the channels of the $K$ UEs in cell $j$ are required. The background noise $\mathbf{n}_j^t$ is assumed to be zero-mean complex Gaussian, but with an unknown spatial covariance matrix $\overline{\mathbf{R}}_{\mathbf{n}_j\mathbf{n}_j} = E \{ \mathbf{n}_j^t \mathbf{n}_j^{t,H}\}$. It is often assumed to be white noise with a diagonal covariance matrix, but here a more general noise model is used, where the noise can be correlated between the antennas of a BS. This happens for example when there are located jamming signals nearby the BS.


\section{Optimal MMSE Channel Estimation}
\label{cemMIMO:sect:optimal_MMSE}

The signal obtained from \eqref{cemMIMO:eq:y_l-training} after despreading with pilot sequence $b^t_{jk}$ can be represented as
\begin{align}
\label{cemMIMO:eq:y}
\mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}] =& \sum_{u=1}^{\tau_p} \mathbf{y}^t_{j}(p) s^*_{b_{jk}^t}(p)  \\
= & \sum_{l=1}^L \sum_{i=1}^K \delta^{t}_{jk,li} \sqrt{p_{li}} \tau_p \mathbf{h}^t_{jli} + \underbrace{\sum_{u=1}^{\tau_p} \mathbf{n}^t_{j}(p) s^*_{b_{jk}^t}(p)}_{\mathbf{n}_{j}^{t,\text{pilot}}[b^t_{jk}]} \nonumber
\end{align}
where the random variable $\delta^{t}_{jk,li}$ is 1 with probability $1/\tau_p$ and 0 otherwise since all other UEs randomly pick 1 of the $\tau_p$ pilot sequences and  $\delta^{t}_{jk,jk}$ is deterministic and equal to 1.

Based on the despread signal in \eqref{cemMIMO:eq:y} a linear MMSE channel estimation criterion, similar to \cite[Sect. 4.2]{Demir;2021}, can be defined to obtain an estimate $\overline{\mathbf{h}}^t_{jjk} = \overline{\mathbf{W}}^H_{jk}\mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}]$ of $\mathbf{h}_{jjk}^t$  in the current coherence block $t$:
\begin{equation}
\label{cemMIMO:eq:W_k-optimal}
\overline{\mathbf{W}}_{jk} = \argmin_{\boldsymbol{W}_{jk}} E \{ ||\mathbf{h}^t_{jjk} - \boldsymbol{W}^H_{jk} \mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}]||^2 \}.
\end{equation}
The optimal solution of \eqref{cemMIMO:eq:W_k-optimal} is given as 
\begin{equation}
\label{cemMIMO:eq:W}
\overline{\mathbf{W}}_{jk} = \sqrt{p_{jk}} \left( \overline{\mathbf{R}}^{\text{pilot}}_{jjk} \right)^{-1} \overline{\mathbf{R}}_{jjk}
\end{equation}
where
\begin{align}
\label{cemMIMO:eq:R_k_dot}
\overline{\mathbf{R}}^{\text{pilot}}_{jjk} &= \frac{1}{\tau_p} E \{\mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}] \mathbf{y}_{j}^{t,\text{pilot},H}[b^t_{jk}] \}  \\
&= p_{jk} \tau_p \overline{\mathbf{R}}_{jjk} + \sum_{i \ne k} p_{ji} \overline{\mathbf{R}}_{jji} + \sum_{l\ne j} \sum_{i=1}^K p_{li} \overline{\mathbf{R}}_{jli} + \overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}} \nonumber
\end{align}
since the channels between different UEs are uncorrelated and are also uncorrelated with the background noise. The fact that $E \{ {\delta^{t^2}_{jk,li}}\} = \frac{1}{\tau_p}$ is also used.  For the background noise term, the following derivation is used:
\begin{align}
E\{&\mathbf{n}_{j}^{t,\text{pilot}}[b^t_{jk}] \mathbf{n}_{j}^{t,\text{pilot},H}[b^t_{jk}]\} \\
&= E\left\lbrace \left(\sum_{u=1}^{\tau_p} \mathbf{n}_{j}^t(p) s^*_{b_{jk}^t}(p)\right) \left(\sum_{u=1}^{\tau_p} \mathbf{n}_{j}^{t,H}(p) s_{b_{jk}^t}(p)\right)^H \right\rbrace \nonumber \\
&=  \sum_{u=1}^{\tau_p} E\{ |s_{b_{jk}^t}(p)|^2 \}  E\{\mathbf{n}_{j}^t \mathbf{n}_{j}^{t,H}\} = \tau_p \overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}} \nonumber
\end{align} 
where the cross terms in the second step are zero since the noise is uncorrelated between different samples.

An important difference with the MMSE channel estimation method proposed in \cite[Sect. 4.1,]{Demir;2021} is that now all UEs appear as interfering UEs in \eqref{cemMIMO:eq:R_k_dot}, while in \cite{Demir;2021} it is assumed that the pilot sequences for all the UEs are predetermined and fixed, so that only UEs sharing the same pilot sequence appear as interfering UEs. However, this assumption requires prior knowledge of the channel statistics of all the UEs and the background noise, while the proposed method allows for a more efficient estimation method which will be discussed next. 

\section{Low-rank Channel Estimation}
\subsection{Low-rank covariance matrix estimator}
\label{cemMIMO:subsect:covariance_estimation}
In practice the covariance matrices $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ and $\overline{\mathbf{R}}_{jjk}$ have to be estimated from the available data. $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ can be estimated using time averaging over the last $T$ coherence blocks, possibly using a diagonal regularizer as in \cite{Ledoit2004,Bjornson2017a}, since perfect synchronization and knowledge of $b_{jk}^t$ at the BS $j$ is assumed. The estimation will be denoted as $\mathbf{R}^{\text{pilot}}_{jjk}$ (without the overline). $\overline{\mathbf{R}}_{jjk}$ as in \eqref{cemMIMO:eq:R_rank} can not immediately be estimated, since the channel is unknown. However, the combined covariance matrix $\mathbf{R}^{\text{all}}_{j}$ can be estimated from all available antenna signals during both the channel estimation phase and uplink data transmission phase, which is an estimate of
\begin{align}
\label{cemMIMO:eq:R}
\overline{\mathbf{R}}^{\text{all}}_{j} &=  E \left\lbrace \frac{1}{\tau_p} \sum_{p=1}^{\tau_p} \mathbf{y}_{j}^{t}(p) \mathbf{y}_{j}^{t,H}(p) + \frac{1}{\tau_u} \sum_{u=1}^{\tau_u} \mathbf{y}_{j}^{t}(u) \mathbf{y}_{j}^{t,H}(u) \right\rbrace  \nonumber \\
&= \sum_{l=1}^L \sum_{i=1}^K p_{li} \overline{\mathbf{R}}_{jli} + \overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}} .
\end{align}
When comparing \eqref{cemMIMO:eq:R_k_dot} and \eqref{cemMIMO:eq:R}, it can be stated that $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ is the covariance matrix estimated after despreading and $\overline{\mathbf{R}}^{\text{all}}_{j}$ is the covariance matrix estimated before despreading. They differ only in a contribution depending on $\overline{\mathbf{R}}_{jjk}$. The following relation can thus be derived:
\begin{equation}
\label{cemMIMO:eq:R_estimate}
p_{jk} \overline{\mathbf{R}}_{jjk} = \frac{1}{\tau_p-1} \left( \overline{\mathbf{R}}^{\text{pilot}}_{jjk} - \overline{\mathbf{R}}^{\text{all}}_{j} \right).
\end{equation}
However, simply replacing $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ and $\overline{\mathbf{R}}^{\text{all}}_{j}$ with the estimated quantities $\mathbf{R}^{\text{pilot}}_{jjk}$ and $\mathbf{R}^{\text{all}}_{j}$ in \eqref{cemMIMO:eq:R_estimate} will not result in a good estimate $\mathbf{R}_{jjk}$ of $\overline{\mathbf{R}}_{jjk}$ since it will generally not satisfy the assumed rank assumption in \eqref{cemMIMO:eq:R_rank} and might even result in an indefinite estimate $\mathbf{R}_{jjk}$. One known remedy consists in using a rank-approximation based on the GEVD \cite{Serizel2014, Dendrinos1991} of the matrix pencil $\{{\mathbf{R}}^{\text{pilot}}_{jjk}, {\mathbf{R}}^{\text{all}}_{j} \}$ given as
\begin{align}
{\mathbf{R}}^{\text{pilot}}_{jjk} &= \mathbf{Q}_{jjk} \bm{\Sigma}_{jjk} \mathbf{Q}^H_{jjk} = \sum_{r=1}^{N} \sigma_{jjk,r} \mathbf{q}_{jjk,r} \mathbf{q}^H_{jjk,r}\\
{\mathbf{R}}^{\text{all}}_{j} &= \mathbf{Q}_{jjk} \mathbf{Q}^H_{jjk} = \sum_{r=1}^{N} \mathbf{q}_{jjk,r} \mathbf{q}^H_{jjk,r}
\end{align}
where $\bm{\Sigma}_{jjk} = \text{diag}\{\sigma_{jjk,1}, ..., \sigma_{jjk,N}\}$ is a diagonal matrix with the generalized eigenvalues sorted from large to small and the columns of $\mathbf{X}_{jjk} = \mathbf{Q}^{-H}_{jjk}$ contain the generalized eigenvectors. The eigenvectors are normalized such that $\mathbf{X}^H_{jjk} {\mathbf{R}}^{\text{all}}_{j} \mathbf{X}_{jjk} = \mathbf{I}$.
The optimal rank $R$ estimate $\mathbf{R}_{jjk}$ is given by keeping only the rank-1 terms belonging to the $R$ largest eigenvalues greater than one in \eqref{cemMIMO:eq:R_estimate}
\begin{equation}
\label{cemMIMO:eq:R_estimate_gevd}
p_{jk} \mathbf{R}_{jjk} = \sum_{r=1}^{R}  \frac{\sigma_{jjk,r}-1}{\tau_p-1} \mathbf{q}_{jjk,r} \mathbf{q}^H_{jjk,r} = \mathds{Q}_{jjk} \mathbb{\Lambda}_{jjk} \mathds{Q}^H_{jjk}
\end{equation}
where  $\mathbb{\Lambda}_{jjk} = \text{diag}\{\frac{\sigma_{jjk,1} -1}{\tau_p-1}, ..., $ $ \frac{\sigma_{jjk,R} -1}{\tau_p-1}\}$.

As stated in \cite{Serizel2014}, the GEVD-based approximation effectively selects the modes with the highest SINR, since approximation errors are weighted relative to the interference and noise aggregated in ${\mathbf{R}}^{\text{all}}_{j}$. Also note that unlike other rank-approximations of \eqref{cemMIMO:eq:R_estimate}, the GEVD-based approximation is immune to scaling and linear combining of the signals, i.e. the output signal and output SINR is independent of such scaling and combining, which is a desirable property.

The proposed covariance matrix estimator is data-driven and requires only synchronization and knowledge of the chosen pilot sequence $b^t_{jk}$ of the UEs in cell $j$ in each coherence block. There is no need for prior knowledge of $\overline{\mathbf{R}}_{jjk}$ (except for the rank $R$ of $\overline{\mathbf{R}}_{jjk}$), $\overline{\mathbf{R}}_{\mathbf{n}_j\mathbf{n}_j}$, or any prior knowledge on the statistics or pilot sequences of the UEs in other cells $l \ne j$.


\subsection{Low-rank approximate MMSE estimator}
\label{cemMIMO:subsect:channel_estimation}
When the estimated quantities are used to calculate the filter in \eqref{cemMIMO:eq:W_k-optimal}, the following filter is obtained:
\begin{align}
\mathbf{W}_{jk} &= \sqrt{p_{jk}} \left( \mathbf{Q}_{jjk} \bm{\Sigma}_{jjk} \mathbf{Q}^H_{jjk} \right)^{-1} \frac{1}{p_{jk}} \mathds{Q}_{jjk} \mathbb{\Lambda}_{jjk} \mathds{Q}^H_{jjk} \\
&= \frac{1}{\sqrt{p_{jk}}} \mathds{X}_{jjk} \mathbb{\Upsilon}_{jjk} \mathds{Q}^H_{jjk} \nonumber
\end{align}
with $\mathbb{\Upsilon}_{jjk} = \text{diag}\{\upsilon_{jjk,1}, ...,  \upsilon_{jjk,R}\} =\frac{1}{\tau_p-1} \text{diag}\{\frac{\sigma_{jjk,1} - 1}{\sigma_{jjk,1}},...,$ $   \frac{\sigma_{jjk,R} - 1}{\sigma_{jjk,R}}\}$ and $\mathds{X}_{jjk} = [\mathbf{x}_{jjk,1} \ ... \ \mathbf{x}_{jjk,R}]$ contains the first $R$ columns of $\mathbf{X}_{jjk}$. The final low-rank approximate MMSE estimate for the channel of UE $k$ in coherence block $t$ is then given as
\begin{align}
\label{cemMIMO:eq:h_estimate}
\hat{\mathbf{h}}^t_{jjk} &= \mathbf{W}^H_{jk} \mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}]
= \frac{1}{\sqrt{p_{jk}}} \mathds{Q}_{jjk} \mathbb{\Upsilon}_{jjk} \mathds{X}^H_{jjk} \mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}] \nonumber \\
&= \frac{1}{\sqrt{p_{jk}}} \sum_{r=1}^{R} \mathbf{q}_{jjk,r} \underbrace{\upsilon_{jjk,r} \left( \mathbf{x}^H_{jjk,r} \mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}] \right)}_{z^t_{jjk,r}}. 
\end{align}
Comparing this expression for $\hat{\mathbf{h}}^t_{jjk}$ with the model used for ${\mathbf{h}}^t_{jjk}$ in \eqref{cemMIMO:eq:model_h}, it can be observed that the estimation method produces an estimate that exactly fits in this model. The random variable $\bar{\mathbf{z}}^t_{jjk}$ of \eqref{cemMIMO:eq:model_h} is estimated from the despread signal $\mathbf{y}_{j}^{t,\text{pilot}}[b^t_{jk}]$ using a simple inner-product with the corresponding eigenvector $\mathbf{x}_{jjk,r}$ and some scalars depending on the generalized eigenvalues.

\subsection{Improved low-rank approximate MMSE estimator}
An improved low-rank approximate MMSE estimator can be derived based on \eqref{cemMIMO:eq:y}. When computing the expectation in \eqref{cemMIMO:eq:W_k-optimal}, the activity of all the UEs is represented by the random variable $\delta^{t}_{jk,li}$. Since BS $j$ is aware of the used pilot sequences of the UEs in its cell, this information can be used to improve the MMSE estimator:
\begin{align}
\label{cemMIMO:eq:y_2}
\mathbf{y}_{j}^{t,\text{pilot,impr}}[b^t_{jk}] =& \sum_{i=1}^K \bar{\delta}^{t}_{jk,ji} \sqrt{p_{ji}} \tau_p \mathbf{h}^t_{jji} \\
&+ \sum_{l \ne j} \sum_{i=1}^K \delta^{t}_{jk,li} \sqrt{p_{li}} \tau_p \mathbf{h}^t_{jli} + \mathbf{n}_{j}^{t,\text{pilot}}[b^t_{jk}] \nonumber
\end{align}
where $\bar{\delta}^{t}_{jk,ji}$ is now deterministic and either equal to 0 or 1 and $\bar{\delta}^{t}_{jk,jk} = 1$. The corresponding cell-specific covariance matrix that now depends on the coherence block index $t$ is given as (alternative to \eqref{cemMIMO:eq:R_k_dot})
\begin{align}
\label{cemMIMO:eq:R_k_dot_2}
\overline{\mathbf{R}}^{\text{pilot,impr,t}}_{jjk} &= \frac{1}{\tau_p} E \{\mathbf{y}_{j}^{t,\text{pilot,impr}}[b^t_{jk}] \mathbf{y}_{j}^{t,\text{pilot,impr},H}[b^t_{jk}] \}  \\
&= \sum_{i=1}^K \bar{\delta}^{t^2}_{jk,ji} p_{ji} \tau_p \overline{\mathbf{R}}_{jji} + \sum_{l\ne j} \sum_{i=1}^K p_{li} \overline{\mathbf{R}}_{jli} + \overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}} \nonumber \\
&= \overline{\mathbf{R}}^{\text{pilot}}_{jjk} + \sum_{i \ne k} (\bar{\delta}^{t^2}_{jk,ji} \tau_p - 1) p_{ji} \overline{\mathbf{R}}_{jji}. \nonumber
\end{align}
where $ \overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ is defined as in \eqref{cemMIMO:eq:R_k_dot}. By replacing the covariance matrices $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ and $\overline{\mathbf{R}}_{jji}$ with their estimated quantities as explained in Subsection~\ref{cemMIMO:subsect:covariance_estimation}, an improved low-rank approximate MMSE estimate 
$\hat{\mathbf{h}}^{t,\text{impr}}_{jjk} = \mathbf{W}^{\text{impr},t,H}_{jk} \mathbf{y}_{j}^{t,\text{pilot,impr}}[b^t_{jk}] $ can be obtained where
\begin{equation}
    \label{cemMIMO:eq:W_2}
    \mathbf{W}^{\text{impr},t}_{jk} = \Big( \mathbf{R}^{\text{pilot}}_{jjk} + \sum_{i \ne k} (\bar{\delta}^{t^2}_{jk,ji} \tau_p - 1) p_{ji} \mathbf{R}_{jji} \Big)^{-1} \mathbf{R}_{jjk}.
\end{equation}

Although this estimator uses more of the information available at BS $j$, it is important to note that  the provided improved low-rank approximate MMSE estimator requires an extra matrix inversion for each UE in its cell, while this is not required for the low-rank approximate MMSE estimator of Subsection~\ref{cemMIMO:subsect:channel_estimation}. 


\section{Numerical Simulations}
\label{cemMIMO:sect:sim}
The same multicell setup with $L = 7$ hexagonal cells as in \cite{Bjornson2017a, Neumann2018} is considered. There are $K = 10$ UEs per cell in a ring around the BS (there are 70 UEs in total) and each BS has a uniform linear array with $N = 100$ antennas. The multipath components from a UE arrive uniformly distributed in an angular interval ($10^\circ$) centered around the geographical angle to the UE \cite{Adhikary2013}. This results in channel covariance matrices with approximately 25 dominant eigenvalues.

The normalized MSE (NMSE) (averaged over 20 Monte-Carlo runs), $E\{ ||\mathbf{h}_{jjk} - \hat{\mathbf{h}}_{jjk}||^2\}/\text{tr}\{\overline{\mathbf{R}}_{jjk}\}$, for an average UE in the center cell is used as a performance metric. The proposed low-rank approximate MMSE estimator is calculated for different ranks: $R=30$ and $R=60$, denoted with \textit{GEVD 30} and \textit{GEVD 60} respectively, and also when the subtraction in \eqref{cemMIMO:eq:R_estimate} is used instead of the GEVD to estimate the channel covariance matrices, denoted with \textit{SUBT}.  The performance of the MMSE estimator with perfect channel covariance matrix knowledge (\textit{MMSE random}) in \eqref{cemMIMO:eq:W} is also shown. Similar results are shown for the proposed improved approximate MMSE estimator, denoted with the word \textit{impr}. As baseline, the simple least-square estimation (\textit{LS fixed}) that does not require covariance matrix information is shown in a scenario where each UE is allocated with a fixed pilot sequence for all coherence blocks \cite{Bjornson2017a}. Also the MMSE estimate with perfect channel covariance matrix knowledge (\textit{MMSE fixed}) for the fixed pilot allocation regime is provided as absolute lower bound.

Figure~\ref{cemMIMO:fig:T_run} shows that a relative small number of coherence blocks (300) suffices to obtain an accurate channel estimate for the proposed methods. The estimators with a rank constraint on the channel covariance matrix outperform the estimators without the rank constraint (SUBT) and $R=30$ works better for the current setup than $R=60$, since it can already capture the dominant eigenspace. The proposed GEVD-based procedure focuses on the modes with the highest power in the estimated covariance matrices and neglect the other modes, that typically contain more estimation error noise than useful signal power. The improved approximate MMSE estimator is also outperforming the approximate MMSE estimator when there are enough coherence blocks available for the estimation procedure. Since all the estimates of the channel covariance matrices are used to calculate the filter in \eqref{cemMIMO:eq:W_2}, this will only perform well when they are estimated using enough samples. Although not shown in Figure~\ref{cemMIMO:fig:T_run}, it should be noted that the \textit{SUBT}-curve will converge to the \textit{MMSE fixed}-curve when $T \rightarrow \infty$ while the \textit{GEVD}-curves perform worse, since the rank assumption on the channel covariance matrices is only approximately valid for this setup.

In Figure~\ref{cemMIMO:fig:tau_run}, the performance with respect to the number of available pilot sequences $\tau_p$ is shown for $T=1500$ coherence blocks. The seemingly arbitrary shapes of the curves with a fixed pilot allocation emerge from the fact that the pilots are allocated in a cyclic way, starting from the leftmost UE in a cell and the cycle is not restarted in the neighboring cell. It is again observed that the proposed low-rank estimators outperform the estimates obtained with subtraction of the covariance matrices and the baseline LS-method. \textit{GEVD 30} provides the best performance, independently of the number of available pilot sequences. Increasing $\tau_p$, increases the dominance of $\overline{\mathbf{R}}_{jjk}$ in \eqref{cemMIMO:eq:R_k_dot}, and has thus a positive impact on the channel estimation performance.


\begin{figure}[]
    \centering
    \inputTikZ{0.46}{NLMS_T_run.tikz}  
    \caption{NMSE for an average intra-cell UE for different estimators with $\tau_p=10$ when the number of coherence block is varied.}
    \label{cemMIMO:fig:T_run}
\end{figure}

\begin{figure}[]
    \centering
    \inputTikZ{0.46}{NLMS_tau_p_run.tikz}  
    \caption{NMSE for an average intra-cell UE for different estimators with $T=1500$ when the number of pilot sequences is varied.}
    \label{cemMIMO:fig:tau_run}
\end{figure}

\section{Conclusion}
\label{cemMIMO:sect:con}
In this paper, a new uplink channel covariance matrix estimator has been presented for low-rank channel covariance matrices, using a GEVD of two covariance matrices that are estimated from the available uplink data. The requirements for the systems are minimal, i.e. UEs in a cell can randomly choose a pilot sequence in each coherence block and only the BS in its cell should be aware of its choice. Except for synchronization, there is no need for communication between the different BSs and no prior knowledge on the background noise is required. The derived approximate MMSE and improved approximate MMSE estimators have been shown to provide a good estimate of the true channel, even when only a small number of coherence blocks is available, while a simple subtraction of the two covariance matrices provides inferior performance due to estimation errors.


\bibliographystyle{IEEEbib}
