Analysis of Joint Angular Distribution for Nonreciprocal Beams via the Mixture of Gaussian Distribution Based on Ray-Tracing

Multiple-input and multiple-output (MIMO) technology can not only provide huge data rates but also overcome the severe propagation attenuation effect, especially in millimeter-wave (mmWave) bands by utilizing beamforming. The nonreciprocal beam is a novel transmission pattern, which indicates that transceivers adopt asymmetrical beamwidths. Such a special pattern can achieve fast beam alignment and alleviate equipment costs. Thorough knowledge of the corresponding wireless channel is pivotal to the system design and optimization, which remains to be investigated. In this paper, we first propose a 3-dimensional (3-D) channel model based on ray-tracing, which is capable of reflection simulation. Based on this model, the ray-based beamforming mechanism is illustrated. The angular distribution is pivotal to beam channel modeling and characterization since transceiver beams filter rays in the angular domain. Then, we conduct an omnidirectional antenna-based channel simulation in an urban macro-cell scenario via the ray-tracing platform. On this basis, we focus on the distribution of the quasiangles, i.e., angles between departure/arrival reflected rays and the line-of-sight (LoS) path. We find that gamma distribution is a better option to fit the quasiarrival angular distribution than the von Mises distribution. Furthermore, to characterize the relationship between quasiangles of departure (AoD) and quasiangles of arrival (AoA), the Gaussian mixture model (GMM) is adopted and the ex-pectation-maximization (EM) algorithm is used to estimate the unknown parameters of GMM. Our findings provide useful insights to beam channel modeling, which should take the joint angular distribution into consideration.


