A Time Variant Outdoor-to-Indoor Channel Model for Mobile Radio Based Navigation Applications

Time based positioning with terrestrial mobile radio signals has gained remarkable attention. To develop and validate positioning algorithms under realistic conditions, an accurate knowledge of the propagation channel is significant. However, there is still a lack of outdoor-to-indoor channel models suitable for positioning applications. To be applicable for positioning, the channel model has to fulfill three requirements that have not been accurately considered so far: the non-line-of-sight bias (affecting ranging accuracy), nondiscrete valued channel parameters (affecting algorithm performance), and the evolution of individual multipath components (MPCs) with time (affecting tracking performance). In this paper, an outdoor-to-indoor channel model is proposed based on an extension of the geometry-based stochastic modeling approach to fulfill the requirements. We consider MPCs occurring due to reflection, scattering, and combinations of both. In the model, three different types of MPCs are modeled separately according to their characteristics. Each MPC is represented by a fixed scatterer, which has a fixed position while the receiver antenna is moving. The parameters of the outdoor-to-indoor channel model are extracted from two channel measurement campaigns. The proposed outdoor-to-indoor channel model is capable of accurately simulating the time variant channel. A comparison of the channel model with the channel measurement data is performed by comparing statistics.


Introduction
Global navigation satellite systems (GNSSs) provide accurate positioning as long as line-of-sight (LoS) conditions between satellites and receiver prevail.However, in indoor scenarios, the position accuracy using GNSSs deteriorates due to signal blockage and multipath propagation [1].Enhancing GNSS based positioning by incorporating terrestrial mobile radio signals in indoor environments can significantly improve the position accuracy compared to a GNSS-only solution [2][3][4].In particular, mobile radio signals transmitted from outdoor base stations (BSs) have generated interest due to the availability of mobile radio [5].
The performance of positioning algorithms strongly depends on the multipath propagation channel.To test and develop suitable signal processing algorithms, accurate channel models are needed to provide a realistic representation of wave propagation.Channel models in the context of positioning have to fulfill the following requirements: (i) Positioning algorithms depend on the delay estimation of the LoS path between the transmitter and the receiver.In outdoor-to-indoor scenarios, the LoS path is often blocked.Therefore, in line-of-sight (NLoS) scenarios, the NLoS bias, defined as the delay offset of the first detectable multipath component (MPC) to the delay of the geometrical line-of-sight (GLoS) path, introduces an additional ranging error [6].Therefore, the NLoS bias is of high importance and needs to be taken into account.
(ii) Multipath mitigation methods like [7,8] rely on the estimation of the ray-based channel impulse response.Therefore, simulations with nondiscrete valued channel parameters, for example, MPCs with non-grid-based delays, are necessary.
(iii) Advanced multipath mitigation strategies based on sequential Bayesian filtering like [9,10] are able to estimate time variant MPCs' parameters.Therefore, International Journal of Antennas and Propagation the channel model has to provide a smooth and continuous evolution of channel parameters with time.
Channel models developed for mobile radio signals normally consider communication applications and do not fulfill all requirements in the context of positioning.Standard channel models based on the wide sense stationary uncorrelated scattering (WSSUS) assumption do not fulfill the requirements as the assumption does not describe the channel behavior accurately in many scenarios [11,12].While the receiver moves, the delay of individual MPCs varies and the number of coexisting MPCs changes [13,14] which deteriorates the WSSUS assumption.Because of its advantages to simulate the time variant behavior of the channel, the geometrybased stochastic channel model (GSCM) approach is widely used [15,16] like in WINNER II, QuaDRiGa, COST-2100, or other models like [17].The well known WINNER II channel model [18,19] is dedicated to both link and system level evaluations for communication applications.Its usability for outdoor-to-indoor positioning applications is limited as it supports neither an accurate representation of the NLoS bias nor a continuous evolution of MPCs due to its drop concept.The extended version of the WINNER channel model, that is, QuaDRiGa [20], is capable of simulating time evolution of MPCs.However, the QuaDRiGa model does not support the simulation of the NLoS bias.Similar to the WINNER model, only single-bounce scattering is considered.Meanwhile, the birth-death process of MPCs is simulated by the so-called drop concept, treating the path lives of all MPCs the same.The recently proposed COST-2100 channel model [21][22][23][24] simulates the birth-death process, that is, the visibility, of individual MPCs by a circular region of receiver positions.Using a circular region is not physically inspired and, therefore, less realistic.Furthermore, the COST-2100 channel model considers only MPCs occurring due to scattering.Reflections at large smooth surfaces [25,26], which have distinct time variant characteristics [27], are not considered.The channel model in [17], parameterized by ultra-wideband (UWB) measurement in gas stations, uses a more realistic approach, that is, assigning a directive beam pattern for each MPC to simulate the path's birth-death process as in [28].However, this model considers only paths occurring due to unique scattering.The distinct time variant characteristics for multiple wave interactions of a propagation path are not considered, which may cause inaccuracies in the simulation as compared to reality.
In this contribution, we propose a GSCM considering paths due to reflection, scattering, and combinations of both.A former version of the outdoor-to-indoor channel model has been published in [29], which uses a different concept for virtual scatterers and is less accurate in modeling the distinct propagation phenomena.The proposed GSCM models three different types of MPCs individually.Each MPC is characterized by a so-called fixed scatterer (FS), which is a virtual scatterer whose position does not change during receiver movement.Using the FS allows to simulate a birth-death process based on a directive beam pattern similar to [17].The proposed channel model in this paper is parameterized by statistics from channel measurement data and fulfills the requirements for channel models in the context of positioning.As one of the most essential parameters relevant to positioning applications, the time variant absolute delay (i.e., the time of a signal propagating from the transmitter to the receiver via a certain path) of an MPC is accurately modeled.Thus, the NLoS bias is implicitly included.Time variant channel parameters like the number of MPCs and the amplitude are also modeled based on the spatial movement of the receiver.
An overview of the modeling concept has been published in [30,31].In this contribution, we provide a detailed description, parameterization, and implementation of the channel model.
The structure of the paper is given as follows: in Section 2, the channel measurements are addressed.Section 3 describes the proposed outdoor-to-indoor channel model for mobile radio based positioning, and Section 4 presents the data extraction for the model.A detailed implementation of the proposed channel model is given in Section 5 and a comparison with the measurement data is presented in Section 6. Section 7 concludes the paper.
Throughout the paper, the following notations are used: (i) Vectors and matrices are denoted by lowercase and capital-case bold letters, respectively.
(ii) [⋅]  represents the transpose of a vector or matrix.
(iii) ‖ ⋅ ‖ is the Frobenius norm of a vector or matrix.
(iv)   ( | ) and   ( | ) represent the conditional probability density function (PDF) and probability mass function (PMF) of  given , respectively.  (, ) represents the joint PDF of  and .
(v) x represents the estimate of the parameters .

