# Geometry-Based Stochastic Model and Statistical Characteristic Analysis of Cell-Free Massive MIMO Channels

Authors: Chengjian Liao; Kui Xu; Xiaochen Xia; Wei Xie; Nan Sha; Lihua Chen; Ningsong Liu
Journal: Security and Communication Networks
 (2022)

---

## Abstract

Cell-free massive multi-input multi-output (MIMO) systems exhibit many characteristics that differ from those of traditional centralized massive MIMO systems, and there are still research gaps in the modeling of cell-free massive MIMO channels. In this paper, a geometry-based stochastic model (GBSM) that combines the double-ring model and hemisphere model to comprehensively consider the distribution of scatterers in the environment was proposed for cell-free massive MIMO channels. Combined with the line-of-sight (LoS) path, single scattering path, and double scattering path components, the channel matrix between the access points (APs) and the user was derived. The proposed model fully considers geometric parameters such as the arrival/departure direction, elevation angle, time delay, and distance, which can accurately characterize the channel. Then, we proved that the traditional channel model, standard block-fading model, and spatial basis expansion model (SBEM) adopted in a cell-free massive MIMO system could not describe the nonstationarity in the space, time, and frequency domains, whereas the proposed GBSM could. Statistical characteristics of the channel were analyzed, including the space cross-correlation function (CCF), time autocorrelation function (ACF), Doppler power spectral density (PSD), level crossing rate (LCR), and average fade duration (AFD). Then, we investigated the proposed model by simulating the space CCF, time ACF, Doppler PSD, LCR, and AFD under the conditions of two different scatterer densities. Through simulations and analyses, some new features of cell-free massive MIMO channels were identified, providing a theoretical basis for in-depth research on cell-free massive MIMO systems. Finally, the measurement-based scenario and the WINNER II channel model are compared to demonstrate that the GBSM is more practical to characterize real cell-free massive MIMO channels.


---

## Body

## 1. Introduction

As one of the key technologies of sixth-generation (6G) mobile communication systems, cell-free massive multi-input multi-output (MIMO) technology is an ideal mobile communication method to further improve spectral efficiency [1, 2]. Cell-free massive MIMO technology introduces a “user-centric” architecture to achieve continuous coverage of a large area, thus breaking through the traditional cellular coverage idea of “cell-centric”, which effectively eliminates the limitations caused by the cell boundary and the static domination of frequency resources. In theory, when the number of antennas on the access point (AP) side is large enough, the total capacity of the system can increase linearly with the increase of the total number of users in the coverage area, which is an ideal mobile communication method to further improve the spectral efficiency.At present, research on cell-free massive MIMO systems has mainly focused on spectral efficiency and energy efficiency [3–9], system performance analyses [10–14], power and resource optimization [15–18], pilot allocation and channel estimation [2, 19, 20], user achievable rate [21–23], and conjugate beamforming [24–26]. In all the abovementioned analyses, the cell-free massive MIMO channel was mainly modeled using two model types. The first type is the standard block-fading channel model [2–8, 10–26], which consists of channel coefficient, large-scale fading including shadow fading and path loss, and small-scale Rayleigh fading obeying an independent and identically standardized normal distribution. It is assumed that small-scale fading is static during each coherence interval and that large-scale fading changes slowly. The second type is the spatial basis expansion model (SBEM) [9], which consists of the complex gain of the multiple paths that represent the small-scale Rayleigh fading and large-scale fading, as well as the array steering vector with the assumption of parallel waves and ignoring the elevation angle parameter. However, the standard block-fading channel model cannot reflect the nonstationarity of the channel. The assumption of parallel waves in the SBEM is inaccurate since the condition that the distance from the transmitting antenna to the scatterer cluster be far greater than the Rayleigh distance will no longer hold in cell-free massive MIMO systems. Moreover, neither of the abovementioned two models can accurately describe the statistical characteristics of cell-free massive MIMO channels. There are still research gaps in the modeling of cell-free massive MIMO channels, and accurate and efficient channel models are urgently needed for further system design and performance evaluation.In a cell-free massive MIMO system, multiple APs are randomly deployed in a coverage area and the distance and scattering environment between users and each AP are different, as shown in Figure1. The channel modeling assumptions of parallel wave incidence in the far-field and the common parameters between antennas [27, 28] (direction of arrival, direction of wave departure, number of scatterers/clusters, and statistical characteristics of random fading) of traditional massive MIMO systems no longer hold. Compared with a traditional centralized massive MIMO system, cell-free massive MIMO system channel modeling and information theoretical analysis face many practical problems, which are listed as follows.The random layout of APs increases the difficulty of accurately characterizing a cell-free massive MIMO channel. The assumptions of far-field parallel wave incidence, generalized channel stability, and common-parameter characteristics no longer hold. The phase of the channel steering vector no longer has a linear relationship with the antenna number. The angle extension is significantly increased due to the shortening of the user-AP distance, and the number of co-scattering clusters is significantly reduced. Additionally, it is necessary to accurately characterize and model nonstationary processes such as the time variation in the channel caused by the movement of users/scatterers.The near-field propagation and environmental evolution characteristics caused by distributed array antennas make the statistical characteristics of cell-free massive MIMO channels significantly different from those of traditional centralized massive MIMO systems. It is necessary to analyze the channel model of a cell-free massive MIMO system from the perspective of multidimensional integration of a time-space-delay antenna to reveal the sparsity and stability characteristics of the channel in a specific domain (space, time, delay, or Doppler) and to explore the common structural characteristics (time and frequency) and spatial distribution structural characteristics of the channel in typical scenarios.Figure 1
Signal propagation environment of a cell-free massive MIMO system.Considering the practical problems of cell-free massive MIMO systems, accurate and efficient channel models urgently need to be established. According to the modeling approach, channel models can be classified into deterministic channel models and stochastic channel models. Deterministic channel models rely on precise information of the propagation environment, while stochastic channel models describe channel parameters using certain probability distributions. Generally, stochastic channel models are mathematically tractable and can be adapted to various scenarios with relatively low accuracy. Among all the stochastic channel models, the geometry-based stochastic model (GBSM) assumes that scatterers are stochastically distributed according to a particular underlying geometry with a regular shape, such as twin-cluster, ellipse, multiring, twin-multiring, hemisphere, sphere, cylinder, ellipsoid, and semiellipsoid. GBSM has been widely used in a variety of scenarios, such as massive MIMO communications [29–31], vehicle-to-vehicle (V2V) communications [32–34], high-speed train (HST) communications [35, 36], and millimeter wave (mmWave) communications [37–39]. Channel modeling research on cell-free massive MIMO systems is still in its infancy, and it is not known which channel model can accurately describe the channel characteristics of cell-free massive MIMO channels. Since the GBSM can accurately describe the characteristics of the channels in space, time, and delay domains through various parameters, such as geographic location, elevation angle, azimuth angle, and time delay, this paper intends to use the GBSM to establish cell-free massive MIMO channels. Note that the references mentioned above considered the perfect CSI; thus, the perfect CSI was also considered in ours paper. The major contributions and novelties of this paper are summarized as follows:Based on the assumption of a distribution of double rings and hemispheres of scatterers, a GBSM is established for cell-free massive MIMO channels. Compared with the previous models for cell-free massive MIMO channels, the innovation of the GBSM is that it combines both two-dimensional (2D) and three-dimensional (3D) characteristics, which gives a more accurate description of the channels.The proposed model fully considers the actual channel propagation environment of a cell-free massive MIMO by combining geometric parameters such as the arrival/departure direction, elevation angle, and distance. Based on the assumption of spherical wavefronts, each link from different APs to a specific scatterer cannot be simply approximated as parallel. However, the proposed GBSM can calculate the transmission distance, delay, and Doppler shift of each link. Thus, we can accurately characterize the channel.For the proposed GBSM, we have proved that the traditional channel model, standard block-fading model, and SBEM adopted in a cell-free massive MIMO system could not describe the nonstationarity in the space, time, and frequency domains, whereas the proposed GBSM could. Statistical characteristics are derived and thoroughly investigated, such as the space cross-correlation function (CCF), time autocorrelation function (ACF), Doppler power spectral density (PSD), level crossing rate (LCR), and average fade duration (AFD).We investigate the proposed model by simulating the space CCF, time ACF, Doppler PSD, LCR, and AFD under the conditions of two different scatterer densities, as well as the influence of antenna spacing and user azimuth angle on the space CCF and time ACF. Through the simulation results and analysis, some new features of cell-free massive MIMO channels are identified.To validate the correctness of the proposed GBSM, the measurement-based scenario and the WINNER II channel model are compared. Simulation results demonstrate that the GBSM was more practical to characterize real cell-free massive MIMO channels.The rest of this paper is organized as follows. Section2 proposes the GBSM for cell-free massive MIMO channels. The statistical characteristic analysis including the space CCF, time ACF, Doppler PSD, LCR, and AFD are carried out in Section 3. Section 4 presents statistical characteristic simulation investigations, and Section 5 verifies the correctness and practicality of the model. The conclusions are drawn in Section 6.


## 2. GBSM for Cell-Free Massive MIMO Channels

### 2.1. GBSM

