1. Introduction

In the last years, indoor wireless optical communications have drawn the interest of researchers due to the need of high-speed and inexpensive data links, overall in environments where radio frequency links are not always viable. However, the communication capacity depends on the characteristics of the room where this communication is established [1]. The multipath penalty, which limits the maximum baud rate, can be characterized by the impulse response. In order to evaluate the impulse response on indoor wireless optical channels, several deterministic methods were firstly proposed [2]. These methods can only be implemented to determine the impulse response until the third reflection due to their computational complexity. Later on, a modified Monte Carlo-based ray-tracing algorithm (MMCA) was introduced, which presents a lower computational cost and without limit in the number of reflections to consider [3, 4]. The algorithm allows easily for the introduction of new models for the devices and materials, such as the Phong’s reflection model [5,6], or the increase of the complexity of the simulation environment (including obstacles and rooms with complex geometries), and it also allows to estimate the error committed by the method during the computation of the impulse response [7].

Recently, a new ray-tracing method, which has been denominated photon-tracing algorithm (PTA), has been proposed [8]. This method is based on the MMCA but now the rays are not always propagated after reflection. Instead, they are only propagated in a certain proportion according to the reflection coefficient, which makes the number of considered rays decrease rapidly after each reflection. Zhang et al. proved the reduction in computational cost of the new re-designed algorithm, but they did not comment on the accuracy of the results. In this paper, we compare in more detail the differences between the two algorithms, MMCA and PTA, according to error analysis and computational cost.

The paper is organized as follows: Section 2 describes the Monte Carlo ray-tracing algorithm, highlighting the differences between MMCA and PTA. The computational cost is evaluated in section 3. The equations that allow us to determine the error in computing the impulse response are presented in section 4. Finally, several computer-simulated results are reported in section 5 in order to compare both methods. The conclusions are summarized in section 6.

2. Algorithms description

In the modified Monte Carlo ray-tracing algorithm, ray directions are randomly generated according to the radiation pattern from the emitter. The contribution of each ray from the source or after a bounce to the receiver is computed deterministically. Consequently, the discretization error is due to the number of random rays. The line-of-sight (LOS) and multiple-bounce impulse responses are considered when calculating the total impulse.

2.1. LOS impulse response

Given an emitter E and receiver R in an environment without reflectors, with a large distance d between both, the received power is approximately

Here n is the mode number of the radiation lobe which specifies the directionality of the emitter, P_E the power radiated by the emitter, A_r the physical area of the receiver, and FOV the receiver field of view (semi-angle from the surface normal).

2.2. Multiple-bounce impulse response

We consider now an emitter E and receiver R in a room with reflectors. The radiation from the emitter can reach the receiver after any number of reflections (see Fig. 1). In the algorithm, many rays are generated at the emitter position with a probability distribution equal to its normalized radiation pattern R_E (ϕ,n)/P_E. The power of each generated ray is initially P_E/N, where N is the number of rays used to discretize the source. In MMCA, when a ray impinges on a surface, the reflection point is converted into a new optical source, thus, a new ray is generated with a probability distribution provided by the reflection pattern of that surface, R_S(θ,θ′). The process continues during the simulation time. After each reflection, the power is reduced by the reflection coefficient of the surface, and the reflected power reaching the receiver is computed by

 

Fig. 1 Geometry of emitter and receiver with reflectors. Reflection pattern of the surface is described by Phong’s model.

Download Full Size | PDF

In PTA, when a ray impinges on a surface, the probability that a new ray is generated is given by the reflection coefficient ρ. Moreover, when a new ray is generated, its direction of propagation is determined according to the probability distribution of the reflection pattern as in MMCA, but now the power of the new ray keeps unalterable and equal to P_E/N. In the same way, after each reflection, the reflected power reaching the receiver is computed by using Eq. (5), but with P_i = P_E/N in Eq. (6). Therefore, both algorithms perform similarly because in PTA only ρN_k rays of power P_ray = P_E/N are reflected (being N_k the number of rays remained after the previous kth bounce) whereas in MMCA all the rays are reflected but with output power P_ray = ρP_i. As we will see, the fact that, in PTA, not all the incident rays are reflected leads the number of computational operations to decrease rapidly with each new bounce with respect to MMCA.

3. Computational complexity