Channel Measurement Campaign
Two outdoor-to-indoor measurement campaigns have been performed at two different locations at the premises of the German Aerospace Center (DLR), denoted as Scenario-1 and Scenario-2.The MEDAV RUSK-DLR broadband channel sounder is used to sound the channel at carrier frequency 5.2 GHz in the single-input single-output (SISO) manner [32].
In the measurement for Scenario-1, the receiver was located on the ground floor of the building TE01 as shown in Figure 1.A sounding signal was transmitted with a bandwidth of 90 MHz.The receiver antenna was omnidirectional and mounted on a model train moving with a speed of about ‖k‖ = 0.075 m/s as described in [6].For all six transmitter positions, denoted as  1 , . . .,  6 , the model train ran on tracks  1 and  2 from the right to the left side of Figure 1.The transmit antenna was also omnidirectional and positioned at four different locations referred to as  1 to  4 on the rooftop of building TE02 at a height 12 m above ground.Propagation conditions from these locations can be regarded as rather mild in terms of wave propagation loss as the signal propagates directly into the part of the building where the receiver antenna is located.Two other transmitter locations  5 and  6 in front of building TE01 were used with an antenna height 2 m above ground.Propagation conditions from these two locations are considered as severe since the signal absorption is larger due to several traversed indoor walls.The LoS is always absent in the severe condition.Both buildings, TE01 and TE02, can be characterized as standard three-story office buildings of concrete with metalized glass windows.
In the measurement for Scenario-2, the receiver antenna was located on the second (top) floor of the building, which was characterized as a standard three-story office building of concrete with nonmetallic window glass as described in [34].A sounding signal was transmitted with a bandwidth of 120 MHz.The rooms located on the second floor (9 m above ground) are usually used for presentations and contain less furniture than normal office rooms.As shown in Figure 2, the rooms include two presentation rooms, a lobby, and a corridor without windows.Similar to Scenario-1, the receiving antenna was mounted on a model train running along tracks  3 ,  4 , and  5 with a speed of ‖k‖ ≈ 0.1 m/s.The transmit antenna was mounted on an aerial work platform positioned at  7 ,  8 ,  9 , and  10 . 7 and  9 are located at the same azimuth point, but at different heights  1 = 12 m and  2 = 18 m, respectively.Similarly,  8 and  10 are located at the same azimuth point, but at different heights  1 and  2 , respectively.
For all measurements of Scenario-1 and Scenario-2, the same Rubidium standard frequency normal for both transmitter and receiver clocks was used for synchronization.The positions of the transmitter and the receiver were precisely determined using a Leica tachymeter giving a nominal accuracy in the sub-cm domain.
Based on the accurate position information, different closely spaced measured channel impulse responses (CIRs) can be combined to a virtual antenna array (VAA) in postprocessing.Thus, joint estimation of the channel parameters in terms of delay, complex amplitude, and angle of arrival (AoA) is possible and is described in more detail in Section 4.1.

