1. Introduction

The finite-difference time-domain (FDTD) [1] method for electromagnetic simulation is today one of the most efficient computational approximations to the Maxwell equations. Its accuracy has motivated several attempts to apply it to the prediction of radio coverage [2, 3], though one of the main limitations is still the fact that FDTD needs the implementation of a highly time-consuming algorithm. Furthermore, the deployment of metropolitan wireless networks in the last years has recently triggered the need for radio network planning tools that aid operators to design and optimize their wireless infrastructure. These tools rely on accurate descriptions of the underlying physical channel in order to perform trustworthy link- and system-level simulations with which to study the network performance. To increase the reliability of these tools, accurate radiowave propagation models are thus necessary.

Propagation models like ray tracing [4, 5] have been around already for some time. They have shown to be very accurate, as well as efficient from the computational point of view, except in environments like indoor where too many reflections need to be computed. In [6], a discrete model called Parflow has been proposed in the frequency domain, reducing a lot the complexity of the problem but bypassing the time-related information such as the delays of the different rays.

The FDTD model, which solves the Maxwell equations on a discrete spatial and temporal grid, can be also considered as a feasible alternative for this purpose. This method is attractive because all the propagation phenomena (reflections, diffractions, refractions, and transmission through different materials) are implicitly taken into account throughout its formulation. In [7], a hybridization of FDTD with a geometric model is proposed. In this approach, FDTD is applied only in small complex areas and combined with ray tracing for the more open space regions. Yet, the running time of such an approach is still too large to consider it for practical wireless networks planning and optimization. The evaluation of the FDTD equations at the frequencies of the current and future wireless networks (UMTS, WiMAX, etc.) requires the use of extremely small spatial steps compared to the size of the obstacles within the scenario. In femtocell environments such as residential areas, this would lead to the use of matrices that require extremely large memory spaces, making infeasible its computation on standard off-the-shelf computers. In order to solve this issue, a reformulation of the problem at a lower frequency [8] is possible and necessary.

The main contribution of this paper is thus the introduction of a heuristics-based calibration approach that solves the lower-frequency approximation by directly matching the FDTD prediction to real WiMAX femtocell measurements. The outcome of this calibration procedure will be the properties of the materials that best resemble the recorded propagation conditions. These can be later reused for further simulations in similar scenarios and at the same frequency. Nevertheless, propagation models always perform better if a measurements-based calibration is carried out in situ [9]. Hence, the approach presented here can also be implemented in a coverage prediction tool and be subject to calibration with new measurements for increased accuracy of the FDTD model in a given scenario.

Over the last few years, the traditional central processing units (CPUs) have started to face the physical limits of their achievable processing speed. This has lead to the design of new processor architectures such as multicore and the specialization of the different parts of computers. On the other hand, programmable graphics hardware has shown an increase on its parallel computing capability of several orders of magnitude, leading to novel solutions to compute electromagnetics [10]. Graphics chipsets are becoming cheaper and more powerful, being their architecture well suited for the implementation of parallel algorithms. In [11], for instance, a ray-tracing GPU implementation has been proposed. FDTD is an iterative and parallel algorithm, being all the pixels updated simultaneously at each time iteration. This fact makes FDTD an extremely suitable method to be implemented on a parallel architecture [12]. By following the recently released compute unified device architecture (CUDA) [13], this paper presents an efficient GPU implementation of an FDTD model able to reduce further the computing time.

One final problem to address when dealing with FDTD is the proper configuration of the absorbing boundary condition (ABC). For efficiency reasons, the convolutional perfectly matched layer (CPML) is to be used. In order to provide the highest absorption coefficient for the problem of interest, adequate parameters must be chosen so a method for the calibration of the CPML parameters is presented.

2. WiMAX Femtocells

Due to the flexibility of its MAC and PHY layers and to the capability of supporting high data rates and quality of service (QoS) [14], wireless interoperability for microwave access (WiMAX) is considered one of the most suitable technologies for the future deployment of cellular networks.

On the other hand, femtocell access points (FAPs) are pointed out as the emerging solution, not only to solve indoor coverage problems, but also to reduce network cost and improve network capacity [15].

Femtocells are low-power base stations designed for indoor usage that have the objective of allowing cellular network service providers to extend indoor coverage where it is limited or unavailable. Femtocells provide radio coverage of a certain cellular network standard (GSM, UMTS, WiMAX, LTE, etc.) and they are connected to the service provider via a broadband connection, for example, digital subscriber line (DSL). These devices can also offer other advantages such as new applications or high indoor data rates, and thus reduced indoor call costs and savings of phone battery.

