#### Abstract

This paper focuses on the modeling, simulation, and experimental verification of wideband single-input single-output (SISO) mobile fading channels for indoor propagation environments. The indoor reference channel model is derived from a geometrical rectangle scattering model, which consists of an infinite number of scatterers. It is assumed that the scatterers are exponentially distributed over the two-dimensional (2D) horizontal plane of a rectangular room. Analytical expressions are derived for the probability density function (PDF) of the angle of arrival (AOA), the PDF of the propagation path length, the power delay profile (PDP), and the frequency correlation function (FCF). An efficient sum-of-cisoids (SOC) channel simulator is derived from the nonrealizable reference model by employing the SOC principle. It is shown that the SOC channel simulator approximates closely the reference model with respect to the FCF. The SOC channel simulator enables the performance evaluation of wideband indoor wireless communication systems with reduced realization expenditure. Moreover, the rationality and usefulness of the derived indoor channel model is confirmed by various measurements at 2.4, 5, and 60 GHz.

#### 1. Introduction

Indoor wireless communications has attracted considerable attention in recent years due to a broad set of emerging indoor services offered by personal communication systems [1, 2], wireless local area networks [3, 4], and wireless private branch exchange networks [5]. The interest in indoor wireless communications is growing especially when optical wireless communication techniques [6, 7] are applied to indoor deployment scenarios. Numerous researches on modulation techniques [8], multiple access techniques [9], and propagation modeling [10] have been carried out for indoor optical wireless communications. In order to design an efficient indoor wireless communication system and to predict accurately its system performance, a prerequisite is to develop an appropriate indoor channel model that describes realistically the underlying propagation conditions.

To improve the indoor channel characterization and modeling, numerous measurement campaigns have been conducted in a variety of indoor environments, such as offices, corridors, buildings, and factory halls. Measurement results have been reported in the literature for various frequency bands, including the frequency bands at 900 MHz [11–13], 1.5 GHz [14, 15], 2.4 GHz [16], 4–5.5 GHz [17–19], 17-18 GHz [18, 20], and 60 GHz [21–24]. Based on measurement results, several empirical statistical channel models [13, 24–27] have been developed for the simulation and performance analysis of various indoor communication systems. The advantage of measurement-based channel models is that they are characterizing the fading behavior utterly realistic. By changing the model parameters, the developed statistical channel models can be employed to simulate indoor propagation channels for a broad range of different scenarios. However, to determine appropriate parameters for different propagation scenarios, a large number of measurement data are required, which is time- and material-consuming and leads thus to high development costs. Alternatively, approaches based on ray-tracing techniques [28–31] have been developed to enable the simulation of indoor propagation channels. Ray-tracing channel models [29, 32–34] can efficiently capture the fading behavior of specified indoor environments. Many developed ray-tracing channel models have shown good agreement with measured channels [35–37]. However, ray-tracing models greatly depend on the physical layouts and materials through which electromagnetic waves propagate, for example, walls, floors, and ceilings. In addition, the main drawback of ray-tracing models comes from their computational costs, which are widely determined by the size and complexity of the geographic database as well as the strategy used for the ray search. Therefore, a trade-off has to be found between the prediction accuracy and the simulation efficiency when modeling indoor propagation channels using ray-tracing techniques.

The two main types of indoor channel models, namely, the empirical statistical channel models and the ray-tracing channel models, have their own strengths and limitations. To cope with the drawbacks mentioned above, a geometrical channel model has been proposed in [38] for indoor environments, where it is assumed that the scatterers are randomly distributed inside of a circle with a base station (BS) in its centre. More specifically, this model is based on the assumption that the distance between the BS and the scatterers follows an exponential distribution. Such an assumption imposes severe restrictions on the applicability of the model. In addition, from the physical point of view, it is questionable whether it is reasonable to represent the geometry of indoor scattering areas, such as offices and corridors, by a circle. On the contrary, a rectangle is a more appropriate geometric figure to model typical indoor propagation environments. A geometrical rectangle scattering model has been proposed in [39] to characterize narrowband indoor radio propagation channels. However, to the best of the authors’ knowledge, the modeling of wideband indoor radio channels by employing an appropriate geometrical scattering model that takes into consideration the unique indoor propagation conditions is still an open problem.

