LIDAR-Assisted Channel Modelling for LiFi in Realistic Indoor Scenarios

We present a new, accurate, low complexity channel modelling methodology for LiFi in realistic indoor scenarios. A LIDAR scanner is used to capture the 3D environment in which the LiFi system is to be deployed. Next, the generated 3D point cloud data is pre-processed to determine the reflectance parameters of the walls and objects in the room. This is easier and more realistic than the manual definition of the environment, which is the current state of the art. As an additional innovation, the complexity of the channel modeling is reduced by i) modeling the line-of-sight and initial reflections precisely in the frequency-domain and ii) using a well-established analytical model based on the integrating sphere for all higher-order diffuse reflections. All steps together yield a substantially simplified channel modelling approach and model the links between multiple optical frontends and multiple mobile devices realistically. As a validation of our new approach, we compare measurements and simulations in two indoor scenarios: an empty room and a conference room with furniture. Simulations and measurements show excellent agreement with a mean square error below 3 percent. Moreover, we evaluate the performance of a distributed multiuser multiple-input multiple-output (MIMO) link and found excellent agreement between the model and measurements. Finally, we discuss the fundamental trade-off between complexity and model error, which depends on the scenario.

geometrical and optical properties of the walls and all other objects in the room.
Several channel modelling methods were already introduced for LiFi, which can be mainly classified into deterministic [7]- [12] and stochastic approaches [13]- [15]. Deterministic models depend on architectural plans of the environment and on the position and orientation of transmitter (Tx) and receiver (Rx). Stochastic models are generic, non-site specific and parameterized according to the measurements taken in specific Tx, Rx, and scattered environments. The stochastic approach is much simpler and useful because of its reduced computational complexity. Deterministic models, on the other hand, are more realistic and include real-world effects such as blockages, partial shading of a link and reflections from nearby objects. In scenarios where the LOS is usually free, the generic approach may be enough for modelling LiFi channels. However, in scenarios with a non-negligible probability of LOS blockage and reflections from the non-line-of-sight (NLOS) path, a deterministic approach is useful. The major challenge is to model sufficient number of channel responses realistically so that the impact of the environment is fully included. The existing deterministic channel modeling techniques for LiFi are not efficient enough for these requirements. Another drawback is that only few LiFi channel models have been validated through measurements in real scenarios [7], [11], [12].
In the recent literature, there are first approaches to include objects inside the room [12], [15], where each object (besides all reflecting surfaces) is currently defined manually. This is very time consuming and lacks the precision required to determine the availability of the LOS and to correctly include blockage in a complex environment. Moreover, the recent channel modelling approach uses commercial optical design software, e.g., Zemax®, which is based on ray tracing [12]. This framework makes it possible to include variable lightemitting-diode (LED) models and to consider both specular as well as mixed specular-diffuse reflections. Moreover, it can include the wavelength-dependent reflectivity of materials. However, ray tracing operates in the time domain, and follows each path until it is either absorbed or received. Modelling multiple-input multiple-output (MIMO) channels and mobility becomes very complex when using this recent approach. The major challenge is to reduce complexity, and it can be decomposed into two parts. First, it is possible to reduce the effort for the 3D scenario generation. Second, one can reduce the computational complexity for the channel modelling algorithm. The present paper introduces new techniques addressing both steps to reduce complexity and compares the new techniques with measurements for validation.
LIDAR scanners are widely used to scan the 3D environment for planning the deployment of radio frequency (RF) based mobile communication systems [16]- [18]. This approach is also applicable for LiFi. In this paper, for the first time, we use a LIDAR scanner to generate the three dimensional (3D) environment for the LiFi scenario. Additionally, one needs to know the reflectance of the surfaces in the room to predict the LiFi channels. Therefore, we have developed a technique to estimate the reflectance parameters of all surfaces directly from the raw LIDAR data and to convert the point cloud data into a file, which can be directly used as input for the channel modelling algorithm. To further reduce complexity, we consider that the required precision of the 3D indoor environment is related to the intended bandwidth of the LiFi system. If the intended bandwidth is low, less time resolution is needed and hence, the modelling of the 3D environment can use lower spatial resolution, while more precise 3D modelling is needed if more bandwidth is targeted. Overall, this approach allows us to replace the manual definition of objects and surfaces in 3D. Preliminary results were reported in [19].
Previous studies show that only the LOS and the first one or two reflections can be identified in the impulse response, while all later reflections merge into a long exponential decay due to the diffuse nature of the environment [11]. Therefore, it is important to model LOS as well as initial NLOS reflections precisely. In our channel modelling algorithm, we use the frequency-domain (FD) method to model the LOS and the initial NLOS reflections [10], which is inherently faster than comparable ray tracing tools operating in the time domain. For higher-order diffuse reflections, it creates a nearly homogeneous illumination inside the room and the contribution to the channel can be explained as a first-order low-pass response as given in [11]. Considering the geometry and the reflection of all objects and walls in the room, we can calculate the contribution of the higher-order diffuse reflections using an analytical formula proposed in [11]. This formula reduces the overall complexity considerably, without losing accuracy. Therefore, we only model the LOS and the initial NLOS reflections using FD method [9] and higher-order diffuse reflections using the integrating sphere model [11]. The new algorithm and initial results were reported in [20].
To validate our approach, we perform channel modelling and measurements in two realistic scenarios, an empty room and a conference room. The empty room scenario has few reflecting surfaces, which is useful for initial testing. However, in the conference room scenario, the room is filled with many objects such as chairs and tables. At first, we perform LIDAR measurements by placing the LIDAR scanner at multiple locations in these rooms. Each scanned dataset is processed and optimized to generate the channel responses for all Tx-Rx links. Finally, channel measurements are performed by keeping Tx and Rx in the same locations as in the modelling. The results demonstrate that the simulated channels are in very excellent agreement with the measured channels. Moreover, we evaluated the performance by estimating the throughput when operating the LiFi system as a distributed multiuser MIMO system and found a good agreement between the model and measurements. Finally, we discuss the trade-off between complexity and the required precision in practical use cases.
Altogether, this paper shows that the complexity of channel modelling for dense LiFi networks can be reduced and computationally efficient, precise and highly realistic results can be achieved. The main innovations in this paper are i) the LIDAR based 3D modelling, ii) an efficient channel modelling approach for LiFi and iii) the validation by measurements in different scenarios.
The remainder of the paper is organized as follows. Section II describes the 3D modelling approach using a LIDAR scanner and discusses the pre-processing for reflectance parameter estimation from the LIDAR data. Section III introduces an efficient LiFi channel modelling methodology. Section IV describes our LIDAR measurement campaigns for capturing 3D environments and the LiFi channel measurements. Modeling and measurement results are compared in Section V. Finally, conclusions are given in section VI.