Description of Channel Model for Positioning Applications
3.1.Overview of the Proposed Channel Model.In the wireless channel, a signal propagates from the transmitter to the receiver along certain geometrical paths.Interactions between the signal and physical objects may occur along each path.According to the number of interactions that occurred, each propagation path is defined as single-bounce (only one interaction) or multibounce (more than one interaction).As described in Section 1, the channel model distinguishes two kinds of interactions: reflection and scattering [25,26].
The transmitter location and the environment are assumed to be static; that is, the physical objects where scattering occurs do not change their positions.It is worth noting that a moving physical scatterer may significantly impact the channel coherence as found in vehicle-to-vehicle channels [28,35].In case the scatterer is moving, further modifications should be done.Based on the interactions, the MPCs are classified into three types: (i) M1: reflected MPCs.
(ii) M2: scattered MPCs with scattering as the last interaction.
(iii) M3: scattered MPCs with reflection as the last interaction.
An MPC of type M1 contains only reflections as interactions, either single-bounce or multibounce.MPCs of types M2 and M3 contain at least one scattering interaction [27,30].It can be shown that all types of MPCs can be geometrically represented by a single FS thereby maintaining the delay and the AoA of the MPC to model.An FS is a virtual scatterer whose position does not change while the receiver moves.Hence, using the FS to represent the MPC is advantageous as it does not require additional movement models like in [29,36,37].It is worth noting that defining an FS for a path is similar to the studies in [27] for types M2 and M3, and to the image method [38] for the type M1.To calculate the position of the FS for a specific type of MPC, the concept of an equivalent scatterer (ES) is used in an initialization step.
An ES is a virtual single-bounce scatterer representing the path by a single-bounce interaction, which results in the same geometrical path length and the same AoA as the path to be modeled [12,15,16,27,29].The three different types of MPCs are explained in detail in Section 3.2.The general representation of the modeled CIR can be written as the sum of a finite number of dominant MPCs as with amplitude   (), absolute delay   (), phase   (), and AoA   () of the th MPC varying with respect to time .a(  ()) is the steering vector for the antenna array depending on the AoA   ().The environment and the transmitter are assumed to be stationary and only the receiver moves.
Compared to the carrier frequency, the bandwidth of around 100 MHz is small; the wavelength difference between band edge and central frequency is 1%; therefore, a narrowband assumption to the steering vector is assumed in this paper; meanwhile, the plane wave assumption can be made as objects in the environment are far enough in comparison to the array aperture size; more details can be found in [39].Thus, the time variation of the channel is caused by the movement of the receiver.For convenience, the space variant CIR is described in dependence on the traveled distance  instead of the time , although the term "time variant" is still used.In this paper, the spatial resolution Δ is used to describe the spatial distance between two consecutive CIR snapshots.
The geometrical representation of the proposed outdoorto-indoor channel model is visualized in Figure 3, where the transmitter and the receiver are denoted by Tx and Rx, respectively.The time variant channel parameters   () and   () are deterministically calculated based on the type of MPC and the positions of Rx and FS.The amplitude of a path is a stochastic process parameterized by a Rician amplitude-fading and a Gaussian spectrum (see Section 4.7).Furthermore, a time variant number of MPCs is simulated in a two-step approach: First, the number of disappearing MPCs from the last to the current time instant is determined according to the visible region (VR) as shown in Figure 3.The VR is determined by the position of the FS and an opening angle   , which can be transferred to a path life for a linear receiver trajectory.Second, the number of newly appearing MPCs is stochastically generated.A detailed implementation of the simulation of the time variant channel parameters is given in Section 5.

Reflected MPC (M1).
For a single-bounce MPC of type M1, the corresponding ES is the interaction point of the propagated wave with the surface, that is, the reflecting point as shown in Figure 4.While the receiver changes its position, the ES moves along the wall [29].For a multibounce MPC of type M1, the corresponding ES is a virtual point, which normally is not identical with a physical interaction point.It can be easily shown that if the receiver moves along a line, the ES of a multibounce MPC of type M1 has a linear trajectory as long as the reflections originate from planar surfaces.The FS of an MPC of type M1 is defined as a virtual transmitter whose position does not change.The distance between the FS and the receiver is the same as the propagation path length to be modeled.As visualized in Figure 4 for a single-bounce reflected MPC, the FS locates at the image point of the real transmitter [38].The positions of the FS and the real transmitter are symmetrical with respect to the wall surface on which the reflection occurs.The FS lies on the line determined by the receiver and the ES positions.It can be seen that the distance between the FS for the th MPC and the receiver ‖l  () − l FS, ‖ equals the geometrical path length   () ⋅ , where  is the speed of light.
The transmitter position l  () and the FS position l FS, () are fixed.Therefore, in this paper, the notations l  () and l FS, () are simplified as l  and l FS, , respectively.The position of the receiver is denoted by l  ().

Scattered MPC, Scattering as the Last Interaction (M2).
For this kind of MPC, the last interaction is scattering.It can be seen that the ES of a single-bounce scattered MPC is a physical object where scattering occurs.While the receiver moves, the ES of a single-bounce scattered MPC does not change its position.For a multibounce MPC of type M2, the ES is not a physical object but a virtual point as visualized in Figure 5. Since the last interaction point (i.e., the physical object or "real scatterer" in Figure 5) does not change its position, the propagation time from the transmitter to the last interaction point is fixed.As a result, while the receiver changes its position, the ES moves along an elliptical trajectory [27].As visualized in Figure 5, the FS for an MPC of type M2 is equivalent to the last real scatterer of the propagation path.It lies on the line through the receiver and the ES.

Scattered MPC, Reflection as the Last Interaction (M3).
For this kind of MPC, the last interaction is reflection as depicted in Figure 6.It has been discussed in [27] that the ES may move along a hyperbolical trajectory according to the geometry.As visualized in Figure 6, the FS is defined as the image point of the last scatterer (i.e., the physical object or "real scatterer" in Figure 6) with respect to the reflecting surface.We denote the distance between the transmitter and the FS as  TM and the geometrical path length between the transmitter and the last scatterer (i.e., the "real scatterer" in Figure 6) as  TS .The ES moves along different types of trajectories according to the relation between  TM and  TS : (1) If  TM is smaller than  TS , the ES has an elliptical trajectory while the receiver moves.Thus, this kind of MPC is equivalent to the type M2.
(2) If  TM equals  TS , the ES and the FS are located at the same position.(3) If  TM is larger than  TS , the ES moves along a hyperbolical trajectory.For simplicity, we will only consider case (3) for the type M3.It can be seen that the FS lies on the line determined by the receiver and the ES.

