Millimeter Wave and Sub-Terahertz Spatial Statistical Channel Model for an Indoor Office Building

Millimeter-wave (mmWave) and sub-Terahertz (THz) frequencies are expected to play a vital role in 6G wireless systems and beyond due to the vast available bandwidth of many tens of GHz. This paper presents an indoor 3-D spatial statistical channel model for mmWave and sub-THz frequencies based on extensive radio propagation measurements at 28 and 140 GHz conducted in an indoor office environment from 2014 to 2020. Omnidirectional and directional path loss models and channel statistics such as the number of time clusters, cluster delays, and cluster powers were derived from over 15,000 measured power delay profiles. The resulting channel statistics show that the number of time clusters follows a Poisson distribution and the number of subpaths within each cluster follows a composite exponential distribution for both LOS and NLOS environments at 28 and 140 GHz. This paper proposes a unified indoor statistical channel model for mmWave and sub-Terahertz frequencies following the mathematical framework of the previous outdoor NYUSIM channel models. A corresponding indoor channel simulator is developed, which can recreate 3-D omnidirectional, directional, and multiple input multiple output (MIMO) channels for arbitrary mmWave and sub-THz carrier frequency up to 150 GHz, signal bandwidth, and antenna beamwidth. The presented statistical channel model and simulator will guide future air-interface, beamforming, and transceiver designs for 6G and beyond.


I. INTRODUCTION
Mobile data traffic is increasing rapidly throughout the world and is predicted to reach 77 exabytes per month by 2022 [1].A large proportion of the data traffic increase comes from emerging indoor wireless applications such as 8K ultra high definition streaming, wireless cognition, and centimeter-level position location, which will be enabled by millimeter-wave (mmWave) and sub-Terahertz (THz) wireless systems due to the vast bandwidths in 6G and beyond [2], [3].
Severe outdoor-to-indoor (O2I) penetration loss of up to 60 dB at mmWave frequencies is beneficial for deploying isolated indoor mmWave systems from outdoor co-channel cellular systems [4].The 60 GHz band has been well studied in the literature [5]- [8] and used in the standards IEEE 802.11ad/ay for wireless local area network (WLAN) [9], [10].However, only a few indoor channel measurements and modeling works at other emerging frequencies or across a vast swath of spectra, such as 28, 73, and 142 GHz, have been published [11], [12].Accurate channel models over mmWave and sub-THz frequencies are needed for the design and evaluation of 6G wireless communications and beyond [2].
MmWave and THz (i.e., 30 GHz -3 THz) have distinct propagation characteristics from sub-6 GHz [13].MmWaves do not diffract well and become more sensitive to the dynamic blockage by humans due to the short wavelength [14], [15].Directional, steerable high gain antennas with beamforming techniques are required to compensate for additional path loss within the first meter of propagation distance as the carrier frequency increases [16].Thus, time-variant directional channel models are vital for efficient beam tracking and selection algorithms and proper system design and deployment guidelines.
The remainder of the paper is organized as follows.Section II provides a brief review of the existing works on channel modeling in indoor environments at mmWave and THz frequencies.Section III describes the 28 and 142 GHz measurement systems used in this work and the indoor office environment, as well as the step-by-step measurement procedure.Section IV presents the directional and omnidirectional path loss data and resulting models, showing that similar path loss exponents were observed at 28 and 142 GHz in the NLOS environment.Section V introduces the 3-D spatial statistical channel impulse response (CIR) model for indoor office scenarios, and Section VI provides empirical statistics and distribution fitting of channel parameters derived from the 28 and 140 GHz measurement datasets in both lineof-sight (LOS) and non-line-of-sight (NLOS) environments.Simulated secondary channel statistics (i.e., root mean square (RMS) delay spread (DS) and RMS angular spread (AS)) are generated from the NYUSIM indoor channel simulator and compared with the measured values, which yield good agreements.Finally, concluding marks in Section VIII show that the number of time clusters follows a Poisson distribution and the number of subpaths within each cluster follows a composite exponential distribution for both LOS and NLOS environments at 28 and 140 GHz, but the total number of observed subpaths at 140 GHz is much fewer than the number at 28 GHz.

