Grid-Based Channel Modeling Technique for Scenario-Specific Wireless Channel Emulator Based on Path Parameters Interpolation

Channel model is one of the key components of a wireless channel emulator (WCE) to realistically examine and assert a wireless communication system. The deterministic channel model (DM) enables the simulation of site-specific propagation with high accuracy. However, it is less preferable for WCE due to computational complexity, especially in a non-stationary environment. In this paper, a grid-based channel modeling technique is formulated to realize a site- and scenario-specific channel response in WCE. To address DM’s computational complexity, path parameters are pre-computed at reference locations (e.g., grid nodes) and stored in WCE database. Then, the time-variant channel response is synthesized via path parameters interpolating from those of the adjacent grid nodes. The applicability of the path interpolation technique is investigated analytically. The proposed channel modeling technique is rigorously validated with the ray tracing simulation result at sub-6 GHz and millimeter wave bands. The effects of antenna types, path clustering, and association algorithm to the interpolated channel are investigated. The performance in the rich multipath environment is evaluated in the office space at sub-6 GHz band.


I. INTRODUCTION
R ECENT technological advancement in communica- tion has driven many industries toward automation systems, such as cooperative intelligent transport systems [1], drone-assisted vehicle communication system [2], smart factory [3], smart farm [4], and smart healthcare [5].In these systems, many wireless communication links between sensors are constantly used.Evaluating the ability of wireless communication links in these systems is necessary to optimize cost and performance [6].The field-measured exhaustive drive test method can be implemented for system validation, but is time and cost inefficient [7].An alternative approach uses a wireless channel emulator (WCE) [8], [9] to design and evaluate a large-scale wireless system based on a synthesized channel model.
In general, WCE reproduces real-world wireless communication in a virtual space based on a radio wave propagation characteristic synthesized from a channel model [10].The 3rd generation partnership project (3GPP) also adopts the use of emulators for the performance testing of the mobile network system [7].The performance of the emulator is significantly dependent on the accuracy of the channel models applied.For example, the work in [1] showed that the block error rate curve of a vehicle-to-vehicle wireless system can vary depending on the channel model.As demonstrated in previous studies [11], [12], the statistical properties, including the K-factor, delay spread, path gain, and spatial correlation function, can exhibit significant variations when different channel models are employed in dynamic scenarios.Generally, a channel model is classified into non-geometrical stochastic (NGSM), geometric-based stochastic (GBSM) and deterministic channel (DM) models [13].NGSM generates a wireless channel based on purely stochastic distributions of radio propagation parameters derived from a site-general measurement campaign [14].GBSM randomly places fictitious interacting objects (IO) based on the stochastic distribution of spatial propagation parameters without the need for an actual environment map [11], [12].With this spatial information, GBSM is capable of simulating a dynamic channel by tracking the movements of Tx, Rx, or IOs and updating its geometric information accordingly [10], [15].Due to the ability to perform link-and system-level simulations with relatively low computational cost, GBSM are widely applied in many commercial WCEs for performance evaluation of wireless communication systems such as Keysight's Propsim channel emulation, Anritsu's MT8000A radio communication test station, and Rohde & Schwarz's CMW500 wideband radio communication tester.Recent work in [10], [13] has also proposed new channel models for WCE based on GBSM to optimize the computation time of channel emulation in 3D non-stationary scenarios.However, these stochastic models can only simulate a site-general channel such as urban, suburban, and rural areas.Therefore, it may not be preferable for the automation systems mentioned above, especially in drone-assisted vehicle communication [2] where precise action or manipulation of sensors in the physical world becomes essential.Therefore, WCE with site-specific emulation capability is needed in these wireless systems.In fact, a concept of cyber-physical wireless emulator was described in [8], [16] to allow high-precision prediction of site-specific propagation and system evaluation.
On the other hand, DM computes the channel based on the deterministic interaction between the radio wave and objects in the site-specific environment [17], as shown in Fig. 1a.Compared to its stochastic counterpart, DM has better prediction accuracy and is able to simulate a site-and scenario-specific channel.However, the high computational complexity becomes the main drawback, especially in a large-scale and dynamic environment.Generally, DM is referred to radio simulation, which applies electromagnetic theory (EM), together with morphological and geographical information about the environment, to estimate a channel impulse response (CIR).Among them, the ray tracing (RT) simulation technique is one of the commonly used DM due to its relatively light computation [18].When dealing with a non-stationary scenario, RT may divide the motion into multiple stationary snapshots.With sufficient sample, the time-variant CIR (TV-CIR) can be estimated as a sequence of each channel sampled [19], but with a huge increase in computational cost and hardware.Hence, it is challenging to practically implement the DM into WCE with real-time computation.Therefore, this paper focuses on exploring a novel channel modeling technique for a DM-based WCE that enables large-scale emulation in a static and dynamic environment with real-time computation.The contributions of this work are summarized as follows: • The authors propose the grid-based channel modeling technique allowing the feasibility of the DM-based WCE with real-time computation.The proposed technique interpolates path parameters between a base station (BS) and a mobile station (MS) moving along an arbitrary trajectory from a set of pre-computed DM parameters at reference positions (i.e., grid node) that are stored in the WCE database, as illustrated in Fig. 1b.Therefore, its time-varying channel characteristic can be synthesized from these interpolated path parameters.Hence, DM's high computation is replaced by the light computation of the interpolation process during channel emulation.The authors briefly described the overview architecture and practical challenges such as path clustering and path association in [20].• The authors analytically derive the interpolated parameters to assert the limitations of the proposed channel modeling technique.Through analytical investigation, the source of the interpolation error can be approximately quantified.
• The authors performed a rigorous simulation to validate the performance of the proposed technique by comparing it with the result obtained directly from the RT simulation.The simulation is conducted in sub-6 GHz and millimeter wave (mmWave) bands, with halfwave dipole and horn antennas to investigate the robustness of the proposed technique in different frequency bands and antenna types.The effect of path clustering and association is also investigated.
As an example of practical implementation, another experiment was conducted in an office space to observe the performance of the interpolated channel in a rich multipath environment in the sub-6 GHz band.The remainder of the paper is organized as follows.In Section II, the proposed grid-based channel modeling technique is first mentioned, which includes the formulation of TV-CIR and the path parameter interpolation and smoothing techniques.Analytical investigation of the interpolation of path parameters is addressed in Section III to discuss the limitation and applicability of the proposed modeling technique.The performance of the proposed modeling is numerically validated and discussed in Section IV.Finally, the key findings are summarized in Section V.