Motivated by the scarcity of proper wideband geometry-based indoor channel models, such a model was developed in [40] by extending the geometrical rectangle scattering model in [39] with respect to frequency selectivity. In this paper, we improve the statistical models in [39–41] for the characterization of the scatterer distribution. We concentrate on the statistical characterization of a wideband reference channel model, which is based on the assumption that an infinite number of scatterers are randomly distributed over the two-dimensional (2D) horizontal plane of a rectangular room. To be more specific, we assume that the probability for the occurrence of a scatterer decays exponentially with the distance from the walls. This assumption includes the uniform distribution of the scatterers as a specific case. In contrast to [42–44], we apply a different procedure for studying the statistical properties of the resulting reference channel model, which requires no knowledge of the complex channel gain. We derive the analytical expressions for the probability density function (PDF) of the angle of arrival (AOA), the PDF of the propagation path length, the power delay profile (PDP), and the frequency correlation function (FCF). Moreover, we derive a wideband sum-of-cisoids (SOC) channel simulator from the reference model. It is shown that the designed wideband SOC channel simulator matches the underlying reference model accurately with respect to the FCF. The obtained SOC channel simulator enables the simulation of indoor mobile fading channels with reduced realization expenditure. At the end, we demonstrate the closeness to reality of the proposed reference channel model by comparing its mean excess delay and the root mean square (RMS) delay spread with the corresponding empirical quantities obtained from measured channels in laboratories at 2.4 GHz [16] and a conference room at 5 GHz [19] as well as corridors and offices at 60 GHz [23].

The remainder of this paper is structured as follows. In Section 2, we introduce the 2D geometrical indoor scattering model, which serves as a starting point for the derivation of the reference channel model. Section 3 analyzes the statistical characteristics of the wideband reference channel model with emphasis on the PDF of the AOA, the PDF of the propagation path length, the PDP, and the FCF. Section 4 introduces a design procedure for wideband indoor SOC channel simulators. Numerical and experimental results are presented in Section 5 to confirm the correctness of the obtained theoretical results. Section 6 validates the usefulness of the proposed indoor channel model by matching its statistical properties to those of measured (real-world) channels. Finally, Section 7 draws the conclusion.

#### 2. The Geometrical Indoor Scattering Model

This section briefly describes the geometrical scattering model, which was first proposed in [40] to characterize indoor propagation scenarios. Although a room has three dimensions, we only consider the 2D horizontal plane in which all local scatterers as well as the BS and the mobile station (MS) are located. As shown in Figure 1, the rectangle represents the 2D horizontal plane of a room. Its length and width are denoted by and , respectively. The BS is considered as the transmitter, leaving the MS to play the role as the receiver. For convenience, it is also assumed that the MS moves along the -axis. Moreover, we consider single bounce scattering, which means that the plane waves emitted from the BS are only bounced once by scatterers before reaching the MS. The black star in Figure 1 represents a single local scatterer. In real-world environments, the number of scatterers is limited and their locations differ from one propagation scenario to another. To avoid studying the fading characteristics for a specific indoor propagation scenario characterized by a specific realization of a finite number of scatterers, we focus on a general statistical model that results from averaging over all possible propagation scenarios. Therefore, we assume that an infinite number of scatterers are randomly distributed inside the room. Such a model acts here as a nonrealizable stochastic reference model. From the reference channel model, an efficient channel simulator with low realization expenditure can be derived by making use of the SOC principle [45].

A straightforward assumption would be that the scatterers are uniformly distributed inside the room. For the modeling of indoor channels, however, it is more realistic to assume that the probability of the occurrence of a scatterer decays exponentially with the distance from the walls. This implies that it is more likely that the received scattered components are coming from the walls rather than from objects located in the room’s centre. Instead of the uniform distribution, we use therefore a mixture of exponential functions to model the distribution of the scatterers inside the 2D horizontal plane of the room. The exponential functions are more general and include the uniform distribution as a special case, as explained in the next section.

#### 3. Statistical Characterization of the Reference Channel Model

This section studies the statistical properties of the proposed indoor wideband reference channel model. In this regard, the PDF of the AOA, the PDP, and the FCF are of special interest.

