Fast Ray Tracing for Microcellular and Indoor Environments
Fernando Aguado Agelet, Fernando P? rez Font? n e a
ETSI Telecomunicaci? n, Universidad de Vigo, 36200 Vigo, Spain o
Universit¨ t des Saarlandes, FB 14 Informatik, 66041 Saarbr¨ cken, Germany a u
—We present a very ef?cient algorithm which calculates deterministically the path losses, delay pro?les and other characteristics of a mobile channel. The algorithm can be used both for indoor and outdoor environments. We apply the technique of tracing rays and compute re?ections, transmissions and diffractions with the help of geometrical optics, heuristical approximations and the uniform theory of diffraction (UTD). I. I NTRODUCTION In principle, the characterization of the mobile channel could be calculated with the help of the electromagnetic ?eld theory based on the equations of M AXWELL while considering the bounding conditions present as physical properties of the walls and other structures in the environment. However, this approach is practically unfeasible and numerical methods approximating the real characteristics have become popular. One of this methods characterizes the mobile channel using the ray tracing technique and approximating the ?eld conditions by applying simpli?ed models, e. g., uniform theory of diffraction UTD. Several publications report suf?cient similarities between calculated results and measured values , . They state that the number of re?ections, transmissions and diffractions is an important factor for the accuracy of the simulations, because especially digital systems are more sensitive to channel degradation in form of multipath distortion and fading. We present a fast algorithm which uses ef?cient data structures and programming paradigms from computational geometry to ?nd all ray paths between two arbitrary points both in microcellular and indoor environments. We discuss the run time behavior of the algorithm based on a theoretical analysis as well as on experimental results. We apply the 2D-3D model as described, for instance, in , that is, all ray paths are calculated in the 2-dimensional plane while assuming in?nite vertical edges. The run times of the algorithm even with high order re?ections, transmissions and diffractions are acceptable for practical usage.
II. C OMPUTATIONAL G EOMETRY PART The most important interactions of an electromagnetic wave take place at obstacles close to the point of interest, e. g., at the walls and edges of the same street or crossing where the emitter is placed. Hence, the geometrical data structure of a quadtree is an appropriate means to exploit the locality. In order to calculate deterministically the multipath visibility between emitting and receiving point, we apply both recursively and iteratively the paradigm of a sweep line taken from computational geometry. A further improvement of the run time is achieved by storing intermediate results in an ef?cient dictionary. Thus, multiple re?ections between diffraction points are calculated only once. In the sequel, we ?rst present some basic de?nitions and notations. Then, we introduce the quadtree. The following part reviews a priority queue which is the basic data structure to maintain the visibility property during the sweep. The polar sweep is explained next and the section concludes with the overall algorithm. At the beginning let us exclude transmissions. We will see at the end that they can be included easily. The program has been implemented with the library LEDA which provides a collection of very ef?cient data structures and algorithms . A. De?nitions and Notations We de?ne a city or an indoor environment by a set of closed polygons. Each polygon is described by an ordered list of segments. Two segments have at most one corner as a common point. Looking from the outside, possible diffraction points are all those visible corners of the polygon which have an inner angle of less than . Looking from inside, it is vice versa. The radiation area of a point P is the part of the plane which possibly can be reached by a wave emitted from P . Thus, the radiation area of the point where the antenna is placed consists of the entire plane. Without loss of generality we assume an isotropic antenna. The radiation area of a virtual point V is given by the area which does not contain V and which is enclosed by the re?ecting segment and the two rays starting from V and going through the end points of the re?ecting segment. The radiation area of a diffraction point D is the larger one of the two regions con?ned by the two rays starting at D and following the direction of the adjacent segments. In order to de?ne the visibility between two points, we con-
Manuscript received March 19, 1996. F. Aguado Agelet, e-mail firstname.lastname@example.org, phone +34-86-812122, F. Per? z Font? n, e-mail email@example.com, fax +34-86-812116, A. Formella, e a e-mail firstname.lastname@example.org, phone +49-681-302-5538, world-wide-web http://www.mpi-sb.mpg.de/LEDA/www/esmi.html This work was supported in part by the Xunta de Galicia and by the German Science Foundation (DFG) through SFB 124, TP D4.
city of Vigo
area for coverage
emitter zoom quadtree
Fig. 1. City of Vigo with quadtree
sider only the segments which lie at least partially within the radiation area of the point serving as emitter. The visible segments of a point P are all those segments which lie within the radiation area of P and which are at least partially visible from P . Note, that for closed polygons the backsides can be detected as not being visible by a simple test. A point P is visible from point Q if P lies in the radiation area of Q and in front of the visible segments of Q. A point P is visible with one re?ection from point Q if there exists a visible segment S of Q such that P is visible from the virtual point Q0 as image of Q in S . A point P is visible with one diffraction from point Q if there exists a diffraction point D such that both P and Q are visible from D. It is straight forward to de?ne the visibility of a point with r re?ections or d diffractions respectively. A ray path between two points Q and P is called an (r; d)path, if P is visible with at most r re?ections and at most d diffractions from Q. Note, that due to symmetry an (r; d)-path between P and Q is an (r; d)-path between Q and P as well. Given two numbers r and d and two points, where one serves as emitter E and the other serves as receiver R, the aim of the algorithm is to ?nd all (r; d)-paths between E and R. B. Data Structures The term quadtree is used to describe a class of hierarchical data structures whose common property is that they are based on the principle of recursive decomposition of 2-dimensional space. The decomposition is governed by properties of the input data. Figure 1 shows how the city of Vigo is divided recursively into quadrants of decreasing size until the number of segments is less than a certain constant. The root node of the tree like data structure corresponds to the entire city. Each of the four descendants of a node represents a quadrant. In order to ?nd the leaf where a certain segment belongs to one has to follow the path from the root to the leaf in the tree. The quadtree is built once at the start of the program. Because segments are usually small and do not intersect, the building time of the quadtree is O(n log n), where n is the number of segments in the city. A priority queue is a data structure for maintaining a set S of elements, each with an associated value called a key. We use
a priority queue that supports the following operations. I N SERT(S,x) inserts element x into set S. M INIMUM(S) returns the element of S with the smallest key. D ELETE(S,x) removes the element x from the set S. The run time of I NSERT() and D ELETE() is O(log n) where n is the actual size of S. The run time of M INIMUM() is constant. The priority queue holds the list of segments being intersected by the sweep line starting at a point O (see next section). In our case the key will represent the distance on the sweep line of a segment from O. The operation M INIMUM() needs a <-relation for elements in S . Therefore, we de?ne a partial ordering of segments. For two segments si and sj the relation si < sj indicates that there exists a ray starting at O that ?rst hits si then sj ; if no such ray hits both segments, the segments are unrelated. We insert a segment s into the queue when we ?nd the ?rst handling point of the segment. We delete the segment s if we ?nd the last handling point of the segment. Hence, segments maintained in the queue are always related. The function M INIMUM() returns the segment actually visible from O. A dictionary is a data structure to maintain a linearly ordered set. Insertions are performed by simple assignments to the appropriate array index i. The function I S D EFINED(i) returns as a boolean value if previously to array index i an insertion has been performed. Operations on the dictionary take run time O(log n), where n is the current size of the dictionary. We use such a dictionary to store the all paths with at most r re?ections already computed between two points of interest. C. Polar Sweep Algorithm The term plane sweep is used to characterize a paradigm employed to solve geometric problems by sweeping a line across the plane and halting at points where the line makes the ?rst or last intersection with any of the objects being processed. At these points, the solution is partially computed, so that at the end of the sweep, a ?nal solution is available. We have modi?ed the classical plane sweep algorithm where the line is swept in the horizontal direction. Our line is rotated clock wise around a point O sweeping the radiation area of the point resulting in a polar sweep. Each time a start or end point of a segment is found, the priority queue is updated. By carefully managing the bounding conditions while beginning and ?nishing the sweep, the sequence of minimum segments during the sweep provides an ordered list of partly or entirely visible segments of point O. To apply the polar sweep technique, we need to sort the halting points of the sweep line, i. e., all start and end points of the segments lying in the radiation area of point O are sorted according to their polar angles. In Fig. 2 we have illustrated the method. The overall run time of the polar sweep is O(n log n), where n is the number of segments falling into the radiation area, because both sorting of 2n points and updating 2n times the priority queue is done in this time. In order to further speed up the algorithm, we combined the use of the polar sweep method and the quadtree introduced above. In a ?rst step only
sweep line visible segments
building emitter leaf quadtree
Fig. 2. Example for Polar Sweep
those segments are swept which fall both into the radiation area and into the leaf of the tree where point O is located. If we already have found a closed visibility line, the search for visible segments can be stopped. If the visibility line is not closed, as it is the case in Fig. 2, the remaining radiation areas must be swept in neighboring leaves. In general, the number of leaves to be visited is small compared to the overall number of segments in the environment. The run time within a leaf is constant, it depends only on the number of segments within the leaf. Searching a neighboring leaf is done in O(log n), where n is the number of segments in the environment. D. Overall Algorithm The overall algorithm can be sketched as follows. The input data consist of a description of the environment, the locations of emitter and receiver, the electromagnetic data, and the number of re?ections r, diffractions d, (and transmissions t). After the input data is read, the segments are sorted into the quadtree. This provides an adaptive subdivision of the plane. The algorithm ?nds all (r; d)-paths works in a recursive, thus depth ?rst search manner for re?ections and in a iterative, thus breath ?rst search manner for diffractions. First, the leaf of the quadtree is found, where the emitter is located in. Then, the polar sweep is performed within the leaf and possibly in neighboring leaves (see Fig. 2). If the receiver is visible, a direct ray path has been found. Now, we start the recursive search for second and higher order re?ections. As long as the limit for the number of re?ections has not been reached, we calculate for each visible segment the virtual image of the emitter. This virtual point serves as new emitter and we repeat the process described above. Note, that during the recursion we ?nd all diffraction points reachable with up to r re?ections. This information is stored in the dictionary. Once the recursion has completed, we start processing the diffraction points. Here, we perform an iteration over all diffraction points being visible with zero re?ections from the emitter ?rst. Each diffraction point serves as a new emitter and all (r; 0)-paths towards the receiver are calculated. Note, that this calculation ?lls the dictionary with information about paths from a diffraction point to the receiver. When the iteration over the diffraction points reachable from the emitter has been ?nished, we start processing with those visible with one re?ection and so forth. At this point, the entries of the dictionary become valuable. If we en-
counter a diffraction point which has already been considered in a previous step, the information about all paths consisting of at most r re?ections towards the receiver is already stored in the dictionary. Hence, almost no additional computation must be performed. In other words, the dictionary is ?lled dynamically with all those parts of the intervisibility graph considering up to r re?ections between diffraction points, which are necessary to compute the solution of the given problem. After all 2-dimensional ray paths have been found, the height of the buildings is checked against the line of sight between emitter and receiver and the re?ections on the ground are incorporated. In the end, the required electromagnetical calculations are performed (see Section III.). We have seen, that diffraction points are handled by treating them as virtual emitters in the iteration phase of the algorithm. Similarly, virtual points as images of emitters in re?ecting walls are treated as new sources of radiation. In the case of transmissions, the same method can be used as for re?ections. Instead of calculating a virtual point, the original source is maintained. The radiation area is the area beyond the segment serving as transmitting surface. In the algorithm, the transmitting surface is added twice to the visibility list, once as re?ector and once as transmitter. If more than one receiver is placed in the environment, we start calculating from the receivers towards the emitter. This can be done because of the symmetry properties of the (r; d)paths. The advantage of this reverse calculation lies again in the dictionary. For the second receiver, a lot of diffraction points which have been treated for the ?rst one, are already available in the dictionary. When we calculate the coverage in a certain area, each grid point serves as a receiver and the run time is dominated by calculating simple re?ections, because most of the paths between diffraction points are stored in the dictionary. Other kinds of antennas with different spatial radiation patterns decrease the range of the sweeping angle at the start of the algorithm, thus we expect even smaller run times. III. P ROPAGATION M ODEL In  the equation h(t) = i=1 Ai (t ? i ) exp(?j i ) is proposed to describe the mobile channel as a form of a bandlimited complex impulse response. The transmitted impulse is mathematically described by a Dirac function. The received signal h(t) is formed from the vector addition of n time delayed rays. Each one is represented by an attenuated and phaseshifted version of the original Dirac waveform. Re?ected rays are calculated through the use of geometrical optics. Diffracted rays are calculated using the geometrical theory of diffraction , . According to the objects encountered on the ray path, we compute the individual receiver ray power Pr as
2 Y P r = P t Gt Gr Rj 2 (4 d) j
#2 " Y
As (s ; s)Dl
4000 3500 run time in seconds 3000 2500 2000 1500 1000 500 0 1 2 3 4 5 6 7 number of point 8 9 10 (5,1) path 20 seg (5,1) path 40 seg (5,1) path 60 seg (5,2) path 20 seg (5,2) path 40 seg (5,2) path 60 seg
Fig. 3. Run times for trajectory of 10 points (see Fig. 1) for (5; 1)- and (5; 2)paths with different maximum number of segments in a leaf of the quadtree on a Sparc2 workstation
where Pt represents the transmitting power of the ray, d is the total length of the ray, is the wavelength, Gt and Gr are the transmitting and the receiving antenna gains in direction of the ray, Rj is the re?ection coef?cient for the j th re?ector, Tk is the wall transmission coef?cient for the k th transmission and Dl is the diffraction coef?cient for the lth diffracting wedge. The diffraction wedges are also multiplied by a spatial attenuation function As (s0 ; s) which ?nds the correct multiplicative diffraction coef?cient given the 1=d2 dependence in the ?rst term. All coef?cients are functions of the angle of incidence of the path and the object characteristics. Although this equation represents the instantaneous received power, all rays are stored as complex ?eld vectors which are multiplied by the complex re?ection, transmission and diffraction coef?cients. This approach is used to model the complex baseband impulse response without loss of phase information. The coverage map in Fig. 4 has been computed in that sense. IV. S IMULATION R ESULTS The propagation model allows to calculate the channel response throughout an indoor or microcellular environment. The signal level, RMS delay spread and fading statistics at each point can be obtained . As results, we present some run time data including high order re?ections and diffractions and a coverage map for the complex power received in an area of the city. Clearly, there exists a trade-off between the number of segments in a leaf of the quadtree and the run time. If the number is too small, too many leaves must be examined; if the number is too large, the sweep within a leaf takes too much time. Simulations have shown that approx. 40 segments in a leaf deliver the best run times. The curves in Fig. 3 demonstrate the tradeoff and present absolute run times on a SunSparc 2 workstation. With approx. 20 Mbyte, the memory requirements have been moderate even in experiments with (5; 3)-paths. The receivers have been placed along the line shown in Fig. 1. The coverage map shown in Fig. 4 has been calculated for an antenna as a dipole of length =2. The resolution was set to a grid of 40 by 40 points. The receiving antenna was assumed to be isotropical and the working frequency was set to 1 GHz. Both antennas have been placed one meter above ground. The buildings are modeled with a relative permeability r = 1, a
Fig. 4. Coverage of a region in the City of Vigo. The region corresponds to the highlighted rectangle in Fig. 1. We calculated approx. 160000 3D ray paths. The grid was set to 40 40 points. The grey scale re?ects the relative power respective to the maximum value, which can be observed close to the emitter at the lower left border of the image.
relative permittivity r ? 0m?1 .
= 2 1,
and a conductivity
= 2 43
V. C ONCLUSION The presented polar sweep algorithm exploiting ray coherence by the use of a quadtree data structure shows good performance for an arbitrary number of re?ections, transmissions and diffractions. The maintenance of previously calculated ray paths between points of diffractions allows for fast table look up when many different receiving points are incorporated in the de?nition of the problem. Coverage maps and probability density functions for the spatially dependent ?at fading envelope amplitude can be calculated in high resolution with run times acceptable for practical usage. High order re?ections and diffractions provide an accurate power delay pro?le which is fundamental to determine the inter-symbol interference and bit error rate for digital systems. The method can be extended to include horizontal diffracting edges. In contrast to previous methods (e. g., ) the sweep line technique is suitable for 3-dimensional calculations as well. R EFERENCES
 G. E. Athanasiadou, A. R. Nix, and J. P. McGeehan, “A ray tracing algorithm for microcellular and indoor propagation modelling”, in Proc. of IEE Conf. Antennas and Propagation no. 407, pp. 231–235, April 1995.  S. J. Fortune, D. M. Gay, B. W. Kernighan, O. Landron, R. A. Valenzuela, and M. .H. Wright, “WISE design of indoor wireless system: practical computation and optimization”, IEEE Computational Science and Engineering. pp. 58–68, Spring 1995.  K. Mehlhorn and S. N¨ her, “LEDA, A platform for combinatorial and a geometric computing”, Communications of the ACM. vol. 38, no. 1, pp. 196–106, 1995.  G. L. Turin, F. D. Clapp, T. L. Johnston, S. B. Fine, and D. Lavry, “A statistical model for urban multipath propagation”, IEEE Trans. Veh. Technol. vol. 21, no. 1, pp. 1–9, February 1972.  R.G. Kouyoumjian and P. H. Pathak, “A uniform geometric theory of diffraction for an edge on a perfectly conducting surface”, Proc. IEEE. vol. 62, no. 11, pp. 1448–1461, November 1974.  R. Luebbers, “Finite conductivity uniform GTD versus knife edge diffraction in prediction of propagation path loss”, IEEE Trans. Anten. Propag. vol. 32, no. 1, January 1984.  H. R. Anderson and P. Eng, “A ray-tracing propagation model for digital broadcast systems in urban areas”, IEEE Trans. on Broadcasting. vol. 39. no. 3, September 1993.