Amplitude of MPCs. The time variant complex amplitude
()   () of the th MPC at the instant  is calculated according to a Rice process [36,40,41] where  0, is the instant when the th MPC appears, ]  the initial phase, and  the wavelength.   is the Rice factor of the th MPC for simulating the stochastic variation of International Journal of Antennas and Propagation the complex amplitude.The initial constant power   ( 0, ) is given as where ( 0, ) is the averaged power,   is the shadowing factor, and the initial phase ]  of the th MPC is randomly drawn from a uniform distribution over [0, 2).In this paper, the shadowing factor is defined as the power difference between the MPC and the average power over all MPCs.Due to different AoAs for unresolvable paths, a spread of the Doppler spectrum has been reported [29,36].Without loss of generality, we define an MPC to consist of   subpaths to simulate the spread in the spectrum.In (2), the phase of each subpath is defined as where V denotes the speed of the receiver.The phase term    of the   th subpath is a random variable that is uniformly distributed over [0, 2).The Doppler frequency offset of the   th subpath to the central Doppler shift is denoted as Δ ,  .

Delay and AoA of MPCs.
In the channel model, we consider the absolute delay   () of the MPC.To exclude the dependence of   () on the delay of the GLoS path   (), the normalized delay is defined as For each MPC, the initial τ ( 0, ) is randomly generated according to a PDF   (τ  ), together with an initial AoA   ( 0, ) drawn from   (  ).
According to the geometry, while the receiver moves, the absolute delays of all types of MPC with either single-bounce or multibounce propagation can be generically presented as where  , is the excess distance defined as The position of the ES at  0, is denoted as l ES, ( 0, ).For type M1, the FS is a virtual transmitter and, therefore,  , = 0.
For type M2, the excess distance equals the major axis of the ellipse.As a result, the excess distance is constant while the receiver moves.A similar conclusion can be made for the type M3 according to the geometrical properties of the hyperbola.Therefore, for all types of MPCs,  , is a constant that is independent of the receiver movement.
While the receiver moves, the AoA   () can be calculated in a straightforward manner according to the positions of the FS and the receiver.

Position of FS.
In many channel models like the COST-2100 channel model and [17], the positions of scatterers are randomly generated.In our approach, we generate the position of the FS for an MPC based on the position of the receiver and a random ES position at the initialization instant  0, .The ES point l ES, ( 0, ) for the MPC is generated based on the transmitter and receiver positions, and the randomly generated values τ ( 0, ) and   ( 0, ).
As visualized in Figures 4, 5, and 6, the FS, the ES, and the receiver position l  () are located on the same line.Therefore, according to geometry, the calculation of the FS position l FS, can be represented generically for all kinds of MPCs by where l  ( 0, ) is the Rx position when the th MPC appears.
The matrix C  reads as where the entity   depends on the type of MPC: where  tr ( 0, ) stands for the distance between the Tx and the Rx,  rs ( 0, ) for the distance between the Rx and the ES, and  ts ( 0, ) for the distance between the BS and the ES at the instant  0, .The so-called FS-ratio used in (10) is defined for type M2 and type M3 individually as International Journal of Antennas and Propagation 7 3.6.Path Life and Visible Region of an MPC.While the receiver moves, the number of coexisting MPCs changes.On the one hand, new MPCs appear randomly.On the other hand, existing MPCs disappear after traveling a certain distance.Path life is a measure of the visibility duration of individual MPCs [36].It is defined as the receiver traveled distance   for which the th MPC is visible to the receiver.A more robust approach is to use the concept of opening angle or directive beam pattern as in [17,42,43], which is capable of simulating a scenario where a receiver is moving along arbitrary trajectories.As visualized in Figure 3, the opening angle   defines the VR of the th MPC.Using the concept of an opening angle requires a static point.As the FS associated with an MPC is fixed in its position, we apply the opening angle to the FS in the proposed channel model.The opening angle can be calculated based on a randomly drawn path life   as described in the appendix.
where Δ is the spatial resolution, that is, the spatial distance of two consecutive CIR snapshots.
To generate (), the terms on the right hand side in (13) must be determined.First, the number of MPCs at the previous instant ( − Δ) is a known number.Second, the number of disappearing MPCs   () is determined according to the beam pattern (i.e., the opening angle) of MPCs.Hence, the problem of generating the number of MPCs at current instant  becomes generating the number of appearing MPCs   () according to the conditional PMF where L() = ( − Δ) −   ().The PMF   (  () | L()) is characterized by the channel measurement campaigns.