II. GRID-BASED CHANNEL MODELING TECHNIQUE A. TIME-VARIANT INTERPOLATED CHANNEL MODEL
In general, the number of propagation paths used in WCE is typically limited according to its hardware limitation [7], [10], [15].Therefore, the sum-of-sinusoid tapped delay line (SoS-TDL) model [10], [17] is usually chosen to represent TV-CIR due to its ability to express multipath propagation in the time domain with a few discrete delay taps.Each tap displays the temporal fading characteristic from a coherent combination of a discrete set of rays (i.e., subpaths) per tap.For simplicity, let us consider a dynamic scenario between fixed BS and moving MS. Continuous time can be discretized in WCE with a combination of the sampling interval of the wireless signal t s and the sampling interval of the emulator orchestrator t o .The former has a relatively shorter interval to adequately approximate the continuous fading characteristic.The latter generally has a longer interval to update the position and velocity of MS by the orchestrator.Also known as a scenario generator in [9], the orchestrator is the emulator unit that coordinates and manages all elements that move in the cyberspace of the WCE according to user input.Hence, a discrete-time instant can be defined as where i ∈ N + and j ∈ N + | j ≤ t o / t s are time sampling indices of t o and t s , respectively.These time indices are valid for all time samples.The determination of both t o and t s depends on the limitation of the WCE hardware.However, the latter should satisfy the Nyquist sampling criteria such that t s ≤ 1 2ν max where ν max is the maximum Doppler frequency [17].
Suppose that N i ≤ N max denotes the number of tapped delays updated at the ith time index.Each nth tapped delay consists of M ni ≤ M max number of rays.Here, the maximum number of tap delay and rays per tap are defined as N max and M max , respectively.Therefore, the SoS-TDL channel model can be expressed in terms of TV-CIR as, where τ n (t ij ), C mn (t ij ), and ν mn (t ij ) are SoS-TDL parameters representing the nth tapped delay with α n (t ij ) weight, magnitude of the mth ray's weight, and the mth ray's Doppler frequency that induced phase rotation mn (t ij ), respectively.These three SoS-TDL parameters are time-varying w.r.t. the temporal MS position.Unlike GBSM and NGSM, SoS-TDL parameters with the grid-based channel modeling technique are interpolated from those reference channel parameters at neighboring grid nodes, as shown in Fig. 1b, via the procedure described in Section II-B.Furthermore, the range and resolution of the SoS-TDL parameters used in WCE are finite due to hardware limitations [7].To determine the appropriate range and resolution, limitations on the emulation scenario should be taken into account, such as the maximum speed of MS, frequency band, and area of the emulated scenario, since these physical parameters directly influence the boundary of the SoS-TDL parameters.