II. EXISTING WORKS ON INDOOR CHANNEL MODELS
Numerous indoor channel measurements and studies have been focused on sub-6 GHz and 60 GHz [6], [8], [17]- [28].Saleh and Valenzuela conducted propagation measurements in an office building using radar-like pulses with 10 ns width at 1.5 GHz, and observed that multipath components (MPCs) arrived in clusters [17].A cluster-based statistical channel model was proposed, where the cluster arrival time and the subpath arrival time within each cluster were Poisson distributed, and the expected cluster power and subpath power were modeled as exponentially decaying functions of cluster arrival time and subpath arrival time within each cluster, respectively.This modeling approach has been extensively used in the past few decades.Rappaport conducted propagation measurements at 1.3 GHz in factories and showed that MPCs arrived independently rather than in clusters for factory and open plan building which contain reflecting objects spread throughout the workspace [19].Other indoor channel measurements and modeling efforts at mmWave frequencies started from the early 1990s, a majority of which were conducted at 60 GHz [6], [8], [20]- [28].
Standard documents such as IEEE 802.11 ad/ay and 3GPP TR 38.901 presented statistical channel models up to 100 GHz for indoor scenarios such as home, office, shopping mall, and factory [9], [10], [29].IEEE 802.11 ad/ay channel models adopted a double-directional CIR model for 60 GHz with dual polarizations based on field measurements and complimentary ray-tracing simulations, which provided detailed temporal and angular channel statistics for conference room, cubical environment, and living room [9], [10].3GPP TR 38.901 proposed a unified geometry-based statistical channel model for indoor and outdoor scenarios for frequencies from 0.5 to 100 GHz, where different scenarios have different values of large-scale parameters (i.e., DS, AS, Rician K factor, and shadow fading) which are required in the channel generation procedure [29].
THz communication systems will most likely be deployed in indoor environments to support extremely high data rates of over 100 Gbps [2].Considering the particular characteristics of THz signals such as high free space path loss (FSPL), high partition loss, dynamic shadowing loss due to human and vehicle blockages, deep understanding of the THz radio propagation channels is critical for 6G and beyond [30]- [36].First of all, atmospheric or molecular absorption induces non-negligible path loss at THz frequencies, especially at several absorption peaks such as 170 and 325 GHz, which causes about 100 dB/km attenuation [37].Thus, a frequencydependent atmospheric absorption term e −k(f )d was introduced in Friis formula [30], [31], where k(f ) is the atmospheric attenuation factor, f and d are the carrier frequency and transmission distance, respectively.A temporal-spatial stochastic channel model for 275-325 GHz was established based on ray-tracing channel simulations having LOS, first- and second-order reflected paths in an office room [32].This stochastic channel model can generate channel transfer function, power delay profile (PDP), and angular power spectrum (APS).The adopted ray tracer was calibrated using vector network analyzer (VNA)-based measurements at 275-325 GHz in an office [32].A generic multi-ray CIR model based on ray tracing consisted of LOS, reflected, diffracted, and scattered paths was proposed and used in channel capacity analysis [33].The reflection, diffraction, and scattering coefficients used in this multi-ray channel model were calibrated by measurements conducted at 0.06-1 THz.Most of the existing channel models for THz frequencies were built upon free space, reflection, and scattering measurements for various materials and constructed as a superposition of LOS, reflected and scattered paths in a ray tracing manner.Most propagation measurements were short-distance within a few meters and confined to a single room [30]- [32].
This paper derives empirical channel statistics based on extensive radio propagation measurements at 28 and 140 GHz conducted on the entire floor of an office building, and proposes a 3GPP-like indoor spatial statistical channel model following the mathematical framework of the NYUSIM outdoor channel models [38], which can generate directional and omnidirectional wideband CIRs from 28 to 140 GHz.

III. 28 GHZ AND 140 GHZ WIDEBAND INDOOR CHANNEL MEASUREMENTS
The 28 and 140 GHz measurement campaigns were conducted in the identical environment, NYU WIRELESS research center on the 9th floor of 2 MetroTech Center in downtown Brooklyn, New York at 2014 and 2019.A wideband sliding correlation-based channel sounder system was used in both measurement campaigns, providing a broad dynamic range of measurable path loss (152 dB at 28 GHz and 145 dB at 140 GHz) [11], [39].A wideband pseudorandom noise (PN) sequence of length 2047 was generated at baseband, then upconverted to a center frequency of 28 and 142 GHz, and transmitted through a directional and steerable horn antenna at the transmitter (TX).The receiver (RX) captured the RF signal via an identical steerable horn antenna and downconverted and demodulated the RF signal into its baseband I and Q signals [40].The demodulated signal was then correlated with a local copy of the transmitted signal with a slightly lower rate, which allowed the received signal to "slide" past the slower sequence [40].An average PDP over 20 instantaneous PDPs was sampled by a high-speed oscilloscope and recorded for further analysis.TX and RX antennas were mechanically steered by two electrically-controlled gimbals with sub-degree accuracy in azimuth and elevation planes and were switched between vertical-and horizontal-polarization modes by a 90degree waveguide twist for co-and cross-polarization studies.The 28 and 140 GHz channel sounder specifications are summarized in Table II.Null-to-null RF bandwidths of 800 MHz and 1 GHz were adopted in the 28 and 140 GHz measurement campaigns, resulting in a time resolution of MPC equal to 2.5 ns and 2 ns, respectively [11], [39].
Omnidirectional channel statistics are often preferred in channel models and channel simulations since arbitrary antenna patterns can be added [38].Thus, omnidirectional PDPs should be recovered from measured directional PDPs by aligning these measured PDPs with absolute time delays [38], [41], [42].However, the 28 GHz channel sounder did not have precise synchronization between TX and RX and cannot provide absolute timing information of measured PDPs since the PDP recording was triggered at the first MPC arrival and only had excess time delay information.A 3-D ray tracer NYURay was employed to provide the time of flight (i.e., absolute time delay) of the first arriving MPC in a measured PDP, which will be explained in Section III-B.The 140 GHz channel sounder was equipped with rubidium standard references at both TX and RX sides for frequency/timing synchronization [40]; however, we used the absolute time delays obtained from the 3-D ray tracer to synthesize omnidirectional PDPs for both 28 and 140 GHz data for processing consistency.