Time Variant Channel Parameter Estimation.
A feature of our measurement campaigns is that the position of the receiver antenna for each measured CIR snapshot is precisely known.Therefore, it is possible to form a linear VAA from the time variant measurement as described in [34].For both scenarios, we construct VAAs consisting of  = 8 elements.The interelement spacings between the antennas of the VAA are 0.3 and 0.45 for Scenario-1 and Scenario-2, respectively.Using a VAA, the AoA can be jointly estimated with the delay and the complex amplitude.Due to the nature of linear arrays, the AoA can be only estimated from 0 ∘ to 180 ∘ , that is, angular ambiguity between left and right sides appears.In the measurement scenarios, the receive antenna was at least 1 m and 1.5 m away from other objects, for Scenario-1 and Scenario-2, respectively.Therefore, the far field distance condition is fulfilled.The time variant parameters of the th MPC are defined as a vector To accurately estimate and track the changing parameters Φ  () of individual MPCs, a recursive Bayesian estimator named Kalman enhanced super resolution tracking (KEST) algorithm [44] is used.The KEST algorithm consists of two stages: the inner estimator using the SAGE algorithm [45] and the outer estimator using a Kalman filter.Figure 8 shows an example of estimated time variant CIRs while the receiver was moving.In this paper, the path life is obtained because the receiver traveled the distance for which the path was observed.The NLoS bias   () is calculated as the additional propagation time of the first estimated path compared to the GLoS path.
As described in Section 3.1, the ES defines a virtual single-bounce scatterer.Therefore, at each time instant , the positions of the ESs for all types of MPCs (either with singlebounce or multiple-bounce propagation) can be calculated as in [34].Based on the trajectory of consecutive positions of the ES, the type of MPC and the corresponding position of the FS can be obtained [27,30].For all types of MPCs, the FS is always located at the intersection point between the lines determined by the receiver and the ES for different instants  as visualized in Figures 4, 5, and 6.For the MPCs of type M1, the distance ‖l FS, − l  ()‖ equals τ ⋅ , whereas ‖l FS, − l  ()‖ is less than τ ⋅  for the MPCs of type M3.An example of an estimated MPC of type M2 and the corresponding position of the FS is visualized in Figure 7.The black line indicates the trajectory of the receiver for which the corresponding MPC could be detected.The obtained positions of the ES and the FS are shown by the blue circle and the red solid dot, respectively.The grey line represents the modeled movement of the ES as in [29].It can be seen that the estimated FS lies on the concrete pillar, which may be the object where the wave is scattered.
After identifying the type of each MPC, the positions of the ES and the FS are used to calculate the FS-ratio according to its definition in (11) and (12).The path life of the MPC is determined as the receiver traveled distance when the corresponding path is tracked by the KEST algorithm.At each time instant, the estimated amplitude α () of the th MPC is normalized according to the Friis free space loss as where   is the transmitted power;   and   are the antenna gains for transmitter and receiver, respectively.Thus, the obtained time variant amplitude â = [ α ( 0, ) , α ( 0, + Δ) , . . ., α ( 0, +   )] is used to calculate the Rice factor    according to the method in [46].The Doppler bandwidth is obtained by estimating the Doppler spectrum based on â as in [36].As discussed in [47], to ensure an accurate simulation of the power delay profile (PDP), the delay is generated according to its PDF with a constant mean power.The mean powers for different types of MPCs are calculated based on â .The shadowing factors are obtained as the residuals after averaging out the fast fading and the mean power.Corresponding parameters are listed in Table 1.It can be seen that the power differences between Scenario-1 and Scenario-2 are similar for all three types of MPCs, which may be attributed to additional loss from different building surface materials.

Type of MPCs.
The obtained PMF for the type   of an MPC is listed in Table 1.For MPCs of type M2, the probability of occurrence is higher in Scenario-1 than in Scenario-2.The environmental difference may cause a different frequency of occurrence for MPCs of type M2.In both scenarios, the frequency of occurrence of MPCs with scattering as the last interaction (M2) is significantly higher compared to other types.This might result from the fact that, indoors, many objects like pillars, metal frames, computers, and rough surfaces of objects scatter the signal.Moreover, reflected MPCs of type M1 require large regular surfaces like walls.Therefore, the frequency of occurrence of scatterer MPCs (M2) can be expected to be higher than other types of MPCs.Although the MPCs of types M1 and M3 have smaller probabilities of occurrence, it is important to include these types.First, due to different trajectories of the ES, the calculation of the delay is different as shown in ( 6) and (7).Therefore, to simulate the delay changes accurately, it is important to distinguish the three types of MPCs.Also, path lives for MPCs of types M1 and M3 are larger in average compared to the path life for type M2 (6 m ninety percent of the time).MPCs with long path lives are advantageous for tracking algorithms.For instance, multipath-aided positioning algorithms based on signal tracking like [48] increase their performance if they are able to track paths for longer times.Furthermore, the mean powers of the MPCs of types M1 and M3 are higher (4 dB) than the mean power of the MPCs of type M2 as summarized in Table 1.

