Modeling the Random Orientation of Mobile Devices: Measurement, Analysis and LiFi Use Case

Light-fidelity (LiFi) is a networked optical wireless communication (OWC) solution for high-speed indoor connectivity for fixed and mobile optical communications. Unlike conventional radio frequency wireless systems, the OWC channel is not isotropic, meaning that the device orientation affects the channel gain significantly, particularly for mobile users. However, due to the lack of a proper model for device orientation, many studies have assumed that the receiver is vertically upward and fixed. In this paper, a novel model for device orientation based on experimental measurements of forty participants has been proposed. It is shown that the probability density function (PDF) of the polar angle can be modeled either based on a Laplace (for static users) or a Gaussian (for mobile users) distribution. In addition, a closed-form expression is obtained for the PDF of the cosine of the incidence angle based on which line-of-sight (LOS) channel gain is described in OWC channels. An approximation of this PDF based on the truncated Laplace is proposed and the accuracy of this approximation is confirmed by the Kolmogorov-Smirnov distance (KSD). Moreover, the statistics of the LOS channel gain are calculated and the random orientation of a user equipment (UE) is modeled as a random process. The influence of the random orientation on signal-to-noise-ratio (SNR) performance of OWC systems has been evaluated. Finally, an orientation-based random waypoint (ORWP) mobility model is proposed by considering the random orientation of the UE during the user's movement. The performance of ORWP is assessed on the handover rate and it is shown that it is important to take the random orientation into account.


I. Introduction
T HE increasing request for wireless data, which is expected to be 49 exabytes per month by 2021 [1], motivates both academia and industry to invest in alternative solutions. These include mmWave, massive multiple-input multiple-output (MIMO), free space optical communication and Light-Fidelity (LiFi) to support the data traffic growth and next-generation high-speed wireless communication systems. Among these technologies, LiFi is a novel bidirectional, high-speed and fully networked wireless communication technology. LiFi uses visible light as the propagation medium in the downlink for the purposes of illumination and communication. It can use infrared in the uplink so that the illumination constraint of a room remains unaffected, and also to avoid interference with the visible light in the downlink [2], [3]. LiFi offers a number of important benefits that have made it favorable for recent and future research.
These include the very large, unregulated bandwidth available in the visible light spectrum (more than 10 3 times greater than the whole RF spectrum), high energy efficiency [4], the rather straightforward deployment which uses off-the-shelf light emitting diodes (LED) and photodiode (PD) devices at the transmitter and receiver ends respectively, and enhanced security as light does not penetrate through opaque objects [5]. One of the key shortcomings of the current research literature on LiFi is the lack of appropriate statistics of device orientation and rotation modeling for system design and handover management purposes.
Smartphones are the most significant and indispensable part of the wireless network generating about 86% of the mobile data traffic [1]. LiFi as part of future 5G can handle this immense data traffic thanks to future LiFi-enabled smartphones. Generally, users tend to work with their smartphones in a comfortable manner which is not necessarily vertically upward. Smartphones are equipped with a gyroscope that can measure the device orientation. This orientation information can be fedback to the access point (AP) via limited-feedback methods [3], [6], [7]. Then, the AP can use the orientation information for resource allocation or handover management. Many previous studies assumed that the receiver is vertically upward and fixed for simplicity purposes and also due to the lack of a proper model for device orientation in LiFi networks. However, there are few studies that have considered the effect of device orientation in their analysis [8]- [19]. Nevertheless, none of these studies have considered the actual statistics of device orientation and have mainly assumed uniform or Gaussian distribution with hypothetical moments for device orientation. In this study, based on practical measurements from forty participants, a new statistical model for device orientation is proposed.

A. Literature Review and Motivation
In [8], the authors consider three standard angles similar to those used in mobile devices to model the device orientation; namely yaw, pitch and roll. Based on this model, the effect of arbitrary orientation on users' throughput and network load balancing is investigated. The problem of handover due to device rotation for downlink in an indoor optical attocell network was first proposed in [9]. The handover probability has been obtained for both sitting and mobile users while considering device orientation. In [10], the handover probability in hybrid LiFi/RF-based networks is evaluated assuming randomly-oriented devices. The impact of the receiver's tilted angle on the channel capacity of visible light communications (VLCs) is investigated in [11]. The lower and upper bounds of the channel capacity for the VLC arXiv:1805.07999v4 [cs.IT] 28 Sep 2018 are presented and by considering an optimization problem the channel capacity has been improved by tilting the receiver plane properly. In [12], a theoretical expression of the bit error ratio (BER) for input-dependent noise of VLC using on-off keying has been derived. Then, a convex optimization problem is formulated based on the derived BER expression to minimize the BER performance by tilting the receiver plane. The impact of device orientation on BER performance of DC biased optical orthogonal frequency division multiplexing (DCO-OFDM) has been evaluated in [13]. A closed form approximation for BER of randomly-orientated UEs is derived.
The effect of device orientation on positioning has been investigated in several studies [14]- [17]. By using the accelerometer sensor, the positioning technique proposed in [14] can be used for any arbitrary device orientation. Downlink and uplink indoor positioning techniques based on VLC while considering the tilting of the device have been developed in [15] and [16], respectively. It is shown that the tilting angle can affect the positioning error significantly. Therefore, device orientation should be considered in the positioning analysis [17]. The signal-to-noise-ratio (SNR) and spectral efficiency improvement of OFDM signals for indoor visible light communication by optimally tilting the receiver plane is proposed in [18]. In [19], the effect of random orientation on the line-of-sight (LOS) channel gain for a randomly located user equipment (UE) has been investigated. The statistical distribution of the channel gain has been derived for a single LED and extended to a scenario with double LEDs. We note that none of these studies are supported by any experimental data. A measurement of the random orientations of mobile devices has been made in [20], but the authors only measure the statistics of the pace of change of the device orientation. Their results, therefore, do not describe the statistical model of the randomly-oriented devices in general. Our recent work in [21] reports some initial results based on the experimental data from 40 participants. The effect of randomly-oriented devices has been studied and a statistical model of the LiFi channel considering the random orientation is proposed.