We consider a cell-free massive MIMO system withQ APs. Suppose that Q APs are randomly deployed in a certain city area, all of which are equipped with a uniform linear array (ULA). Each AP is equipped with P antennas, and the distance between each antenna element is δ. The far-field plane wavefront assumption in conventional massive MIMO systems is not fulfilled in cell-free massive MIMO systems, as the condition that the distance from the transmitting antenna to the scatterer cluster being far smaller than the Rayleigh distance cannot be upheld. As a result, the wavefront emitted from the n-th cluster to the receive array is assumed to be spherical. Then, the angle of arrival (AoA) is no longer linear along the array, and it needs to be computed based on geometrical relationships. Considering the channel between any AP and the user, the signal sent from the AP transmitter is scattered through the line-of-sight (LoS) path and the non-line-of-sight (NLoS) path to the user. The types of scatterers in the NLoS path can be divided into three categories: scatterers around APs, scatterers around users, and scatterers such as tall buildings and trees. We fully consider the geographic location of the APs and users and the geometric relationship of the scatterer distribution and construct a GBSM with the spherical wave assumption for cell-free massive MIMO channels, as shown in Figure 2. Some of the main parameters and definitions are shown in Table 1. The GBSM needs to be extended for different scenarios to effectively improve the accuracy of channel modeling, that is, the combination of the double-ring and hemisphere models. The reason for adopting the double-ring and hemisphere models is as follows. This paper considers a cell-free massive MIMO system model in an urban rich scattering scenario. In order to distinguish between near-end scattering and far-end scattering, this paper adopts a GBSM combined with a double-ring model and a hemisphere model. When the user faces a single AP with a short distance, there are fewer channels reflected by the AP through the higher scatterers, which can be ignored, and the channel with the AP is modeled as a single ring. This is similar to the user and the single ring model is adopted. For the channels between the user and other APs that are far away, the reflection effect of the higher scatterers needs to be considered to form a hemispherical model. Since APs are randomly distributed in the system, some APs are close to the user and the direct path is strong at this time. The scattering of high scatterers can be ignored and only the single ring model is considered. Some APs far away from the user have a weak direct path and can be ignored, while the scattering of high scatterers cannot be ignored, so only the hemispherical model is considered.Figure 2
GBSM for cell-free massive MIMO channels.Table 1
Definition of parameters.
ParametersDefinitionλCarrier wavelengthfcCarrier frequencyδAntenna spacingvuSpeed of the userβuAzimuth angle of the user’s movement directionKRician factorfDDoppler frequencyfRmaxMaximum Doppler frequency with respect to the userRq, RDistribution radii of scatterers around the AP and user, respectivelyPii=1,2,3Normalized power coefficient of the single scattering pathPDNormalized power coefficient of the double scattering pathφni, φn1,n2Initial random phase of the LoS for single and double scattering paths, respectivelyDqpDistance between theq-th AP p-th antenna and userDqp,niTi=1,2,3Distance between theq-th AP p-th antenna and the ni cluster of scatterers through the single scattering pathDqp,niRi=1,2,3Distance between theni cluster of scatterers and the user through the single scattering pathDqp,n1n2Distance between then1 and the n2 clusters of scatterers through the double scattering pathαqpT, βqpTAAoD/EAoD between theq-th AP p-th antenna and the user through the LoS pathαqpR, βqpRAAoA/EAoA between theq-th AP p-th antenna and the user through the LoS pathαqp,niT, βqp,niTAAoD/EAoD between theq-th AP p-th antenna and the ni cluster of scatterers through the single scattering pathαqp,niR, βqp,niRAAoA/EAoA between theni cluster of scatterers and user through the single scattering pathαqp,n3, βqp,n3Elevation angle and azimuth angle between the hemisphere center and then3 cluster of scatterersτqpTime delay between theq-th AP p-th antenna and user through the LoS pathτqpniTime delay between theq-th AP p-th antenna and the ni cluster of scatterers through the single scattering pathτqp,n1,n2Time delay between then1 and the n2 cluster of scatterers through the double scattering pathFor scatterers around APs and users where the elevation angle/depression angle changes are relatively small, we assume that both categories are distributed on a single ring. Since the scatterers around tall buildings and trees are relatively concentrated, the 3D distribution characteristics of APs are particularly prominent. Hemisphere models need to be used to model cell-free massive MIMO channels. We take the first antenna element of the AP transmitter and the user as the centers to establish a double-ring model with radii ofRq and R and establish a hemisphere model with the two centers as diameter Dq. We set the layout direction of the ULA to the Y axis. Suppose that N1 scatterers are distributed on a circle centered on the AP and that N2 scatterers are distributed on a circle centered on the user. The moving speed of the user is vu, and the moving azimuth angle is βu. Suppose that N3 stationary scatterers are distributed in the hemisphere area. Note that the GBSM proposed in this paper includes 2D and 3D characteristics. Among them, the double-ring model represents the 2D distribution of scatterers and the hemispherical model represents the 3D distribution of scatterers. The fading channel matrix from all the APs to the user is described by an Q×P matrix of complex fading envelopes, which can be represented as follows.(1)Ht=h11th12t…h1Pth21th22t…h2Pt…………hQ1thQ2t…hQPtQ×P,where hqpt denotes the channel from the q-th AP p-th antenna to the user. Note that hqpt includes the LoS path component and the NLoS path component, where the NLoS path component consists of the single scattering path component passing through T−ni−Ri=1,2,3 and the double scattering path component passing through T−n1−n2−R. The channel of hqpt can be represented as follows:(2)hqpt=hqpLoSt+∑i=13hqpSit+hqpDt.The channel of LoS path component is written as follows:(3)hqpLoSt=KΥqpK+1e−j2πfcτqp×ej2πtfRmaxcosβqpR−βucosαqpR.The channel of the single scattering component is written as follows:(4)hqpSit=PiΥqpK+1limNi⟶∞∑ni=1Ni1Niejφnie−j2πfcτqpni×ej2πtfRmaxcosβR/qp,ni−βucosαR/qp,ni.The channel of the double scattering component is written as follows:(5)hqpDt=PDΥqpK+1limN1,N2⟶∞∑n1=1N1∑n2=1N21N1N2×ejφn1,n2e−j2πfcτqp,n1,n2×ej2πtfRmaxcosβqp,n2R−βucosαqp,n2R,where K denotes the Rice factor. Υqp=Ωqp+Xqp denotes the large-scale fading that consists of path loss and shadowed fading. Ωqp denotes distance-dependent path loss, which can be expressed as Ωqp=Amin+10nlogdqp/dmin, where dmin denotes the minimum propagation distance, Amin denotes the path loss at the minimum propagation distance, and n denotes the path loss index and depends on the specific propagation environment. Xqp denotes the shadowed fading. The measured results show that the log-normal distribution can be used to describe shadowed fading, and the corresponding probability density function can be expressed as fXx=1/2πσxexp−lnx−lnμ2/2σ2, where μ and σ denote the mean power and variance of shadowed fading, respectively. fRmax denotes the maximum Doppler frequency caused by the user, and Pi and PD denote the normalized power coefficient that satisfies ∑i=13Pi+PD=1, specifying how much the single- and double-bounced rays contribute to the total scattered power. The phases φni and φn1,n2 are independent and identically distributed (i.i.d.) random variables that satisfy uniform distributions over −π,π. Note that Sii=1,2,3 represents the type of single scattering, i.e., the single scattering path T−ni−R from the transmitter AP through the scatterer nini=1,2,…,Ni to the user. For the double scattering path, only the path through the two single-ring scatterers to the user is considered, i.e., the T−n1−n2−R path. Different paths have different delays to users. To accurately characterize the channel model, it is necessary to calculate the delays of each path. The time delay of the LoS path, the single scattering path, and the double scattering path can be represented as τqp=Dqp/c, τqpni=Dqp,niT+Dqp,niR/c, and τqp,n1,n2=Dqp,n1T+Dqp,n1n2+Dqp,n2R/c, respectively. According to the geometric relationship, the abovementioned distance parameters can be obtained as Dqp=Dq2+δp−12, Dqp,n1T=Rq−p−1δsinβqp,n1T, and Dqp,n2R=Rq.(6)Dqp,n1R=Rq−p−1δsinβqp,n1Tsinβqp,n1T+p−1δ2+Dq−Rq−p−1δsinβqp,n1Tcosβqp,n1T2.(7)Dqp,n1n2=Rsinβqp,n2R−p−1δ−Dq−p−1δsinβqp,n1Tsinβqp,n1T2+Dq−Dq−p−1δsinβqp,n1Tcosβqp,n1T+Rcosβqp,n2R2.The parameter definitions that appear above are shown in Table1, where Dqp,n2T=Dq+Rcosβqp,n2R2+Rcosβqp,n2R−p−1δ2. The expressions of distance parameters Dqp,n1R and Dqp,n1n2 are given by (6) and (7), respectively.Proof.
See Appendix A.
For a hemisphere, ifDqp,n3TR is expressed directly in terms of the departure of the AP and the arrival angle of the user, the formula is very complicated. Therefore, the position of the scattering cluster, i.e., the intersection of the departure angle and the arrival angle, is used to express Dqp,n3TR. We take the center of the sphere as the coordinate origin to construct a Cartesian coordinate system. The 3D coordinates of the scattering clusters distributed in the hemisphere can be expressed as Dq/2cosαqp,n3cosβqp,n3cosαqp,n3sinβqp,n3sinαqp,n3T, the coordinates of the q-th AP p-th antenna can be expressed as −Dq/2,−δp−1,0, and the coordinates of the user can be expressed as Dq/2,0,0. According to the above coordinates, the expressions of Dqp,n3T and Dqp,n3R are given by (8) and (9), respectively.Proof.
See Appendix B.(8)Dqp,n3T=Dq2cosαqp,n3cosβqp,n3+Dq22+Dq2cosαqp,n3sinβqp,n3+δp−12+Dq2sinαqp,n32.(9)Dqp,n3R=Dq2cosαqp,n3sinβqp,n32+Dq2sinαqp,n32+Dq2cosαqp,n3cosβqp,n3−Dq22.
Note that the azimuth angle of departure (AAoD) of each single scattering path corresponds to the azimuth angle of arrival (AAoA) with respect to the geometric relationship, and the elevation angle of departure (EAoD) and the elevation angle of arrival (EAoA) have a similar relationship. Due to this correspondence, the unknown variable groupsαqp,niT,βqp,niT and αqp,niR,βqp,niR are reduced by two, i.e., we need to know only αqp,niT,βqp,niT or αqp,niR,βqp,niR to know the other parameters. According to the geometric relationship, the available angle parameters are as βqp,n1T=arctanDq−Rsinβqp,n1R/Rcosβqp,n1R, βqp,n2R=arctanDq−Rsinβqp,n2T/Rcosβqp,n2T, αqp,n3T=arcsinDqsinαqp,n3/2/Dqp,n3T, βqp,n3T=arctanDqcosαqp,n3sinβqp,n3+2δp−1/Dqcosαqp,n3cosβqp,n3+Dq, αqp,n3R=arcsinDqsinαqp,n3/2/Dqp,n3R, and βqp,n3R=arctancosαqp,n3sinβqp,n3/cosαqp,n3cosβqp,n3−1.
The double scattering path can be randomly generated and the abovementioned correspondence does not exist. At this point, the channel model between the AP and user is completed. The specific process to establish a GBSM is given in Algorithm1. Note that the GBSM proposed for cell-free massive MIMO channels can be degenerated into a centralized massive MIMO channel model. When there exists only 1 AP, i.e., Q=1, the channel matrix in (1) can be re-expressed as H′t=h11t,h12t,…,h1Pt1×P, which denotes the channel matrix of the transmitter with P antennas and h1pt denotes the channel from each antenna to the user. When P⟶∞, H′t is equivalent to the SBEM of a centralized massive MIMO as shown in [9], h=1/L∑l=1Lβlαlaϕl, where αl denotes the complex gain of the lth path that represents the small-scale Rayleigh fading, βl denotes the large-scale fading, and aϕl denotes the array steering vector and aϕl=1/N1,ej2πδ/λsinϕl,…,ejN−12πδ/λsinϕl. Although the degeneration model H′t can characterize centralized massive MIMO channels, there are still some differences compared with the multipath model h. The multipath model assumes parallel waves and the phase differences of the subantennas are all 2πδ/λsinϕl. However, each subantenna of the degeneration model has an independent and accurate expression by substituting into a specific position. It takes into account the phase difference caused by the distance and the Doppler frequency shift, which is a more accurate channel model.
Algorithm 1: The specific process to establish a GBSM.
Step 1: Take theq-th AP and user as the centers of circles to establish a double-ring model. Take the center of the connection line as the center of the sphere to establish a hemisphere model and Cartesian coordinate system.Step 2: Randomly generate the distribution of scatterers, generate AAoD/EAoD/AAoA/EAoA for each single scattering path and double scattering path, and calculate the distance of each path.Step 3: Obtain the time delaysτqp, τqpni, and τqp,n1n2 for each LoS path, single scattering path, and double scattering path.Step 4: Allocate the normalized power coefficientsPi and PD for the single scattering path and double scattering path, respectively.Step 5: Combine Step 2–Step 4 to obtain the normalized channel expression, i.e.,hqpt, between the q-th AP p-th antenna and user.Step 6: Repeat Step 1–Step 5 to build the channel matrixHt.