B. SOS-TDL PARAMETER INTERPOLATION AND SMOOTHING
In the grid-based channel modeling technique, the emulated 3D space is uniformly divided into a small grid cube structure (or a square in the 2D space), as shown in Fig. 2. Each grid has a unique set of nodes (e.g., grid's vertices) that contain the reference spatial propagation characteristic between BS and the grid node in the same tap/ray structure of the SoS-TDL model.Such a structure is similar to the clustered delay line model in the 3GPP standard [21].Specifically, each nth tapped delay of the node's reference parameter contains the mth ray's dual-polarized path weight, angle-of-arrival (AoA), and angle-of-departure (AoD).
Consider a simple scenario in which the MS position at t ij is located within a particular grid and, therefore, surrounded by a set of D number of neighboring nodes with grid resolution g as shown in Fig. 2. A set of neighboring nodes is represented by G(t ij ) = {g 1 , g 2 , . . ., g D }, with each g d being the position vector of the dth neighboring node in a Cartesian coordinate system.With the MS position provided by WCE orchestrator, WCE can determine this MSassociated grid by inquiring its database to search for the grid whose all nodes surround the MS position.As conceptualized in [20], when the WCE is running, only the reference spatial parameters of these neighbors are converted to SoS-TDL parameters given the updated information on the MS velocity and antenna patterns at every interval t o .Generally, C mn (t ij ) is computed by applying the antenna radiation patterns of MS and BS to the reference dual-polarized ray's weight parameters [21], while ν mn (t ij ) is the projection of the AoA parameters in the MS velocity direction normalized by the wavelength [17], [21], respectively.Note that the contribution of AoD is omitted in this case, since BS is presumably fixed.On the other hand, τ n (t ij ) remains the same during conversion.
Suppose that the reference spatial parameters at the dth grid node have already been converted to the SoS-TDL parameters and are denoted by τ nd (t ij ), C mnd (t ij ), and ν mnd (t ij ).Since the SoS-TDL parameters can be updated every time instant t i1 according to the latest information from the orchestrator, the temporal change can be described with a piecewise function such that Consequently, the interpolated SoS-TDL parameters of MS located at p i1 will also change in a piecewise manner as shown in Fig. 3.With this model, the interpolated channel of any jth time index within each t o interval can be considered as a local stationary region of a deterministic process, as depicted in Fig. 4.This structure allows the WCE to estimate the characteristic of statistical fading in terms of the local scattering function (LCF) [22].Therefore, channel statistical parameters such as delay and Doppler spreads and the Rician K-factor can be obtained to evaluate the performance of the communication link at each t o interval.However, this abrupt transition unavoidably results in a discontinuity of the TV-CIR.This discontinuity is visualized by the dashed lines in Fig. 4 in the delay domain and the channel envelope at the beginning of the t o interval.To synthesize the continuous TV-CIR characteristic, an additional modification is needed to smoothen the interpolated SoS-TDL parameters.
According to (2c), a piecewise change of ν mnd (t) in (3c) fortunately can produce a continuous linear rotation of the TV-CIR phase in the time domain, as seen in the middle graph in Fig. 3, and thus does not require modification.Nevertheless, as depicted with the gray line on the right graph of Fig. 3, non-negligible discontinuity may occur from an abrupt change of C mnd (t ij ) in (3b) given a longer t o interval.Therefore, a linear weight transition is introduced to smooth out the abrupt change at each t i1 .Note that the phase term of C mnd (t ij ) is omitted under the assumption that the fast fading is strongly dominated and continues to rotate by the ray's Doppler.For simplicity, let us assume that these rays with the same mth and nth indices are already associated among D neighboring grid nodes.Hence, the set of the associated weight and Doppler can be defined by L c,mn 2, . . ., D , respectively.Both rays' associated Doppler and weight corresponding to the MS position at t ij can be calculated by where is a linear slope of ray's weight and F (•) is an interpolation function.The latter function computes interpolated SoS-TDL parameters at an arbitrarily MS position p i1 from a set of reference parameters located at the neighboring node positions in the set G. It should be noted that the weight of the interpolated ray C mn (t ij ) becomes a real value due to the aforementioned interpolation assumption.As shown in Fig. 3, a gray line with a piecewise transition is replaced by a linear slope of C mn (t ij ) within β t s duration given a length coefficient β that can be arbitrary set between However, its drawback is the time lag by β t s within a transition region shown in Fig. 4.
Similarly, the interpolated tapped delay at the MS position can be calculated by where . ., D is a set of associated delay taps among grid nodes.At first glance, interpolating the tapped delay seems to only change the location of the tapped delay in the TV-CIR in (2a) without affecting either magnitude or phase.However, when considering a received signal waveform, its phase will rotate due to the delay bin.As a consequence, a sudden change in the  tapped delay depicted in Fig. 4 introduces discontinuity in the received signal.To mitigate an abrupt phase rotation due to the movement of the tapped delay, an additional phase term is introduced at α n (t ij ) as, where f 0 is a carrier frequency.By substituting (2b) with (7), this manipulation will effectively cancel out delay-induced phase rotation, thus preserving the continuous fading characteristic of the received signal due to the interpolated mth ray's Doppler component.

III. ANALYTICAL INVESTIGATION OF THE SOS-TDL PATH PARAMETERS INTERPOLATION
The performance of the SoS-TDL channel model in (2a) depends on the accuracy of the interpolated SoS-TDL parameters.Therefore, the applicability of the proposed modeling technique and the source of error can be identified by investigating the interpolated parameters.Without loss of generality, the following investigation focuses only on the analytical derivation of a single-ray propagation geometry with a single reflection mechanism from BS to MS and to all D neighboring grid nodes, as depicted in Fig. 5.In practical implementation, the path clustering algorithm is required to group a set of deterministic paths in the same tap delay, while the parameter interpolation scheme also requires spatial association among D nodes [20].Since these algorithms introduce additional artifacts to the interpolated parameters, a perfect association of rays among D nodes and a single ray per tap condition are assumed.It should be noted that these processes are non-trivial in the practical implementation of the proposed modeling technique, and thus require further in-depth investigation.Although it is beyond the scope of the following analysis, its effect will be numerically discussed through simulation in Section IV.Suppose that the interpolation function in this analysis is formulated in the weight summation form to interpolate a path parameter x intp such that where x d and w(g d , p) are the interpolating parameters, and the spatial weighted function such that g d ∈G w(g d , p) = 1, respectively.For brevity, the weight notation is shortened to w d = w(g d , p), and the subscripts m, n, and the notation t ij will be omitted in this analysis.

A. INTERPOLATED SOS-TDL TAP DELAY
Due to the reflection mechanism, the reflected rays of the BS-node and BS-MS path trajectories (e.g., link) as depicted in Fig. 5 can be unfolded into virtual links of the line-ofsight (LoS) path [23].Let us define τ MS as the ideal tap delay of the BS-MS link.From the geometry of MS and grid positions, the tap delay of link between BS and the dth node can be expressed in terms of τ MS as, where ψ BS,d and τ d are the angle between BS-MS and BS-node links, and the delay offset from τ MS to τ d cos ψ BS,d , respectively.Substituting ( 9) into ( 8) yields a relation between the interpolated delay τ and τ MS , by From the above equation, the interpolation condition with τ = τ MS may ideally be achieved through designing w d such that the second term is minimized while the inverse cosine factor in the first term becomes unity.
It is intuitively obvious that 1 cos ψ BS,d and τ d degrade the interpolation performance, which depend greatly on the distance between BS-MS and g, respectively.For example, suppose that the BS-MS distance is larger relative to the grid size, cτ MS g where c is the light speed.The angle ψ BSd → 0 and Taylor expansion [24] can be applied to the inverse cosine factor such that the approximation becomes where ε τ is the relative residual of the interpolated delay which is expressed by With a large relative distance assumption, the accuracy of the interpolated delay will gradually improve as τ τ MS , and thus ε τ will become negligible.Therefore, this example scenario can be considered as the optimal condition for the delay interpolation.

B. INTERPOLATED SOS-TDL DOPPLER
In this scenario, the Doppler frequency is the projection of the AoA direction along the direction of motion [17], [21], which is described by the MS velocity vector v.This projection is expressed by a factor cos ϑ where ϑ is the Doppler angle between the AoA direction and v.As depicted in Fig. 5, the AoA directions of the BS-MS and BS-dth node rays are different due to the rays' propagation geometry, and thus their Doppler angles defined as ϑ MS and ϑ d , respectively.Similar to the delay case, the ray's Doppler frequency of the dth node is formulated in terms of ϑ MS as where λ and ϑ d = ϑ MS − ϑ d are the wavelength, and the Doppler angle offset, respectively.The interpolated Doppler can be expressed as the weight summation of D Doppler angles by substituting ( 13) into (8) as where ν MS = ( v /λ) cos ϑ MS is the Doppler frequency of BS-MS link.This relation shows that the interpolated Doppler becomes more precise as ϑ d approaches zero.This could happen when the reflected points are further away from the MS compared to the size of g.Under this condition, it may be possible to truncate the Doppler angle offset with the first order Taylor expansion such that cos ϑ d ≈ 1 and sin ϑ d ≈ ϑ d .Hence, ( 14) can be simplified into where ε ν is the relative residual of the interpolated Doppler formulated by Here, the approximated form in this scenario confirms the performance of the interpolated Doppler as ε ν becomes negligible with a smaller ϑ d .

C. INTERPOLATED SOS-TDL WEIGHT
According to [20], [21], the weight of SoS-TDL ray can be expressed as the multiplication of the dual-polarized complex ray's weight matrix to both MS and BS antenna radiation pattern in its respective polarization.However, for simplicity, the analysis of the weight interpolation will consider only a single polarization.Since only the magnitude term is considered for the interpolated SoS-TDL weight as mentioned in Section II-A, the phase term will not be are the azimuth and zenith angles of the BS-dth node link.Superscripts A and D associate those angles with the AoA and AoD quantities, respectively.Therefore, the weight of the SoS-TDL ray at the dth node becomes where R(ϕ d ) and s d are a magnitude of a single polarized Fresnel reflection coefficient and the propagation distance, respectively, of BS-dth node link.Here, ϕ d is the incident angle relative to the normal vector of the reflected surface w.r.t. the BS-dth node link.Similarly, the interpolated ray weight can be expressed by Similar to the delay and Doppler interpolation, the propagation distance, reflected angle, AoA, and AoD of BS-dth node link can be expressed by where the first term and the second term on the right-hand side of the above equations are actual propagation parameters of the BS-MS link, and the parameter offsets to those of the BS-dth node link, respectively.Therefore, the actual weight of the BS-MS link, C MS , can be formulated similarly to (17) but with the BS-MS propagation parameters indicated above.Let us assume the scenario used in the delay and Doppler interpolation: larger distances of the BS-MS link and the MS-reflected point link.The first order Taylor expansion can be applied to (18) d and s d , to express the approximated form of the interpolated weight as, where ε C is the relative residual of the interpolated weight expressed as The superscript prime represents the first order derivative of the function, while the additional subscripts θ and φ on the antenna radiation functions indicate the partial derivation direction along the azimuth and zenith angular variables, respectively.Due to the highly nonlinearity and multiple variables in the path weight model, the linear approximation of the weight error term may be less accurate than other interpolated SoS-TDL parameters.However, the approximation shows the interpolation capability as s MS s d and θ and ϕ d approach zero.

A. RIGOROUS VALIDATION OF THE INTERPOLATED CHANNEL 1) SIMULATION SETUP
In this experiment, the empty hallway space with a partition wall and a doorway was constructed as a simulation scenario to numerically validate the applicability of the interpolated SoS-TDL parameters and the channel.The geometry of the hallway space and the positions of the BS and two MS trajectories are shown in Fig. 6.The materials in the hallway and the partition wall were set to be perfect electrical conductors (PEC), with infinite permittivity and conductivity [25].
The channel response along the A-trajectory represented the LoS scenario, whereas the B-trajectory exhibited the transition of channel response between LoS and non-line-ofsight (NLoS) scenarios due to the presence of the partitioner and doorway.The reference grid nodes were assigned at the same height as the MS and covered the entire two trajectories.MS along the A-trajectory were initially placed in the middle of the grid system on the y-axis and moved along the x-axis, and vice versa for the B-trajectory.With this 2D grid structure, there were four nodes surrounding MS at any given time.The reference spatial path parameters of the links between the BS and D nodes were simulated using Wireless InSite [26], the RT propagation simulation software.These reference parameters were generated at each grid node and used to synthesize the interpolated channel.For validation, the reference SoS-TDL parameters and channel along the same trajectories were simulated using RT with a 1-mm spatial resolution.The interpolated channel is quantitatively validated in terms of normalized Doppler and delay spread because these quantities are essential for the evaluation of wireless communication performance [17].
The maximum Doppler frequency and the symbol duration were used as normalization terms for Doppler and delay spread, respectively.The latter is inversely proportional to the signal bandwidth.To avoid the effect of the window function, the interpolated SoS-TDL parameters were used to calculate both spread quantities at every t o interval.The relevant simulation parameters are listed in Table 1.
The first simulation in Section IV-A2 assumes a perfect path association among the nodes, and the emulator has no limitation in terms of N max .Therefore, the performance of the interpolated channel and its parameters can be rigorously validated without applying any path clustering and association algorithm.In other words, the SoS-TDL channel has only M max = 1 ray per tap.The halfwave dipole antenna of both BS and MS has a half-power beamwidth (HPBW) of 78 • and a directivity gain of 2.156 dBi [27].The perfect path association condition was achieved by comparing the positions of the IOs along each propagation path trajectory observed from the links of the BS and D nodes.The A-trajectory was considered to avoid the birthand-death path behavior [13] during the transition of the LoS-NLoS scenario.Only 10 RT path parameters with the highest path gain were selected.Initially, the performance of the interpolated channel was evaluated at 5.25 GHz.
The second simulation in Section IV-A3 examines the mmWave channel at 28 GHz with both halfwave dipole and horn antennas to assess the robustness of the interpolated channel across different frequency bands and antenna types.Therefore, another set of pre-computed spatial parameters was generated at each node for the mmWave channel.The horn antenna radiation pattern was generated using Wireless InSite, based on the physical geometry of the commercially available Pasternack PE9851/2F-20 antenna This horn antenna has the HPBW of 17.4 • on the horizontal plane and 17 • on the vertical plane, with a directivity gain of 20 dBi.
The third simulation condition in Section IV-A4 investigates the effect of the path clustering and association algorithm on the interpolated channel.The same RT path parameters as in the previous simulation at 5.25 GHz were used for the A-trajectory.The B-trajectory was also considered to investigate artifacts from path clustering and association algorithms.The K-powermean clustering algorithm [29] was applied to the RT path parameters to reorganize in the SoS-TDL tap-ray structure explained in Section II.The clustering data was based only on the RT delay parameter, with the RT path gain being the power weight.Therefore, the RT path in the same cluster became a set of M ni rays sharing the same tapped delay, which is determined as the clustering delay centroid.The sortingbased path association method was applied to the reference SoS-TDL parameters of the D neighboring nodes prior to interpolation.Specifically, a set of reference SoS-TDL tapped delays at each node was sorted ascendingly and the tapped delays with the same sorted order were assumed to be spatially associated.The set of SoS-TDL rays from D nodes within the same associated tapped delay was associated in the same manner but based on the Doppler value of the rays.

