Massive MIMO Extensions to the COST 2100 Channel Model: Modeling and Validation

TL;DR

This paper extends the COST 2100 channel model for massive MIMO systems by introducing visibility regions at both the base station and mobile station sides, supported by measurement-based statistical analysis and validation.

Contribution

It proposes new models for cluster and MPC visibility regions, including their statistical distributions, and validates these extensions through analysis and simulations.

Findings

BS-VR lifetimes are exponentially distributed.

Number of BS-VRs follows a Poisson distribution.

MPC lifetimes and radii are lognormally distributed.

Abstract

To enable realistic studies of massive multiple-input multiple-output systems, the COST 2100 channel model is extended based on measurements. First, the concept of a base station-side visibility region (BS-VR) is proposed to model the appearance and disappearance of clusters when using a physically-large array. We find that BS-VR lifetimes are exponentially distributed, and that the number of BS-VRs is Poisson distributed with intensity proportional to the sum of the array length and the mean lifetime. Simulations suggest that under certain conditions longer lifetimes can help decorrelating closely-located users. Second, the concept of a multipath component visibility region (MPC-VR) is proposed to model birth-death processes of individual MPCs at the mobile station side. We find that both MPC lifetimes and MPC-VR radii are lognormally distributed. Simulations suggest that unless MPC-VRs are applied the channel condition number is overestimated. Key statistical properties of the proposed extensions, e.g., autocorrelation functions, maximum likelihood estimators, and Cramer-Rao bounds, are derived and analyzed.

Tables

Table 1. TABLE I: Estimated parameters and means of cases 1, 2, and 3.

Case	$\hat{\lambda }$	${\widehat{𝜽}}_Y$	$E(Y;{\widehat{𝜽}}_Y)$
1	52.30	${\hat{L}}_{\mathrm{MPC}}=0.81$	0.81
2	171.60	${\hat{\mu }}_{\mathrm{MPC}}=-16.92,$ ${\hat{\sigma }}_{\mathrm{MPC}}^2=94.60$	0.25
3	52.13	${\hat{\mu }}_{\mathrm{MPC}}=-11.09,$ ${\hat{\sigma }}_{\mathrm{MPC}}^2=89.91$	0.84

Table 2. TABLE II: Parametrization of the COST 2100 model extension for closely-located users with PLA and CLA at 2.6 Ghz.

Parameter	Outdoor VLA (NLOS)1	Indoor CLA (LOS)
Length of BS-VRs
$L_{\mathrm{BS}}$ [m]	$3.2$	-
Slope of BS-VR gain
${\mu }_{\mathrm{BS}}$ [dB/m]	0	-
${\sigma }_{\mathrm{BS}}$ [dB/m]	0.9	-
MPC gain function
${\mu }_{R_{\mathrm{MPC}}}$ [dB]	-	-19.8
${\sigma }_{R_{\mathrm{MPC}}}$ [dB]	-	10.1
Average number of visible far clusters
$N_C$	2.9$\times \left(L_{\mathrm{BS}}+L\right)$	15
Radius of the cluster visibility region
$R_C$ [m]	10	5
Radius of cluster transition region
$T_C$ [m]	2	0.5
Number of MPCs per cluster
$N_{\mathrm{MPC}}$	31	10002
Cluster power decay factor
$k_{\tau }$ [dB/$\mu$s]	43	31
Cluster cut-off delay
${\tau }_B$ [$\mu$s]	0.91	0.25
Cluster shadowing
${\sigma }_S$ [dB]	7.6	2.7
Cluster delay spread
$m_{\tau }$ [$\mu$s]	0.14	0.005
$S_{\tau }$ [dB]	2.85	1.5
Cluster angular spread in azimuth (at BS)
$m_{{\psi }_{\mathrm{BS}}}$ [deg]	7.0	4.6
$S_{{\psi }_{\mathrm{BS}}}$ [dB]	2.4	2.1
Cluster angular spread in elevation (at BS)
$m_{{\theta }_{\mathrm{BS}}}$ [deg]	0	3.7
$S_{{\theta }_{\mathrm{BS}}}$ [dB]	0	2.6
Cluster angular spread in azimuth (at MS)
$m_{{\psi }_{\mathrm{MS}}}$ [deg]	19	3.63
$S_{{\psi }_{\mathrm{MS}}}$ [dB]	2.0	2.13
Cluster angular spread in elevation (at MS)
$m_{{\theta }_{\mathrm{MS}}}$ [deg]	0	0.73
$S_{{\theta }_{\mathrm{MS}}}$ [dB]	0	3.63
Cluster spread cross-correlation
${\rho }_{{\sigma }_S\tau }$	-0.09	-0.45
${\rho }_{{\sigma }_S{\psi }_{\mathrm{BS}}}$	0.04	-0.56
${\rho }_{{\sigma }_S{\theta }_{\mathrm{BS}}}$	0	-0.50
${\rho }_{\tau {\psi }_{\mathrm{BS}}}$	0.42	0.70
${\rho }_{\tau {\theta }_{\mathrm{BS}}}$	0	0.34
${\rho }_{{\psi }_{\mathrm{BS}}{\theta }_{\mathrm{BS}}}$	0	0.50
Radius of LOS visibility region
$R_L$ [m]	-	303
Radius of LOS transition region
$T_L$ [m]	-	03
LOS power factor
${\mu }_{K_{\mathrm{LOS}}}$ [dB]	-	-5.2
${\sigma }_{K_{\mathrm{LOS}}}$ [dB]	-	2.9
XPR
${\mu }_{\mathrm{XPR}}$ [dB]	0	94
${\sigma }_{\mathrm{XPR}}$ [dB]	0	34

Equations