A. Measurement Environment and Procedure
The measurements were conducted in a typical indoor office environment (65.5 m × 35 m × 2.7 m) with offices, conference rooms, classrooms, long hallways, open-plan cubicles, and elevators, as shown in Fig. 1.Common obstructions are desks, chairs, cubicle partitions, glass doors, and walls made of drywall with metal studs.
Five TX locations and 33 RX locations were selected in the 28 GHz measurement campaign in 2014.Overall, measurements were conducted at nine LOS location pairs (i.e., a pair of TX and RX locations) and 35 NLOS location pairs, where the 3-D TX-RX (T-R) separation distances ranged from 3.9 m to 45.9 m.The identical five TX locations and a subset of the RX locations were measured in the 140 GHz measurement campaign due to the limit of maximum transmit power in 2019 and 2020, resulting in nine LOS location pairs and 13 NLOS location pairs.The T-R separation distance ranged from 3.9 m to 39.2 m.In Fig. 1, TX and RX locations measured at both 28 GHz and 140 GHz are denoted as stars and circles  [11], [43].was measured with TX2, TX3, and TX4).
For each T-R location pair, eight unique antenna azimuth sweeps were measured to investigate the spatial statistics of arrival and departure, where six RX antenna azimuth sweeps and two TX antenna azimuth sweeps were performed.During each sweep, TX (RX) horn antenna was rotated in step increments of the antenna half-power beamwidth (HPBW) [11] so that the directional measurements can emulate channel measurements using omnidirectional antennas.The detailed description of each measurement sweep is listed in Table III.
The equivalent omnidirectional received power can be synthesized by summing the received powers from all measured unique pointing angles obtained at antenna HPBW step increments in both planes [41].The sweeping step was equal to the antenna HPBW (30°for 28 GHz and 8°for 140 GHz), which corresponded to 12 and 45 rotation steps over the complete azimuth plane, respectively.At each rotation step, an averaged PDP over 20 instantaneous PDPs with accurate excess timing information with time resolutions of 2.5 and 2 ns was recorded for 28 GHz and 140 GHz, respectively.Note that two antenna polarization configurations, vertical-tovertical (V-V) and vertical-to-horizontal (V-H), were measured using the identical procedure described above, resulting in 16 measurement sweeps in total at each unique T-R location pair.This paper mainly focuses on the co-polarized (V-V) polarization to develop the omnidirectional and directional indoor channel models.For each T-R location pair, at most 96 (= 8×12) at 28 GHz and 360 (= 8×45) at 140 GHz directional PDPs were acquired with V-V polarization configuration.[11].TX and RX locations measured at both 28 GHz and 140 GHz are denoted as stars and circles with checkerboard texture, respectively.RX locations only measured at 28 GHz are denoted as solid circles.Each of the five TX locations is denoted in a different color, and the RX locations paired with a TX location is denoted in the same color.

B. Synthesizing Omnidirectional PDPs
A 3-D mmWave ray-tracing software, NYURay [44], was used to predict possible propagating rays between the TX and RX and provided the time of flight (i.e., absolute time delay) of the first arriving MPC of a measured directional PDP.Since the horn antennas had beamwidths of 30°and 8°for 28 GHz and 140 GHz, the exact angle of departure and angle of arrival of MPCs were unknown.Each measured directional PDP was assigned to a predicted ray which was closest to this directional PDP in space, then the absolute time delay of the first arriving MPC of the PDP was set to be the time of flight of the corresponding predicted ray.Directional PDPs were aligned in the temporal domain and summed to generate an omnidirectional PDP.Being closest in space between a measured PDP and the set of predicted rays means the antenna gain in the direction of a predicted ray when the antenna is pointing to the direction of the measured PDP is the highest among all predicted rays: (1) where P is the set of predicted rays from NYURay.∆φ AOD , ∆θ ZOD , ∆φ AOA , ∆θ ZOA denote the absolute difference of the azimuth angle of departure (AOD), zenith angle of departure (ZOD), azimuth angle of arrival (AOA), zenith angle of arrival (ZOA) between the measured directional PDP and a predicted ray, respectively.G φ and G θ represent the antenna pattern in the azimuth and elevation planes, respectively.The gain at the peak of the antenna main lobe is normalized to 0 dB.Thus, the antenna gain is -3 dB when the angle difference is 1/2 HPBW from the peak of the main lobe.
MPCs recorded in different directional PDPs may have originated from the same predicted ray, in which case the measured MPC was an antenna-gain weighted version of the true MPC.Thus, the directional PDPs assigned to the same predicted ray were summed in powers and generate a partialomnidirectional PDP for MPC extraction to avoid double counting.The direction of the extracted MPC was assumed to be the direction of the measured directional PDP which was closest to the predicted ray in space.An omnidirectional PDP was recovered, and MPCs were extracted by applying this procedure to all directional PDPs measured at each T-R location pair.Due to the mismatching between the measured PDPs and predicted rays at a few locations, we recovered omnidirectional PDPs for 37 of 44 location pairs at 28 GHz and 20 of 22 location pairs at 140 GHz.

A. Directional Path Loss Modeling
Antenna arrays with many elements will enable steerable and narrow beams to compensate for the large free space path loss in the first meter at mmWave and sub-THz frequencies.Directional path loss modeling is increasingly critical for future 6G and beyond communication system design.Thus, in this work, rotatable high-gain horn antennas were used at both the TX and RX during the 28 GHz and 140 GHz measurements, as shown in Table II, to study double-directional channels.
We use the close-in free space reference distance (CI) path loss model with 1 m reference distance [13], as this has The TX and RX antennas were pointed directly towards each other on boresight in both the azimuth and elevation planes (for LOS or NLOS environments).The RX antenna was then swept in the azimuth plane in steps of HPBW, for a fixed TX antenna at the boresight azimuth and elevation angles.

RX sweep
With respect to the boresight angle in elevation, the RX antenna was uptilted by HPBW and then swept in the azimuth plane in steps of HPBW, for a fixed TX antenna at the boresight azimuth and elevation angles.

RX sweep
With respect to the boresight angle in elevation, the RX antenna was downtilted by HPBW and then swept in the azimuth plane in steps of HPBW, for a fixed TX antenna at the boresight azimuth and elevation angles.

RX sweep
With respect to the boresight angle in elevation, the TX antenna was uptilted by HPBW.The RX antenna was fixed at the boresight elevation angle, and then swept in the azimuth plane in steps of HPBW.