Introduction
With the burst use of video streaming, immersive VR/AR, and the Internet-of-ings (IoT), the proliferation of the fth-generation (5G) mobile communication system boosts wireless mobile connectivity, enhances data rate, and yields great user satisfaction [1]. Enhanced mobile broadband (eMBB) communication is one of three primary 5G New Radio (NR) use cases de ned by the 3GPP, which aims to achieve a high peak user data rate of 1 Gbps. Moreover, this target is set to 1 Tbps in the next sixthgeneration (6G) communication system, which brings new challenges and receives extensive attention both from academia and industry. Under such a circumstance, the trend of exponentially growing demand in mobile data drives wireless networks to migrate to higher frequency bands.
Millimeter-wave (mmWave), spanning the spectrum from 30 GHz to 300 GHz, is envisioned as the key technology of the next-generation mobile communication system, including both beyond 5G (B5G) and 6G communication systems [2]. However, the mmWave band leads to unbearable propagation loss due to a variety of factors including atmospheric absorption, basic Friis transmission-e ect, and low penetration. To address this issue, it is widely believed that mmWave propagation together with the multiple-input and multiple-output (MIMO)-based beamforming technology can handle the drastic increase of mobile data tra c [3].
Beamforming employs facilities to vary the beams' power and their directions with the aid of amplitude/phase variations. Modern wireless system uses antenna arrays instead of a single antenna to synthesize a beam. Ideally, transceivers' beams point to each other all the time, i.e., beam alignment can provide a high propagation gain. However, perfect beam alignment incurs high delay overhead, especially over fast fading channel scenarios like highspeed railways (HSR), vehicle-to-vehicle/infrastructure (V2V/I) communications, and unmanned aerial vehicles (UAV). Conventional beam pattern adopts the symmetric beamwidth, reciprocal beams in brief, which makes nding a perfect match pair of beams time-consuming by exhaustive search method. Furthermore, the computational complexity explodes when using narrow-width beams. Besides, such a beam pattern requires the same channels of both transmitter (Tx) and receiver (Rx), which incurs substantial hardware complexity and expensive hardware costs. To cope with these issues, some prevailing technologies, such as hybrid beamforming (the combined analog and digital beamforming), low-resolution analog-to-digital converters (ADCs), and digital-to-analog converter (DACs) technologies, have been thoroughly investigated. Nevertheless, the former requires a heavy beam training overhead for channel estimation, whereas the latter uses complex iterative algorithms to recover the lost information at the receiver. e authors in [4,5] rst proposed a concept of the fulldigital arrays with nonreciprocal transmitting and receiving beam patterns. e novelty of this strategy lies in the fact that the nonreciprocal beam patterns use an unequal number of Tx and Rx channels, or beams with di erent widths, which makes the transmitting and receiving patterns asymmetrical. Figure 1 presents a diagram of this pattern. As seen in this gure, the base station (BS) uses a narrow-width beam to transmit signals in the downlink, whereas the mobile station (MS) adopts a wide-width beam to receive. In the uplink, the usage of beams is quite opposite. Bene ting from this exible architecture, this pattern owns merits like wide beam scanning range and low hardware complexity. erefore, it can achieve higher energy and cost e ciencies. is nonreciprocal pattern could not only release the burden of data processing as well as reduce the hardware cost and power consumption of the entire system but also can realize ecient beam alignment quickly. However, such a pattern leads to the nonreciprocity for uplink and downlink channels, indicating that time division duplex (TDD)-based systems cannot implement this strategy directly.
A thorough understanding of the wireless channel is a prerequisite for further system design, algorithm optimization, and network deployment. Generally, current prevailing channel models mainly contain the geometry-based stochastic models (GBSMs) and the ray-tracing at the mmWave bands. e former model is widely used due to its mathematical tractability, especially for the regular-shaped GBSM (RS-GBSM) [6]. Ray-tracing is validated by extensive comparisons with actual measured data [7]. Based on this, a huge volume of MIMO-related ray-tracing works are also performed. Nonetheless, some of them use the horn antenna to mimic the beamforming [8], whereas others consider the MIMO array response but fail to investigate the beam ltering e ect [9]. As for the beam channel measurement works, most of them are performed on symmetrical beam patterns, like [10,11]. To our best knowledge, few related works on the asymmetrical pattern have been studied until now, which remains a gap to be investigated thoroughly and we will address it herein.
To this end, we aim to propose a beam channel model based on the ray-tracing and investigate the stochastic channel properties of high re ection bounces. e motivations and our contributions are presented as follows: (1) We propose a 3-dimensional (3-D) MIMO channel model based on the ray-tracing, which incorporates the multiple-bounce re ection propagations. Our proposed model is elaborated from the phase expression of propagation distance and Doppler shift. en, we prove the feasibility of beamforming and enable our model capable of beamforming with di erent beamwidths and steering directions.
(2) To get the full knowledge of channel angular information of multiple-bounced paths, we analyze the angular distribution relationship based on a ray-tracing platform in an urban street scenario. Speci cally, we focus on transceiver angular distributions on the assumption of perfect beam alignment, i.e., the angles between Tx/Rx and departure/arrival rays, denoted as quasi-angles of arrival (AoAs)/departure (AoDs) in brief. (3) We focus on the arrival side and nd that the gamma distribution is a better option than the von Mises distribution to characterize the distribution of quasi-AoAs of multiple re ection bounces. Moreover, to obtain the angular distribution at the departure side, we employ the Gaussian mixture model (GMM) to depict the relationship between quasi-AoAs and quasi-AoDs. Besides, transceiver beams are considered as two constraints of the GMM. e remainder of this paper proceeds as follows. In Section 2, we present an overview of the mmWave beam channel modeling. Section 3 introduces the proposed channel model in detail. In Section 4, we analyze the angular distribution based on a ray-tracing platform, and the nal relevant conclusions are drawn in Section 5.