B. Contributions and Outcomes
The lack of a proper model for device orientation along with its analysis was a motivation to perform a set of experimental measurements, in which participants use their smartphones in landscape or portrait mode in both sitting and walking conditions. It is worth mentioning that the initial statistical model proposed for device orientation can be also used in mmWave communications as it depends on the user behavior and not the access technology. The main contributions of this paper are summarized as follows: • Performing a set of practical measurements, proposing a new model for device orientation based on the measurement results and deriving the probability density function (PDF) and statistics of the model.
• Deriving the PDF and statistics of both LOS channel gain and received SNR.
• Proposing an orientation-based random waypoint (ORWP) mobility model that considers the random orientation of the UE while the user is moving. The ORWP mobility model can also be used in mmWave networks due to the importance of the angle of arrival in those systems. Notations: f X and F X denote the PDF and the cumulative distribution function (CDF) of the random variable (RV) X, respectively. E [·] denotes the expected value over the RV inside the argument, and E X [·] denotes the expected value with respect to the RV, X. [] T stands for transpose operator. The inner product operator is shown by · and · expresses the norm of a vector. Further, tan −1 (y/x) is the four-quadrant inverse tangent.

II. Rotation Geometry
According to Euler's rotation theorem [22], any rotation in R 3 space can be uniquely achieved by composing three elemental rotations, i.e., the rotations about the axes of a coordinate system. Depending on whether the device (local) or Earth (global) coordinate system is chosen, there are two types of rotations. Intrinsic rotation corresponds to a rotation about the device coordinate system and extrinsic rotation, which conforms to a rotation about the Earth coordinate system. Throughout this paper, we will show the device and Earth coordinate system with xyz and XYZ, respectively. The Earth and device coordinates are shown in Fig. 1(a).
Since the matrix multiplication is not commutative, the elemental rotation order matters. Therefore, there are six possible choices of rotation axes for proper Euler angles. World wide web consortium (W3C) has chosen the intrinsic rotation orders of (z → x → y ) as standard for device orientation where x y z and x y z are respectively the device coordinate systems after rotation about the z-axis and subsequently after rotation about x -axis [23]. These rotation orders will be also adopted in this study. According to W3C specification, the angles α ∈ [0, 360), β ∈ [−180, 180) and γ ∈ [−90, 90) correspond to the rotation about z, x and y axes, respectively. The elemental rotation angles, α, β and γ are called yaw, pitch and roll, respectively. The rotation about the axes are depicted in Fig. 1 Now we derive the concatenated rotation matrix with respect to the Earth coordinate system. Let n u = [n 1 , n 2 , n 3 ] T and n u = [n 1 , n 2 , n 3 ] T be the normal vectors before and after device rotation. Based on the Euler's theorem, we have n u = R α R β R γ n u , where R α , R β , and R γ denote the rotation matrices with respect to the z, x and y axes, respectively. Let's R = R α R β R γ be the ultimate rotation matrix. Assuming that the Earth and device coordinate systems initially coincide and n u = [0, 0, 1] T , the corresponding rotated normal vector, after applying the rotation matrices R α , R β , and R γ , is given by (1) at the top of this page. Note that since the initial device coordinates are aligned with the Earth coordinates, the rotation matrix, R, can be represented in the Earth coordinates, XYZ. As can be seen, the rotated normal vector is a function of the three elemental Euler's angles.
The rotated normal vector, n u , can be represented in the spherical coordinate system (corresponding to XYZ) with the polar angle θ and azimuth angle ω. Thus, θ is the angle between n u and the positive direction of the Z-axis, while ω denotes the angle between the projection of n u in the XY-plane and positive direction of the X-axis. Note that cos θ = n u ·Ẑ/ n u , whereẐ = [0, 0, 1] T is the unit vector of the Z-axis. Then, from (1), the polar angle θ can be obtained as: It is clear from (2) that the polar angle only depends on the pitch and roll rotation angles which are further associated with the movements of human's wrists. Referring to Fig. 3, the azimuth angle, ω, can be obtained as follows: ω = tan −1 n 2 n 1 = tan −1 sin α sin γ − cos α cos γ sin β cos γ sin α sin β + cos α sin γ From a mobility point of view, let's define the angle Ω which represents the angle between the direction of movement and the X-axis in the Earth coordinate system. The angle Ω can be described as the angle between the projection of the y-axis on the XY-plane and X-axis when the user is working with the cellphone in portrait mode, or the angle between the projection of the x-axis on the XY-plane and X-axis when the user is working with the cellphone in landscape mode. Thus, whereŷ andx are the unit vectors of y and x axes in the device coordinate system andX is the unit vector of X-axis in the Earth coordinate system. When the user is working with the cellphone in portrait mode, the angle between the y-axis and North is defined as α, and Ω is specified based on the angle between the y-axis and East. Hence, the relationship between Ω and α can be expressed as: Since the difference between portrait mode and landscape mode is just π/2 which follows the right-hand rule, for landscape mode, we have:

III. Mobile Device Orientation Statistics
An experiment is designed to study mobile users behavior and to develop a statistical model for the orientation of mobile devices that act as the receiver for wireless communication systems. During the experiment, 40 participants were asked to use their cellphones normally that create 222 datasets for orientation. They were asked to use the cellphone in both portrait and landscape modes for one minute. The orientation data is measured for both sitting and mobile users. In the experimental measurement, the application Physics Toolbox Sensor Suite has been used as it can provide instantaneous rotation angles, α, β and γ [24]. This application can be running in the background while the participants can perform activities that require data connection, e.g., browsing or watching streaming videos. Below is a summary of the experimental setup: • Activities while sitting: 1) Browsing twice in portrait mode, 2) Watching streaming videos twice in landscape mode, • Activities while walking following a certain path: 1) Browsing in portrait mode, 2) Watching streaming videos in landscape mode. The path that the participants took was a straight corridor with dimensions of 40 m ×1.5 m. The participants were asked to walk down the corridor once. We note that the shape of the test area should not affect the experimental results and the model for the elevation angle as it mostly depends on the posture and physical attributes of typical users rather than the environment. This has been confirmed with sets of uncontrolled data collections from participants using their device in different environments [25].