Delay and AoA.
In [18,19,47,49], excess delays of MPCs are considered by shifting   () according to the delay of the first detectable path  0 ().An exponential distribution is proposed to model the excess delay.As described in Section 3.4, we consider the normalized delay τ (), which relates to  the absolute delay   () as in (5).Reference [50] evaluates the delay τ () based on ray tracing in an indoor-to-indoor scenario and proposes the Log-normal distribution for the normalized delay.However, a different PDF is proposed based on the measurement data described in this paper.Figure 9(a) visualizes the PDF of the obtained normalized delay τ () for Scenario-1.A clustered structure is visible in Figure 9(a).
The first cluster may be caused by direct signal propagation from the transmitter to the receiver through the window front, while the second cluster may have resulted from a strong double reflection between both buildings.The distance between the two buildings is around 21 m.And the additional path length of the second cluster is around 42 m compared to the first cluster, which translates to a delay of about 140 ns.It is worth mentioning that the similar clustering effect in delay is observed for all transmitter locations.A mixture distribution consisting of two Weibull distributions is used to characterize τ () for Scenario-1 as where  0,1 ,  0,2 ,  1,1 , and  1,2 are the parameters of the Weibull mixture distributions;  0 and  1 are the weights of the corresponding Weibull distribution with  0 +  1 = 1.It is worth noting that the difference between the distributions for τ () for the different types of MPCs is negligible.
Figure 9(b) visualizes the PDF of estimated AoAs for different types of MPCs based on the measurement for Scenario-1.For visualization convenience, the PDF is estimated using a Gaussian kernel estimator [51].The reflected MPCs are mainly caused by the walls of the corridor, which results in a high peak at 16 ∘ in the PDF.AoAs of MPCs of type M2 are almost uniformly distributed.For both scenarios pictured in the channel model, the AoAs of the type M2 MPCs are modeled as uniformly distributed random variables.For the other two types of MPCs, a Gaussian mixture model is used to model the AoA.Table 1 shows the parameters of the Gaussian mixture for both scenarios.  with  ∈ {0, 1, 2} represents the weight for the Gaussian component, and   and   with  ∈ {0, 1, 2} represent the mean and the standard deviation of the th component, respectively.Furthermore, there is no significant dependency between the AoA and the delay.

FS-Ratio.
For MPCs of types M2 and M3, the FS-ratio is used to determine the position of the FS as described in Section 3.5.Based on the position of the ES at instant  0, , the FS-ratios  M2, and  M3, are calculated according to (11) and (12), respectively.
The statistics of the obtained FS-ratio  M2, from the measurement for Scenario-1 are visualized in Figure 10, where the Log-normal distribution is found to fit well.It is worth mentioning that most of the values for the ratio are small.Based on the estimation results, most of the indoor scattering is caused by concrete pillars, metal frames, or internal metallic objects like computers.The statistics of the FSratio  M3, obtained from the measurement for Scenario-1 are visualized by the cumulative distribution function (CDF) in Figure 10, where the Weibull distribution is found to fit well.The  values based on the Kolmogorov-Smirnov (KS) test [33] are 0.47 and 0.76 for types M2 and M3, respectively.

Path Life.
According to the definition in Section 3.6, the path life   of the th MPC is estimated as the traveled distance over which the MPC was detected by the KEST algorithm.Figure 11 visualizes the CDF of the path life of MPCs measured in Scenario-1.A significant dependence between the path life and the type of MPC is noticed.The path life of MPCs of type M1 has a higher probability for larger values compared with the path life of MPCs of type M2.A similar conclusion can be made for Scenario-2 based on the results.To compare the path life to literature, we first introduce the PDF of the path life   (  ), which can be represented by the marginalization over the MPC type   ∈ {M1, M2, M3} Since the type M2 occurs most often among all three types of MPCs according to Table 1, the probability of the overall path life is dominated by the term   (  | M2) ⋅   (M2).Hence, the averaged overall path life is short similar to [14,36].However, as described in Section 4.2, MPCs with a long path life significantly influence tracking algorithms' performances.Therefore, MPCs of types M1 and M3 need to be taken into account.
It has been reported in [14] that the path life is approximately exponentially distributed.However, a Log-normal distribution fits best to path lives obtained from our measurement data.The difference in the model of the PDF may be caused by a different spatial resolution (i.e., the increment Δ in the receiver traveled distance ) during the measurement.
If the spatial resolution is large, MPCs with smaller path lives cannot be estimated.Furthermore, different measurement environments may result in different statistics.Nevertheless, a similar shape of the PDF for the path life has been reported in [36].It is worth mentioning that the receiver moves along straight lines during the measurement process.As a result, the obtained path life characterizes the visibility duration of the MPCs while the receiver moves along a linear trajectory.

Number of MPCs.
The maximum number of appearing MPCs and/or disappearing MPCs obtained is two for Scenario-1 when Δ = 3 mm and one for Scenario-2 when Δ = 1 mm.These small values are mainly caused by the high spatial resolution in the estimation of channel parameters.For larger Δ, the change in the number of MPCs is larger.Figure 12 visualizes the obtained conditional PMF   (  () | L()) for the number of appearing MPCs obtained from the measurement with Δ = 3 mm, where L() = ( − Δ) −   ().When L() reaches its minimum value, that is, L() = 4, the probability of adding one MPC is increased to one.When L() reaches its maximum value, that is, L() = 26, the probability of adding one MPC is decreased to zero.The values of conditional PMFs,   (  () = 0 | L()) and   (  () = 1 | L()), obtained from the measurement data are listed in Tables 2 and 3,  respectively.The corresponding conditional PMF for the observed maximum number of appearing paths for a spatial resolution of Δ = 3 mm is given as (20) 4.7.Rice Factor and Doppler Bandwidth.As described in Section 3.3, the time variant amplitude of an individual MPC experiences small scale fading due to the limited bandwidth.Figure 13 visualizes the estimated time variant amplitude of an MPC that can be well modeled by a Rice distribution.Figure 14 visualizes the CDFs of    for both Scenario-1 and Scenario-2.The estimated Rice factors fit well to a Gaussian distribution.Another finding is that the Rice factor depends only slightly on the scenario.There is no significant dependency between the Rice factor and the types of MPCs.The statistical values of the Rice factor for Scenario-1 are almost the same as those for Scenario-2.As another important factor to characterize the fast fading process, the estimated 3 dB bandwidth Δ 3 for Scenario-1 is depicted in Figure 15.A Weibull distribution is used to model the distribution of Δ 3 .Due to the tail of the Weibull distribution, a truncated distribution according to the maximum Doppler is applied to limit the generated Doppler spectrum similar as in [36].It is worth mentioning that the Doppler bandwidth can be translated into the angular spread of the MPC like in the WINNER channel model.In this work, we focus on the Doppler bandwidth approach as in [36].