$$N_{\mathrm{birth}}(x)=\text{ \# births events in }(0,x],$$

$$N_{\mathrm{new}}(x_1,x_2)=N_{\mathrm{birth}}(x_2)-N_{\mathrm{birth}}(x_1).$$

$$N_{\mathrm{new}}(x_1,x_2)\in \mathrm{Po}(\lambda \cdot (x_2-x_1)),$$

$$p(n)=e^{-m}\frac{m^n}{n!},n=0,1,2,\dots ,$$

$$N(x_1,x_2)\in \mathrm{Po}(\lambda \cdot (x_2-x_1)+\lambda \cdot E(Y)).$$

$$P(Y\ge x_1-u)=1-F_Y(x_1-u).$$

$$\lambda (x)=\lambda \cdot P(Y\ge x_1-x).$$

$$N_{\mathrm{alive}}(x_1)=N_{\mathrm{birth}}^{\prime }(x_1)-N_{\mathrm{birth}}^{\prime }(0)=\mathrm{Po}(\lambda \cdot {\int }_0^{x_1}P(Y\ge x_1-u)du).$$

$$N_{\mathrm{alive}}(x_1)=\underset{x_1\to \infty }{\lim }\mathrm{Po}(\lambda \cdot {\int }_0^{x_1}P(Y\ge u)du)=\mathrm{Po}(\lambda \cdot E(Y)),$$

$$N(x_1,x_2)\to N_{\mathrm{alive}}(x)\text{ as }{x}_1,x_2\to x,$$

$$N(x)\in \mathrm{Po}(\lambda \cdot E(Y)).$$

$${\mathcal{N}}_{00}=$$

$${\mathcal{N}}_{10}=$$

$$n=n_{00}+n_{01}+n_{10}+n_{11}$$

$$\mathcal{L}(\mathbf{\theta };\mathbf{x})$$

$$\times \prod _{i\in {\mathcal{N}}_{00}}{f}_Y({\upsilon }_i;{\mathbf{\theta }}_Y)\times \prod _{j\in {\mathcal{N}}_{01}\cup {\mathcal{N}}_{10}}{\int }_{{\upsilon }_j}^{\infty }{f}_Y(t;{\mathbf{\theta }}_Y)dt$$

$$\times {\left({\int }_L^{\infty }(t-L)f_Y(t;{\mathbf{\theta }}_Y)dt\right)}^{n_{11}},$$

$$(\hat{\lambda },{\hat{\mathbf{\theta }}}_Y)$$

$$f_Y(y)=\{\begin{array}{ll}\frac{1}{L_{\mathrm{BS}}}{e}^{-\frac{y}{L_{\mathrm{BS}}}},& \text{for }y\ge 0;\\ 0,& \text{otherwise,}\end{array}$$

$$\hat{\lambda }$$

$${\hat{L}}_{\mathrm{BS}}$$

$$\mathcal{L}(\mathbf{\theta };\mathbf{x})\propto {\lambda }^n{e}^{-\lambda (L-{\Delta }_0+L_{\mathrm{BS}})e^{-\frac{{\Delta }_0}{L_{\mathrm{BS}}}}}{L}_{\mathrm{BS}}^{\nu }{e}^{-{\sum }_{i=1}^n\frac{{\upsilon }_i}{L_{\mathrm{BS}}}},$$

$$\frac{\partial \log \mathcal{L}(\mathbf{\theta };\mathbf{x})}{\partial \lambda }$$

$$\frac{\partial \log \mathcal{L}(\mathbf{\theta };\mathbf{x})}{\partial L_{\mathrm{BS}}}$$

$$0$$

$$-\frac{n}{L+L_{\mathrm{BS}}}+\frac{\nu }{L_{\mathrm{BS}}}+\frac{{\sum }_{i=1}^n{\upsilon }_i}{L_{\mathrm{BS}}^2}=0,$$

$$a{L}_{\mathrm{BS}}^2+b{L}_{\mathrm{BS}}+c=0$$

$${\hat{L}}_{\mathrm{BS}}=\frac{-b\pm \sqrt{b^2-4ac}}{2a}.$$

$${\hat{L}}_{\mathrm{BS}}=\frac{-b+\sqrt{b^2-4ac}}{2a}$$

$$\mathbf{I}(\mathbf{\theta })={\lambda }_0\left[\begin{array}{cc}(L_{\mathrm{BS}}+L_0)/{\lambda }^2& 1\\ 1& (L_{\mathrm{BS}}+L_0)/L_{\mathrm{BS}}^2\end{array}\right].$$

$${\mathrm{CRLB}}_{\lambda }=[\mathbf{I}(\mathbf{\theta }{)}^{-1}{]}_{11},{\mathrm{CRLB}}_{L_{\mathrm{BS}}}=[\mathbf{I}(\mathbf{\theta }{)}^{-1}{]}_{22}.$$

$$\mathrm{Var}\left(\frac{\hat{\lambda }}{\lambda }\right)\ge \frac{{\mathrm{CRLB}}_{\lambda }}{{\lambda }^2}=\left(\frac{1}{{\lambda }_0{L}_0}\right)\frac{1+L_{\mathrm{BS}}/L_0}{1+2L_{\mathrm{BS}}/L_0},$$

Full text

I Introduction

Massive multiple-input multiple-output (MIMO), MaMi for short, has since its inception [3] attracted considerable attention in the wireless communication community [4, 5, 6]. Within the last few years, an abundant body of theoretical [5, 6] and experimental [7, 8] research has shown that MaMi systems can improve the energy and spectral efficiencies of today’s wireless communication systems by one to two orders of magnitude. For this reason, MaMi is considered a crucial component of the new radio (NR) air interface of the fifth generation (5G) wireless communication standard [9].

As is well known, the radio propagation channel ultimately dictates the performance of any wireless communication system. The availability of sufficiently accurate propagation channel models is therefore of critical importance to the design and evaluation of new wireless systems. The 3GPP SCM-3D [10], the WINNER II/WINNER+ [11, 12], and the COST 2100 [13, 14] are examples of channel models customarily used in the development and validation of 5G networks. Among them, the Quasi Deterministic Radio Channel Generator (QuaDRiGa) [15, 16], an extension of the 3GPP-3D and WINNER II/WINNER+ channel models, is especially popular because of its enhanced spatiotemporal consistency. Moreover, it possesses certain features that enable MaMi simulations, namely support for spherical wavefronts, dispersion of clusters in elevation, and independence of arrival/departure angles [16].

Experiment shows [17, 18, 19] that at least two additional aspects of the MaMi propagation channel need to be considered. The first of them [17, 18] is the presence at the base station (BS) of non-stationarities caused by the appearance and disappearance of clusters along physically-large arrays (PLAs). Such BS-side non-stationarities have been addressed in theoretical MaMi channel models [20, 21] in which the appearances and disappearances of clusters are typically modeled as Markov birth-death processes obeying certain survival probabilities. However, multiuser consistency becomes problematic in such models, because the Markov birth-death processes are not linked to the geometry of the environment. Whereas the original QuaDRiGa model neither considers BS-side non-stationarities, nor ensures multiuser consistency, some progress has been made by Oliveras et al. [22] to remedy these difficulties. However, the novel ideas presented in [22] are yet to be validated by experiments, and integrated in the mainstream QuaDRiGa model.

Secondly, it appears [23, 24, 25, 26, 27, 28, 29, 30, 31] that, just as clusters, individual multipath components (MPCs) undergo birth-death processes. The MPC lifetimes when tracked along a straight-line segment at the mobile station (MS) side have been found to obey certain probability distributions, e.g., the exponential distribution [23, 28, 31], the lognormal distribution [26, 29, 30, 31], and the Birnbaum-Sanders distribution [30, 31]. More important, closely-located users only a few wavelengths apart can encounter different MPCs [19]. Channel models that incorporate birth-death processes at the level of individual MPCs are of interest to emerging applications such as MaMi communications [32, 33], vehicle-to-vehicle communications [28, 30], or multipath-assisted localization [26, 27, 31, 34], which can use this information to improve their performance. To the best of our knowledge, none of the earlier referred 5G channel models currently considers birth-death processes for individual MPCs.

In this paper, MaMi extensions of the COST 2100 channel model are presented that describe BS-side non-stationarities of the clusters and MS-side birth-death processes of individual MPCs, introduced above. The main contributions of the paper are the following:

• •

The notion of BS-side visibility region (BS-VR) is proposed to model BS-side non-stationarities. Measured outdoor MaMi channels show that the number of BS-VRs is Poisson distributed and their lifetimes obey an exponential law.

• •

The maximum likelihood estimator and the Cramer-Rao bound of the parameters of the BS-VRs are rigorously derived. A closed-form expression of the autocorrelation function (ACF) of the number of the number of observed BS-VRs is provided.

• •

The influence of the choice of BS-VR parameters on the MaMi propagation channel is investigated through simulations. We find that under certain conditions longer BS-VRs can help decorrelating closely-located users.

• •

The notions of MPC-VR and MPC gain function are proposed to model birth-death processes of individual MPCs at the MS side. Measured indoor MaMi channels reveal that MPC lifetimes and MPC-VR radii are lognormally distributed.

• •

The small-scale fading ACF with the MPC gain function is always smaller than without it. When the gain function is applied to a group of closely-located users their condition number decreases by a few dB and comes to agree with measurements.

The present work extends our preliminary studies in [1] with further analyses, revised insights, and new measurements to characterize the MPC gain function. A MATLAB implementation of the COST 2100 model with the proposed extensions is freely available at [35].

II Notation

Throughout the paper, boldface lowercase letters represent column vectors, boldface uppercase letters matrices, and calligraphic letters sets. Thus, $\boldsymbol{I}$ denotes the identity matrix, $\boldsymbol{A}^{\operatorname{T}}$ the transpose, $\boldsymbol{A}^{{\dagger}}$ the Hermitian transpose, $i$ the imaginary unit, $\operatorname{\mathbb{R}}^{+}$ the set $[0,\infty),$ $\operatorname{Pr}(A)$ the probability of event $A,$ $1_{A}(\cdot)$ the indicator function of $A$, $\operatorname{E}(\cdot)$ the expectation operator, and $\Phi(x)$ the CDF of $\mathcal{N}(0,1).$

III The COST 2100 MIMO Channel Model

We briefly review some concepts and terminology of the COST 2100 model. For further details, the reader is referred to [13, 14]. The COST 2100 model falls into the category of geometry-based stochastic channel models (GSCM) [36], which are built around the notion of a cluster, i.e., a group of MPCs with similar delay and angular parameters [37]. Clusters model interactions of the transmit signal with scattering objects in the environment such as building facades, trees and street furniture, in outdoor environments, or inner walls, pillars and office equipment, in indoor settings. Interactions happen at so-called scattering points. In GSCMs, the locations of clusters and scattering points are drawn randomly from prescribed multivariate distributions depending on the simulated environment. MPCs are then obtained by mere geometric ray tracing of the transmit signal through the scattering points [13].

The distinguishing feature of the COST 2100 model is perhaps the use of so-called visibility regions (VRs) to model non-stationarities in a spatiotemporally consistent way. When a MS moves inside a VR, the associated cluster is active and its MPCs contribute to the MS-BS double-directional impulse response. Similarly, when the MS moves out of a VR, the associated cluster is no longer visible and it does not contribute. The situation is depicted in Fig. 1a. The two MaMi channel aspects alluded to in the introduction are also shown. The first aspect, illustrated in Fig. 1b, has to do with the presence of non-stationarities at the BS side in MaMi settings with PLAs [17, 18]; the second one, illustrated in Fig. 1c, concerns the realization that individual MPCs undergo birth-death processes at the MS side [23, 24, 25, 26, 27, 28, 29, 30, 31]. To model these phenomena, the COST 2100 channel model is extended with the concepts of BS-VRs, in Sec. IV, and MPC-VRs and MPC gain functions, in Sec. V.

IV Extension for Physically-Large Arrays

IV-A Measurements and Data Processing

Channels used in this section were measured at the parking lot outside the E-building of the Faculty of Engineering of Lund University, Lund, Sweden. At the BS side, a 7.5 m-long, 128-element virtual uniform linear array (ULA) was mounted on a rooftop, three floors above the ground level. Eight MS sites, four in line-of-sight (LOS) propagation conditions to the BS and four in non-LOS (NLOS), were considered. At each site, the radio channel between the ULA and a single-antenna MS was measured at five locations interspaced by 0.5 m. Measurements were acquired at a carrier frequency of 2.6 GHz, and over a bandwidth of 50 MHz. All antennas were vertically-polarized and omnidirectional on the horizontal plane.

To extract clusters of MPCs, the measured channels were subsequently processed in the following way. First, the ULA was partitioned into sets of neighboring antennas by using a 10 element-long sliding window. Then, MPC parameters (i.e., delay, azimuth angle-of-arrival, and complex amplitudes) were estimated by applying the SAGE [38] super-resolution channel estimation algorithm to each window. Finally, clusters were identified using the KPowerMeans [37] joint clustering and tracking algorithm. Clusters containing less than 2.5% of the total power in any window, or surviving less than 5 windows, were discarded. See [7] for further details.

IV-B System Model and Assumptions

Because of the nature of the measurements available, the scope of this section is limited to ULAs. However, the results obtained below can be extended to the case of two-dimensional arrays by applying, e.g., the techniques demonstrated in Sec. V-B. Such an extension is best addressed with complementary measurements and therefore we leave it for future work.

Let the ULA be located at $0\leq x_{1}\leq x\leq x_{2}$, and assume that BS-VRs can be modeled as arising from a birth-death process having the following properties:

1. A1.

Birth events (i) are independent from each other; (ii) for positive increments $h$, $\operatorname{Pr}($exactly one birth event during$\ (x,x+h])=\lambda h+o(h)$ as $h\to 0$ for some $\lambda>0$ called the birth rate; and (iii) $\operatorname{Pr}($at least two birth events during$\ (x,x+h])=o(h)$ as $h\to 0$. By [39, Chapter 8, Definition II], this is equivalent to the assertion that the birth count process $\{N_{\operatorname{birth}}(x),x\geq 0\}$, where

[TABLE]

is a Poisson process with intensity $\lambda>0$. 2. A2.

The true lifetimes of BS-VRs are (i) independent copies of some nonnegative random variable $Y$ with probability density function (PDF) $f_{Y}:\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+}$ with finite first moment; (ii) independent of birth events. 3. A3.

In experiments one can only observe (perhaps partially) BS-VRs that overlap with $[x_{1},x_{2}]$. Hence, the outcome of any experiment is a sequence of intervals $[a_{1},b_{1}],\ldots,[a_{n},b_{n}]\subseteq[x_{1},x_{2}]$, where $a_{1}\leq\ldots\leq a_{n}$ are the observed birth positions, $\upsilon_{1}=b_{1}-a_{1},\ldots,\upsilon_{n}=b_{n}-a_{n}$ the observed lifetimes, and $n$ the number of observed BS-VRs.

IV-C Distribution of the Number of BS-VRs

Let the random variable $N(x_{1},x_{2})$ denote the number of BS-VRs observed in the interval $[x_{1},x_{2}].$ Then, some natural questions to ask are: What is the expected number of observed BS-VRs, $\operatorname{E}(N(x_{1},x_{2}))$? And, how does this number depend on the interval length $x_{2}-x_{1}\geq 0$? To answer these questions, in this section we derive the distribution function of $N(x_{1},x_{2})$. We start by noting that the total number of observed BS-VRs can be written as the sum of two independent contributions:

• •

The number, $N_{\operatorname{new}}(x_{1},x_{2})$, of BS-VRs “born” in the interval $(x_{1},x_{2}]$. Clearly, this quantity can be obtained by the difference

[TABLE]

Since $N_{\operatorname{birth}}(x)$ is, by hypothesis, a Poisson process, it follows that

[TABLE]

where $\operatorname{Po}(m)$ denotes the Poisson distribution with parameter $m\geq 0$ and probability mass function (PMF)

[TABLE]

where one should assign the value 1 to $0^{0}$.

• •

The number, $N_{\operatorname{alive}}(x_{1})$, of BS-VRs born before, and still alive at, $x=x_{1}.$ The distribution of $N_{\operatorname{alive}}(x)$ for $x_{1}$ sufficiently far away from $x=0$ is given as part of the proof of Theorem 1.

The following result is a key ingredient of our proposed channel model extension for PLAs.

Theorem 1** (Number of observed BS-VRs for PLAs).**

Let the assumptions in Sec. IV-B hold. Then, the number of observed BS-VRs in the interval $[x_{1},x_{2}]$ has the distribution

[TABLE]

Proof.

By (A2) and the fact that the sum of two independent Poisson-distributed random variables is also Poisson-distributed, it suffices to show that $N_{\operatorname{alive}}(x_{1})\in\operatorname{Po}(\lambda\cdot\operatorname{E}(Y)).$

Clearly, only birth events $N_{\operatorname{birth}}(x_{1})$ that lead to BS-VRs extending past $x=x_{1}$ contribute to $N_{\operatorname{alive}}(x_{1}).$ Let $0\leq u<x_{1}.$ Then, BS-VRs emanating from a neighborhood $(u,u+h]$ of $u$ are observed at $x=x_{1}$ with probability

[TABLE]

Since $\{N_{\operatorname{birth}}(x),x\geq 0\}$ has independent increments and BS-VR lifetimes are generated independently of each other, it follows that the surviving BS-VRs may be interpreted as arising from a non-homogeneous Poisson process $\{N^{\prime}_{\operatorname{birth}}(x),x\geq 0\}$ with position-dependent intensity $\lambda(x)$ given by

[TABLE]

From the theory of non-homogeneous Poisson processes [39] we obtain

[TABLE]

Finally, assuming that the starting point of the array can be chosen as far from $x=0$ as desired produces

[TABLE]

which completes the proof. ∎

Remark 1*.*

Note that Theorem 1 gives us the answers we seek, namely that for PLAs, the expected number of observed BS-VRs, $\operatorname{E}(N(x_{1},x_{2}))$, depends on the environment through (i) the intensity $\lambda>0$ of the underlying Poisson process, and (ii) the mean $\operatorname{E}(Y)$ of the true BS-VR lifetimes; and on the system through (iii) the length $L=x_{2}-x_{1}$ of the PLA. In particular, note that the specific shape of the distribution $f_{Y}:\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+}$ does not matter.