Measurement and Analysis
The statistics of the azimuth and polar angles can be measured based on the collected experimental data set. In order to evaluate the similarity of measurement data with a particular distribution, we consider Kolmogorov-Smirnov distance (KSD), skewness and kurtosis. The two-sample Kolmogorov-Smirnov distance (KSD) is the maximum absolute distance between the CDFs of two data vectors, which can be obtained as [26], whereF 1 (x) andF 2 (x) are the CDFs of the first and second data vectors, respectively. Smaller values of KSD correspond cos γ sin α sin β + cos α sin γ sin α sin γ − cos α cos γ sin β cos β cos γ to more similarity of the distributions. Skewness is described as a measure of the symmetry or asymmetry of a probability distribution. A perfectly symmetrical distribution will have a skewness of 0. For example, the Gaussian and Laplace distributions both have a skewness of 0. Mathematically, the skewness of a RV X is defined as [27]: where µ and σ are the mean and variance of the distribution X. Kurtosis is a measure of the tailedness of a probability distribution which is given as: The kurtosis of Laplace, Gaussian and Uniform distributions are 6, 3 and 1.8, respectively. The sample PDF of the azimuth angle for sitting and walking activities are represented in Fig. 2-(a) and Fig. 2-(b), respectively. As can be seen, the azimuth angle closely follows a uniform distribution, i.e., ω ∼ U[−π, π), with the skewness of −0.03, kurtosis of 1.68 and Kolmogorov-Smirnov distance (KSD) of 0.034 for sitting activities. Also it is shown that walking activities have a skewness of −0.0045, kurtosis of 1.85 and KSD of 0.019. Note that the peaks in the PDFs (more visible for sitting) are due to the limited number of users, therefore, the azimuth angle is quite correlated for each sitting user along with the direction of the chair they used.
For the polar angle, θ, we model it with Laplace and Gaussian distributions taking into account the first and second laws of error. According to the first law of error proposed by Laplace in 1774, the frequency of an error can be represented as an exponential function of the numerical magnitude of the error, regardless of the sign as f X (x) = k 2 e −k|x−x q | , where x and x q are the measured data and actual value, respectively and k is a constant [28]. The second law of error was also proposed by Laplace four years later and it states that the frequency of error is an exponential function of the square of the error, with σ as the variance of measured data. This is also called the normal distribution or Gauss law of error [28].
The PDF of the polar angle, θ, is estimated based on a set of almost uncorrelated data samples taken from the measurement data. Since the acquired data from the application are unevenly-spaced in time, we first generate a set of sufficiently separated data in time to ensure that samples taken from an individual user are almost uncorrelated. To define the proper separation between uncorrelated samples, we need to calculate the coherence lag between the samples that can be acquired by the autocorrelation function (ACF). However, due to uneven distribution of samples in time, the classical Fourier analysis is no longer valid. In fact, interpolation can be used to generate evenly spaced data samples but this would affect the autocorrelation by removing higher frequencies [29], [30]. Another solution would be the use of least-squares spectral analysis (LSSA) [31], [32]. Based on the LSSA method, we first calculate the power spectral density (PSD) of the unevenly spaced samples and then, the ACF can be estimated by calculating the inverse Fourier transform of the PSD. This technique has been explained in detail in our recent paper [33] that reports the coherence time of the polar angle for sitting and walking activities as 373 ms and 134 ms, respectively.
The histogram of θ using (2), where β and γ are collected from the experimental measurements is shown in  Table I, the kurtosis of the polar angle, θ, for sitting activities is 6.36 which is closer to the kurtosis of the Laplace distribution than the Gaussian distribution, whereas for walking activities the kurtosis of empirical data is 3.77 which confirms a greater similarity to the Gaussian distribution compared to the Laplace distribution. Hence, for sitting activities, the Laplace distribution closely matches the distribution of the experimental measurements in comparison with the Gaussian distribution when considering both KSD and kurtosis metrics. For walking activities, however, the Gaussian distribution matches the experimental data more closely. For the rest of the paper up to section VI, we mainly focus on the sitting activities and therefore the derivations are based on the Laplace distribution for the polar angle. However, one can readily apply a similar methodology to obtain the analytical results for the Gaussian distribution. Note that, in section VI, we use the Gaussian model because the UE is assumed to be mobile. The other important observation here is in regards to the moments of the measured data. The mean of the sitting and walking activities are about 41 • and 30 • , respectively and both with a standard deviation of less than 9 • . Based on the experimental results, the PDF of f θ (θ) can be properly fitted with the truncated Laplace distribution. Mathematically, θ ∈ [0, π], however, as shown by the experimental measurements, the samples of the angle θ are restricted to the range [0, π 2 ]. Therefore, the PDF of θ can be conveniently denoted as: where b θ = σ 2 L /2 > 0. The mean and scale parameters are set to the values from Table I Note that with the parameters for µ θ and b θ given in Table I, we have G( π 2 ) ≈ 1 and G(0) ≈ 0. Thus, (10) can be simplified as:f where π 2 0f θ (θ) dθ ≈ 1. The CDF of θ is also given as: IV. Analysis of Random Orientation Effect on Channel Gain The downlink geometry of a generic optical wireless communication system is shown in Fig. 3. The locations of a UE and an AP are denoted by position vectors (x u , y u , z u ) and (x a , y a , z a ), respectively. The vertical distance of the AP and the UE is denoted by h, where h = z a − z u , and the Euclidean distance of the AP and the UE is denoted by d. The half-power semiangle of the LED is denoted as φ 1/2 and the field-of-view (FOV) of the PD is Ψ c . The angle of incidence at the PD is denoted by ψ. The LOS channel gain of the optical wireless channel is given as [34]: where m =− ln(2)/ln cos φ 1/2 , rect( ψ Ψ c ) = 1 for 0 ≤ ψ ≤ Ψ c and 0 otherwise; A is the physical area of the PD. From Fig. 3 and (13), it is clear that for a particular UE position, the statistics of the LOS highly depend on the device orientation. In order to develop the statistics of the LOS channel gain, the statistics of cos ψ needs to be first determined based on the statistics of the device orientation discussed in section III. We note that the radiance angle, φ, is not affected by the random orientation and we always have cos φ = −n a · d. It should be mentioned that in this paper, we focus on the effect of UE's random orientation and do not consider the possible movements in the user's hands that may change the distance d. Note that such random changes in the distance would be small compared to the average distance between the AP and the UE. Nevertheless, if an analysis of the random movement of the UE is of interest, our results can be readily used as the conditional statistics of the channel gain for a given location of the UE to develop the joint statistics assuming that the statistics of the random movement is available. This is similar to how we analyze the effect of random orientation in a mobile scenario in section VI where the random movement is modeled based on the random waypoint model.