It is shown in Figure 1 that the BS is located at the position , while the MS is located at the origin of the coordinate system. The position of the local scatterer is described by , where and are independent random variables. As mentioned in Section 2, it is assumed that the scatterers are exponentially distributed within the 2D horizontal plane of the room. This implies that the random variables and must be exponentially distributed over the intervals and , respectively. Thus, the PDF of and the PDF of are given bywhere are real-valued parameters andIt should be mentioned that the distributions in (1a) and (1b) include the uniform distribution as a special case if for . With the rule of de l’Hospital, it can be shown that in this case and tend to and , respectively. Owing to the assumption that the random variables and are statistically independent, we can express the joint PDF of and as the product of the marginal PDFs and ; that is,

Next, we describe the position of the scatterers by using polar coordinates , where denotes the path length from the position of the scatterer to the MS, and is the AOA (see Figure 1). Notice that and are random variables. The joint PDF of the position of the scatterers in polar coordinates can be written by using (3) asLet be the propagation path length, which is the plane wave’s overall traveling distance from the BS via the scatterer to the MS; then is given byOur next focus is on finding the joint PDF of the path length and the AOA using (4). For this purpose, we introduce the auxiliary random variable . Solving the system of equations determined by and gives us the following two solutions:By applying the fundamental theorem of transformation of random variables [46, p. 201], the joint PDF of the path length and the auxiliary random variable is given bywheredenotes the Jacobian determinant of the inverse transformation. Substituting (8) into (7) and using (4), (6a), and (6b), we find the expression for the joint PDF . Owing to , the joint PDF can be expressed directly using the joint PDF by replacing in (7) by . Thus, we obtainwhere and are the densities given by (1a) and (1b), respectively.

##### 3.1. The PDF of the AOA

In this subsection, we are concerned with the derivation of the PDF of the AOA . This PDF can be obtained from the joint PDF in (9) by integrating over the range of ; that is,where represents the interval in which and are the minimum and the maximum propagation path length, respectively. The minimum propagation path length follows from (5) for ; that is, . Notice that the maximum propagation path length depends on . It is shown in Appendix A that is given as in (11), where

With and , as given in (11), the PDF can readily be obtained by solving the integral in (10) numerically.

##### 3.2. The PDF of the Propagation Path Length

The PDF of the propagation path length , denoted by , can be calculated by means of the relationIt is shown in Appendix B [see (B.9)] that the integral in (13) can be expressed by (14), which is presented as follows:where the path lengths are given in Appendix B, and denotes the integralin which

##### 3.3. The PDF and PDP of the Propagation Delays

The objective of this subsection is to find the PDF of the propagation delays . Here, the propagation delay is defined aswhere denotes the speed of light. Taking this relation into account and applying the concept of transformation of random variables [46, p. 130], we can express the PDF of the propagation delays asHere, can be obtained directly from the PDF in (14) by replacing the variable by the term . For brevity, the expression for is not presented here. Let the total power of the received multipath components be represented by , and let us denote the PDP by . Then, it follows that holds. By taking the relationship and the property into account, we can conclude that the following equation must hold for the PDP:Hence, the PDP can be obtained directly from the PDF in (14).

##### 3.4. The FCF

According to the Wiener-Khinchin theorem, the PDP and the FCF form a Fourier transform pair. The FCF is defined as the Fourier transform of the PDP ; that is,Substituting (19) in (20) results in the following expression for the FCF:Since no closed-form solution exists, the integral above has to be solved numerically.

#### 4. Design of an SOC Wideband Indoor Channel Simulator

This section deals with the design of a stochastic SOC channel simulator for the proposed wideband indoor channel model. The time-variant impulse response of a wideband SOC channel simulator consisting of discrete propagation paths can be expressed as follows [47, Eq. ]:Here, the path gains are determined by the square root of the PDP assigned to the th discrete propagation delay. By applying the procedure described in [47, p. 374], we havewhere are the intervals defined in [47, p. 374]. The stochastic complex process in (22) models the sum of all received scattered components having the same propagation delay . According to the SOC principle, such a stochastic process can be represented by a sum of cisoids of the form [45]Here, the quantities , , and represent the gain, the Doppler frequency, and the Doppler phase of the th component of the th discrete path, respectively. The gains and the Doppler frequencies are constant. They can be determined by a proper parameter computation method, for example, the modified method of equal areas (MMEA) [48], according to which the gains are given by