The main advantage of Monte Carlo ray-tracing algorithms is that they allow for the evaluation of the impulse response for rooms of complex geometries [7], in contrast to other methods [2,5], without a meaningful increase in the computational cost. This can be explained by the number of elementary calculations that are performed: $N_{\text{op}}^{\text{MMCA}}=KN{N}_S$, where N is the number of rays, K is the number of reflections that are considered, and N_S is the number of surfaces that define the room. An elementary calculation is defined as the calculation of power contribution and delay from a point source (emitter or reflection point of a ray) to the receiver, as in Eq. (5), or the assessment of the propagation of the new generated ray to determine a new point source. The previous value of $N_{\text{op}}^{\text{MMCA}}$ is an upper bound, because not always, after a reflection, the N_S − 1 remaining surfaces have to be considered as possible future reflecting surfaces (in rooms with complex geometries, some surfaces can be placed on other greater surfaces, e.g. windows or doors over walls; therefore, when a ray sets up from a surface with others superimposed to it, all of them can be discarded during the searching of the new collision point). Moreover, not always the ray contributes in the received power, because sometimes the emitting point is out of the FOV of the receiver and its contribution does not have to be computed. Similarly, not all the rays reach the maximum number of reflections K during the simulation time t_max either. However, this value represents a good approximation to the required number of elementary calculations of MMCA.

In PTA, if we consider that N rays (photons) are initially launched, after the kth bounce, only ρ̃_kN_k−1 rays (photons) continue their path, where N_k−1 is the number of photons remained after the (k − 1)-th bounce (with N₀ = N), and ρ̃_k is an average parameter of the reflection coefficient at the kth bounce which depends on the reflection coefficients of the surfaces, but also on the radiation and reflection patterns, and the position and other characteristics of the emitter. The number of elementary calculations can be computed as:

It must be noticed that the last rays considered to compute their power contributions are those after the (K − 1)-th reflection. These last rays are propagated to determine the reflecting surfaces (and the reflecting points) and then their power contributions are computed. This is the reason why the previous equation has been truncated in the N_K−1 term.

Evidently, ρ̃_k depends on which bounce is being considered. However, we can assume an average reflection coefficient ρ̃ which allows us to represent any bounce obtained as an average over all the ρ̃_k. Then, the previous series is reduced to N_S × (N + ρ̃N + ρ̃²N + ... + ρ̃^K−1N). This numeric series converges to:

As an example, we are going to consider the room studied in Section 5. The average reflection coefficient is ρ̃ = 0.69 (only taking into account the areas of the surfaces), the number of rays is N = 500000, the number of surfaces is N_S = 6 and the number of considered reflections K = 10. With MMCA, the number of elementary calculations is upper bounded by $N_{\text{op}}^{\text{MMCA}}=30\times {10}^6$, whereas in PTA, according to Eq. (8), we will only need $N_{\text{op}}^{\text{PTA}}=9.3\times {10}^6$ elementary calculations, that is, the 31.1% of those required by MMCA. It must be noticed that the previous results are only crude approximations, but they allow us to state that, from a computational point of view, PTA is much more efficient than MMCA.

4. Error estimation

Both PTA and MMCA are Monte Carlo based ray-tracing algorithms. Therefore the error analysis for MMCA given in [7] is directly applicable to PTA. There, it was demonstrated that the error in computing the power P′_j reaching the receiver during a small time interval Δt can be estimated from the variance of P′_j, var (P′_j) The biggest admissible Δt is defined as the largest interval which ensures that the same ray does not contribute twice to a receiver near the walls, being var (P′_j) given by [7]:

In the previous equation, N_f,j is the total number of flights during the jth time interval (rays flying during that interval), N_j the number of rays that contribute in the power reaching the receiver during Δt, and p_i,j is the power contribution of the ith ray (i = 1, 2,...,N_j) arriving during that interval obtained by using Eq. (5). Therefore, P′_j is the power value in the jth histogram interval computed with the Monte Carlo ray-tracing algorithm. In MMCA, N_f,j coincides with N, because during a certain interval j there are exactly N rays flying, those generated by the emitter, which are never discarded. However, in PTA, after each reflection a certain number of rays are removed according to the reflection coefficient of surfaces, then N_f,j decreases exponentially along the time, being different for each time interval (which has been indicated by the sub-index j in N_f,j : N_f,j1 ≥ N_f,j2, j₁ < j₂).