Remark 2*.*

In this paper we are interested mainly in extensions of the COST 2100 model supporting MaMi, but it is also instructive to compare these with the “conventional” COST 2100 model, which only considers physically-compact arrays (PCAs). Observe that

[TABLE]

so that the number of VRs observed by a PCA placed anywhere along the line of the PLA is

[TABLE]

As a matter of fact, the conventional COST 2100 models the number of visible far clusters (and thus of VRs) as being Poisson-distributed with a certain prescribed intensity which depends on the simulated environment; see [13, Chapter 3]. Thus, our proposed channel model extension for PLAs naturally contains and generalizes the conventional COST 2100 model.

Fig. 2 verifies Theorem 1 empirically. In particular, the empirical cumulative density functions (ECDFs) of $N_{\operatorname{alive}}(x_{1})$, $N_{\operatorname{new}}(x_{1},x_{2})$, and $N(x_{1},x_{2})$ for the measured LOS sites described in Sec. IV-A are presented, along with Poisson fits. Empirical curves are based on 20 samples, five from each measured LOS site. A quick inspection reveals that these quantities are well approximated by Poisson-distributed random variables; the Chi-square goodness-of-fit test at significance level $\alpha=0.05$ statistically confirms this visual impression. Empirical evidence from measured NLOS sites (not shown due to lack of space) also supports the hypothesis that the number of observed BS-VRs is well modeled by (2), (4), and (7).

IV-D Maximum-Likelihood Estimation of BS-VR Lifetimes

Let $\mathcal{I}=\{1,\ldots,n\},$ and define

[TABLE]

to be the subsets of indices whose BS-VRs are either observed completely, extend into $x>x_{2}$, into $x<x_{1}$, or into both $x<x_{1},$ $x>x_{2}$, respectively. Further, let $n_{00}$, $n_{01}$, $n_{10}$, $n_{11}$ denote the number of elements in $\mathcal{N}_{00}$, $\mathcal{N}_{01}$, $\mathcal{N}_{10},$ $\mathcal{N}_{11}$ so that

[TABLE]

holds. It is shown in [40] that, under the assumptions Sec. IV-B, the likelihood function that makes use of all the available information has the form

[TABLE]

where $\boldsymbol{\theta}=[\lambda,\boldsymbol{\theta}_{Y}^{\operatorname{T}}]^{\operatorname{T}}$ is the $(p+1)$-dimensional vector of unknown deterministic parameters to be estimated, $\boldsymbol{\theta}_{Y}$ is the $p$-dimensional vector parametrizing the PDF $f_{Y}(y;\boldsymbol{\theta}_{Y})$, $\boldsymbol{x}=[a_{1},b_{1},\ldots,b_{n}]^{{\operatorname{T}}}$ is a $2n$-dimensional vector of data, $\Delta_{0}\geq 0$ is the minimum feature size, and $Y_{>\Delta_{0}}$ is the restriction of $Y$ to $\{\omega:Y(\omega)>\Delta_{0}\}$. The maximum likelihood estimator (MLE) $\hat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$ given $\boldsymbol{x}$ is