RX sweep
With respect to the boresight angle in elevation, the TX antenna was downtilted by HPBW.The RX antenna was fixed at the boresight elevation angle, and then swept in the azimuth plane in steps of HPBW.

TX sweep
The TX and RX antennas were pointed directly towards each other on boresight in both the azimuth and elevation planes.The TX antenna was then swept in the azimuth plane in steps of HPBW, for a fixed RX antenna at the boresight azimuth and elevation angles.

RX sweep
This measurement was an RX sweep with the TX antenna set to the second strongest AOD in the azimuth and elevation plane.The second strongest AOD was determined by comparing the signal level from all the AODs during Measurement 6, except for the angles corresponding to the main angle of arrival.The RX antenna was fixed at the boresight elevation angle and then swept in steps of HPBW in the azimuth plane.

TX sweep
This measurement corresponds to the second TX sweep with TX antenna either uptilted or downtilted by HPBW after determining the elevation plane with the strongest received power from Measurement 4 and Measurement 5 during measurements.The RX antenna was pointed towards the initial boresight azimuth and elevation angles, and the TX was uptilted or downtilted by HPBW, and then swept in steps of HPBW in the azimuth plane.
been proven to be superior for modeling path loss over many environments and frequencies [45].PL CI represents the path loss in dB scale, which is a function of distance and frequency: where n denotes the path loss exponent (PLE), and χ σ is the shadow fading (SF) that is commonly modeled as a lognormal random variable with zero mean and σ standard deviation in dB.d is the 3-D T-R separation distance.d 0 is the reference distance, and FSPL(f, d 0 ) = 20 log 10 (4πd 0 c/f ).The CI path loss model uses the FSPL at d 0 = 1 m as an anchor point and fits the measured path loss data with a straight line controlled by a single parameter n (PLE) obtained via the minimum mean square error (MMSE) method.Throughout this paper, LOS and NLOS locations are defined according to whether the TX and RX can see each other.Here, for the directional path loss modeling, we define the LOS direction for LOS locations as the direction when the TX and RX directional antennas are pointed directly to each other.The LOS direction can be calculated based on the relative position of the TX and RX -the LOS direction is along the line of bearing between the TX and RX.The NLOS-Best direction is only defined for NLOS locations, and represents the best pointing direction for which minimum path loss is measured, which can be found by thoroughly rotating TX and RX directional antennas in the 3D space.The NLOS direction is defined for both LOS and NLOS locations, and represent all pointing directions which received detectable powers other than the LOS and NLOS-Best directions [11], [46].Fig. 2 shows the directional CI path loss model using measured path loss data at 28 GHz and 142 GHz [11].The path loss in the LOS direction is represented by a green circle for the LOS locations, and the path loss in the NLOS-Best direction is represented by a blue diamond for the NLOS locations.Measurements pointed to other directions are denoted by red crosses as NLOS directions for both LOS and NLOS locations.Comparing Fig. 2a and Fig. 2b, the LOS PLEs at 28 and 142 GHz are 1.7 and 2.1 which may be due to the differences in antenna HPBWs (30°and 8°).Wider beamwidths may capture more energy through reflection and scattering in the vicinity of the LOS direction, causing a PLE of less than 2. Furthermore, the NLOS-Best PLE at 28 and 142 GHz are about 3.0, suggesting strong NLOS paths are available to provide a sufficient link margin and can be leveraged by intelligent reflecting surfaces [47].

B. Omnidirectional Path Loss Modeling
Even though the directional path loss model will be widely used in future wireless system deployment, the omnidirectional path loss model is fundamental and serves as a reference model in various standard documents [10], [29].In Fig. 3, we present the omnidirectional path loss data and the fitted CI path loss model.The omnidirectional path loss is synthesized from received powers from all directions measured in the 3-D space [41].Fig. 3 marks LOS and NLOS scenarios in green and blue, respectively.The LOS PLEs at both frequencies are lower than 2.0, where 28 GHz shows a surprisingly low PLE of 1.2, which can be attributed to the waveguide effect in some corridor measurement locations.Note that both 28 and 142 GHz have a comparable PLE of about 2.7 in the NLOS environment, indicating that the signal power drops equally versus distances after the first meter in the mmWave band of 28 GHz and the sub-THz band of 140 GHz [11], [16].