The relative error can be computed as the square root of the variance (one standard deviation) of P′_j given by Eq. (9) divided by the computed contribution power defined as

Then

The above equation allows us to estimate the relative error in computing the received power in the jth time interval.

5. Simulation results

In order to compare PTA and MMCA methods, the simulation room described in [4,8] has been evaluated. The main characteristics of this room are indicated in Table 1. In this example, the reflecting surfaces are purely diffuse Lambertian (r_d = 1), but rooms with materials characterized by the Phong’s model (such as the examples described in [7]) have also been evaluated, obtaining similar results to those presented below. The emitter is placed at the center of the room rested on the ceiling and pointing straight down. The receiver is located on the floor near a corner pointing straight up. The maximum number of considered reflections is K = 10 and the simulation time t_max = 120 ns.

Table 1. Parameters for simulation

View Table | View all tables in this article

In Fig. 2 we show the impulses responses (received power) obtained by using both methods when N = 500000 rays (photons) are generated from the emitter: in Fig. 2(a) we have the response computed by the PTA method, whereas that obtained by MMCA is depicted in Fig. 2(b). We can observe how very similar impulse responses are obtained with both methods, only distinguished by a slight greater ripple in that provided by PTA.

 

Fig. 2 Impulse responses obtained with (a) photon-tracing (PTA) and (b) ray-tracing (MMCA) algorithms.

Download Full Size | PDF

The accuracy of both methods can be compared quantitatively by means of the relative error committed by both algorithms during the calculation of the impulse response. Figure 3 shows the relative error for PTA and MMCA calculated by using Eq. (11). It can be observed how MMCA offers a solution (impulse response) with a relative error roughly constant along the time. However, in PTA the relative error tends to increase with time, what is logical because the number of rays (photons) that contribute in the calculated received power decreases along the time, since many of them begin to be discarded when they collide against the room surfaces. This is understood clearly if we observe the number of rays that contribute in the received power along the time or at a certain bounce k (see Fig. 4). One can see how in MMCA the number of contributions always presents a relative high value (> 1/2 peak value) except at extreme time instants (t > 80%t_max) when the rays are moved apart from each other and their contributions become more and more spread. In addition, no more rays are generated after the kth reflection. On the contrary, in PTA the number of contributions decreases rapidly (in an exponential way) with time or with the bounce index. Therefore, the received power (impulse response) is computed by using a lower and lower number of information samples, leading the algorithm to present a greater relative error at longer time instants.

 

Fig. 3 Relative error of both algorithms: PTA (blue) and MMCA (green).

Download Full Size | PDF

 

Fig. 4 Number of contributions (a) along the time and (b) at the kth reflection.

Download Full Size | PDF

In Table 2 we compare PTA and MMCA according to mean relative error and number of elementary calculations. The mean relative error has been weighted by the value of P′_j:

Table 2. Comparison between PTA and MMCA

View Table | View all tables in this article

The number of elementary calculations required for both algorithms (see Table 2) demonstrates that PTA is much more efficient than MMCA (∼28% out of the total of operations required by the latter), as we had stated in section 3. We can observe how the predicted values in that section are very close to those obtained during the simulations. Therefore these expressions can be used to determine in advance the number of required calculations approximately. In comparison with MMCA using N = 500000, the simulations performed with the PTA for N = 1000000 and N = 1500000 needed the 56% and 84%, respectively, of elementary calculations. However, for the latter we obtained a relative error inferior to that of MMCA when t < 50 ns in spite of requiring 16% less computation time.

Finally, we show in Table 2 the power distribution according to the bounce index. We can observe how both algorithms present very similar values, demonstrating a good behavior of the PTA during the distribution of power in the remained rays after each reflection.

6. Conclusions

In this paper, we have compared the conventional modified Monte Carlo ray-tracing algorithm with a recently proposed one, which has been called photon-tracing algorithm. We have established quantitative parameters to carry out this comparison according to computational cost and accuracy of the provided solution. We have stated analytically and by means of simulation results that PTA presents a lower computational cost than conventional MMCA. However, regarding the error committed by both algorithms, MMCA is more reliable than PTA, although the mean relative error of the latter can be considered acceptable taking into account the reduction in computing time. In addition, more rays can be used by the new method still requiring lower simulation run-time in order to obtain more accurate results. Therefore, we can conclude that PTA begins to appear as a good substitute to MMCA with superior performance.