A. PDF of cos ψ
Throughout this paper, we calculate the PDF of cos ψ conditioned on ω. Let the vector d be the distance vector from the UE to the AP as shown in Fig. 3. Noting that the vector n u denotes the normal vector of the UE receiver after rotation, the cosine of the incidence angle, ψ, can be obtained as: where ω is given in (3) and depends on α, β and γ. According to the experimental data, the standard deviation of γ in portrait mode for sitting and walking activities are 0.072 and 0.113 radian, respectively. The standard deviation of β in landscape mode for sitting and walking activities are 0.044 and 0.091 radian, respectively. These small values of standard deviation confirm that when a user is working with the cellphone in portrait mode, the roll angle, γ, is almost zero. Whereas, when the user works with the cellphone in landscape mode, the pitch angle, β, is close to zero. Hence, substituting γ = 0 in (3), it can be simplified asω = tan −1 − cos α sin α . Therefore, for portrait mode, we have: Similarly for landscape mode, we have: Comparing Ω given in (5) and (6), with the above equations forω, one can easily deduce Ω =ω + π. Note that, based on the measurements, the mean absolute error of approximating the angles ω asω is about 0.09 and 0.14 radian respectively for sitting and mobile users, which confirms a relatively good accuracy of the approximation. For the rest of the paper, we consider the angle Ω in the equations since it has a better physical interpretation compared to the anglesω or ω. As noted before, the angle Ω represents the angle of the direction the user is facing and the X-axis of Earth coordinates. Therefore, (14) can be approximated as cos θ. For simplicity of notation, this equation can be rewritten as: where The coefficients a and b play a prominent role in the analysis of cos ψ so it is worth investigating them in detail to help readers intuitively understand them. Throughout this paper, the AP is always located above the UE, i.e., b > 0. The coefficient a can be rewritten as: Fixing the AP location, the peak of the coefficient a only depends on the position of the UE. Meanwhile, the sign of coefficient a is negative if the UE does not face the AP (opposite direction to the AP). For example, if the AP is located in the position (0, 0, 2), and the UE is located in (−1, 0, 0), then a will be negative if −π

2
< Ω < π 2 . In addition, the coefficients a and b should satisfy the following inequalities: Based on the sign of a, the term cos ψ may be a monotonically decreasing function of θ (Case 1) or it may be a concave downward function with one peak (Case 2). These two cases are explained next and the PDF of cos ψ is derived for each case. The detailed derivations are provided in Appendix-A. 1) For a < 0 (Case 1): If a < 0, b > 0 and 0 < θ < π 2 , it can be seen from (17) that cos ψ is a monotonically decreasing function of θ. Using the fundamental theorem to calculate the PDF of a function of an RV given in [35], we get the PDF of cos ψ as: where τ denotes the realization of the RV cos ψ. Let's define: as the supremum of the support of f cos ψ . This metric is presented to emphasize that ss f is not always 1. For a < 0, ss f = b which is strictly less than 1.
2) For a ≥ 0 (Case 2): If a ≥ 0, b > 0, and 0 < θ < π 2 , (17) is a concave downward function. Then, the maximum value of (17) is at the point θ * given as follows: Therefore, the PDF of cos ψ can be expressed as in (22a) for g(0) ≤ g( π 2 ) or in (22b) for g( π 2 ) < g(0), which are given at the top of the next page. Note that for a ≥ 0, ss f = √ a 2 + b 2 ≤ 1. It is worth noting that the PDF expressions given for the two cases do not particularly depend on whether θ is Laplace or Gaussian distributed so one can apply either of the distributions into (19) and (22) and calculate the distribution of cos ψ.