Literature Review
In terms of mmWave beam channel measurement, numerous works have been performed with respect to different frequency carriers. ese works can be generally classified into two types based on beam synthesis method, i.e., hornbased measurement and phased-array-based measurement.
e first method emulates a beam with a certain beam width by using a horn antenna, which is widely used due to its simple implementation. Related works are like [12][13][14][15]. Indeed, this method can work to a large extent by designing an antenna pattern similar to an array response. However, the beam alignment needs to be achieved via the real-time adjustment of antenna rotation.
Other works were conducted based on the phased-array antennas, which form real beams in the desired direction with appropriate power. A list of related mmWave channel measurement campaigns is presented in Table 1. However, due to the high equipment cost and expenditure, few works have been performed, not to mention the complete channel characterization of different beamwidths. Hence, there is a great demand for an accurate theoretical model with the ability of beamforming simulation. e measurement-based channel modeling indeed provides some insightful mmWave beam models, but it cannot cover a wide range of beamwidth, whereas theoretical models can fill this research gap. Currently, the theoretical methods of mmWave channel modeling can be classified into two types, i.e., the deterministic and stochastic methods. As for the former one, ray-tracing is widely used due to its high accuracy, especially considering that electronic wave exhibits an evident mechanism of reflection at high-frequency bands. Some works have combined ray-tracing with beamforming, like [19]. Moreover, some commercial simulation platforms, e.g., Wireless Insite, has already incorporated beamforming into ray-tracing [20]. However, RT is time-consuming when considering multiple bounces and cannot efficiently emulate the scattering mechanism.
In terms of the stochastic method, a prevailing scheme is the geometry-based stochastic method (GBSM) [21], which is widely adopted in some 5G organization standards, e.g., COST2100 and COST IC1004 channel models, QuaDRiGa, etc. [22]. Much work has been performed on the GBSM, representative works like [23][24][25]. Moreover, the application of GBSM involves rich scenarios, like V2V and UAV [26,27]. Conventional GBSM assumes that channel parameters follow a certain stochastic process, which imposes little computational burdens and can be quickly adapted to various scenarios by changing parameter values. Other works aim to apply beamforming into GBSM, like [28]. However, the AoD distributions of these works are obtained via geometric manipulations based on the singlebounced or double-bounced assumption, which faces exceptional geometric difficulties for higher bounces.
Another stochastic method is the propagation-graph (PG) theory, which abstracts the propagation environment into a directed graph [29]. is method requires much less computational complexity compared to RT. It can yield a satisfying simulation accuracy since that it mainly considers the scattering mechanism, enriching the "tails" of multipath taps in channel impulse response (CIR). Moreover, some works applied scattering and diffraction models into PG to enhance modeling accuracy [29,30]. Nevertheless, these works all fail to incorporate the beamforming effect while simulating the multiple-input and multiple-output (MIMO) system.
As for works on the nonreciprocal beam patterns, current efforts are mainly devoted to antenna design and implementation. A concept of a full-duplex nonreciprocalbeam-steering model is proposed based on the metasurface in [31], which uses nonreciprocal radiation beams to realize point-to-point telecommunications. Furthermore, the authors in [32] proposed a nonreciprocal-beam phased-array based on the transistor-based nonreciprocal phase shifters, which can exhibit different beams for transmission and reception states. Notably, the above works are performed on the metasurface, which is beyond the scope of commercial mobile communications. Since being investigated in [4], nonreciprocal beam patterns attract much attention due to the low hardware complexity of radio frequency (RF) chains. In [33], the authors designed an asymmetrical transceiver architecture for massive MIMO systems, which consists of a different number of RF chains for Tx and Rx. Reference [34] reports a high-performance wideband frequency synthesizer for nonreciprocal beams. However, to our best knowledge, little research on channel modeling and characterization over nonreciprocal beam patterns has been performed until now.  Parameter Definition T p , R q , r n e p-th Tx, q-th Rx, n-th ray D n (t) Path distances of the n-th ray Azimuth/elevation angles of Tx beam-steering direction Azimuth/elevation angles of MS velocity

Ray-Tracing-Based Channel Model.
For further derivations and descriptions, we list some important variables together with definitions in Table 2. e channel transfer function (CTF) obtained by ray-tracing is the sum of all possible reflected paths and line-of-sight (LoS) path.
Different from the conventional expression, we take the channel gain as well as the phase information into account. Without loss of generality, we use the word "ray" instead of path in the following text since that our proposal is inspired by ray-tracing. Specifically, the CTF of the n-th ray from the the p-th Tx (T p ) to the q-th Rx (R q ) can be formulated as where F R,q,V and F R,q,H are the radiation patterns of the q-th receiver antenna element for vertical and horizontal polarizations, respectively. F T,p,V and F T,p,H denote the radiation patterns of the p-th transmit antenna element for vertical and horizontal polarizations, respectively. Parameters Θ VV n , Θ VH n , Θ HV n , and Θ HH n refer to the initial phases of the n-th ray for 4 different polarization combinations, i.e., vertical-vertical (VV), vertical-horizontal (VH), horizontalvertical (HV), and horizontal-horizontal (HH). κ n is the corresponding cross polarization power ratio between H and V, and λ is the wavelength. D n means the propagation distance from Tx center to Rx center. D pq,n is the propagation distance of the n-th reflected path from T p to R q , and τ p,q,n is the corresponding propagation delay. g n is the path gain after M reflection bounces, which can be expressed as where ϑ m is the incident angle of the m-th reflection and Γ m is the corresponding Fresnel reflection coefficient concerning ϑ m . In equation (1), η n is the additional phase caused by reflection interaction and initial phase. Φ →D n is the departure angle unit vector of the n-th ray with the elevation angle of departure (EAOD) β V n and the azimuth angle of departure (AAOD) β H n . Φ →D bm is the beam-steering angle unit vector of Tx, which consists of EAOD β V bm and AAOD β H bm .
Similarly, the arrival angle unit vector Φ →A n of the n-th ray with the elevation angle of arrival (EAOA) α V n and the azimuth angle of arrival (AAOA) α H n can be marked as Parameters r → p � (x p ′ , y p ′ , z p ′ ) and r → q � (x q ″ , y q ″ , z q ″ ) denote the corresponding coordinates of T p and R q in the local coordinate systems (LCS), i.e., O ′ (x ′ , y ′ , z ′ ) and O ″ (x ″ , y ″ , z ″ ), respectively. ] n is the Doppler shift caused by the movement of mobile station (MS), which is calculated as where where variable ] MS means the MS speed.