V. 3-D SPATIAL STATISTICAL CHANNEL MODEL
A received signal can be viewed as a superposition of multiple replicas of the transmitted signal with different delays and angles for any wireless propagation channel [48].An extended S-V channel model [17] was commonly used to represent the double directional channel in the 3-D space [8], [29].MPCs were observed to arrive in clusters in delay and angular domains from 28 GHz and 140 GHz indoor channel measurements, which agreed with many early works [8], [17].Current standards document such as the 3GPP TR 38.901 channel model defined a cluster as a group of MPCs closely spaced in the joint temporal-spatial domain, where each cluster represented a reflector or a scatterer in the environment [8], [9], [29].
We observed in the measurements that MPCs traveling close in time may arrive from very different directions due to the symmetric structure of the environment like hallways [13], [40].Conversely, MPCs arriving from a similar direction may have very different propagation times.A time cluster  spatial lobe (TCSL) approach was introduced to characterize temporal and angular domains separately [38].A time cluster (TC) comprises MPCs traveling close in time and arriving from potentially different directions.A spatial lobe (SL) represents a main direction of arrival or departure where MPCs can arrive over hundreds of nanoseconds [38].
Both modeling methodologies are valid, where the 3GPP model is more widely used and the NYUSIM model using TCSL has a more straightforward and physically-based structure [49], [50].Performance evaluation with respect to spectrum efficiency, coverage, and hardware/signal processing requirements between the 3GPP and NYUSIM channel models were provided in [51].
The cluster-based omnidirectional CIR h omni (t, where t is the absolute propagation time, − → Θ = (φ AOD , θ ZOD ) is the AOD vector, and − → Φ = (φ AOA , θ ZOA ) is the AOA vector.N and M n denote the number of TCs and the number of subpaths within each TC, respectively.For the mth subpath in the nth TC, a m,n , ϕ m,n , τ m,n , − −− → Θ m,n , and − −− → Φ m,n represent the magnitude, phase, absolute time delay, AOD vector and AOA vector, respectively.Note that MPC and subpath are used interchangeably.The PDP and APS can be obtained by integrating the square of the CIR in space and time domains, respectively.
The PDP and APS can be easily partitioned based on TCs and SLs, respectively (see Fig. 10 and Fig. 11 in [38]).The partition in the time domain is realized by defining a minimum inter-cluster time void interval (MTI).Two sequentially recorded MPCs belong to two distinct TCs if the difference of the excess time delays of these two MPCs is beyond MTI.These two MPCs are considered as the last MPC of the former TC and the first MPC of the latter TC, respectively.For example, 25 ns was used as MTI for an outdoor urban microcell (UMi) environment [38], while 6 ns is used as MTI in this paper for an indoor office (InO) environment since the width of a typical hallway in the measured indoor office environment is about 1.8 m (i.e., ∼6 ns propagation delay).
The partition in the space domain is realized by defining a spatial lobe threshold (SLT) [38].The angular resolution of the measured APS depends on the antenna HPBW (30°and 8°for 28 GHz and 140 GHz, respectively).A linear interpolation of the directional received powers in azimuth and elevation planes with 1°resolution was used to reconstruct the 3-D spatial distribution of the received power.A power segment is generated for every 1°direction in the 3-D space.Neighboring power segments above the SLT form an SL.The SLT was -15 dB below the maximum directional power in the APS.
As defined in [38], the primary statistics such as the number of TCs and SLs, cluster delays, and cluster powers are used in the channel generation procedure given in Section VI.The secondary statistics such as RMS DS and RMS AS are not required in the channel generation but necessary in the channel validation.The presented channel model will be validated in Section VII by showing that the simulated and measured secondary statistics yield good agreements.

VI. STATISTICS OF CHANNEL GENERATION PARAMETERS
As described in Section V, temporal and spatial channel parameters are extracted from the measured PDP and APS.Temporal parameters are the number of TCs (N ) and SPs in a TC (M n ), TC excess delay (τ n ) and intra-cluster subpath excess delay (ρ m,n ), TC power (P n ) and subpath power (Π m,n ).Spatial parameters are the number of SLs (L), the mean azimuth and elevation angle of an SL (φ and θ), and the azimuth and elevation angular offset of a subpath (∆φ and ∆θ) with respect to the mean angle of the SL.
Since the 140 GHz measurement locations is a subset of the 28 GHz measurement locations, the 28 GHz common set was created out of the 28 GHz all set to have a fair comparison with the 140 GHz dataset (referred to the 140 GHz common set below).The 28 GHz common set and the  140 GHz common set have identical TX-RX location pairs.Table IV presents the channel parameters required for channel generation procedure.Table V provides statistics of channel parameters derived from the 28 GHz all set, 28 GHz common set, and 140 GHz common set, for LOS and NLOS scenarios.Step 2 # Cluster subpaths Mn Mn ∼ DU(1, Ms) Mn ∼ (1 − β)δ(Mn) + DE(µs) Step 3 Cluster delay τn (ns) Step 4 Intra-cluster delay ρm,n (ns) 1+Xn , m = 1, 2, ..., Mn, n = 1, 2, ..., N ρm,n ∼ Exp(µρ) Step 5 Cluster power Pn (mW) Zn ∼ N (0, σ Z ), n = 1, 2, ..., N Step 6 Subpath power Πm,n(mW) , Um,n ∼ N (0, σ U ), m = 1, 2, ..., Mn Step 7 Subpath phase ϕ (rad) Uniform(0, 2π) Step 9 Spatial lobe mean angle φ i , θ i (°) Step 10 Subpath angle offset ∆φ i , number of TCs, which is given by: The 28 GHz channel has about three more TCs than the 140 GHz channel in both NLOS and LOS scenarios, which can be attributed to the higher partition loss at 140 GHz (e.g., 4-8 dB higher than 28 GHz for different materials [43]).The channel sparsity at 140 GHz should be considered in the channel estimation and beamforming algorithms for sub-THz frequencies.The NLOS scenario has about one more TC than the LOS scenario.Note that the Poisson distribution of the number of TCs for the indoor scenario is different from the uniform distribution used for the outdoor scenario, as given in Table IV [38].
2) Number of Cluster Subpaths: The number of cluster subpaths M n is negatively correlated to the number of TCs depending on the MTI.A larger MTI causes fewer TCs and more subpaths within each TC, and vice versa.Fig. 5 presents the empirical histograms for three datasets, indicating that the number of cluster subpaths is close to exponentially distributed.The exponential distribution is continuous and starts from zero, while the value of the number of subpaths is discrete and starts from one.Thus, a discrete exponential (DE) distribution is applied to fit the empirical histogram of M n = M n − 1. Fig. 5a shows that about half of the measured TCs only have one subpath at 28 GHz, making a simple DE distribution unsuitable.We proposed a composite distribution with a δ-function at M n = 0 and a DE distribution, which is given by where µ s is the mean of the DE distribution, and β is the weight of the DE distribution in the composite distribution.By maximizing the joint probability mass function (PMF) of all data samples over β and µ s simultaneously, the MLE of µ s and β is 5.3 and 0.7 for 28 GHz NLOS all set, respectively.The identical composite distribution for 28 GHz LOS all set shows that µ s = 3.7, β = 0.7, suggesting that the NLOS scenario forms larger clusters than the LOS scenario.The composite distribution yields a good agreement with empirical histograms of 28 GHz LOS and NLOS scenarios.Moreover, the large TCs with more than 25 subpaths were mainly from locations in the corridor environment (e.g., TX4 and RX15).for both LOS and NLOS scenarios is about 1, suggesting that LOS and NLOS scenarios have similar sizes of clusters which contain about two subpaths on average.The clusters at 140 GHz are much smaller than the clusters at 28 GHz, indicating that the 140 GHz channel is much sparser than the channel at 28 GHz.The detailed comparison of channel parameters between 28 GHz and 140 GHz common sets can be found in Table V.
3) Inter-cluster Excess Delay: The cluster excess delay τ n is defined as the time difference between the first arriving subpath in the PDP and the first arriving subpath in a cluster, as given in Table IV, where τ n−1 is the cluster excess delay of the former cluster, ρ Mn−1,n−1 is the intra-cluster excess delay of the last subpath in the former cluster.∆τ n is the intercluster excess delay without MTI (i.e., 6 ns).The empirical cumulative distribution function (CDF) of the inter-cluster delay for LOS and NLOS scenarios of 28 GHz all set are shown in Fig. 6.An exponential distribution with the mean 10.9 ns fits the 28 GHz NLOS scenario, while a lognormal distribution with the mean 2.1 ns and standard deviation 1.6 ns fits the 28 GHz LOS scenario well since a few clusters with long cluster delays were observed in the LOS corridor environment.
Inter-cluster delays at 140 GHz can be well fitted using an exponential distribution for both LOS and NLOS scenarios, where the mean values are 14.6 ns and 21.0 ns, respectively.The distributions for the 28 GHz and 140 GHz LOS scenarios are different (lognormal and exponential), likely due to the higher partition loss and the smaller measurable range at 140 GHz.In addition, clusters with large inter-cluster delays were mainly observed in the corridor environment due to the waveguide effect, indicating that the corridor scenario may be considered a distinct indoor scenario and requires more channel measurements for accurate characterization.
4) Intra-cluster Excess Delay: The intra-cluster excess delay is defined as the time difference between the first arriving subpath and the targeted arriving subpath within the same TC.As shown in Fig. 7, an exponential distribution shows a good agreement with the empirical CDF for 28 GHz LOS and NLOS scenarios, where the mean intra-cluster excess delay is 3.4 ns and 22.7 ns for 28 GHz LOS and NLOS all set, suggesting a larger intra-cluster delay is usually observed in the NLOS scenario.
5) Cluster Power and Subpath Power: Cluster power is defined as the sum of the subpath powers in the cluster, and the normalized cluster power over the total received power in the PDP can be well modeled by an exponentially decaying function of cluster excess delay with a lognormal-distributed shadowing term, as given in Table IV.P0 is the mean power in the first arriving TC, Γ is the cluster decay time constant, and Z n is a lognormal distributed (normal in dB scale) shadowing term for the cluster power with zero-dB mean and standard deviation σ Z .P n represents the cluster power so that the sum of P n is equal to the total omnidirectional received power P r .The normalized cluster powers measured in the 28 GHz NLOS scenario is shown in Fig. 8, where P 0 is 0.68, and Γ is 23.6 ns, indicating that the expected first cluster occupies about 68% of the total received power and the expected cluster power is less than 34% of the total received power when the cluster excess time delay is over 23.6 ns.
Similarly, the normalized subpath power over the cluster  power can be modeled as an exponentially decaying function over the intra-cluster excess delay, as given in Table IV.Π0 is the mean power in the first arriving subpath in a TC.γ is the subpath decay time constant, and U m,n is a lognormal distributed shadowing term for the subpath power with zero-dB mean and standard deviation σ U .Fig. 9 shows that the first subpath in the cluster is about 42% of the cluster power on average, suggesting a relatively large RMS intra-cluster DS.
The expected subpath power is less than 21% of the cluster power when the intra-cluster excess time delay is over 9.2 ns.