B. Approximate PDF of cos ψ
In order to make the performance analysis of OWC systems with random orientation tractable, one would be interested in approximating the closed-form PDF equations (19) and (22) with simpler expressions. In this section, we show that the truncated Laplace distribution can be used to approximate the PDF of cosine of the incidence angle, f cos ψ for the sitting scenario. Assuming a walking scenario, an approximation of the exact PDF as a truncated Gaussian distribution can be made similarly but was not included here.
Note that the non-zero mean random variable θ can be expressed as θ = µ L + θ , where θ follows a zero mean Laplace distribution with the same variance as θ. Substituting θ = µ L + θ in (17), and expanding the cosine and sine functions, we have: cos ψ = a sin θ cos µ θ + a sin µ θ cos θ + b cos θ cos µ θ −b sin θ sin µ θ .
Noting thatĝ(θ ) is a linear function of θ andτ min ≤ cos ψ ≤τ max , the PDF of cosine of the incidence angle can be expressed based on the truncated Laplace distribution given as follows: Fig. 4. Behavior of the sample PDF of experimental data and the closed-form PDFs with different locations of UE and fixed (x a , y a , z a ) = (0, 0, 2), z u = 0, Ω = π, and for the sitting activity.
whereτ min = −1 and taking into account the supremum of the support of f cos ψ , then,τ max = b for a < 0 andτ max = √ a 2 + b 2 for a ≥ 0. In (25), θ follows a Laplace distribution with the following parameters: Substitutingτ min andτ max in (25), the truncated Laplace distribution of cos ψ conditioned on Ω is given as: where and the support range off cos ψ is −1 ≤τ ≤τ max . Furthermore, the corresponding CDF can be obtained by simplifying thẽ F cos ψ (τ) = τ −1f cos ψ (x)dx, as follows: In order to evaluate the accuracy of our proposed truncated Laplace model based on the first-order Taylor series, we have calculated the average KSD with respect to the experimental data over different positions in a room of size 10 × 10 m 2 . The average KSD for sitting activities of the proposed model given in (28) and the exact PDF given in (19) and (22) are 0.055 and 0.026, respectively. These values are relatively small so that our exact derived PDF and the proposed approximation based on the first order Taylor series show a close similarity with the sample PDF from the experimental measurements. Note that, as indicated before, an approximation of the exact PDF for the walking scenario as truncated Gaussian distribution can be made similarly using the linearity ofĝ(θ ) in (24). Fig. 4 presents a comparison between the PDF of cos ψ given in (19) and (22) with an approximate PDF expressed in 28 for different positions of the UE with θ following the Laplace distribution. The UE is assumed to be stationary with Ω = π and located in the XY-plane, i.e., z u = 0. The AP location is (x a , y a , z a ) = (0, 0, 2). As can be seen, the first order Taylor approximate PDF provides a well-matched approximation with the experimental measurements. Taking a close look at the results, it can be inferred that for a ≤ 0, f cos ψ is still closely Laplace distributed. Although, for a > 0 and especially when it is close to 1, the shape of f cos ψ does not resemble a Laplace distribution, however, it still can be approximated with a truncated Laplace distribution.

C. Behavior of the PDF of cos ψ
In this subsection, we give a deeper analysis on the behavior of the exact PDF of cos ψ to provide more justifications to support the approximation proposed in the previous subsection. Our analytical approach is to identify the critical points of f cos ψ and analyze them. For both Laplace and Gaussian distributions, we expect a peak's location (referred as τ * ) to strictly fall between the boundaries of the support of f cos ψ . To further justify the approximation of the exact PDF as either Laplace or Gaussian distributions, it will be shown that the rate of change of f cos ψ is exponentially increasing and decreasing respectively for τ < τ * and τ > τ * .
We classify the analysis of the behavior of PDF of cos ψ into two discussions, i.e., for case 1 (a < 0) and for case 2 (a ≥ 0). 1) For a < 0 (Case 1): For case 1, we need to focus on (19) whose support is the interval (a, b).
A proof of Proposition 1 is provided in Appendix-B. To have a physical meaning of the conditions in Proposition 1 in terms of room dimension, they can be translated into following. In an indoor room with dimensions of L × L m 2 , 1 h = z a − z u , and with the AP located in the center of the room, the following condition needs to be met: To get the above condition, we use (18) and the definition of b. For L = 20 m and h = 2 m, b θ should be less than 0.1414 rad which is always true based on our experimental data. In addition, as the PDF f cos ψ increases exponentially in the left tail, i.e., a < τ < τ * , and decreases exponentially in the right tail, i.e., τ * < τ < b, the truncated Laplace distribution can be used to approximate (19).
Proposition 2 For a ≥ 0, f cos ψ has following characteristics: 1) τ * does not always exist in the support min{a, b}, √ a 2 + b 2 ; 2) f cos ψ is well defined for min{a, b} < τ < √ a 2 + b 2 ; 3) when τ * exists in the support min{a, b}, √ a 2 + b 2 , f cos ψ (τ) increases in an exponential rate for min{a, b} < τ < τ * , decreases in an exponential rate for τ * < τ < a 2 +b 2 1+b 2 θ , and increases in an exponential function for 4) when τ * exists in the support min{a, b}, √ a 2 + b 2 , f cos ψ (τ) increases in an exponential function for min{a, b} < τ < √ a 2 + b 2 ; and 5) f cos ψ is continuous at τ = τ * when τ * exists in the support min{a, b}, A proof of Proposition 2 is provided in Appendix-C. Based on Proposition 2, we have an exponentially increasing function of f cos ψ in the lower tail, but it is not always an exponentially decreasing function of f cos ψ in the higher tail as in case 1. It can be inferred from (27) and (53) that when (x u , y u ) ∈ C w and x δ = y δ = 0, then,b θ = 0. 2 To give readers a context, the set C w is a line (a hyperline in R 2 ) at which the orientation of the UE faces the AP. In other words, when (x u , y u ) ∈ C w the channel attenuation is small. Under this condition, the PDF of the approximation method turns into a probability mass function, i.e., it becomes a discrete RV. In this case, we observe a sudden jump in the KSD value [21]. In fact, this anomaly only occurs when (x u , y u ) ∈ C w . However, the area of the line C w in R 2 is zero and so this is a negligible part of the typical indoor room area. Furthermore, in our recent study [13], we have shown that if (x u , y u ) ∈ C w , the BER is less sensitive to the random orientation.

V. The Statistics of Channel Gain
In this section, we discuss the statistics of the channel gain. First, the PDF of the channel gain given the location and 1 These dimensions can also be interpreted as the domain of the locations of the UE. 2 Refer to the discussion in Appendix-C for the definition of C w . direction of the user is obtained in closed form. Then, the effect of random orientation change is described as a random process over time.