### 2.2. Conversion of a Discrete Angle into a Continuous Angle

Since it is assumed that the number of scatterers is infinite, i.e.,Ni⟶∞i=1,2,3, the model proposed in (2) is actually a mathematical reference model. Due to infinite complexity, these conditions cannot be achieved in practice. However, as described in Reference [40], the reference model can be used for theoretical analysis, the design of communication systems, and the realizable simulation models with reasonable complexity. For our reference model, we assume that Nii=1,2,3 is a finite value and convert the discrete angle value into a continuous angle value with a certain probability density distribution. Many different distributions have been proposed to characterize the angular distribution, such as the uniform, Gaussian, wrapped Gaussian, and heart-shaped probability density function (PDF) distributions [41]. Here, we use the von Mises–Fisher (vMF) PDF distribution [42], which can effectively simulate the distribution of scatterers in 3D space, and the expression is as follows.(10)fα,β=kcosα4πsinhkekcosα0cosαcosβ−β0+sinα0sinα,where sinh⋅ denotes the hyperbolic sine function and α and β denote the random variables of the elevation angle and azimuth angle, respectively. α0 and β0 denote the average values of the elevation angle and azimuth angle, respectively. When the scatterers are distributed on a double ring, the elevation angle tends to be ignored. At this time, only the distribution of the azimuth angle is considered. The vMF distribution should simplify to the von Mises distribution by integrating over the angle α, and the probability density function is fβ=ekcosβ−β0/2πI0k. Note that in the two single scattering paths, due to the low height of the surrounding scatterers, the influence of the elevation angle is ignored, i.e., αqp,niT≈0, αqp,niR≈0. The calculation of the available angle parameters indicates that due to the corresponding relationship between the azimuth angles βqp,niT and βqp,niR, the 2D variable βqp,niT,βqp,niR|i=1,2 can be converted into a probability density function that needs to consider only the one-dimensional variable βqp,niR. In the double scattering path T−n1−n2−R, the azimuth angles of the AP and user are independent of each other. Therefore, the probability density function products of βqp,n1T and βqp,n1R should be considered. In the single scattering path, since the scatterers are distributed on the hemisphere, the joint probability density distribution function of the elevation angle and the azimuth angle αqp,n3T,βqp,n3R should be considered.


## 2.1. GBSM

