# A Three-Dimensional Geometry-Based Statistical Model of2×2 Dual-Polarized MIMO Mobile-to-Mobile Wideband Channels

Authors: Jun Chen; Thomas G. Pratt
Journal: Modelling and Simulation in Engineering
 (2012)

---

## Abstract

A three-dimensional (3D) model for wide-band dual-polarized (DP) multiple-input-multiple-output (MIMO) mobile-to-mobile (M2M) channels is proposed. Using geometrical scattering based on concentric spheres at the transmitter (Tx) and at the receiver (Rx), a 3D parametric reference model for2×2 M2M DP multipath fading channels is developed. The channel model assumes the use of colocated half-wavelength dipole antennas for vertical and horizontal polarizations at both transmit and receive stations. Model parameters include the velocities of the Tx and Rx nodes, the distance between the nodes, the 3D antenna pattern gains, the azimuth and elevation angles of arrival and departure, the geometrical distribution of the scatterers, the Rician K-factors defining the fading envelope distributions, the maximum Doppler frequency, the scattering loss factors, the cross-polar power discrimination ratio (XPD), and the copolarization power ratio (CPR). Using the proposed model, expressions for joint time-frequency correlation functions (TFCFs) are derived which are used to investigate system behavior over different wide sense stationary uncorrelated scattering (WSSUS) channel realizations. The numerical results illustrate the sensitivities of the TFCF to simultaneous time and frequency offsets for the 2×2 DP-MIMO architectures.


---

## Body

## 1. Introduction

For nearly two decades, multiple-input-multiple-output (MIMO) communications systems [1–3] have been rigorously studied, leading to the adoption of MIMO in a number of wireless communications standards such as Long Term Evolution (LTE) and 802.11n. The LTE provides high-speed wireless communications for mobile phones and data terminals, while 802.11n supports MIMO, frame aggregation, and security improvements for wireless local area networks (WLAN). In these systems, multiple antennas can be used at both the transmitter (Tx) and receiver (Rx) to exploit the diversity and channel capacity afforded by MIMO architectures. In array-based MIMO systems, CP antennas at mobile stations and access points are ideally separated by at least one-half wavelength and at cellular base stations by at least ten wavelengths to achieve significant diversity or multiplexing gains [4]. In space-constrained implementations such as with handheld devices, use of antenna spacings less than one-half wavelength can lead to antenna correlations that lend themselves to beamforming gains, but ohmic coupling effects can deplete any gains that would otherwise be achieved with these strategies [5]. An alternative architecture for space-constrained deployments involves the use of dual-polarized (DP) antennas. This strategy has the advantage of offering colocated antennas that provide largely uncorrelated fading responses, although usually with asymmetrical average powers among the MIMO channel matrix entries. Hence, the channel behaviors for DP systems are substantively different from conventional copolarized (CP) channels, the latter which exhibit similar average powers. The relative performance of DP-MIMO in comparison to traditional MIMO with CP arrays will be dictated by the disparities in these responses. While DP architectures are receiving increasing attention, they have thus far been studied with much less rigor than CP systems.To characterize DP systems, especially in wideband signaling applications, modeling of dispersive input-to-output polarimetric channel behavior is important. Many theoretical models [6–11] provide statistical representations of polarimetric channels based on cross-polarization coupling (XPC), but do not characterize polarization-frequency behavior that is needed to represent wider-band, frequency-selective channels for the analysis of DP communications systems. For example, XPC is most often analyzed in the narrowband sense and yet it has been shown that the XPC can vary with frequency. It also exhibits correlation in the frequency dimension in a manner that depends on the channel’s temporal dispersion properties and polarization coupling [12]. Statistical channel models in literature that address time-selective and frequency-selective behavior often assume that the time dispersion (due to time delays) and the frequency dispersion (due to Doppler spread) are statistically independent, where the time delays depend on the relative locations of the random scatterers and the Doppler spreads depend on the motion of the Tx and Rx antennas [13, 14]. However, this assumption is broken if both the time delays and the Doppler spreads depend on the relative location of the random scatterers (e.g., the angles of departure and the angles of arrivals). Because these phenomena would not be expected to be statistically independent, it is of interest to examine joint correlation functions in time and frequency, which is an application of the model that we consider in our paper.Analytical approaches such as Kronecker or eigenbeam models have also been used for model simulations, and these have dealt with narrowband channels. An analytical framework based on a Kronecker model is presented in [15] to model narrowband DP Rayleigh and Rician fading channels for arbitrary array sizes. The framework uses a relatively small number of physical parameters to analyze the benefits of multiple polarization architectures. In other work, a 3D polarized narrowband spatial channel model is presented in [7], and the impact of elevation angle on capacity is reported in [16]. These models assume ideal dipole antennas and unity cross-polar discrimination ratios.Another approach to characterize channels for DP systems involves the use of geometric scattering models. Geometrical scattering models [17–19] have the advantage that they are able to represent important channel behaviors, including mutual effects that arise from polarization coupling, time-dispersion, angle-dispersion, and frequency-dispersion. However, they also are more complex. Early two-dimensional (2D) geometric scattering models were developed for narrowband CP single-input single-output (SISO) mobile-to-mobile (M2M) Rayleigh fading channels [20–22]. In [23], the polarisation-sensitive geometric modelling is developed with a direction-of-arrival (DOA) distribution that depends on the polarisation states of the transmitting antennas, the receiving antennas and the polarisation properties of the scatterers. These models were 2D in the sense that the models treated electromagnetic propagation only in a fixed-elevation plane and did not model elevation antenna pattern dependencies or distributions of scatterers in the elevation dimension.Some three-dimensional (3D) models have been proposed to overcome shortcomings of 2D M2M channel models, particularly in environments where deployed antenna heights are lower than surrounding buildings and obstacles. In [18], a 3D wideband M2M mathematical reference model was proposed based on a concentric-cylinders geometry using a superposition of line-of-sight (LOS), single bounce at the transmit side (SBT), single bounce at the receive side (SBR), and double-bounce (DB) rays in a variety of urban environments. The analysis approach proved useful to characterize M2 M correlations, but to date has only been applied to the analysis of CP MIMO channels with unity-gain idealized dipole antennas. Recently, a 3D polarized channel model has been proposed to treat spatially-separated orthogonally polarized elements [19]. The model deals exclusively with the double-bounce ray and narrowband channels.In this paper, we develop a 3D geometric scattering model for2×2 M2M DP wideband channels based on concentric spheres to evaluate joint time frequency correlation functions associated with the subchannel fading envelopes. The use of concentric spheres is motivated by the ease with which scatterer locations may be identified for a given radius, azimuth AOD, and elevation AOD, and by its capability to support the modeling of overhead reflectors. The DP-MIMO channel is constructed assuming a DP antenna at both the Tx and the Rx, where the polarization basis at the Rx is matched to the polarization basis used at the Tx. The received signals can be translated to an arbitrary orthogonally-polarized basis by applying a unitary transformation to the received signal vectors [24] without impacting theoretical performance measures such as capacity and diversity.We assume a channel with wide sense stationary uncorrelated scattering (WSSUS) where the channel correlation function is invariant over time, and the scatterers with different path delays are uncorrelated [25, 26]. Channel transfer functions for each transmit antenna/receive antenna pair are derived as a superposition of LOS, SBT, SBR, and DB component rays. The transfer functions are used to form time-frequency correlation functions (TFCF) for the assumed 3D nonisotropic scattering environment, where scattering distributions are characterized in the azimuth dimension through the von Mises distribution and in the elevation dimension by a cosine distribution. The 3D scattering propagation model used to evaluated the TFCFs simulates the effects of the antenna pattern gains, the geometrical distribution of scatterers and the associated azimuth/elevation angles of arrival and departure, the K-factor of the fading distributions, the maximum Doppler frequency, the scattering loss factors, the cross-polar discrimination (XPD), and the copolarization ratio (CPR).The remainder of the paper is organized as follows. The 3D DP2×2 MIMO M2M model is presented in Section 2. In Section 3, we derive transfer functions and TFCFs for each subchannel of the DP-MIMO channel and for each ray path type for 3D nonisotropic scattering environments. Numerical results of the joint correlation functions associated with the 2×2 DP channel models are presented in Section 4. We conclude with a summary of our findings in Section 5.


## 2.2×2 DP-MIMO Channel Model Based on Concentric Spheres