A. PDF of Channel Gain
After justifying the proposed approximation for the PDF of cos ψ, we can use the simple equation given in (28) instead of the more complicated equations given in (19) and (22) to calculate the LOS channel gain and further related derivations. Note that the PDF of the channel gain derived in this section is the conditional PDF given the location and direction of UE. Therefore, having the statistics of the user location, the joint PDF of the channel gain with respect to both UE orientation and location can be readily obtained. In the next section, we use this approach to develop an orientation-based random waypoint model to analyze mobile wireless systems.
The DC gain of the LOS optical channel can be expressed . Thus, for any given position (x u , y u , z u ), the PDF of H can be expressed as follows: where h n = H 0 d m+2 is the normalizing factor. The support range of f H ( ) is h min ≤ ≤ h max , where h min and h max are given as: The proof of (31) is provided in Appendix-D. After some manipulations, (31) can be rewritten for the sitting scenario (i.e., Laplace distribution assumption for θ) as: (34) with the support range of h min ≤ ≤ h max ; Note that in (34), H denotes a Laplace distribution with the following parameters: Moreover, for the special case of b H ≈ 0, we have f H ( ) δ( − h max ) and for the case that Pr(ψ < Ψ c ) ≈ 0, the PDF would be f H ( ) δ( ). Fig. 5 represents the PDF of the LOS channel gain obtained from analytical results given in (34) and the measurement-based simulations. The results are provided for different positions in the room with Ω = π 4 . The simulation parameters are provided in Table II. The results show the accuracy of the derived analytical PDF. The magnitude of the Dirac delta term is almost zero in Fig. 5-(a) and Fig. 5-(c) while it is 0.0336 and 0.006 in Fig. 5-(b) and Fig. 5-(d), respectively. Fig. 5-(a) and Fig. 5-(b) illustrate two positions that correspond to a = 0 and a < 0, respectively while Fig. 5-(c) and Fig. 5-(d) present the positions associated to a > 0. As can be seen, for a ≤ 0, the analytical derivation of LOS channel gain and the measurement-based simulations match very well. For a > 0, we still observe a good accommodation between analytical and simulation results and in fact the difference happens at = h max as shown in Fig. 5-(c) and Fig. 5-(d). We also observe that the distribution of the channel gain significantly changes from the Laplacian shape as a −→ 1 or equivalently (x u , y u ) approaches to the C w region as shown in Fig. 5-(d).
SNR Statistics: The received electrical SNR in an optical wireless channel can be obtained as: where R PD is the PD responsivity; P opt is the transmitted optical power; N 0 is the single sided power spectral density of noise; B is the modulation bandwidth. Let's define S 0 R 2 PD P 2 opt N 0 B , then we have S = S 0 H 2 . Using the fundamental theorem of determining the distribution of a random variable [35] and noting that H ≥ 0, we have: with the support range of s ∈ (s min , s max ), where s min = S 0 h 2 min and s max = S 0 h 2 max , with h min and h max given in (32) and (33), respectively.

B. Random Process Model for Device Orientation
If a fixed orientation is assumed for the UE with an angle ψ 0 ∈ [0, Ψ c ], the LOS channel gain remains constant and equals to H = H 0 cos ψ 0 d m+2 . However, due to the random orientation of the UE, the incidence angle indeed fluctuates around the angle ψ 0 as observed in the experimental measurements. Thus, the LOS channel gain also varies and it can be observed as a stochastic random process (RP). Note that in a realistic scenario whether the user is in a walking or sitting position, it is fair to assume that the angle Ω (corresponding to the direction of the movement) is fairly stable and does not fluctuate with the same rate as the polar angle. We therefore first focus on the conditional probability of the channel gain given a particular direction of the movement. We then in the next section introduce a random waypoint model to consider the random effect of Ω (random change of direction) on the channel gain and thus the communication system performance. Fig. 6 illustrates one ensemble of the stochastic processes β, γ and θ. As can be seen, the pitch angle, β, varies around the mean value of E[β] = −35.81 • and the roll angle, γ, fluctuates around its mean value which is zero. This means that the user tends to hold its cellphone with the pitch angle of −35.81 • and roll angle of zero. However, due to the random nature of users behavior, fluctuations are observed in pitch, roll and polar angles. The variations in pitch and roll angles also affect the LOS channel gain. Fig. 6 shows three ensembles of the normalized LOS channel gain, H H 0 , for UE's locations of (x u , y u ) = (−2, −2), (x u , y u ) = (3, 3) and (x u , y u ) = (0, 0) with Ω = π 4 . To characterize the temporal behavior of the device random orientation as a random process, the coherence time of the channel needs to be known. Autocorrelation is one common way to calculate the coherence time of a stochastic random process. Denoting the autocorrelation of θ as R θ (t), the coherence time of the polar angle stochastic process denoted by T c,θ represents the time it takes for the random process to become uncorrelated with its initial value and can be defined as R θ (T c,θ ) = 0.05R θ (0) = 0.05 [36]. As expressed in section III, our experimental results show that the coherence time of the angle θ is in the order of a few hundreds of milliseconds [33]. However, the coherence time of the LOS channel gain highly depends on the UE's location and direction (determined by the angle Ω), which is in the order of a few tens of milliseconds. Accordingly, the effect of random orientation may be modeled as a slow large-scale fading effect as the coherence time of the channel variation due to random orientation, T c,H , is very large compared to the typical symbol period in OWC. This implies that the channel gain can be presumed as roughly constant over a large number of transmitted symbols [36].

VI. Orientation-based Random Waypoint
A commonly used mobility model in simulation based studies of wireless networks is the random waypoint (RWP) model. According to the RWP mobility model, i) users choose their destinations randomly in the room area, ii) they move with a constant speed on a straight line between a source and a destination. Note that destinations are uniformly randomly distributed in the room area [37]. RWP mobility model is described as a discrete-time stochastic process. Here, we consider the RWP movement in two-dimensional space. The angle between the direction of movement and the positive direction of the X-axis is defined as the angle of direction. This angle is the same as Ω given in (5) and (6). The RWP can be mathematically denoted as an infinite sequence of triples, (P n−1 , P n , v n ) n ∈ N , where n stands for the nth movement period. The UE moves from the random waypoint P n−1 = (x n−1 , y n−1 ) to the destination point P n = (x n , y n ) with speed of v n chosen from a random distribution of f v . The Euclidean distance between two consecutive waypoints, P n−1 and P n is defined as the transition length and is given by D n = P n − P n−1 . The sequence of transition lengths {D 1 , D 2 , ...} are independent identically distributed (i.i.d.) RVs with the PDF described in [37], and the expected transition length of E[D] = 0.5214L in a room of size L × L m 2 . The elapsed time between two successive movements for the nth period can be obtained as T e,n = D n /v n . The expected value of the elapsed time is given as E[T e ] = E[D]E 1 v . In the context of LiFi as well as mmWave cellular networks, the effect of UE's orientation on the performance of the system is significant. In fact, a significant change of UE's orientation, whether individually or combined with the mobility of the user, can lead to a handover that would not normally happen for UEs with a constant orientation. So in order to provide a framework to analyze the performance of mobile wireless networks more realistically, we need to combine the conventional RWP with the random orientation model. Therefore, the orientation-based random waypoint (ORWP) can be modeled as an infinite Algorithm 1 Orientation-based random waypoint (ORWP) 1: Initialization: n ←− 1; k ←− 1; denote P n = (x n , y n ) as the nth location of UE and P 0 = (x 0 , y 0 ) as the initial UE's position; N r as the number of runs; v as the speed of UE; T c,θ as the coherence time of the polar angle; E[θ] and σ 2 θ as the mean and variance of Gaussian RP; 2: for k = 1 : N r do 3: Choose a random position P k = (x k , y k ) 4: Compute D k = P k − P k−1