We consider a cell-free massive MIMO system withQ APs. Suppose that Q APs are randomly deployed in a certain city area, all of which are equipped with a uniform linear array (ULA). Each AP is equipped with P antennas, and the distance between each antenna element is δ. The far-field plane wavefront assumption in conventional massive MIMO systems is not fulfilled in cell-free massive MIMO systems, as the condition that the distance from the transmitting antenna to the scatterer cluster being far smaller than the Rayleigh distance cannot be upheld. As a result, the wavefront emitted from the n-th cluster to the receive array is assumed to be spherical. Then, the angle of arrival (AoA) is no longer linear along the array, and it needs to be computed based on geometrical relationships. Considering the channel between any AP and the user, the signal sent from the AP transmitter is scattered through the line-of-sight (LoS) path and the non-line-of-sight (NLoS) path to the user. The types of scatterers in the NLoS path can be divided into three categories: scatterers around APs, scatterers around users, and scatterers such as tall buildings and trees. We fully consider the geographic location of the APs and users and the geometric relationship of the scatterer distribution and construct a GBSM with the spherical wave assumption for cell-free massive MIMO channels, as shown in Figure 2. Some of the main parameters and definitions are shown in Table 1. The GBSM needs to be extended for different scenarios to effectively improve the accuracy of channel modeling, that is, the combination of the double-ring and hemisphere models. The reason for adopting the double-ring and hemisphere models is as follows. This paper considers a cell-free massive MIMO system model in an urban rich scattering scenario. In order to distinguish between near-end scattering and far-end scattering, this paper adopts a GBSM combined with a double-ring model and a hemisphere model. When the user faces a single AP with a short distance, there are fewer channels reflected by the AP through the higher scatterers, which can be ignored, and the channel with the AP is modeled as a single ring. This is similar to the user and the single ring model is adopted. For the channels between the user and other APs that are far away, the reflection effect of the higher scatterers needs to be considered to form a hemispherical model. Since APs are randomly distributed in the system, some APs are close to the user and the direct path is strong at this time. The scattering of high scatterers can be ignored and only the single ring model is considered. Some APs far away from the user have a weak direct path and can be ignored, while the scattering of high scatterers cannot be ignored, so only the hemispherical model is considered.Figure 2
GBSM for cell-free massive MIMO channels.Table 1
Definition of parameters.
ParametersDefinitionλCarrier wavelengthfcCarrier frequencyδAntenna spacingvuSpeed of the userβuAzimuth angle of the user’s movement directionKRician factorfDDoppler frequencyfRmaxMaximum Doppler frequency with respect to the userRq, RDistribution radii of scatterers around the AP and user, respectivelyPii=1,2,3Normalized power coefficient of the single scattering pathPDNormalized power coefficient of the double scattering pathφni, φn1,n2Initial random phase of the LoS for single and double scattering paths, respectivelyDqpDistance between theq-th AP p-th antenna and userDqp,niTi=1,2,3Distance between theq-th AP p-th antenna and the ni cluster of scatterers through the single scattering pathDqp,niRi=1,2,3Distance between theni cluster of scatterers and the user through the single scattering pathDqp,n1n2Distance between then1 and the n2 clusters of scatterers through the double scattering pathαqpT, βqpTAAoD/EAoD between theq-th AP p-th antenna and the user through the LoS pathαqpR, βqpRAAoA/EAoA between theq-th AP p-th antenna and the user through the LoS pathαqp,niT, βqp,niTAAoD/EAoD between theq-th AP p-th antenna and the ni cluster of scatterers through the single scattering pathαqp,niR, βqp,niRAAoA/EAoA between theni cluster of scatterers and user through the single scattering pathαqp,n3, βqp,n3Elevation angle and azimuth angle between the hemisphere center and then3 cluster of scatterersτqpTime delay between theq-th AP p-th antenna and user through the LoS pathτqpniTime delay between theq-th AP p-th antenna and the ni cluster of scatterers through the single scattering pathτqp,n1,n2Time delay between then1 and the n2 cluster of scatterers through the double scattering pathFor scatterers around APs and users where the elevation angle/depression angle changes are relatively small, we assume that both categories are distributed on a single ring. Since the scatterers around tall buildings and trees are relatively concentrated, the 3D distribution characteristics of APs are particularly prominent. Hemisphere models need to be used to model cell-free massive MIMO channels. We take the first antenna element of the AP transmitter and the user as the centers to establish a double-ring model with radii ofRq and R and establish a hemisphere model with the two centers as diameter Dq. We set the layout direction of the ULA to the Y axis. Suppose that N1 scatterers are distributed on a circle centered on the AP and that N2 scatterers are distributed on a circle centered on the user. The moving speed of the user is vu, and the moving azimuth angle is βu. Suppose that N3 stationary scatterers are distributed in the hemisphere area. Note that the GBSM proposed in this paper includes 2D and 3D characteristics. Among them, the double-ring model represents the 2D distribution of scatterers and the hemispherical model represents the 3D distribution of scatterers. The fading channel matrix from all the APs to the user is described by an Q×P matrix of complex fading envelopes, which can be represented as follows.(1)Ht=h11th12t…h1Pth21th22t…h2Pt…………hQ1thQ2t…hQPtQ×P,where hqpt denotes the channel from the q-th AP p-th antenna to the user. Note that hqpt includes the LoS path component and the NLoS path component, where the NLoS path component consists of the single scattering path component passing through T−ni−Ri=1,2,3 and the double scattering path component passing through T−n1−n2−R. The channel of hqpt can be represented as follows:(2)hqpt=hqpLoSt+∑i=13hqpSit+hqpDt.The channel of LoS path component is written as follows:(3)hqpLoSt=KΥqpK+1e−j2πfcτqp×ej2πtfRmaxcosβqpR−βucosαqpR.The channel of the single scattering component is written as follows:(4)hqpSit=PiΥqpK+1limNi⟶∞∑ni=1Ni1Niejφnie−j2πfcτqpni×ej2πtfRmaxcosβR/qp,ni−βucosαR/qp,ni.The channel of the double scattering component is written as follows:(5)hqpDt=PDΥqpK+1limN1,N2⟶∞∑n1=1N1∑n2=1N21N1N2×ejφn1,n2e−j2πfcτqp,n1,n2×ej2πtfRmaxcosβqp,n2R−βucosαqp,n2R,where K denotes the Rice factor. Υqp=Ωqp+Xqp denotes the large-scale fading that consists of path loss and shadowed fading. Ωqp denotes distance-dependent path loss, which can be expressed as Ωqp=Amin+10nlogdqp/dmin, where dmin denotes the minimum propagation distance, Amin denotes the path loss at the minimum propagation distance, and n denotes the path loss index and depends on the specific propagation environment. Xqp denotes the shadowed fading. The measured results show that the log-normal distribution can be used to describe shadowed fading, and the corresponding probability density function can be expressed as fXx=1/2πσxexp−lnx−lnμ2/2σ2, where μ and σ denote the mean power and variance of shadowed fading, respectively. fRmax denotes the maximum Doppler frequency caused by the user, and Pi and PD denote the normalized power coefficient that satisfies ∑i=13Pi+PD=1, specifying how much the single- and double-bounced rays contribute to the total scattered power. The phases φni and φn1,n2 are independent and identically distributed (i.i.d.) random variables that satisfy uniform distributions over −π,π. Note that Sii=1,2,3 represents the type of single scattering, i.e., the single scattering path T−ni−R from the transmitter AP through the scatterer nini=1,2,…,Ni to the user. For the double scattering path, only the path through the two single-ring scatterers to the user is considered, i.e., the T−n1−n2−R path. Different paths have different delays to users. To accurately characterize the channel model, it is necessary to calculate the delays of each path. The time delay of the LoS path, the single scattering path, and the double scattering path can be represented as τqp=Dqp/c, τqpni=Dqp,niT+Dqp,niR/c, and τqp,n1,n2=Dqp,n1T+Dqp,n1n2+Dqp,n2R/c, respectively. According to the geometric relationship, the abovementioned distance parameters can be obtained as Dqp=Dq2+δp−12, Dqp,n1T=Rq−p−1δsinβqp,n1T, and Dqp,n2R=Rq.(6)Dqp,n1R=Rq−p−1δsinβqp,n1Tsinβqp,n1T+p−1δ2+Dq−Rq−p−1δsinβqp,n1Tcosβqp,n1T2.(7)Dqp,n1n2=Rsinβqp,n2R−p−1δ−Dq−p−1δsinβqp,n1Tsinβqp,n1T2+Dq−Dq−p−1δsinβqp,n1Tcosβqp,n1T+Rcosβqp,n2R2.The parameter definitions that appear above are shown in Table1, where Dqp,n2T=Dq+Rcosβqp,n2R2+Rcosβqp,n2R−p−1δ2. The expressions of distance parameters Dqp,n1R and Dqp,n1n2 are given by (6) and (7), respectively.Proof.
See Appendix A.
For a hemisphere, ifDqp,n3TR is expressed directly in terms of the departure of the AP and the arrival angle of the user, the formula is very complicated. Therefore, the position of the scattering cluster, i.e., the intersection of the departure angle and the arrival angle, is used to express Dqp,n3TR. We take the center of the sphere as the coordinate origin to construct a Cartesian coordinate system. The 3D coordinates of the scattering clusters distributed in the hemisphere can be expressed as Dq/2cosαqp,n3cosβqp,n3cosαqp,n3sinβqp,n3sinαqp,n3T, the coordinates of the q-th AP p-th antenna can be expressed as −Dq/2,−δp−1,0, and the coordinates of the user can be expressed as Dq/2,0,0. According to the above coordinates, the expressions of Dqp,n3T and Dqp,n3R are given by (8) and (9), respectively.Proof.
See Appendix B.(8)Dqp,n3T=Dq2cosαqp,n3cosβqp,n3+Dq22+Dq2cosαqp,n3sinβqp,n3+δp−12+Dq2sinαqp,n32.(9)Dqp,n3R=Dq2cosαqp,n3sinβqp,n32+Dq2sinαqp,n32+Dq2cosαqp,n3cosβqp,n3−Dq22.
Note that the azimuth angle of departure (AAoD) of each single scattering path corresponds to the azimuth angle of arrival (AAoA) with respect to the geometric relationship, and the elevation angle of departure (EAoD) and the elevation angle of arrival (EAoA) have a similar relationship. Due to this correspondence, the unknown variable groupsαqp,niT,βqp,niT and αqp,niR,βqp,niR are reduced by two, i.e., we need to know only αqp,niT,βqp,niT or αqp,niR,βqp,niR to know the other parameters. According to the geometric relationship, the available angle parameters are as βqp,n1T=arctanDq−Rsinβqp,n1R/Rcosβqp,n1R, βqp,n2R=arctanDq−Rsinβqp,n2T/Rcosβqp,n2T, αqp,n3T=arcsinDqsinαqp,n3/2/Dqp,n3T, βqp,n3T=arctanDqcosαqp,n3sinβqp,n3+2δp−1/Dqcosαqp,n3cosβqp,n3+Dq, αqp,n3R=arcsinDqsinαqp,n3/2/Dqp,n3R, and βqp,n3R=arctancosαqp,n3sinβqp,n3/cosαqp,n3cosβqp,n3−1.
The double scattering path can be randomly generated and the abovementioned correspondence does not exist. At this point, the channel model between the AP and user is completed. The specific process to establish a GBSM is given in Algorithm1. Note that the GBSM proposed for cell-free massive MIMO channels can be degenerated into a centralized massive MIMO channel model. When there exists only 1 AP, i.e., Q=1, the channel matrix in (1) can be re-expressed as H′t=h11t,h12t,…,h1Pt1×P, which denotes the channel matrix of the transmitter with P antennas and h1pt denotes the channel from each antenna to the user. When P⟶∞, H′t is equivalent to the SBEM of a centralized massive MIMO as shown in [9], h=1/L∑l=1Lβlαlaϕl, where αl denotes the complex gain of the lth path that represents the small-scale Rayleigh fading, βl denotes the large-scale fading, and aϕl denotes the array steering vector and aϕl=1/N1,ej2πδ/λsinϕl,…,ejN−12πδ/λsinϕl. Although the degeneration model H′t can characterize centralized massive MIMO channels, there are still some differences compared with the multipath model h. The multipath model assumes parallel waves and the phase differences of the subantennas are all 2πδ/λsinϕl. However, each subantenna of the degeneration model has an independent and accurate expression by substituting into a specific position. It takes into account the phase difference caused by the distance and the Doppler frequency shift, which is a more accurate channel model.
Algorithm 1: The specific process to establish a GBSM.
Step 1: Take theq-th AP and user as the centers of circles to establish a double-ring model. Take the center of the connection line as the center of the sphere to establish a hemisphere model and Cartesian coordinate system.Step 2: Randomly generate the distribution of scatterers, generate AAoD/EAoD/AAoA/EAoA for each single scattering path and double scattering path, and calculate the distance of each path.Step 3: Obtain the time delaysτqp, τqpni, and τqp,n1n2 for each LoS path, single scattering path, and double scattering path.Step 4: Allocate the normalized power coefficientsPi and PD for the single scattering path and double scattering path, respectively.Step 5: Combine Step 2–Step 4 to obtain the normalized channel expression, i.e.,hqpt, between the q-th AP p-th antenna and user.Step 6: Repeat Step 1–Step 5 to build the channel matrixHt.