MIMO-Based Beamforming Effects.
MIMO-based beamforming technology can overcome the severe propagation attenuation for mmWave communications. In this subsection, we aim to illustrate the mechanism of MIMOenabled beamforming. We use the phased array to synthesize the digital beams, which means that each antenna element is assigned with a RF converter. Besides, each channel is assigned with a weight and additional phase to synthesize a desirable beamwidth and beam-steering direction. On the assumption of far-field condition, we adopt the 2-dimensional uniform rectangular array (URA) and the phase parts A1 in equation (1) can be characterized by the Tx array steering vector, i.e., where m � 1, . . . , N H t , m ′ � 1, . . . , N V t . We assume that the origin of the LCS is located at the array center. k � 2π/λ is 4 International Journal of Antennas and Propagation the wave number. N H t , N V t mean the Tx and Rx antenna element number along the horizontal and vertical directions, respectively. Furthermore, each element channel is assigned with a weight to formulate an idea beam shape. Consequently, the synthesis of digital beamforming can be regarded as the design of a finite impulse response (FIR) filter. Herein, we adopt the Dolph-Chebyshev window function ω d � [ω n ] N×1 , which can achieve the lowest sidelobe level (SLL) for a given beamwidth among all weightbased beam designing methods [35]. It can be expressed as equation (8) In equation (8), function (·)! means the factorial function, and z is an auxiliary variable with an expression of e first case of z means that the main lobe width α W is a constant while designing, whereas the second case refers that the SLL is determined in decibel scale regardless the main lobe width. As for the desired beam-steering direction , ω d should be modified by p(α H s ) and p(α V s ), and the final corresponding expression of the beam response where i ∈ H, V { }, and symbol ⊙ means the element-wise operator. e major merits of the Dolph-Chebyshev window function consist of two aspects. Importantly, it can achieve the minimum average SLL while guaranteeing a constant main lobe width. Meanwhile, it can provide the narrowest beamwidth while guaranteeing a constant SLL. As such, we can roughly regard the filter gain as 1 within the passband frequency range and 0 for the stopband case, i.e., Finally, the final beam-enabled CTF is expressed as where Ω D n � α H n , α V n , Ω A n � β H n , β V n , and ⎦g n · e j2πD n /λ · e j2π] n t · e jη n δ τ − τ n . (13) It can be seen that the combined effect of all transceivers' elements is equivalent to the Dolph-Chebyshev filter. erefore, the antenna element index is omitted in this formula.

Ray-Tracing-Based
Channel Characterization