5:
Compute Ω = tan −1 y k −y k−1 Compute P n = (x n , y n ) with x n = x n−1 + vT c,θ cos Ω and y n = y n−1 + vT c,θ sin Ω 9: Generate θ n based on the AR(1) model 10: Return (P n−1 , P n , v, θ n ) as ORWP specifications 11: n ←− n + 1 12: t move ←− t move + T c,θ 13: end while 14: if t move Generate θ n based on the AR(1) model 16: P n ←− P k 17: Return (P n−1 , P n , v, θ n ) as ORWP specifications 18: n ←− n + 1 19: end if 20: k ←− k + 1 21: end for sequence of quadruples, (P n−1 , P n , v n , θ n (t)) n ∈ N , where θ n (t) is a random process describing the UE's polar angle during the movement from waypoint P n−1 to waypoint P n . More discussions about how to generate these sequences are presented next and the ORWP is summarized in the Algorithm 1.

A. Correlated Gaussian Random Process
As shown in section III, the polar angle for walking activities follows a Gaussian distribution. The experimental measurements also illustrate that the adjacent samples of the RP, θ, are correlated. Therefore, in order to incorporate the orientation with the RWP mobility model, it is required to generate a correlated Gaussian RP that statistically matches the experimental measurements. Possible ways of generating a correlated Gaussian RP can be found in [38], [39] and references therein. A simple method to generate a correlated Gaussian RP is to pass a white noise process through a linear time-invariant (LTI) filter, e.g., using a linear autoregressive (AR) model. Let w[n] denote the white noise process, then, after passing it through the LTI filter, the nth time sample of the correlated Gaussian RP, θ[n], is given as: where c 0 determines the biased level and c i for i = 1, . . . , p are constant factors of the AR order of p, AR(p). Note that AR(p) contains p + 2 unknown parameters including: c 0 , c 1 , . . . , c p , σ 2 w , where σ 2 w is the variance of white noise RP, w. In this study, we focus on matching the generated random process to the moments and the coherence time of the polar angle measured experimentally rather than its exact ACF. Therefore, it is sufficient to assume p = 1 and consider first-order AR model to generate the correlated Gaussian RP. Thus, the nth sample of the AR(1) model can be expressed as: where c 1 should meet the condition |c 1 | < 1 to guarantee the RP, θ, is wide-sense stationary. Noting that for AR(1), the mean, variance and ACF are given as [40]: Using the above equations and noting that R θ ( θ = T c,θ T s ) = 0.05 where T s is the sample time, we have: Once the parameters of the AR(1) model are determined, the nth time sample of the correlated Gaussian RP, θ, can be specified according to (40). Using the method described above, the ORWP is presented in the Algorithm 1.

B. Use Case: Handover Rate
Here, we investigate the effect of the UE's orientation on the handover rate as a case study. In fact, one of the key metrics in cellular network design is the handover rate which is defined as [9]: . For a room of length L the expected number of handover can be obtained as [9]: where the UE is assumed to move from the initial RWP (x 0 , y 0 ) to the destination point (x, y). Fig. 7 shows the Monte-Carlo simulations of handover rate as a function of room length for a UE moving at a speed of v = 1 m/s, v = 1.4 m/s and v = 2 m/s. These three speed values are chosen around the average human walking speed [41]. In this simulation setup, four APs are assumed that they divide the network area into four separate quadrants (attocells) with one AP at the center of each attocell. The UE is assumed to be initially connected to the corresponding AP denoted as AP j for j ∈ {1, 2, 3, 4}. We compared a vertically upward UE (shown by the solid lines in Fig. 7) with a UE that follows the orientation-based RWP mobility model in which the random orientation of the UE is generated based on the correlated Gaussian RP (shown by the dotted line).  The mean and variance of the generated samples are chosen from Table I for walking activities, which are E[θ] = 29.67 and σ 2 θ = 7.78. Furthermore, the coherence time of T c,θ = 130 ms as reported in [33] for walking activities is chosen. As observed from Fig. 7, the handover rate is decreasing overall with an increase in the network dimensions. On the other hand, as the UE's speed increases, the handover rate increases. This is more remarkable when the network dimensions are smaller. More importantly, we observe that the effect of random orientation significantly increases the value of handover rate in all conditions while smaller network sizes are more influenced by this effect. That is, as the room length is smaller, the gap between the random orientation and vertically upward scenarios is greater.