[TABLE]

For arbitrary PDFs $f_{Y}(y;\boldsymbol{\theta}_{Y}):\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+}$, closed-form solutions to (13) are in general not available and, as a result, one has to resort to numerical methods. It turns out, however, that the true BS-VR lifetimes $y_{1},y_{2},\ldots$ in our channel measurements can be well approximated by an exponential law $\operatorname{Exp}(L_{\operatorname{BS}})$, as indicated by Fig. 3.111Incidentally, an exponential law is assumed in [20, 21] to model the appearance and disappearance of clusters on both the array and time axes, but no empirical evidence is given to support this assumption. That is, the random variable $Y$ has the distribution

[TABLE]

where $L_{\operatorname{BS}}$ is the mean BS-VR lifetime. In this case, the MLE takes on a relatively simple form.

Theorem 2** (MLE for PLAs).**

Consider a birth-death process with intensity $\lambda>0$ obeying the assumptions in Sec.IV-B. Let $Y\in\operatorname{Exp}(L_{\operatorname{BS}})$ be the distribution of the true BS-VR lifetimes. Then

1. (i)

the quantities $n$, $\nu=n_{11}-n_{00}$, and $\Lambda_{0}=\sum_{i=1}^{n}(\upsilon_{i}-\Delta_{0})$ are sufficient statistics; 2. (ii)

the MLEs of the unknown parameters $\lambda,$ $L_{\operatorname{BS}}$ are given by

[TABLE]

where $L_{0}=L-\Delta_{0}$ is the array length shortened by the minimum feature size, $\Delta_{0}.$

We will need the following result.

Lemma 1**.**

Let the assumptions in Theorem 2 hold, let $\Delta_{0}\geq 0,$ and define $\lambda_{0}=\lambda e^{-\frac{\Delta_{0}}{L_{\operatorname{BS}}}}.$ Then, two PLAs given by the 4-tuples $(L,L_{\operatorname{BS}},\lambda,\Delta_{0})$ and $(L_{0},L_{\operatorname{BS}},\lambda_{0},0)$ have the same likelihood function, provided that observations $\upsilon_{1},\ldots,\upsilon_{n}$ in the former one are shortened by $\Delta_{0}.$

Proof.

Assume $(L,L_{\operatorname{BS}},\lambda,\Delta_{0})$, apply (14) to (IV-D) to obtain

[TABLE]

which is unchanged under the mapping $L\mapsto L_{0},$ $\lambda\mapsto\lambda_{0},$ $\Delta_{0}\mapsto 0$ $\upsilon_{i}\mapsto\upsilon_{i}-\Delta_{0}.$ ∎

We are now ready to prove Theorem 2.

Proof.

By Lemma 1 we can assume without loss of generality that $\Delta_{0}=0.$ Claim (i) follows directly from an application of the Neyman-Fisher factorization theorem. We now prove (ii). Taking the partial derivatives of the logarithm of the likelihood function (17), we have

[TABLE]

and setting them equal to zero produces

[TABLE]

Solving (18a) for $\lambda$ and inserting into (18b) produces

[TABLE]

which, by grouping terms of the same degree, can be rewritten as

[TABLE]

with $a=\nu-n,$ $b=L\nu+\sum_{i=1}^{n}\upsilon_{i},$ and $c=L(\sum_{i=1}^{n}\upsilon_{i}).$ Solving for $L_{\operatorname{BS}}$ yields the two solutions

[TABLE]

Clearly, we have $c\geq 0$ and, by (11), we see that $a\leq 0$. It follows that $b^{2}-4ac\geq b^{2},$ and so only the solution

[TABLE]

pertains to the permissible range $\hat{L}_{\operatorname{BS}}\geq 0$ of the parameter $L_{\operatorname{BS}}.$ Since $a\leq 0$, it is clear that $\hat{L}_{\operatorname{BS}}$ indeed maximizes the log-likelihood function. Finally, we see from (18a) that the MLE of the parameter $\lambda$ must be (15). ∎

Remark 3*.*

Note in Theorem 2 that if $n_{00}=n_{11}+\frac{\Lambda_{0}}{L_{0}}$ holds (and so $n_{11},\frac{\Lambda_{0}}{L_{0}}\in\operatorname{\mathbb{N}}^{+}$), then (16) reduces to $\hat{L}_{\operatorname{BS}}=(\frac{L_{0}\Lambda_{0}}{n-\nu})^{\frac{1}{2}}$. On the other hand, if $n_{00}=n_{01}=n_{10}=0$, then $\hat{L}_{\operatorname{BS}}=\infty.$

It is a well-known result [41, Chapter 7] that under certain regularity assumptions on the likelihood function, the MLE is asymptotically unbiased (and thus $\operatorname{E}(\hat{\boldsymbol{\theta}})\to\boldsymbol{\theta}$ as $n\to\infty$) and efficient (and thus $\operatorname{Var}(\hat{\boldsymbol{\theta}})\to\text{CRLB}$ as $n\to\infty$ as well, where CRLB means the Cramer-Rao lower bound of the PDF $f_{\boldsymbol{\upsilon}}(\boldsymbol{\upsilon};\boldsymbol{\theta}).$) It is therefore of interest to obtain the CRLB for the estimates of the intensity $\lambda$ and lifetimes $L_{\operatorname{BS}}.$

Theorem 3** (Cramer-Rao lower bound).**

Consider a birth-death process with intensity $\lambda>0$ satisfying the assumptions in Sec.IV-B, and let $Y\in\operatorname{Exp}(L_{\operatorname{BS}})$ be the distribution of the true lifetimes. Then, the Fisher information matrix (FIM) is

[TABLE]

The CRLB is found as

[TABLE]

The proof is a routine exercise in view of the equalities $\operatorname{E}(n)=\lambda_{0}(L_{\operatorname{BS}}+L_{0}),$ $\operatorname{E}(n_{11}-n_{00})=\lambda_{0}(L_{\operatorname{BS}}-L_{0}),$ $\operatorname{E}(\Lambda_{0})=\lambda_{0}L_{\operatorname{BS}}L_{0},$ and is thus omitted.222Even more is true. The equalities $\operatorname{E}(n)=\lambda(L_{\operatorname{BS}}+L),$ $\operatorname{E}(n_{11}-n_{00})=\lambda(L_{\operatorname{BS}}-L),$ $\operatorname{E}(\Lambda)=\lambda L_{\operatorname{BS}}L$ hold for arbitrary $f_{Y}:\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+},$ and $\Delta_{0}=0.$ Numerous examples can be found in, e.g., [41]. Note that the expectation is with respect to $f_{\boldsymbol{x}}(\boldsymbol{x};\boldsymbol{\theta}).$

The well-known information inequality [41] asserts that $\operatorname{E}((\hat{\boldsymbol{\theta}}-\boldsymbol{\theta})(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta})^{\operatorname{T}})\succeq\boldsymbol{I}(\boldsymbol{\theta})^{-1}$ for $\hat{\boldsymbol{\theta}}$ unbiased, from which it follows that

[TABLE]

and

[TABLE]

For fixed birth rate $\lambda>0$, the normalized Cramer-Rao bounds satisfy

[TABLE]

and so the influence of the parameter $L_{\operatorname{BS}}$ is limited. In fact, (23) reveals that the normalized CRLBs are dominated by the inverse of the intensity-aperture product $\lambda_{0}L_{0}.$ Increasing the length $L$ of the array can therefore reduce the variance of $\hat{\boldsymbol{\theta}}$ significantly. Likewise, better performance of $\hat{\boldsymbol{\theta}}$ should be expected in environments with many scatterers (clusters) such that $\lambda$ is large. We also note from (23) that the effect of the minimum feature size $\Delta_{0}$ on the CRLBs can be mostly neglected so long as $\frac{\Delta_{0}}{L_{\operatorname{BS}}}\ll 1,$ $\frac{\Delta_{0}}{L}\ll 1$ hold.

To evaluate the extent to which the estimators $\hat{\lambda},$ $\hat{L}_{\operatorname{BS}}$ given by (15), (16) approach the CRLBs (21), (22) we run Monte Carlo simulations. Fig. 4 shows the root-mean-square error (RMSE) relative to the magnitude of the estimated parameter. As baseline for comparison we use the method-of-moments estimator (MoME)

[TABLE]