## 2.2. Conversion of a Discrete Angle into a Continuous Angle

Since it is assumed that the number of scatterers is infinite, i.e.,Ni⟶∞i=1,2,3, the model proposed in (2) is actually a mathematical reference model. Due to infinite complexity, these conditions cannot be achieved in practice. However, as described in Reference [40], the reference model can be used for theoretical analysis, the design of communication systems, and the realizable simulation models with reasonable complexity. For our reference model, we assume that Nii=1,2,3 is a finite value and convert the discrete angle value into a continuous angle value with a certain probability density distribution. Many different distributions have been proposed to characterize the angular distribution, such as the uniform, Gaussian, wrapped Gaussian, and heart-shaped probability density function (PDF) distributions [41]. Here, we use the von Mises–Fisher (vMF) PDF distribution [42], which can effectively simulate the distribution of scatterers in 3D space, and the expression is as follows.(10)fα,β=kcosα4πsinhkekcosα0cosαcosβ−β0+sinα0sinα,where sinh⋅ denotes the hyperbolic sine function and α and β denote the random variables of the elevation angle and azimuth angle, respectively. α0 and β0 denote the average values of the elevation angle and azimuth angle, respectively. When the scatterers are distributed on a double ring, the elevation angle tends to be ignored. At this time, only the distribution of the azimuth angle is considered. The vMF distribution should simplify to the von Mises distribution by integrating over the angle α, and the probability density function is fβ=ekcosβ−β0/2πI0k. Note that in the two single scattering paths, due to the low height of the surrounding scatterers, the influence of the elevation angle is ignored, i.e., αqp,niT≈0, αqp,niR≈0. The calculation of the available angle parameters indicates that due to the corresponding relationship between the azimuth angles βqp,niT and βqp,niR, the 2D variable βqp,niT,βqp,niR|i=1,2 can be converted into a probability density function that needs to consider only the one-dimensional variable βqp,niR. In the double scattering path T−n1−n2−R, the azimuth angles of the AP and user are independent of each other. Therefore, the probability density function products of βqp,n1T and βqp,n1R should be considered. In the single scattering path, since the scatterers are distributed on the hemisphere, the joint probability density distribution function of the elevation angle and the azimuth angle αqp,n3T,βqp,n3R should be considered.


## 3. Statistical Characteristic Analysis

Since the small-scale fading of the wireless channel causes the received signal to fluctuate randomly, it is necessary to use statistical characteristics to describe the fading characteristics of the wireless channel. Statistical characteristics are used to better understand the nonstationary channel behavior. Statistical feature parameters mainly include first-order statistical properties probability density function (PDF), and second-order statistical properties, including space CCF, time ACF, DPSD, LCR, and AFD. The first-order statistical parameters cannot fully describe the characteristics of the wireless channel, so it is necessary to use the second-order statistical parameters to supplement the characteristics of the wireless channel.


### 3.1. Space-Time Correlation Function (ST CF)

The ST CF can describe the variation in the channel in time and space and indicate the characteristics of the channel [43]. Therefore, to further verify the performance of the proposed GBSM, it is necessary to investigate the ST CF. The ST CF between any two channels hqpt and hq′p′t can be represented as follows:(11)Hq′p′,qpΔr,τ=Ehqpt×hq′p′∗t−τ,where ⋅∗ denotes the complex conjugate operation, E⋅ is the statistical expectation operator, p,p′∈1,2,…,P, and q,q′∈1,2,…,Q. Note that (11) is a function of space separation Δr caused by different APs and antennas and time separation τ. With the assumption that the LoS path and NLoS path components are independent [43], we can rewrite formula (11) as the sum of the correlation of the LoS path, single scattering path, and double scattering path components, which can be given by the following equation:(12)Hq′p′,qpΔr,τ=HhqpLoShq′p′LoSΔr,τ+∑i=13HhqpSihq′p′SiΔr,τ+HhqpDhq′p′DΔr,τ.For the LoS path component, the azimuth angleβqp,RLoS≈π and the elevation angle αqp,RLoS≈0 can be considered. Incorporating formula (12), the ST CF of the LOS path component is represented as follows:(13)HhqpLoShq′p′LoSΔr,τ=KΥqpK+1ej2πfcτq′p′−τqpe−j2πτfRmaxcosβu.Here,τqp denotes the time delay and τqp=Dq2+δp−12/c, τq′p′=Dq′2+δp′−12/c. The ST CF of a single scattering component through n1 and n2 scatterers, i.e., the paths T−n1−R and T−n2−R are represented as follows:(14)HhqpSihq′p′SiΔr,τ=PiΥqpK+1∫−ππe−j2πτfRmaxcosβR/qp,ni−βu×ej2πfcτq′p′ni−τqpnifβqp,niRdβqp,niR,where fβqp,niR=ekcosβqp,niR−βqpR/2πI0k and nii=1,2. The ST CF of the single scattering component through n3 scatterers, i.e., the path T−n3−R, is represented as (15):(15)HhqpS3hq′p′S3Δr,τ=P3ΥqpK+1∫0π/2∫−ππej2πfcτq′p′n3−τqpn3×e−j2πτfRmaxcosβqp,n3R−βucosαqp,n3Rfαqp,n3R,βqp,n3Rdαqp,n3Rdβqp,n3R.where(16)fαqp,n3R,βqp,n3R=kcosαqp,n3R4πsinhk×ekcosαqp,n3cosαqp,n3Rcosβqp,n3R−βqp,n3+sinαqp,n3sinαqp,n3R.The ST CF of the double scattering component throughn1,n2 scatterers, i.e., the path T−n1−n2−R, is represented as follows:(17)HhqpDhq′p′DΔr,τ=PDΥqpK+1∫−ππ∫−ππej2πfcτq′p′,n1,n2−τqp,n1,n2×e−j2πτfRmaxcosβqp,n2R−βufβqp,n1Tfβqp,n2Rdβqp,n1Tdβqp,n2R.wherefβqp,n1T=ekcosβqp,n1T−βqpT/2πI0k and fβqp,n2R=ekcosβqp,n2R−βqpR/2πI0k.


### 3.2. Space CCF

For formula (12), we set the delay factor τ=0, and we can obtain the space CCF.(18)Hq′p′,qpΔr=HhqpLoShq′p′LoSΔr+∑i=13HhqpSihq′p′SiΔr+HhqpDhq′p′DΔr.The space CCF of the LoS path can be obtained as follows:(19)HhqpLoShq′p′LoSΔr=KΥqpK+1ej2πfc/cDq′2+δp′−12/−Dq2+δp−12.The space CCF of a single scattering component throughn1 and n2 scatterers, i.e., the paths T−n1−R and T−n2−R are represented as follows:(20)HhqpSihq′p′SiΔr=PiΥqpK+1∫−ππekcosβqp,niR−β02πI0k×ej2πfcτq′p′ni−τqpnidβqp,niR.The space CCF of a single scattering component throughn3 scatterers, i.e., the path T−n3−R is represented as (21), and the double scattering component through n1,n2 scatterers is represented as formula (22).(21)HhqpS3hq′p′S3Δr=P3ΥqpK+1∫0π2∫−ππkcosαqp,n3R4πsinhkej2πfcτq′p′n3−τqpn3×ekcosαqp,n3cosαqp,n3Rcosβqp,n3R−βqp,n3+sinαqp,n3sinαqp,n3Rdαqp,n3Rdβqp,n3R.(22)HhqpDhq′p′DΔr=PDΥqpK+1∫−ππ∫−ππ12πI0k2ekcosβqp,n1T−βqpT×ekcosβqp,n2R−βqpRej2πfcτq′p′,n1,n2−τqp,n1,n2dβqp,n1Tdβqp,n2R.


### 3.3. Time ACF