VII. Conclusions
A new model for device orientation based on the experimental measurements is proposed. The experimental measurements are taken from forty participants creating 222 datasets for orientation. It is shown that the PDF of the polar angle follows a Laplace distribution for static users while it is better fitted to a Gaussian distribution for mobile users. The exact PDF of the cosine of the incidence angle is derived with analytical discussions. We also proposed an approximation PDF using the truncated Laplace distribution and the accuracy of this approximated method was confirmed by using the KSD test. The LOS channel gain statistics are calculated and it is described as a random process for which the coherence time was discussed. The influence of the random orientation on received SNR of OWCs is evaluated. By means of AR(1) model, an orientation-based random waypoint mobility model is proposed by considering the random orientation of the UE during the user's movement. The performance of the mobility model is assessed on the handover rate and it is shown that it is important to take the random orientation into account.
2) For a ≥ 0 (Case 2): For case 2, we need the intervals of the domain of cos ψ such that the transformation (44) is one-to-one. Using (17) and (21), we have two intervals which are g(0) < cos ψ < g(θ * ) and g( π 2 ) < cos ψ ≤ g(θ * ). The inverse functions are denoted by: Note that both the inverse functions have the same the Jacobian of the transformation which is |J| = 1/ a 2 + b 2 − cos 2 φ.

B. Proof of Proposition 1
Proof: To prove all characteristics in Proposition 1, let's first define some critical points that make f cos ψ undefined. The possible points are ± √ a 2 + b 2 . For a < 0 and b > 0, it is straightforward to see the following inequality: Therefore, these critical points are outside the support of f cos ψ . Considering the absolute value term in (19), the other critical point called τ * is given by: This point leads to a peak for f cos ψ similar to the peak of the Laplace or Gaussian distributions. It is straightforward to see that a < τ * < b, if 0 < µ θ < π 2 . Using (47), we have: Hence, this proves that f cos Ψ is well defined for a < τ < b and continuous at τ = τ * . To prove the first and the third characteristics, let's denote the following: where Note that we base this on the Laplace distribution to derive f 1 and f 2 for simplicity. For the Gaussian case, they have similar forms, but the argument inside the exponential function is squared and both f 1 and f 2 are multiplied by a function of: Note that this does not change the fact that the function f cos ψ is still a monotonically increasing or decreasing function in an exponential manner.
Since f 1 (τ) and f 2 (τ) are both positive and monotonic functions in their domains of interest, we need to focus on the coefficients C 1 (τ) and C 2 (τ) as they determine the monotonicity of f cos ψ .
Let's focus on the scenario where x a x u or y a y u . Since the denominators are always positive, we only need to observe the numerators, so it is straightforward to have: Satisfying the above inequalities guarantees that f cos ψ (τ) is a monotonically increasing function for a < τ < τ * and a monotonically decreasing function for τ * < τ < b. The consequence of satisfying (50) is that τ * is a global maximum of f cos ψ in the interval (a, b). This completes the proof of Proposition 1.

C. Proof of Proposition 2
Proof: As in Appendix--B, we have following inequality: It means that for the lower tail, the support does not include the term − √ a 2 + b 2 , which makes f cos Ψ undefined. To see the higher tail, it is not straightforward and we need to introduce a region C w so that for given (x a , y a , z a ), z u , and Ω, the CDF of cos ψ is almost 1 in the vicinity of maximum value of cos ψ. Thus, C w can be formalized as: = (x u , y u ) x u = x u , y u = y u and tan Ω = − x δ y δ , for small > 0. Here, x δ and y δ are two auxiliary variables to move the UE position along the line C w in the room. Note that C w for x δ = y δ = 0 can be expressed in terms of the coefficients a and b as: Considering the absolute value term in (22), the critical point τ * is given as: It is obvious to see that if 0 < µ θ < π 2 , min{a, b} < τ * ≤ √ a 2 + b 2 , and it is equal when (x u , y u ) ∈ C w . Therefore, unlike the case 1 (a < 0), τ * may not exist in the interval min{a, b}, √ a 2 + b 2 . It is when τ * = √ a 2 + b 2 which is outside of the support of f cos ψ , and the case τ * = √ a 2 + b 2 only happens when (x u , y u ) ∈ C w . Hence, this proves the first characteristic in Proposition 2.
If (x u , y u ) C w , we have min{a, b} < τ * < √ a 2 + b 2 and a continuity at τ * , i.e.,: This proves the fifth characteristic in Proposition 2. In addition, we have the following identity: which says that for case 2, the end tail of f cos ψ always tends to infinity. Therefore, discussion about the rate of change in case 2 is as not straightforward as that in case 1.
The function f cos ψ in (22) can be rewritten as a combination of the following two expressions: Since the derivatives of f a (τ) and f b (τ) give the same coefficients as those in the Appendix--B, we are interested only in the coefficients as expressed in (49). However, in case 2, C 1 (τ) is used for the support min{a, b} < τ < τ * and C 2 (τ) is used for the support τ * < τ < √ a 2 + b 2 . For min{a, b} < τ < τ * , C 1 (τ) is always positive. It means that we always have a monotonic increasing function of f cos ψ with the support of min{a, b} < τ < τ * . Note that for the case (x u , y u ) ∈ C w , we have τ * = √ a 2 + b 2 by definition of C w in (53). In this case, we have only a monotonic increasing function that goes to infinity at √ a 2 + b 2 . For τ * < τ < √ a 2 + b 2 , let's first define which is the point where the coefficient C 2 (τ) starts changing from negative to positive values. This expression is obtained as in (50) by only observing the sign of the numerator of C 2 (τ). Then, we have the following result for C 2 (τ): τ * < τ d < √ a 2 + b 2 =⇒ C 2 (τ) < 0 for τ * < τ < τ d , and If τ d < τ * , C 2 (τ) is always positive meaning that f cos ψ is a monotonically increasing function. This proves the second to the fourth characteristics in Proposition 2 and completes the proof of Proposition 2.

D. Proof of (31)
Proof: The CDF of the LOS channel gain, H = H 0 cos ψ/d m+2 rect ψ Ψ c is given as: where U( ) is the unit step function; the support range of F H ( ) is h min ≤ ≤ h max . Then, the corresponding PDF can be obtained as follows: c N h n f cos ψ h n + F cos ψ (cos Ψ c ) δ( ).
(55) where in the last equation, we define h n H 0 /d m+2 . The Dirac delta function comes out due to the discontinuity of CDF given in (54) at = 0. To ensure that the h max h min f H ( )d = 1, the normalizing factor should be c N = 1. This completes the proof of (31).