According to recent surveys [16], around 90% of the data services and 60% of the mobile phone calls take place in indoor environments. Scenarios such as homes or offices are the favorite locations of the users, and these areas will support most of the traffic in the following years. WiMAX femtocells appear thus as a good solution to improve indoor coverage and support higher data rates and QoS. Furthermore, there are already several companies involved in the manufacture [17] and deployment [18] of these OFDMA-based devices.

Since a massive deployment of femtocells is expected to occur as soon as of 2010, the impact of adding a new femtocell layer to the existing macrocell layer stills needs to be investigated. The number and position of the femtocells will be unknown, and hence a controlled deployment of macrocells throughout traditional network planning can no longer be a solution used by the operator to enhance the network performance. Therefore, a detailed analysis of the interference between both layers, femto and macro, and the development of self-configuring and self-healing algorithms and techniques for femtocells are needed. Due to this, accurate network link-level and system-level simulations will play an important role to study these scenarios before femtocells are widely deployed.

Since femto-macrocell deployments will take place in hybrid indoor/outdoor scenarios, propagation models able to perform well in both environments are required. On the one hand, empirical methods [19] such as Xia-Bertoni or COST231 Walfish-Ikegami are not suitable for this task because they are based on macrocell measurements and are specifically designed for outdoor environments. Ray tracing has shown excellent performance in outdoor scenarios but its computational requirements become too large [20] when they come to compute diffraction- and reflections-intense scenarios. For instance, in indoor environments this results in long computation times [21], forcing ray-based approaches to restrict the amount of reflections that are computed. The same happens in cases where the simulation of street canyons requires a large number of reflections. On the other hand, finite-difference methods such as FDTD are able of accounting for all of the field interactions as long as the simulation is run until the steady state and the grid resolution describes accurately the environment. Therefore, these methods appear as an appealing and accurate alternative [22] for the modeling of hybrid indoor/outdoor scenarios.

3. Optimal FDTD Implementation

Since femtocells are designed to be located indoors and have an effect only in the equipment premises and a small surrounding area, in the case of low-buildings residential areas, properly tuned bidimensional propagation models should be able to precisely predict the channel behavior. The problem under consideration (femtocells coverage prediction) can be thus restricted to the two-dimensional case. Considering typical femtocells antennas with a vertical polarization and following the terminology given in [23], the FDTD equations can be written in the mode as follows:where is the magnetic field and is the electrical field in a discrete grid sampled with a spatial step of . , , and are the update coefficients that depend on the properties of the different materials inside the environment. , , , and are discrete variables with nonzero values only in some CPML regions and are necessary to implement the absorbing boundary.

However, the propagation of cylindrical waves in 2D FDTD simulations is by nature different from the 3D case. In order to minimize the error caused by this approximation, the current model is calibrated using femtocell measurements recorded in a real environment (see Section 5). This guarantees that the final simulation result resembles the real propagation conditions as faithfully as possible. It is also to be noticed that femtocell antennas are omnidirectional in the horizontal plane, emitting thus much less energy in the vertical direction. Moreover, in residential environments containing houses with a maximum of two floors, the main propagation phenomena occur in the horizontal plane. That is why restricting the prediction to the 2D case is only acceptable for this or similar cases, and not appropriate for constructions with bigger open spaces such as airports, train stations, or shopping centers.

From the computational point of view, restricting the problem to the 2D case is still not enough to achieve timely results for the study of femtocells deployments and their influence into the macrocell network. FDTD is very computationally demanding and therefore a specific implementation must be developed. The main purpose of this section is thus to present two techniques that aid to solve the scenario within reasonable execution times. The first technique reduces the complexity of the problem by increasing the spatial step used to sample the scenario, that is, it chooses a simulation frequency lower than that of the real system. The second technique presents a programming model that optimizes memory access for implementations in standard graphics cards.

3.1. Lower-Frequency Approximation and Model Calibration

The running time of the FDTD method depends, among other things, on the number of time iterations required to reach the steady state, that is, the stable state of the coverage simulation. To summarize, this number of iterations depends on the following.

(i)The number of obstacles inside the environment under consideration: the more the walls are, the more reflective and diffractive effects that will occur.(ii)The size of the environment in FDTD cells: a larger environment will need more iterations for the signal to reach all the cells of the scenario. In order to accurately describe the environment, the number of obstacles should not be reduced. It is thus interesting to try to reduce the size of the problem, which can be achieved by using a larger spatial step . To describe the simulation scenario, must also be small compared to the size of the obstacles. Furthermore, to avoid dispersion of the numerical waves within the Yee lattice, the spatial step also needs to be several times smaller than the smallest wavelength to be simulated [24]. For example, an WiMAX simulation would require a spatial step smaller than according to