The second equation in (24) follows from the identity $\operatorname{E}(\upsilon)=\frac{L_{\operatorname{BS}}L_{0}}{L_{\operatorname{BS}}+L_{0}},$ and by approximating $\operatorname{E}(\upsilon)$ by $T=\frac{\Lambda_{0}}{n}.$ MoM-based estimators are sometimes used in the literature [42].333The derivation of (24) uses (i) $Y\in\operatorname{Exp}(L_{\operatorname{BS}})$, and (ii) $\operatorname{E}(\upsilon)=\frac{L_{\operatorname{BS}}L_{0}}{L_{\operatorname{BS}}+L_{0}},$ none of which is obvious. Less sophisticated estimators are possible. One can, for instance, assume $Y$ constant (and thus knowledge of type (i) is not needed), or ignore the fact that observed BS-VR lifetimes may be truncated (and thus knowledge of type (ii) is not needed). Such simpler estimators abound in the literature, but their performance is generally poor (neglecting (i) leads to data-model mismatch, and (ii), to biased estimators). In this work we ignore these simpler estimators and focus on the MoME and, especially, the MLE. According to the simulations, the MLE performs better than the MoME. As the intensity-aperture product $\lambda_{0}L_{0}$ increases, the MoME lags behind the MLE, which steadily approaches the CRLB and becomes efficient.444We have numerically verified that both the MLE and the MoME are asymptotically unbiased as $\lambda_{0}L_{0}\to\infty.$ As can be seen, estimating $\lambda$ is easier than estimating $L_{\operatorname{BS}}$, where “easier” means here that the RMSE approaches $\sqrt{CRLB}$ more rapidly. For instance, the MLE of $\lambda$ can be considered efficient as early as $\lambda_{0}L_{0}=5$, while that of $L_{\operatorname{BS}}$ requires $\lambda_{0}L_{0}=15$ or larger to be considered within range of the CLRB. A somewhat surprising observation is that the performance of the MLE improves as the ratio $\frac{L_{\operatorname{BS}}}{L_{0}}$ increases. (By contrast, the performance of the MoME appears to be rather insensitive to the value of $L_{\operatorname{BS}}.$) This can be explained by the fact that as $L_{\operatorname{BS}}$ grows larger for fixed $\lambda_{0}L_{0}$, more BS-VR lifetimes are (partially) observed by the PLA, contributing to improve estimation accuracy. Another interesting remark is that the amount of fully observed BS-VR lifetimes is only a fraction of the total; thus, developing estimators that can reliably cope with truncated observations appears to be of great importance.

Applying Theorems 2 and 3 to the data from the LOS sites, we obtain the estimates

[TABLE]

Compared to Fig. 4, the relative RMSEs have been divided by $\sqrt{4}.$ This is because each of the four measurement sites is assumed to produce uncorrelated observations. In view of the above, when designing measurement campaigns with PLAs one should perhaps ensure

[TABLE]

or so, where $N_{0}$ is here the number of uncorrelated experiments. In the next section we discuss when two experiments can be considered uncorrelated.

IV-E Autocorrelation Function of Observed BS-VRs

Of considerable importance to characterize an stochastic process is the autocorrelation function (ACF) of the process.

Theorem 4** (Autocorrelation function).**

Let $N(x,\boldsymbol{r}),$ with $(x,\boldsymbol{r})\in\operatorname{\mathbb{R}}^{+}\times\operatorname{\mathbb{R}}^{2},$ denote the number of far clusters visible from BS array location $x$ and MS location $\boldsymbol{r},$ and define the covariance function

[TABLE]

where $m=\operatorname{E}(N(x,\boldsymbol{r}))$ is assumed constant. Then, the ACF $R_{N}(\Delta x,\Delta\boldsymbol{r})=C(\Delta x,\Delta\boldsymbol{r})/C(\boldsymbol{0})$ of the number of far clusters in the COST 2100 model is given by

[TABLE]

where $R_{N}(\Delta x)=e^{-\frac{|\Delta x|}{L_{\operatorname{BS}}}}$ is the BS-side ACF, and

[TABLE]

is the MS-side ACF, where $\chi=\cos^{-1}(\frac{|\Delta\boldsymbol{r}|/2}{R_{\operatorname{C}}})$, and $R_{\operatorname{C}}\geq 0$ is the radius of the MS-VRs.

Proof.

The factorization of the ACF into BS- and MS-side ACFs comes from a separability assumption. Let us first consider the BS-side ACF

[TABLE]

Since $R_{N}(-\Delta x)=R_{N}(\Delta x)$, we need only to consider the case $\Delta x>0$. Using the notation introduced in Sec. IV-D with $x_{1}=x$, and $x_{2}=x+\Delta x$, we can write

[TABLE]

A moment’s reflection will convince the reader that $R_{N}(\Delta x)$ can only depend upon the number $N_{11}$ of BS-VRs visible from both $x$ and $x+\Delta x$, that is

[TABLE]

Using that $N(0)\in\operatorname{Po}(\lambda\cdot\operatorname{E}(Y))$, and $N_{11}\in\operatorname{Po}(\lambda\cdot\int_{\Delta x}^{\infty}(t-\Delta x)f_{Y}(t)\,dt)$, and after some manipulations, we obtain $R_{N}(\Delta x)=1-\frac{\int_{0}^{\infty}\min(t,\Delta x)f_{Y}(t)\,dt}{\operatorname{E}(Y)},$ which for $Y\in\operatorname{Exp}(L_{\operatorname{BS}})$ becomes $R_{N}(\Delta x)=e^{-\frac{\Delta x}{L_{\operatorname{BS}}}}$. The computation of the MS-side factor $R_{N}(\Delta\boldsymbol{r})$ proceeds along the same lines and is omitted. ∎

If we consider two measurements at locations $(x,\boldsymbol{r})$, $(x+\Delta x,\Delta\boldsymbol{r})$ to be uncorrelated whenever $R_{N}(\Delta x,\Delta\boldsymbol{r})\leq{1}/{e},$ then experiments become uncorrelated if conducted at $|\Delta x|>L_{\operatorname{BS}}$ apart, or $|\Delta\boldsymbol{r}|>(1-{1}/{e})2R_{\operatorname{C}}$ apart, or a combination of both. We should also note the important fact that the ACF of the large-scale fading (LSF), $R_{\operatorname{F}}(\Delta x,\Delta\boldsymbol{r})$ is also given by Theorem 4. The reason is, intuitively, that if the clusters’ powers are assumed independent and identically distributed (which is a reasonable assumption for $|\Delta x|,$ $|\Delta\boldsymbol{r}|$ in the order of a few tens of wavelengths) and independent of the VRs, then only the proportion of common clusters affects the ACF.

IV-F Results and Validation

To gain more insight into the behavior of the proposed BS-VR extension of the model and the influence of the mean lifetime and intensity, simulation results of a PLA and $K=9$ single-antenna users in NLOS propagation conditions are presented in Fig. 5. Users are located in an outdoor environment at distance $R$ from the BS; the angular separation between adjacent MSs is one degree. At the BS side, a ULA with $M=128$ $\lambda/2$-spaced antennas spanning $L=7.5$ meters is used. All antennas are vertically-polarized and omnidirectional in the azimuth plane. Simulation results have been obtained from 100 runs using the parameters in Table II. The carrier frequency is 2.6 GHz and the bandwidth, sampled at 257 equispaced points, 50 MHz.

We study the condition number, $\kappa(\boldsymbol{H})=\frac{\sigma_{1}(\boldsymbol{H})}{\sigma_{K}(\boldsymbol{H})},$ of the $K\times M$ MU-MIMO channel matrix $\boldsymbol{H}.$ In particular, CDFs of

[TABLE]

are shown, of which several interesting remarks can be made. Firstly, $\kappa_{\operatorname{dB}}$ can decrease as the mean BS-VR lifetime, $L_{\operatorname{BS}},$ increases (for fixed $\lambda$), which is in contrast to common wisdom stipulating that spatial separability improves as more array sections are exposed to different environments. Secondly, $\kappa_{\operatorname{dB}}$ decreases as the BS-VR process birth-rate, $\lambda,$ increases for fixed $L_{\operatorname{BS}}>0.$ In this case, the reason is that larger values of $\lambda$ entail richer environments with more scatterers (clusters), which facilitates spatial separability. Thirdly, $\kappa_{\operatorname{dB}}$ decreases as the BS-MS distance, $R,$ increases for fixed MS angular separation. By the discussion in the preceding section, the MS-side LSF ACF, $R_{N}(\Delta\boldsymbol{r})$, is low for $\frac{R}{R_{\operatorname{C}}}>>1,$ thereby driving toward zero the condition number, $\kappa_{\operatorname{dB}}.$ Conversely, $R_{N}(\Delta\boldsymbol{r})$ goes to one as $\frac{R}{R_{\operatorname{C}}}\to 0$, the overall effect being a noticeable increase of $\kappa_{\operatorname{dB}}$. (This is despite the fact that spherical effects become more noticeable at close range, pulling the BS-side small-scale fading (SSF) ACF, $R_{N}(x_{1},x_{2}),$ toward zero.)

V Extension for Closely-Located Users