II. 3D MODELLING USING LIDAR
Nowadays, LIDAR scanners are commonly used to capture 3D models of indoor environments for applications like redesign, visualization, monitoring and simulations [21]- [23]. To accurately capture the indoor environment, the LIDAR scanner is placed at multiple locations [24]. In this paper, we used this common approach to capture the 3D environment with a Leica RTC 360 scanner. The captured point cloud data from each scanning location are then pre-processed. In this section, we describe the LIDAR data processing and the reflectance estimation for the surfaces.

A. LIDAR DATA PROCESSING
The original output data of the LIDAR scanner consist of (x, y, z) spatial coordinates of the reflected points and the intensity of the light reflected from each point. Our goal is to use these point cloud data as input into the channel modelling algorithm.
Before using this data as an input to the channel modelling algorithm, the point cloud data needs to be optimized. There are three major steps after the scanning.
1. Most of the LIDAR scanners have difficulties to detect transparent and specular reflecting surfaces (e.g. mirrors and glass surfaces) [25]. In addition, due to dynamic changes in the environment during the scanning, some of the points can be scattered. Hence, such noisy points should be removed from the point cloud data.
2. LIDAR scanners on the market have a resolution down to below 2 millimeter and can record millions of points within a few seconds. Due to the reduced bandwidth of LiFi systems, we only need the point resolution in the centimeter range. Hence, we can significantly reduce the resolution and thus the size of the point cloud data set.
3. In our channel-modelling algorithm, we need the surface normal of points as additional information besides the Cartesian coordinates. The normal is required to calculate the angle of incidence of light at each point.
All three steps are well covered research topics in the literature. There are many methods to remove the noisy points [26]- [28], calculate the surface normal [29]- [30], and reduce the resolution of point cloud data [31]- [32]. However, as a proof of concept, we use the open-source software CloudCompare to post-process the point cloud data [33]. In this software, at first all noisy data are removed, then the point cloud data set is down-sampled, and finally, surface normals are calculated. After processing, each data set is exported into the polygon file format (ply).