B. Spatial Channel Parameters
1) The Number of Spatial Lobes: An SL represents a main direction of arrival or departure.The angular resolution of the measured APS depends on the antenna HPBW, which are 30°and 8°in 28 and 140 GHz measurements, respectively.A linear interpolation of the measured directional powers with 1°angular resolution in the azimuth and elevation planes was implemented to model the 3-D spatial distribution of the received power.The SLT is -15 dB below the peak power.Measurement results show that there are at most two main directions of arrival in the azimuth plane, except that there are a few NLOS locations measured at 28 GHz which observed three main directions of arrival, as shown in Fig. 10.Thus, a simple DU distribution is used to characterize the number of spatial lobes, which is given in Table IV.2) Mean Direction of Spatial Lobes: Each SL has a mean direction in the azimuth and elevation planes.A simple partition can be applied to generate the azimuth mean direction of an SL by equally dividing the azimuth plane into several sectors, each of which corresponds to an SL.The elevation mean direction of an SL is modeled as a normal random variable N (µ l , σ l ), as given in Table IV, where θ i is defined with respect to the horizontal plane.Considering that the TX height is usually higher than the RX height in a downlink of base station to mobile device setting, µ l of ZOD is typically negative, and µ l of ZOA is typically positive, which represents that ZOD and ZOA are below and above the horizon, respectively.
3) Subpath Angular Offset: For each spatial lobe, the RMS lobe AS is extracted from the partitioned AOA and AOD APS.A generated subpath is randomly assigned to one of the generated SLs, and the angles of this subpath (i.e., AOD, ZOD, AOA, and ZOA) are calculated by adding angular offsets with respect to the mean angle of the SL.The angular offset follows a normal distribution with zero mean and a standard deviation of the median of the measured RMS lobe AS, as given in Table IV.Such angular offset generation deviates from the 3GPP TR 38.901 channel model where angular offsets of 20 MPCs in a cluster are constant [29].