Using indoor measured channels, we next characterize, analyze, and validate an extension of the COST 2100 model that describes individual MPCs at the MS side in terms of birth-death processes, and which is of particular relevance for MaMi systems with closely-located users.

V-A Distribution and Intensity of MPC Lifetimes

The results of this section are based on channel measurements in a large sports hall (20 m$\times$36 m$\times$7.5 m) in the 2.6 GHz band. The radio channel between a conical monopole omnidirectional antenna moving on a 15-m-long straight-line track and a fixed cylindrical array with 128 antennas was sampled at an average rate of 400 snapshots/m. The snapshots obtained in this way were then processed using a phase-based extended Kalman filter (EKF) that was able to track the MPC parameters MPC of the propagation channel along a 15 m route. The reader is referred to [31] for further discussion of the measurements and their processing.

To analyze the lifetimes of MPCs, we adopt the modeling framework introduced in Sec. IV-B. Let $\mathcal{I}=\{1,\ldots,n\},$ where $n$ is the number of detected MPCs. The output of the EKF consists then of a sequence $a_{1}\leq\ldots\leq a_{n}$ of observed MPC birth locations such that $x_{1}\leq a_{i}\leq x_{2}$ for all $i\in\mathcal{I}$, and a sequence $\upsilon_{1},\ldots,\upsilon_{n}$ of observed MPC lifetimes555Technically, the observed MPC lifetimes $\upsilon_{1},\ldots,\upsilon_{n}$ are outcomes of the random variables $\Upsilon_{1},\ldots,\Upsilon_{n}.$ These are assumed identically distributed, and are denoted generically by $\Upsilon.$ such that $\Delta_{0}\leq\upsilon_{i}\leq L$ for all $i\in\mathcal{I}$ and for some minimum feature size $\Delta_{0}>0.$ For the particular track that we will study we have $L=15$ m, $\Delta_{0}=7.5$ cm, and $n=752.$

As before, let $\boldsymbol{\theta}=[\lambda,\boldsymbol{\theta}_{Y}^{\operatorname{T}}]^{\operatorname{T}}$ be a $(p+1)$-dimensional vector of parameters to be estimated, and let $\boldsymbol{x}=[a_{1},b_{1},a_{2},\dots,b_{n}]$ be a $2n$-dimensional vector of data. To model the MPC lifetimes we restrict attention to the following three cases:

$Y\in\operatorname{Exp}(L_{\operatorname{MPC}})$. Thus, $f_{Y}(y;\boldsymbol{\theta}_{Y})$ takes the form (14) with $\boldsymbol{\theta}_{Y}=[L_{\operatorname{MPC}}]$ for some $L_{\operatorname{MPC}}\in\operatorname{\mathbb{R}}^{+}$, and the likelihood function $\mathcal{L}(\boldsymbol{\theta};\boldsymbol{x})$ is given by (17). As this and the BS-VR problem addressed in Sec. IV-D are formally identical, it follows that the vector $\hat{\boldsymbol{\theta}}=[\hat{\lambda},\hat{L}_{\operatorname{MPC}}]^{\operatorname{T}}$ of MLEs is available in closed-form from (15) and (16) in Theorem 2. 2. 2)

$Y\in\operatorname{LN}_{10}(\mu_{\operatorname{MPC}},\sigma^{2}_{\operatorname{MPC}})$ is lognormally distributed. Thus, $f_{Y}(y;\boldsymbol{\theta}_{Y})$ takes the form

[TABLE]

with $\boldsymbol{\theta}_{Y}=[\mu_{\operatorname{MPC}},\sigma^{2}_{\operatorname{MPC}}]^{\operatorname{T}}$ for some $\mu_{\operatorname{MPC}}\in\operatorname{\mathbb{R}},$ $\sigma^{2}_{\operatorname{MPC}}\in\operatorname{\mathbb{R}}^{+}$ such that

[TABLE]

and $\mathcal{L}(\boldsymbol{\theta};\boldsymbol{x})$ is given by (IV-D). As no closed-form solution is available in this case, the vector $\hat{\boldsymbol{\theta}}=[\hat{\lambda},\hat{\mu}_{\operatorname{MPC}},\hat{\sigma}^{2}_{\operatorname{MPC}}]^{\operatorname{T}}$ of MLEs has been obtained by solving (13) numerically.666To facilitate the numerical evaluations in (13), closed-form expressions, as a function of $\Phi(x),$ of some of the terms involved are given in [40]. 3. 3)

$\Upsilon\in\operatorname{TLN}_{10}(\mu_{\operatorname{MPC}},\sigma^{2}_{\operatorname{MPC}},\Delta_{0},L),$ the observed MPC lifetimes, follows a truncated lognormal distribution. In this case $f_{\Upsilon}(\upsilon;\boldsymbol{\theta}_{\Upsilon})$ is given by

[TABLE]

with $\boldsymbol{\theta}_{\Upsilon}\!=\![\mu_{\operatorname{MPC}},\!\sigma^{2}_{\operatorname{MPC}}]^{\operatorname{T}}$, and $\mu_{\operatorname{MPC}},\sigma^{2}_{\operatorname{MPC}},f_{Y}$ as in case 2. The vector $\hat{\boldsymbol{\theta}}_{\Upsilon}\!=\![\hat{\mu}_{\operatorname{MPC}},\!\hat{\sigma}^{2}_{\operatorname{MPC}}]^{\operatorname{T}}$ has been selected to minimize the vertical distance between the empirical CDF, obtained from channel measurements, and $F_{\Upsilon}(\upsilon;\boldsymbol{\theta}_{\Upsilon})$, the reference CDF, with $\boldsymbol{\theta}_{\Upsilon}$ in some range of interest; $\hat{\lambda}$ has simply been set to $\hat{\lambda}=n/L.$ We note that this case seems to be the most common approach in the literature [30, 31].

Fig. 6 shows a comparison of the ECDF obtained from measurements, and the CDFs corresponding to cases 1), 2), and 3). According to the figure, the best fit is the lognormal one (case 2). An exponential fit (case 1) matches measurements poorly, underestimating the cumulative density function for short- and medium-lived MPCs ($\upsilon\leq 2$ m) and overestimating it for long-lived ones ($\upsilon>2$ m). The lognormal fit of case 2, on the other hand, shows an excellent agreement over the entire range of the observations. Even the proportion of lifetimes exceeding $L$, namely $\operatorname{Pr}(Y>L)$, is best explained by the lognormal distribution, as shown in the inset. Turning to the truncated lognormal case, the figure shows that when observed ($\Upsilon$) rather than true ($Y$) MPC lifetimes are modeled, the resulting fit may not be satisfactory even when the “correct” (lognormal) distribution is being used (see also [30, 31]). The problem is that the observation model of case 3) takes observations $\upsilon_{1},\ldots,\upsilon_{n}$ at face value, and neglects the fact that MPCs may extend beyond the measurement track, at one or both ends. The data shown in Table I further emphasizes this point.

V-B Distribution of MPC-VR Radii

In order to model the lifetime of individual MPCs within a cluster we introduce the concept of MPC visibility region (MPC-VR), which is schematically illustrated in Fig. 1c in Sec. III. In the figure, MPC-VRs are represented by discs contained within a certain MS-VR. There is a one-to-one correspondence between the MPC-VRs of a certain MS-VR, and the scattering points of the associated cluster. When the MS enters a MPC-VR, the corresponding scattering point becomes active and contributes to the propagation channel. To give an example, active MPC-VRs in Fig. 1c have been marked gray. It is a simple but crucial observation that different MSs do not necessarily “see” the same set of active MPCs, even when they do have the same set of active MS-VRs; we will apply this later to establish that the SSF ACF decays more rapidly than when the concept of MPC-VR is not applied.

Let a MS move along the $X$-axis, and let MPC-VRs be distributed uniformly on the $XY$-plane. Then, the distribution of the MPC lifetimes encountered during the movement (i.e., the length of the intersecting chords) is

[TABLE]

where $R_{\operatorname{MPC}}\in\operatorname{\mathbb{R}}^{+}$ is the radius of the MPC-VRs. (For simplicity, we have assumed that the same set of MS-VRs is always visible during the movement. This assumption is often made in indoor environments [43, 14].) The CDF (33) has been plotted in Fig. 7, top, for several values of $R_{\operatorname{MPC}}.$

To reconcile (33) with our previous findings, namely that true MPC lifetimes appear to be lognormally distributed, we need to find $R_{\operatorname{MPC}}\in\operatorname{\mathbb{R}}^{+}$ such that $\breve{F}_{Y}(y;R_{\operatorname{MPC}})$ approaches

[TABLE]

with $\hat{\boldsymbol{\theta}}_{Y}$ given by case 2 in Table I, and $m,\psi$ as before. However, a quick examination of Fig. 7, top, shows that no such number exists: The desired lognormal curve (34) is concave for $y\geq 0$, but for every fixed $R_{\operatorname{MPC}}\in\operatorname{\mathbb{R}}^{+}$, (33) turns out to be convex in the interval $0\leq y\leq 2R_{\operatorname{MPC}}.$