2) PERFORMANCE OF INTERPOLATED SOS-TDL PARAMETERS AND CHANNEL
The performance of the interpolated SoS-TDL parameters was numerically investigated.The relative residual of the SoS-TDL parameters of the nth tap delay and the mth ray in each ith time index, ˜ ζ,mni = ζmni ζMS,mni − 1, was calculated as the evaluation metric where ζ ∈ {τ, ν, C} and ζMS ∈ {τ MS , ν MS , C MS } represent the interpolated SoS-TDL parameters and the actual SoS-TDL parameters between BS and MS obtained from the RT simulation, respectively.The precision of the interpolated SoS-TDL parameters gradually improves as ˜ ζ,mni → 0, and so does the interpolated channel.Linear approximation forms of the relative residual, ζ,mni , in ( 12), ( 16), and (20a) were calculated to validate the applicability of the analytical model to predict the performance of the interpolation.
Relative residual profiles of the first five paths with the strongest path gain along the A-trajectory with a grid size of 10 cm are shown in Fig. 7.The first path corresponds to LoS propagation, while others are reflected paths bouncing between ceiling and floor.Upon examination of τ,mni (solid line with a circle marker) in Fig. 7a, the approximate linear residual model in (12) was able to accurately predict the residual profiles of the interpolated delay, as the difference between the numerical residual (solid line with a cross marker) was negligible.Since the grid size was significantly small relative to the propagation path, the assumption of a smaller φ BS,d in ( 11) is highly valid.Increasing the grid size in this scenario would significantly show the observable discrepancy.In the case of interpolated Doppler and weight, a larger discrepancy was clearly seen between the approximate linear model and the numerical result of the relative residual, as shown in Figs.7b -7c.Because the reflected point was closer to MS than the total propagation length, ϑ d became slightly larger than φ BS,d , resulting in a larger discrepancy than in the delay case.A similar reason could be applied to the relative residual model of the weight.Additionally, highly non-linearity of the antenna pattern directly contributed to the prediction discrepancy.However, as MS moved farther away from BS, the linearly approximated residual model became more accurate.This investigation has validated the applicability of the relative residual model to predict the performance of the interpolated SoS-TDL parameters, which can be beneficial in determining the optimal grid size for a given emulation environment.
The largest residual was observed at the initial position of the MS, where the propagation distance was the shortest.Even in this scenario, the interpolated residual was still less than 1% given the smaller 10-cm grid size.As MS moved away from BS, the relative residual was significantly reduced.Obviously, the level of relative residual would increase proportionally to the size of the grid.The first path (LoS propagation) had the largest relative residual for the interpolated delay and Doppler, followed by the reflected paths with the least residual at the fifth path.Since the total propagation path and the distance between the last reflected point and the MS increased proportionally in this trajectory, the relative residuals of these paths also decreased in the same way according to (12) and ( 16).On the other hand, the fifth and second paths have the largest and smallest relative residuals of the weight.As the residual weight depends on the antenna pattern and reflection coefficient rather than the propagation geometry according to (20a), its residual levels were different from delay and Doppler.Therefore, the reduction rate of the relative weight residual may vary according to the material and antenna pattern, in addition to the propagation geometry.
The interpolated channel along the A-trajectory was synthesized by substituting the interpolated SoS-TDL parameters into the channel model in (2a).The bandlimited time-variant power delay profile (PDP) was then synthesized to validate the channel of the wideband signal by convoluting the interpolated channel with the signal filter.Here, the signal was shaped using the raised cosine filter with 0.5 roll-off factor, assuming a 400 MHz signal bandwidth.The interpolated PDP depicted in Fig. 8a produced a pattern similar to that computed directly from the RT depicted in Fig. 8b because of a small discrepancy of the interpolated SoS-TDL parameters.The local stationarity period within the t o interval introduced a discrete delay change in the interpolated PDP compared to the continuous change in the PDP obtained directly from RT.As illustrated in Figs.8c -8d, two PDPs at the 1-s and 2.5-s time instances were observed to closely match the RT results.However, there were minor differences in the path gain level at specific peaks in both PDPs, especially the peaks around 60 -80 ns of the 2.5-s PDP.Additionally, a larger difference was observed at the interpolated PDP around the delay with a lower path gain level.In fact, the fluctuation of the path gain level in the interpolated time-variant PDP in Fig. 8a was found to have a slight difference from that of the PDP produced by RT depicted in Fig. 8b.Such a discrepancy is expected since the absolute phase was not considered during the parameter interpolation, and thus a different temporal fading pattern.However, the fading statistic, which is governed by the Doppler spectrum, should still remain correlated to that of the RT channel.The comparison in terms of the narrowband channel in Fig. 8e visibly indicated the difference in the fast fading pattern of the interpolated channel.The random phase approach (RPA) [23], [30] was applied to validate the similarity of the fading statistic based on the confidence interval of the fading statistic.Specifically, the RPA estimates the CDF of the fading statistic describing the probable range of path gain fluctuation at a given time, assuming a coherent summation of all rays with a random phase condition.Here, the path gain level corresponding to the CDF from 2.5% to 97.5% was chosen as the 95% confidence interval of the fading statistic.The fact that most of the narrowband fading channel obtained from RT was bounded inside the confidence interval line, as shown in Fig. 8e, has evidently indicated a high correlation in terms of the fading statistic with the interpolated channel.Furthermore, weight smoothing during the transition region mitigated the transient of the interpolated channel at the border between two consecutive i-th time indices, resulting in the continuously smooth fading channel.
The statistical performance of the interpolated channels generated from different grid sizes was compared with that of the RT channel in terms of the normalized Doppler and delay spread profiles in Fig. 8f and 8g, respectively.Overall, the three interpolated channels could statistically reproduce results similar to the reference RT channel due to small discrepancies in both spread profiles.When the discrepancy was carefully investigated up close, the performance of the interpolated channel in the same scenario slightly deteriorated as the grid size increased.The RMSE of the normalized Doppler spread of the interpolated channel with a 10-cm grid size was 8.41 × 10 −6 .This was significantly lower than the RMSEs of the channels with grid sizes of 50 cm and 100 cm, which were 23.875 and 97.875 times higher, respectively.The channel with a 10-cm grid size interpolation had the lowest RMSE of the normalized delay spread, which was approximately 2.63 × 10 −4 .In comparison, the RMSEs for the channels with 50-cm and 100-cm grid sizes were 19.407 and 68.255 times higher, respectively.In Fig. 8f, a slightly larger normalized Doppler spread of the interpolated channel was observed at 0 -1.5 s compared to the RT result, but the discrepancy was reduced as MS moved away from BS.However, the trend suddenly reversed afterward, with the normalized Doppler spread of the interpolated channel significantly lagged behind the RT result.This is because the Doppler spectrum experienced a considerable alteration over the 1.5 to 2 s period, as evidenced by the accelerated rate of change in the normalized Doppler spread.In fact, such an effect is also seen in Fig. 8e as the fading patterns abruptly shifted during this period from slower to more rapid fluctuations in the path gain level.Nevertheless, the discrepancy gradually decreased throughout the trajectory.In terms of the normalized delay spread in Fig. 8g, its discrepancy gradually reduced along the MS trajectory proportional to the trend of the relative residual of the interpolated delay.Given the ideal path association and the absence of path clustering, the validation results showed that the grid-based channel modeling technique is capable of estimating the time-variant channel.Its accuracy highly depends on the interpolated SoS-TDL parameters, which can be approximately predicted by the proposed relative residual model.