B. REFLECTANCE ESTIMATION
For the channel modelling, reflectance of the surfaces inside the room is needed. Most LIDAR scanners are designed for range measurements rather than intensity measurements [34]. The intensity information of the point cloud data is based on a measurement of the electrical signal at the receiver side which is obtained by converting and amplifying the backscattered optical power of the emitted signal. This intensity information depends on the reflectivity of the surfaces, the distance between the laser source and the target, the angle of incidence of the laser beam with respect to the surface normal and the amplifiers inside the LIDAR device [35]. Hence, intensity needs to be corrected by removing those dependencies and evaluate the reflectance for each target point [36]. Implementing and testing an appropriate method for reflectance calculation is a major contribution in this paper.
Radiometric correction is the term for the process of converting the raw intensity data to a value proportional to the surface reflectance at each scanned point [35]. Recently, different correction methods have been introduced to evaluate the reflectance from the raw intensity data [35], [37]. Since the LIDAR scanners internal variables are unknown, e.g. laser and receiver temperatures, most of the existing methods only correct the effects due to geometrical parameters such as the distance and the angle of incidence. Therefore, it is only possible to obtain a pseudo-reflectance, which is an increasing function depending on the reflectance of the surface points [37].
Radiometric correction methods are mainly divided into two approaches: model-driven and data-driven. The modeldriven approach is based on a theoretical model, the so-called LIDAR range equation [37]. The data-driven approach is based on quantified observations of the intensity measurement, which is more adaptive and suitable for the correction. Here, we follow the data-driven approach to calculate the pseudo-reflectance of the scanned points.
We follow the approach in [38] to evaluate the pseudoreflectance from the intensity data. In this method, the measured intensity is expressed as three separable functions, which cover the effects of reflectance, angle of incidence and distance. Therefore, the intensity recorded by the LIDAR can be expressed as: where is the reflectance varying from 0 to 1, is the incidence angle of the laser beam with respect to the target point, and is the distance from the laser source to the target point. The functions 1 , 2 , and 3 is given in more detail in equation (10) in [38]. From the measured intensity value (1), the effects of attenuation due to the distance and the angle of incidence should be corrected. Therefore, functions 2 and 3 need to be determined. To do so, additional intensity measurements were carried out at different distances and angles with commercially available reflectance targets. We follow the measurement procedure similar to [37]. The objective is to determine the parameters and coefficients, which control the functions 2 and 3 . Therefore, we consider four Lambertian diffuse surfaces with known reflectance values 99%, 80%, 50%, and 20%. The 99% reflectance target is of 5 x 5 cm 2 size while the 80%, 50%, 20% targets are of 10 x 10 cm 2 size. The reflectance value of these targets is calibrated by the manufacturers in the range between 400 nm -2000 nm. The measurement has been done by keeping the targets at distances between 1 m to 16 m and varying the angle between 0 degree to 70 degree. From the measured data, we generate a best-fitting polynomial and determine all the required parameters for 2 and 3 . From these determined polynomial functions, we can calculate the corrected intensity of the targets as: where is the measured intensity and = (5, 0) is a normalization factor respective to the intensity calculated for the target kept at 5m at an angle of 0 degree. Using equation (2), the corrected intensities acquired at the other incidence angles and distances have been corrected with respect to 0 degree and 5 m, which has been chosen as the reference scanning geometry. The corrected intensity can be expressed as: Using (3), we correct the intensity of every scanned point and nullify the effects due to distance and angle similar as mentioned in [37]. Note that this correction method is independent of the reflectance of the targets. Fig. 1 shows the experimental setup for LIDAR intensity correction measurements using the Lambertian reflectance targets. The Lambertian targets used in our experiment are shown in Fig. 1(a) and Fig. 1(c). At first, the targets are kept at different distances starting from 1 m (see Fig. 1 Fig. 1(d)) to 16 m with an increment of 1 m, in order to take a LIDAR scan. In the same environment, we performed further measurements by varying the angle, while the Lambertian targets are kept at 5m distance and tilted from 0 degree to 70 degree with an increment of 5 degree.
From the measured 3D point cloud data, the intensity profile received from each target was extracted by using the software CloudCompare. Then we calculated the average intensity for each target at its given distance and angle. Fig. 2 shows the results after applying the intensity correction for different Lambertian targets. Intensity as a function of distance is shown in Fig. 2(a) and as a function of angle in Fig. 2(b). From these results, all parameters were estimated the same way as described in [37]. Polynomial regression orders were chosen in regard of the minimal mean error, i.e. polynomial order N2 = 2 for angle of incidence, and polynomial order N3 = 6 for the distance. The polynomial coefficients were normalized by leading term chosen as 1.
The normalization constant was chosen at the 50% reflectance curve at 5 m, at an incidence angle of 0•. Intensity correction was performed using equation (3) and finally scaled to the range 0-1. Final corrected intensities for different angles and distances are shown in Fig. 2(c) and Fig. 2(d), respectively. When comparing the corrected intensities to the ideal reflectance values of Lambertian targets, the mean square error (MSE) is below 5% for angles less than 50 degrees and MSE is more than 10% for angles larger than 60 degrees. With respect to distances, MSE is less than 2%. These results validate that we can evaluate the reflectance with minor error by using the described intensity correction method [37, [38].
In a similar way, using equation (3), the reflectance values were estimated from the raw intensity data of a 3D room. Example results for two considered room scenarios are reported in the Appendix A. Note that, this method can only detect the diffuse reflectance targets since specular reflections would saturate the LIDAR receiver.
Finally, all scanned data are merged to get the complete point cloud data of the room. The data set consists of Cartesian coordinates of the scanned points, the surface normal at each scanned point, and the reflectance value belonging to each point. This pre-processed LIDAR data set is used as an input into our channel modelling algorithm.

III. CHANNEL MODELLING METHODOLOGY
In this section, we describe our channel modelling approach for LiFi by using 3D point cloud data. Modelling is performed in the frequency domain rather than in the time domain. Hence, this method delivers channel transfer functions instead of impulse responses. Note that LiFi links are modelled separately for LOS and NLOS links.

A. LOS CHANNEL MODEL
The LOS channel model is based on the orientation and optical parameters of Tx and Rx as described in [7]. The generalized LOS channel model is given as where , is the transfer function coefficient between Tx and Rx which can be described as where 0 ( ) is the radiation intensity of the optical transmitter for each elevation angle , which is given with respect to the surface normal. Here, ( ) is the angulardependent effective detector area at the Rx side. Total optical power emitted from the Tx is given as . The coefficient , depends on the FOV of the photodiode and the radiation pattern of the optical transmitter [9]. The variable , is the propagation time which depends on , and being the distance between Tx and Rx and the speed of the light, respectively [10]. The visibility factor , is equal to one if there is a free LOS between Tx and Rx and zero otherwise [19].

B. INITIAL NLOS REFLECTIONS
We consider the reflecting surfaces as surface elements. From the LIDAR data, we consider that each point in the point cloud data represents the center position of the corresponding surface element. The average resolution of the point cloud data defines the size of each surface element. The FD method allows assembling all mutual LOS links between all surface elements as well as the links between the surface elements and the Rx and Tx in a matrix form and to compute NLOS reflections by consecutive matrix multiplications [10].
For a single Tx to Rx scenario, as described in section III in [10], the entire NLOS channel model can be represented as where is the reflection order, ( ) is the matrix consists of LOS transfer functions for the links from Tx to all surface elements, is a diagonal reflectivity matrix, where each diagonal element represents the reflectivity of the ℎ surface element, T ( ) is the matrix consist of LOS transfer functions for the links from each surface element in the room to the Rx. The intrinsic function ( ) of the room is a matrix where each element , ( ) in this matrix represents the LOS transfer function from the ℎ surface element to the ℎ surface element.
Our main objective is to simplify the modelling and achieve high precision with limited complexity. Therefore, we model only the first reflections precisely, so that we consider = 1, 2, 3 accordingly. This means that the matrix multiplications in equation (6) are calculated in a recursive manner as described in section III. C in [10].
When compared to the mathematical descriptions given in section III in [10], we modify the definition of the matrices ( ), ( ), and ( ). Each element in these matrices contains an additional visibility factor (similar to equation (4)) which defines the visibility in the link. The visibility factor is equal to one if the link is free and zero otherwise.
For example, an element in the matrix ( ) can be represented as where , , , and , are the visibility factor, the transfer function coefficient and the delay time between ℎ and ℎ surface elements, respectively.

C. VISIBILITY FACTOR
Visibility analysis in the point cloud data is a vital research problem in the field of computer graphics, computer vision, robotics and photogrammetry [39]- [44]. Usually, visibility is estimated from a certain viewpoint e.g. the camera position [40]. The estimation of visibility of a point cloud consists of assigning a label to each point of the scene. The label marks the point as visible if it lies on an object that is directly visible from a given viewpoint and non-visible otherwise [41]. This is a fundamental step in many computer graphics applications. There are many methods to estimate the visibility in the point cloud data such as surface reconstruction methods, convex hull methods, and likelihood methods [40]. Among these, use of the hidden point removal (HPR) operator, which is an example for the convex hull method, is known as a simple and fast solution with sufficient accuracy [42]. The advantage of this technique is to avoid creating surfaces from the point cloud data and to analyze visibility efficiently both in sparse and dense point clouds. By replacing the view point with Tx, Rx and surface point locations, we can estimate which points are visible from each other. Detailed discussion of HPR based visibility analysis is given in the Appendix B where we describe the selection of parameter R for the HPR based visibility analysis. From our study, we observe that, the HPR method estimates visibility for our LiFi channel modelling application with the lowest computational time and error.

D. HIGHER-ORDER REFLECTIONS
From our previous channel measurements, we observe that the higher-order diffuse reflections have all very similar characteristics and could be considered jointly [45]. It is intuitive that higher-order diffuse reflections create a nearly homogeneous illumination in the room and depend more on the environment than on the orientation of Tx and Rx. This motivates the use of a heuristic model for the higher-order reflections. They can be modelled jointly by using Ulbricht's integrating sphere model, which has been adapted to regular room dimensions [11].
In contrast to the microscopic approach described in section III B, here the macroscopic approach does not include any details of the room. Rather, an analytic formula for the channel transfer function is used that depends on basic parameters such as total surface area, volume and average reflectance of the room. These parameters can be obtained from the point cloud data of the room.
The generalized diffuse channel model for the considered room is given by equation (9) in [11]. Since initial reflections are calculated using the FD method, we need to calculate only the higher-order diffuse reflections, which can be calculated by removing initial diffuse reflections from equation (9) in [11]. The higher-order diffuse reflections are expressed as where is the exponential decay time which is related to the room parameters and ∆ is the delay time, compared to the LOS, after which the diffuse components arrive at the receiver side [11]. The variable is the diffuse channel gain. After excluding the initial diffuse reflections, it can be expressed as where is the area of the room surface, is the area of the receiver, 1 is the reflectivity of the region initially illuminated by the Tx and 〈 〉 is the average reflectivity of the room [11]. The variable is the reflection order until which reflections are precisely modeled by using the FD method.
We consider the delay time Δ is approximately equal to 〈 〉 which is the average time between two reflections as given in equation (13) in [11]. We observe that in most of the typical indoor scenarios Δ is less than or equal to 10 ns. This model provides an approximate result for higher order reflections in a given room.
Finally, the complete transfer function of the LiFi channel can be represented as Using Equation (10), we calculate the channel transfer function at each sampling frequency .

E. MIMO AND MOBILITY
For a mobile MIMO scenario, at first, the LOS corresponding to each Tx-Rx link is calculated using equation (4). Then, for the first NLOS reflections, we follow the same approach by using the FD method as given in [10] (see section III. G), where T ( ) and ( ) are calculated for each Tx and Rx configuration in the room. In the FD method, the NLOS channelfor a given MIMO scenario is represented by equation (47) in [10]. It is a general advantage of the FD method that, for a given scenario, the intrinsic transfer function matrix ( ) of the room is computed once and can be used for all Tx-Rx links in a MIMO. This helps to further reduce the computational time. Furthermore, in case of a mobile MIMO scenario with fixed transmitters, we only need to update T ( ) corresponding to the varying position of the receivers in the room.
Due to the integrating sphere method, all results for the higher-order diffuse reflections do not depend on the Rx orientation. Therefore, in the downlink, where the Tx is static, we can calculate the higher-order contributions only once. Otherwise, one would need to calculate the higherorder diffuse reflections for each Tx position.

F. LIMITATIONS
For simplicity, we considered the reflecting surfaces inside the room as Lambertian for estimating the reflectance values from raw LIDAR data using equation (3). The channel model is not limited regarding the surface reflection characteristics. It is possible to define the surfaces with any reflectance parameter and characteristics (like Phong reflection model [46]). Moreover, we did not consider wavelength dependency. Wavelength-dependent channel models for LiFi include three main effects: i) the ℎ ℎ ( ) = 1+ 2 2 ∆ , spectral characteristics of the light source, ii) wavelengthdependence of the surface reflectance and iii) responsivity of photo detectors [14], [47]. LIDARs use a monochromatic source, and thus it is difficult to measure wavelength dependent reflectance. Note that, the wavelength dependency could be included in the channel modelling algorithm if available.