Implementation of the Channel Model
The flow chart of the channel model implementation is depicted in Figure 16.At the beginning of the simulation, the positions of the Tx and the Rx are given by the user.An averaged number of MPCs (0) is used to initialize the CIR.Modeling the time variant CIR consists of two steps: initialization and update.The initialization step for each new MPC is indicated by a separate block, "Initialization of MPCs," in the lower right corner of Figure 16.Each MPC is initialized according to statistics obtained from the measurement data.The blocks within "Updating of MPCs" illustrate the update steps of the time variant CIR for each instant.Based on the initialization, the time variant CIR is deterministically updated.In the following parts, each block is described in detail.

Initialization of MPCs.
In the initialization step, each new MPC is stochastically parameterized.The parameters include the type of MPC, the normalized delay, the AoA, the complex amplitude, the fading parameters, the position of the FS, and the opening angle.The initialization of the parameters is implemented using the following steps.

Type of MPCs.
The type of the th MPC   ∈ {M1, M2, M3} is randomly chosen according to Table 1.In the following study, the type of MPC influences the initialization of several parameters.For instance, the determination of FS position, the calculation of excess distance, and the path life depend on the type of MPC.

Amplitude Parameters of MPCs.
The mean power ( 0, ) without shadowing is assigned according to Table 1.Furthermore, an MPC shadowing factor   in dB is randomly drawn from a Gaussian distribution N(0,  2  ) based on the measurement data.Thus, the amplitude   ( 0, ) in ( 2) is initialized as The initial phase ]  of the th MPC is randomly drawn from a uniform distribution over [0, 2).The Rice factor    of the th MPC is randomly generated from a Gaussian distribution.  = 20 subpaths are used to model the fast fading process.A Doppler offset Δ ,  is randomly drawn from a Gaussian distribution and assigned to the   th subpath.The 3 dB Doppler bandwidth Δ 3 is drawn from a Weibull distribution as described in Section 4.7.  ( 0, ), the normalized delay τ ( 0, ), and the AoA   ( 0, ), the ES position l ES, ( 0, ) is uniquely determined and can be calculated as in [34].It is worth noting that l ES, () is assumed to be on the same horizontal plane as the receiver.

FS Position
The FS-ratios  M2, and  M3, are needed to determine the FS positions of MPCs of types M2 and M3, respectively.The ratio  M2, is randomly drawn from a Log-normal distribution whereas the ratio  M3, is drawn from a Weibull distribution as described in Section 4.4.
As described in Section 3.5, the position of the FS l FS, is calculated according to (8) based on the FS-ratio and the positions of the ES, the Tx, and the Rx at  0, .

Excess Distance.
The FS position does not change while the Rx moves.Therefore, the excess distance from the Tx to the FS is constant.The excess distance  , is initialized according to (7).
5.1.6.Opening Angle.According to the type of MPC, the path life is randomly drawn from the Log-normal distribution with   (  |   ).The generated path life is transferred to the opening angle   according to (A.3).
For each instant  in the simulation, the steps above are applied to each appearing MPC.After the initialization of the MPC, the time variant changes are deterministically calculated according to the Rx movement.A summary of the parameters needed for the initialization of individual MPCs is listed in Table 1.

Update Position of Rx.
As a first step, the position of the Rx is updated.It is worth noting that the calculation of the ES position is only required in the initialization step of an MPC to determine the position of the FS.Although the ES position changes while the Rx moves, a new calculation of the ES position at each instant is not needed.

Remove MPC.
Based on the new position of the Rx, each MPC is checked to see whether the Rx is in the VR.If the updated Rx position is outside the VR of an MPC, the MPC is removed from the CIR.When the Rx moves along a straight line, removing the MPC according to the opening angle is equivalent to removing the MPC according to the path life as defined in Section 3.6.

5.2.3.
Update Delay, AoA, and Amplitude.For each of the remaining MPCs, the delay   (), the AoA   (), and the amplitude   () are updated.According to the geometry at instant , the MPC delay is calculated based on ( 6) and (7) for the different types of MPCs.The AoA is calculated based on the FS position, the Rx position, and the Rx moving direction assuming a coordinate system that is aligned with the moving direction.The complex amplitude is calculated according to (2).

Add New MPC(s). If the number of appearing MPCs
() is larger than 0, one or more MPCs are added to the CIR.The initialization step described in Section 5.1 is performed for each added MPC.