3) EFFECT OF THE FREQUENCY AND ANTENNA TYPES TO THE INTERPOLATED SOS-TDL CHANNEL
With the same simulation setup, the time-variant PDP of the interpolated channel with the 10-cm grid at 28 GHz exhibited a weaker path gain of around 20 dB compared to the channel at the sub-6 GHz band, but with more frequent fading.The agreement on the time-variant PDP with that of the RT result shown in Fig. 9b indicated the robustness of the interpolated channel in different frequency bands.When the horn antenna was applied to the BS, the overall path gain level of the interpolated time-variant PDP increased to a level similar to the sub-6 GHz channel, as shown in Fig. 9c.However, the path gain level of those components with a larger delay (NLoS paths) was significantly attenuated due to the smaller HPBW.The interpolated channel was still capable of generating a time-variant PDP that resembled the result obtained from the RT depicted in Fig. 9d.However, a slight discrepancy in terms of path gain level was observed, especially at the NLoS components around 70 ns.The effect of the antenna pattern on the path gain level can be explained by the relative residual C in (20a) where the antenna pattern contributes to the error during the interpolation of ray weight.To quantify the statistical performance of the interpolated channel across frequency and antenna type, the absolute discrepancies of normalized Doppler and delay spread profiles were chosen as evaluation metrics.Here, the absolute discrepancy is calculated by taking the absolute difference between those two spread profiles of interpolated and RT channels.As depicted in Figs.10a -10b, interpolated channels between two frequency bands have almost the same discrepancies of normalized Doppler and delay spread profiles given the same grid size.On the other hand, changing the type of BS antenna from half-wave dipole to horn antenna drastically increased the absolute discrepancies of delay and Doppler profiles by approximately ten times.Relatively, these discrepancies decreased slightly as MS moved away from BS.Since halfwave dipole is omnidirectional and has wider HPBW, the effect of antenna radiation pattern to the interpolated ray's weight is less prominence than the horn antennas with narrower HPBW.In comparison, enlarging the grid size from 10 cm to 100 cm increased the absolute discrepancies of delay and Doppler spread approximately 100 times in both antenna types.In conclusion, antenna types significantly influence the performance of the interpolated channel than frequency bands.