IV. MEASUREMENT SCENARIOS AND SETUP
To validate our new channel modelling approach, measurements were performed in two scenarios; i) empty room and ii) conference room. In each scenario, at first, we performed LIDAR measurements for capturing the 3D environment, then modelled the LiFi channels using our channel modelling method and finally performed LiFi channel measurements as validation. In this section, we describe the LIDAR measurement setup, the Tx-Rx configurations and the LiFi channel measurement setup for both scenarios.

A. EMPTY ROOM SCENARIO
In this scenario, the room has a size of 4.5 x 5.8 x 3.1 m³ and is empty with only walls and has no furniture (see Fig. 3). The LIDAR scanner scanned the room at three positions as shown in Fig. 3(a), Fig. 3(b) and Fig. 3(c) in less than 3 minutes per scan. After the scan, offline processing was performed to yield reflectance. Finally, the point cloud data were merged as shown in Fig. 4. This 3D data set is further used for the channel modelling algorithm.
To validate our channel modelling approach, we also performed channel measurements in the same room as shown in Fig. 5. The measurements were conducted in a downlink configuration. At first, where we performed a 4x4    Fig. 5. Note that, the measurement of Rx1 and Rx2, whose positions are marked as P1 and P2 in Fig. 5, was done separately at a different time.
In the mobile user scenario, using the same measurement setup, measurements were done at 40 different positions where Rx1 moved at the edge of a 2 m x 2 m square while Rx2 was kept fixed at the center position of the room. Fig. 6 shows the DC optical power distribution in the room at 1 m height. The positions of Rx1 and Rx2 are marked by blue and green color, respectively. Note that, positions of transmitters are the same as shown in Fig. 5 for both scenarios.

B. CONFERENCE ROOM SCENARIO
In the conference room scenario, the room contained 6 tables, 12 chairs, and a mobile table carrying the channel sounder. Here, the room size was 5.8 m x 6.8 m x 3.1 m. Since the scenario contained more reflecting surfaces, we placed the LIDAR scanner at five different positions, where position 1 is shown in Fig. 7 and all others are marked as green crosses in Fig. 8. After the LIDAR data processing, all 3D point cloud data were merged as shown in Fig. 9 and used as input to the channel modelling algorithm.
In the same room, we performed LiFi channel measurements for a 12x6 distributed MIMO downlink scenario, which comprises of two measurements with 6 Txs each, as shown in Fig 10. Four Rxs were kept on the tables at 75 cm height and two Rxs in the middle of the room at 1 m height. All positions of transmitters are marked as red color in Fig. 8. Transmitters were separated at 1.5 m distance to realize a homogeneous coverage. Note that, transmitters were fixed at 2.9 m height pointing downwards and receivers are pointing upwards.
In this room, we considered three different multiuser (MU) MIMO scenarios. In scenario 1, all receivers were kept at the positions as shown in pink color in Fig 8. In scenario 2, we considered two receiver pairs. Rx1 was moved towards Rx2 and then kept nearby. Similarly, Rx4 was moved towards Rx3 and then kept nearby. In this scenario, receiver pair Rx1 and Rx2 had better coverage from Tx1 and the second receiver pair Rx3 and Rx4 were under the coverage of Tx3. Finally, in scenario 3 receivers Rx1, Rx2, Rx3 and Rx4 were kept close to each other under the coverage of Tx1. In all scenarios, the positions of Rx5 and Rx6 were at the same place as shown in Fig 8.