Distribution of Quasi-AoAs.
e above section elaborates on the proposed MIMO-enabled beamforming channel model based on ray-tracing. In this section, we will investigate the angular properties of the beam channel based on a commercial ray-tracer platform. It is noteworthy that the simulation is performed via this platform without considering the beam e ects, i.e., it is conducted based on the omnidirectional antenna.
e nonreciprocal beams are considered as two constraints based on the conclusions in equation (11). e simulation is conducted based on the Wireless Insite, a commercial ray-tracing platform with high simulation e ciency. As shown in Figure 2, the roadside BS communicates with the moving vehicle in an urban street. Furthermore, we simplify this cellular-assisted vehicle communication scenario and present in Figure 3. Note that single-input and single-output (SISO) omnidirectional antenna pattern is used in the simulation and the beam e ect is realized for each ray according to equation (11). e BS is located at the roadside with a height of 10 m, whereas the MS with an Rx mounted on the vehicle moves forwards along the main street at a speed of 60 km/h and a 147-meters-long street is simulated. Tall buildings are around the street with heights larger than that of BS. Consequently, scarce re ected rays come from the top obstacles and most of the received rays di er dramatically in the horizontal plane. Detailed values of the simulation parameters are presented in Table 3.
As mentioned earlier, beams act like a spatial lter that poses e ects on rays with di erent angular information. Figure 4 compares the ltered rays via nonreciprocal beams. Red lines indicate re ected rays with higher received power, whereas brown ones represent lower received power, and green lines mean the lowest power. Obviously, this special pattern yields di erent ltering results on re ected rays. erefore, to obtain the beam-enabled stochastic channel model, the primary pivotal task is to investigate the angular distribution of AoAs and AoDs. Figure 5 presents the  International Journal of Antennas and Propagation statistical histogram of the AoA results for all rays over this process. Besides, the widely used von Mises probability distribution function (PDF) is utilized to depict the AoA distribution and the tting result is plotted in this gure. Clearly, this curve shows a good tting trend to the simulated data. Considering that perfect beam alignment is achieved regardless of MS position and velocity, we focus on the angle between the LoS path and the arrival path from the last re ection bounce to the MS, abbreviated as quasi-AoA. For instance, as shown in Figure 3, α 1 and β 1 denote the AoD and AoA of the rst ray path, whereas α 1 and β 1 mean the corresponding quasi-AoD and quasi-AoA, respectively. Based on this novel de nition, the corresponding expression α n of the n-th path is where Φ →D s is the angle unit vector of Tx beam-steering direction with respect to α H s . Similarly, the expressions of Φ →A , Φ →A s , and β can also be derived. Figure 5(b) plots the statistics results of the quasi-AoAs over the MS movement process. It can be observed that the majority appears at lower values and a tail comes into view at higher values. Remember that all values of quasi-AoAs are larger than 0. Taking these two factors into account, the gamma distribution is a promising candidate. Hence, we attempt to t the results by using the gamma distribution function with a detailed expression as where a is the shape parameter, b is the scale parameter, and Γ(·) means the gamma function. As compared in Figure 5(b), the gamma distribution is a better option with a lower root mean squared error (RMSE) result. is angular distribution model provides insights into the angular generation for current prevailing channel models like GBSM.
e Doppler frequency shift is another pivotal parameter for a ray. It can be calculated as where β 0 means the AoA of the LoS link, symbol " ± " depends on the relative position of BS and MS, and c MS is the movement direction of MS.