4) EFFECT OF THE PATH ASSOCIATION AND CLUSTERING TO THE INTERPOLATED SOS-TDL CHANNEL
Despite the minimal distortion caused by interpolation of the SoS-TDL parameters, the implementation of the interpolated channel necessitates a path clustering and association algorithm due to the hardware restrictions of the emulator and the requirement for interpolation with multiple sets of data, respectively.Therefore, the path clustering and association algorithms mentioned in Section IV-A1 were used as an implementation example to numerically investigate additional artifacts in the interpolated channel.When using a condition of N max = 3 tap delays and M max = 5 rays, the interpolated channel experienced more rapid fading than the reference RT, as seen in the time-variant PDP in Fig. 11 for g of 100 cm.This large path gain fluctuation is the result of coherent summation of M n rays at each tap delay as opposed to a single ray per tap in the reference RT channel.Overall time-varying path trajectories could be preserved with a group of multiple paths within 5 -40 ns were collapsed into two taps, while another two paths around 60 -70 ns were merged into one.The perseverance of the path trajectories can be viewed as a certain level of success for the heuristic path association in this ideal LoS scenario.By increasing N max to 5, the time-variant PDP became more similar to that of the RT as depicted in Fig. 11b where the number of paths within 5 -40 ns increased to four with delay trajectories similar to the reference one.However, the path gain level of the path around delay of 35 ns within 0 -0.8 s was incorrectly synthesized in terms of path gain level.This artifact may occur when the path clustering algorithm incorrectly assigned a certain ray among the D grid nodes to different tap delay.Note that the effect of path clustering incorrectness may be decreased when dealing with a narrower signal bandwidth relative to the delay spread.To investigate only the effect of path association, another interpolated channel was synthesized without applying the clustering algorithm (M max = 1 and N max = 10).Observing the normalized Doppler spread in Fig. 11c, the performance was quite similar between the interpolated channels with and without path clustering.However, a larger discrepancy was observed during 1 -2 s when the channel experienced rapid change in terms of Doppler spectrum.This is possibly due to the fact that the grid size g = 100 cm was sparse compared to the actual change in the Doppler spectrum.Therefore, these adjacent grid nodes may contain paths with slight spatial inconsistencies.In such a scenario, incorrect path association may be unavoidable.However, a slightly larger discrepancy was observed during 0 -1 s for the normalized Doppler spread of the interpolated channel with the N max = 5 taps condition.In fact, a similar larger discrepancy was also observed in the normalized delay spread profile, as depicted in Fig. 11d.Since some paths were incorrectly clustered as previously discussed, inappropriate path association within the same grid was expected, resulting in the overall discrepancy of these spread profiles.Despite this interval, the performance of the normalized delay spread increases proportionally to the number of N max .However, a larger number of taps also come with a high possibility of incorrect path clustering.Therefore, the artifacts caused by these processes are highly dependent on each other, especially in this LoS scenario.
Through simulation of the interpolated channel along the B-trajectory, the performance of the interpolated channel in the NLoS and the transition between the LoS-NLoS scenarios can be investigated.The actual time-variant PDP of the first 10 paths with the strongest path gain is visualized in Fig. 12. Considering a confined room behind the partitioner with a relatively longer distance from BS, some propagation paths arrived MS at similar delay and AoA angles, resulting in fewer distinguishable paths given 400 MHz with a lesser fading pattern.The LoS-NLoS transition occurred at around 1.8 and 4.5 s, where the majority of existing paths disappeared and the new paths appeared immediately.A similar transition was also observed within the NLoS region around 0.2 and 5 s, respectively.The interpolated channel was synthesized given the clustering condition of N max = 3 taps and M max = 5 rays with the same path association algorithm.In the event of a smaller grid size ( g = 10 cm), artificial paths were introduced during the abrupt transition, as shown in Fig. 12b.For example, an artificial delay of 175 ns was created to connect the movement of tap delay from/to 190 ns to/from 140 ns during a time instance of 0.25 and 5.3 s, respectively.Although the paths at adjacent grid nodes before and after the transition were spatially inconsistency, they were considered associated by the path association algorithm.Therefore, through interpolation, the artificial SoS-TDL parameters were synthesized as a weighted average based on the relative position of the MS within the grid.When the grid size was sufficiently fined, such an artifact may be mitigated, as seen during the transition at 1.8 s.When the grid size increased to 100 cm, the artificial transition period became wider, but with a more continuous transition, as clearly observed in Fig. 12c.Consequently, the shorter transient characteristics during 0.25 and 5.3 s were completely lost.Similarly, the statistical performance of the interpolated channel during the transition became severed as a larger discrepancy in the normalized Doppler and delay spread profiles observed in Figs.12d and 12e, respectively.As the artificial transition period is proportional to the grid size, selecting the grid size within such a region should be than in other areas.In other words, the grid size could be considered as a scenariodependent parameter rather than uniformly determined for optimizing the interpolated channel performance.