C. LIFI CHANNEL MEASUREMENT SETUP
The characterization of MIMO LiFi channels has been performed using the LiFi channel sounder developed at   Fraunhofer HHI [48]. It can perform broadband 8×8 channel measurements at frequencies of up to 250 MHz. Our optical frontends are described and characterized in [45]. The transmitter side consisted of a 16-bit, 8-channel arbitrary waveform generator (AWG) (Spectrum DN 2.662-08), in which the first channels were connected to optical Tx frontends. At the receiver side, the optical signals were received by optical Rx frontends and connected to the first channels of the 8-channel digitizer (Spectrum DN 2.445-08). algorithm.
A multi-carrier approach, also denoted as direct current (DC) biased orthogonal frequency-division multiplexing (OFDM), was used for simultaneous measurement of MIMO LiFi channels at all modulation frequencies. The OFDM waveforms were generated using Matlab and transferred into the memory of the AWG. The corresponding signals were simultaneously sent from the parallel AWG outputs to the optical frontends in which a DC is added to modulate the LED around a certain bias resulting in the DC-OFDM waveform.
The transmitted waveform has a similar structure like a physical layer frame in the ITU-T G.9991 standard [49]. It consists of three parts; (a) framing sequence (FS) for detection of the frame start, (b) training sequence (TS) for the estimation of the channel frequency response and (c) OFDM symbols corresponding to the payload data.
After transmitting these packets through the wireless channel, the optical signals were detected at the optical receivers where they are converted into the electrical domain. Then, the DC bias was removed and the signal was digitized for the purpose of channel estimation. Based on FSs defined in [49], we used frame detection by a modified autocorrelation method [50]. Subsequently, the OFDM blocks were processed individually, followed by the de-multiplexing via a DFT. The first OFDM symbols are TSs containing the known complex-valued pilot symbols on a grid of supported subcarriers only, where the used grid depends on the transmitter index. Based on these known symbols, channel estimation on the supported sub-carriers and channel interpolation between them was performed for each pair of Tx and Rx [51]. Accordingly, the frequency response of the MIMO LiFi channel was recorded through simultaneous transmission and detection [52].

V. RESULTS
In this section, we compare channel modelling and channel measurement results for various Tx-Rx configurations in the empty room and in the conference room scenarios. Moreover, we compute the throughput as a metric in exemplary mobile scenarios and show how useful the model is for performance evaluation in complex MIMO scenarios intended for LiFi.

A. EMPTY ROOM SCENARIO
At first, we consider two simple SISO channel scenarios to check the accuracy of the channel modelling methodology. In these SISO scenarios, we model the channels based on varying simulation parameters and compare them with experimental results. In this way, we obtain the parameter values used for the rest of the simulations. Finally, we model the channels and report the results for a 4x4 MIMO scenario and a mobility scenario.

1) SISO
Most of the optical channel modelling parameters have more significant influence for NLOS channels than for LOS channels. There are two major simulation parameters that must be optimized to obtain nearly perfect modelling results. These are i) the number of first reflections , and ii) the spatial resolution of the point cloud data ∆ . By varying each parameter, all simulations were performed for two NLOS scenarios with dominant first-order reflection, for which measurements have already been reported in [45]. Note that the indoor scenario reported in [45] has been the same room where the LIDAR scan was now performed. We compared the simulated channels with regard to the number of first reflections modelled. We increased from 1 to 3 and calculated the channel responses as shown in Fig. 11. Note that here we considered ∆ = 5 cm. For the 3 m scenario in Fig. 11(a), we achieved an MSE of 10%, 4% and 3% for = 1, 2, and 3, respectively, compared to the experimental result. Similarly, for the 2.5 m scenario in Fig. 11(b), the MSE was 2.6%, 1.1% and 0.8% for = 1, 2 and 3. As expected, as the reflection order increases, simulation and experimental results agreed better. However, the simulation time increased significantly as increased. There is a fundamental trade-off between complexity and accuracy of the results. Finally, we considered the impact of varying the average resolution of the point cloud data ∆ for the values 1 cm, 3 cm, and 5 cm. To limit computational complexity, we modeled only first-order reflections precisely using the FD method, i.e. = 1. We observed for the 3 m scenario (see Fig. 12(a)), the MSE = 10%, 8% and 7% for ∆ = 1 cm, 3 cm, and 5 cm, and for the 2.5 m scenario (see Fig. 12(b)), MSE = 8%, 7% and 5%, respectively. Disregarding the insignificant differences between the results, we observed a possibility for overestimating the reflected signals and hence a more frequency selective channel for decreasing resolution of the point cloud data. Additionally, when the resolution of the point cloud is small, the matrix dimensions in the equation (6) are larger and, hence, the computational effort. For very dense point cloud data, it is very complex to compute the channels with more than one reflection. Overall, from these results, we chose the parameters = 3, and ∆ = 5 cm to model the channels for the rest of our simulations, which we consider the best trade-off between complexity and accuracy of the results.