Numerical dispersion in 2D FDTD simulations causes anisotropy of the propagation in the spatial grid. However, these effects can be reduced if a fine enough spatial grid is used. It is shown in [25] that with , the velocity-anisotropy error is , introducing thus a distortion of about cells for every propagated cells. However, these errors become meaningless after the calibration procedure introduced in Section 5.3, which corrects the power distribution so that it resembles the real propagation case according to the recorded measurements.

A scenario for the study of femto-to-macro interference has a typical size of around 100 × 100 meters so sampling the scenario with is not feasible in terms of computer implementation. A frequency reduction is thus necessary [26] to cope with memory and computational restrictions. This frequency reduction comes obviously at a cost because the reflections, refractions, transmissions, and diffractions behave differently depending on the frequency. Since the physical properties of the different materials are frequency dependent, reflections, refractions and transmissions through materials will vary. To overcome this problem, the approach presented here consists on performing a calibration of such parameters. This calibration, based on real measurements, will find values for the materials parameters in order to model, at a lower frequency, their behavior at the real frequency. This search is performed by minimizing the root mean square error (RMSE) between simulation and measurements, and the details of such a method are described in Section 5.3.

The effects of simulating with a lower frequency for WiMAX at have been already studied in [8], where it was shown that even after calibration, the predictions are still subject to an error due to diffractive effects. Nevertheless, it is well known that reflections dominate over diffractions in indoor environments, and the main power leakage of the femtocell from indoor to outdoor occurs by means of transmissions through wooden doors and glass windows (see Figure 1). Furthermore, in streets like the one shown in the current scenario, canyon effects caused by reflections are the main propagation phenomenon so it is clear that diffraction is not a significant propagation means in femtocell environments.

Figure 1: Example of a calibrated femtocell coverage prediction subject to diffraction errors due to lower-frequency FDTD simulation.

Additionally, it was shown in [8] that the absolute value of the error due to diffraction is limited and that the overall error of the simulation will depend only on the number of diffractive obstacles. In Section 5.4 a postprocessing filter is proposed as a means to reduce the fading errors due to this phenomenon. For comparative purposes, an unfiltered lower-frequency prediction is shown in Figure 1. The more accurate postprocessed prediction is explained later and displayed in Figure 9.

3.2. Parallel Implementation on GPU

If the previously described simplification reduces the size of the environment to simulate, the focus of this section is to present an implementation of the algorithm that reduces further the simulation time. In wireless networks planning and optimization, the aim is to run several system-level simulations to test hundreds of combinations of parameters for each base station. This implies that several base stations (emitting sources) must be simulated. It is thus necessary to reach simulation times on the order of seconds for each source. In order to reach this objective and since each cell of an FDTD environment performs similar computation (update of the cell own field values taking into account the neighboring cells), an approach is the use of parallel multithreaded computing.

The implementation of finite-difference algorithms on parallel architectures such as field-programmable gate-arrays (FPGAs) [27] and graphics processing units [28] has been recently highly regarded by the FDTD community. For instance, speeds of up to 75 Mcps (mega cells per second) have been claimed [29] for a 2D implementation in an FPGA. However, FPGAs are costly devices whose use is not as common as that of GPUs, which are present today in almost every personal computer. Therefore, the interest on programmable graphics hardware has increased and some solutions are already being proposed [10] as a feasible means of achieving shorter computation times for this type of algorithms.

By programming an NVIDIA GPU device with the new CUDA architecture [13], a 2D-FDTD algorithm has been implemented. With this technology, it is not necessary to be familiarized with the graphics pipeline and only some parallel programming and C language knowledge are necessary to exploit the properties of the GPU. This reduces the learning curve for scientists interested in quickly testing their parallel algorithms, while eliminating the redundancy of general purpose computing on GPU (GPGPU) code based on graphics libraries such as OpenGL.