By setting the indexesq′=q and p′=p in formula (12), the time ACF is represented as follows:(23)Hqp,qpτ=HhqpLoShqpLoSτ+∑i=13HhqpSihqpSiτ+HhqpDhqpDτ.The time ACF of the LoS path is represented as follows:(24)HhqpLoShqpLoSτ=KΥqpK+1ej2πfc/cDq2+δp−12−2πτfRmaxcosβu.The time ACF of a single scattering component throughn1 and n2 scatterers, i.e., the paths T−n1−R and T−n2−R, is represented as formula (25). The single scattering component through n3 scatterers, i.e., the path T−n3−R is represented as formula (26), and the double scattering component through T−n1−n2−R path is represented as formula (27).(25)HhqpSihqpSiτ=PiΥqpK+1∫−ππekcosβqp,niR−β02πI0ke−j2πτfRmaxcosβqp,niR−βudβqp,niR.(26)HhqpS3hqpS3τ=P3ΥqpK+1∫0π/2∫−ππkcosαqp,n3R4πsinhkekcosαqp,n3cosαqp,n3Rcosβqp,n3R−βqp,n3+sinαqp,n3sinαqp,n3R×ej2πfcτqpn3e−j2πτfRmaxcosβqp,n3R−βucosαqp,n3Rdαqp,n3Rdβqp,n3R.(27)HhqpDhqpDτ=PDΥqpK+1∫−ππ∫−ππ12πI0k2ekcosβqp,n1T−βqpTekcosβqp,n2R−βqpR×ej2πfcτqp,n1,n2e−j2πτfRmaxcosβqp,n2R−βudβqp,n1Tdβqp,n2R.


### 3.4. Doppler PSD Function

Applying the Fourier transform to the space-time function in (12) in terms of τ, we can obtain the corresponding Doppler PSD function as follows:(28)Sq′p′,qpfD=∫−∞∞Hq′p′,qpΔr,τe−j2πfDτdτ=ShqpLoShq′p′LoSfD+∑i=13ShqpSihq′p′SifD+ShqpDhq′p′DfD,where fD denotes the Doppler frequency, which satisfies fD≤fRmax.


### 3.5. Envelope LCR and AFD

The LCR and AFD are second-order statistical parameters that describe how fast the fading of a wireless channel changes. The LCR refers to the rate at which the amplitude of a fading signal crosses a specified levelr with a positive or negative slope in a unit time. Using the traditional PDF-based method [44], we derive the expression of the LCR for cell-free massive MIMO channels of the LoS path as follows:(29)Lr=2rK+1π3/2b2b0−b12b02×e−K−K+1r2×∫0π2cosh2KK+1⋅rcosθ×e−χsinθ2+πχsinθ⋅erfχsinθdθ,where cosh⋅ is the hyperbolic cosine function, erf⋅ is the error function,χ=Kb12/b0b2−b12, b0=1/2K+1, bm=∑i=13bmSi+bmD,m∈1,2, bmSi=Pi/2K+12πm×1/Ni∑ni=1NifRmaxcosβqp,niR−βucosαqp,niRm, and bmD=PD/K+12πm×1/N2∑n2=1N2fRmaxcosβqp,n2R−βucosαqp,n2Rm.The LCR for NLoS conditions can be obtained from (8) by setting K=0. The AFD is defined as the mathematical expectation of a duration lower than the specified level in a unit time [45], and the AFD of the LoS path can be written as follows:(30)Ar=1−Q2K,2K+1r2Lr,where Qa,b denotes the Marcum Q function. Qa,b=∫b∞x⋅e−x2+a2/2I0axdx, where I0⋅ is the zeroth-order modified Bessel function of the first kind. We can obtain the closed-form solution of the Marcum Q function as follows:(31)Qa,b=e−a2+b2/2I0ab+a2π∫0πcosθ⋅e−asinθ2/2erfcb−acosθ2dθ.Proof.
See Appendix C
For NLoS conditions, the AFD in (30) simplifies to Ar=1−e−r2/Lr.


## 3.1. Space-Time Correlation Function (ST CF)

The ST CF can describe the variation in the channel in time and space and indicate the characteristics of the channel [43]. Therefore, to further verify the performance of the proposed GBSM, it is necessary to investigate the ST CF. The ST CF between any two channels hqpt and hq′p′t can be represented as follows:(11)Hq′p′,qpΔr,τ=Ehqpt×hq′p′∗t−τ,where ⋅∗ denotes the complex conjugate operation, E⋅ is the statistical expectation operator, p,p′∈1,2,…,P, and q,q′∈1,2,…,Q. Note that (11) is a function of space separation Δr caused by different APs and antennas and time separation τ. With the assumption that the LoS path and NLoS path components are independent [43], we can rewrite formula (11) as the sum of the correlation of the LoS path, single scattering path, and double scattering path components, which can be given by the following equation:(12)Hq′p′,qpΔr,τ=HhqpLoShq′p′LoSΔr,τ+∑i=13HhqpSihq′p′SiΔr,τ+HhqpDhq′p′DΔr,τ.For the LoS path component, the azimuth angleβqp,RLoS≈π and the elevation angle αqp,RLoS≈0 can be considered. Incorporating formula (12), the ST CF of the LOS path component is represented as follows:(13)HhqpLoShq′p′LoSΔr,τ=KΥqpK+1ej2πfcτq′p′−τqpe−j2πτfRmaxcosβu.Here,τqp denotes the time delay and τqp=Dq2+δp−12/c, τq′p′=Dq′2+δp′−12/c. The ST CF of a single scattering component through n1 and n2 scatterers, i.e., the paths T−n1−R and T−n2−R are represented as follows:(14)HhqpSihq′p′SiΔr,τ=PiΥqpK+1∫−ππe−j2πτfRmaxcosβR/qp,ni−βu×ej2πfcτq′p′ni−τqpnifβqp,niRdβqp,niR,where fβqp,niR=ekcosβqp,niR−βqpR/2πI0k and nii=1,2. The ST CF of the single scattering component through n3 scatterers, i.e., the path T−n3−R, is represented as (15):(15)HhqpS3hq′p′S3Δr,τ=P3ΥqpK+1∫0π/2∫−ππej2πfcτq′p′n3−τqpn3×e−j2πτfRmaxcosβqp,n3R−βucosαqp,n3Rfαqp,n3R,βqp,n3Rdαqp,n3Rdβqp,n3R.where(16)fαqp,n3R,βqp,n3R=kcosαqp,n3R4πsinhk×ekcosαqp,n3cosαqp,n3Rcosβqp,n3R−βqp,n3+sinαqp,n3sinαqp,n3R.The ST CF of the double scattering component throughn1,n2 scatterers, i.e., the path T−n1−n2−R, is represented as follows:(17)HhqpDhq′p′DΔr,τ=PDΥqpK+1∫−ππ∫−ππej2πfcτq′p′,n1,n2−τqp,n1,n2×e−j2πτfRmaxcosβqp,n2R−βufβqp,n1Tfβqp,n2Rdβqp,n1Tdβqp,n2R.wherefβqp,n1T=ekcosβqp,n1T−βqpT/2πI0k and fβqp,n2R=ekcosβqp,n2R−βqpR/2πI0k.


## 3.2. Space CCF

For formula (12), we set the delay factor τ=0, and we can obtain the space CCF.(18)Hq′p′,qpΔr=HhqpLoShq′p′LoSΔr+∑i=13HhqpSihq′p′SiΔr+HhqpDhq′p′DΔr.The space CCF of the LoS path can be obtained as follows:(19)HhqpLoShq′p′LoSΔr=KΥqpK+1ej2πfc/cDq′2+δp′−12/−Dq2+δp−12.The space CCF of a single scattering component throughn1 and n2 scatterers, i.e., the paths T−n1−R and T−n2−R are represented as follows:(20)HhqpSihq′p′SiΔr=PiΥqpK+1∫−ππekcosβqp,niR−β02πI0k×ej2πfcτq′p′ni−τqpnidβqp,niR.The space CCF of a single scattering component throughn3 scatterers, i.e., the path T−n3−R is represented as (21), and the double scattering component through n1,n2 scatterers is represented as formula (22).(21)HhqpS3hq′p′S3Δr=P3ΥqpK+1∫0π2∫−ππkcosαqp,n3R4πsinhkej2πfcτq′p′n3−τqpn3×ekcosαqp,n3cosαqp,n3Rcosβqp,n3R−βqp,n3+sinαqp,n3sinαqp,n3Rdαqp,n3Rdβqp,n3R.(22)HhqpDhq′p′DΔr=PDΥqpK+1∫−ππ∫−ππ12πI0k2ekcosβqp,n1T−βqpT×ekcosβqp,n2R−βqpRej2πfcτq′p′,n1,n2−τqp,n1,n2dβqp,n1Tdβqp,n2R.


## 3.3. Time ACF