2) MIMO
The measurement and simulation results for a 4x4 MIMO channel measurement are shown in Fig. 13. Here, the bold and dashed lines denote the measured and simulated responses, respectively. Fig. 13 (a) shows the channel at Rx1, placed at position P1 in the middle between Tx3 and Tx4. There are strong LOS signals from Tx3 as well as from Tx4 and weak signals from Tx1 and Tx2. Since all transmitters were kept in a 2 m x 2 m grid, the simulated channel responses for Tx3 to Rx1 and Tx4 to Rx1 are the same. In the experiment, however, due to small differences in the optical frontends, wires and connectors, measured responses differ slightly from each other and do not overlap perfectly like in the simulation. The channels with respect to Tx1 and Tx2 have lower signal strength and measurements were more affected by the noise. Therefore, we observed measured results fluctuate more than in the simulation. As a measure of accuracy, we calculated the relative MSE between measurement and simulation for all links between all Rx and Tx. The MSE of links from Tx1, Tx2, Tx3, Tx4 to Rx1 are 22%, 25%, 1.1%, 2.7%, respectively. Fig. 13 (b) shows the channel responses for Rx2, placed in the middle between Tx1 and Tx2. It can be observed that Rx2 has strong LOS signals from Tx1 and Tx2 and weak signals from Tx3 and Tx4. In the simulated channel responses, the links between Rx2 to Tx1 and Rx2 to Tx3 have the same response. Similarly, links between Rx2 to Tx2 and Rx2 to Tx4 are similar. As explained before, due to mismatch in the optical frontends and other connectors, there were minor deviations in the measurement results, which are not identical to those in the simulations. Here the MSE of links from Rx2 to Tx1, followed by links to Tx2, Tx3, and then Tx4 are 1.5%, 1.3%, 26%, and 29%, respectively. Fig. 13 (c) shows the channel responses for Rx3, placed in the middle between Tx2 and Tx3. For Rx3, signals from Tx2 and Tx3 had strong LOS and signals from Tx1 and Tx4 were weak. Here the MSE of links from Rx3 to Tx1, followed by links to Tx2, Tx3, and then Tx4 are 25%, 2.6%, 1.7%, and 29%, respectively. Finally, the channel responses for Rx4 are shown in Fig. 13 (d). Rx4 placed in the middle between Tx1 and Tx4 had strong LOS signals from Tx1 and Tx4 and weak signals from Tx2 and Tx3. The MSE of the links from Rx4 to Tx1, followed by links to Tx2, Tx3, and then Tx4 are 2.3%, 17%, 21%, and 1.6%, respectively.  From these results, we observed that the MSE is less than 3% for channels with strong signal strength and 10% -30% for weak signals.

3) MOBILITY
In this scenario, we considered a 4x2 MIMO scenario where Rx1 moved at 40 positions around in the room while another receiver Rx2 was kept at a fixed position. Positions of all transmitters and receivers is shown in Fig. 6. We performed MIMO channel measurements and simulations corresponding to all 40 positions of Rx1.
To plot all channels for all frequencies at 40 positions it require larger graph area in the paper. Therefore, to compare the experimental data with the simulations, we calculated the channel gain at a lower frequency of 5 MHz corresponding to each position. Fig. 14 shows the variation of the channel gain versus the distance between Tx and Rx1. The channel gain is plotted for each Tx to the Rx1 link for all positions. When the receiver move away from the center of illumination of one transmitter, the corresponding channel gains are reduced. In the experiment, we observed that the channel gain variation is between -13 dB to -50 dB for distances from 1.85 m to 3.33 m. Due to mismatches in the optical frontends, there are negligible differences in the channel gains at lower distances. When Rx1 is far from the transmitters, the corresponding channel gain is lower and there is a random variation due to noisy data. In the simulations, since all transmitters are placed in a 2 m x 2 m grid, channel gain variations for all Tx links with respect to Rx1 is nearly the same. Therefore, we plotted the average channel gain variation for Tx-Rx1 link. As shown in Fig. 14, the channel gain variations for the simulated channel at 5 MHz varied from -15 dB to -45 dB. When compared with experimental data, the deviation of simulated channel gain was smaller and higher at shorter and longer distances, respectively, due to noise. This study shows that approximate results for channel gain variations can be estimated for a mobile device with low error. From the numerically calculated channel, in the next section we extend our study to estimate channel throughput.

4) THROUGHPUT ANALYSIS
In this section, we extend the study on mobile scenario to calculate the achievable data rate. From the MIMO channel matrices corresponding to all Rx positions we compute further metrics such as singular values and channel throughput. For estimating the throughput, we followed the same approach as given in [52].
Based on the modified Foschini formula [53], [52], the throughput at each position is given as where is the number of considered frequency points, is the MIMO channel matrix, Δ is the bandwidth covered by each subcarrier, = 10 is an empirical scaling factor taking into account impairments like non-linear distortions (clipping) and imperfect constellation shaping, is the total Tx power at each Tx, is the noise power. As explained in [52], the Tx power-to-noise ratio (PNR) can be expressed as where is the average path loss of the MIMO matrix , and is the received power. Finally, the PNR has substituted in equation (12) and the channel throughput corresponding to each user position is calculated.
The total throughput for both receivers along the track is shown in Fig. 15. Note that by substituting equation (12) in (11) the average path loss is normalized out, and thereby any changes in the SNR, which is assumed fixed and equal to 20 dB. Because of (12), only the singular values of the MIMO channel matrix have an impact on the throughput results. Rx1 is moving as follows: first from A to B, then B to C, next C to D and finally D to A: Here A, B, C, and D are the 2D positions of transmitters Tx1 to Tx4 as shown in Fig.  6. Data rate is lowest when both receivers are close to each other and Rx1 receives the minimum power. Intermediate Minima occur always if the Rx1 is at points A, B, C or D. Data rates are the highest if Rx1 is in between two Tx positions, i.e. between B and C or between C and D. In the measurement, due to noise, we observe slight overestimation of the throughput at some points. Overall, the measured and simulated throughput results agree very well. From this results, we can say that, our channel modelling methodology can be used also for the performance analysis of LiFi systems.