The Doppler frequencies are determined by the AOAs according to the well-known relationwhere stands for the maximum Doppler frequency. The AOAs are constant parameters, which can be computed by means of the MMEA. This method requires solving the following equation:for all . The phases of the stochastic SOC channel simulator are independent and identically distributed (i.i.d.) random variables, each following a uniform distribution over . Thus, the stochastic SOC channel simulator can be interpreted as a family of sample functions determined by the realizations of the phases . A sample function can be obtained from the stochastic channel simulator by fixing all phases. Such a sample function, which can be interpreted as a waveform generated by a deterministic channel simulator, can be used to simulate the proposed wideband indoor channel.

#### 5. Numerical and Simulation Results

This section illustrates the theoretical results given by (10) and (21). The correctness of the theoretical results will be verified by simulations. We will also show that the SOC channel simulator matches the reference channel model with respect to the FCF. The validity of the proposed indoor channel model is also confirmed by indoor channel measurements by studying the mean delay and the RMS delay spread. In all simulations, we consider a room with length m and width m as our indoor environment. The parameters and have been chosen to be equal to 2 and 1, respectively. The SOC channel simulator is designed by applying the MMEA [48] using cisoids for each discrete propagation path.

The theoretical result of the PDF of the AOA [see (10)] of the wideband reference channel model is presented in Figure 2 for different BS locations ( and ). It can be observed from this figure that the shape of the PDF of the AOA is independent of the position of the BS. The theoretical results illustrated in Figure 2 have been verified by experiments. In our experiments, we determined the horizontal (and vertical) locations of all scatterers in the Cartesian coordinate system as outcomes of a random generator, which generates uniformly distributed scatterers over and according to the distributions in (1a) and (1b), respectively. It should be mentioned that the results in Figure 2 correspond to the PDF of the AOA of the narrowband channel model (see [39, Eq. ]). In Section 6, we will graphically present the PDF of the AOA in case that the locations of the scatterers are exponentially distributed. The impact of the parameters and on the shape of the PDF of the AOA is illustrated in Figures 3 and 4, respectively. It can be observed from Figure 3 that an increase in the value of from to results in a decrease in the PDF at small values of the AOAs, where they are confined to the interval . In contrast, the PDF increases with increasing values of at large values of the AOA; that is, .

Figure 4 reveals that an increase in leads to an inverse effect on the shape of the PDF compared to an increase in . Figures 3 and 4 show that both the parameters and , which control the symmetry of the room, have impacts on the shape of the PDF of AOA. The main reason is that when or changes, it means the location of the mobile station relative to the room area moves (see the geometrical indoor scattering model in Figure 1). While a scatterer is fixed, the position of the mobile station determines the value of the AOA. If the mobile station moves closer to a wall, the AOA bounced by the scatterers in the triangle formed by the mobile station and two corners of the wall takes a larger range. Therefore, the PDF of one angle from that range decreases.

Figure 5 shows the absolute value of the FCF that has been calculated by using (21) for different values of the room length . For comparison purposes, we also plot the theoretical curve for the absolute value of the FCF of the designed SOC simulation model by using Eq. in [42]. We consider discrete propagation paths for the SOC channel simulator, where the corresponding delays are equally spaced between 0 and the maximum propagation delay . The value of can be obtained from (18) as the maximum value of such that . The path powers assigned to different propagation paths have been calculated according to the method described in [47, pp. 374–375]. The AOAs of the SOC channel simulator are computed by employing the MMEA [48]. We observe from Figure 5 that the FCF of the SOC channel simulator can be brought into extremely good agreement with that of the reference model. The FCF decays faster with increasing the frequency separation if the room length increases (here from 10 m to 30 m). This means that the coherence bandwidth becomes smaller as increases (the coherence bandwidth of a wideband mobile radio channel is the smallest positive value of the frequency separation variable for which the condition is fulfilled).

The influence of the room width on the FCF is presented in Figure 6 from which similar conclusions can be drawn as from Figure 5. When the dimensions of the room increase, the frequency correlation of time-variant transfer functions reduces at the same frequency separation. In other words, the time-variant transfer functions of the developed SOC channel simulator become independent at a smaller frequency separation, which leads to a smaller coherence bandwidth.

#### 6. Experimental Verification

To demonstrate the usefulness of the proposed geometrical-based indoor reference channel model, we will show how the statistics of the reference channel model can be fitted to the statistics of measured real-world channels by optimizing the key parameters of the reference model. Here, two of the most important characteristic quantities of wideband mobile radio channels are considered, namely, the mean excess delay and the RMS delay spread.