By setting the indexesq′=q and p′=p in formula (12), the time ACF is represented as follows:(23)Hqp,qpτ=HhqpLoShqpLoSτ+∑i=13HhqpSihqpSiτ+HhqpDhqpDτ.The time ACF of the LoS path is represented as follows:(24)HhqpLoShqpLoSτ=KΥqpK+1ej2πfc/cDq2+δp−12−2πτfRmaxcosβu.The time ACF of a single scattering component throughn1 and n2 scatterers, i.e., the paths T−n1−R and T−n2−R, is represented as formula (25). The single scattering component through n3 scatterers, i.e., the path T−n3−R is represented as formula (26), and the double scattering component through T−n1−n2−R path is represented as formula (27).(25)HhqpSihqpSiτ=PiΥqpK+1∫−ππekcosβqp,niR−β02πI0ke−j2πτfRmaxcosβqp,niR−βudβqp,niR.(26)HhqpS3hqpS3τ=P3ΥqpK+1∫0π/2∫−ππkcosαqp,n3R4πsinhkekcosαqp,n3cosαqp,n3Rcosβqp,n3R−βqp,n3+sinαqp,n3sinαqp,n3R×ej2πfcτqpn3e−j2πτfRmaxcosβqp,n3R−βucosαqp,n3Rdαqp,n3Rdβqp,n3R.(27)HhqpDhqpDτ=PDΥqpK+1∫−ππ∫−ππ12πI0k2ekcosβqp,n1T−βqpTekcosβqp,n2R−βqpR×ej2πfcτqp,n1,n2e−j2πτfRmaxcosβqp,n2R−βudβqp,n1Tdβqp,n2R.


## 3.4. Doppler PSD Function

Applying the Fourier transform to the space-time function in (12) in terms of τ, we can obtain the corresponding Doppler PSD function as follows:(28)Sq′p′,qpfD=∫−∞∞Hq′p′,qpΔr,τe−j2πfDτdτ=ShqpLoShq′p′LoSfD+∑i=13ShqpSihq′p′SifD+ShqpDhq′p′DfD,where fD denotes the Doppler frequency, which satisfies fD≤fRmax.


## 3.5. Envelope LCR and AFD

The LCR and AFD are second-order statistical parameters that describe how fast the fading of a wireless channel changes. The LCR refers to the rate at which the amplitude of a fading signal crosses a specified levelr with a positive or negative slope in a unit time. Using the traditional PDF-based method [44], we derive the expression of the LCR for cell-free massive MIMO channels of the LoS path as follows:(29)Lr=2rK+1π3/2b2b0−b12b02×e−K−K+1r2×∫0π2cosh2KK+1⋅rcosθ×e−χsinθ2+πχsinθ⋅erfχsinθdθ,where cosh⋅ is the hyperbolic cosine function, erf⋅ is the error function,χ=Kb12/b0b2−b12, b0=1/2K+1, bm=∑i=13bmSi+bmD,m∈1,2, bmSi=Pi/2K+12πm×1/Ni∑ni=1NifRmaxcosβqp,niR−βucosαqp,niRm, and bmD=PD/K+12πm×1/N2∑n2=1N2fRmaxcosβqp,n2R−βucosαqp,n2Rm.The LCR for NLoS conditions can be obtained from (8) by setting K=0. The AFD is defined as the mathematical expectation of a duration lower than the specified level in a unit time [45], and the AFD of the LoS path can be written as follows:(30)Ar=1−Q2K,2K+1r2Lr,where Qa,b denotes the Marcum Q function. Qa,b=∫b∞x⋅e−x2+a2/2I0axdx, where I0⋅ is the zeroth-order modified Bessel function of the first kind. We can obtain the closed-form solution of the Marcum Q function as follows:(31)Qa,b=e−a2+b2/2I0ab+a2π∫0πcosθ⋅e−asinθ2/2erfcb−acosθ2dθ.Proof.
See Appendix C
For NLoS conditions, the AFD in (30) simplifies to Ar=1−e−r2/Lr.


## 4. Statistical Characteristic Simulation and Discussion

In this section, we first compare the GBSM with the standard block-fading channel model and the SBEM as shown in Figure3. We can conclude from the simulation result that both the standard block-fading channel model and the SBEM could not describe the nonstationarity in the frequency domain since the normalized Doppler PSD does not change with the frequency domain but is always a discrete point. The same result appeared in the space and time domains. However, the proposed GBSM could describe the nonstationarity of cell-free massive MIMO channels. Moreover, based on the statistical characteristic analysis of the channel, the proposed GBSM is simulated by the space CCF, time ACF, Doppler PSD function, envelope LCR, and AFD. In addition, the effects of the antenna space and azimuth angle on the channel are analyzed in detail. The following parameters [46] are chosen for our simulations. Dq=200, Rq=20, R=20, Q=10, P=5, fc=5.9GHz, fRmax=570Hz, vu=6m/s, βqp,n1R=135.4°, βqp,n2T=−35.7°, αqp,n3=50.2°, βqp,n3=70.6°. βqp,n1T, βqp,n2R, αqp,n3T, αqp,n3R, βqp,n3T, and βqp,n3R can be obtained by geometric relationship in Section 2.Figure 3
Comparison between the GBSM, the standard block-fading channel model, and the SBEM.Note that the value ofk is related to the distribution of scatterers. Generally, the smaller the value of k is, the greater the density of scatterers distributed around the user and AP. The value of K is related to the power intensity of the LoS component. It can be considered that the higher the density of the scatterer is, the smaller the value of K because under this condition. The NLoS component accounts for a larger proportion, while the LoS component accounts for a smaller proportion. This paper simulates two scattering environments, namely, low-density and high-density scenarios, and analyzes the proposed model. Nii=1,2,3 represents the number of finite scatterers, which have different values in different density scenarios. The values of each parameter are shown in Table 2.Table 2
Key parameters of different scenarios.
ParametersLow densityHigh densityK3.7860.156k3.61.3P10.3350.176P20.2030.176P30.4110.163PD0.0510.485N11540N21540N31540


### 4.1. Space CCF

Figure4 analyzes the relationship between the space CCF and the index of antennas with different antenna spacings δ and scatterer densities, i.e., low density and high density. It is obvious that the space CCF varies significantly in different scattering environments. In the high-density scatterer scenario, that is, in a rich scattering environment, the space correlation between antennas is relatively small. The reason is that the higher the scatterer density, the more spatial diversity the massive MIMO channel has. In addition, the antenna spacing has an impact on the space CCF. The space CCF is always lower when δ=λ than when δ=1/2λ. Increasing the antenna spacing causes the CCF to weaken.Figure 4
The relationship between the space CCF and the index of antennas with different antenna spacingsδ and scatterer densities.


### 4.2. Time ACF

Figure5 analyzes the variation between the time ACF and τ with different antenna spacings δ and scatterer densities. It is obvious that the time ACF varies significantly in different scattering environments. Combined with Figure 4, we can see that both the space CCF and time ACF are relatively small and that the values remain approximately unchanged in a rich scattering environment. It can be considered to conform to wide-sense stationarity (WSS) in a certain period of time under high-density conditions. Similar to the space CCF, increasing the antenna spacing causes the ACF to weaken. The time ACF is always lower when δ=λ than when δ=1/2λ, which indicates that time ACF can be reduced by increasing the antenna spacing.Figure 5
Variation between the time ACF andτ at different times with different antenna spacings δ and scatterer densities.Since the value of Rice factorK is related to the power intensity of the LoS component, Figure 6 analyzes the impact of K on the statistics of ACFs for the LoS scenario. The simulation results indicate that when Rice factor K is smaller, the LoS path occupies a smaller space and the ACF of the channel is weaker. The simulation results demonstrate that the cell-free massive MIMO system has more favorable propagation in a rich scattering environment.Figure 6
Impact of Rice factorK on the statistics of ACFs.


### 4.3. Doppler PSD

Figure7 analyzes the relationship between the PSD and the Doppler frequency in the LoS path and NLoS path under different scattering scenarios. The simulation results show that in both the LoS and NLoS paths, the Doppler PSD has a peak and tends to be evenly distributed on the sides. In a rich scattering environment, the PSD in the high-density scenario is larger than that in a low-density scenario. The underlying physical reason is that in the high-density scenario, the received power is mainly obtained by the NLoS path, including all directions reflected by the scatterers distributed on the circle and hemisphere. In contrast, the LoS path obtains a higher PSD in a low-density scattering environment since the received power comes mainly from specific directions identified by LoS components.Figure 7
The relationship between the PSD and the Doppler frequency in the LoS and NLoS paths, whereδ=1/2λ.


### 4.4. Envelope LCR and AFD

Figure8 analyzes the relationship between the envelope LCR and the envelope level for different azimuth angles and scatterer densities. Again, the density of the scatterers significantly affects the envelope LCR for cell-free MIMO channels. This indicates that the LCR is higher when the channel is in a rich scattering environment. Additionally, the azimuth angle of the user has little influence on the LCR.Figure 8
Relationship between the envelope LCR and the envelope level with different azimuth anglesβu and scatterer densities, where δ=1/2λ.Figure9 analyzes the relationship between the envelope AFD and the envelope level with different azimuth angles and scatterer densities. Unlike the LCR, the AFD tends to be small when the channel is in a rich scattering environment. Additionally, the azimuth angle of the user has little influence on the AFD.Figure 9
Relationship between the envelope AFD and the envelope level with different azimuth anglesβu and scatterer densities, where δ=1/2λ.


## 4.1. Space CCF