B. CONFERENCE ROOM SCENARIO
In this section, we report the results for 12x6 distributed MIMO configuration in a larger conference room. We consider three different MIMO configurations as mentioned in section IV B.
In the previous cases, we focused on the matching between simulated and measured channels. In this section, we expand our study to provide insights into the MIMO channel characteristics and estimate the throughput for different MIMO scenarios.
At first, using the 3D point cloud data, all LiFi channels are modelled for each MIMO configuration. Second, all channels are measured using the channel sounder system [48]. Finally, from both simulated and measured channels, we calculate the singular values and the throughput same as explained in [52]. From a previous study, it is known that the MIMO channel matrix and singular values are dependent on the positions of receivers relative to the transmitters and between each other. Therefore, the positions of optical frontends in the infrastructure have a significant impact on the achievable throughput [52]. Fig. 16 shows the normalized singular values for each MIMO configuration. Here, the bold and dashed lines denote the measured and simulated singular values. Note that for all considered MIMO configurations, positions of receivers Rx5 and Rx6 are the same as shown in Fig. 8 and receivers Rx1, Rx2, Rx3 and Rx4 are always on the table. In all scenarios, receivers Rx6 and Rx5 are closer to the transmitters and have higher signal strength compared to other receivers.
In scenario 1, where all receivers were well separated from each other (see Fig. 8), there are six dominant singular values as shown in Fig 16(a). This indicates the potential of six parallel data links for wireless communications. Some of the singular values have almost the same values indicating nearly orthogonal channels with similar strength, between the ceiling infrastructure and mobile devices. Since Rx5 and Rx6 are located closer to the infrastructure, the first two singular values are higher compared to the other four singular values.
Singular values in the scenario 2 are shown in Fig. 16(b). In this scenario, we considered two receiver pairs, i) Rx1 and Rx2 and ii) Rx3 and Rx4. Here, receivers Rx5 and Rx6 were in the same position as in the previous scenario. Therefore, there are four significant singular values available in this scenario, which show that up to four parallel data streams are useful. At the same time, the first data stream is received by Rx1 and Rx2, the second data stream by Rx3 and Rx4 and the final two streams separately by Rx5 and Rx6. Hence, the channel throughput may be lower as compared to the previous scenario.  Finally, in scenario 3, all receivers on the tables, i.e., Rx1-Rx4, were kept close together and receivers Rx5 and Rx6 are at the same position as shown in Fig 8. Singular values are plotted as shown in Fig. 16(c). In the simulation, there are three dominant singular values, which indicate up to three data streams could be used to communicate in parallel. In the measurement, the fourth singular value is higher than expected, possibly due to orientation error of the optical frontend. The dominant singular values differ more from each other than in the previous scenarios, because the four receivers are kept closer together, so that the correlation of the channels may be high. Only one data stream can be transmitted in the same time slot from the infrastructure to the receivers Rx1 to Rx4, i.e., users can only be served sequentially by using time-division multiple access. The other two data streams are attributed to the strong links between Tx8 and Rx5, and Tx5 and Rx6, respectively. Thus, in this scenario, the throughput becomes rather low compared to previous scenarios.
In Fig. 16(b) and (c), there are singular values that are nonzero but very small i.e. below one, which do not contribute to the throughput at low SNR and maybe used for data communication only at very high SNR. This is visible from comparing the throughput results in the scenarios, which is shown in Fig. 17. Overall, the throughput is reduced from scenario 1 to 3. At low SNR, only the largest singular values matters, and the curves converge. However, due to the differences in the singular values, the curves diverge at higher SNR. Also in scenario 3, the slope increases at high SNR, which is due to the contribution from the lower singular values.
We evaluated the throughput versus the SNR varying from 5 dB to 25 dB. The throughput is calculated according to the equation (4) in [52]. In Fig. 17, bold lines show the experimental results and dashed lines are simulation results. In the simulation results, the channel throughput for scenario 1 varies from around 344 Mb/s to 4.3 Gb/s, for scenario 2 from 332 Mb/s to 3.6 Gb/s and for scenario 3 from 310 Mb/s to 3.1 Gb/s. From the experimental data, the calculated throughput for scenario 1 varies from 346 Mb/s to 4.4 Gb/s, for scenario 2 from 331 Mb/s to 3.4 Gb/s and for scenario 3 from 305 Mb/s to 2.8 Gb/s. Overall, from the results, we observe that the channel throughput is high when all users are spatially well separated from each other and low when all users closer together. As a measure of accuracy, we calculate the relative MSETP between measured and simulated throughput results. We observe that, the MSETP results for scenario 1, scenario 2 and scenario 3 are 0.02%, 0.08%, and 0.4% respectively. The estimated throughput results show a complex relation between Rx and Tx configuration, singular values and SNR, similar to our previous study [52]. In all scenarios, the distance between Rx and Tx has an impact onto the channel throughput.
Note that due to mismatch in the optical frontends and noise in the experimental setup, there are always slight differences between experimental results compared to simulations. However, the overall throughput results show that the simulated and measured throughput are well matched and yield MSETP below 1%.

C. COMPLEXITY RESULTS
The computational complexity of our new channel modelling approach depends mainly on i) the number of points in the point cloud data, ii) the number of first reflections m modelled precisely. As the number of points in the point cloud data gets larger, computation time increases. However, by keeping the resolution ∆S in the range of 5 cm (which is equivalent to a time resolution of 0.167 ns), we can realize a point cloud data set which is sufficient for the currently discussed LiFi systems with several 100 MHz bandwidth and compatible with currently available computing power. By considering ∆S = 5 cm, there are around 35000 and 50000 points in the point cloud data of the empty room scenario and in conference room scenario, respectively.
Similarly, as the number of first reflections m increases, computational time increases significantly. To compare the computational time with respect to first reflections m we consider a NLOS scenario in the empty room (see Fig. 5). We consider a receiver Rx2 kept at P2 position and looking towards the wall. Note that in this case, there exist strong NLOS links from Tx1 and Tx2 and weak NLOS links from other transmitters Tx3 and Tx4. Table I shows the comparison of computational time, average MSE and data rate as a function of . Note that, here we model reflections precisely using FD method and include all other reflections with the integrating sphere model. To compare the variations of the results, we calculate the MSE when compared to the experimental result. We observe that, as number of reflections increases, MSE reduces slightly. We calculated the data rate at the Rx2 using the equation (10). Even though data rate are nearby the same, computational time increases significantly as the number of first reflections increases.
In the conference room, LOS links are available in most of the scenarios in the room. Additionally, the walls and objects in the room are far away from the receiver. Therefore, we did not observe significant differences in MSE as well as for the data rate when increasing m. However, computation time is 50 seconds, 3.2 hours and 6.1 hours for = 1, 2 and 3. The steepest increase in computational time is between = 1 and = 2. This is caused by the intrinsic room function ( ) in equation (6) is considered from = 2 onwards. Therefore, by choosing a reasonable tradeoff between accuracy of the modelling results and number of first reflections , our approach can be used to model mobile LiFi channels in a timely manner. In most parts of this paper, a minor modeling error was targeted, for which we have used = 2 as a reasonable compromise. For mobile LiFi network simulations with large numbers of access points and mobile devices, another trade-off at = 1 may be reasonable. This approach would model the LOS links and the first diffuse reflection precisely and include the residual reflections by the integrating sphere model. In this way, computational time can be significantly reduced. For example, in the empty room scenario, computational time is reduced by factor around 180 for a single point on the track (see Fig. 6). The computation time could be well spent to model a dense grid of points in the room. Some very relevant effects, such as blockage of the LOS and the possibility of NLOS communications via a nearby object would then be included, which is a major advantage over the common assumption to model the LOS links with no blockage at all.
Since calculations are based on matrix or vector multiplications, it is possible to accelerate the modelling by using high performance Graphics processing unit (GPU) or Tensor Processing Unit (TPU).