The number of single instruction, multiple thread (SIMT) multiprocessors in each GPU varies between different cards, and each multiprocessor is able to execute a block of parallel threads by dividing them into groups named warps. Depending on the features (memory and processing capability) of a given multiprocessor, a certain number of threads will be executed parallely. It is thus important to balance the amount of memory that a thread will use, otherwise the memory could be fully occupied by less threads than the maximum allowed by the multiprocessor. It is in the programmer best interest to maximize the number of threads to be executed simultaneously [30]. Therefore and to maximize the multiprocessor occupancy, five different types of kernels (GPU programs) have been designed to compute different parts of the scenario as shown in Figure 2. The central area is the computational domain containing the scenario that needs to be simulated. Meanwhile, the other four areas represent the four absorbing boundary regions at the limits of the environment.

Figure 2: Fragmentation of the simulation scenario for independent kernels execution.

To compare the performance of such an implementation with traditional nonparallel approaches, the simulation of a pixels scenario has been tested under three different platforms. iterations were required to reach the steady state in this environment. MATLAB, which makes use of the AMD core math library (ACML) and is thus very optimized for matrices computation, is used as the nonparallel reference. Then a standard laptop graphics card (GeForce 8600M GT) and a high-performance computing card (TESLA C870) are tested. The main differences between these two cards are the number of multiprocessors (4 and 32) and the card memory (256 MB and 1.5 GB). The different performance results can be checked in Table 1.

Table 1: Performance of the algorithm running on different platforms when computing three thousand iterations of a scenario of size .

The achieved running time (8 seconds) for a complete radio coverage can be considered as a reasonably quick propagation prediction, fulfilling thus the requirements in terms of speed for wireless network planning in the presence of randomly distributed femtocells. This way, a high number of network configurations can be tested within acceptable times by the operator.

4. Calibration of the Absorbing Boundary Condition

FDTD is a precise method for performing field predictions in small environments and it has been widely applied in several areas of the industry, such as the simulation of microwave circuits or antennas design. But during many years, the computation of precise solutions in unbounded scenarios remained a great challenge.

In 1994 Berenger introduced the perfectly matched layer (PML) [31], an efficient numerical absorbing material matched to waves of whatever angle of incidence. The next improvement of this method occurred in 2000, when Roden and Gedney presented a more efficient implementation called convolutional perfectly matched layer (CPML) [32], which has since been one of the better regarded choices for this purpose.

However, the CPML must be carefully configured in order to properly exploit its full potential. The absorptive properties of the CPML depend mainly on the wave -vector, which is a function of the type of source being used, and it will therefore present different reflection coefficients for simulations performed at different frequencies. A proper selection of parameters is thus necessary.

An error function based on the reflection error of the CPML is presented next, as well as a continuous optimization approach to find its minimum in the solutions space formed by the CPML parameters.

4.1. The CPML Error Function

4.1.1. The Optimization Parameters

The CPML method maps the Maxwell equations into a complex stretched-coordinate space by making use of the complex frequency-shifted (CFS) tensorwhere, following the notation of [24], indicates the direction of the tensor coefficient.

In order to avoid reflections between the computational domain (CD) and the CPML boundary due to the discontinuity of , the losses of the CPML must be zero at the CD interface. These losses are then gradually increased [31] in an orthogonal direction from the CD interface to the outer perfect electric conductor (PEC) boundary. A polynomial grading of , , and has shown [24] to be quite efficient for this task:where is the depth of the CPML, and are the grading orders. An approximate optimal can also be estimated to outcome a given reflection error withwhere is the impedance of the background material [24].

However, which precise values of , , and to choose for a specific FDTD simulation remains an open question. The solution to this problem is thus the combination of parameters that configures the most absorbing CPML for a given source and number of FDTD time steps. Since the optimal value of is close to (5), the factor can be defined for notation simplicity and be subject to the optimization process. The intervals to search for the optimal solution when using a continuous soft source are presented in Table 2 and can be defined as

Table 2: Typical properties of the search parameters.

4.1.2. The Error Function

This section presents CPML calibration results for 2D simulations where the electrical field is the output magnitude from each FDTD simulation. In order to evaluate a given solution we compare it to a reference simulation that is free of reflections at the border. This reference simulation must be computed [24] using a grid large enough to avoid that reflections bounce back into the computational domain. As long as the FDTD simulation is implemented with first-order derivatives, a wavefront can only advance one cell per time step. In order to construct the extended grid in this case, the number of cells that must be added to the original grid in each direction can be thus calculated by simply considering the number of FDTD steps and the position of the source (see Figure 3).

Figure 3: Sounding points in a 2D grid of size . The depth of the extended grid in each direction varies depending on the position of the source.