Joint Distribution of Quasi-AoAs and Quasi-AoDs.
In this subsection, we aim to characterize the distribution of quasi-AoDs. We try to obtain its distribution based on the prior distribution of quasi-AoAs. To characterize the joint distribution relationship between quasi-AoAs and quasi-AoDs, we make statistics of those rays up to 3 re ection bounces since that those high-bounce rays contribute little power to the received signal. e statistical histogram results of these rays are scattered in Figure 6(a), where the abscissa and ordinate mean quasi-AoA and quasi-AoD, respectively. For a certain bin, the lighter color it is, the higher probability it means. From this gure, we can draw two conclusions. First, the movement of vehicles leads to discernible trajectories with di erent lifetimes, which demonstrates the channel evolution. Usually, it can be characterized from a perspective of the birth-death process [36]. Second, the distribution result indicates that quasi-AoAs and quasi-AoDs are not irrelevant to each other. Two dense clusters with irregular shapes can be recognized roughly. Herein, we adopt the GMM to depict this joint probability distribution, which can be expressed as where (18) and (19) denote the beamwidth constraints under the condition of perfect beam alignment. As illustrated in equation (11), beams act like spatial lters regarding rays with di erent AoAs and AoDs. erefore, only a few rays survive when considering beamforming as presented in Figure 4. Under this case, the statistical results will be less accurate as a result of few samples. Consequently, the GMMbased joint distribution of multiple rays is obtained based on the omnidirectional antenna, and the nonreciprocal beams are considered as constraints since that perfect beam alignment is assumed in this paper.
A data point consists of two variables x → ( α, β), K is the cluster number, which can be determined by the Akaike information criterion (AIC). e k-th component has a mixing proportion ρ k , a mean of μ k → , and a covariance matrix of Σ k for this multivariate case. Parameter θ k can be expressed as θ k μ k → , Σ k and ϕ( x → |θ k ) refers to the k-th component, which can be calculated by where symbol | · | means the determinant operator, and Σ −1 k indicates the inverse matrix of Σ k . To obtain these 3 key parameters (ρ k , μ → k , Σ k ), the maximum likelihood estimation (MLE) is adopted to estimate from the mixture PDF. Assume that X → 1 , . . . , X → n are observed random variables obeying Gaussian distribution with unknow parameters, then the likelihood function is given as  International Journal of Antennas and Propagation and the corresponding log likelihood function of equation (21) can be calculated as By taking partial derivatives of equation (22) to maximize the likelihood, these 3 parameters can be obtained. When estimating, the expectation-maximum (EM) algorithm is usually adopted. Initially, the parameters are randomly chosen for the mixture model parameters en, the parameters are updated in each iteration, including E-stimation and M-aximization steps, until the convergence criteria are satis ed, whereas the E-step calculates the membership coe cients for all data point (i 1, . . . , N) and mixture components (k 1, . . . , K) by utilizing the current parameters θ 1: K , i.e., where c ik denotes the possibility of data X → i belonging to the k-th Gaussian model with an obvious constraint K k 1 c ik 1. In terms of the M-step, it aims to calculate the coe cient and nd the parameters, which can be expressed as where E[·] is the mathematical expectation operator, and Var[·] means the calculation of variance. By adopting the GMM, Figure 6(b) presents the joint probability distribution histogram of quasi-AoAs and quasi-AoDs together with the GMM tting results expressed in the contour style. e cluster number is set to 2 and the corresponding estimated parameters are listed in Table 4.
Once the quasi-AoA value β of a ray is determined, the distribution of 2-dimensional GMM is reduced to 1-dimensional GMM. Consequently, the corresponding quasi-DoD value can be obtained by the Bayes estimation, i.e., f( α| β) P α| β, θ 1: K , where seed is a random number distributed evenly within the interval (0, 1), and CDF(·) means the cumulative distribution function (CDF) of f(·). Once the β is determined, the two-dimensional GMM is reduced to one-dimensional function. Hence, the CDF(·) is a sum of several error functions. Let ϕ( α|μ k| β , σ k| β ) denote the k-th Gaussian distribution component. en, the corresponding total CDF can be expressed as where erf(·) means the error function. Obviously, it is hard to nd the explicit expression of the inverse CDF; hence, numerical solution is used herein.

Joint Distribution of Quasi-AoAs and Propagation Delay.
Another crucial parameter for a path is the propagation delay/distance. In this section, we focus on characterizing the joint distribution between quasi-AoAs and the propagation distance. To better illustrate the scattered datapoints, we use the propagation delay in decibel instead of distance, and the corresponding results are presented in Figure 7. Generally speaking, these datapoints can be divided into 3 regions roughly by two dashed lines. e rst part is the single-bounce trajectory at the bottom, which stems from the ground re ection. is part is a deterministic process that can be easily derived by some geometric manipulations. ose datapoints located in the middle region present an ascending trend, which is characterized by a power function with a time-varying variance. It is noteworthy that these datapoints mainly stem from single-bounce and double-   International Journal of Antennas and Propagation bounce reflections. In terms of the rest datapoints on the top, they exhibit a drifting status and we use the Gaussian distribution to depict them. Table 5 presents a summary of these 3 regions.

Conclusions
In this paper, we propose a 3-D geometric beam channel model based on the ray-tracing and conduct the angular distribution analysis based on the simulation results for multiple reflection bounces. Results reveal that the arrival angles between reflected rays and the LoS path, abbreviated as quasi-AoA, can be characterized by the gamma distribution with high accuracy. en, the GMM is used to characterize the joint distribution between quasi-AoAs and quasi-AoDs. Furthermore, we analyze the relationship between quasi-AoAs and propagation delays. Consequently, nonreciprocal beam pattern channel modeling should take the joint angular distribution relationship into account seriously and our findings are important to beam channel models for next-generation mobile communications.

Data Availability
e simulation data used to support the findings of this study were supplied by Beijing Jiaotong University under license and so cannot be made freely available. Requests for access to these data should be made to Liu Liu, liuliu@bjtu.edu.cn.

Conflicts of Interest
e authors declare that they have no conflicts of interest.