The concentric sphere geometry-based scattering model assumes a mobile Tx and a mobile Rx, both equipped with DP antennas using a matched polarization basis. Radio propagation between the Tx and the Rx is characterized by 3D WSSUS under channels that can include line-of-sight (LOS) and non-line-of-sight (NLOS) components, where scattering centers in the latter case reside on concentric spheres about either the Tx, the Rx, or both. The linearly time-variant MIMO channel can be represented by a2×2 impulse response matrix in terms of time t and delay τ(1)H(t,τ)=[hqp(t,τ)hqp~(t,τ)hq~p(t,τ)hq~p~(t,τ)].Figure1 illustrates the concentric sphere model for a MIMO M2M channel with DP antennas. One DP antenna with vertically-polarized and horizontally-polarized components denoted by At(p) and At(p~), respectively, is located at the center of the Tx sphere (Ot). A second DP antenna with corresponding components denoted by Ar(q) and Ar(q~) is located at the center of the Rx sphere (Or). At the Tx, M fixed scatterers reside within the volume of a sphere defined by a radius Rt. The mth transmit scatterer is denoted by St(m), where 1≤m≤M and resides on the surface of a sphere with radius Rm, where Rm<Rt. Similarly at the Rx, N scatterers occupy the spherical volume with a radius Rr. The nth receive scatterer is denoted by Sr(n), where 1≤n≤N and resides on the surface of a sphere with radius Rn, where Rn<Rr. The set of scatterers is comprised of H scatterers (reflecting horizontally polarized waves) and V scatterers (reflecting vertical polarized waves). We assume that the distribution of the H scatterers and the V scatterers are identical and that the number of scatterers for each are the same, although this may not be true in general.Figure 1
Concentric-spheres for 3D colocated DP channel model. The signals from the Tx antennas arrive in parallel at the Rx antenna array.The center of the Tx sphere serves as the global origin of a rectangular coordinate system. At timet=0, the Rx is a distance D from the Tx with XYZ coordinates denoted by (Δx, Δy, Δz). The height difference between the Tx and the Rx antennas is included in the offset Δz. The symbols ϵp,m, ϵm,q, ϵp,n, ϵn,q, ϵm,n and ϵp,q denote distances d(At(p),St(m)), d(St(m),Ar(q)), d(At(p),Sr(n)), d(Sr(n),Ar(q)), d(St(m),Sr(n)), and d(At(p),Ar(q)), respectively, where d(·) denotes the distance between the two coordinates. The symbols θt(m),θt(n) are the elevation angles of departure (EAoD, relative to the XY plane) to the scatterers St(m) and Sr(n), respectively, whereas ϕt(m),ϕt(n) are the azimuth angles of departure (AAoD, in the XY plane relative to the Z-axis) to the scatterers St(m) and Sr(n) respectively. Similarly, the symbols θr(m),θr(n) denote the elevation angles of arrival (EAoA, relative to the x-y plane) reflected from the scatterers St(m) and Sr(n), respectively, whereas ϕr(m),ϕr(n) denote the azimuth angles of arrival (AAoA, in the XY plane relative to the z-axis) reflected from the scatterers St(m) and Sr(n) respectively. For the LOS component between the Tx and the Rx, the symbols θtLOS and ϕtLOS are the elevation angle of departure (relative to the XY plane) and the azimuth angle of departure in the XY plane relative to the Z-axis, respectively. Similarly, the symbols θrLOS and ϕrLOS are the elevation angle of arrival and the azimuth angle of arrival, respectively. The Tx and Rx are moving with speeds vt and vr in directions described by the elevation angles θvt and θvr relative to the XY plane, respectively, and by the azimuth angles ϕvt and ϕvr in the XY plane, respectively.


### 2.1. Scatterer Distributions