To assess the optimal CPML configuration, it is necessary to analyze the time evolution of the simulated grid. For the sake of efficiency and to provide a reasonable estimation of the behavior of the CPML, the grid will be sounded only at certain key points. The highest reflection error occurs typically near the borders and corners of the CD so a homogeneous selection of sounding points is that shown in Figure 3.

The output of the reference simulation will therefore be the value of the electrical field at each sounding point with coordinates and at time step . Defining similarly the output of each optimization simulation as , the relative error for the same sounding point and at the same time step is

Each optimization simulation performs FDTD time steps. Therefore to obtain an indicator of the relative error value over the time, the RMS relative error is computed for each sounding point:

Finally, and in order to obtain a general indicator of the error for the whole scenario, the average value of (8) for all the sounding points is to be computed. The error function for a given combination of parameters can be thus defined as

Numerical experiments have shown that (9) does not vary much by adding more sounding points. represents therefore a good compromise between sounding efficiency and reliability of the error function.

4.2. The Calibration Process

4.2.1. The Optimization Algorithm

The objective of this section is to present a method to compute the combination of that minimizes (9). Several tests indicate that (9) is unimodal along the , , and dimensions, that is, (9) has only one minimum in the region given by (6). In order to find the optimum without evaluating the error function at a large number of candidate solutions, a smarter approach can be applied by minimizing (9) along each dimension sequentially and independently. Algorithm 1 presents this approach, being the stop condition the location of a satisfactory minimum lower than or the evaluation of a maximum number of iterations.

Algorithm 1: Minimization of the error function by means of coordinatewise minimization subroutines.

In order to find the minimum of the error function for each dimension of the space of solutions, it is necessary to evaluate (9) at several positions within the search intervals (6). Each of these evaluations needs to perform an FDTD simulation, which is the most time-consuming part of the algorithm. To minimize these, a Fibonacci search algorithm [33] is to be used. This algorithm narrows down the search interval by sequentially evaluating the error function at two positions within the interval and reusing one of these evaluations in the next step. Therefore only one function evaluation is necessary at each step. Table 2 contains the precision achieved for the example intervals and the required length of the Fibonacci sequence for each parameter.

4.3. ABC Calibration Results

Figure 4 presents a contour plot of the error function described by (9). The function values were obtained by computing the error at 2500 different locations of the 2D space of solutions given by for the optimal value of . The size of the FDTD scenario for this example is of 256 × 256 cells with the source located at the coordinates and being the spatial and time steps and respectively. The CPML has a depth of 16 cells and a total of FDTD time steps were performed to compute each value of the error function. The applied source was a Gaussian pulse with a temporal width of and modulated with a sinusoidal frequency of , which is the frequency of WiMAX in Europe.

Figure 4: Contour plot of the error function with for a modulated Gaussian pulse of width 0.4 nanosecond and an oscillating frequency of 3.5 GHz. The graph also shows the solutions found by Algorithm 1 and the evolution until the optimum.

Figure 4 also displays the error points found at each iteration of Algorithm 1 after minimizing in the and dimensions. In this example, is initialized with a random value within its range and the optimal solution is reached in just 3 iterations. Without fixing and optimizing in all three dimensions, the minimum is reached in only 4 iterations. But clearly the number of required FDTD simulations is much greater and can be calculated byTo obtain, for instance, the precision shown in Table 2, accounts for a total of simulations. Using the previously mentioned parallel computing architecture, these can be computed in less than 2 minutes on a laptop graphics card.

Once the algorithm has converged, the quality of the solution can be tested by computing an FDTD simulation using the obtained CPML calibration parameters. Figure 5 presents the change over time of the relative error at a corner point in the scenario described by Figure 3. It is clear in this example that the relative error never exceeds , yielding thus an excellent absorption coefficient.

Figure 5: Time evolution of the relative error (solid line) at the upper left point (see Figure 3). The dash-dotted line is the value of the electrical field at the same sounding point.

5. Calibration of the Propagation Model

In FDTD, the parameters that define each material and therefore play an active role in the final simulation result are three: (i)relative electrical permittivity ;(ii)relative magnetic permeability ;(iii)electrical conductivity .

Due to the 2D and lower-frequency simplifications applied to this model, it should not be expected that the values of the materials parameters at the real frequency perform the same as at the simulation frequency. The correct values of these parameters must be therefore chosen carefully in order for the simulation result to resemble faithfully the reality. As advanced in Section 3.1, this can be achieved by using real coverage measurements to find the proper combination of parameters that better match the prediction to the measurements.

5.1. Coverage Measurements