C. Discussions
Each temporal and spatial channel parameter discussed above is generally fitted well by an identical distribution for 28 and 140 GHz, but the values of each parameter for these two frequencies are quite different.The channel at 140 GHz has fewer time clusters and fewer subpaths within each cluster than the channel at 28 GHz.Greater partition loss and higher path loss in the first meter of propagation distance at 140 GHz cause a smaller signal propagation range (the difference of maximum measurable path loss between two frequencies has been considered); thus, some of RX locations which could receive signals at 28 GHz were in outage at 140 GHz.

VII. SIMULATION RESULTS
The statistical channel model presented in Section VI was implemented in an indoor channel simulator based on NYUSIM outdoor channel simulator to investigate the accuracy of the simulated temporal and spatial statistics by comparing with the measured statistics.Note that the parameters listed in Table V are primary statistics used in the channel parameter generation procedure.The metrics used in this section for channel validation are secondary statistics such as RMS DS and RMS AS, which are not explicitly used in the channel generation, but the simulated and measured secondary statistics should yield good agreements.10,000 simulations were carried out for each of four frequency scenarios (i.e., 28 GHz LOS, 28 GHz NLOS, 140 GHz LOS, and 140 GHz NLOS) presented in this work by generating 10,000 omnidirectional and directional PDPs, and 3-D AOD and AOA PASs as sample functions of (3) using the NYUSIM indoor channel simulator.

A. Simulated RMS Delay Spreads
The RMS DS describes channel temporal dispersion, a critical metric to validate a statistical channel model.Fig. 11 shows the simulated and measured omnidirectional RMS DS at 28 GHz and 140 GHz in LOS and NLOS scenarios.As shown in Fig. 11, we obtained the empirical and simulated medians as 10.8 and 10.8 ns for the 28 GHz LOS scenario, 17.0 and 16.7 ns for the 28 GHz NLOS scenario, 3.0 and 2.6 ns for the 140 GHz LOS scenario, and 9.2 and 6.7 ns for the 140 GHz NLOS scenario.The simulated CDFs yield good agreements with the empirical CDFs for four frequency scenarios.By applying the antenna pattern, the directional CIR based on (3) can be given by where g TX ( − → Θ ) and g RX ( − → Θ ) can be arbitrary 3-D TX and RX complex amplitude antenna patterns.The antenna pattern of horn antennas were used in directional CIR simulations in NYUSIM to compare with measured directional RMS DS from 28 GHz and 140 GHz measurements.The antenna gain of a horn antenna can be calculated using the given antenna HPBW, which was given by ( 45)-( 46) in [38].For each omnidirectional channel realization (i.e.PDP), the simulated horn antenna was pointing to the direction of each generated MPC, which output the same number of directional RMS DSs as the number of MPCs.The comparison between the measured and simulated directional RMS DS is shown in Fig. 12.The simulated TX and RX antennas have 15 dBi gain with

B. Simulated RMS Angular Spreads
The omnidirectional azimuth and elevation AS describe the angular dispersion at a TX or RX over the entire 4π steradian sphere, also termed global AS.The AOA and AOD global AS were computed using the total (integrated over delay) received power over all measured azimuth/elevation pointing angles.The measured and simulated global AOA RMS AS were calculated using Appendix A-1,2 in [29].Fig. 13 shows that the simulated and measured median global ASs match well for 140 GHz but not well for 28 GHz due to the difference in the measured and simulated statistics of spatial lobes and the limited number of data samples.The simulated number of spatial lobes was uniformly distributed, which cannot perfectly recreate the specific statistics of spatial lobes measured in this environment, but may be well generalized to measurement data from various indoor office environments.The sheer increasing of the cumulative probability from 0 to 0.6 within 10°at 140 GHz indicates only one SL (one main direction of arrival) observed at the RX, where the global RMS AS would be close to the lobe RMS AS as shown in Fig. 14.
The directional azimuth and elevation AS describe the degree of angular dispersion in a certain direction, which can be regarded as the lobe AS due to the definition of SLs.A -15 dB SLT was applied to obtain SLs before calculating the lobe RMS AS.The simulated and measured AOA RMS lobe AS for 28 and 140 GHz NLOS scenarios are compared in Fig. 14, where the median lobe AS of measured and simulated channels show an excellent agreement (within 0.5°).The measured lobe AS for 28 GHz is larger than the measured

VIII. CONCLUSION
The paper presented a 3-D spatial statistical channel model for mmWave and sub-THz frequencies in LOS and NLOS scenarios based on the extensive measurements at 28 and 140 GHz in an indoor office building.The omnidirectional and directional CI path loss models were derived from measurements, suggesting that NLOS propagation at both frequencies experience similar path loss over distance after removing the effect of the first meter of free space propagation loss.The extracted channel statistics showed that the number of TCs and the number of subpaths within each TC decrease as frequency increases.The channel generation procedure was listed step by step in Table IV, and the values for required parameters obtained from 28 and 140 GHz LOS and NLOS measurements were given in Table V.The indoor channel simulator NYUSIM 3.0 based on the presented statistical model was used to generate tens of thousands of PDP and APS samples.The simulated secondary channel statistics (i.e., omnidirectional and directional RMS DS, global and lobe RMS AS) yielded good agreements with the measured channel statistics.The empirical channel statistics and corresponding unified statistical channel models across mmWave and sub-THz frequencies will provide insights for future propagation measurement and modeling in such frequency range and support analysis and design of 6G indoor wireless systems and beyond.