VI. CONCLUSION
In this paper, we propose a new channel modelling method for LiFi. We used a LIDAR scanner to automatically capture the 3D environment and use the point cloud data as an input to our channel modelling algorithm. In the channel modelling algorithm, we used a frequency domain technique to model the LOS and the first NLOS reflections and an integrating sphere model for the rest of the higher order diffuse reflections.
In the 3D modelling using LIDAR, at first, the LIDAR scanner was kept at multiple locations in the room to accurately scan the room. As a pre-processing, at first, we remove the noisy points, secondly calculated the surface normal, and then reduced the resolution of point cloud data. After that, we estimated the reflectance of the surface points based on the method was described and tested in this paper. Finally, all scanned data were merged to get the complete point cloud data of the room, based on which the LiFi channels for the given Tx-Rx configuration can be modeled as described above.
To validate our simulation method, we compared modelling and measurements in two realistic indoor scenarios, an empty room and a conference room. In an empty room scenario, we first obtained the best parameter values used for the rest of our simulations based on the channel modelling results. For a realistic LiFi scenario in a 4x4 distributed MIMO configuration, measurements and modelling results are in a good agreement with less than 3% error for the channels with high signal strength. Moreover, we demonstrated the possibility to model mobile scenarios in the same room, in a 4x2 MIMO scenario where one receiver is moving while another is at a fixed position. Measurement and modelling results are in good agreement for channels with high gain. For the mobile MIMO channel, we calculated the achievable data rate and observed that the modelled and measured throughput results agree very well.
In a further step, we performed channel modelling and measurements in a larger conference room scenario, where the LIDAR scanner was kept at five different locations to capture all objects including blockages from all necessary points of view. We modeled the LiFi channels for a 12x6 distributed MIMO scenario with three different MIMO configurations. From the channel data, we evaluated the singular values and channel throughput. The overall throughput results show that the simulated and measured throughput results agree very well and yield an error below 1%.
Overall, our approach can be used to model LiFi channels in complex environments including multi-story buildings, apartments, auditoriums, hospitals etc. In the future, it is possible to integrate our channel modelling algorithm with the LIDAR pre-processing software to model the LiFi channels in many indoor scenarios. As a future vision, it is possible to create the 3D map of a coverage area, determine how many access points are required and evaluate the SNR and throughput of the LiFi systems. Additionally, using advanced algorithms, it may be possible to predict the channels for dynamic movements and evaluate the performance with minor error. One of the major applications will be the performance analysis of PHY and MAC layer algorithms in research and standardization. Moreover, our technique can be used for LiFi network planning in real scenarios, including cost-versus-performance optimization.

A. INTENSITY CORRECTION RESULTS
Using the intensity correction method described in the section II B, we estimate the reflectance of walls and objects in a 3D room from the raw intensity data measured by the LIDAR. In this section, we report the results for two considered rooms, an empty room and a conference room.
Based on the LIDAR measurements in the empty room (see Fig. 3), intensity data are obtained along with the point coordinates. Fig. 18 (a) shows the raw intensity data of the room obtained from our measurement. To obtain the reflectance parameter of each point, the raw intensity data need to be corrected. From this raw intensity data along with considered 3D coordinate positions of LIDAR scanner in the room, we correct the intensity and estimate the values as shown in Fig. 18(b). These corrected intensity values of each point in the point cloud data are considered as reflectance values for rest of the simulations in this empty room.
Similarly, intensity correction method is applied for the conference room scenario. The LIDAR measurement setup in the room are mentioned in section IV. Scanning locations of LIDAR scanner in this room is given in Fig. 8. The raw intensity data of the conference room scenario obtained from our measurement is shown in Fig. 19 (a). After the correction method, the corrected intensity is shown as Fig. 19(b). These corrected intensity values of each point in the point cloud data are considered as reflectance values for rest of the simulations in this conference room.

B. VISIBILITY ANALYSIS USING HPR METHOD
The HPR method consist of two main steps: inversion and convex hull construction. At first, an inversion of the point cloud data is performed using the so-called spherical flipping method. In the next step, for the given view point, the convex hull construction of the inverted point cloud data set is performed. We follow all these steps as mentioned in [42]. To calculate the visibility from the given view point, we use the same MATLAB code as given in Section 5 in [42]. Besides the simplicity of the implementation, major advantages of this method are to i) determine the visibility without reconstructing a surface and ii) calculate the visibility for dense as well as sparse point clouds. By using this approach, we calculate the visibility factor in each optical link between Tx, Rx and surface elements.
In the HPR-based visibility analysis, it is important to select a suitable radius parameter, which can be defined as 10 [42]. As increases, more points are marked visible and when is small fewer points are marked visible [42]. Therefore, there are always some false points being marked wrongly as visible or not [42]. Fig. 20 shows an example plot for visibility analysis for varying values for . We consider the point cloud data of an empty room scenario (see Fig. 18(b)). In this example, we calculate the visibility between Tx location to all the points in the point cloud data. Note that, Tx is kept nearby one corner of the room (marked as red color) and is looking towards the next corner of the room. Due to the mobile table, which is kept in the middle of the room (see Fig. 5 or Fig. 6), there will be blockage between Tx to some of the points in the point cloud data. Therefore, there is a possibility that some of the points should not be visible while looking from the Tx point. As shown in Fig. 20(a) many points are not visible from the Tx location when = 1. In the Fig. 20(b), when is increases to 2 many more points are visible and only few points are not visible. However, when = 3 (Fig. 20(c)) and = 4 ( Fig. 20(d)) most of the points are visible even there is a blockage present due to the mobile table. We observe that, large lead to more points that are visible and lower values of lead to very less visible points, hence the wrong estimation of channel. From this study, we observe that results are better when lies in between 1.5 to 2. Although the minimum is not the same, indicating the best choice, we consider the value for = 1.5 for our rest of the simulations. Even though the value of R is considered as 1.5, still there will be fewer false points, which do not affect significantly on the LiFi channel modelling.
When the point cloud is very noisy or non-uniformly sampled, a robust HPR operator (RHPR) can also be used to evaluate the visibility as reported in [43]. Note that, instead of HPR method, we can also use any other method to evaluate the visibility of the point cloud data [44]. From our study, we observe that, the HPR method estimates visibility for our LiFi channel modelling application with the lowest computational time and reasonable error.