In order to measure the accuracy of the presented model, a measurements campaign has been performed. The chosen scenario was a residential area with two-floor houses in a medium-size British town. The femtocell excitation is an oscillatory source implemented on a vector signal generator and configured as shown in Table 3. The emitting antennae are omnidirectional in the azimuth plane with a gain of 11 dBi in the direction of maximum radiation.

Table 3: Main parameters of experimental femtocell.

Since one of the main objectives of this work is to introduce a propagation model for the study of interference scenarios in femtocells deployments, the measurements have been performed mainly outdoors. This way, the indoor-to-outdoor propagation case, proper of femto-to-macro interference scenarios, is characterized. Figure 6 shows the collected power data laid over a map of the scenario under study.

Figure 6: Power measurements and simulation scenario. The location of the transmitter is marked with a magenta square.

5.1.1. Measurements Postprocessing

Raw power measurements are not yet useful for the calibration of a finite-difference propagation model. The data must first undergo a postprocessing phase during which outliers will be removed. Such postprocessing is detailed next.

Removal of Location Outliers
The location of the outdoor measurement points has been obtained using GPS coordinates but these coordinates are sometimes subject to errors. At this stage every measurement matching the next properties must be removed: out of range GPS coordinates, coordinates inside of a building, no GPS coverage or coordinates outside of the scenario.

Removal of Noise Bins
In areas of low coverage, it is possible that the measured signal becomes indistinguishable from the background noise. Those measurements are thus also classified as outliers. In order to clearly distinguish signal from noise, a large recording of the noise in the examined frequency band and location area has been performed. This way, the noise has been found to follow an approximate normal distribution with mean of and a standard deviation of . Any measurement value that falls within a range of is thus considered to be an outlier.

Spatial Filtering
The used source is a narrowband frequency pulse. Therefore, the collected measurements are also subject to narrowband fading effects which are usually modelled using random processes. In order for these measurements to be useful for the calibration of deterministic models, the randomness due to fading needs to be removed. Hence, a spatial filtering of the measurements has been applied by following the 40-Lambda averaging criterion [34]. The final state of the measurements is shown in Figure 7.

Figure 7: Power measurements after postprocessing.

5.2. The Materials Error Function

The objective of the model tuning is to configure the materials involved in the FDTD simulation so that they show in the computational domain a similar behavior to the reality. If represents the properties of material , a solution to a problem involving materials is thus :

Each measurement point (with and the number of points) is assigned a measured power value . Similarly and for an FDTD prediction calibrated with the same point can be assigned a predicted power value . The error of the prediction at point can be then expressed asbeing the mean error of all points, which can also be interpreted as the offset between the measurements and the predictions. Once the model is calibrated, the tuned mean error is computed. Then the of any other prediction can be greatly reduced by simply adding to the predictions.

The root mean square error is often used as a good estimate of the accuracy of a propagation model. The RMSE will hence be the error function subject to minimization. For an FDTD configuration , the RMSE can be thus computed as

5.3. Metaheuristics-Based Calibration

Once the error function has been defined, a brute-force approach to find an optimal solution to the problem could be, for instance, to test all possible until a solution that minimizes (13) is found. Since the properties of the materials are all real, the space of solutions for is infinite and a smarter approach is needed. In this work, a meta-heuristics optimization algorithm is proposed as a feasible way of searching the space of solutions. The algorithm applied here is simulated annealing, though the same concept also applies to other heuristic algorithms, as long as they are properly configured.

Simulated Annealing (SA) [35] is an optimization algorithm based on the physical technique of annealing, widely used in metallurgy. From the point of view of the minimization of an error function, SA works by setting the state of the system to a solution , and evaluating neighbor solutions to try to find a better one (). If a better solution is found, then the current state of the system is updated to the new solution . If, however, a worst solution is found, the state of the system is set to this new neighbor solution with probability . is called the acceptance probability function (APF) and it is a function of , and a variable called the temperature that is progressively decreased as the calibration progresses. The acceptance probability function must meet certain requirements in order to accept better solutions than the current state and worst solutions when the temperature is high, that is, at the beginning of the calibration process. A simple APF that follows these criteria isbut the user of SA is free to choose any APF to its convenience.

The way the temperature is decreased is also subject to many implementations. In this paper, the value of the temperature at each stage is obtained as follows:with and is the number of different temperature levels. is called the annealing factor and it is related to the rate with which the temperature decreases from one stage to the next one.