The mean excess delay is defined as the first moment of the PDP [26, 49]; that is,The RMS delay spread is defined as the square root of the second central moment of the PDP [26, 49]; that is,

The mean excess delay and the RMS delay spread of the reference model can easily be calculated by substituting the PDP [see (19)] into (28) and (29), respectively.

The measurement results of the mean excess delay and the RMS delay spread considered in this section are taken from [16, 19, 23]. The measurement campaigns were carried out in laboratories, a conference room, and corridors. The measurements in [16] have been obtained at 2.4 GHz in a laboratory with dimensions 7.8 m × 9.95 m. Furthermore, the measurements reported in [19] have been obtained at 5 GHz in a conference room with dimensions 6.6 m × 5.9 m. Finally, the measurements in [23] have been obtained at 60 GHz in two 30 m 1.75 m × 2.80 m corridors as well as inside a laboratory with dimensions 19.5 m × 7.5 m.

We combine all the relevant model parameters (except the parameters and ) which control the mean excess delay and the RMS delay spread of the reference channel model into a parameter vector denoted and defined by . We also introduce the following error function:for measuring the deviations between the mean excess delays and as well as between the RMS delay spreads and . In (30), and represent proper weighting factors.

The optimization of the parameter vector has been carried out by minimizing numerically the above error function by means of the quasi-Newton optimization procedure [50]. The optimization results are presented in Tables 1, 2, and 3, for the laboratory, the conference room, and corridor scenarios. For comparison purpose, the measured mean assess delay and the measured RMS delay spread , presented in [16, 19, 23], are also shown in Tables 1, 2, and 3. Note that the dimensions of the different rooms, which are controlled by the parameters and , have been chosen such that they are equal to the lengths and widths of the rooms in which the corresponding measurement campaigns have been conducted. Hence, and are fixed, and thus they have not been included in the optimization procedure. In Table 1, the measured mean assess delay is in the range from ns to ns, while the RMS delay spread varies from ns to ns. Similarly, in Table 2, the measured mean assess delay is in the range from ns to ns, while the RMS delay spread varies from ns to ns. In Table 3, finally, the measured mean assess delay covers the range from ns to ns and holds for hallways, while the range from ns to ns refers to offices. In addition, the measured RMS delay spread varies from ns to ns for corridors and from ns to ns for laboratories. In our optimization procedure, the weighting factors and were chosen to be and , respectively. For reasons of brevity, we have only presented the values of the optimized parameters in all tables with two decimals of precision. As shown in Tables 1, 2, and 3, the mean excess delay and the RMS delay spread of the proposed wideband indoor channel model are very close to the corresponding measured quantities for all considered indoor propagation scenarios, which demonstrates that the reference channel model is useful for characterizing real-world wideband indoor mobile radio channels.

For completeness, we present in Figure 7 the nonuniform PDF of the AOA for three of the indoor propagation scenarios, namely, the Average Corridor, the Average Lab (NLOS), and the overall channel scenario. In this figure, we see that the PDF of the AOA has 4 peaks corresponding to the characteristics of the exponential distribution within the 2D horizontal plan of the room. The results in this figure have been obtained by evaluating (10) with the optimized channel parameters listed in Table 1.

#### 7. Conclusion

In this paper, we have developed a wideband mobile radio channel model for indoor propagation environments. The starting point for the derivation of the reference channel model was a geometrical scattering model, where we have assumed that an infinite number of scatterers are exponentially distributed within the 2D horizontal plane of a rectangular room. Analytical expressions have been derived for the PDF of the AOA, the PDF of the propagation delays, the PDP, and the FCF. We have shown that the shape of the PDF of the AOA is independent of the position of the transmitter (BS). Both the room length and width have a strong influence on the PDF of the AOA and the FCF. If the room length or width increases, the FCF decays faster if the frequency separation increases. The coherence bandwidth decreases with increasing the room size. The usefulness of the proposed reference channel model has been demonstrated by a close match between the channel statistics of the reference model and measured channels.

An efficient channel simulator has been derived from the reference channel model by applying the SOC principle. It has been shown that the SOC channel simulator matches closely the wideband reference model with respect to the FCF. The designed SOC channel simulator can be used for the performance evaluation of wideband indoor wireless communication systems under realistic propagation conditions.