B. PERFORMANCE OF THE INTERPOLATED CHANNEL IN THE ACTUAL INDOOR ENVIRONMENT 1) SIMULATION SETUP
In this experiment, a simulation scenario with surrounding furniture based on actual office space in Tokyo, Japan, was constructed to investigate the performance of interpolated   channels in a rich multipath environment.BS was placed inside the meeting room where the MS was horizontally moving in the common area for 10 m, as visualized in Fig. 13.The furniture in the rooms on the left and right was discarded due to a negligible contribution to the propagation channel between BS and MS.A set of grid nodes was generated along the MS trajectory, where each node consists of 30 DM paths.The positions of BS, MS, and the grid nodes in the simulation environment are visualized in Fig. 13.The grid nodes were generated using Wireless Insite [26] with the same height as MS.Other relevant simulation parameters are summarized in Table 2, and the dielectric properties of the material used in the RT simulation are listed in Table 3.The path clustering and association algorithms described in Section IV-A1 were applied for the generation of the interpolated channel.Similarly to the previous simulation, the RT channel was synthesized as a reference.For comparison, the 3GPP standard channel was stochastically synthesized assuming the NLoS Indoor-Office scenario based on the GBSM channel generation procedure in [21].Procedure B of spatially-consistent UT/BS mobility model was applied to update the propagation parameters and synthesized the 3GPP TV-CIR along MS trajectory.

2) PERFORMANCE EVALUATION OF THE INTERPOLATED CHANNEL
The time-variant PDP depicted in Fig. 14a was simulated using RT with a spatial resolution of 1 mm to represent the reference site-specific channel along the MS trajectory in this experiment.Due to the large number of IOs, many paths in this time-variant channel have a shorter lifetime with a propagation delay ranging from 25 to 200 ns.Based on a default setting of the emulator hardware prototype used in the authors' project, the maximum number of tap delay and rays per tap in the interpolated channel were set to N max = 6 and M max = 10, respectively.With these limitations, the interpolated PDPs produced from grid sizes of 10 cm and 100 cm are comparatively shown in Figs.14b -14c, respectively.With its finer grid size, the 10-cm interpolated PDP still showed some degree of similarity to the reference RT channel, particularly in the time duration of 0 -0.5 s, 3 -4 s, and 6 -7 s, where most paths remained alive.However, a rapid change in tap delays was introduced in the remaining interval.This characteristic was similar to the B-trajectory in Fig. 12b as the adjacent spatially inconsistent paths of birth and death tried to connect with each other according to the applied path association algorithm.Obviously, increasing N max may reduce this artifact, as it adds a degree of freedom to associate the path.Additionally, due to the delay sorting condition of the applied path association algorithm, the path-crossing behavior at delays around 70 ns during a time instance of 1 s could not be reproduced by the interpolated channel.With a coarse grid size of 100 cm, the artifact characteristic of rapid change in tap delays was replaced by a slow transition between adjacent regions that are spatially inconsistency.In this region, the interpolated channel was expected to be less accurate because the transition region was created artificially.This characteristic was previously discussed in a similar situation but in a more rigorous setting in Fig. 12c.A single realization of the standard 3GPP channel with the same number of taps and rays is depicted in Fig. 12d for comparison.Due to the stochastic approach, the 3GPP time-variant PDP failed to reproduce a site-specific characteristic of the reference RT channel in terms of path gain level and delay line distribution.In contrast, the interpolated channel is superior in the aspect of reproducing the site-specific characteristic to some extent, despite having an artificial pattern during the transition of adjacent spatial consistency regions.
The interpolated channels with a higher number of tap delays (N max = 16, M max = 10) and without applying path clustering (N max = 30, M max = 1) were synthesized to relatively compare the performance of the interpolated channel in Figs.14b -14c.The level of discrepancy in terms of normalized Doppler spread was quite similar between the three interpolated channels when g = 10 cm, as shown in Fig. 14e.A slightly larger discrepancy was found during the period when the path delay was rapidly fluctuating (i.e., 4 -5 s).With a coarse-grained grid size of 100 cm, the normalized Doppler spread discrepancy increased significantly for three cases of interpolated channels, as shown in Fig. 12f.Nevertheless, the normalized Doppler spreads from those three interpolated channels were still time-varying with a pattern similar to that from the reference RT result.The interpolated channel with a higher number of tap delays produced a smaller discrepancy in terms of normalized delay spread in both grid sizes, as clearly seen in Figs.14g -14h.The magnitude of the discrepancy was found to be around 5 -10 for the interpolated channel with N max = 6 in both grid sizes.Increasing N max to 16 reduced the discrepancy of the normalized delay spread, as clearly observed during 0 -2 s when g = 10 cm.The normalized delay spread of the interpolated channel without path clustering certainly agreed with the reference result when g = 10 cm, but with a larger discrepancy when g = 100 cm.The results of the normalized delay spread profiles indicated that the effect of path clustering and association was more significant for the wideband channel characteristic.In case of 3GPP channel, delay and Doppler spread profiles are continuously smooth along MS trajectory owing to the mobility modeling with spatial consistency assumption.However, their profiles could not resemble those produced by the RT channel.Although the overall trend of the 3GPP Doppler spread profile in this case showed a certain degree of similarity to the RT channel in Fig. 14e -14f, different realization may induce totally different characteristic due to its stochastic parameters.In comparison, the delay and Doppler spread of the interpolated channel showed better performance, especially with a smaller grid size.Nevertheless, it should be emphasized that inappropriate path clustering and association may lead to larger discrepancies comparable to that of the 3GPP stochastic counterpart.
Since the discrepancy of the interpolated channel given a perfect path association and no path clustering is negligible small as rigorously investigated in Section IV-A2, it can be summarized that the practical performance of the interpolated channel is highly dependent on the path clustering and association algorithm.Therefore, an intensive investigation of these algorithms is needed to optimize the performance of the interpolated channel.Additionally, non-uniform resolution grid size may be preferable to mitigate discrepancy during the transition region with smaller grid size, while utilizing larger grid size elsewhere to optimize number nodes stored in the emulator database.