The evolution of the state of the system by means of SA is displayed in Figure 8, as well as the evolution of the temperature. For this calibration, different levels of temperature have been defined and the system is let free to test different neighbors at each temperature level. This way, the physical process of annealing is resembled much more faithfully than if the temperature was decreased at each SA iteration. The idea behind this is to allow the system to perform a deeper search at each temperature level before decreasing its chances of escaping local minima.

Figure 8: Evolution of the RMSE of the FDTD prediction when choosing the materials parameters using simulated annealing. The temperature is expressed in natural units, and .

Figure 9: Filtered coverage prediction of a WiMAX femtocell with a 3.5 GHz measurements-based calibrated FDTD model.

The way neighbor solutions are chosen can also be decided freely by the user. Since the purpose here is to find the optimal values of different materials, only one material is modified at each stage. Furthermore, only one parameter of this material is modified. This way, a local search in the very neighborhood of the current state is guaranteed.

The calibration displayed in Figure 8 is performed using the measurements and scenario shown in Figure 6. For this scenario and according to the most commonly used construction materials in the United Kingdom, five different materials have been assumed: air as the background material, plaster for the inner walls, wood for the doors and furniture, glass for the windows, and brick for the houses outer walls. The final values of the parameters for these materials after the calibration are shown in Table 4. The electrical conductivity is expressed in and the refraction index , computed as , is provided as reference.

Table 4: Calibrated parameters of the materials at 3.5 GHz.

5.4. Fading Removal Filter

The spatial step for this calibration is with for good isotropic propagation, yielding thus a wavelength of . This means that the simulation frequency is approximately , while the real frequency of the WiMAX measurements is . Following the terminology presented in [8], the frequency reduction factor is defined aswhich has in this case a value of . Due to the reasons expressed in Section 3.1, a prediction performed with the final calibration results of Table 4 is still subject to errors at diffracting obstacles. Such an error is limited and can be easily evaluated for each obstacle withwhere is a geometrical parameter that depends on the specific disposition of the scenario (see [36] for details).

Since diffraction introduces wrong fading effects, a spatial (2D) average moving filter has been applied as a postprocessing method to reduce the impact of the frequency reduction. A decrease of up to has been observed in the value of the , and up to in macrocell-calibrated models. A coverage prediction performed by the calibrated FDTD model and postprocessing filter is shown in Figure 9 along with the measurements used for the calibration.

After the postprocessing filter, the final obtained is of and a comparison between the FDTD predictions and the measurements is displayed in Figure 10.

Figure 10: Comparison between the FDTD predictions and the measurements at 3.5 GHz. .

5.5. Accuracy Validation

Finally, in order to assess the accuracy of the FDTD propagation model, calibrations have been performed at the real and several lower frequencies. The analyzed range of simulated frequencies comprises values of between and , being displayed in Figure 11 the errors of the simulations after calibration. From this figure it is also clear how the filtering introduced in the previous section contributes to the reduction of the .

Figure 11: Evolution of the after calibration, with respect to the frequency reduction factor .

Furthermore, the data also shows that proper lower-frequency calibrations of the model are able to reach performances close to that of the true frequency. However, the simulation frequency should not be reduced indefinitely. This is because of the increase in the size of the spatial step as decreases. If becomes too large, the spatial resolution might not be enough to accurately describe the simulation scenario and the diffraction phenomena, bypassing some features of the environment. As a consequence of this, the error grows quickly and reaches values that could be achieved with simpler propagation models. Therefore, a compromise between the computational complexity and the model accuracy must be achieved. For the scenario under consideration, Figure 11 shows that a value of has been chosen. This value, located in the elbow of the curve guarantees a low error without compromising the execution time and is used in Section 6 to perform system-level simulations of WiMAX femtocells.

In order to examine the achievable accuracy in the overall scenario, a different measurements route has also been used to test the coverage prediction. For this purpose, the transmitter was placed in a different room within the femtocell premises and new measurements were recorded along the street. When compared to the FDTD prediction performed with the total error was which differs from the originally calibrated error. This indicates that the accuracy of the model calibration can be improved by taking more points into consideration. Nevertheless it also indicates that the results presented in Table 4 can still be used in similar scenarios while keeping reasonable levels.

The reduction of the simulation frequency also has an effect on the interference patterns that arise in the simulation as a result of phase differences in the propagated waves. In order to analyze the phase behavior of the simulation, the received power distribution is illustrated in Figure 12 as a box plot. The lower and upper limits of the boxes represent the first and third quartiles, while the red horizontal line is the median of the data. The mean received power is indicated by a dark dot and the extremes of the whiskers are located one standard deviation below and above the mean. Due to the calibration of the received power, it is clear from this figure that the overall power distribution remains approximately invariant for those values of where a low-prediction error can be achieved.