APPENDIX
The processed data used to generate and calibrate the omnidirectional channel models in this paper are given in Table VI and VII as follows.
TABLE VI: 28 GHz omnidirectional channel statistics with corresponding environment (Env.),TX IDs, RX IDs, T-R separation distance (T-R) in meters, path loss (PL) [11] in dB, the number of TCs (#TC), the number of SPs (#SP), and the omnidirectional RMS DS in ns.

Fig. 1 :
Fig.1: Floor plan of the 9th floor, 2 MetroTech Center[11].TX and RX locations measured at both 28 GHz and 140 GHz are denoted as stars and circles with checkerboard texture, respectively.RX locations only measured at 28 GHz are denoted as solid circles.Each of the five TX locations is denoted in a different color, and the RX locations paired with a TX location is denoted in the same color.

(a) 28
GHz indoor directional CI path loss model and data[11].(b)142 GHz indoor directional CI path loss model and data.

Fig. 2 :
Fig. 2: 28 GHz and 142 GHz indoor directional CI path loss models and scatter plots with TX antenna height of 2.5 m and RX antenna height of 1.5 m for V-V polarization.

(a) 28
GHz indoor omnidirectional CI path loss model and data[11].(b)142 GHz indoor omnidirectional CI path loss model and data.

Fig. 3 :
Fig. 3: 28 GHz and 142 GHz indoor omnidirectional CI path loss models and scatter plots with TX antenna height of 2.5 m and RX antenna height of 1.5 m for V-V polarization.
(a) The number of TCs in the NLOS scenario.(b) The number of TCs in the LOS scenario.

Fig. 4 :
Fig. 4: Histograms and Poisson distribution fittings of the number of TCs of 28 GHz all set, 28 GHz common set, and 140 GHz common set in the (a) NLOS scenario (b) LOS scenario.

A. Temporal Channel Parameters 1 )
The Number of Time Clusters: TCs are obtained by partitioning the measured PDPs based on the MTI.In Fig.4, the empirical histograms of the number of TCs (N ) of three datasets (i.e., 28 GHz all set, 28 GHz common set, and 140 GHz common set) for the LOS and NLOS scenarios with a 6 ns MTI are shown to be well fitted by the Poisson distribution.Since the Poisson distribution starts from zero while the number of TCs is at least one.Thus, N = N − 1 is used for distribution fitting, and the maximum likelihood estimator (MLE) of the parameter λ c of the Poisson distribution is the sample mean of N .The simulated number of TCs from the Poisson distribution is added by one to obtain the actual Fig. 5b shows that the a simple DE distribution matches the empirical histogram of 140 GHz LOS and NLOS common sets since the optimal β NLOS = 1.The mean number of M n (a) The number of subpaths for 28 GHz all set.(b) The number of subpaths for 140 GHz common set.

Fig. 5 :
Fig. 5: Histograms and composite distribution fittings of the number of subpaths of (a) 28 GHz all set and (b) 140 GHz common set.

Fig. 10 :
Fig. 10: Histograms of the number of AOA spatial lobes for 28 GHz and 140 GHz common sets.

Fig. 11 :
Fig. 11: Omnidirectional RMS DS for 28 GHz and 140 GHz LOS and NLOS scenarios.Meas stands for measurement, and Sims stands for simulations.

Fig. 12 :
Fig. 12: Directional RMS DS for 28 GHz and 140 GHz LOS and NLOS scenarios.The simulated TX and RX antenna HPBWs for 28 GHz and 140 GHz in the azimuth and elevation plane are 30°and 8°, respectively.

Fig. 13 :
Fig. 13: Simulated and measured RMS global AOA AS for 28 and 140 GHz LOS and NLOS scenarios.values for 140 GHz, which may be partly attributed to the difference in antenna HPBW (30°and 8°HPBW in 28 GHz and 140 GHz measurements, respectively).

Fig. 14 :
Fig. 14: Simulated and measured RMS lobe AOA AS for 28 and 140 GHz LOS and NLOS scenarios.

TABLE I :
ACRONYMS

TABLE II :
SPECIFICATIONS FOR THE 28 GHZ AND 142 GHZ SLIDING CORRELATOR CHANNEL SOUNDING SYSTEMS Table VI and Table VII in the Appendix give TX-RX location pairs used in this paper.

TABLE III :
TX/RX ANTENNA SWEEP DESCRIPTION

TABLE IV :
INPUT PARAMETERS FOR CHANNEL COEFFICIENT GENERATION PROCEDURE

TABLE V :
VALUES OF REQUIRED PARAMETERS IN THE CHANNEL GENERATION PROCEDURE DERIVED FROM 28 GHZ ALL SET, 28 GHZ COMMON SET, AND 140 GHZ COMMON SET FOR LOS AND NLOS SCENARIOS.
30°HPBW and 27 dBi gain with 8°HPBW in both azimuth and elevation planes in 28 and 140 GHz channel simulations.The measured directional RMS DSs in the LOS and NLOS scenarios are close with respect to both 28 GHz and 140 GHz.Furthermore, the median values of the measured and simulated directional RMS DS yield good agreements in the 28 GHz and 140 GHz LOS and NLOS scenarios.

TABLE VII :
140 GHz omnidirectional channel statistics with corresponding environment (Env.),TX IDs, RX IDs, T-R separation distance (T-R) in meters, path loss (PL) in dB, the number of TCs (#TC), the number of SPs (#SP), and the omnidirectional RMS DS in ns.