The directions of scatterers are described by azimuth and elevation angle distributions. Several different distributions, such as uniform, Gaussian, Laplacian, and von Mises, have been used in the literature to characterize the azimuth angles of departure and arrival. von Mises Fisher (VMF) distributions are effective to model both AAoD and AAoA for spatial fading correlation models [27] and have the advantage that the probability density function approximates many of distributions and admits closed-form solutions [18]. To simplify the modeling, it is assumed that the azimuth and elevation angles are independent, enabling the use of a product of distribution functions instead of a joint distribution. The scatterer distributions are synthesized using different angle distributions for the azimuth and the elevation dimensions. The Von Mises Fisher probability density function (pdf) distribution is used for azimuth dimensions, including for angles of departure (ϕt(m) and ϕr(n)) and angles of arrival (ϕt(n) and ϕr(m)). The von Mises pdf is defined as [28]
(2)f(ϕ)=exp[κcos(ϕ-μ)]2πI0(κ),ϕ∈[-π,π),
where μ∈[-π,π) is the mean value of the scatterer directions in the azimuth plane, κ controls the spread of scatterers around the mean μ, and I0(·) is the zero-order modified Bessel function of the first kind. When κ=0, f(ϕ)=1/2π is a uniform distribution, corresponding to isotropic scattering in azimuth. As κ increases, the scatterers become more clustered about the mean angle μ, and scattering is nonisotropic.The random elevation angles of departure and arrival can be characterized by a uniform, cosine or Gaussian distribution. In [18, 29], a cosine pdf is used as it may fit with the typical propagation in M2M communications, where the Tx and Rx are in motion and equipped with low elevation antennas. A cosine distribution employed for the elevation angles of departure and arrival is given by [29]:
(3)f(θ)={π4|θmax|cos(π2θθmax),|θ|≤|θmax|≤π2,0,otherwise.
where θmax is the maximum elevation angle which we assume to have a value near 20°. This maximum elevation angle is typical of M2M wireless communications where both the Tx and Rx are equipped with low elevation antennas.The above-mentioned distributions are used to represent the response from a cluster, which describes a group of scatterers located within an isolated solid angle, and the response from multiple clusters compose the aggregate response from the channel. In [27], a mixture of Von Mises distributions is proposed for modeling the 3D direction of scatterers in the presence of multiple clusters of scatterers over the propagation channel. Mathematically, each distribution in the mixture of M VMFs at the Tx or N VMFs at the Rx can be described as f(ϕt(m)|ϕ-t(m),κm)(1≤m≤M) or f(ϕr(n)|ϕ-r(n),κn)(1≤n≤N). Note that ϕ-t(m) or ϕ-r(n) is the mean azimuth angle of the mth or nth cluster; κm or κn is the concentration of the mth or nth cluster. Hence, the overall density functions of the mixture model consisting of M VMF distributions for the cluster azimuth angles ϕt(m) at the transmitter and N VMF distributions for the cluster azimuth angles ϕr(n) at the Rx can be described as [27]
(4)f-(ϕt)=∑m=1Mνm(fϕt(m)∣ϕ-t(m),κm),f-(ϕr)=∑n=1Nνn(fϕr(n)∣ϕ-r(n),κn),
where M and N the number of clusters at the Tx and Rx, respectively; νm or νn is defined as the a priori probability that mth or nth cluster was generated. Similarly, the overall density functions of the mixture model consisting of the cosine distributions for the elevation angles θt(m) at the Tx and θr(n) at the Rx can be expressed as
(5)f-(θt)=∑m=1Mνm(θt(m)∣fθ-t(m)),f-(θr)=∑n=1Nνn(θr(n)∣θ-r(n)),
where θ-t(m) and θ-r(n) are the mean elevation angle of the mth or nth cluster at the Tx and Rx, respectively.Signals from the Tx antenna elements that propagate directly to the Rx antenna elements form the LOS component. Signals reflected exclusively from the scatterers located around the Tx before arriving at the Rx antenna elements are collectively called the SBT component. Similarly, transmit signals reflected only by scatterers located around the Rx before arriving at the Rx antenna elements form the SBR component. The DB component is formed from the signals that are reflected from scatterers about both the Tx and the Rx before arriving at the Rx antenna elements. For each realization of the WSSUS channel, the channel impulse responses and TFCFs can be written as superpositions of the LOS, SBT, SBR, and DB signal components.The channel impulse responses for the subchannelsAt(p)→Ar(q)→,At(p~)→Ar(q), At(p)→Ar(q~) and At(p~)→Ar(q~) can be written as a superposition of the LOS, SBT, SBR and DB signals,(6)hab(t,τ)=habSBR(t,τ)+habSBR(t,τ)+habDB(t,τ)+habLOS(t,τ),
where the habl(t,τ)(a∈{p,p~},b∈{q,q~},l∈{LOS,  SBT,  SBR,  DB}) represent the impulse response functions between antenna elements a and b along the l path respectively, and hab(t,τ) represents the total response between antenna elements a and b.


### 2.2. Single Bounce Channel Impulse Response

The time varying impulse responses of the single-bounce components are given by(7)habSBT(t,τ)=ηt,pqK+11N∑n=1Mξa,n,bga,n,b(t)×δ(τ-τa,m,b)×[Ga(v)(ϕt(m),θt(m))Ga(h)(ϕt(m),θt(m))]×[1rpq~(m)rp~q(m)1][Gb(v)(ϕr(m),θr(m))Gb(h)(ϕr(m),θr(m))]×exp(jϕab(m))ρmSBT,habSBR(t,τ)=ηt,pqK+11N∑n=1Mξa,n,bga,n,b(t)×δ(τ-τa,n,b)×[Ga(v)(ϕt(n),θt(n))Ga(h)(ϕt(n),θt(n))]×[1rpq~(n)rp~q(n)1][Gb(v)(ϕr(n),θr(n))Gb(h)(ϕr(n),θr(n))]×exp(jϕab(n))ρnSBR,
where the parameter K is the Rician K factor which is a ratio of the power through the LOS path relative to the power through the scattered paths. The power ratio coefficients ηt,ab and ηr,ab (a∈{p,p~}, b∈{q,q~}) in (7) specify the fraction of the power of the single-bounce components along the subchannels At(a)→Ar(b) with respect to the average transmit power Pab between the Tx antenna element a and the Rx antenna element b. The relationships among these, and related power ratio coefficients are shown in (29). ξa,m,b and ξa,n,b denote the signal amplitudes and τa,m,b and τa,n,b denote time delays of the multipath components propagating from antenna element a to antenna element b via transmit scatterer m and from antenna element a to antenna element b via receive scatterer n (a∈{p,p~}, b∈{q,q~}), respectively.G
a
(
v
)
(
·
,
·
), Ga(h)(·,·), Gb(v)(·,·), and Gb(h)(·,·)(a∈{p,p~},b∈{q,q~}) denote the antenna patterns of the pth, p~th, qth, and q~th antenna element with vertical(v) or horizontal(h) polarizations, respectively. We assume that a half-wavelength dipole is used for each element where the spherical coordinate system defining antenna orientations are illustrated in Figure 2. The antenna feed is situated at the origin of the coordinate system and the antenna elements are inclined at the angle α from the Z axis in the vertical ZX plane. A thin dipole is assumed and the element radius is ignored. The radiation power gain patterns of the half-wavelength dipole for vertical and horizontal polarizations are given by [30]
(8)G(v)(θ,ϕ)=1.641(cosθcosϕsinα-sinθcosα)2cos2(πζ/2)(1-ζ2)2,G(h)(θ,ϕ)=1.641sin2ϕsin2αcos2(πζ/2)(1-ζ2)2,
where
(9)ζ=sinθcosϕsinα+cosθcosα.Figure 2
Half-wavelength dipole antenna and its coordinate system.For illustration, antenna gain patterns of the half-wavelength dipole with inclination angles of0 and π/2 are shown in Figures 3 and 4 respectively.Figure 3
Gain patterns of half-wavelength dipole antenna with inclination angle of 0:G(v)(θ,ϕ). The horizontal component G(h)(θ,ϕ) is zero.Gain patterns of half-wavelength dipole antenna with inclination angle ofπ/2: (a) G(v)(θ,ϕ), and (b) G(h)(θ,ϕ).
(a)
(b)τ
a
,
m
,
b and τa,n,b denote time delays of the multipath components propagating from antenna element At(a) to antenna element Ar(b) via the scatterer St(m) and from antenna element At(a) to antenna element Ar(b) via the scatterer Sr(n), respectively. ga,m,b(t) and ga,n,b(t) are the time-varying phase functions along the path At(a)→St(m)→Ar(b) and a→Sr(n)→b, respectively. Since the antenna elements are collocated and share a common phase center, the time delays τa,m,b or τa,n,b are equivalent between antennas elements At(a) and Ar(b), that is,
(10)τp,m,q=τp,m,q~=τp~,m,q=τp~,m,q~,τp,n,q=τp,n,q~=τp~,n,q=τp~,n,q~.
Similarly, ga,m,b(t) and ga,n,b(t) satisfy
(11)gp,m,q(t)=gp,m,q~(t)=gp~,m,q(t)=gp~,m,q~(t),gp,n,q(t)=gp,n,q~(t)=gp~,n,q(t)=gp~,n,q~(t).The random variableϕab(m) represents the uniformly-distributed random phase offset associated with the path to each scatterer St(m) between the V or H component of the Tx antenna elements and the V or H component of the Rx antenna elements. Similarly, the random variable ϕab(n) is the uniformly-distributed phase offset associated with the path to each scatterer Sr(n) between the V or H component of the Tx antenna elements and the V or H component of the Rx antenna elements.r
p
~
q
(
m
), rpq~(n), rp~q(n), and rp~q(n) represent the inverse of the channel cross-polar discrimination ratios along the single bounce path to scatterer m or scatterer n, and have statistical means that follow(12)E{rpq~(m)}=E{rpq~(n)}=1XPDp,E{rp~q(m)}=E{rp~q(n)}=1XPDp~.Here, the channel cross-polar discrimination for each transmit polarization is defined as(13)XPDp=E|hqp|2E{|hq~p|2},XPDp~=E{|hq~p~|2}E{|hq~p~|2},
where E{·} denotes the expectation operator.The XPD depends on channel parameters and the environment, such as the distance between the Tx and the Rx, the angles of arrival and departure (both azimuth and elevation), the delay spread of the multipath components, and the transmit and receive antenna polarization basis.We define a parameterCPRqpq~p~  as the copolar power ratio between the average powers transmitted through the vertical-vertical subchannel and the average powers transmitted through the horizontal-horizontal subchannel:
(14)CPRqp,q~p~=E{|hq~p~|2}E{|hq~p~|2}.
The CPR depends on the Brewster angle phenomenon [31].The XPD and CPR, when expressed in decibel (dB), are often observed as having the normal distribution withN(μ,σ). In [7], the mean of XPD varies from 0 to 18 dB, with the standard deviation in order of 3~8 dB. Normally, the received power in the vertical-to-vertical transmission is reported to be greater than that in the horizontal-to-horizontal transmission (CPR>0 dB). For example, the mean of the CPR is reported to vary between 0 and 6 dB [15]. Depending on the propagation environment and transmission configuration, the CPR may be less than 0 dB when the amplitude of vertically polarized waves is degraded more than that of horizontally polarized waves or the transmission power of horizontally polarized waves is greater than that of vertically polarized waves. When the XPD or CPR in dB is denoted as X, the expectations in (12) can be computed as
(15)E{rpq~(m)}=E{rp~q(m)}=E{10-(X/20)},E{1CPRqp,q~p~}=E{10-(X/20)}.WhenX has the normal distribution with N(μ,σ), then the expectation E{10-X/20} over an interval [x1, x2] can be expressed as
(16)E{10-(X/20)}=∫x1x210-(x/20)12πσ2exp(-(x-μ)22σ2)dx.For simplicity, we adopt the normal XPD model to characterize cross-polarized components in the 3D modeling of2×2 DP-MIMO M2M channels.When the elevation angles of arrival/departure are small, the distancesd(At(p),St(m)), d(St(m),Ar(q)), d(At(p),Sr(n)), d(Sr(n),Ar(q)), d(St(m),Sr(n)) and d(At(p),Ar(q)) can be approximated as below
(17)|d(At(p),St(m))|=ϵp,m≈Rt,|d(St(m),Ar(q))|=ϵm,q≈D,|d(At(p),Sr(n))|=ϵp,n≈D,|d(Sr(n),Ar(q))|=ϵn,q≈Rr,|d(At(p),Ar(q))|=ϵpq≈D,|d(St(m),Sr(n))|=ϵm,n≈(Rt-Rr)2+D2.Furthermore, using the Taylor series expansion, the amplitudes of the multipath components,ξp,m,q or ξp,m,q~ is approximately given by
(18)ξp,m,q=ξp,m,q~=[|d(At(p),St(m))|+|d(St(m),Ar(q))|]-γ/2≈Am(1-γ2RtD),
where γ is the path loss exponent, whose value is normally in the range of 2 to 4, depending upon the propagation environment. ξp~,m,q~ or ξp~,m,q is approximated as
(19)ξp~,m,q~=ξp~,m,q=1CPRqp,q~p~[|d(At(p),St(m))|+|d(St(m),Ar(q))|]-γ/2≈Am1CPRqp,q~p~(1-γ2RtD),
where Am is defined as
(20)Am=D-γ/2.Similarly,ξp,n,q or ξp,n,q~ is approximated as
(21)ξp,n,q=ξp,n,q~=[|d(Ar(p),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am(1-γ2RrD),ξp~,n,q~ or ξp~,n,q is approximated as
(22)ξp~,n,q~=ξp~,n,q=1CPRqp,q~p~[|d(Ar(p),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am1CPRqp,q~p~(1-γ2RrD).The time delaysτp,m,q and τp,n,q are the travel times of the signals scattered from the Tx scatterer St(m), and from the Rx scatterer Sr(n), respectively
(23)τp,m,q=ϵp,m+ϵm,qc,(24)τp,n,p=ϵp,n+ϵn,qc,
where c is the speed of the light.The time-varying phase functionga,m,b(t) along the path At(p)→St(m)→Ar(q) is given by
(25)gp,m,q(t)=exp[-j2πλc(ϵp,m+ϵm,q)]×exp[j2πtftmaxcos(ϕt(m)-ϕvt)cosθt(m)cosθvt]×exp(j2πtftmaxsinθt(m)sinθvt)×exp[j2πtfrmaxcos(ϕr(m)-ϕvr)cosθr(m)cosθvr]×exp(j2πtfrmaxsinθr(m)sinθvr),
where ftmax=vt/λc and frmax=vr/λc are the maximum.Doppler frequencies associated with the Tx and Rx respectively andλc denotes the carrier wavelength.The time-varying phase functiongp,m,b(t) along the path At(p)→Sr(n)→Ar(q) is given by
(26)gp,n,q(t)=exp[-j2πλc(ϵp,n+ϵn,q)]×exp[j2πtftmaxcos(ϕt(n)-ϕvt)cosθt(n)cosθvt]×exp(j2πtftmaxsinθt(n)sinθvt)×exp[j2πtfrmaxcos(ϕr(n)-ϕvr)cosθr(n)cosθvr]×exp(j2πtfrmaxsinθr(n)sinθvr).The parametersρmSBT and ρnSBR represent scattering loss factors, and are governed by [32]
(27)ρ=exp[-8(πσhsinθiλc)2]I0[8(πσhsinθiλc)2],
where σh is the standard deviation of the surface height of the scatterer, θi is a given angle of incidence, and I0(·) is the zero-order modified Bessel function of the first kind.


### 2.3. Double-Bounce Channel Impulse Response

The double-bounce components of the channel impulse responses are(28)hpqDB(t,τ)=ηtr,pqK+11MN∑m=1M∑n=1Nξa,m,n,bga,m,n,b(t)×δ(τ-τa,m,n,b)×[Ga(v)(ϕt(m),θt(m))Ga(h)(ϕt(m),θt(m))]×[1rpq~(m)rp~q(m)1][1rpq~(n)rp~q(n)1]×[Gb(v)(ϕr(n),θr(n))Gb(h)(ϕr(n),θr(n))]×exp(jϕab(m,n))ρmSBTρnSBR,
where the power ratio coefficients ηtr,ab(a∈{p,p~},b∈{q,q~}) in (28), specify the ratio of the power of the double-bounce ray along the subchannels At(a)→Ar(b) with respect to the averaged transmit power Pab between the Tx antenna element a and the Rx antenna element b. The power ratio coefficients of the single-bounce and the double-bounce rays satisfy
(29)ηt,ab+ηr,ab+ηtr,ab=1,(a∈{p,p~},b∈{q,q~}).ξa,m,n,b denotes the amplitude and τa,m,n,b denotes the time delay of the multipath components along the path At(a)→St(m)→Sr(n)→Ar(b) (a∈{p,p~}, b∈{q,q~}). The amplitudes of the multipath components ξp,m,n,q and ξp,m,n,q~ are given by
(30)ξp,m,n,q=ξp,m,n,q~×[|d(At(p),St(m))|+|d(St(m),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am(1-γ2Rt+RrD).The amplitudes of the multipath componentsξp~,m,n,q~ and ξp~,m,n,q are given by
(31)ξp~,m,n,q~=ξp~,m,n,q=1CPRqp,q~p~×[|d(At(p),St(m))|+|d(St(m),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am1CPRqp,q~p~(1-γ2Rt+RrD).As DP antenna pair are colocated at the center of Tx and the center of Rx without considering the size of antennas, the time delaysτa,m,n,b are equivalent between antennas At(a) and Ar(b)(32)τp,m,n,q=τp,m,n,q~=τp~,m,n,q=τp~,m,n,q~.The time delayτp,m,n,q is the travel time of the signals scattered from both the Tx scatterer St(m) and the Rx scatterer Sr(n)(33)τp,m,n,q=ϵp,m+ϵm,n+ϵn,qc.Similarly,ga,m,n,b(t) satisfy
(34)gp,m,n,q(t)=gp,m,n,q~(t)=gp~,m,n,q(t)=gp~,m,n,q~(t).The time-varying functiongp,m,n,q(t) is given by
(35)gp,m,n,q(t)=exp[-j2πλc(ϵp,m+ϵm,n+ϵn,q)]×exp[j2πtftmaxcos(ϕt(m)-ϕvt)cosθt(m)cosθvt]×exp(j2πtftmaxsinθt(m)sinθvt)×exp[j2πtfrmaxcos(ϕr(n)-ϕvr)cosθr(n)cosθvr]×exp(j2πtfrmaxsinθr(n)sinθvr).The random variableϕab(m,n) is the phase offset of the double-bounce path scattered by both the Tx scatterer St(m) and the Rx scatterer Sr(n) between the V or H component of the Tx antenna elements and the V or H component of the Rx antenna elements.


### 2.4. LOS Channel Impulse Response

The LOS components of channel impulse responses are(36)habLOS(t,τ)=KK+1ξa,bLOSga,bLOS(t)δ(τ-τa,bLOS)×Ga(ϕtLOS,θtLOS)Gb(ϕrLOS,θrLOS),
where ξa,bLOS (a∈{p,p~}, b∈{q,q~}) denotes the amplitude and τa,bLOS denotes the time delay of the LOS components between the antenna element At(a) at the Tx and the antenna element Ar(b) at the Rx. The time-varying function ga,bLOS(t) is defined as
(37)gp,qLOS(t)=gp,q~LOS(t)=gp~,qLOS(t)=gp~,q~LOS(t)=exp(-j2πλcϵp,q)×exp[j2πtftmaxcos(ϕtLOS-ϕvt)cosθtLOScosθvt]×exp(j2πtftmaxsinθtLOSsinθvt)×exp[j2πtfrmaxcos(ϕrLOS-ϕvr)cosθrLOScosθvr]×exp(j2πtfrmaxsinθrLOSsinθvr).The amplitude of the LOS componentsξa,bLOS is given by
(38)ξp,qLOS=ξp,q~LOS≈D-γ2=Am,ξp~,q~LOS=ξp~,qLOS≈1CPRqp,q~p~D-γ2=Am1CPRqp,q~p~.The time delayτa,bLOS is the travel time of the signal from the Tx antenna element At(a) and the Rx antenna element Ar(b).  Consider the  following:
(39)τp,qLOS=τp,q~LOS=τp~,qLOS=τp~,q~LOS=ϵp,qc.It is assumed that the elevation angles (θtLOS, θrLOS, θt(m), θt(n), θr(m), and θr(n)) and the azimuth angles (ϕtLOS, ϕrLOS, ϕt(m), ϕt(n), ϕr(m) and ϕr(n)) are independent random variables. The radii Rt of the Tx sphere and Rr of the Rx sphere are also independent. The phase offsets ϕab(m), ϕab(n), and ϕab(m,n) are assumed to be uniformly random variables on the interval [-π,π) that are independent from the elevation angles, azimuth angles, and radii of the scattering spheres. Using the Central Limit Theorem [33, 34], we posit that the delay-spread functions habSBT(t,τ), habSBR(t,τ), habLOS(t,τ), and habDB(t,τ) are zero-mean complex Gaussian random processes.


### 2.5. Time-Variant Transfer Functions

The time-variant transfer function is the Fourier transform of the channel impulse response with respect to the delayτ. Consider the following:
(40)Tab(t,f)=Fτ{hab(t,τ)}=TabSBT(t,f)+TabSBR(t,f)+TabDB(t,f)+TabLOS(t,f),
where TabSBT(t,f), TabSBR(t,f), TabDB(t,f), and TabLOS(t,f) (a∈{p,p~},b∈{q,q~}) are time-variant transfer functions for the SBT, SBR, DB and LOS components, respectively.


## 2.1. Scatterer Distributions

The directions of scatterers are described by azimuth and elevation angle distributions. Several different distributions, such as uniform, Gaussian, Laplacian, and von Mises, have been used in the literature to characterize the azimuth angles of departure and arrival. von Mises Fisher (VMF) distributions are effective to model both AAoD and AAoA for spatial fading correlation models [27] and have the advantage that the probability density function approximates many of distributions and admits closed-form solutions [18]. To simplify the modeling, it is assumed that the azimuth and elevation angles are independent, enabling the use of a product of distribution functions instead of a joint distribution. The scatterer distributions are synthesized using different angle distributions for the azimuth and the elevation dimensions. The Von Mises Fisher probability density function (pdf) distribution is used for azimuth dimensions, including for angles of departure (ϕt(m) and ϕr(n)) and angles of arrival (ϕt(n) and ϕr(m)). The von Mises pdf is defined as [28]
(2)f(ϕ)=exp[κcos(ϕ-μ)]2πI0(κ),ϕ∈[-π,π),
where μ∈[-π,π) is the mean value of the scatterer directions in the azimuth plane, κ controls the spread of scatterers around the mean μ, and I0(·) is the zero-order modified Bessel function of the first kind. When κ=0, f(ϕ)=1/2π is a uniform distribution, corresponding to isotropic scattering in azimuth. As κ increases, the scatterers become more clustered about the mean angle μ, and scattering is nonisotropic.The random elevation angles of departure and arrival can be characterized by a uniform, cosine or Gaussian distribution. In [18, 29], a cosine pdf is used as it may fit with the typical propagation in M2M communications, where the Tx and Rx are in motion and equipped with low elevation antennas. A cosine distribution employed for the elevation angles of departure and arrival is given by [29]:
(3)f(θ)={π4|θmax|cos(π2θθmax),|θ|≤|θmax|≤π2,0,otherwise.
where θmax is the maximum elevation angle which we assume to have a value near 20°. This maximum elevation angle is typical of M2M wireless communications where both the Tx and Rx are equipped with low elevation antennas.The above-mentioned distributions are used to represent the response from a cluster, which describes a group of scatterers located within an isolated solid angle, and the response from multiple clusters compose the aggregate response from the channel. In [27], a mixture of Von Mises distributions is proposed for modeling the 3D direction of scatterers in the presence of multiple clusters of scatterers over the propagation channel. Mathematically, each distribution in the mixture of M VMFs at the Tx or N VMFs at the Rx can be described as f(ϕt(m)|ϕ-t(m),κm)(1≤m≤M) or f(ϕr(n)|ϕ-r(n),κn)(1≤n≤N). Note that ϕ-t(m) or ϕ-r(n) is the mean azimuth angle of the mth or nth cluster; κm or κn is the concentration of the mth or nth cluster. Hence, the overall density functions of the mixture model consisting of M VMF distributions for the cluster azimuth angles ϕt(m) at the transmitter and N VMF distributions for the cluster azimuth angles ϕr(n) at the Rx can be described as [27]
(4)f-(ϕt)=∑m=1Mνm(fϕt(m)∣ϕ-t(m),κm),f-(ϕr)=∑n=1Nνn(fϕr(n)∣ϕ-r(n),κn),
where M and N the number of clusters at the Tx and Rx, respectively; νm or νn is defined as the a priori probability that mth or nth cluster was generated. Similarly, the overall density functions of the mixture model consisting of the cosine distributions for the elevation angles θt(m) at the Tx and θr(n) at the Rx can be expressed as
(5)f-(θt)=∑m=1Mνm(θt(m)∣fθ-t(m)),f-(θr)=∑n=1Nνn(θr(n)∣θ-r(n)),
where θ-t(m) and θ-r(n) are the mean elevation angle of the mth or nth cluster at the Tx and Rx, respectively.Signals from the Tx antenna elements that propagate directly to the Rx antenna elements form the LOS component. Signals reflected exclusively from the scatterers located around the Tx before arriving at the Rx antenna elements are collectively called the SBT component. Similarly, transmit signals reflected only by scatterers located around the Rx before arriving at the Rx antenna elements form the SBR component. The DB component is formed from the signals that are reflected from scatterers about both the Tx and the Rx before arriving at the Rx antenna elements. For each realization of the WSSUS channel, the channel impulse responses and TFCFs can be written as superpositions of the LOS, SBT, SBR, and DB signal components.The channel impulse responses for the subchannelsAt(p)→Ar(q)→,At(p~)→Ar(q), At(p)→Ar(q~) and At(p~)→Ar(q~) can be written as a superposition of the LOS, SBT, SBR and DB signals,(6)hab(t,τ)=habSBR(t,τ)+habSBR(t,τ)+habDB(t,τ)+habLOS(t,τ),
where the habl(t,τ)(a∈{p,p~},b∈{q,q~},l∈{LOS,  SBT,  SBR,  DB}) represent the impulse response functions between antenna elements a and b along the l path respectively, and hab(t,τ) represents the total response between antenna elements a and b.


## 2.2. Single Bounce Channel Impulse Response

The time varying impulse responses of the single-bounce components are given by(7)habSBT(t,τ)=ηt,pqK+11N∑n=1Mξa,n,bga,n,b(t)×δ(τ-τa,m,b)×[Ga(v)(ϕt(m),θt(m))Ga(h)(ϕt(m),θt(m))]×[1rpq~(m)rp~q(m)1][Gb(v)(ϕr(m),θr(m))Gb(h)(ϕr(m),θr(m))]×exp(jϕab(m))ρmSBT,habSBR(t,τ)=ηt,pqK+11N∑n=1Mξa,n,bga,n,b(t)×δ(τ-τa,n,b)×[Ga(v)(ϕt(n),θt(n))Ga(h)(ϕt(n),θt(n))]×[1rpq~(n)rp~q(n)1][Gb(v)(ϕr(n),θr(n))Gb(h)(ϕr(n),θr(n))]×exp(jϕab(n))ρnSBR,
where the parameter K is the Rician K factor which is a ratio of the power through the LOS path relative to the power through the scattered paths. The power ratio coefficients ηt,ab and ηr,ab (a∈{p,p~}, b∈{q,q~}) in (7) specify the fraction of the power of the single-bounce components along the subchannels At(a)→Ar(b) with respect to the average transmit power Pab between the Tx antenna element a and the Rx antenna element b. The relationships among these, and related power ratio coefficients are shown in (29). ξa,m,b and ξa,n,b denote the signal amplitudes and τa,m,b and τa,n,b denote time delays of the multipath components propagating from antenna element a to antenna element b via transmit scatterer m and from antenna element a to antenna element b via receive scatterer n (a∈{p,p~}, b∈{q,q~}), respectively.G
a
(
v
)
(
·
,
·
), Ga(h)(·,·), Gb(v)(·,·), and Gb(h)(·,·)(a∈{p,p~},b∈{q,q~}) denote the antenna patterns of the pth, p~th, qth, and q~th antenna element with vertical(v) or horizontal(h) polarizations, respectively. We assume that a half-wavelength dipole is used for each element where the spherical coordinate system defining antenna orientations are illustrated in Figure 2. The antenna feed is situated at the origin of the coordinate system and the antenna elements are inclined at the angle α from the Z axis in the vertical ZX plane. A thin dipole is assumed and the element radius is ignored. The radiation power gain patterns of the half-wavelength dipole for vertical and horizontal polarizations are given by [30]
(8)G(v)(θ,ϕ)=1.641(cosθcosϕsinα-sinθcosα)2cos2(πζ/2)(1-ζ2)2,G(h)(θ,ϕ)=1.641sin2ϕsin2αcos2(πζ/2)(1-ζ2)2,
where
(9)ζ=sinθcosϕsinα+cosθcosα.Figure 2
Half-wavelength dipole antenna and its coordinate system.For illustration, antenna gain patterns of the half-wavelength dipole with inclination angles of0 and π/2 are shown in Figures 3 and 4 respectively.Figure 3
Gain patterns of half-wavelength dipole antenna with inclination angle of 0:G(v)(θ,ϕ). The horizontal component G(h)(θ,ϕ) is zero.Gain patterns of half-wavelength dipole antenna with inclination angle ofπ/2: (a) G(v)(θ,ϕ), and (b) G(h)(θ,ϕ).
(a)
(b)τ
a
,
m
,
b and τa,n,b denote time delays of the multipath components propagating from antenna element At(a) to antenna element Ar(b) via the scatterer St(m) and from antenna element At(a) to antenna element Ar(b) via the scatterer Sr(n), respectively. ga,m,b(t) and ga,n,b(t) are the time-varying phase functions along the path At(a)→St(m)→Ar(b) and a→Sr(n)→b, respectively. Since the antenna elements are collocated and share a common phase center, the time delays τa,m,b or τa,n,b are equivalent between antennas elements At(a) and Ar(b), that is,
(10)τp,m,q=τp,m,q~=τp~,m,q=τp~,m,q~,τp,n,q=τp,n,q~=τp~,n,q=τp~,n,q~.
Similarly, ga,m,b(t) and ga,n,b(t) satisfy
(11)gp,m,q(t)=gp,m,q~(t)=gp~,m,q(t)=gp~,m,q~(t),gp,n,q(t)=gp,n,q~(t)=gp~,n,q(t)=gp~,n,q~(t).The random variableϕab(m) represents the uniformly-distributed random phase offset associated with the path to each scatterer St(m) between the V or H component of the Tx antenna elements and the V or H component of the Rx antenna elements. Similarly, the random variable ϕab(n) is the uniformly-distributed phase offset associated with the path to each scatterer Sr(n) between the V or H component of the Tx antenna elements and the V or H component of the Rx antenna elements.r
p
~
q
(
m
), rpq~(n), rp~q(n), and rp~q(n) represent the inverse of the channel cross-polar discrimination ratios along the single bounce path to scatterer m or scatterer n, and have statistical means that follow(12)E{rpq~(m)}=E{rpq~(n)}=1XPDp,E{rp~q(m)}=E{rp~q(n)}=1XPDp~.Here, the channel cross-polar discrimination for each transmit polarization is defined as(13)XPDp=E|hqp|2E{|hq~p|2},XPDp~=E{|hq~p~|2}E{|hq~p~|2},
where E{·} denotes the expectation operator.The XPD depends on channel parameters and the environment, such as the distance between the Tx and the Rx, the angles of arrival and departure (both azimuth and elevation), the delay spread of the multipath components, and the transmit and receive antenna polarization basis.We define a parameterCPRqpq~p~  as the copolar power ratio between the average powers transmitted through the vertical-vertical subchannel and the average powers transmitted through the horizontal-horizontal subchannel:
(14)CPRqp,q~p~=E{|hq~p~|2}E{|hq~p~|2}.
The CPR depends on the Brewster angle phenomenon [31].The XPD and CPR, when expressed in decibel (dB), are often observed as having the normal distribution withN(μ,σ). In [7], the mean of XPD varies from 0 to 18 dB, with the standard deviation in order of 3~8 dB. Normally, the received power in the vertical-to-vertical transmission is reported to be greater than that in the horizontal-to-horizontal transmission (CPR>0 dB). For example, the mean of the CPR is reported to vary between 0 and 6 dB [15]. Depending on the propagation environment and transmission configuration, the CPR may be less than 0 dB when the amplitude of vertically polarized waves is degraded more than that of horizontally polarized waves or the transmission power of horizontally polarized waves is greater than that of vertically polarized waves. When the XPD or CPR in dB is denoted as X, the expectations in (12) can be computed as
(15)E{rpq~(m)}=E{rp~q(m)}=E{10-(X/20)},E{1CPRqp,q~p~}=E{10-(X/20)}.WhenX has the normal distribution with N(μ,σ), then the expectation E{10-X/20} over an interval [x1, x2] can be expressed as
(16)E{10-(X/20)}=∫x1x210-(x/20)12πσ2exp(-(x-μ)22σ2)dx.For simplicity, we adopt the normal XPD model to characterize cross-polarized components in the 3D modeling of2×2 DP-MIMO M2M channels.When the elevation angles of arrival/departure are small, the distancesd(At(p),St(m)), d(St(m),Ar(q)), d(At(p),Sr(n)), d(Sr(n),Ar(q)), d(St(m),Sr(n)) and d(At(p),Ar(q)) can be approximated as below
(17)|d(At(p),St(m))|=ϵp,m≈Rt,|d(St(m),Ar(q))|=ϵm,q≈D,|d(At(p),Sr(n))|=ϵp,n≈D,|d(Sr(n),Ar(q))|=ϵn,q≈Rr,|d(At(p),Ar(q))|=ϵpq≈D,|d(St(m),Sr(n))|=ϵm,n≈(Rt-Rr)2+D2.Furthermore, using the Taylor series expansion, the amplitudes of the multipath components,ξp,m,q or ξp,m,q~ is approximately given by
(18)ξp,m,q=ξp,m,q~=[|d(At(p),St(m))|+|d(St(m),Ar(q))|]-γ/2≈Am(1-γ2RtD),
where γ is the path loss exponent, whose value is normally in the range of 2 to 4, depending upon the propagation environment. ξp~,m,q~ or ξp~,m,q is approximated as
(19)ξp~,m,q~=ξp~,m,q=1CPRqp,q~p~[|d(At(p),St(m))|+|d(St(m),Ar(q))|]-γ/2≈Am1CPRqp,q~p~(1-γ2RtD),
where Am is defined as
(20)Am=D-γ/2.Similarly,ξp,n,q or ξp,n,q~ is approximated as
(21)ξp,n,q=ξp,n,q~=[|d(Ar(p),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am(1-γ2RrD),ξp~,n,q~ or ξp~,n,q is approximated as
(22)ξp~,n,q~=ξp~,n,q=1CPRqp,q~p~[|d(Ar(p),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am1CPRqp,q~p~(1-γ2RrD).The time delaysτp,m,q and τp,n,q are the travel times of the signals scattered from the Tx scatterer St(m), and from the Rx scatterer Sr(n), respectively
(23)τp,m,q=ϵp,m+ϵm,qc,(24)τp,n,p=ϵp,n+ϵn,qc,
where c is the speed of the light.The time-varying phase functionga,m,b(t) along the path At(p)→St(m)→Ar(q) is given by
(25)gp,m,q(t)=exp[-j2πλc(ϵp,m+ϵm,q)]×exp[j2πtftmaxcos(ϕt(m)-ϕvt)cosθt(m)cosθvt]×exp(j2πtftmaxsinθt(m)sinθvt)×exp[j2πtfrmaxcos(ϕr(m)-ϕvr)cosθr(m)cosθvr]×exp(j2πtfrmaxsinθr(m)sinθvr),
where ftmax=vt/λc and frmax=vr/λc are the maximum.Doppler frequencies associated with the Tx and Rx respectively andλc denotes the carrier wavelength.The time-varying phase functiongp,m,b(t) along the path At(p)→Sr(n)→Ar(q) is given by
(26)gp,n,q(t)=exp[-j2πλc(ϵp,n+ϵn,q)]×exp[j2πtftmaxcos(ϕt(n)-ϕvt)cosθt(n)cosθvt]×exp(j2πtftmaxsinθt(n)sinθvt)×exp[j2πtfrmaxcos(ϕr(n)-ϕvr)cosθr(n)cosθvr]×exp(j2πtfrmaxsinθr(n)sinθvr).The parametersρmSBT and ρnSBR represent scattering loss factors, and are governed by [32]
(27)ρ=exp[-8(πσhsinθiλc)2]I0[8(πσhsinθiλc)2],
where σh is the standard deviation of the surface height of the scatterer, θi is a given angle of incidence, and I0(·) is the zero-order modified Bessel function of the first kind.


## 2.3. Double-Bounce Channel Impulse Response

The double-bounce components of the channel impulse responses are(28)hpqDB(t,τ)=ηtr,pqK+11MN∑m=1M∑n=1Nξa,m,n,bga,m,n,b(t)×δ(τ-τa,m,n,b)×[Ga(v)(ϕt(m),θt(m))Ga(h)(ϕt(m),θt(m))]×[1rpq~(m)rp~q(m)1][1rpq~(n)rp~q(n)1]×[Gb(v)(ϕr(n),θr(n))Gb(h)(ϕr(n),θr(n))]×exp(jϕab(m,n))ρmSBTρnSBR,
where the power ratio coefficients ηtr,ab(a∈{p,p~},b∈{q,q~}) in (28), specify the ratio of the power of the double-bounce ray along the subchannels At(a)→Ar(b) with respect to the averaged transmit power Pab between the Tx antenna element a and the Rx antenna element b. The power ratio coefficients of the single-bounce and the double-bounce rays satisfy
(29)ηt,ab+ηr,ab+ηtr,ab=1,(a∈{p,p~},b∈{q,q~}).ξa,m,n,b denotes the amplitude and τa,m,n,b denotes the time delay of the multipath components along the path At(a)→St(m)→Sr(n)→Ar(b) (a∈{p,p~}, b∈{q,q~}). The amplitudes of the multipath components ξp,m,n,q and ξp,m,n,q~ are given by
(30)ξp,m,n,q=ξp,m,n,q~×[|d(At(p),St(m))|+|d(St(m),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am(1-γ2Rt+RrD).The amplitudes of the multipath componentsξp~,m,n,q~ and ξp~,m,n,q are given by
(31)ξp~,m,n,q~=ξp~,m,n,q=1CPRqp,q~p~×[|d(At(p),St(m))|+|d(St(m),Sr(n))|+|d(Sr(n),Ar(q))|]-γ/2≈Am1CPRqp,q~p~(1-γ2Rt+RrD).As DP antenna pair are colocated at the center of Tx and the center of Rx without considering the size of antennas, the time delaysτa,m,n,b are equivalent between antennas At(a) and Ar(b)(32)τp,m,n,q=τp,m,n,q~=τp~,m,n,q=τp~,m,n,q~.The time delayτp,m,n,q is the travel time of the signals scattered from both the Tx scatterer St(m) and the Rx scatterer Sr(n)(33)τp,m,n,q=ϵp,m+ϵm,n+ϵn,qc.Similarly,ga,m,n,b(t) satisfy
(34)gp,m,n,q(t)=gp,m,n,q~(t)=gp~,m,n,q(t)=gp~,m,n,q~(t).The time-varying functiongp,m,n,q(t) is given by
(35)gp,m,n,q(t)=exp[-j2πλc(ϵp,m+ϵm,n+ϵn,q)]×exp[j2πtftmaxcos(ϕt(m)-ϕvt)cosθt(m)cosθvt]×exp(j2πtftmaxsinθt(m)sinθvt)×exp[j2πtfrmaxcos(ϕr(n)-ϕvr)cosθr(n)cosθvr]×exp(j2πtfrmaxsinθr(n)sinθvr).The random variableϕab(m,n) is the phase offset of the double-bounce path scattered by both the Tx scatterer St(m) and the Rx scatterer Sr(n) between the V or H component of the Tx antenna elements and the V or H component of the Rx antenna elements.


## 2.4. LOS Channel Impulse Response

The LOS components of channel impulse responses are(36)habLOS(t,τ)=KK+1ξa,bLOSga,bLOS(t)δ(τ-τa,bLOS)×Ga(ϕtLOS,θtLOS)Gb(ϕrLOS,θrLOS),
where ξa,bLOS (a∈{p,p~}, b∈{q,q~}) denotes the amplitude and τa,bLOS denotes the time delay of the LOS components between the antenna element At(a) at the Tx and the antenna element Ar(b) at the Rx. The time-varying function ga,bLOS(t) is defined as
(37)gp,qLOS(t)=gp,q~LOS(t)=gp~,qLOS(t)=gp~,q~LOS(t)=exp(-j2πλcϵp,q)×exp[j2πtftmaxcos(ϕtLOS-ϕvt)cosθtLOScosθvt]×exp(j2πtftmaxsinθtLOSsinθvt)×exp[j2πtfrmaxcos(ϕrLOS-ϕvr)cosθrLOScosθvr]×exp(j2πtfrmaxsinθrLOSsinθvr).The amplitude of the LOS componentsξa,bLOS is given by
(38)ξp,qLOS=ξp,q~LOS≈D-γ2=Am,ξp~,q~LOS=ξp~,qLOS≈1CPRqp,q~p~D-γ2=Am1CPRqp,q~p~.The time delayτa,bLOS is the travel time of the signal from the Tx antenna element At(a) and the Rx antenna element Ar(b).  Consider the  following:
(39)τp,qLOS=τp,q~LOS=τp~,qLOS=τp~,q~LOS=ϵp,qc.It is assumed that the elevation angles (θtLOS, θrLOS, θt(m), θt(n), θr(m), and θr(n)) and the azimuth angles (ϕtLOS, ϕrLOS, ϕt(m), ϕt(n), ϕr(m) and ϕr(n)) are independent random variables. The radii Rt of the Tx sphere and Rr of the Rx sphere are also independent. The phase offsets ϕab(m), ϕab(n), and ϕab(m,n) are assumed to be uniformly random variables on the interval [-π,π) that are independent from the elevation angles, azimuth angles, and radii of the scattering spheres. Using the Central Limit Theorem [33, 34], we posit that the delay-spread functions habSBT(t,τ), habSBR(t,τ), habLOS(t,τ), and habDB(t,τ) are zero-mean complex Gaussian random processes.


## 2.5. Time-Variant Transfer Functions

The time-variant transfer function is the Fourier transform of the channel impulse response with respect to the delayτ. Consider the following:
(40)Tab(t,f)=Fτ{hab(t,τ)}=TabSBT(t,f)+TabSBR(t,f)+TabDB(t,f)+TabLOS(t,f),
where TabSBT(t,f), TabSBR(t,f), TabDB(t,f), and TabLOS(t,f) (a∈{p,p~},b∈{q,q~}) are time-variant transfer functions for the SBT, SBR, DB and LOS components, respectively.


## 3. Polarization Matched Time-Frequency Correlation Functions

Wide-sense stationarity and uncorrelated scatterers are often assumed to be valid for mobile radio channels [18]. In this paper, time and frequency dispersion are modeled dependently over a wide sense stationary uncorrelated scattering channel for a 3D nonisotropic scattering environment. For such channels, the time-frequency correlation function is an effective way of characterizing the statistical dependencies in the temporal and frequency domains associated with the mobile-to-mobile channels. Since the polarization states of antennas at both the Tx and Rx are matched, we derive the matched polarization-basis time-frequency correlation function to show the relationship between the frequency and the time of the 3D statistical model for 2×2 DP antennas.For a WSSUS channel, the TFCF of the time-varying transfer functionTab(a∈{p,p~},b∈{q,q~}) is defined in terms of the time difference Δt and the frequency separation Δf(41)ρab(Δt,Δf)=E[Tab(t,f)*Tab(t+Δt,f+Δf)],
where a∈{p,p~},b∈{q,q~} and p, p~, q, q~ correspond to vertical and horizontal antenna polarizations. E[·] is the expectation operator, and * denotes the complex conjugate operation. The TFCFs in (41) can be rewritten as the superposition of the TFCFs of the SBT, SBR, DB and LOS components. Consider
(42)ρab(Δt,Δf)=ρabSBT(Δt,Δf)+ρabSBR(Δt,Δf)+ρabDB(Δt,Δf)+ρabLOS(Δt,Δf).In order to estimate the time frequency correlations, we make the assumption that the number of scatterers in the 3D reference model is infinite so that the discrete azimuth angles (ϕt(m), ϕr(m), ϕt(n) and ϕr(n)), the discrete elevation angles (θt(m), θr(m), θt(n), and θr(n)) and the radii Rt, Rr can be represented as continuous random variables with probability density functions of f(ϕr(s)), f(ϕt(s)), f(θr(s)), f(θt(s)), f(Rt), and f(Rr), respectively. The form of these pdf’s are as follows.The  Tx  azimuth angles are characterized by the von Mises pdf in (2) as
(43)f(ϕt(S))=exp[κt(S)cos(ϕt(S)-μt(S))]2πI0(κt(S)),
where μt(S)∈[-π,π) is the mean value of angle at which the Tx scatterers are distributed in the x-y plane and κt(S) controls the spread of Tx scatterers around the mean μt(S). Similarly, the Rx azimuth angles are characterized by
(44)f(ϕr(S))=exp[κr(S)cos(ϕr(S)-μr(S))]2πI0(κr(S)),
where μr(S)∈[-π,π) is the mean value of angle at which the Rx scatterers are distributed in the x-y plane and κr(S) controls the spread of Rx scatterers around the mean μr(S).The Tx and Rx elevation angles are characterized by the cosine pdf  in (3) as
(45)f(θt(S))={π4|θtm|cos(π2θt(S)θtm),|θt(S)|≤|θtm|≤π2,0,otherwise,f(θr(S))={π4|θrm|cos(π2θr(S)θrm),|θr(S)|≤|θrm|≤π2,0,otherwise,
where θtm is the maximum elevation angle for the Tx scatterers and θrm is the maximum elevation angle for the Rx scatterers.The pdfs used to characterize the radiiRt and Rr are given by
(46)f(Rt)=3Rt2Rt23-Rt13Rt1≤Rt≤Rt2,f(Rr)=3Rr2Rr23-Rr13Rr1≤Rr≤Rr2.The pdf of phase-offset random variablesϕm(pol), ϕab(m), ϕab(n), and ϕab(m,n) (a∈{p,p~} and b∈{q,q~}) is given by a uniform distribution on the interval [-π,π). Using these distributions, the TFCFs ρSBT(Δt,Δf), ρSBR(Δt,Δf), ρDB(Δt,Δf), and ρLOS(Δt,Δf) of the time-variant transfer functions for the corresponding SBT, SBR, DB and LOS components are derived below.Using these pdfs it is possible to formulate the TFCFs for each of the bounce path components. The TFCFs along the SBT path are given by(47)ρabSBT(Δt,Δf)=ηt,pqAm2K+1IabSBT(Δt,Δf),
where IabSBT(Δt,Δf) is defined as
(48)IabSBT(Δt,Δf)=18|θtm|I0(κt(s))∫Rt1Rt2∫-θtmθtm∫-ππ(1-γ2RtD)2×exp[-j2πcΔf(ϵp,m+ϵm,q)]3Rt2Rt23-Rt13×cos(π2θt(s)θtm)exp[κt(s)cos(ϕt(s)-μt(s))]×exp[j2πΔtftmaxcos(ϕt(s)-ϕvt)cosθt(s)cosθvt]×exp(j2πΔtftmaxsinθt(s)sinθvt)×exp[j2πΔtfrmaxcos(ϕr(s)-ϕvr)cos(θr(s))cosθvr]×exp(j2πΔtfrmaxsinθr(s)sinθvr)×exp[-8(πσhsinθt(i)λc)2]I0[8(πσhsinθt(i)λc)2]×E[1CPRqp,q~p~][Ga(v)(ϕt(s),θt(s))Ga(h)(ϕt(s),θt(s))]×[1E[rpq~(m)]E[rp~q(m)]1][Gb(v)(ϕr(s),θr(s))Gb(h)(ϕr(s),θr(s))]×dϕt(s)dθt(s)dRt.E[rpq~(m)] and E[rp~q(m)] represent the expectations of the random variables rpq~(m) and rp~q(m), respectively, over the interval [−20 dB, 20 dB]. They are calculated using the pdf of the normal XPD model in (16). E[CPRqp,q~p~] denotes the expectation of random variable of copolar ratio over the interval [0, 6 dB].At the Rx, the azimuth angleϕr(s) and elevation angle θr(s) induced by the scatterer St(m) are approximated as
(49)ϕr(s)≈3π2+RtDsinϕt(s),θr(s)≈RtDθt(s)+ΔzD.The angle of incidenceθt(i) with the surface of Tx scatterer can be determined by using trigonometric identities in the triangle Ot-St(m)-Or as follows:
(50)θt(i)≈12arccos(Rt2+dm,r2-D22Rtdm,r),
where dm,r is given by
(51)dm,r2=(Δx-Rtcosθt(s)cosϕt(s))2+(Δy-Rtcosθt(s)sinϕt(s))2+(Δz-Rtsinθt(s))2.The TFCFs along the SBR path are given by(52)ρabSBR(Δt,Δf)=ηr,abAm2K+1IabSBR(Δt,Δf),
where IabSBR(Δt,Δf) is defined as
(53)IabSBR(Δt,Δf)=18|θrm|I0(κr(s))∫Rr1Rr2∫-θrmθrm∫-ππ(1-γ2RrD2)2×exp[-j2πcΔf(ϵp,n+ϵn,q)]3Rr2Rr23-Rr13×cos(π2θr(s)θrm)exp[κr(s)cos(ϕr(s)-μr(s))]×exp[j2πΔtftmaxcos(ϕt(s)-ϕvt)cosθt(s)cosθvt]×exp(j2πΔtftmaxsinθt(s)sinθvt)×exp[j2πΔtfrmaxcos(ϕr(s)-ϕvr)cosθr(s)cosθvr]×exp(j2πΔtfrmaxsinθr(s)sinθvr)×exp[-8(πσhsinθr(i)λc)2]I0[8(πσhsinθr(i)λc)2]×E[1CPRqp,q~p~][Ga(v)(ϕt(s),θt(s))Ga(h)(ϕt(s),θt(s))]×[1E[rpq~(n)]E[rp~q(n)]1][Gb(v)(ϕr(s),θr(s))Gb(h)(ϕr(s),θr(s))]×dϕr(s)dθr(s)dRr,
where E[rpq~(n)] and E[rp~q(n)] represent the expectations of the random variables rpq~(n) and rp~q(n), respectively. They are calculated over the interval [−20 dB, 20 dB] using the pdf of the normal XPD model in (16).At the Tx, the azimuth angleϕt(s) and elevation angle θt(s) induced by the scatterer Sr(n) are approximated as
(54)ϕt(s)≈RrDsinϕr(s),θt(s)≈RrDθr(s)-ΔzD.The angle of incidenceθr(i) with the surface of Rx scatterer can be determined by using trigonometric identities in the triangle Ot-Sr(n)-Or as follows:
(55)θr(i)≈12arccos(Rr2+dt,n2-D22Rrdt,n),
where dt,n is given by
(56)dt,n2=(Δx+Rrcosθr(s)cosϕr(s))2+(Δy+Rrcosθr(s)sinϕr(s))2+(Δz+Rrsinθr(s))2.The TFCFs along the DB path are given by(57)ρabDB(Δt,Δf)=ηtr,abAm2K+1IabDB(Δt,Δf),
where IabDB(Δt,Δf) is defined as
(58)IabDB(Δt,Δf)=164|θtm||θrm|I0(κt(s))I0(κr(s))×∫Rt1Rt2∫Rr1Rr2∫-θtmθtm∫-θrmθrm∫-ππ∫-ππ(1-γ2Rr+RtD)2×exp[-j2πcΔf(ϵp,m+ϵm,n+ϵn,q)]×3Rr2Rr23-Rr133Rt2Rt23-Rt13cos(π2θt(s)θtm)cos(π2θr(s)θrm)×exp[κt(s)cos(ϕt(s)-μt(s))]exp[κr(s)cos(ϕr(s)-μr(s))]×exp[j2πΔtftmaxcos(ϕt(s)-ϕvt)cosθt(s)cosθvt]×exp(j2πΔtftmaxsinθt(s)sinθvt)×exp[j2πΔtfrmaxcos(ϕr(s)-ϕvr)cosθr(s)cosθvr]×exp(j2πΔtfrmaxsinθr(s)sinθvr)×exp[-8(πσhsinθt(i)λc)2]I0[8(πσhsinθt(i)λc)2]×exp[-8(πσhsinθr(i)λc)2]I0[8(πσhsinθr(i)λc)2]×E[1CPRqp,q~p~][Ga(v)(ϕt(s),θt(s))Ga(h)(ϕt(s),θt(s))]×[1E[rpq~(m)]E[rp~q(m)]1][1E[rpq~(n)]E[rp~q(n)]1]×[Gb(v)(ϕr(s),θr(s))Gb(h)(ϕr(s),θr(s))]dϕr(s)dϕt(s)dθr(s)dθt(s)dRrdRt.The TFCFs along the LOS path are written in forms as(59)ρpqLOS(Δt,Δf)=Am2K+1ILOS(Δt,Δf),whereILOS(Δt,Δf) is defined as
(60)ILOS(Δt,Δf)=Kexp(-j2πcΔfϵpq)×exp[j2πΔtftmaxcos(π-ϕr(LOS)-ϕvt)sinθvt]×exp[j2πΔtfrmaxcos(ϕr(LOS)-ϕvr)sinθvr].Numerical integration methods are used to evaluate these TFCFs in the simulation analysis section.


## 4. Numerical Results and Analysis

In this section, we present numerical results associated with the analysis of the2×2 DP-MIMO system. The matched polarization-basis time-frequency correlation functions are evaluated using the models described in Section 3. In the analysis, the propagation paths include LOS, SBT, SBR, and DB bounce paths, where the DB paths predominate over the SBT and SBR paths. The temporal dispersion is on the order of 600 ns, and the Doppler spread is approximately 200 Hz. The mean CP power ratio CPRqp,q~p~ is set to 0 dB so that the average powers of the V-V link and the H-H link are approximately identical, a condition that we recognize will not always hold in operational channels. Further, we also assume that the von Mises distributions for the azimuth angles and the cosine distributions for the elevation angles are equivalent for both H scatterers (reflecting horizontally polarized waves) and V scatterers (reflecting vertical polarized waves), and that the XPD’s probability distribution functions are identical for the H-V and V-H components. As before, these are assumed for expediency. The specific simulation parameters are listed in Table 1. The center frequency is assumed to be 1 GHz.Table 1
The list of parameters used in the simulation.
Parameter
K
γ
f
t
max [Hz]
f
r
max[Hz]
σ
h [mm]
Value
0
3
100
100
0.2
Parameter
Δ
x [m]
Δ
y [m]
Δ
z [m]
f
c [GHz]
η
t
Value
1500
0
0
1.0
0.2
Parameter
θ
v
t
θ
v
r
ϕ
v
t
ϕ
v
r
η
r
Value
pi/4
pi/4
π
/
4
π
/
4
0.2
Parameter
XPD-μ
XPD-σ
CPR-μ
CPR-σ
η
t
r
Value
10 dB
6 dB
0 dB
3 dB
0.6
V-V/H-H link
Parameter
R
t
1[m]
R
t
2[m]
R
r
1 [m]
R
r
2 [m]
κ
t
(
s
)
Value
10
100
10
100
3
Parameter
κ
r
(
s
)
μ
t
(
s
)
μ
r
(
s
)
θ
t
m [deg]
θ
r
m [deg]
Value
3
π
/
4
7
π
/
4
20
20
V-H/H-V link
Parameter
R
t
1 [m]
R
t
2 [m]
R
r
1 [m]
R
r
2 [m]
κ
t
(
s
)
Value
10
100
10
100
3
Parameter
κ
r
(
s
)
μ
t
(
s
)
μ
r
(
s
)
θ
t
m [deg]
θ
r
m [deg]
Value
3
π
/
4
7
π
/
4
20
20Figures5, 6, 7, 8, and 9 provide examples of joint TFCFs that were computed for specific subchannels and bounce path components. ρpqSBT(Δt,Δf), ρpqDB(Δt,Δf), and ρpq(Δt,Δf) are shown for the V-V CP components in Figures 5, 6, and 7, respectively. Similarly, ρp~qSBT(Δt,Δf) and ρp~qDB(Δt,Δf) are shown in Figures 8, and 9 for the H-V link. The correlations along the Δt axis are driven by the Doppler frequency spread [35] and correlations along the Δf axis are driven by the delay spread associated with each component. The results indicate that the TFCFs are not separable into a product of a time correlation function and a frequency correlation function, a conclusion that we anticipated since the scatter locations contribute to both the time dispersion and the Doppler spread. The peak correlation value at Δt=Δf=0 in each plot is determined by the relative power of the paths being considered. The physical model parameters that play a role in generating temporal and Doppler spread include the radii Rr1, Rr2, Rt1, and Rt2 associated with the Rx and Tx spheres of scatterers, the maximum elevation angles θrm and θtm, the mean value of the azimuth angles μr(s) and μt(s), the spread factor of scatterers κr(s) and κt(s), the angles defining the direction of the velocity, θvr, ϕvr, θvt, and ϕvt, and the maximum Doppler frequencies frmax and ftmax.Figure 5
Normalized time-frequency correlation function of V-V link for SBT ray.Figure 6
Normalized time frequency correlation function of V-V link for DB ray.Figure 7
Normalized time frequency correlation function of V-V link for the superposition of SBT, SBR and DB rays.Figure 8
Normalized time-frequency correlation function of V-H link for SBT ray.Figure 9
Normalized time-frequency correlation function of V-H link for DB ray.Frequency correlation functions (for a fixedΔt) are driven primarily by the delay profile of the channel. Small delay profiles yield frequency correlation functions that decorrelate at higher frequencies. The DB component will exhibit the lowest decorrelation frequency since this component has the most significant delay spread. In the concentric-sphere based scattering model shown in Figure 1, the maximum delay spread is determined by the maximum distance of signal traveling along the double-bounce path induced by the Tx scatterer or the Rx scatterer. It can be roughly calculated as (2Rt+2Rr)/c, where c is the speed of light. The coherence bandwidth Wc in Hz is given approximately by the inverse proportion of the time delay spread τd. The minimum coherence bandwidths versus the radii of the Tx and Rx spheres used in the channel model are estimated and presented in Figures 11 and 12. The coherence bandwidth decreases when the radii of the Tx or Rx spheres increase. In Figure 7, if we approximate the coherence bandwidth as the frequency separation over which the magnitude of correlation coefficients is larger than a half of the its maximum at the zero frequency separation, the coherence bandwidth for the test profile used in our analysis is approximately 2.5 MHz.The time correlation function dependencies on the maximum Doppler frequency are illustrated for the double-bounce ray withΔf=0, in Figure 10. Large Doppler frequencies lead to small decorrelation times. We anticipate that the double-bounce components will decorrelate at least as quickly as the single-bounce components due to the fact that the Doppler spread of this component will be greater than or equal to the Doppler spread of the SBT component, a trend that is confirmed in Figures 8 and 9.Figure 10
Normalized time-correlation functions varying with the time delay and the maximum Doppler Tx/Rx frequencies for DB ray. The configuration is defined in Table1 with ηt=ηr=ηtr=1/3.Figure 11
The minimum coherence bandwidth of2×2 DP system varies with upper radii of spheres of Tx and Rx scatterers respectively.Figure 12
The minimum coherence bandwidth of2×2 DP system varies with upper radii of spheres of scatterers in both Tx and Rx.The use of non-zero XPD values leads to cross-polarized components (V-H and H-V) with power contributions lower than the CP components (V-V or H-H). Due to the similarity of parameters used to define the co- and crosspolarized channels, the cross-polarized response is approximately related to the CP ones by a scale factor. This result can be observed from cross-polarized parts of TFCFs in Section3. The normalized correlation values for the SBT and DB components at Δf=Δt=0 are plotted as a function of the mean XPD value μ for both matched polarization and cross polarization components in Figure 13. The normalized correlation values for the matched polarizations are constant since the correlations and powers are independent of μ. For the crosspolarization components, the normalized correlation is seen to decrease as μ increases, largely due to changes in relative power. At μ=0, when the scatterers yield the same average power at both matched and cross-polarized components, the normalized correlation is identical for both the matched and cross-polarization links.Figure 13
Normalized time-correlation functions of SBT and DB rays vary with the mean XPDμ and the configuration is set up in Table 1.


## 5. Conclusion

In this paper, a 3D geometrical propagation model has been proposed for DP-MIMO mobile-to-mobile Rayleigh and Rician fading channels. The parametric channel model incorporates a number of physical parameters to characterize the channel, including antenna patterns, azimuth/elevation angles of arrival and departure, geometrical distribution of scatterers, K-factor, the maximum Doppler frequency, scattering loss factor, crosspolar discrimination, and copolarization power ratio. Using the key parameters of the model, matched polarization-basis time-frequency correlation functions were formulated and numerically computed for WSSUS 3D non-isotropic scattering environments. The numerical results show that the joint TFCFs are not separable into independent time-correlation and frequency-correlation functions. The normalized TFCFs can be parsed into bounce path components for matched polarization and copolarization links, leading to power-normalized marginal TFCFs for the channel realization. The flexibility of the model enables control of channel parameters to achieve a variety of multipath fading environments to investigate2×2 DP architectures.Our intention in future research is to validate the proposed channel model using measurements. A series of experiments are planned to verify and justify our model.



---

*Source: 756508-2012-11-08.xml*