Figure 12: Dependence of the power distribution with respect to the frequency reduction factor .

6. System-Level Simulations (SLSS)

The applicability of the presented propagation predictions to the study of a macro-femtocell hybrid scenario is presented here by means of mobile WiMAX (IEEE 802.16e-2005) system-level simulations with private access femtocells. The target of this experimental evaluation is to show how a measurements-based calibrated FDTD model can help the operator to predict common interference problems between users of the macrocell and the femtocell.

The scenario used for this experimental evaluation was the same residential street presented in Figure 6. A nonuniformly deployed WiMAX hybrid network formed by one macrocell and five femtocells was used for this case of study. The femtocells were located in five different households along the street, while the macrocell is positioned in an area located further away from the street under consideration. This is realistic, since femtocells are mainly aimed at users with poor indoor macrocell coverage. To perform the system-level simulation, different traffic maps were used for both the indoor and outdoor environments. There is one indoor traffic map per femtocell and house, which contains two randomly positioned users, and there is one outdoor traffic map in the street, containing five users.

The static system-level simulator functions by recording hundreds of snapshots with random positions of the macrocell and femtocell users. As the power distribution remains constant with the reduction of the (see Figure 12) and the location of the users varies between different snapshots, particular phase errors at given sites in the coverage prediction do not affect the final SLS statistics. Furthermore, it has been experimentally confirmed that the probability of outage, as well as the average throughput of the different cells in the system-level simulations, is not altered by the reduction of the simulation frequency.

This case of study makes use of private access femtocells, which means that indoor users are allowed to connect, depending on the signal quality, to their own femtocell or to the macrocell. On the other hand, outdoor users are only allowed to connect to the macrocell. For illustration purposes of the applicability of the presented propagation model, only downlink is considered.

It is illustrated in Figure 13 that an outdoor user connected to a distant macrocell is jammed due to the interference coming from nearby femtocells. In this case, the green users are successful, while the blue users suffer outage in downlink. A user will be successful or in outage depending on whether they are able or not of obtaining their requested throughputs and QoS from the network in order to use their services (video). In the example shown here, it occurs that there are three users on the street connected to the macrocell, who are using the same WiMAX subchannel as another femtocell user during the same time interval (symbol). In this case and as predicted by the FDTD model, the signal level of the carrier is smaller than the signal power of the interference, resulting thus in a poor signal quality. Due to this, the macrocell user is jammed and the communication cannot be supported by the network.

Figure 13: WiMAX system-level simulation in a hybrid femtocell/macrocell scenario.

7. Conclusion

In this paper, the coverage prediction of WiMAX femtocells by means of a calibrated FDTD model is studied. The reduction of the simulation frequency is proposed as a simplification of the problem which is required for computational reasons. The use of a parallel architecture such as the computation on a graphics card is also proposed as a feasible mean of reducing the computation time.

An optimal method to obtain an acceptable combination of parameters, which maximizes the absorbing properties of the CPML boundary condition for FDTD electromagnetic simulations, is also proposed. Furthermore, an error function that measures the relative error of the electrical field prediction near the CPML has been modelled. In addition to this, a Fibonacci search-based method is presented as a fast way to explore the solutions space and reach the minimum point without falling in the need to compute the error function at thousands of different solutions.

A method for the calibration of the materials involved in the FDTD simulation has also been presented. This model tuning approach, based on simulated annealing, is introduced as a way to match the propagation predictions to the reality. Then, a spatial averaging filter has been used as a mean to solve prediction errors at diffractive obstacles due to the lower-frequency simplification. The accuracy of the method has been validated by performing calibrations at a wide range of simulation frequencies, analyzing the power distribution and evaluating the predictions with a different measurements route.

Finally, system-level mobile WiMAX simulations that use this FDTD propagation model have been presented. This exemplifies the interference caused by indoors-located WiMAX femtocells to outdoor users of the macrocell. This way, the need for hybrid indoor/outdoor propagation models is evinced.

Acknowledgments

This work is supported by the EPSRC-funded research Project EP/F067364/1 with title “The feasibility study of WiMAX based femtocell for indoor coverage.” It is also partially supported by two EU FP6 projects on 3G/4G Wireless Network Design: “RANPLAN-HEC” with Grant no. MEST-CT-2005-020958 and EU FP6 “GAWIND” with Grant no. MTKD-CT-2006-042783.