The solution to this difficulty is to use discs of various sizes, an idea which is supported by the measurements. Let $R_{\operatorname{MPC}}$ be a random variable with PDF $f_{R_{\operatorname{MPC}}}:\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+}$. Then by (the continuous version of) the law of total probability [39] we have

[TABLE]

Thus, $\breve{F}_{Y}(y)$ can be thought of as a “weighted sum” of MPC-VR functions (33), and our task is to find nonnegative weights $f_{R_{\operatorname{MPC}}}(r)$ for all $r\in\operatorname{\mathbb{R}}^{+}$ such that $\breve{F}_{Y}(y)$ is close to lognormal according to some prescribed criterion of closeness. Our approach here is to adopt the least-squares criterion. As we have no a priori reason to expect $f_{R_{\operatorname{MPC}}}(r)$ to be easily obtainable by analytical methods, we solve a discretized version of problem (35), namely

[TABLE]

where $\boldsymbol{A}$ is a $p\times q$ matrix with $[\boldsymbol{A}]_{ij}=\breve{F}_{Y}(y_{i};r_{j})$ as the $ij$-th entry, $\boldsymbol{w}$ is a $q\times 1$ vector of nonnegative weights $w_{1},\ldots,w_{q}$, and $\boldsymbol{b}$ is a $p\times 1$ vector with $F_{Y}(y_{i};\boldsymbol{\theta})$ as the $i$-th entry. Further, we have defined the sets $\mathcal{Y}=\{y_{1},\ldots,y_{p}\}$ of sampling lifetimes, and $\mathcal{R}=\{r_{1},\dots,r_{q}\}$ of sampling radii. Note that the condition $\sum_{i=1}^{q}w_{i}=1$ is needed to ensure that the mapping

[TABLE]

is a proper PMF. Crucially, the optimization problem (V-B) is a quadratic program [44], and can be readily solved by using a number of available numerical packages, such as MATLAB or CVX [45].777A different version of (V-B) sets $[\boldsymbol{A}]_{ij}=r_{i}\,\breve{F}_{Y}(y_{i};r_{j}),$ $1\leq i\leq p,$ $1\leq j\leq q,$ and uses the mapping $p_{R_{\operatorname{MPC}}}(r_{i})=(\frac{w_{i}}{r_{i}})/(\sum\limits_{k=1}^{p}\frac{w_{k}}{r_{k}})$ instead of (37). While closer to the original, continuous formulation, this alternative way of discretizing (35) is prone to numerical issues and instability for small values of $r_{i}.$ Because of this, in this work we restrict attention to (V-B), (37).

Fig. 7, bottom left, shows a comparison of the lognormal CDF (34) and its approximation $\tilde{F}_{Y}(y_{i})$ obtained by solving (V-B) with $\mathcal{Y}=\{.0025,\allowbreak.0525,\allowbreak\ldots,\allowbreak 14.9525\}$ and $\mathcal{R}=\{.000,\allowbreak.025,\allowbreak\ldots,\allowbreak 23.000\}$, so that $\boldsymbol{A}$ is of size $300\times 921.$ The approximation is excellent, with a root-mean-square error

[TABLE]

about $1.1\times 10^{-14}$ in the range $0\leq y\leq L.$ The goodness of the fit suggests that the proposed method, namely to model MPC birth-death processes using MPC-VRs with radii $R_{\operatorname{MPC}}$ drawn from a certain distribution $f_{R_{\operatorname{MPC}}}:\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+}$, is rather flexible and, as such, applicable to a greater range of environments having different propagation characteristics. For each such environment, the appropriate distribution of $R_{\operatorname{MPC}}$ can be found by solving an optimization problem like (V-B).

In our case, however, a further generalization is possible: Fig. 7, bottom right, shows the cumulative sum of the weights $w_{1},\ldots,w_{921}$, and a lognormal approximating curve which has the form (34) with parameters

[TABLE]

The radii $R_{\operatorname{MPC}}$ are now drawn from a two-parameter continuous distribution which conveniently summarizes all $q=921$ weights. Again, the agreement between the two curves is excellent, with a $\operatorname{RMSE}\approx 2.2\times 10^{-4}$ whenever $R_{\operatorname{MPC}}>.025.$ Motivated by this, we introduce in the COST 2100 model circular MPC-VRs with lognormally distributed radii according to the parameters $\mu_{R_{\operatorname{MPC}}},$ $\sigma^{2}_{R_{\operatorname{MPC}}}$, given in Table II.

V-C Distribution of MPC-VR Amplitudes

We model the relative contribution of each MPC to the total channel gain by means of the so-called gain function [7, 1]. More precisely, a Gaussian profile

[TABLE]

multiplies the complex amplitude of each MPC, where the weight $g_{{\operatorname{MPC}},\ell}(\boldsymbol{r}_{\operatorname{MS}})$ depends on the Euclidean distance between the user position $\boldsymbol{r}_{\operatorname{MS}}=[r_{{\operatorname{MS}},x},\allowbreak r_{{\operatorname{MS}},y},\allowbreak r_{{\operatorname{MS}},z}]^{\operatorname{T}}$ and the center $\boldsymbol{r}_{{\operatorname{g}},\ell}=[r_{{\operatorname{g}},\ell,x},r_{{\operatorname{g}},\ell,y}]^{\operatorname{T}}$ of the $\ell$-th MPC-VR measured on the $XY$-plane:

[TABLE]

Thus, at distance, e..g., $\sigma_{{\operatorname{g}},\ell}$ from its center, the MPC gain function has decayed by 4.3 dB. The centers $\boldsymbol{r}_{{\operatorname{g}},\ell}$ are generated uniformly inside the corresponding MS-VRs, which are assumed to have radius $R_{\operatorname{C}}\in\operatorname{\mathbb{R}}^{+}$. The profile widths $\sigma_{{\operatorname{g}},\ell}$ are identified with the radii $R_{\operatorname{MPC}}$, and follow a lognormal distribution with parameters see (39). Thanks to the MPC gain function feature, variations observed [46, 26, 31] in the gain of individually tracked MPCs can be adequately modeled, helping to smoothly ramp up (fade out) the contributions of new-born (dead) MPCs, as users navigate in an out of the MPC-VRs.

Let $N_{\operatorname{MPC}}^{\operatorname{eff}}$ denote the average effective number of MPCs per cluster, defined as the average number of MPCs whos amplitude at any given point within the MPC-VR exceeds a certain value. Then, the required number of MPCs per cluster is

[TABLE]

where

[TABLE]

and $\log$ denotes the natural logarithm. We remark that if backwards compatibility is desired, one can simply set $N_{\operatorname{MPC}}=N_{\operatorname{MPC}}^{\operatorname{eff}},$ $\mu_{R_{\operatorname{MPC}}}=\infty,$ $\sigma_{R_{\operatorname{MPC}}}=0,$ and the MPC-VRs concentric with the MS-VRs.

V-D Results and Validation

To validate the proposed extension, we here study the autocorrelation function and the channel condition number.

V-D1 Autocorrelation Function

One important result in this paper is the following:

Theorem 5** (ACF of SSF).**

The space-frequency ACF of the SSF is asymptotically given by

[TABLE]

where $R_{\operatorname{H}}(\Delta\boldsymbol{r},\Delta f)$ is the ACF of the SSF without the gain function, and

[TABLE]

In particular, $R(\Delta\boldsymbol{r},\Delta f)$ obeys the bound

[TABLE]

Proof.

We start by writing the propagation channel of the COST2100 model as

[TABLE]

where $\boldsymbol{k}_{\ell}=\frac{2\pi}{\lambda_{0}}\hat{\boldsymbol{k}}_{\ell}$ is the (vector) wavenumber, $\hat{\boldsymbol{k}}_{\ell}$ is a unit vector pointing in the direction of propagation of the $\ell$-th MPC, and we assume that the complex amplitudes $a_{\ell}$ (and thus $H(\boldsymbol{r},f))$ have zero mean. Note that the subindices $\ell$ can possibly extend over several MPC-VRs. For any $\boldsymbol{r}\in\operatorname{\mathbb{R}}^{3},f\in\operatorname{\mathbb{R}}$, the covariance function of the SSF is given by

[TABLE]

Now, by defining $\mathcal{N}_{11}:=\mathcal{N}(\boldsymbol{r})\cap\mathcal{N}(\boldsymbol{r}+\Delta\boldsymbol{r})$, and assuming uncorrelated scattering, (V-D1) reduces to

[TABLE]

where $C_{H}(\boldsymbol{r},f)$ is the covariance function of the SSF without the gain function. Finally, normalize $C(\Delta\boldsymbol{r},\Delta f)$ by dividing by $C(0,0)=C_{H}(0,0)$ and the claim follows. ∎

Example 1**.**