Figure4 analyzes the relationship between the space CCF and the index of antennas with different antenna spacings δ and scatterer densities, i.e., low density and high density. It is obvious that the space CCF varies significantly in different scattering environments. In the high-density scatterer scenario, that is, in a rich scattering environment, the space correlation between antennas is relatively small. The reason is that the higher the scatterer density, the more spatial diversity the massive MIMO channel has. In addition, the antenna spacing has an impact on the space CCF. The space CCF is always lower when δ=λ than when δ=1/2λ. Increasing the antenna spacing causes the CCF to weaken.Figure 4
The relationship between the space CCF and the index of antennas with different antenna spacingsδ and scatterer densities.


## 4.2. Time ACF

Figure5 analyzes the variation between the time ACF and τ with different antenna spacings δ and scatterer densities. It is obvious that the time ACF varies significantly in different scattering environments. Combined with Figure 4, we can see that both the space CCF and time ACF are relatively small and that the values remain approximately unchanged in a rich scattering environment. It can be considered to conform to wide-sense stationarity (WSS) in a certain period of time under high-density conditions. Similar to the space CCF, increasing the antenna spacing causes the ACF to weaken. The time ACF is always lower when δ=λ than when δ=1/2λ, which indicates that time ACF can be reduced by increasing the antenna spacing.Figure 5
Variation between the time ACF andτ at different times with different antenna spacings δ and scatterer densities.Since the value of Rice factorK is related to the power intensity of the LoS component, Figure 6 analyzes the impact of K on the statistics of ACFs for the LoS scenario. The simulation results indicate that when Rice factor K is smaller, the LoS path occupies a smaller space and the ACF of the channel is weaker. The simulation results demonstrate that the cell-free massive MIMO system has more favorable propagation in a rich scattering environment.Figure 6
Impact of Rice factorK on the statistics of ACFs.


## 4.3. Doppler PSD

Figure7 analyzes the relationship between the PSD and the Doppler frequency in the LoS path and NLoS path under different scattering scenarios. The simulation results show that in both the LoS and NLoS paths, the Doppler PSD has a peak and tends to be evenly distributed on the sides. In a rich scattering environment, the PSD in the high-density scenario is larger than that in a low-density scenario. The underlying physical reason is that in the high-density scenario, the received power is mainly obtained by the NLoS path, including all directions reflected by the scatterers distributed on the circle and hemisphere. In contrast, the LoS path obtains a higher PSD in a low-density scattering environment since the received power comes mainly from specific directions identified by LoS components.Figure 7
The relationship between the PSD and the Doppler frequency in the LoS and NLoS paths, whereδ=1/2λ.


## 4.4. Envelope LCR and AFD

Figure8 analyzes the relationship between the envelope LCR and the envelope level for different azimuth angles and scatterer densities. Again, the density of the scatterers significantly affects the envelope LCR for cell-free MIMO channels. This indicates that the LCR is higher when the channel is in a rich scattering environment. Additionally, the azimuth angle of the user has little influence on the LCR.Figure 8
Relationship between the envelope LCR and the envelope level with different azimuth anglesβu and scatterer densities, where δ=1/2λ.Figure9 analyzes the relationship between the envelope AFD and the envelope level with different azimuth angles and scatterer densities. Unlike the LCR, the AFD tends to be small when the channel is in a rich scattering environment. Additionally, the azimuth angle of the user has little influence on the AFD.Figure 9
Relationship between the envelope AFD and the envelope level with different azimuth anglesβu and scatterer densities, where δ=1/2λ.


## 5. Validation of the GBSM

Note that Section4 considered several features of cell-massive MIMO channel based on reasonable assumptions. We would utilize some relevant channel measurement data [47] about cell-free massive MIMO channel to claim that the proposed GBSM model is accurate and efficient in this Section. Figure 10 indicates the practical measurement scenario where the user trajectory consists of regions R4 and R5. The AP is blocked by the tall building, dividing the user trajectory into NLOS regions. Three different AP array configurations are implemented with a significant variation in the array aperture and antenna element spacing. Specifically, the AP array I denotes a single AP with a linear antenna array of 32 elements aligned horizontally and spaced 0.64λ with frequency of f=3.2 GHz. The AP array II denotes two APs with two antenna arrays of 16 elements, and the two APs are distributed with a distance of 44.6 m between them. The AP array III denotes 32 APs with a single antenna element horizontally aligned but spaced with a distance of 16λ from one AP to the next.Figure 10
Top view of the measurement scenario where three different AP array configurations are implemented.To validate the correctness of the proposed GBSM, the measurement-based scenario [47] and the WINNER II channel model [48] are introduced in this section. The specific details of the two comparation scenarios are described as follows.The measurement-based scenario for cell-free massive MIMO channel in reference [47] is performed in an outdoor environment with a 32-element polarized virtual linear array covering both LoS and NLoS scenarios. As the reference described, it is the world-first fully parallel and coherent widely distributed array channel sounding measurements in a high mobility scenario. Here, we utilize the user’s trajectory of R4 in [47], which can be equivalent to the NLoS path in ours paper.The WINNER II channel model [48] follows a geometry-based stochastic channel modeling approach, which is based on statistical distributions extracted from channel measurement. Here, we utilize the urban microcell scenario in the WINNER II model. In the urban microcell scenario, the height of both the antenna at the BS and at the MS is assumed to be well below the tops of surrounding buildings. The propagation environment of the urban microcell scenario is similar to the cell-free channel in ours paper.We evaluate the delay spread and the Doppler spread as defined in reference [47] on each individual AP antenna from a transmitter for the three considered AP array configurations along the user trajectory.


### 5.1. The Delay Spread

The time delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model are compared in Figure11, Figure 12, and Figure 13 with configurations AP array I, AP array II, and AP array III, respectively. It is clear that the delay spread of the proposed GBSM aligns well with the practical measured data, which reflects the advantages of the GBSM proposed in this paper. However, the WINNER II channel model deviates greatly from the practical measurement data with the user trajectory. The reason may be the WINNER II channel model cannot accurately describe the propagation environment of cell-free massive MIMO channels while the proposed GBSM could.Figure 11
The delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array I.Figure 12
The delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array II.Figure 13
The delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array III.


### 5.2. The Doppler Spread

The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model are compared in Figures14–16 with configurations AP array I, AP array II, and AP array III, respectively. It is clear that the Doppler spread of the proposed GBSM aligns well with the practical measured data, which demonstrates that the GBSM was more practical to characterize real cell-free massive MIMO channels. However, the WINNER II channel model deviates greatly from the practical measurement data with the user trajectory.Figure 14
The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array I.Figure 15
The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array II.Figure 16
The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array III.


## 5.1. The Delay Spread

The time delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model are compared in Figure11, Figure 12, and Figure 13 with configurations AP array I, AP array II, and AP array III, respectively. It is clear that the delay spread of the proposed GBSM aligns well with the practical measured data, which reflects the advantages of the GBSM proposed in this paper. However, the WINNER II channel model deviates greatly from the practical measurement data with the user trajectory. The reason may be the WINNER II channel model cannot accurately describe the propagation environment of cell-free massive MIMO channels while the proposed GBSM could.Figure 11
The delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array I.Figure 12
The delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array II.Figure 13
The delay spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array III.


## 5.2. The Doppler Spread

The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model are compared in Figures14–16 with configurations AP array I, AP array II, and AP array III, respectively. It is clear that the Doppler spread of the proposed GBSM aligns well with the practical measured data, which demonstrates that the GBSM was more practical to characterize real cell-free massive MIMO channels. However, the WINNER II channel model deviates greatly from the practical measurement data with the user trajectory.Figure 14
The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array I.Figure 15
The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array II.Figure 16
The Doppler spread of the proposed GBSM, the practical measurement data, and the WINNER II channel model with AP array III.


## 6. Conclusion

In this paper, we proposed a GBSM for cell-free massive MIMO channels that combined the double-ring model and hemisphere model to comprehensively consider the distribution of scatterers in the environment. Combined with the LoS path, single scattering path, and double scattering path component, the channel matrix between APs and the user is derived. The proposed model fully considers geometric parameters such as the arrival/departure direction, elevation angle, and distance and can accurately characterize the channel. We have proved that the traditional channel model, standard block-fading model, and SBEM adopted in a cell-free massive MIMO system could not describe the nonstationarity in the space, time, and frequency domains, whereas the proposed GBSM could. The statistical characteristics of the channel, including the space CCF, time ACF, Doppler PSD, LCR, and AFD, are analyzed and simulated. By comparing these simulation results, we can see that the scattering environment, antenna space, and azimuth angle have an impact on all channel statistical properties. Specifically, the space CCF, time ACF, Doppler PSD, and AFD are relatively small in a rich scattering environment, while the LCR is relatively large. Since the space CCF and time ACF are relatively small and the values remain approximately unchanged in a rich scattering environment, it can be considered to conform to the WSS condition in a certain period of time. Moreover, as the antenna spacing and user azimuth angle decrease, the space CCF and time ACF gradually weaken. In summary, the simulations and analyses clearly showed that a rich scattering environment always obtains better channel performance than a low-density scattering case. To validate the correctness of the proposed GBSM, the measurement-based scenario and the WINNER II channel model were compared. Simulation results demonstrated that the GBSM was more practical to characterize real cell-free massive MIMO channels. Our research work can be considered theoretical guidance for establishing more purposeful cell-free massive MIMO measurement campaigns in the future [49].



---

*Source: 4730044-2022-10-17.xml*