Channel Measurement and Ray-Tracing-Statistical Hybrid Modeling for Low-Terahertz Indoor Communications
11 Channel Measurement andRay-Tracing-Statistical Hybrid Modeling forLow-Terahertz Indoor Communications
Yi Chen,
Student Member, IEEE , Yuanbo Li, Chong Han,
Member, IEEE ,Ziming Yu, and Guangjian Wang
Abstract
TeraHertz (THz) communications are envisioned as a promising technology, owing to its un-precedented multi-GHz bandwidth. One fundamental challenge when moving to new spectrum is tounderstand the science of radio propagation and develop an accurate channel model. In this paper, awideband channel measurement campaign between 130 GHz and 143 GHz is investigated in a typicalmeeting room. Directional antennas are utilized and rotated for resolving the multi-path components(MPCs) in the angular domain. With careful system calibration that eliminates system errors andantenna effects, a realistic power delay profile is developed. Furthermore, a combined MPC clusteringand matching procedure with ray-tracing techniques is proposed to investigate the cluster behaviorand wave propagation of THz signals. In light of the measurement results, physical parameters andinsights in the THz indoor channel are comprehensively analyzed, including the line-of-sight path loss,power distributions, temporal and spatial features, and correlations among THz multi-path characteristics.Finally, a hybrid channel model that combines ray-tracing and statistical methods is developed for THzindoor communications. Numerical results demonstrate that the proposed hybrid channel model showsgood agreement with the measurement and outperforms the conventional statistical and geometric-basedstochastic channel model in terms of the temporal-spatial characteristics.
Index Terms
Y. Chen, Y. Li and C. Han are with the Terahertz Wireless Communications Laboratory, University of Michigan-Shanghai Jiao Tong University Joint Institute, Shanghai Jiao Tong University, Shanghai, China (email: { yidoucy, yuanbo.li,chong.han } @sjtu.edu.cn).Z. Yu and G. Wang are with Huawei Technologies Co., Ltd, China (email: { yuziming, wangguangjian } @huawei.com). a r X i v : . [ c s . I T ] J a n Terahertz communications, Channel measurement, Temporal-Spatial analysis, Ray-Tracing-Statistichybrid modeling.
I. I
NTRODUCTION
In the past few decades, wireless data traffic has witnessed exponential growth and bringsan increasing demand for high-speed wireless communication. New spectral bands are requiredto support Terabit-per-second (Tbps) data rates for future wireless applications [1]. Currently,wireless local area networks (WLAN) techniques, i.e., 802.11ad protocol and fifth-generation(5G) mobile networks, have opened up the millimeter-wave (mmWave) spectrum (10-100 GHz)to seek for broader bandwidth and higher data rates [2], [3]. However, still limited to severalGHz bandwidth, the mmWave cannot support Tbps requirements. To further move up the carrierfrequency, the Terahertz (THz) band spanning over 0.1 and 10 THz, is envisioned as one ofthe promising spectrum bands to enable ultra-broadband 6G communications. With tens andhundreds of GHz contiguous spectrum resources, the use of THz band can address the spectrumscarcity and capacity limitations of current wireless systems. To make THz communicationscome into reality in the very near future, the first spectral window of the THz band is centeredat 140 GHz band, which can be implemented by mature electronic-based device technologies [4].When moving to the THz spectrum, fundamental challenges to designing THz wireless systemsinclude understanding the science of radio propagation, development of efficient yet accuratechannel models, and thorough analysis of channel characteristics in the THz band. To understandthe science of radio propagation, and develop channel models to describe major propagationprocesses, three main approaches are adopted, including physical measurement using a channelsounder [5], geometric-optic method [6], and full-wave EM field analysis [7]. In this paper, weadopt physical channel measurement to acquire the realistic THz propagation features betweenthe transmitter and receiver in a meeting room.From the literature, a number of channel measurement campaigns at THz frequencies have beenreported for short-distance indoor scenarios [8]–[14], outdoor vehicular communications [15]–[17], and aircraft cabin communications [18]. Most of these studies are above 300 GHz, which ishigher than 140 GHz that we concern. The channel measurement at 140 GHz for indoor scenariosincludes short-range on a desk [9], and in-building [12]. In particular, these studies mainly focuson the analysis of path loss, reflection loss, and penetration loss [9], [12]. However, an extensivechannel measurement with various transmitter (Tx) and receiver (Rx) positions in a meeting room environment for 140 GHz is still missing. Furthermore, thorough channel characterizationon the temporal-spatial THz propagation features as well as large-scale fading path loss is notwell captured, which includes K-factor, sparsity, temporal and spatial dispersion, correlation ofmulti-path features, among others.In this paper, we present a wideband channel measurement campaign based on the vectornetwork analyzer (VNA) at 130-143 GHz in a typical meeting room. With two fixed Tx, multi-path measurements are conducted at ten different receiver positions. Directional antennas areequipped at Tx and Rx to resolve multi-path components (MPCs) in the angular domain. Toassist post-processing of the measured data, we invoke ray-tracing (RT) techniques to simulatethe MPCs with detailed propagation trace in the same geometry where our channel measurementis implemented. Then, a combined MPCs clustering and matching procedure with a multi-pathcomponent distance (MCD) based Density-Based Spatial Clustering of Applications with Noise(DBSCAN) is proposed to cluster the measured MPCs and label the MPC clusters based on theRT results.In light of the measurement results, we characterize the physical parameters and provideinsights of the THz indoor channel. Specifically, the measured path loss of the line-of-sight (LoS)path is validated with theoretical computation results. Moreover, the power-delay-angular profiles(PDAPs) are presented to reveal the multi-path channel sparsity in the THz indoor environmentand the significance of reflection over walls. Furthermore, temporal and spatial features, includingthe K-factor, delay, and power spreads, are analyzed. A correlation matrix is formulated toanalyze the dependence among the multi-path characteristics of the THz indoor channel, e.g.,Tx-Rx separation distance, number of clusters, K-factor, delay spread, and angular spread. Atlast, we develop a hybrid semi-deterministic channel model that combines the geometric-opticmethod (ray tracing) and statistical modeling method, which suggests that the non-line-of-sight(NLoS) path reflects over the walls are dominating in the channel. Numerical results showthat the proposed hybrid model outperforms the conventional statistical and the 3rd GenerationPartnership Project (3GPP) geometric-based stochastic channel models (GSCM) models with themetric of PDAP.Compared to our preliminary and shorter version [19], this full-fledged work includes multi-path clustering effects and development of a RT-statistical hybrid channel model for the indoorTHz channel, with substantially more performance analysis. The distinctive contributions of thiswork are summarized as follows. • We establish a wideband THz channel measurement campaign in a meeting room, between130 GHz and 143 GHz via VNA. Especially, Rx with a directional antenna scans in bothazimuth ( [0 ◦ , ◦ ] ) and elevation domains ( [ − ◦ , ◦ ] ) with the rotation step of ◦ forresolving the MPCs in the spatial domain. 27 Tx-Rx positions are measured, including 21for channel analysis and model parameter extraction and 6 for model validation. • An MPCs clustering and matching procedure with an MCD-based DBSCAN algorithm isproposed to cluster the MPCs and match the measured MPCs with the simulated MPCswith the RT simulator. The proposed procedure outperforms conventional K-Power-Means(KPM) algorithm in our channel measurement data. • A comprehensive analysis of the channel characteristics is provided, including the path loss,reflection loss, cluster behavior, K-factor, delay spread, and angular spread. The correlationmatrix of the channel parameters is evaluated to discuss the temporal-spatial properties. • An RT-statistical hybrid channel model for the indoor THz channel at 140 GHz is developedon the basis of the analysis of the channel characteristics. The proposed hybrid model haslow complexity without requiring detailed geometry information of the environment. Delayand angular spreads of the proposed model are validated and shows good agreement withthe channel measurement. • A measure of figure structure similarity, Structural Similarity Index Measure (SSIM), isintroduced to measure the accuracy of the channel model along with the root mean squareerror (RMSE). SSIM and RMSE between the generated PDAP by the proposed hybrid modeland measured PDAP are evaluated, respectively. Compared to the conventional statisticalchannel model and the 3GPP GSCM [20], the proposed hybrid model shows improvedperformance in characterizing PDAP of the THz multi-path channel.The remainder of this paper is organized as follows. In Sec. II, we describe the channelmeasurement campaign, post-processing of measured data, and an MPCs clustering and matchingprocedure. Then, we present the channel measurement results and, more importantly, thoroughlyanalyze the channel characteristics in Sec. III. Furthermore, an RT-statistical hybrid channelmodel is proposed and validated in Sec. IV. Finally, the paper is concluded in Sec. V.II. C
HANNEL M EASUREMENT C AMPAIGN
In this section, we describe the THz measurement campaign, including the specification of thehardware system, meeting room environment, and measurement deployment. Moreover, system
Fig. 1: Channel measurement platform.calibration is carried out for eliminating the effects of the cables and radio-frequency (RF)fonts on the measured channel, as detailed in [1]. Ray-tracing techniques are used to trace thepropagation of measured MPCs. We further propose a clustering and MPC matching procedurefor analyzing MPCs clustering and match the measured data with the RT simulator output.
A. Channel Measurement System
The measurement platform mainly consists of 140 GHz RF fronts with horn antennas and aVNA. The local oscillator (LO) signal of 10.667 GHz is multiplied by a factor of 12 to 128GHz, and the intermediate frequency (IF) signals from 2 GHz to 15 GHz is up-converted by themultiplied LO signal to the frequency band from 130 to 143 GHz. As a result, the measuredbandwidth B w is 13 GHz, and the time domain resolution, ∆ τ = 1 /B w , is 76.9 ps, whichsuggests we can resolve two paths if the distance difference between them is larger than 2.3 cm.The recorded points in the frequency domain or equivalently, the sweeping frequency pointsare 1301, which corresponds to the sweeping interval at ∆ f = 10 MHz. The maximum excessdelay, τ m = 1 / ∆ f , is calculated as 100 ns, and accordingly, the largest traveling distance of adetectable path to be measured is L m = 30 m, which is sufficient for channel measurement ina meeting room.Both Tx and Rx are equipped with horn antennas, which are shown in Fig. 1. The antennaat Tx produces an antenna gain of 15 dBi with the half-power beamwidth (HPBW) of ◦ at140 GHz, to enable a sufficient angular coverage. The Rx antenna gain is 25 dBi, and the HPBWis ◦ , which is one-third of that at Tx for high spatial resolution. The two antennas are mountedon two rotation units, which can be rotated by step motors. In addition, the power of the testsignal is 1 mW, and the noise floor of our THz measurement platform is − dBm. In our measurement, the maximum received signal strength is − . dBm. The detailed parameters ofthe measurement system are summarized in Table I.TABLE I: Parameters of the Measurement System. Parameter Symbol ValueStart frequency f start
130 GHzEnd frequency f end
143 GHzBandwidth B w
13 GHzSweeping points N ∆ f
10 MHzAverage noise floor P N -160 dBmTest signal power P in HP BW Tx ◦ HPBW of receiver
HP BW Rx ◦ Antenna gain at Tx G t
25 dBiAntenna gain at Rx G r
15 dBiTime domain resolution ∆ τ ∆ L τ m
100 nsMaximum path length L m
30 mAzimuth rotation range [0 ◦ : 10 ◦ : 360 ◦ ] Elevation rotation range [ − ◦ : 10 ◦ : 20 ◦ ] B. Meeting Room Environment and Measurement Deployment
We carry out the channel measurement in a typical meeting room with an area of 10.15 m × × ◦ is fixed as an AP ata corner of the meeting room. The Tx height is 2 m, which is higher than the Rx or UE with aheight of 1.4 m. The wide HPBW for Tx guarantees the coverage of target UEs, while the small (a) Meeting room for THz channel measurements. (b) Overview of the meeting room and measurement deployment. Fig. 2: The meeting room layout and measurement deployment.HPBW for Rx of ◦ provides high spatial resolution for spatial-domain channel measurement.In our measurement deployment, 2 positions of Tx as well as 12 positions of Rx are set in themeeting room, as depicted in the top view of the meeting room in Fig. 2b. Tx1 is at the corner ofthe meeting room while Tx2 is behind the table and near wall 1. In this measurement campaign, two measurement sets are conducted. In the measurement set 1, Tx is placed on the position Tx1,the Rx is placed on the positions Rx1-4 and Rx6-10. In the measurement set 2, Tx is placed onthe position Tx2 while Rx is placed on the positions Rx1-12. For the measurement of each Rx,the main beam of Tx is directed to the Rx. By contrast, Rx with the spatial resolution of ◦ scans the receiving beam in the azimuth domain from ◦ to ◦ to detect sufficient multi-paths.As a directional beam of transmitter points to the Rx, the antenna gain of the reflected paths fromthe ceiling and the floor are 16 dB lower than that of the LoS path. Therefore, the consideredreflected paths collected in our experiment are mainly from the desk, chairs and walls, whoseelevation angles are sufficiently confined within [- ◦ , ◦ ]. C. Measured Power Delay Angular Profile
We perform Inverse Discrete Fourier Transform (IDFT) on the transfer function measured fromVNA to compute the channel impulse response (CIR) for each azimuth and elevation angle. As aresult, we can obtain the three-dimensional (3D) power-delay-angular profile (PDAP), p r ( i, j, k ) ,denoting the received power with time-of-arrival of i ∆ τ , azimuth angle-of-arrival of j ∆ θ , and elevation angle-of-arrival k ∆ φ . Then, we compose all the received power with the same azimuthangle and different elevation angles to obtain the two-dimensional (2D) PDAP of each Tx-Rxpair, calculated as PDAP ( i, j ) = (cid:88) k p r ( i, j, k ) . (1)In particular, we do not use the space-alternating generalized expectation-maximization (SAGE)algorithm to estimate the parameters of MPCs for eliminating the effect of antenna pattern andsimply regard each point on the PDAP is an MPC. The reason is that the SAGE algorithm in thechannel measurement with rotated directional-antenna is highly affected by phase inaccuracy,and requires accurate 3D antenna radiation pattern, which is challenging in practice and timeconsuming [21]. We take Rx6 as an example and present its path loss over the spectrum andthe composed 2D PDAP in Fig. 3(a) and 3(b), respectively.Since Tx and Rx are well synchronized by the VNA, the delay of the peak in the channelimpulse function refers to the absolute time of arrival of an MPC. To verify this, we observe thatthe delay of the first peak is 19.15 ns, which matches well with the theoretical time-of-arrival(ToA) of the LoS path of 18.4 ns, with the traveling distance of 5.52 m. The minor deviation iscaused by a small displacement of the Tx and Rx antennas. The reason is that the Rx antennais mounted on a rotator and the rotator is placed on a trolley. When we move the trolley to themarked position, the head of the antenna might not be precisely located at the desired position. D. MPCs Clustering and Matching with Ray-Tracing Simulator
MPCs clustering is to group the MPCs with similar ToA and angle-of-arrival (AoA) into acluster. Physically, a cluster consists of one reflected path along with some scattered paths sur-rounding the reflection point. Therefore, the cluster behavior is one of the most significant channelcharacteristics, based on which many statistical channel models and GSCMs are developed [22]–[25]. To cluster MPCs, the unsupervised learning techniques, e.g., K-means and Gaussian MixtureModel (GMM), K-Power-Means (KPM) [26], Kernel-power-density (KPD) [27] method areutilized. However, those methods except for KPD require manual selection of cluster numberand are sensitive to the initial cluster centers. Especially, K-means is limited at handling withthe clusters with a ball shape and not flexible for others. Moreover, KPD cannot identify theoutliers which are common in our channel measurement results.As the channel measurement only tells the amplitude, angle, and delay information of theMPCs, the physical propagation information of each path is still unknown to us. To solve this
130 132 134 136 138 140 142Spectrum [GHz]89.59090.59191.59292.593 P a t h l o ss [ d B ] (a) Path loss over the spectrum of Rx6 in measurement set1. (b) Composed 2D PDAP of Rx6. Fig. 3: Measured path loss and composed power delay angular profile with antenna gain of Rx6in measurement set 1.problem, we adopt the RT technique to help understand the physical propagation of THz signalsin the study of THz channel measurement results. First, RT simulation is conducted in themeasurement environment, and then the simulated MPCs from the RT simulator are matched withthe MPCs extracted from the measurement. From the literature, one of the most common MPCsmatching methods is the threshold-based method [28]. If the multi-path component distance(MCD) between two MPCs, is a measure to quantify the distance between MPCs, is lower thana threshold, the two MPCs are regarded to be matched. Nevertheless, the effectiveness of theMCD matching method is sensitive to the pre-defined threshold.
Ray Tracing
SimulatorChannel Measurement
MPCs MCD-based
DBSCANSystem Calibration and Post processing
Measured
MPCs
Simulated MPCsMeasurement results Matched
Clusters
Outliers
Non-matched
ClustersEnvironment geometry
Fig. 4: Combined MPCs clustering and matching procedure.Instead, we propose a clustering and matching procedure that combines the MPCs clustering and MPCs matching with RT simulation, with reduced complexity and good clustering perfor-mance. The core of clustering and matching in our analysis is MCD-based DBSCAN, whichis a density-based clustering algorithm [29], for the following benefits. First, DBSCAN canautomatically select the number of clusters, without pre-defining. Second, DBSCAN can findan arbitrarily shaped cluster, which is not restricted to a certain shape like K-means. Third,DBSCAN is robust to outliers, which classifies the MPCs that do not belong to any clusteras noise. By contrast, K-means, GMM, KPM and KPD assign every MPC to a cluster, whichcauses over-clustering in practice. Fourth, DBSCAN is not sensitive to the initial cluster centers.On the downside, DBSCAN is not good at addressing the dataset where the clusters are closeto each other. Fortunately, the THz channel is sparse in both temporal and angular domains,which is suitable to invoke DBSCAN. Furthermore, DBSCAN requires a minimum number ofpoints in a cluster and measures the separation between two points the Euclidean distance. TheEuclidean distance is not suitable for measuring the distance in a circular domain (AOA in ourcase). Therefore, we replace the Euclidean distance in the original DBSCAN algorithm by MCD,which is given by, MCD = (cid:113) || MCD AoA || + ζ MCD τ , (2)where MCD AOA and MCD τ are the MCD of AoA and MCD of ToA, respectively. In addition, ζ is a suitable delay scaling factor. In particular, the MCD of AoA is obtained asMCD AoA = 12 | cos( θ )sin( θ ) − cos( θ )sin( θ ) | , (3)where θ and θ are the AoAs of two MPCs, respectively. MCD of ToA is calculated asMCD τ = τ std (∆ τ max ) | τ − τ | , (4)where τ and τ are the ToAs of two MPCs, respectively. Moreover, τ std = 20 ns denotes thestandard deviation of ToAs of all the MPCs, and ∆ τ max = max ( τ l ) − min ( τ l ) = 100 ns.The flow of the combined MPCs clustering and the matching procedure is depicted in Fig. 4.In order to filter out the noise in the measured channel as much as possible, we first filter outthe noise whose power is lower than a constant threshold of -140 dBm, and then set the minimalnumber of subpaths in a cluster as 6 in the DBSCAN algorithm. Therefore, if an MPC does nothave at least other 5 surrounding MPCs, it does not belong to a cluster and instead, is regardedas noise. The simulated MPCs from ray tracing and measured MPCs are clustered by the MCD-based DBSCAN algorithm. To be concrete, the output of MCD-based DBSCAN is divided into three parts, including matched clusters, non-matched clusters, and outliers. First, the matchedclusters are the clusters consisting of at least one simulated MPC. Second, the non-matchedclusters are the clusters consisting of no simulated MPC. The reasons for the appearance of non-matched clusters are twofold. On one hand, the environment geometry fed into the RT simulatoris simplified and differs from the real room. On the other hand, the types of traced MPCs dependon the configurations in the RT simulator. In our simulation, since we trace the reflected pathsfrom walls and the non-matched clusters are the paths that reflect from obstacles like chairs anddesks. Third, the outliers refer to the MPCs with negligible received power, or the non-matchedMPCs from the RT simulator, which are excluded from the channel analysis and modeling. ToA [ns] A o A [ ° ] (a) Clustering by the KPM algorithm.
21 3 4 5 67 8 (b) Clustering by the MCD-based DBSCAN procedure.
Fig. 5: Comparison between the K-Power-Means (KPM) algorithm and the proposed MPCsclustering and matching procedure for Rx1. A cross represents an MPC while a circle is noise.Clusters are marked with different colors. Subpaths in a cluster are with an identical color. Someof the matched clusters in (b) are circled with numeric labels.Fig. 5 illustrates the clustering outcomes of the KPM and the proposed MCD-based DBSCANalgorithms for Rx1 in the measurement set 1. In Fig. 5, markers with different colors denotedifferent clusters. A cross represents a measured MPCs while a circle represents a noise inFig. 5(b). The number of clusters is set to be the same for a fair comparison. We observethat KPM groups many MPCs with sparse separation into one cluster in Fig. 5(a), due to thedrawback that the KPM cannot find outliers. Moreover, KPM splits the LoS cluster into severalsmall clusters, although they are very close. This is due to the fact the number of clusters should meet the pre-determined parameter, and the ball-shape clusters are favored by the KPMalgorithm. By contrast, Fig. 5(b) shows that the proposed clustering and matching procedure withMCD-based DBSCAN algorithm outperforms the K-means algorithm. In particular, the matched-clusters are labeled with red circles, and their propagation in the environment is provided by theRT simulator. To be concrete, cluster 1 represents the LoS path with some subpaths caused bythe channel measurement platform. Cluster 2 denotes the paths reflected from the desks closeto wall 2 in Fig. 2b. Cluster 3 refers to the reflected path from the chairs around the desk inthe room center. Cluster 4 and cluster 7 incorporate the first-order reflected paths from wall 1and wall 2, respectively. Cluster 5 contains the second-order reflected paths traveling throughwall 1 and wall 4. Finally, the paths in cluster 6 undergo triple reflections on the wall. Byaccurately capturing the clustering effects and matching with the measured MPCs in the indoorenvironment, the proposed MCD-based DBSCAN algorithm performs well.III. TH Z C HANNEL C HARACTERIZATION AND A NALYSIS
Based on the processed measurement data, we shed light on the THz indoor channel charac-teristics, from the perspectives of path loss of the LoS path, MPCs reflection properties, and thetemporal and spatial multi-path features in the THz band. The properties and physical parametersrevealed in this section are useful as guidelines for THz communication system design.
Tx-Rx separation [m] P a t h l o ss [ d B ] MeasurementCI model (PLE=1.75)FSPL
Fig. 6: Comparison of measured path loss, the fitted CI model and the FSPL results at 140 GHz.
A. Path Loss Model
The path loss at each Rx in the measurement sets 1 and 2 is calculated by dividing the receivedpower by the transmitted power as depicted in Fig. 6, where the received power is the sum of the power of multi-paths above the noise floor. In particular, a close-in free space reference distance(CI) path loss model is developed based on the measured path loss, asPL CI [dB] = 10 × PLE × log ( dd ) + FSPL ( d ) + X σ , (5)where PLE is the path loss exponent, d denotes the distance between Tx and Rx, and d representsthe reference distance which is 1m in this work. X σ is a zero-mean Gaussian random variablewith standard deviation σ SF in dB, which represents the fluctuation caused by shadow fading.Moreover, we compute the free-space path loss (FSPL) by invoking the Friis’ law, given by,FSPL ( d , f ) = −
20 log ( c πf d ) , (6)where c denotes the speed of light, f represents the frequency. In addition, PLE in (5) isdetermined by minimizing σ SF via a minimum mean square error (MMSE) approach. Themeasurement results yield that the PLE value is 1.75 and σ D equals to 3.44 dB for the LoSchannel at 140 GHz in a meeting room.
100 110 120 130 140 150 160 170
Path loss [dB] C u m u l a ti v e p r ob a b ilit y Fig. 7: CDF of path loss of the reflected paths with different reflection orders at 140 GHz.
B. Reflection Loss
The cumulative density functions (CDFs) of the path loss of first-order, second-order, andthird-order reflected MPCs in our measurement campaign are computed in Fig. 7. By includingthe propagation and reflection attenuation, the average path loss values of the three kinds are124.2 dB, 129.5 dB, and 138.6 dB, respectively. These results suggest that one reflection leadsto additional loss between 6 dB and 9 dB, which agrees with the relative dielectric constant of drywall ( (cid:15) r = 6 . ) measured at 140 GHz [12]. Being consistent with intuition, a higher-order reflected path suffers from a higher path loss since it travels through a longer propagationdistance and bears excessive attenuation from reflection. Indeed, the reflection loss depends onthe incident angle, material parameters of objects, and roughness of the reflection surface [30]. C. Wall-reflection Paths
In the indoor environment, we further classify the reflected paths into wall-reflection pathsand obstacle-reflection paths. A wall-reflection path is one that interacts only with the wallsinside the meeting room. Otherwise, the reflection belongs to the obstacle-reflection category.The motivation of this classification is that in an extreme case where there are no obstacleslike furniture, only wall-reflection paths exist. Note that reflected MPCs from the ceiling andground are omitted as their elevation angles are out of the scope of our measurement equipment.On the other hand, the existence of furniture and other obstacles in a room can block somewall-reflection paths and produce multiple obstacle-reflection paths.To investigate the power distribution among the MPCs, we calculate the power ratio valuesbetween the wall-reflection clusters and the obstacle-reflection clusters, R w , which are listed inTable. II and Table. III, for the two measurement sets respectively. The power ratio values atall the receivers exceed one, which indicates that the energy carried by wall-reflection MPCs ishigher than that carried by obstacle-reflection MPCs, and consequently, the contribution fromwall-reflection rays is important to the received signal strength in the THz WLAN environment.The only exception appears at Rx2 in measurement set 2 since there exist some strong NLoSpaths that reflect over the desk and chairs near wall 2. We notice that in the measurement set1, Rx1, Rx2, and Rx4 which are on the same side of the desk as Tx receive multiple first-orderobstacle-reflection MPCs that come from the chairs around the desk. The ToAs of those MPCsare much shorter than the wall-reflection MPCs, which results in lower free-space propagationloss. Therefore, the values of R w at these receivers are small. By contrast, Rx3 and Rx6 are moredistant from the Tx, and the wall-reflection multi-paths undergo shorter propagation, which leadsto very large R w values. This observation is well aligned with the fact that the obstacle-reflectionMPCs are decided by the relative positions of Tx, Rx, and obstacles. In the measurement set 2, R w values at Rx2, Rx3, Rx4, Rx5, and Rx6 are much lower than others. The reason is that thescreen near wall 4 blocks the first-order wall-reflection path and cause another strong first-orderobstacle-reflection path for Rx2-6. TABLE II: Temporal and spatial characteristics for each Rx in measurement set 1. Rx d [m] N K
DS [ns] AS [ ◦ ] R w TABLE III: Temporal and spatial characteristics for each Rx in measurement set 2. Rx d [m] N K
DS [ns] AS [ ◦ ] R w D. Temporal and Spatial Features of THz Multi-path Propagation
The temporal and spatial distributions are important characteristics of multi-path propagation.In Fig. 5b, we clearly observe that the MPCs are resolvable in both temporal and angulardomains. Given in Table II and Table III, the number of clusters, N , is around ten in bothmeasurement sets, varying with the receiver position. Therefore, the indoor channel at 140 GHzpresents sparsity . Since the cost of the highly directional beam at Tx is that the LoS blockagemight cause a severe service interruption of the THz link, it is inappropriate to neglect multi-patheffects although the THz channel is sparse. According to our measurement results, the multi- path effect is indeed significant when the LoS path is blocked, where considerable energy iscontributed by the NLoS MPCs.To characterize the significance of the LoS path, we evaluate the K-factor , K , which is definedas the ratio between the power of the LoS cluster and the power of the remaining NLoS clusters,as computed in Table II and Table III for measurement sets 1 and set 2, respectively. A largeK-factor suggests that the LoS path dominates in the channel. With a good LoS condition in themeeting room, we compute that the K-factor value ranges from 8.59 to 56.43 in the measurementset 1, and ranges from 6.04 to 326.94 in the measurement set 2, which change greatly over variousreceiver positions in the room. The highest K-factor appears at those positions with shorter Tx-Rx separation, where the propagation distances of reflected MPCs are substantially longer thanthat of the LoS. Taking Rx1 as an example, the strongest reflected path is the one that reflectson wall 1. However, its ToA is 6 times of the Los ToA, and the incident angle is almost normalto the plane, which causes very high reflection loss. Thus, its received power is 25 dB lowerthan that of the LoS path. The second strongest NLoS paths are the first-order obstacle-reflectionpaths coming from the desk and chairs. Additional 6 dB free-space path loss is introduced sincetheir ToAs are double of the LoS ToA. Meanwhile, the angles-of-departure (AoDs) of the tworeflected MPCs are ◦ away from the LoS path, which leads to a noticeable beam misalignmentand introduces further attenuation in antenna gains. We observe that the received power of thetwo obstacle-reflection paths is at least 30 dB lower than that of the LoS path. To sum up, alonger propagation distance, high reflection loss, and beam misalignment of the reflected MPCsare the three critical factors that cause a high K-factor value.The power dispersion of MPCs in the temporal domain is quantified by root-mean-square(RMS) delay spread (DS) , τ rms , which is calculated as, τ rms = (cid:115) (cid:80) N τ i =0 ( i ∆ τ − ¯ τ ) A c ( i )∆ τ (cid:80) N τ i =0 A c ( i )∆ τ (7)with ¯ τ = (cid:80) N τ i =0 i ∆ τ A c ( i )∆ τ (cid:80) N τ i =0 A c ( i )∆ τ (8)where N τ is the number of sampled ToAs, and A c ( i ) = (cid:80) j (cid:80) k p r ( i, j, k ) denotes the powerdelay profile. Similarly, we compute angular spread (AS) to characterize the power dispersionof MPCs in the spatial domain. The DS and AS at the measured positions in the measurementset 1 and set 2 are listed in Table II and Table III, respectively. On the one hand, the delay and angular spreads differ among the receivers. However, the deviation is insignificant compared tothe metric of K-factor. A small DS value suggests that the received power is confined in thetemporal domain, while a small AS value indicates the received power comes from a narrowspatial region, which degrades the spatial degree of freedom.Even with very similar geometry properties, the temporal and spatial features could be distinct.For example, we notice that Rx6 and Rx7 are both at a corner of the meeting room and residesymmetrically about the diagonal of the meeting room. Besides, the distance difference from Tx1to the two receivers is only 0.38 m. However, Rx6 and Rx7 have substantially different temporaland spatial characteristics, in terms of K-factor, delay spread, and R w . For example, Rx7 hasthe largest delay spread of . ns and the smallest K-factor. By contrast, the smallest DSand the largest R w are measured at Rx6. The reason is that there exists a very strong first-orderwall-reflection MPC (i.e., reflecting on the wall 1 in Fig. 2b) that arrives at Rx6 shortly afterthe LoS path. Therefore, Rx6 has a high K-factor, large R w , as well as small DS. By contrast,at Rx7, the wall-reflection MPC (i.e., reflecting on wall 4), which is expected to be strong,is partially obstructed by a metal cabinet. Therefore, all the reflected MPCs have comparablepower, which is significantly weaker than the LoS path. As a result, one can observe a lowK-factor, large power dispersion in the time domain, and small R w at Rx7. In a nutshell, delayspread is severely influenced by the critical NLoS MPCs with low path loss and short ToA.Furthermore, we use log-normal distribution to fit the THz temporal and spatial features,including N , K , DS, AS and R w . The log-mean and log-standard-deviation values for thosechannel parameters in measurement set 1, set 2 and combined set 1&2 are computed in Table. IV.There is no noticeable difference in the mean and standard deviation values for the characteristicparameters between the measurement set 1 and set 2, which suggests that the two measurementsets are statistically consistent. E. Intra-cluster Characteristics
For intra-cluster characteristics, we analyze cluster delay spread (CDS) and cluster angularspread (CAS) for each cluster. The distributions of CDS and CAS in measurement set 1, set 2 areshown and compared in Fig. 8(a) and 8(b). The log-mean and log-standard-deviation values forCDS and CAS in measurement set 1, set 2 and combined set 1&2 are computed in Table. IV.We notice that the mean values of CDS and CAS are both smaller than DS and AS. To be TABLE IV: Parameters of log-normal distribution for the THz temporal and spatial features.
Measurement set1 2 1&2 µ ln N [ln( · )] 2.30 2.16 2.22 σ ln N [ln( · )] 0.40 0.26 0.33 µ ln K [ln( · )] 3.01 3.12 3.07 σ ln K [ln( · )] 0.63 1.00 0.86 µ ln DS [ln(ns)] 1.55 1.48 1.51 σ ln DS [ln(ns)] 0.48 0.39 0.43 µ ln AS [ln( ◦ )] 3.56 3.51 3.53 σ ln AS [ln( ◦ )] 0.26 0.43 0.37 µ ln Rw [ln( · )] 2.21 3.10 2.72 σ ln Rw [ln( · )] 1.66 1.64 1.71 µ ln CDS [ln(ns)] -1.01 -1.11 -1.06 σ ln CDS [ln(ns)] 1.18 0.95 1.07 µ ln CAS [ln( ◦ )] 1.77 1.79 1.78 σ ln CAS [ln( ◦ )] 0.54 0.50 0.52 C u m u l a ti v e p r ob a b ilit y Measurement set 1&2Measurement set 1Measurement set 2 (a) CDS. ° ]00.20.40.60.81 C u m u l a ti v e p r ob a b ilit y Measurement set 1&2Measurement set 1Measurement set 2 (b) CAS.
Fig. 8: CDS and CAS distribution of the measured channel at 140 GHz.concrete, the mean values of DS and AS are 4.23 ns and . ◦ , while the mean values of CDSand CAS are 0.35 ns and . ◦ . TABLE V: Correlation matrix among THz multi-path characteristics. d N K
DS AS R w d N -0.58 1 - - - - K -0.39 0.12 1 - - -DS 0.12 0.04 -0.35 1 - -AS 0.50 0.06 -0.63 0.44 1 - R w F. Correlation among THz Multi-path Characteristics
To obtain deep understanding and insights of the THz indoor channel properties, we proceed toevaluate the correlation among the temporal and spatial multi-path characteristics . In particular,a correlation matrix containing the separation distance among Tx and Rx, the number of clusters,K-factor, delay spread, angular spread, and the power ratio between wall-reflection and obstacle-reflection MPCs, R w , is presented in Table V. The element in the correlation matrix is thecorrelation coefficient between two channel characteristic parameters, which is calculated as, ρ X,Y = Cov ( X, Y ) (cid:112) Var [ X ] Var [ Y ] (9)where X and Y are the two sequences of measured channel characteristics, Cov ( X, Y ) denotesthe covariance of X and Y , and Var [ · ] is the variance. The value of the correlation coefficientranges from -1 to 1. In particular, ρ X,Y = 1 denotes a positively linear correlation, 0 denotesnon-correlation, and − denotes a negatively linear correlation.The main observations are drawn as follows. First, the distance is negatively correlated to thenumber of clusters and K-factor, which suggests that a long separation distance results in a higherchannel sparsity and a lower K-factor. This can be explained that if Tx and Rx separate far, thedifference in propagation distance between the LoS and NLoS MPCs decreases. Therefore, forthose receiver points far from Tx, the power of LoS is weakened compared with NLoS MPCs.Moreover, the number of clusters and K-factor demonstrate weak correlation, as their corre-lation coefficient is equal to 0.12. Furthermore, the number of clusters is neither correlated withdelay spread nor to angular spread. This suggests that the number of clusters has no impacton power distribution among the MPCs as well as power dispersion in the temporal and spatialdomains. Third, the correlation coefficient between delay spread and angular spread is 0.44, which implies that the temporal and spatial spreads are related. This result coincides with theintuition that the power is likely to disperse in the angular domain if it is dispersive in thetime domain. This is further aligned with our discussions about Rx6 and Rx7 in Sec. III-D,where the two receivers with almost the same separation distance have significantly differentcharacteristics. Fourth, K is negatively related to DS and AS, which is consistent with the factthat LoS power is strong with high K and in this case, the power is confined in the temporaland angular domains. Fifth, R w illustrates a positive correlation with the distance. Thus, if areceiver is far away from the transmitter, the wall-reflection power needs to be carefully treated,since it contributes significantly to the received power and dominates over other NLoS MPCs.IV. H YBRID C HANNEL M ODELING FOR TH Z I NDOOR P ROPAGATION
The temporal and spatial analysis of the channel measurement results suggest that the criticalMPCs, i.e., the LoS path and wall-reflection paths, dominate in the THz indoor channel. Inspiredby this, we develop a semi-deterministic or, more precisely, RT-statistical hybrid channel model.The proposed channel model is further validated with the measurement data. Moreover, theperformance of the proposed channel model is evaluated and compared with the conventionalstatistical model and 3GPP channel model [20].
A. Hybrid Channel Model
The RT-statistical hybrid channel model generates the channel by combining an RT and astatistical approach. First, we use the RT geometric-optic approach to produce the dominantMPCs, including one LoS path and several wall-reflection paths. Each of these MPCs is thecenter of a cluster with subpaths. In the RT method, only the dimensions of the room and thepositions of Tx and Rx are required. Second, a statistical approach supplements intra-clustersubpaths in the wall-reflection clusters as well as the additional obstacle-reflection clusters.Therefore, the spatial channel impulse response, h ( τ, θ, f ) , consists of the deterministic CIRgenerated by the RT approach, h RT ( τ, θ, f ) , and the statistical CIR, h s ( τ, θ, f ) , given by h ( τ, θ, f ) = h RT ( τ, θ, f ) + h s ( τ, θ, f ) . (10)where τ and θ denote the delay and azimuth AoA. We emphasize that the AoD is not consideredsince AoD is not scanned in the measurement campaign, as described in Sec. 2. Also, as elevationAoA is scanned from − ◦ to ◦ for each measurement point, which doesn’t occupy the whole elevation domain, we compose the CIRs in the elevation domain. Therefore, elevation AoA iseliminated in the model.
1) RT channel modeling:
The CIR generated by the RT approach, h RT ( τ, θ, f ) , is representedas, h RT ( τ, θ, f ) = A t ( (cid:126)φ LoS ) α LoS ( f ) δ ( τ − τ LoS ) δ ( θ − θ LoS )+ L RT (cid:88) l =0 A t ( (cid:126)φ l, ) α l, ( f ) δ ( τ − τ l, ) δ ( θ − θ l, ) , (11)where A t ( · ) represents the antenna pattern at Tx. α LoS ( f ) , τ LoS , θ LoS and (cid:126)φ
LoS denote theamplitude gain, ToA, azimuth AoA and AoD vector of the LoS path, respectively. L RT statesthe number of RT clusters. As we only consider the wall reflection with up to 3 reflection times,the total number of RT clusters is L RT = 20 (4 first-order wall-reflection paths, 8 second-orderwall-reflection paths, and 8 third-order wall-reflection paths). α l, ( f ) , τ l, , θ l, and (cid:126)φ l, denotethe amplitude gain, ToA, azimuth AoA and AoD vector of the central path in the l th RT cluster,respectively.In particular, the amplitude gain α l, is calculated as, α l, ( f ) = | Γ (1) l ( f ) || Γ (2) l ( f ) | · · · | Γ ( R l ) l | πf τ l, ( f ) (12)where Γ ( n ) l is the n th -order reflection coefficient of the l th RT cluster, R l is the reflection timesof the l th RT cluster, and f is the frequency. Γ l ( f ) is frequency-dependent and its detailedcalculation can be found in [30]. B. Statistical Modeling
In this subsection, we elaborate the statistical part of the RT-statistical hybrid model. In ourhybrid channel model, the subpaths of the RT clusters except the central path and all the non-RTclusters are statistically generated. The statistical CIR, h s ( τ, θ, f ) , is calculated as, h s ( τ, θ, f ) = L RT (cid:88) l =1 P l (cid:88) p = − Q l ,p (cid:54) =0 A t ( (cid:126)φ l,p ) α l,p δ ( τ − τ l,p ) δ ( θ − θ l,p ) , + L s (cid:88) q =1 S q (cid:88) s = − T q α q,s δ ( τ − τ q,s ) δ ( θ − θ q,s ) , (13)where P l denotes the number of pre-cursor intra-cluster subpaths in the l th RT cluster, while Q l stands for the number of the post-cursor intra-cluster subpaths in the l th RT cluster. α l,p , τ l,p , TABLE VI: Statistical parameters of the RT-statistical hybrid model.
Variable Model PrameterRT-cluster Non-RT clusterPre-cursor Post-cursor Pre-cursor Post-cursorNumber of clusters Given in (13) - -Number of subpaths ln( · ) ∼ N ( µ, σ ) µ = 2 . µ = 4 µ = 1 . µ = 2 . σ = 1 . σ = 1 . σ = 1 . σ = 1 . Inter-cluster time of arrival Given in (14) - λ inter = 13 . Intra-cluster time of arrival Given in (15) λ intra = 0 . λ intra = 0 . λ intra = 0 . λ intra = 0 . Inter-cluster angle of arrival Given in (16) µ v = − . µ v = − . κ v = 0 κ v = 0 Intra-cluster angle of arrival Given in (16) µ v = 0 µ v = 0 . µ v = 0 µ v = 0 κ v = 33 κ v = 3 . κ v = 16 κ v = 16 Inter-cluster amplitude Given in (17) a inter = 0 . a inter = 0 . b inter = − . b inter = − . Intra-cluster amplitude Given in (18) a intra = 0 . a intra = 0 . a intra = 0 . a intra = 0 . b intra = − . b intra = − . b intra = − . b intra = − . θ l,p and (cid:126)φ l,p denote the amplitude gain, ToA, azimuth AoA and AoD vector of the p th subpathin the l th RT cluster, respectively. In particular, the subpath with p = 0 refers to the central pathof the RT cluster, which has been characterized by RT. In addition, L s is the number of non-RTclusters. T q denotes the number of pre-cursor intra-cluster subpaths in the q th non-RT cluster,while S q describes the number of the post-cursor intra-cluster subpaths in the q th non-RT cluster. α q,s , τ q,s , θ q,s denote the amplitude gain, ToA, azimuth AoA and AoD vector of the s th subpathin the q th non-RT cluster, respectively.One should note that the antenna pattern of Tx is not involved in the statistical CIR in (13),which is different from the RT CIR in (11). The reason is that all the statistical parameters ofthe hybrid model are extracted from the channel measurement with a directional antenna at Tx.Therefore, the effect of antenna pattern at Tx is contained in the generated amplitudes. As anoverview, the main statistical parameters of the RT-statistical hybrid model are summarized inTable. IV, which are detailed as follows.
1) Number of clusters:
Fig. 9 shows that the number of non-RT clusters is linearly related tothe distance between Tx and Rx. In particular, the channel of far-separated Tx-Rx pair consists L s Measurementy=-1.056x+9.776
Fig. 9: Number of non-RT clusters versus the Tx-Rx separation.of fewer non-RT clusters. Therefore, the number of clusters can be fitted by a linear model as, L s = (cid:100)− . d + 9 . (cid:101) , (14)where (cid:100)·(cid:101) denotes the ceiling function to produce an integer number of clusters. P r ob a b ilit y d e n s it y f un c ti on PostNum dataLog-normal distribution
Fig. 10: Log-normal fitting of the number of post-cursor subpaths in non-RT clusters.
2) Number of subpaths:
The numbers of pre-cursor and post-cursor subpaths in both RT andnon-RT clusters follow a log-normal distribution. In Fig. 10, we present the CDF of the numberof post-cursor subpaths in a non-RT cluster and the log-normal fitting, for which the parametersare given in Table IV.
3) ToA:
The arrivals of inter-cluster MPCs and intra-cluster MPCs are modeled by Poissonprocesses. Therefore, the ToA between two adjacent inter-cluster MPCs, ∆ τ inter , is independentlyexponentially distributed, expressed as τ l +1 , = τ l, + ∆ τ inter , ∆ τ inter ∼ E ( λ inter ) . (15) where E ( · ) denotes the exponential distribution and λ inter is the inter-cluster rate of arrival.Similarly, the ToA between two adjacent intra-cluster subpaths is independently exponentiallydistributed as, τ l,p +1 = τ l,p + ∆ τ intra , ∆ τ intra ∼ E ( λ intra ) . (16)where λ intra is the intra-cluster rate of arrival.
4) AoA:
The azimuth AoA of inter-cluster MPCs and the azimuth AoA of the intra-clustersubpaths both follow a Von Mises distribution. The probability density function (PDF) of theVon Mises distribution is given by f v ( θ ) = exp( k v cos( θ − µ v )) / (2 πI ( κ v )) . (17)where the parameter µ v is a measure of location, κ v is a measure of concentration, and I ( · ) isthe modified Bessel function of order 0. It should be noted that µ v and κ v for AoA follow aVon Mises distribution in the unit of rad, as given in Table VI.
5) Amplitude:
The amplitude of the inter-cluster MPCs is calculated as α q, = α LoS a inter | τ l, − τ LoS | b inter . (18)where a inter and b inter are the coefficient and the exponent for the inter non-RT clusters.The amplitudes of the intra-cluster MPCs have the same representation, α l,p = α LoS a intra | τ l,p − τ l, | b intra . (19)where a intra and b intra are the coefficient and exponent of the intra cluster subpaths.In addition, the phase of each MPC is assumed to be independently uniformly distributed. C. Model Validation and Evaluation
We validate the developed RT-statistical hybrid channel model with measured data, and com-pare it with a conventional statistical model, and the 3GPP TR 38.901 model [20]. In the con-ventional statistical channel model, the clusters generated by RT are replaced by the statistically-generated ones in which the parameters are adopted from Table VI. The 3GPP TR 38.901 modelis GSCM, which first stochastically generates a virtual map with a certain number of scatteringclusters and then constructs the multi-paths based on a simplified RT technique. The GSCMchannel model is designed to characterize the angular channel characteristics, and address theproblems of smooth time evolution and spatial consistency of the conventional statistical channel Delay spread [ns] C u m u l a ti v e p r ob a b ilit y f un c ti on MeasurementRT-Statistical modelStatistical model3GPP GSCM (a) DS.
Angular spread [ ° ] C u m u l a ti v e p r ob a b ilit y f un c ti on MeasurementRT-Statistical modelStatistical model3GPP GSCM (b) AS.
Fig. 11: DS and AS validation of the proposed RT-statistical channel model and comparisonwith conventional statistical model and 3GPP GSCM.model. To generate the GSCM channel, the QuaDRiGa simulator [31] is utilized, which is anopen-source reference implementation of the baseline 3GPP TR 38.901 model. Specifically,indoor scenario in 3GPP TR 38.901 model [20] is chosen for the simulation. To feed thesimulator, the parameters of the 3GPP TR 38.901 model, including path loss exponent, themean and standard deviation of K-factor, the delay spread, and angular spread as well as thecorrelation matrix, are taken from our channel measurement for a fair comparison.The measured channel data for channel model validation and performance comparison isobtained from a separate channel measurement campaign, where Tx is placed at the position ofRx3 while Rx is placed at the position of Rx7-12 (i.e., 6 measured positions) in Fig. 2b.
1) Delay Spread and Angular Spread:
The delay spread and angular spread are widelyaccepted measures of small-scale fading of a channel model. For validation, Fig. 11(a) and 11(b)compare the results of the delay spread and angular spread based on the proposed hybridchannel model, the conventional statistical model, the 3GPP TR 38.901 model, and the channelmeasurement, respectively. Good agreement is observed, which suggests that the three channelmodels well characterize the small-scale fading of the THz indoor channel in the temporaldomain and angular domains, respectively.
2) Power-delay-angular Profile:
The power-delay-angular profile (PDAP) contains the jointtemporal-spatial characteristics of the multi-paths of THz channel. We regard that a good channelmodel should be consistent with the measured PDAP. The RMSE of the simulated PDAP is calculated by, RMSE = (cid:114) N τ N θ (cid:88) ( PDAP m ( i, j ) − PDAP s ( i, j )) (20)where PDAP m ( i, j ) denotes the measured power with time-of-arrival (ToA) of i ∆ τ and azimuthangle-of-arrival (AoA) of j ∆ θ . Similarly PDAP s ( i, j ) denotes the simulated power with thetime-of-arrival of i ∆ τ and azimuth angle-of-arrival of j ∆ θ . In addition, N τ is the number ofsampled ToAs while N θ is the number of sampled AoAs. As RMSE calculates the absolute errorof the simulated PDAP compared with the measured PDAP, a small RMSE value justifies thatthe channel model can well characterize the full channel response. The CDF of the PDAP withrespect to RMSE generated by those three channel models is shown in Fig. 12(a). The proposedRT-statistical model has the lowest RMSE, while the conventional statistical model incurs thehighest RMSE. The mean RMSE values of the proposed RT-statistical model, 3GPP GSCM, andconventional statistical model are 3.65 dB, 4.22 dB and 5.73 dB, respectively.The RMSE value accounts for the absolute error of all the points between two PDAPs.Therefore, we introduce a Structural Similarity Index Measure (SSIM), which is widely utilized inimage compression, image restoration, and pattern recognition as a measure of figure similarity,and has been repeatedly shown to significantly outperform MSE [32]. In particular, SSIM isa perception-based model that considers image degradation as perceived changes in structuralinformation variation [33]. The range of SSIM is from 0 to 1, which indicates two figures aretotally different (i.e., SSIM=0) to exactly the same (i.e., SSIM=1). Fig. 12(b) compares theCDF of PDAP SSIM of the different channel models. The mean SSIM values of the proposedRT-statistical model, conventional statistical model, and 3GPP GSCM are 0.51, 0.05, and 0.14,respectively. It can be observed that the proposed RT-statistical model has the highest SSIM,which is approximately 10 times higher than the conventional statistical model. Similar to themetric of RMSE, the 3GPP GSCM is a compromised model. As a result, we state that theproposed RT-statistical hybrid model accurately captures the power distribution in both temporaland angular domains, which performs substantially better than the conventional statistical andthe 3GPP GSCM channel models. D. Comparison and Discussion
From the numerical results in Sec. IV-C, we observe that the proposed RT-statistical hybridchannel model, conventional statistical channel model, and 3GPP GSCM can well characterizethe delay spread and angular spread of the measured THz indoor channel. However, this does not C u m u l a ti v e p r ob a b ilit y f un c ti on RT-Statistical modelStatistical model3GPP GSCM (a) RMSE of PDAP. C u m u l a ti v e p r ob a b ilit y f un c ti on RT-Statistical modelStatistical model3GPP GSCM (b) SSIM of PDAP.
Fig. 12: RMSE and SSIM of the proposed RT-statistical channel model, conventional statisticalmodel and 3GPP GSCM.suffice to validate the accuracy of a channel model. When it comes to PDAP, the conventionalstatistical channel model shows the unacceptable performance in terms of RMSE and SSIM. Thereason is that the conventional statistical channel model only takes the Tx-Rx separation intoaccount, without requiring geometry of the propagation environment as prior information. TheToA and AoA are independently and randomly generated and assigned to each MPC, accordingto the parameter statistics pre-defined in the statistical model. Therefore, the generated PDAP isnot realistic.By contrast, the GSCM is semi-deterministic, which considers the positions of Tx and Rx.Therefore, the ToA, AoA, and amplitude of the LoS path in the GSCM are consistent withthe measured LoS path. Nevertheless, the scattering environment is still randomly generatedaccording to the pre-defined statistics, before the simplified RT technique is implemented. Thisleads that the remainder of MPCs is noticeable inconsistent with the measured channel data.As a result, although the 3GPP GBSCM is more deterministic and more accurate than theconventional statistical in terms of PDAP, it still performs worse than our proposed RT-statisticalhybrid channel model clearly.The proposed RT-statistical hybrid channel model requires the positions of Tx and Rx, and thedimensions of the room for RT, which introduces additional computational complexity comparedwith the conventional statistical and the GSCM models. However, the complexity of RT is lowsince we only trace a small number of LoS and wall-reflection rays with up to triple reflections. Moreover, it should be emphasized that the complexity of the statistical part of our hybrid channelmodel is reduced compared to the conventional statistical model and GSCM counterparts, sincethe large-scale fading model, i.e., path loss and shadowing, are not necessary. This is due tothe fact that the power of LoS and wall-reflection multi-paths dominate in the channel, which iswell captured in the RT technique. As a result, both the RMSE and SSIM metrics with respectto PDAP of the proposed model outperform the existing models from the literature.V. C
ONCLUSION
In this paper, we have elaborated extensive channel measurements from 130 GHz to 140 GHzin an indoor meeting room. The directional antenna at Rx is scanned in both azimuth andelevation directions for resolving the multi-paths in the spatial domain. A novel MPC clusteringand matching procedure along with ray tracing techniques is utilized to post-process the receivedmulti-path signals. THz channel characterization and analysis suggest that the LoS and reflectedpaths from walls dominate in the THz channel of the meeting room, resulting in the measuredhigh K-factor and R w values.In the meeting room environment, we develop a hybrid cluster-based THz channel model, bycombining ray-tracing and statistical methods. The numerical results show that compared withthe statistical model and 3GPP GSCM, the hybrid channel model has better agreement with themeasured PDAPs for all Tx-Rx positions. This suggests that the deterministic part of the hybridchannel model introduced by ray tracing presents spatial consistency in nature and captures themost significant paths in the THz signal propagation. The statistical part of the hybrid channelmodel is necessary to complete the multi-path characteristics.R EFERENCES [1] I. F. Akyildiz, J. M. Jornet, and C. Han, “Terahertz band: Next frontier for wireless communications,”
PhysicalCommunication (Elsevier) Journal , vol. 12, pp. 16 – 32, Sep. 2014.[2] M. Shafi, J. Zhang, H. Tataria, A. F. Molisch, S. Sun, T. S. Rappaport, F. Tufvesson, S. Wu, and K. Kitao, “Microwavevs. millimeter-wave propagation channels: Key differences and impact on 5G cellular systems,”
IEEE CommunicationsMagazine , vol. 56, no. 12, pp. 14–20, Dec. 2018.[3] C.-X. Wang, J. Bian, J. Sun, W. Zhang, and M. Zhang, “A survey of 5G channel measurements and models,”
IEEECommunications Surveys & Tutorials , vol. 20, no. 4, pp. 3142–3168, Aug. 2018.[4] Y. Xing and T. S. Rappaport, “Propagation measurement system and approach at 140 GHz-moving to 6G and above 100GHz,” in Proc. of IEEE Global Communications Conference (GLOBECOM) , pp. 1–6, Dec. 2018. [5] K. Guan, B. Peng, D. He, J. M. Eckhardt, H. Yi, S. Rey, B. Ai, Z. Zhong, and T. K¨urner, “Channel Sounding and RayTracing for Intra-Wagon Scenario at mmWave and Sub-mmWave Bands,” IEEE Transactions on Antennas and Propagation ,Aug. 2020.[6] C. Han and Y. Chen, “Propagation modeling for wireless communications in the terahertz band,”
IEEE CommunicationsMagazine , vol. 56, no. 6, pp. 96–101, May 2018.[7] Y. Chen and C. Han, “Channel modeling and characterization for wireless networks-on-chip communications in themillimeter wave and terahertz bands,”
IEEE Transactions on Molecular, Biological and Multi-Scale Communications ,vol. 5, no. 1, pp. 30–43, Nov. 2019.[8] S. Priebe, C. Jastrow, M. Jacob, T. Kleine-Ostmann, T. Schrader, and T. Kurner, “Channel and propagation measurementsat 300 GHz,”
IEEE Transactions on Antennas and Propagation , vol. 59, no. 5, pp. 1688–1698, Mar. 2011.[9] S. Kim, W. T. Khan, A. Zaji´c, and J. Papapolymerou, “D-band channel measurements and characterization for indoorapplications,”
IEEE Transactions on Antennas and Propagation , vol. 63, no. 7, pp. 3198–3207, Apr. 2015.[10] S. Kim and A. Zaji´c, “Characterization of 300-GHz wireless channel on a computer motherboard,”
IEEE Transactions onAntennas and Propagation , vol. 64, no. 12, pp. 5411–5423, Oct. 2016.[11] J. M. Eckhardt, T. Doeker, S. Rey, and T. K¨urner, “Measurements in a real data centre at 300 GHz and recent results,” inProc. of 13th European Conference on Antennas and Propagation (EuCAP) , pp. 1–5, Mar. 2019.[12] Y. Xing, O. Kanhere, S. Ju, and T. S. Rappaport, “Indoor wireless channel properties at millimeter wave and sub-terahertzfrequencies,” in Proc. of IEEE Global Communications Conference (GLOBECOM) , pp. 1–6, Dec. 2019.[13] J. Fu, P. Juyal, and A. Zaji´c, “Modeling of 300 GHz chip-to-chip Wireless Channels in metal enclosures,”
IEEE Transactionson Wireless Communications , vol. 19, no. 5, pp. 3214–3227, Feb. 2020.[14] C.-L. Cheng, S. Sangodoyin, and A. Zaji´c, “Thz cluster-based modeling and propagation characterization in a data centerenvironment,”
IEEE Access , vol. 8, pp. 56 544–56 558, Mar. 2020.[15] K. Guan, G. Li, T. K¨urner, A. F. Molisch, B. Peng, R. He, B. Hui, J. Kim, and Z. Zhong, “On millimeter wave and THzmobile radio channel for smart rail mobility,”
IEEE Transactions on Vehicular Technology , vol. 66, no. 7, pp. 5658–5674,Nov. 2016.[16] N. A. Abbasi, A. Hariharan, A. M. Nair, A. S. Almaiman, F. B. Rottenberg, A. E. Willner, and A. F. Molisch, “Doubledirectional channel measurements for THz communications in an Urban environment,” in Proc. of IEEE InternationalConference on Communications (ICC) , pp. 1–6, Jun. 2020.[17] V. Petrov, J. M. Eckhardt, D. Moltchanov, Y. Koucheryavy, and T. Kurner, “Measurements of reflection and penetrationlosses in low terahertz band vehicular communications,” in Proc. of 14th European Conference on Antennas andPropagation (EuCAP) , pp. 1–5, Mar. 2020.[18] J. M. Eckhardt, T. Doeker, and T. K¨urner, “Indoor-to-Outdoor path loss measurements in an aircraft for terahertzcommunications,” in Proc. of IEEE 91st Vehicular Technology Conference (VTC2020-Spring) , pp. 1–5, May 2020.[19] Z. Yu, Y. Chen, G. Wang, W. Gao, and C. Han, “Wideband channel measurements and temporal-spatial analysis for terahertzindoor communications,” in Proc. of IEEE International Conference on Communications Workshops (ICC Workshops) , pp.1–6, Jun. 2020.[20] 3GPP, “TR 38.901: Study on channel model for frequencies from 0.5 to 100 GHz (Release 15),” 3GPP Recommendation,Tech. Rep., 2018.[21] X. Wu, C.-X. Wang, J. Sun, J. Huang, R. Feng, Y. Yang, and X. Ge, “60-ghz millimeter-wave channel measurements andmodeling for indoor office environments,”
IEEE Transactions on Antennas and Propagation , vol. 65, no. 4, pp. 1912–1924,2017. [22] C. Gustafson, K. Haneda, S. Wyne, and F. Tufvesson, “On mm-wave multipath clustering and channel modeling,” IEEETransactions on Antennas and Propagation , vol. 62, no. 3, pp. 1445–1455, Dec. 2013.[23] C. Ling, X. Yin, R. M¨uller, S. H¨afner, D. Dupleich, C. Schneider, J. Luo, H. Yan, and R. Thom¨a, “Double-directionaldual-polarimetric cluster-based characterization of 70–77 GHz indoor channels,”
IEEE Transactions on Antennas andPropagation , vol. 66, no. 2, pp. 857–870, Nov. 2017.[24] J. Li, B. Ai, R. He, M. Yang, Z. Zhong, Y. Hao, and G. Shi, “On 3d cluster-based channel modeling for large-scale arraycommunications,”
IEEE Transactions on Wireless Communications , vol. 18, no. 10, pp. 4902–4914, Jul. 2019.[25] M. Yang, B. Ai, R. He, L. Chen, X. Li, J. Li, B. Zhang, C. Huang, and Z. Zhong, “A cluster-based three-dimensionalchannel model for vehicle-to-vehicle communications,”
IEEE Transactions on Vehicular Technology , vol. 68, no. 6, pp.5208–5220, Apr. 2019.[26] N. Czink, P. Cera, J. Salo, E. Bonek, J.-P. Nuutinen, and J. Ylitalo, “A framework for automatic clustering of parametricmimo channel data including path powers,” in Proc. of IEEE Vehicular Technology Conference (VTC) , pp. 1–5, Sep. 2006.[27] R. He, Q. Li, B. Ai, Y. L.-A. Geng, A. F. Molisch, V. Kristem, Z. Zhong, and J. Yu, “A kernel-power-density-basedalgorithm for channel multipath components clustering,”
IEEE Transactions on Wireless Communications , vol. 16, no. 11,pp. 7138–7151, Aug. 2017.[28] N. Czink and C. Mecklenbrauker, “A novel automatic cluster tracking algorithm,” in Proc. of IEEE International Symposiumon Personal, Indoor and Mobile Radio Communications (PIMRC) , pp. 1–5, Sep. 2006.[29] M. Ester, H.-P. Kriegel, J. Sander, X. Xu et al. , “A density-based algorithm for discovering clusters in large spatial databaseswith noise,” in Proc. of the Second International Conference on Knowledge Discovery and Data Mining , vol. 96, no. 34,pp. 226–231, Aug. 1996.[30] C. Han, A. O. Bicen, and I. F. Akyildiz, “Multi-Ray channel modeling and wideband characterization for wirelesscommunications in the terahertz band,”
IEEE Transactions on Wireless Communications , vol. 14, no. 5, pp. 2402–2412,Dec. 2015.[31] S. Jaeckel, L. Raschkowski, K. B¨orner, and L. Thiele, “QuaDRiGa: A 3-D multi-cell channel model with time evolutionfor enabling virtual field trials,”
IEEE Transactions on Antennas and Propagation , vol. 62, no. 6, pp. 3242–3256, Mar.2014.[32] Z. Wang and Q. Li, “Information content weighting for perceptual image quality assessment,”
IEEE Transactions on ImageProcessing , vol. 20, no. 5, pp. 1185–1198, Nov. 2010.[33] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structuralsimilarity,”