V. CONCLUSION
Applying DM in WCE enables its applicability to evaluate the wireless system in a site-specific environment with high accuracy prediction.To realize such a wireless emulator, the grid-based channel modeling technique was proposed.The channel can be synthesized from the pre-computed deterministic path parameters at the reference grid node.The relative residual model of the interpolated SoS-TDL parameters was analytically derived to quantify the discrepancy of the interpolated channel during the interpolation of the SoS-TDL parameters.
Rigorous simulations were performed in the empty hallway to numerically verify the proposed relative residual model at 5.25 GHz.The simulation also confirmed the ability of the proposed channel modeling technique to reproduce the continuous fading pattern in the interpolated SoS-TDL channel response.The performance of the interpolated channel between a fixed BS and a moving MS was validated in terms of delay and Doppler spread.The prediction accuracy of the interpolated channel is confirmed to be highly dependent on the interpolated SoS-TDL parameters.Furthermore, the performance of the interpolated channel is highly dependent on the antenna types, whereas the effect of the frequency bands was negligible.The effect of path clustering and association on interpolated channel performance was investigated with different grid sizes and different numbers of tap delay and rays per tap.KPowermean path clustering and heuristic path-sorting algorithms were applied to the interpolated channel.The artifact caused by these algorithms, such as a long transition between adjacent spatial inconsistency, produced a significantly larger discrepancy than those introduced by the parameter interpolation residual.Another simulation was conducted based on the office space environment to investigate the performance of the interpolated channel in a rich multipath environment.A similar result was observed where a smaller grid size and an increasing number of tap delays can lead to better performance of the interpolated channel.However, despite a smaller residual of interpolation, artifacts introduced by path clustering and association significantly degraded the accuracy of the interpolated channel.Therefore, further analysis of these algorithms and optimal grid size resolution is essential for the practical implementation of the proposed modeling.

FIGURE 1 .
FIGURE 1. Concept comparison of (a) deterministic channel model and (b) grid-based channel modeling technique for WCE.

FIGURE 2 .
FIGURE 2. Grid-based deterministic channel models during the channel emulation at tij -th time instance.

FIGURE 3 .
FIGURE 3. Time-varying characteristic of interpolated SoS-TDL parameters between BS and MS at nth tapped delay and mth associated ray.

FIGURE 4 .
FIGURE 4. Characteristic of interpolated SoS-TDL channel response with local stationary region and transition period.

FIGURE 5 .
FIGURE 5. Geometric comparisons of a single reflected ray trajectory associated with BS-dth node and BS-MS links.
included in the derivation.Let us denote F 1 (θ D d , φ D d ) and F 2 (θ A d , φ A d ) as the magnitude of a single polarized antenna radiation pattern of BS and MS, respectively, where θ {D|A} d and φ {D|A} d

FIGURE 6 .
FIGURE 6. Simulation hallway environment with a room partition (magenta) for rigorous validation of interpolated channel.

FIGURE 7 .
FIGURE 7. Accuracy of the interpolated SoS-TDL parameters ( g = 100 cm) in terms of interpolation residual for (a) delay, (b) Doppler, and (c) weight along A-trajectory.

FIGURE 8 .
FIGURE 8. Time-variant PDP of (a) interpolated ( g = 10 cm) and (b) RT channels, samples of wideband channel in terms of PDP at time instances of (c) 1 s, and (d) 2.5 s, (e) comparison of narrowband fading channels, and its normalized (f) Doppler and (g) delay spreads profiles along A-trajectory.

FIGURE 10 .FIGURE 11 .
FIGURE 10.Statistical performance comparison of interpolated channels at 5.25 and 28 GHz with BS halfwave dipole and horn antennas in terms of absolute discrepancy of (a) normalized Doppler and (b) spread profiles.

FIGURE 12 .
FIGURE 12. Time-variant PDP of (a) RT channel, and the interpolated channel with path clustering and association given Nmax = 3 taps, and Mmax = 5 rays at (a) g = 10 cm, and (b) g = 100 cm grid size, and its normalized (c) Doppler and (d) delay spread profiles along B-trajectory.

FIGURE 13 .
FIGURE 13.Top view of the simulation scenario in an actual office space.

FIGURE 14 .
FIGURE 14. Time-variant PDP of the MS trajectory in the office floor of (a) RT channel, interpolated channels with (b) g = 10 cm, and (c) g = 100 cm, and (d) 3GPP GBSM channel, with the performance of interpolated channel in terms of normalized Doppler spread profiles with (e) g = 10 cm, and (f) g = 100 cm, and normalized delay spread profiles with (g) g = 10 cm and (h) g = 100 cm.