Let $R_{\operatorname{MPC}}$ have distribution $F_{R_{\operatorname{MPC}}}(r)=\delta_{R_{0}}(r)$ for some $R_{0}\geq 0$, where $\delta_{x_{0}}(x)$ is the Dirac delta function. Then the ACF of the number of MPCs has the form

[TABLE]

where $\chi=\cos^{-1}(\frac{d/2}{R_{0}})$ (and so $0\leq\chi\leq\frac{\pi}{2}$), and $d=|\Delta\boldsymbol{r}|$ (cf. BS-side ACF of Theorem 4). The autocorrelation function (48) is sometimes called the circular correlation function [47, 48]. Note that (48) vanishes beyond $d=R_{0}$, and so the SSF becomes uncorrelated.

Example 2**.**

Let $R_{\operatorname{MPC}}$ have arbitrary PDF $f_{R_{\operatorname{MPC}}}:\operatorname{\mathbb{R}}^{+}\to\operatorname{\mathbb{R}}^{+}$ with finite second-order moment. Then the ACF of the number of MPCs takes the form

[TABLE]

with $\chi(r)=\cos^{-1}(\frac{d/2}{r})$ and, consequently, $\sin(2\chi(r))=2(\frac{d/2}{r})\sqrt{1-(\frac{d/2}{r})^{2}}.$

The examples above show that the MPC gain function can play an important role at reducing multiuser correlation, especially when the channel as seen by the users exhibits small angular spreads, or when the MPC lifetimes are particularly short (see Example 1). For typical angular spreads in outdoor and (especially) indoor environments, however, $R_{\operatorname{H}}(\Delta\boldsymbol{r},\Delta f)$ often decays rapidly with $|\Delta\boldsymbol{r}|$, in effect rendering the contribution of $R_{Y}(\Delta\boldsymbol{r})$ less apparent. In fact, a careful comparison888Due to lack of space, the detailed results of this comparison are omitted, and only a short summary is given. See towards the end of the next section for related comments. of user correlation levels in (i) channels measured in an indoor environment, (ii) synthetic channels emulating the said environment and including the proposed gain function extension, and (iii) synthetic channels without the gain function extension, revealed no statistically significant differences between the three cases ($|\Delta\boldsymbol{r}|$ being in the range between 0.5 and 3.5 m). Notwithstanding the above, the importance of modeling the MPC gain function becomes clear when higher-order statistics, involving more than two users, are considered, as we show next.

V-D2 Channel Condition Number

It turns out that the MPC gain function has quite a noticeable effect on the condition number $\kappa(\boldsymbol{H}(t,f)),$ where $\boldsymbol{H}(t,f)$ is the $K\times M$ channel matrix at time $t$ and frequency $f$ of an MU-MIMO system with $K$ single-antenna users and $M$ BS antennas. In Fig. 8, a comparison of simulated MU-MIMO channels with and without the MPC gain function is presented. In particular, CDFs of the logarithmic channel condition number

[TABLE]

of $K=9$ closely-located, single-antenna users in a LOS, public auditorium [19, indoor site MS 1] are shown. We consider both directive and omnidirectional antenna patterns. Directive antenna patterns have been measured with an upper body phantom [49] and in our simulations are oriented towards the auditorium stage. At the BS side, a compact cylindrical array with measured antenna patterns have been used. We consider $M=32$ and $128$ BS antennas. The positions and spreads of clusters are set as in the actual measured environment, while the rest of the COST 2100 model parameters have been extracted from measured indoor channels [50] and are reported in Table II. Simulation results are based on data from 10 runs. For each run, $B=257$ frequency points over 50 MHz of bandwidth, and $T=10$ positions (i.e., snapshots) over 0.25 m of a straight line, one for each user, have been simulated. For comparison, the CDFs obtained from measured channels are also included (black, dotted line). Prior to computing $\kappa_{\operatorname{dB}}$, the complex gains of both simulated and measured channels have been normalized by the expression

[TABLE]

where $h_{km}^{\operatorname{meas}}(t,f)$ are the unnormalized channels, so that the total channel gain for each user is $MTB.$

As can be seen, when the MPC gain function is included in the simulations (“ON” label), a good agreement is observed between simulated and measured channels for both 32 and 128 BS antennas. On the other hand, not including it (“OFF” label) results in too large values of $\kappa_{\operatorname{dB}}$ by about 10 dB. This gap between the “ON” and “OFF” cases is reduced by 2 dB or so when measured rather than omnidirectional antenna patterns are applied.999This can be explained as follows. When the MPC gain function is turned “OFF” directive antennas can, by selecting different MPCs, help decorrelating the users’ radio channels. When the gain function is turned “ON” users already “see” different MPCs and directive antenna patterns, by filtering away some of the MPCs, may have an adverse effect. This improvement alone, however, is not sufficient to explain the observed low values of $\kappa_{\operatorname{dB}}$ measured, and the introduction of a MPC gain function seems to be required for realistic channel simulations.

To gain further insights into the behavior of the MPC gain function in multiuser systems, we next study the dependence of $\kappa_{\operatorname{dB}}$ on the angular spreads of the environment. For this we adopt a simple propagation model consisting of only one cluster whose MS- and BS-side angular spreads, denoted respectively by $\Omega_{\operatorname{MS}}$ and $\Omega_{\operatorname{BS}}$, can be controlled independently. (This is called a “twin” cluster in the terminology of the COST 2100 model [13]). Fig. 9 shows the average

[TABLE]

of $\kappa_{\operatorname{dB}}$ for various MS- (left) and BS-side (right) angular spread values. Here, $T=300,$ and for each $1\leq t\leq T$ the locations of the users have been randomly assigned inside a $2~{}$m radius circle. Furthermore, we have set the MPC-VR radius $R_{\operatorname{MPC}}=0.5,$ the effective number of MPCs $N_{\operatorname{MPC}}^{\operatorname{eff}}=100,$ the number of users $K=9,$ and the number of BS antennas $M=128.$

In general, the channel condition number, $\bar{\kappa}_{\operatorname{dB}}$, decreases as the angular spreads increase. When the MPC gain function is turned “ON”, $\bar{\kappa}_{\operatorname{dB}}$ always attains a lower value. Moreover, this value mainly depends on the BS angular spread $\Omega_{\operatorname{BS}}$, and not so much on the MS angular spread $\Omega_{\operatorname{MS}}$. When, on the other hand, the MPC gain function is turned “OFF”, we see that $\bar{\kappa}_{\operatorname{dB}}$ depends on both $\Omega_{\operatorname{BS}}$ and $\Omega_{\operatorname{MS}}$. This has interesting implications for the design and validation of MaMi systems. Most notably, the plots in Fig. 9 suggest that users having angular spreads as small as $\Omega_{\operatorname{MS}}=5^{\circ}$ can still be served concurrently provided that $\Omega_{\operatorname{BS}}$ and $M$ are sufficiently large. Conversely, if the decorrelating effects of the MPC gain function are not properly accounted for, results involving closely-located users may be overly pessimistic, even for large $\Omega_{\operatorname{BS}}$ and $M$.

The significance of the MPC gain function will depend on the spatial density of concurrently served users, that is, the number of users communicating over the same time-frequency resource in a certain limited region. In the bottom pane of Fig 9, plots of the gap

[TABLE]

between the “ON” and “OFF” states are presented for $\Omega_{\operatorname{BS}}=60^{\circ},$ $\Omega_{\operatorname{MS}}=15^{\circ},30^{\circ},60^{\circ},$ and $K=2,\ldots,18.$ As can be seen, the gain function has a negligible effect for low values of $K$, but its impact becomes noticeable as $K$ increases, in good agreement with our comments in Sec. V-D1.

VI Summary and Conclusions

Certain important aspects of the MaMi propagation channel are not captured by prevalent channel models. In this paper, we extend the conventional COST 2100 model with the concepts of (1) BS-VRs, to model the appearance and disapperance of clusters along the axis of a physically-large array; and (2) MPC-VRs and MPC gain function, to model, at the MS side, birth-death processes of individual MPCs. Based on measurements of MaMi channels, the parameters of the proposed extensions have been statistically characterized, and their impact on the properties of MaMi propagation channels investigated theoretically and by simulations. In addition, it has been shown that simulation results are consistent with measurements of closely-located users, which demonstrates that the proposed extensions are able to capture the intended MaMi characteristics. A MATLAB implementation of the COST 2100 model with the proposed extensions (as well as full support of other essential, but not discussed in this paper, MaMi characteristics such as spherical wavefronts and cluster dispersion in elevation) is freely available at [35]. Most importantly, this implementation provides a complete simulation framework for studies of MaMi channels for 5G and beyond.

Acknowledgment

The presented investigations are based on data obtained in measurement campaigns performed by Xiang Gao, Fredrik Tufvesson, Ove Edfors, Tommy Hult, and Meifang Zhu, and Sohail Payami. We greatly benefited from discussions with Fredrik Rusek. This research has been conducted under the auspices of the COST-IRACON action.