#### Appendix

#### A. Derivation of the Maximum Propagation Path Length for a Given Value of the AOA

The function , defined bydescribes the propagation path length from the BS to the MS via a single scatterer located at an arbitrary place of the 2D horizontal plane of the room. Owing to the fact that the derivative of with respect to is always positive, that is, we may conclude that is a monotonic increasing function with respect to . That means, for a given value of the AOA , the propagation path length takes a maximum value if the distance from the MS to the scatterer is maximum. The maximum of , denoted by , always occurs if a scatterer is located at the boundary of the rectangle.

For the derivation of the maximum propagation path length, we partition the scattering region inside the room into four subregions , as illustrated in Figure 8. In the Subregion , the AOA is confined to the interval, .

According to the geometrical relationship in Figure 8, we haveSubstituting (A.3) into (A.1) gives an analytical expression for the maximum propagation path length :which depends only on the AOA . The maximum propagation path length for the other subregions , and can be computed analogously. For brevity, we only present the final expression for in (11).

#### B. Determination of the Values for and in (14)

In this appendix, we first determine the maximum value and the minimum value of for the four scattering subregions illustrated in Figure 8.

*Subregion . *. It can easily be seen from (A.4) that is a monotonic increasing function with respect to within the range . We find that takes the minimum value if , whereas equals the maximum value if ; that is,If , then decreases if increases. Thus, the maximum over this range is given by

*Subregion . *. If a given value of the AOA belongs to the range , then we obtain the fixed point by setting the first derivative of [see the second part of the piecewise function in (11)] with respect to to zero. Since the second derivative of with respect to is positive, it follows that has a minimum value at the fixed AOA ; that is,We notice that is a monotonic increasing function if . Therefore, takes a maximum value at ; that is,

If , then decreases if increases. Thus, we have .

*Subregion . *. By studying the first and second derivatives of [see the third part of the piecewise function in (11)], it can be shown that is a monotonic decreasing function within . Thus, we have , whereIf , then , where

*Subregion . *. By conducting similar calculations as for the Subregion , we obtain the range of for the remaining interval of the AOA . For brevity, we present only the final results here. The function has a minimum value at the fixed AOA ; that is,The value of decreases if the AOA changes from to the fixed AOA given above. Within this range, we have . If , then increases and the maximum value occurs at , and thus we have .

An example for the maximum path length is illustrated in Figure 9. It should be mentioned that the path length is always less than or equal to . Since is the minimum value of [see Figure 9], for a given value of the path length , we can conclude that all AOAs satisfy the inequality . However, if the path length , then the inequality can only be guaranteed if or . Here, the AOAs and can be obtained by solving the equationAs shown in Figure 9, . Thus, we select the second piece of in (11) as the expression for the left-hand side of (B.8). Solving (B.8) results finally in the expressions for and as presented in (16a).

For the other given values of , it is necessary to meet the condition . According to Figure 9, we obtained the relations in (B.9), which are presented as follows:The AOAs can be determined similarly by solving (B.8). The expression for on the left-hand side of (B.8) is chosen in accordance with the range of . For example, as , the fourth part of [see (11)] corresponding to the Subregion is selected.

#### Disclosure

This paper is an extended version of “Modeling and Statistical Characterization of Wideband Indoor Radio Propagation Channels,” which was published in the proceedings of the International Congress on Ultra Modern Telecommunications and Control Systems and Workshops (ICUMT), 2010. The starting point in this submitted paper is new. A mixture of exponential functions is used to model the distribution of scatterers. The new mixture exponential functions are more general and include the uniform distribution considered in the ICUMT paper as a special case. This influences all key results of the paper. In addition, the authors in this paper validate the developed channel model with experimental data. The paper shows, by optimizing some of the key parameters with proper methods, that it is possible to use the geometrical channel mode to represent different realistic propagation scenarios.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported in part by the Key Scientific Research Team Cultivation Project founded by College of Information Science and Engineering, Shandong Agricultural University, China. Also, this work was supported in part by the Spanish Ministry of Economy and Competitiveness through the RACHEL Project (TEC2013-47141-C4-2-R), the CARMEN Project (TEC2016-75067-C4-3-R), and the COMONSENS network (TEC2015-69648-REDC).