Abstract
We propose a model of double-directional indoor nonline of sight (NLOS) channels for multiterminal communications. We derive a simple channel matrix that describes input-output relationship for such channels. The multiterminal systems may consist of several terminals that act as amplify-and-forward (AF) relays, where source, relays, and destination have arbitrary numbers of antennas. We complete our model by characterizing the parameters of double-directional channel impulse response of such channels through measurements in indoor environment using 3D synthetic array antenna at 2.5 GHz band. To find out the relation between spatial characteristic of channels in each hop, we observe the direction of arrival (DOA) and direction of departure (DOD) of multipath component signals at the terminal that acts as relay. We find that there are several closely matched azimuths of DOAs and azimuths of DODs which follow uniform distribution in the range of −180° to 180° for elevation around the broadside direction of vertical omnidirectional elements of arrays.
1. Introduction
The significant improvements in spectral efficiency and system reliability through spatial multiplexing-diversity tradeoff promised by the multiple-input multiple-output (MIMO) transmission techniques have triggered their adoption in various broadband wireless systems. Multiterminal communication systems may involve terminals equipped with multiantennas, which are known as multiuser MIMO. Further, these systems grow up by incorporating other terminals as relays between source and destination, in the form of multihop communication. The evaluation of such systems needs a channel model that describes statistically the actual physical phenomenon of radio waves propagation with sufficient accuracy.
The standard MIMO channel models are being developed, as within the European Cooperation in the field of Scientific and Technical Research (COST), and the 3rd Generation Partnership Project (3GPP) which then followed by Wireless world Initiative New Radio (WINNER) project. On the other hand, works on multiple-link indoor channel characterization are reported in some literatures, for example, in [1], where the authors examined the correlation of dual-link MIMO channel matrices between a mobile transmitter and two receivers by studying the similarity in the eigenspace of correlation matrix of those two channel matrices, which might only partially be applicable to the multihop modeling problem at hand. Works reported in [2, 3] were carried out for characterization of indoor multilink large-scale fading parameters. In [2], the authors further investigated the correlation of channel matrices of two-hop MIMO channels. They expressed the channel matrix between source and destination that consists of two MIMO hops as the Kronecker product of the two MIMO channel matrices, so that their expression is restricted to the condition that scatterers around the first and second hops that are uncorrelated. However, the uncorrelated scattering condition is given without further insight into the multipath phenomenon that actually occurs in radio channels in a multihop system. Also, works reported in [1–3] have not considered a general model of multihop channels, that is, channels with arbitrary number of hops and antennas at each node, and joint characterization of spatial-temporal parameters of the channels.
An analytical result in [4] showed that radio wave propagation over MIMO amplify-and-forward relay channel is influenced by the propagation of radio wave in individual hops and the responses of the relays. However, the model assumes that the array responses of transmitter and receiver of a relay are orthogonal, which does not hold for isotropic or omnidirectional elements. Furthermore, no verification was attempted to the assumption by field measurement. Likewise, a recent literature [5] proposed a unified channel model for multilink cooperative MIMO. Based on geometry-based stochastic modeling, the model can be applied to many cooperative scenarios by adjusting its parameters. The paper also derived spatial correlation between any two links. However, the channel gains and their spatial correlations obtained are determined by the type of links in a multilink scenario, where the parameters have been adjusted according to the physical properties of the local scatterers around each end node of the links, which still needs further investigations.
Our contribution includes firstly the derivation of a simple expression of channel matrix that describes input-output relationship of indoor multihop communication involving amplify-and-forward (AF) relays as illustrated in Figure 1, assuming that there is no direct link between the source and the destination. The expression in general is that the number of relays and the numbers of antennas at the source (), relays (), and the destination () are arbitrary. It also incorporates any correlation between model parameters that may exist in the channel. The model is similar to that reported in [4], but in our expression, we include the expression of delays as well as that of DOD and DOA of MPCs in each hop of the channel, and we do not consider any assumptions for the array response at relay.
In a two-hop link, the first and the second hops that share a common node, that is, a relay in Figure 1, which can be a base station or access point in multi user MIMO, might also share the same interacting objects surrounding the relay. This situation might lead to the appearance of multipath component signals (MPCs) that arrive at and depart from relay at the same directions, which might cause an extent of correlation between the hops. The common assumption, as taken in [2], in which the first and second hop channel matrices are uncorrelated, still needs to be verified. Our second contribution in this work is an investigation into the spatial characteristics at a relay. We observe how DOA and DOD at a relay are related to each other to obtain a clear understanding about the correlation between the two hops described previously.
To complete our model, we derive statistical characterization of model parameters based on the results of a measurement campaign conducted in our laboratory building. Although several statistical characteristics of channel parameters have been reported in the literature, for example, in [6–8], according to the best of our knowledge, the joint statistical characterizations of all double-directional indoor channel parameters have not been reported yet. By using a 3D array antenna at the transmitter and at the receiver instead of using a linear or planar array, we could estimate completely the directional parameters at both sides, that is, azimuth and elevation, from −π to π and from 0 to π, respectively, without any ambiguity in elevation. A 3-D cube geometry with isotropic elements theoretically has isotropic array factor desired for better DOD and DOA estimations [9]. Although we use omnidirectional antennas as array elements, better estimation results compared to those obtained by using linear or planar arrays should still be achieved. By applying a high resolution channel parameter estimation first proposed in [10], we could obtain the directions, delays, and powers contained in MPCs simultaneously, so that we could observe the relation among parameters simultaneously too. Thus, our third contribution is a new joint statistical characterization of channel parameters and channel power spectra that would be meaningful to gain more complete and accurate statistical model, particularly for the 2.5 GHz band in NLOS indoor environment.
The paper is organized in the following way. In Section 2, we derive double-directional multihop channel model. Section 3 describes the measurement system and measurement environment, followed by signal processing technique used to estimate channel parameters in Section 4. In Section 5, we discuss the statistical characteristics of model parameters, including a look into the relation between DOA and DOD at relays. In Section 6, we present the power spectra of the channels. Finally, some concluding remarks are given in Section 7.
2. Multihop Channel Model
The multihop channel considered in this paper is illustrated in Figure 1, where source and destination employ and antennas, respectively. The channel consists of parallel relays, each of which consists of serial relays, so that the number of relays is in total. This relay structure is general so that simple arrangement of serial or parallel relays can be derived from it. Each relay has a gain of and is equipped with antennas. Each hop in Figure 1 generally represents an MIMO link. On link (), if represents output discrete time signal vector of , where and represents time index, then following the description in [11], received signal vector at relay R11, can be expressed as follows: In (1), denotes the number of MPCs at link (), and is delay of the th MPC with respect to the first MPC, where the delay of the first MPC is 0 ns. is noise vector at R11 antennas, and is an MIMO channel matrix that consists of complex channel gain between all possible transmit-receive antenna pairs encountered by the th MPC in link () and can be written as where is the number of antenna at R11.
According to [12], the double-directional impulse response at link () can be written as where , , , and are azimuth and elevation variables with respect to R11 receiver and array, respectively, and is delay variable. Also , , , , and are azimuth and elevation variables of DOA at R11 receiver, delay, azimuth, and elevation of DOD at of the th MPC, respectively, while is complex amplitude of the th MPC. The definition of elevation and azimuth variables at any transmitter and receiver array in a single-hop MIMO channel and the example of azimuth and elevation of DOD and DOA of the th MPC related to array orientation are depicted in Figure 3.
The matrix elements in (2) can be obtained by computing double-directional channel impulse response at each pair of and R11 receiver antenna given by assuming that the array at and R11 is narrowband and all MPCs incident to array are plane waves. In (4), denotes radiation pattern of the th element of array, represents the position vector of the th element of array with respect to array mass center, and is the wave number vector given by where is wavelength at the center frequency of radio wave and stands for transpose operation. Meanwhile, and are wave number vectors at R11 receiver and related to the th MPC, respectively.
If the complex steering vector of any array containing elements for the th MPC is denoted by and for , then (2) can be expressed as where and are complex steering vectors of the th MPC at R11 receiver and , respectively.
The complete expression of channel matrix of link () that combines all MPCs and includes all MPCs parameters is simply where and are complex steering matrices at R11 receiver and , respectively. Parameter is diagonal matrix that represents complex amplitude and delay factors of MPCs at link () and is expressed as
Using (7) and by dropping time index to simplify the expression, input-output relationship in (1) can be written as
On the the second hop of the first serial relay in Figure 1, the signal vector radiated by R11 transmitter is . If this second hop channel is denoted by (), then the signal vector received by relay R12 is , where the is noise vector at R12 receiver. Then, can be written as where is the accumulated noise. The combined two-hop channel matrix is
If we suppose that the number of MPCs on link () is and the number of antennas at relay R12 is , then and are complex steering matrices at R12 receiver and R11 transmitter, respectively. Parameter is diagonal matrix that represents complex amplitude and delay factors of MPCs at link ().
Following this derivation, if the number of such serially arranged relays is , then the matrix of ()-hop channel is where and are complex steering matrices at and is a diagonal matrix that represents complex amplitude and delay factors of MPCs at link (), respectively, while is the number of MPCs on that link. The channel matrix in (12) is generic in that every ()-hop channel can be considered as an extension of -hop channel by inserting a matrix factor similar in shape to (7) into (12) to represent the additional relay.
Then, the signal vector received at from these serial relays can be written as where is the accumulated noise in the first serial relays.
The other sets of serial relays have the same expression as (12). Therefore, the matrix of multihop channel with relay structure depicted in Figure 1 can be written as and signal vector received at is
While the expression in (12) and (14) can be applied to multihop links involving AF relays with any number of hops, the statistical characteristics of any single link in the following sections apply also for multihop links with detect and forward (DF) relays.
3. Measurement System and Environment
The measurement system is based on vector network analyzer (VNA) focusing on 2.17–2.5 GHz band and comprising 201 points of frequency [13] with the measurement configuration as depicted in Figure 2. We use a 3-D cube-shaped synthetic array antenna at both transmitter and receiver. Each of the eight elements is a wideband biconical type antenna with 3 dB elevation beamwidth of 60° and mounted vertically at the cube corners. Measurements in anechoic chamber have shown that the antenna pattern in azimuth and elevation variables did not change significantly over the measured frequency band. The standard deviations of normalized power measured in all frequency points averaged over the azimuth and elevation samples are 0.04 and 0.09, respectively, which prove the statement mentioned previously. The distance between array elements is 6 cm and is equal to half wavelength of the upper frequency to avoid spatial aliasing caused by spatial undersampling over the array. The resulting dimension of the array fulfills the requirement of narrowband array, so that there is no significant delay difference experienced by an MPC when it arrives at the array elements.
The cube array orientation and the definition of azimuth and elevation variables for transmitter and receiver in a single-hop link are depicted in Figure 3, where the th MPC departs from transmitter at azimuth and elevation variables of and arrives at receiver at . As shown in Figure 3, the reference direction of 0° azimuth of DOD and DOA is in line to their respective axes when the arrays are facing each other.
The measurement floor plan is depicted in Figure 4. There are two measurement scenarios. The first is conducted to obtain statistical characteristics of double-directional channel parameters. This scenario consists of single-hop double-directional channel measurement, where the transmitter and receiver in each hop are described in Table 1, and the numbers indicate transmitter and receiver locations in Figure 4.
The second scenario is conducted to investigate the relation between DOA and DOD at a relay in a two-hop channel. Source, relay, and destination in each of the two-hop channel are described in Table 2. The two-hop channel is composed of two one-hop channels measured in sequence, where the “relay” serves as a receiver and a transmitter in the first and second hops, respectively. The relation between DOA and DOD at relay in the two-hop channel can be used to describe the relation between DOA and DOD at any relay in a multihop channel. This is because a multihop channel can be decomposed into several two-hop channels, and any additional relay or hop can be represented by a matrix factor as in (7).
This second measurement scenario concerns the spatial characteristics at each relay of a multihop channel with an arbitrary number of hops. For example, in the measurement realizations depicted in Table 2 and referring to Figure 4, the three-hop channel consisting of nodes 1-2-3-4 is decomposed into two two-hop channels, that is, 1-2-3 and 2-3-4, and also measurement realizations of four-hop channel consisting of nodes 5-6-7-8-9 are decomposed into three two-hop channels, that is, 5-6-7, 6-7-8, and 7-8-9.
From the scenario in Table 2, we obtain seven two-hop double-directional channel transfer functions. In all measurement realizations, the distances between transmitter and receiver and those between transmitter or receiver and wall, glass door, or metal cupboard are kept in far field region to fulfill the requirement of plane wave measurement.
The environment under consideration is a laboratory room which represents a small laboratory or office with its furniture and equipments of any shapes, sizes, and materials. The common furniture materials are wood, plastic, and glass, while the common equipment material is metal. There are also some glass windows and doors in the environment. Transmitter and receiver array’s height is 80 cm above the floor to approximately represent users’ equipment or access point’s height. In all measurement scenarios, the transmitter and the receiver are in different rooms separated by brick wall or in the same room but obstructed by metal cupboards in between.
The use of synthetic array antenna requires that the environment should be kept static during one double-directional measurement of each single-hop channel in the first measurement scenario. Meanwhile, in the second scenario, the static condition must be fulfilled during two one-hop channel measurements described previously. This is achieved by doing the measurement activities in the evening or during weekends. In addition, by using synthetic arrays, the measured channel transfer functions are free from mutual coupling effect. To avoid interference from any wireless local area network signals, we turn off access points around measurement locations during the measurement activities and by using a spectrum analyzer to assure that there is only noise in the measurement band when the system is turned off.
Measurement using cube-shaped synthetic array antenna is carried out by placing an array element on a corner of the transmitter cube and by placing another array element on a corner of the receiver cube. The elements of transmitter array and that of receiver array are connected to VNA at port 1 and port 2, respectively, and then the channel transfer function is measured. This measurement is repeated by moving the element of receiver array to the other cube corners sequentially so that we can obtain eight channel transfer functions between the first element of transmitter array and each element of the receiver array. To obtain one double-directional channel transfer function, this procedure is repeated for the other transmitter array element positions. The results are arranged in a matrix of 64 × 201 size that represents 201-frequency point channel transfer function between 8 × 8 pairs of transmitter and receiver antennas.
4. Data Post Processing
There are some methods that can be used to estimate channel parameters from measured data, such as the frequency-domain space generalized expectation-maximization (FD-SAGE) [14]. FD-SAGE was developed for estimating single-input multiple-output (SIMO) channel parameters from data measured in the frequency domain. It was based on previous work reported in [10, 15], where [10] discussed general SAGE algorithm and [15] studied the implementation of [10] for multipath channel parameter estimation.
The computational load in FD-SAGE depends on the number of MPCs assumed to exist in the channel and the number of parameters of each MPC that has to be estimated. In an indoor environment, naturally there are many MPCs. Modeling MIMO channels relates to the larger number of MPC parameters that have to be estimated when compared to that in SIMO channels. As indicated in (3), the set of parameters of the th MPC is . Accordingly, we modify some procedures in FD-SAGE to achieve less computational load. In order to do so while maintaining the accuracy, we divide our procedure into two stages, that is, the IFT stage and FD-SAGE stage. We call this method as IFT-FD-SAGE.
The first step in IFT-FD-SAGE consists of applying window and inverse Fourier transform (IFT). The result of this process is a delay domain impulse response of a channel between the th and mth antennas at receiver and transmitter, respectively. For a pair of transmitter-receiver with and antennas, respectively, we obtain impulse responses. A threshold is then applied to those impulse responses, where the components below it are considered as noise and removed. In this work, the threshold is set to 3 dB above the highest side-lobe level of the window used, which is the Hamming window.
Further, we use the first criterion in [16] to determine if an MPC is not a spike noise; that is, all of impulse responses are examined for consistent emergence of an MPC. Inconsistent spikes are considered as noise and removed from the MPCs set. The emergence of a spike is considered as consistent if it is greater than or equal to % of the impulse responses, where in our procedure we choose 75% of 64 impulse responses. Then, the power of MPCs is calculated as the mean of the corresponding power values in all impulse responses. By completing the first stage, we have estimated the number of MPCs, the delay of each MPC, and its power.
In the second stage, we apply FD-SAGE to estimate four parameters that have not been estimated yet, that is, . Because we have determined the number of MPCS through thresholding in the first stage, unlike in [14], we do not any longer require MPC removal after FD-SAGE is finished.
We verify our method by comparing the results of the method when it is applied to the measured data with their counterparts obtained from a simple 2D ray tracing in one of the rooms in floor plan of Figure 4, that is, the room in the middle that contains measurement points 2, 7, and 10. The results are depicted in Figure 5. There are five paths obtained from IFT-FD-SAGE and those obtained from ray tracing that closely matched. These paths consist of the direct, single-reflection, and double-reflection paths. The other paths which encounter more reflections are not considered in the ray tracing so that the results of the algorithm which have longer delay cannot be confirmed. In Figure 5 we can see that path number 2 and path number 6 have larger error (defined herein as absolute difference between the value of parameter obtained from ray tracing and that obtained from the algorithm) compared to the other paths. This is mainly because the delay difference between path number 2 and path number 6 is below the delay resolution, that is, 3 ns according to the measurement bandwidth, whereas based on ray tracing calculation, the delays of path number 2 and path number 6 are 39.2 ns and 41.8 ns, respectively.
(a)
(b)
On the other hand, we also perform simulations to observe the performance of the method and its capability to estimate parameters of a sufficiently large number of MPCs. For example, in case of estimating the parameters of twenty MPCs, the method works well if the SNR is high enough, that is, 20 dB in our simulations. Following this simulation result, we keep the minimum SNR at 20 dB in all measurement realizations, except one measurement whose the maximum SNR achieved is 18 dB.
By using the SAGE-based method, the angular resolution can be set high enough although the number of array elements is limited. In this work, we set it at 1°.
5. Statistical Characteristic of Channel Parameters
It was introduced in [17] that the channel statistic depends on the product of probability density function (pdf) of the channel parameters and power spectral density of the channel. In this section, we present pdfs of channel parameters, whereas power spectral density of the channel will be discussed in Section 6.
From the first measurement scenario described in Section 3, we obtain 23 single-hop double-directional channel transfer functions in total. The number of MPCs and their parameters in each hop were obtained by applying IFT-FD-SAGE algorithm described in Section 4 to these channel transfer functions. Based on these results, we characterize statistical distribution of channel parameters described in (3).
5.1. Marginal Distribution of Parameters
Marginal distribution of parameters is obtained by applying Kolmogorov-Smirnov (K-S) goodness of fit (GOF) test at significance level of 0.05 to the empirical cumulative distribution functions (cdfs) of parameters. The results are listed in Table 3, where in the table, and are the mean and the standard deviation of parameters, respectively, is the maximum distance between empirical cdf and theoretical cdf, is the value of K-S statistic at significance level of , and is the number of samples within the range of parameter values used to calculate the empirical and theoretical cdfs. The empirical data are well fitted to a certain theoretical distribution if is lower than [18].
We observe that many strong MPCs arrive at the array around the azimuth of 0°. There is also another group of weaker MPCs arriving around 180°. These two groups appear to be separated by the even weaker powers at around 120° and 240°, as can be seen in Figure 6; however, the MPCs coming around 180° have lower powers than those coming around 0°. We also observe that the departure characteristic of MPCs is similar to their arrival characteristic above. These results indicate that in the measured environment, many MPCs whose powers are higher than those coming to or depart from around 120° and 240° seem coming to or departing from the front of and behind the array. At the receiver site, these MPCs are likely to come from the direct path through the wall, diffracted by objects or building structures near the direct path, and then penetrate the wall, and some of them whose powers are lower and arrive behind the array may come then after being reflected by objects behind the arrays. On the other hand, the small number of MPCs around 120° and 240° indicates that the MPCs reflected or diffracted by object or building structures in those directions may experience longer propagation path and have lower power than the others.
Therefore, in order to properly characterize DOA and DOD azimuths, we divide those parameters into intervals of –120° to 120° and 120° to 240° for the array orientation depicted in Figure 3. In those two intervals, DOD and DOA follow a curtailed Gaussian distribution. Each of DOD and DOA azimuths has mean values around 0° and 180° for the two interval described previously. This gives a new understanding that MPCs are actually not distributed uniformly, although transmitter and receiver are surrounded by objects and the radio links are NLOS. These results are different from those reported in [7], where for NLOS environment, DOA azimuth of intercluster MPCs is distributed uniformly in −180° to 180°. We compare this result with intercluster azimuth in [7], although we do not model MPCs as clusters for the reason that will be discussed in Section 6.
We also observed statistics of the number of MPCs contained in the two intervals. We found that on average, the number of MPCs in each hop is 25, while the number of MPCs in front of the array, that is, in the range of –120° to 120°, is almost four times that of behind the array, that is, in the range of 120° to 240°. In these two intervals, they follow the Gaussian distribution.
Statistics of those six normally distributed parameters are listed in Table 3 and their marginal pdf can be written as, where can be one of the above-described parameters, and are the mean and the standard deviation, respectively, and are the lower and upper bounds values of parameter .
On the other hand, DOA and DOD elevations are better fitted to asymmetric Laplace than to Gaussian distribution, although the two are accepted according to their K-S statistics. Their pdfs can be expressed as where can be or , that is, DOD or DOA elevation, and also can be or , that is, mean DOD or DOA elevation. , where can be or , that is, standard deviation of DOD or DOA elevation. Figure 7 depicts the empirical pdf of DOA and DOD elevations superimposed with their theoretical Laplace pdf.
As described in Section 3, we use vertical biconical antennas for array elements and we set the same height for transmitter and receiver antennas. Accordingly, most of the MPCs should depart from or arrive at the antenna from the broadside direction of the array element. The Laplace pdf with the mean value around 75° and the small standard deviation suggests that most of MPCs are scattered by objects whose positions are higher than the array, for example, lamp shades, the top of metal cupboards, equipments in high shelves, and so forth. We see from Figure 7 that MPCs departing to or coming from elevation below the broadside direction of the array element are few. This gives information that MPCs which are resulted from reflections at the floor are few.
Time of arrival or delay of MPCs is measured relative to that of the first MPC, which is 0 ns. Result of K-S GOF test confirms that the delay distribution fits to the exponential distribution.
The pdf is expressed as where is mean of delay and the value can be seen in Table 3. The delay distribution of MPCs cluster reported in [7] for NLOS environment is also exponential with mean 53 ns.
5.2. DOA and DOD Relation at Relay
Our second contribution is the investigation of DOA and DOD relation at relay. It was based on the idea that any relay in multihop communication may be surrounded by objects that will produce the same direction to the wave coming to the relay or leaving it. This gives an insight into the presence of scatterers around relay common to the incoming and outgoing waves and about the physical properties of scatterers around the relay. As channel matrices correlation between incoming and outgoing hop of a relay should depend on the common scatterers of the incoming and outgoing waves, the spatial characteristics discussed below relate to the correlation between those channel matrices. The investigation was performed by measuring the channel with the second measurement scenario described in Section 3.
As was described in Section 5.1, when considered in elevation variable, most of the MPCs arrive at or depart from array near the broadside direction of the array elements. This fact motivates us to focus the investigation of the closely equal DODs and DOAs only on MPCs with elevation of DODs and that of DOAs around this direction. By applying this procedure, we neglect elevation variable and the investigation of DOA and DOD at relays is carried out only in the azimuth variable. The closely matched azimuths are defined as the DODs and DOAs azimuths that have magnitude of difference value below or equal to 2°.
In reference to Figure 3, closely matched DOAs and DODs azimuths at represent the same directions with respect to array at . MPCs in those directions arise as radio wave scattered to by any objects, and those same objects scatter radio wave from .
According to K-S GOF test at the significance level of 0.05, closely matched DOAs and DODs azimuths are uniformly distributed in the range of −180° to 180°, with an average of −6° and standard deviation of 109°. This result indicates that in the environment considered, any objects will act as the same scatterer around for the incoming and outgoing MPCs with the same probability.
From the seven two-hop links considered, we found that the number of closely matched DOAs and DODs azimuths is uniformly distributed in the range of 2 to 6, with an average of 4 and a standard deviation of 2. It suggests that on average 4 of 25 MPCs arrive at relay after scattered by some objects. Then these MPCs depart from the relay to the same objects and scattered again by these objects to the direction of destination. The different pairs of arriving and departing waves with closely matched azimuth DOAs and azimuth DODs observed originate from different scatterers, as indicated by their values. The low average number of these pairs indicates the low correlation between channel matrix of incoming hop and that of outgoing hop of a relay in the multihop channel considered. However, the statistics reported in this section need more measurements to be confirmed because we only observe 7 pairs of dual-hop link.
5.3. Joint Distribution of Parameters
By using the FD-SAGE-based algorithm, we obtain all parameters of MPCs, that is, , , , , , and simultaneously for . There is no assumption that these parameters are independent to each other, and hence, we investigate the dependence between any pair and out of the five parameters above. Since the random parameters evaluated are not necessarily Gaussian, we therefore resort to the least-squares independence test (LSIT) [19], rather than evaluating the correlation coefficients. The test is performed on parameters of 23 single-hop measured channels. The results are depicted in Table 4, where we divide the azimuths of DOD and that of DOA into two intervals as in Section 5.1. However, a DOD azimuth in front of the array may form a pair with a DOA azimuth in front of or behind the array, and also one behind the array may form a pair with another in front of or behind the array, so that the sample spaces of DODs and DOAs in those two intervals are not the same, and we cannot define any joint distributions between those parameters. Accordingly, we suggest considering them independent to each other. We choose a reasonable significance level of 0.1 to reject the null hypothesis, in which and are independent, where the value of LSIT is listed in Table 4.
The joint pdf of DOD and DOA elevations, whose marginal pdf is asymmetric Laplace, is a bivariate asymmetric Laplace expressed as [20] where , , is covariance matrix of and and denotes inverse operation. is the third kind of Bessel function of zeroth order. The values of and are depicted in Table 3. The covariance matrix can be constructed if the correlation coefficient between and are known, and the standard deviation of and are also known. Correlation coefficients between and calculated from all measured channels appear to be random in the range of −0.5 to 0.5 following a Gaussian distribution with a mean of 0.13 and a standard deviation of 0.25. Meanwhile, the values of and can be seen in Table 3.
Following the derivation in [7], the joint distributions of DOD elevation and delay and that of DOA elevation and delay are expressed as conditional pdfs, that is, pdfs of DOD elevation and that of DOA elevation conditioned on a specific delay interval , where for and is the number of delay steps of width. The relation between intervals and is described in Table 5. The pdf of elevation , which can be DOD elevation or DOA elevation conditioned on delay , is expressed as
A set of empirical DOA elevation pdfs conditioned on are depicted in Figure 8, where ns are chosen to obtain a moderate interval in the 0 to 120 ns range of delay. All of the empirical pdfs in Figure 8 are well fitted to Laplace pdf according to K-S GOF test. The average of means is 78°, in agreement with the mean of DOA elevation marginal distribution. The standard deviations tend to decay with increasing , which implies that they tend to decay with increasing delay. The pdfs of DOD elevation conditioned on have similar characteristics. The conditional pdf plots are not shown here, but the results in general are similar to those for the DOA given in Figure 8. Figure 9 shows the variation of standard deviations of DOA and DOD elevation pdfs conditioned on .
(a)
(b)
(c)
(d)
(e)
(f)
Pdfs of DOD and DOA elevation conditioned on can be written as where , , and are, respectively, the mean and the standard deviation of DOD and DOA elevations for delay interval centered at .
The variation of standard deviation on delay can be approximated by the following functions:
The expression in (17) now can be written as, where the value of can be approximated by the mean of in Table 3.
Because DOD and DOA elevations are dependent and are expressed as joint pdf in (19), then the joint pdf of DOA elevation and DOD elevation and delay is expressed as joint DOD and DOA pdfs conditioned on a certain delay interval as follows: where now is covariance matrix of DOD and DOA elevations at , obtained by calculating the standard deviation of DOD and DOA elevations at in (22) and the correlation coefficient as described before.
The decreasing standard deviation of elevation pdfs as a function of delay indicates that the MPCs which arrive with longer delay most probably are not caused by multiple reflections by objects in which the position is higher or lower than the array. The MPCs that depart to or arrive from objects above or below the arrays most probably present in the delay below 40 ns, as indicated by Figure 8.
Then, joint pdf of all parameters can be expressed as where and are described in (16) for each azimuth interval, as in (24), and as in (18). The expression in (25) completely describes joint pdf of all double-directional parameters appearing in (3).
6. Channel Power Spectra
Power spectra of the channel are considered separately at transmitter and receiver sides. This is because in each side, they represent the total power spread in the delay-angular domain. Power spectrum in delay-angular variable considered at transmitter and receiver sides can be expressed as respectively. The operator and denote expectation and absolute values, respectively. Those two power spectra are related to each other by the delay variable and satisfy
Power spectrum considered in any certain variable could be obtained by integrating (26) over other variables. For example, power spectrum considered in delay variable, that is, power delay spectrum (PDS), is calculated by
In the following, we investigate two kinds of power spectra. The first is power spectra of single variables defined in (28), namely, empirical power delay spectra (PDS), departure and arrival power elevation spectrum (DPES and APES), and departure and arrival power azimuth specta (DPAS and APAS).
However, those power spectra do not describe the relation between power spectrum in a certain variable and that in another. Therefore, we need a complete description of power spectrum involving all variables, that is, delay, elevation, and azimuth, defined in (26), which is called power delay-angular spectrum. Recalling the results of LSIT to the double directional parameters in Section 5.3, the only dependencies are between DOA elevation and delay, DOD elevation and delay, and between DOA and DOD elevations. Because the power spreading in the double-directional variables depends on the probability of the emergence of MPCs in those variables, that is, their pdfs, the power spreading should also consider the known dependencies above. Hence, to obtain power delay-angular spectrum, in the following, we only investigate power delay-elevation spectrum and consider that the departure and arrival power azimuth spectra are independent to delay. This is the second power spectra that will be discussed.
In this work, the expectation in (26) is estimated by averaging the 23 single-hop double-directional impulse responses obtained previously. This is because the desired power spectra are those that represent the characteristic of the measured environment, not those that represent a local averaging over a certain area of radius about several wave lengths.
6.1. Power Delay Spectrum (PDS)
Empirical PDS is calculated as follows: where represents the average power of th MPC, is the number of measurements, and is the power of th MPC in the th measurement. The empirical PDS should satisfy as in (27).
The obtained empirical PDS is depicted in Figure 10. It fits to an exponential function expressed as Parameter is delay spread and defined as the square root of the second central moment of . The estimated value of is 25 ns.
In [6, 21], PDS of clusters and rays inside clusters are modeled by decaying exponential function. Although in this work we do not consider MPCs cluster, generally our model agrees with those results, because PDS of clusters is the envelope of PDS of rays inside clusters. In the cluster-based model, the center of cluster is MPC whose power is the highest. If we could observe rays inside a cluster, its power should not be higher than the cluster center. According to this, we believe that the PDS in (30) is more relevant to compare to the PDS of clusters. Also, the comparison of pdf of delay and that of DOA azimuth in Section 5.1 with those reported in [7] is based on the same reason above.
6.2. Departure and Arrival Power Elevation Spectra (DPES and APES)
DPES and APES are power elevation spectra considered at transmitter and receiver sides, respectively. The calculation procedure of empirical DPES and APES is similar to that used for empirical PDS.
The obtained empirical APES is depicted in Figure 11. The empirical DPES is not shown; however, its shape is similar to the empirical APES. The empirical APES and also the empirical DPES follow asymmetric Laplace function expressed as where may denote or , that is, DOD or DOA elevation, , and can be or , that is, departure or arrival elevation spread which is defined as the square root of second central moment of DPES or APES, respectively, and also can be or , that is, the central moment of DPES and APES. The values of and are 75° and 77°, respectively, whereas the values of and are 14° and 17°, respectively. In (31), the factor is normalized constantly to satisfy (27).
These results show that much of the power distribution in elevation is concentrated in the narrow region which is slightly higher above the broadside direction of the array element, in agreement with the result in Section 5.1, where the elevation departure and arrival of MPCs are mostly in the these directions.
6.3. Departure and Arrival Power Azimuth Spectra (DPAS and APAS)
Empirical APAS can be seen in Figure 6, whereas empirical DPAS has similar shape to it. As described in Section 5.1, visually we can see that there are two groups of MPCs, that is, having DOD and DOA around 0° and 180°, which tend to have higher power then MPCs groups having DOD and DOA around 120° and 240°. Based on these results, we do not consider that DPAS and APAS are uniform in −180° to 180°. Instead, we divide the intervals of DPAS and APAS into −120° to 120° and 120° to 240° and fit the empirical DPAS and APAS with the Gaussian function.
The Gaussian fit shows much lower sum of square error (SSE) than uniform fit, and, accordingly, we model DPAS and APAS as follows: where can be or , that is, DOD or DOA azimuth, respectively, whereas and can be or and or , that is, the central moment of DPAS or APAS in −120° to 120° and in 121° to 239°, respectively. Their values are −15°, −15°, 176°, and 177°, respectively. Also, and can be or and or , that is, departure or arrival azimuth spread in the two intervals above, respectively, the values of which are 55°, 56°, 31°, and 30°. The constants and in (32) are required to fulfill .
These results are different from the power azimuth spectrum reported in [7], where intercluster APAS in NLOS environment is uniform in −180° to 180°. Our results can be understood when we recall that in LOS environment; the power azimuth spectrum in several literatures, for example, [4, 5], is modeled by the Laplace function. When the links are obstructed or separated by wall, the Laplace function may not become uniform function directly, but part of the power which is higher than the others may still fill the region in front of the arrays, and one, which is lower than the first, fills the region behind the arrays, where the mechanism of filling those regions has been discussed in Section 5.1.
6.4. Departure and Arrival Power Elevation-Delay Spectra
Following the method in [22], departure and arrival power elevation-delay spectra are expressed as partial departure and arrival power elevation spectra (PDPES and PAPES) at a certain delay interval and are given by where relates to a specific delay interval defined in Table 5. A series of empirical PAPES is shown in Figure 12. Empirical PDPES is not shown here, but its shape is similar to that of PAPES. Figure 12 shows that PAPES on all specified delay intervals could be approximated by asymmetric Laplace function. However, the elevation spreads seem to decrease when increase. As in (31), asymmetric Laplace function is now used to model PAPES as follows: whereas before, might represent or , while denotes the central moment of PDPES or PAPES at interval . does not vary significantly with and can be approximated by the central moment of PDPES or PAPES described in Section 6.2. , where denotes the departure or arrival elevation spread at . In (34), the notation means that is proportional to . The variation of on implies that varies according to and is shown in Figure 13. They can be approximated by exponential functions expressed as for PAPES and PDPES, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
To obtain a complete description of power delay-elevation spectrum, we observe the variation of power contained in each delay interval as a function of . As can be seen in Figure 12, the total power contained in each delay interval seems to vary, which fits an exponential function well in Figure 10.
From these results, we can model power delay-elevation spectrum (PDES) from (30) and (34) as follows: where varies according to the variation of in (35). The values of are given in Section 6.1.
6.5. Channel Power Delay-Angular Spectrum (PDAS)
In this section, the results of the previous power spectra are used to describe the general power delay-angular spectrum of the channel. The power delay-angular spectrum in (26) can be expressed as partial power angular spectrum at a certain delay interval given by where has been described in Section 5.3, and as before, may denote or , and similarly, may be or .
Examined at either receiver or transmitter, channel power angular spectrum is given by
Because and are independent, therefore , and partial power angular spectrum becomes . However, . Hence, partial power angular spectrum can be written as
Since implicitly describes the variation of power elevation spectrum on delay as in (36), therefore, (39) also describes power angular-delay spectrum and can be expressed as where and are given by (32) and (36), respectively. In Table 6, we present the summary of the channel power spectra discussed in Section 6.
7. Conclusion
An expression for channel matrix that describes input-output relationship of double-directional indoor Multiterminal communication involving AF relays has been presented. The expression in general is that the number of relays and the number of antennas at transmitter, relay, and receiver are arbitrary. If there is any correlation between model parameters in each hop or in multihop channel, the correlation can be incorporated in the expression.
The spatial characteristics at a relay of Multiterminal communication have been described by closely matched DOA and DOD azimuths at the relay. Based on the measurement results, we found that there are some closely matched DOA and DOD azimuths which follow uniform distribution in –180° to 180° for elevation around the broadside of the vertical omni-directional array element. The numbers of closely matched DOAs and DODs azimuths are uniformly distributed in the range of 2 to 6 with an average of 4.
We present new statistical characteristics of double-directional channel model parameters and channel power spectra. The statistical characteristics are based on intensive measurements using cube synthetic array antenna so that the double-directional angular parameters obtained consist of elevation without ambiguity. The proposed statistical characteristics of channel parameters and channel power spectra are expressed by including their variations on delay.
The proposed delay distribution is exponential, which agrees with the delay distribution proposed before in the literatures. The marginal distributions of DOD and DOA elevations are asymmetric Laplace with a mean of about 75° and a standard deviation of about 15°. They are dependent to each other and can be modeled by bivariate asymmetric Laplace pdf characterized by their covariance matrix. Further, the DOD and DOA elevations are dependent to delay so that their joint distribution is characterized into partial joint pdf, that is, joint pdf considered in a specific delay interval, where their standard deviations decrease with increasing delay. Meanwhile, DOD and DOA azimuths follow the Gaussian distribution for the direction in front of and behind the array with means near 0° and 180°.
Each of the channel power spectra considered at transmitter and receiver is expressed as the multiplication of its power delay-elevation spectrum and power azimuth spectrum. The power delay-elevation spectrum considered is modeled by the Laplace function whose magnitude and elevation spread vary with delay; meanwhile the power azimuth spectrum is modeled by the Gaussian function in the region in front of and behind the array.