Simulation Example and Comparison
To demonstrate the proposed time variant channel model, we present simulation examples for Scenario-1 in the following.An example of generated CIRs in an NLoS scenario is visualized in Figure 17.The Rx moves over a total distance of 12 m.The normalized powers of MPCs are depicted in colors and range between −70 dB and −35 dB.To visualize the NLoS bias, the time variant delay of the GLoS path is shown by the grey line.Variations in delay and amplitude can be clearly seen from the figure.Some MPCs of types M1 and M3 are marked.It is worth mentioning that the delays are independent of the type of the MPCs.For instance, there are also few MPCs of type M1 with shorter delays but not marked in the figure due to visualization reason.It can be seen that MPCs of types M1 and M3 have relatively long path lives and high power levels.For algorithms based on path tracking like [9,10], long-existing MPCs are of high importance.Therefore, although the probability of the types M1 and M3 MPCs is small, it is necessary to consider these two types of MPCs in the channel model.
The proposed time variant channel model does not rely on a room layout; it relies on the GSCM approach.To verify the modeling approach, we compare the statistics obtained by channel model simulations with the statistics extracted from measurement data.Then a comparison is performed based on the following statistics: the NLoS bias, the PDP, and the number of MPCs.Four independent runs are simulated by the channel model.In each run, the Rx travels 21 m with a speed of 0.03 m/s and the Tx is located at different positions.The simulation sampling rate is 10 Hz resulting in a spatial resolution of Δ = 3 mm.The carrier frequency in the simulation is 5.2 GHz and the simulated CIR is band-limited to 90 MHz.
Figure 18 visualizes the PDF for each delay bin of the PDP obtained from the measurement data and through simulations using the channel model.A normalization is performed in power with respect to free space loss and in delay with respect to delay of the GLoS path.Comparing both results, a similar shape can be observed.Particularly, the Figure 19 visualizes the CDF of   () ⋅  (i.e, corresponding to the first MPC) obtained from the simulation in comparison with the CDF obtained from the measurement data.There is a good match between both CDF curves.To evaluate the spatial characterization of the NLoS bias, the normalized autocovariances of   () over traveled distance obtained from simulation and measurement data are visualized in Figure 20.Considering the decorrelation level of 1/ = 0.37, similar decorrelation distances of the NLoS bias (i.e., around 0.5 m) are obtained for both cases.Another comparison is conducted in terms of the number of MPCs.As visualized in Figure 21, the PMF of the overall number of MPCs obtained from the simulation fits well with the PMF obtained from the measurement data.Thus, we come to the conclusion that the proposed channel model simulates the propagation conditions realistically and fulfills the requirements in the context of positioning.

Conclusion
The contribution of this paper is an outdoor-to-indoor channel model fulfilling the requirements for simulations of mobile radio based positioning algorithms.The non-line-ofsight bias and the continuous evolution of multipath components with time are taken into account.The channel model is based on an extension of the geometry-stochastic modeling  the transmitter, the receiver, and the fixed scatterer, the time variant channel impulse response is deterministically calculated.
The statistical parameters provided in this paper are obtained from two different measurement campaigns for an outdoor located transmitter to a receiver positioned in an indoor office.For different indoor environments like shopping malls, corresponding parameter settings may be different.However, the proposed channel modeling approach is generic and can be applied to other buildings.

Figure 3 :
Figure 3: Visualization of the proposed channel model.For representation convenience, only the MPCs of types M1 and M2 are visualized.

Figure 4 :Figure 5 :
Figure 4: MPC of type M1: example of the ES located at the position of the reflection on the surface.The ES moves along a linear trajectory while the FS is fixed during receiver movement.

Figure 6 :
Figure 6: MPC of type M3: example of an ES moving along a hyperbolical trajectory while the FS is fixed during the receiver movement.

InternationalFigure 7 :Figure 8 :
Figure 7: An example of the estimated positions of ES and FS for Scenario-1.
ϱ M3,l based on measurement ϱ M2,l based on measurement ϱ M2,l Log-normal fit for ϱ M3,l Weibull fit for

Figure 11 :
Figure 11: CDF of the path life   for Scenario-1.The  values based on the KS test are 0.48 for type M1 MPC, 0.74 for M2, and 0.73 for M3.

Figure 13 :
Figure 13: CDF of amplitude for an estimated MPC.

Figure 14 :Figure 15 :
Figure 14: CDF of estimated Rice factors.The  values obtained by the KS test are 0.3 and 0.52 for Scenario-1 and Scenario-2, respectively.

Figure 17 :
Figure 17: Example of generated time variant CIRs.Some examples of modeled MPCs of types M1 and M3 are marked.It is worth mentioning that the delays are independent of the type of the MPCs.For instance, there are also few MPCs of type M1 with shorter delays but not marked in the figure due to visualization reason.

Figure 18 :
Figure18: PDF of the normalized power on each delay bin for channel measurement and channel model.The PDP is normalized in power to the free space loss and shifted in delay to the GLoS.

Figure 19 :Figure 20 :Figure 21 :
Figure 19: NLoS bias comparison between measurement based estimation and channel model based simulation for Scenario-1.
3.7.Number of MPCs.The number of MPCs () at instant  is represented as the number of MPCs at the previous instant  − Δ minus the number of disappearing MPCs   () plus the number of appearing MPCs   () as

Table 1 :
Summary of the parameters used in the channel model.
(a) PDF of the normalized delay τ () for Scenario-1 (b) PDF of the AoA θ () for Scenario-1