text:
Member, IEEE Marco Moretti

Member, IEEE Luca Sanguinetti

Fellow, IEEE Xiaodong Wang

Spatial Multiplexing MIMO-OFDMA Systems

come

I. INTRODUCTION

Dynamic resource allocation in multiple-input multiple-output (MIMO) systems based on orthogonal frequency-division multiple-access (OFDMA) technologies has gained considerable research interest [START_REF] Li | Resource allocation for clustered network MIMO-OFDMA systems[END_REF]. In most cases, subcarriers are assigned to the active users in an exclusive manner without taking advantage of the multi-user diversity offered by the spatial domain. A possible solution to exploit the spatial dimension is to make use of space-division multiple-access (SDMA) schemes, which allow the simultaneous transmission of different users over the same frequency band. The main impairment of SDMA is represented by multiple-access interference (MAI). In downlink transmissions, MAI mitigation can only be accomplished at the BS using pre-filtering techniques. The most common approach for interference mitigation is zero-forcing (ZF) linear beamforming, which relies on the idea of pre-inverting the channel matrix at the transmitter. Another approach is represented by the block-diagonalization ZF (BD-ZF) scheme originally proposed in [START_REF] Spencer | Zero-forcing methods for downlink spatial multiplexing in multiuser MIMO channels[END_REF]. Particular attention has been also devoted to dirty paper M. Moretti and L. Sanguinetti are with the University of Pisa, Dipartimento di Ingegneria dell'Informazione, Italy ( {marco.moretti,luca.sanguinetti}@iet.unipi.it). X. Wang is with the Electrical Engineering Department of Columbia University, USA (wangx@ee.columbia.edu) coding (DPC) techniques [START_REF] Tran | Multiuser diversity for successive zero-forcing dirty paper coding: Greedy scheduling algorithms and asymptotic performance analysis[END_REF] even though their implementation is still much open. A possible solution in this direction is represented by Tomlinson-Harashima precoding (THP), which can be seen as a one dimensional DPC technique [START_REF] Fisher | Precoding and signal shaping for digital transmission[END_REF] and has been widely used in the downlink of single-user and multiuser MIMO systems [START_REF] Stankovic | Successive optimization Tomlinson-Harashima precoding (SO-THP) for multi-user MIMO systems[END_REF]- [START_REF] D'amico | Tomlinson-Harashima precoding in MIMO systems: A unified approach to transceiver optimization based on multiplicative Schur-convexity[END_REF]. In combination with pre-filtering, another way to deal with interference in SDMA-OFDMA systems is user partitioning, which basically consists in properly selecting the set of users transmitting on the same subcarriers. As illustrated in [START_REF] Maciel | On the performance, complexity, and fairness of suboptimal resource allocation for multiuser MIMO-OFDMA systems[END_REF], a common approach is to first group together users whose channels have low spatial cross-correlation and then to assign the subcarriers to the various groups. In [START_REF] Zhang | An efficient resource-allocation scheme for spatial multiuser access in MIMO/OFDM systems[END_REF], the authors follow a completely different approach in which the users are first divided into groups such that the spatial cross-correlations among users in different groups is low as much as possible and then subcarriers are sequentially assigned within each group.

From the above discussion, it follows that the use of SDMA schemes in MIMO-OFDMA systems makes the problem of resource allocation more challenging as it requires the joint optimization of a) channel assignment and user partitioning; b) power allocation over all active links; c) transmit and receive filters. To the best of our knowledge, there exists only a few works dealing with all the above problems together. In [START_REF] Ho | Optimal resource allocation for multiuser MIMO-OFDM systems with user rate constraints[END_REF], the authors employ BD-ZF and Lagrange dual decomposition to derive a resource allocation scheme for minimizing the power consumption when individual user rate constraints are imposed. The main drawback of this approach is that an exhaustive search is required to find the best user allocation on each subchannel. A reduced complexity solution is illustrated in [START_REF] Hassan | Low complexity margin adaptive resource allocation in downlink MIMO-OFDMA system[END_REF], in which a two-step procedure is adopted to decouple BD-ZF beamforming from subcarrier and power allocation.

Although simpler than [START_REF] Ho | Optimal resource allocation for multiuser MIMO-OFDM systems with user rate constraints[END_REF], it still requires an exhaustive search over a subset of users. In [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF], the author exploits a layered architecture in which a user partitioning technique (resembling that discussed in [START_REF] Zhang | An efficient resource-allocation scheme for spatial multiuser access in MIMO/OFDM systems[END_REF]) is first used in conjunction with BD-ZF to partially remove multiuser interference and then carrier assignment is performed jointly with transceiver design using a linear programming (LP) formulation of the allocation problem [START_REF] Kim | Use of linear programming for dynamic subcarrier and bit allocation in multiuser OFDM[END_REF].

In this work, we return to the layered architecture investigated in [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF] and extend it in several directions. First, we reformulate the power minimization problem assuming that the quality-of-service (QoS) constraint of each user is given as a sum of the mean-square-errors (MSEs) over all subcarriers rather than on the sum of the achievable rates. Second, transceiver design is carried out employing a nonlinear THP precoder operating at user level at the transmitter. Third, the choice of the user partitioning strategy is motivated by its combination with the THP precoding technique. This allows us to completely remove the multiuser interference (rather than partially removing it) and to make use of a close-to-optimal partitioning strategy. All this leads to a resource allocation scheme of affordable complexity, which is shown by means of numerical results to outperform the solution presented in [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF].

II. PROBLEM DESCRIPTION

We consider 1 the downlink of an OFDMA network in which a total of N subcarriers is used to communicate with K MTs, each equipped with N R ≥ 2 antennas 2 . The BS is endowed with N T > N R transmit antennas. We denote by s n,k the N T -dimensional vector collecting the data transmitted to user k on subcarrier n and by a n,k ∈ {0, 1} the binary allocation variable, which is equal to one if subchannel n is assigned to user k and zero otherwise. The goal of this work is to minimize the total power consumption given by

P T = N n=1 K k=1 E s H n,k s n,k (1)
while satisfying user QoS requirements given as a function of the sum of the MSEs over all their assigned subcarriers. To be more specific, the expression for the kth user constraint is

N n=1 a n,k L ℓ=1 MSE n,k (ℓ) ≤ γ k ( 2
)
where L denotes the number of streams transmitted to the kth user over the nth subcarrier and MSE n,k (ℓ)

denotes its corresponding MSE. The quantities γ k > 0 are design parameters that specify different QoS requirements for each user. We assume that a maximum number of Q = ⌊N T /N R ⌋ users can be simultaneously allocated over each subcarrier, so that it is K k=1 a n,k ≤ Q for each channel n. To avoid the trivial solution where a user with no allocated subcarrier consumes no power and has a zero MSE, we require that at least n k subcarriers are assigned to each user so that it is N n=1 a n,k ≥ n k ∀k.

III. MULTI-USER INTERFERENCE ELIMINATION AND USER PARTITIONING

Unfortunately, solving the optimization problem described above requires an exhaustive search over all possible subcarrier allocations. Moreover, it needs also the joint optimization of the transmit and receive processing matrices for each allocation. All this makes its complexity extremely large for any practical scenario. To address this issue, we follow the approach of [START_REF] Zhang | An efficient resource-allocation scheme for spatial multiuser access in MIMO/OFDM systems[END_REF] and [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF], in which the population of K users is partitioned into Q different subsets {S (1) , S (2) , . . . , S (Q) }. This allows us to break the original 1 We use A = blkdiag {A1, A2, . . . , AK} to represent a block diagonal matrix whereas A -1 and tr {A} denote the inverse and trace of a square matrix A. We denote IK the identity matrix of order K while we use E {•} for expectation, • for the Euclidean norm of the enclosed vector and the superscript * , T and H for complex conjugation, transposition and Hermitian transposition. The notation [•] k,ℓ indicates the (k, ℓ)th entry of the enclosed matrix. 2 The results can be easily extended to a more versatile system in which a different number of services is required by each MT. In this case, K would simply denote the total number of services.

problem into a sequence of Q lower-complexity optimization sub-problems, each assigning all radio resources to a subset of users. Users within the same subset are transmitted on orthogonal subcarriers and do not interfere with each other. Channels allocation is performed sequentially starting from set S (1) .

From the above discussion, it follows that, after the Q allocation sub-problems are solved, there will be Q users assigned to each subcarrier. Without loss of generality, we focus on subcarrier n. Let us denote by K n the set of users assigned to n and by µ n (i) the user in S (i) associated to subcarrier n. To simplify the notation, in the following derivations the indexes µ n (i) will be relabelled according to the map µ n (i) → i. The signal x n,k ∈ C NR×1 received at the kth MT over the nth subcarrier can be thus written as

x n,k = H n,k Q i=1 s n,i + w n,k (3)
where w n,k ∈ C NR×1 is a Gaussian vector with zero-mean and covariance matrix σ 2 I NR and H n,k ∈ C NR×NT is the channel matrix over the nth subcarrier. From (3), it follows that the interference term is given by two different contributions, namely,

H n,k k-1 i=1 s n,i and H n,k Q i=k+1 s n,i .
The first term represents the interference caused by the active users already allocated before the kth assignment subproblem has been solved (i.e., users belonging to sets S (i) with indexes i < k), while the second term accounts for the users with indexes i > k (i.e., users which have been allocated after user k). In [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF], a BD-ZF scheme is employed to remove the first term while the second one is treated as Gaussian noise.

In the sequel, a THP technique operating at user level is used to remove both terms.

A. Multi-user interference elimination

The L ≤ ⌊N T /Q⌋ symbols transmitted to the kth user over the nth subcarrier are denoted by {d n,k (ℓ); ℓ = 1, 2, . . . , L}. They belong to an M -ary quadrature-amplitude modulation (QAM) alphabet with variance σ 2 d = 2(M -1)/3 and are stacked in the L-dimensional vector d n,k . As depicted in Fig. 1, the QL-dimensional data vector

d n = [d T n,1 , d T n,2 , . . . , d T n,Q ] T is pre-coded in a recursive fashion using a strictly block lower triangular matrix B n ∈ C QL×QL and a non-linear operator MOD M (•) that constrains the entries of b n,i ∈ C L×1 into the square region ℵ = {x (R) + jx (I) |x (R) , x (I) ∈ (- √ M , √ M ]}.
Denoting by

[B n ] i,ℓ ∈ C L×L the (i, ℓ)th block of B n , we have that b n,i ∈ C L×1 can be iteratively computed as [4] b n,i = d n,i - i-1 ℓ=1 [B n ] i,ℓ b n,ℓ + ς n,i i = 1, 2, . . . , Q (4)
where

[B n ] i,ℓ ∈ C L×L is the (i, ℓ)th block of B n , ς n,i is defined as ς n,i = 2 √ M ξ n,i and ξ n,i = [ξ n,i (1), ξ n,i (2)
, . . . , ξ n,i (L)] T with ξ n,i (ℓ) complex-valued quantity, whose real and imaginary parts are suitable integers that reduce b n,i (ℓ) to ℵ. The above equation indicates that the modulo operator is equivalent to adding a vector ς n to the input data d n . This produces the modified data vector

v n = d n + ς n = [v T n,1 , v T n,2 , , . . . , v T n,Q ] T from which b n is obtained as follows b n = C -1 n v n where C n ∈ C
LQ×QL is a block unit-diagonal and lower triangular matrix given by C n = B n + I LQ . The pre-coded vectors b n,i ∈ C L×1 are then linearly processed through the forward transmit matrices F n,i ∈ C NT ×L to produce s n,i = F n,i b n,i . The vectors s n,i for n = 1, 2, . . . , N and i = 1, 2, . . . , Q are finally fed to the OFDMA modulator and transmitted over the channel using the N T antennas of the BS array. As depicted in Fig. 2, at the MT the incoming waveforms are implicitly combined by the receive antennas and passed to an OFDMA demodulator whose outputs take the form in [START_REF] Tran | Multiuser diversity for successive zero-forcing dirty paper coding: Greedy scheduling algorithms and asymptotic performance analysis[END_REF]

with s n,i = F n,i b n,i . The complete elimination of H n,k Q i=k+1 F n,i b n,i at
F n,k = V (0) Hn,k U n,k (5)
where

U n,k ∈ C [NT -(k-1
)NR]×L is an arbitrary matrix and

V (0) Hn,k ∈ C NT ×[NT -(k-1
)NR] is a matrix whose columns form a basis for the null space of Hn,k obtained from its singular value decomposition (SVD). Setting F n,k as in ( 5) into (3) and stacking the received signals of all users into a single vector

x n = [x T n,1 x T n,2 • • • x T n,Q
] T , we may write

x n = T n b n + w n (6)
where

w n = [w T n,1 , w T n,2 , . . . , w T n,Q ] T and T n ∈ C NRQ×QL is a block lower triangular matrix with blocks [T n ] k,i ∈ C NR×L given by [T n ] k,i = H n,k V (0) Hn,i
U n,i for k ≥ i. We are now left with the problem of removing the interference term [START_REF] Tran | Multiuser diversity for successive zero-forcing dirty paper coding: Greedy scheduling algorithms and asymptotic performance analysis[END_REF]. To this end, we decompose T n in [START_REF] Zhou | Joint Tomlinson-Harashima precoding and scheduling for multiuser MIMO with imperfect feedback[END_REF] as

H n,k k-1 i=1 V (0) Hn,i U n,i b n,i in
T n = D n L n where D n = blkdiag{[T n ] 1,1 , [T n ] 2,2 , . . . , [T n ] Q,Q } and L n is a block unit-diagonal
and lower triangular matrix with

[L n ] k,i = [T n ] H k,k [T n ] k,k [T n ] H k,k -1 [T n ] k,i (7)
for

k > i. Substituting T n = D n L n into (6) and recalling that b n = C -1 n v n yields x n = D n L n C -1 n v n + w n from which setting C n = L n we obtain x n = D n v n + w n . Recalling that D n has a block-diagonal structure with blocks given by [T n ] k,k = H n,k V (0) Hn,k
U n,k , it follows that the multi-user MIMO system has been decoupled into |K n | parallel single-user MIMO links given by

x n,k = H ′ n,k U n,k v n,k + w n,k (8)
each of which represented by the equivalent channel transfer matrix

H ′ n,k = H n,k V (0) Hn,k
. This means that each user may operate in its corresponding link independently without affecting the other active users. Henceforth, we denote by

H ′ n,k = Ω H ′ n,k Λ 1/2 H ′ n,k V (1) H H ′ n,k
the SVD of H ′ n,k . As mentioned before, the vectors {x n,k } are processed by the kth mobile terminal for data recovery.

B. User partitioning

As mentioned above, MAI mitigation in SDMA-OFDMA systems is accomplished not only by precoding the users' data but also by partitioning the users and dynamically assigning the radio channels.

Unfortunately, optimal grouping is a problem of combinatorial complexity whose solution can only be found through an exhaustive search. To overcome this problem, a heuristic approach widely used in the literature is to partition users on the basis of their space cross-correlations (see for example [START_REF] Maciel | On the performance, complexity, and fairness of suboptimal resource allocation for multiuser MIMO-OFDMA systems[END_REF]).

Although reasonable, this approach has still a large complexity as it requires the calculation of the crosscorrelations among all users in the system over all available channels. Alternatively, in this work we exploit the fact that THP can be viewed as the transmit counterpart of the vertical Bell Labs layered space-time (V-BLAST) architecture and thus we order the users according to their channel qualities as as originally proposed in [START_REF] Wolniansky | V-BLAST: an architecture for realizing very high data rates over the rich-scattering wireless channel[END_REF] and later extended to THP in [START_REF] Kusume | Efficient Tomlinson-Harashima precoding for spatial multiplexing on flat MIMO channel[END_REF]. In our context, the channel quality of the kth user is measured by the following quantity:

π(k) = 1 N N n=1 tr H H n,k H n,k = 1 N N n=1 L ℓ=1 λ Hn,k (ℓ) (9)
where {λ Hn,k (ℓ)} denote the eigenvalues of H H n,k H n,k . The above quantities are used to partition users according to a worst-first criterion. In doing so, the users with the most attenuated channels are allocated in set S (1) whereas the users with the best channels are grouped in S (Q) . This choice is motivated by the fact that the null-space projection in (5) progressively reduces the available spatial diversity as the group index tends to Q and the number of rows of Hn,k increases up to (Q -1) N R . Therefore, since power consumption is in general dominated by users with the worst channel conditions, we give those users higher priority by placing them in set S (1) . Observe that the MAI arising among users (in different sets) allocated on the same subcarriers is mitigated jointly by THP and dynamic channel assignment. With the objective of minimizing the overall required power, channel assignment will automatically couple users that tend to not interfere with each other. It is worth observing that the same ordering strategy is used in [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF] following a different line of reasoning.

IV. LINEAR PROGRAMMING SUBCARRIER ASSIGNMENT

Without loss of generality, we focus on the resource allocation problem over the K/Q users within the set S (q) . For notational convenience, we denote by a (q) and U (q) the vector and the matrix obtained stacking the allocation variables and the precoding matrices of the users in S (q) , respectively. As before, the user indexes µ n (i) will be relabelled according to the map µ n (i) ← i. To make the problem mathematically tractable, we assume also that the precoded symbols b n,k are statistically independent and with the same power of user data3 , i.e., E{b n,k b H n,k } = σ 2 d I L . In these circumstances, using (5) it follows that the power required by the BS to transmit the signal s n,k is given by E{s

n,k s H n,k } = σ 2 d tr{U H n,k U n,k }.
The optimization problem can be thus mathematically formulated as: min

U (q) ,a (q) N n=1 k∈S (q) a n,k tr U H n,k U n,k (10)
subject to

N n=1 a n,k L ℓ=1 MSE n,k (ℓ) ≤ γ k k ∈ S (q) and N n=1 a n,k ≥ n k k ∈ S (q) (10.1)
which is a mixed-integer non-linear problem and thus not convex and very difficult to solve. A possible way out is to decouple the power allocation and subcarrier assignment problems. This can be achieved by assigning n k subcarriers to the kth user and designing the processing matrices such that the following constraint is satisfied

L ℓ=1 MSE n,k (ℓ) ≤ γ k n k . (11)
In this framework, the power is no longer an optimization variable but simply the cost of using n k subcarriers [START_REF] Moretti | A layered architecture for fair resource allocation in multicellular multicarrier systems[END_REF]. In particular, the cost c n,k of using subcarrier n for user k ∈ S (q) can be computed as

min Un,k tr U H n,k U n,k subject to L ℓ=1 MSE n,k (ℓ) ≤ γ k n k . ( 12
)
Once the solution of ( 12) is obtained, (10) can be recast as a linear integer programming (LIP) problem:

min a (q) N n=1 k∈S (q) a n,k c n,k (13)
subject to

N n=1 a n,k = n k k ∈ S (q) and k∈S (q) a n,k ≤ 1 ∀n
where the objective function and the constraints are linear in {a n,k }. In general, the solution of LIP problems can be found either performing an exhaustive search or relaxing the integrality condition on the allocation variable. In this particular case, the channel assignment in (13) has the advantage that can be modelled as a minimum cost flow problem and as such it is possible to show that the solution obtained by relaxing the integral condition is the optimal integral solution, so that very efficient solvers can be employed with no performance degradation [START_REF] Moretti | A layered architecture for fair resource allocation in multicellular multicarrier systems[END_REF].

A. Receiver design

To keep the complexity of the MTs at a tolerable level, we assume that a linear receiver is used for data recovery. As depicted in Fig. 2, vector x n,k in ( 8) is first processed by G n,k ∈ C L×NR to obtain

y n,k = G n,k H ′ n,k U n,k v n,k + G n,k w n,k (14)
and then passed to the same modulo operator employed at the transmitter so as to remove the effect of

ς n,k . The output z n,k = [z n,k (1), z n,k (2)
, . . . , z n,k (L)] T is finally fed to a threshold unit which delivers an estimate of d n,k . From ( 14), it follows that the received samples depend on G n,k and U n,k . The latter must be designed so as to mitigate co-channel interference while satisfying the QoS constraints. For this purpose, we adopt a ZF approach in which multi-stream interference is completely eliminated and the remaining degrees of freedom are exploited to minimize the power consumption under the constraint on the MSEs. The complete elimination of the multi-stream interference implies that

G n,k H ′ n,k U n,k = I L . ( 15
)
In these circumstances, the output z n,k (ℓ) from the modulo operator takes the form 4

z n,k (ℓ) = d n,k (ℓ) +
n n,k (ℓ) and its corresponding MSE results given by

MSE n,k (ℓ) = σ 2 [G n,k G H n,k ] ℓ,ℓ .
It can be shown that the optimal G n,k satisfying [START_REF] Wolniansky | V-BLAST: an architecture for realizing very high data rates over the rich-scattering wireless channel[END_REF] and minimizing each MSE n,k (ℓ) is the minimum norm solution of ( 15) [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. The latter is found to be

G n,k = (U H n,k H ′ H n,k H ′ n,k U n,k ) -1 U H n,k H ′ H n,k from which it follows that MSE n,k (ℓ) = σ 2 [(U H n,k H ′ H n,k H ′ n,k U n,k ) -1 ] ℓ,ℓ
. We now proceed with the design of the matrix U n,k , which requires to solve the following problem:

min {Un,k} tr U H n,k U n,k subject to L ℓ=1 σ 2 U H n,k H ′ H n,k H ′ n,k U n,k -1 ℓ,ℓ ≤ γ k n k . ( 16
)
The solution can be computed as follows.

Proposition 1: The optimal U n,k in ( 16) takes the form

U n,k = V (1) H ′ n,k Λ 1/2 Un,k S H n,k (17)
where

V (1) H ′ n,k is obtained from the SVD of H ′ n,k , Λ Un,k is diagonal and S n,k ∈ C L×L is a suitable unitary matrix such that MSE n,k (ℓ) = ǫ k for ℓ = 1, 2, . . . , L with ǫ k = 1 L γk nk .
In addition, the diagonal elements of Λ Un,k are given by

λ Un,k (ℓ) = ν n,k σ 2 λ H ′ n,k (ℓ) ℓ = 1, 2, . . . , L (18)
4 In writing z n,k (ℓ) = d n,k (ℓ) + n n,k (ℓ), we have neglected for simplicity the modulo-folding effect on the thermal noise.

Although not rigorous, this assumption is quite reasonable for moderate values of signal-to-noise ratios (see for example the book of Robert F. H. Fisher [START_REF] Fisher | Precoding and signal shaping for digital transmission[END_REF] for a complete treatment of the subject).

where

ν n,k is such that L ℓ=1 σ 2 λU n,k (ℓ)λ H ′ n,k (ℓ) = γk
nk . Proof : The proof is omitted for space limitations but it can be derived using the results illustrated in [START_REF] Sanguinetti | On the design of amplify-and-forward MIMO-OFDM relay systems with QoS requirements specified as Schur-convex functions of the MSEs[END_REF] since the sum of the MSEs is a Schur-convex function.

Using the results of Proposition 1, the cost c n,k in ( 13) is eventually given by

c n,k = L ℓ=1 λ Un,k (ℓ) (19)
with λ Un,k (ℓ) computed as in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF].

B. Complexity analysis

All the operations required by the proposed solution are summarized in Algorithm 1 whose computational load can be assessed in terms of the number of required floating point operations (flops) as follows 5 .

Observe that computing the quantities {π(k)} requires O(N KN T N R ) flops whereas computing the power cost c n,k according to [START_REF] Sanguinetti | On the design of amplify-and-forward MIMO-OFDM relay systems with QoS requirements specified as Schur-convex functions of the MSEs[END_REF] basically requires first to evaluate the SVDs of Hn,k for k = 1, 2, . . . , K/Q and n = 1, 2, . . . , N and then those of H ′ n,k in (8) for k = 1, 2, . . . , K and n = 1, 2, . . . , N . The total number of flops required for these two operations are summarized in the second and third row of Table I. In writing these figures, we have taken into account that evaluating the SVDs of

H ′ n,k requires O(Q/2(Q -1)N T N 2 R + QN R N 2 T ) flops in total since O(QN R N 2 T ) flops are needed to compute H ′ n,k = H n,k V (0) Hn,k whereas O(Q/2(Q -1)N T N 2 R
) flops are required for the SVD. Summing all the above terms it turns out that the overall complexity for computing all costs {c n,k } is approximately given by O(N KQN R N 2 T ). The complexity of solving ( 13) is an open research issue. The latest results (see for example [START_REF] Huang | Fast b-matching via sufficient selection belief propagation[END_REF] references therein) place the complexity of the assignment problem in a range between O κ 2 and O κ 2.5 with κ being the total number of nodes. In our case, the number of nodes is the sum of the number of users per single allocation problem plus the number of subcarriers, i.e., κ = N +K/Q. Since we have Q distinct subproblems to solve, the overall complexity of the LP optimization is approximately

given by O Q(N + K/Q) 2.5 flops. The computation of B n = C n -I LQ in (4) with C n = L n can be assessed as follows. Evaluating each [L n ] k,i in (7) requires O((N 3 R +4LN 2 R )) flops. Since the total number of matrices [L n ] k,i is Q/2(Q -1), it follows that O(N Q/2(Q -1)(N 3 R + 4LN 2 R
)) flops are required to obtain all matrices {C n } and thus all {B n }. The computational load for obtaining {F n,k }, {G n,k } and {U n,k } can be reasonably neglected as it basically require to put together all the unitary matrices computed above with SVDs. The processing requirements of the proposed two-layer architecture are 5 In doing so, we make use of the following results: i) the multiplication of A ∈ C m×n and B ∈ C n×p requires O(2mnp) flops; ii) evaluating the SVD of A ∈ C m×n needs O(mn 2 ) flops; iii) the inverse of A ∈ C n×n requires O(n 3 ) flops.

summarized in Table I from which it follows that the overall number of flops is approximately given by

O Q(N + K/Q) 2.5 + N KQN R N 2 T + N Q 2 N 3 R .
The latter is comparable to the computational load required by the scheme illustrated in [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF] as it is dominated by the computational burden required by the LP approach, especially when the number of subcarriers relatively large. However, as shown in the sequel, the proposed solution provides much better performance in terms of power reduction with respect to [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF] thanks to the underlying THP scheme.

V. NUMERICAL RESULTS

We consider a system with K uniformly distributed users in a cell of radius R = 100 m. The propagation channel is static, frequency-selective and modelled as a Rayleigh fading process with an exponentially decaying power delay profile. The path loss exponent is β = 4. Unless noted differently, the number of users is K = 16.

We compare the proposed architecture, denoted by THP Tx -Lin Rx, with three other algorithms: a) a ZF linear beam-former, denoted as ZF Tx, b) a THP scheme, denoted as THP Tx (see for example [START_REF] Zhou | Joint Tomlinson-Harashima precoding and scheduling for multiuser MIMO with imperfect feedback[END_REF]), and c) the architecture proposed in [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF] that employs linear processing at both the transmitter and the receiver (Lin Tx -Lin Rx). In details, letting

H n = [H T n,1 H T n,2 • • • H T n,Q ] T and F n = [F n,1 F n,2 • • • F n,Q ] T , the precoding matrix for ZF Tx is F n = H H n (H n H H n ) -1 .
The THP Tx architecture is realized by setting F n = Q n and C n = R -H n with Q n and R n being computed as the QR decomposition H H n , i.e., H H n = Q n R n . Both ZF Tx and THP Tx schemes are designed to remove the inter-stream and inter-user interference at the transmitter so that the receive filter is G n,k = I L .

We consider three different scenarios, summarised in Table II, which are designed to observe the behaviour of the proposed algorithms when the total number of available channels per user is fixed and frequency channels are progressively replaced by streams in the spatial domain. More in details, the first scenario, referred to as S (1) , is a 2 × 1 MIMO system with a bandwidth W (1) = 10 MHz and N (1) = 64 orthogonal subchannels. The bandwidth of Scenario S (2) is W (2) = 5 MHz, spanning N (2) = 32 subchannels with a 4 × 2 MIMO configuration. Scenario S (3) transmits over a bandwidth W (3) = 2.5 MHz with N (3) = 16 subchannels and employs a 8 × 4 configuration. For each scenario we assume that the number of allocated subcarriers is n

(i) k = N (i) × Q/K
and the total number of channels per user is n (i) k L (i) = 8 (i = 1, . . . , 3; k = 1, . . . , K) regardless of the scenario considered. Figs. 3 -5 report the total transmit power for the three scenarios as a function of the average target MSE ρ per data stream. By design, for a given value of ρ, the overall MSE is γ (i) k = 8ρ (i = 1, . . . , 3; k = 1, . . . , K). Results show that the gains obtained thanks to the implementation of non-linear processing progressively increase from scenario S (1) to S (3) , as the spatial dimension becomes more important.

In particular, Fig. 3 shows that, with a 2 × 1 configuration and 64 channels, all the schemes, except ZF Tx, tend to have similar performance. The effect of resource allocation is predominant and the users transmitting on the same channel are sufficiently separated regardless of the specific architecture.

As the number of orthogonal frequency channels is reduced, the consequent diminution in frequency diversity is only partially compensated by the larger number of antennas: in facts, even if the total number of channels is the same, the spatial streams tend to be more correlated. In this case, the choice of the transceiver architecture plays a very important role since channel allocation alone is not able to fully exploit all the diversity of the the system. The results plotted in Fig. 4 show that the THP-based schemes largely outperform all other solutions.

The same trend appears in Fig. 5, where THP Tx -Lin Rx effectively exploits the spatial diversity provided by the multiple antennas. Scenario S (1) requires less power when compared to S (2) and S (3) as it occupies a larger bandwidth. In scenarios S (2) and S (3) , the proposed scheme takes advantage of the increased spatial dimension to transmit the same amount of data employing a comparable amount of power and occupying only a fraction of the bandwidth. Fig. 6 shows the total transmit power for an average target MSE ρ = 0.25 as a function of K for S (1) and S (3) . For ease of representation, only the results of THP Tx, Lin Tx -Lin Rx and THP Tx -Lin Rx are reported. As before, the parameters are set so that the number of data stream per user is the same (regardless of the specific scenario). An accurate inspection of the results shows that for scenario S (1) , the performance of the three algorithms tend to be very close for K ≥ 16, when the resource allocation algorithm is able to fully exploit both multi-user and frequency diversity. The situation is remarkably different for scenario S (3) where it appears that resource allocation alone is not sufficient to completely deal with MAI. In fact, all multiuser diversity is already exploited for K = 8 and further increase of the number of users produce only marginal improvements. In this case, the THP Tx -LIN Rx configuration outperforms the other two schemes thanks to its capability to cancel the MAI.

VI. CONCLUSIONS

We have derived a resource allocation scheme for the downlink of SDMA-MIMO-OFDMA systems.

The proposed solution relies on a layered architecture in which MAI is first removed by means of a THP technique operating at user level and then channel assignment and transceiver design are jointly addressed using a ZF-based linear programming approach that aims at minimizing the power consumption while satisfying specific QoS requirements given as the sum of the MSEs over the assigned subcarriers. The proposed approach outperforms the existing solutions, especially when the frequency diversity is small and the number of spatial modes is large. tr H H n,k H n,k . 3: end for 4: Sort users according to π(k) and group them in Q sets {S (1) , . . . , S (Q) }. for subcarrier n = 1 to N do 8:

Compute the power cost c n,k according to [START_REF] Sanguinetti | On the design of amplify-and-forward MIMO-OFDM relay systems with QoS requirements specified as Schur-convex functions of the MSEs[END_REF] Solve the resource allocation problem in [START_REF] Moretti | Efficient margin adaptive scheduling for MIMO-OFDMA systems[END_REF].

12:

for subcarrier n = 1 to N do ZF Tx THP Tx Lin Tx -Lin Rx THP Tx -Lin Rx Fig. 3. Total power consumption for the investigated solutions in scenario S (1) for different MSEs for data stream. ZF Tx THP Tx Lin Tx -Lin Rx THP Tx -Lin Rx Fig. 5. Total power consumption for the investigated solutions in scenario S (3) for different MSEs for data stream.

Fig. 1 .Fig. 2 .Algorithm 1

 121 Fig. 1. Block diagram of the THP technique operating at user level.

5 : 6 :

 56 for group i = 1 to Q do for user k = 1 to |S (i) | do 7:

  .

13 :

 13 Compute B n = C n -I LQ . for user k = 1 to |S (i) | do 16:for subcarrier n = 1 to N do 17:Compute {F n,k , G n,k , U n,k }.

Fig. 4 .

 4 Fig.[START_REF] Fisher | Precoding and signal shaping for digital transmission[END_REF]. Total power consumption for the investigated solutions in scenario S(2) for different MSEs for data stream.

Fig. 6 .

 6 Fig.6. Total power consumption for the investigated solutions for different number of users.

  the transmitter can be achieved by constraining F n,k to lie in the null space of Hn,k = [H T

n,1 , H T n,2 , . . . , H T n,k-1 ] T . Accordingly, this amounts to letting F n,k have the following structure

Although not rigorously true, this assumption is reasonable for large M -QAM constellations with size M ≥ 16[START_REF] Fisher | Precoding and signal shaping for digital transmission[END_REF].


Jean Peyhardi
email: jean.peyhardi@gmail.com

New qualitative choice models incorporating individual and choice characteristics

Keywords: Qualitative choices, Conditional logit model, Link function, Design matrix

In the econometric framework, multinomial and conditional logit models are the most usual regression models for qualitative choices. They differ by their parametrization while sharing the canonical link function. This link function can be decomposed into the reference ratio of probabilities and the logistic cumulative distribution function (cdf). We propose to conserve the reference ratio, appropriate for qualitative choices, but to select the cdf among an enlarged family containing the Student cdf for instance. These new qualitative choice models often outperform logit models in terms of likelihood and error rate of classification and stay easily interpretable. This is illustrated with a benchmark dataset of travel demand between Sydney and Melbourn.

Multinomial and conditional logit models for qualitative choices

Let Y i be the response variable corresponding to the choice of individual i (with alternatives j = 1, . . . , J) and x i be the vector of individual characteristics (e.g. sex, age). In the context of travel demand, some choice characteristics ω i,j are also used, such as the cost of alternative j for individual i. In the following we will suppress the individual subscript i without loss of generality.

Luce's choice axiom [START_REF] Luce | Individual Choice Behavior: A Theoretical Analysis[END_REF] and the principle of random utility maximisation lead to the logit model defined by

π j = exp(η j ) 1 + J-1 k=1 exp(η k )
This paper was published as a part of the proceedings of the 30th International Workshop on Statistical Modelling, Johannes Kepler Universität Linz, 6-10 July 2015. The copyright remains with the author(s). Permission to reproduce or extract any parts of this abstract should be requested from the author(s).

for j = 1, . . . , J -1, where π j = P (Y = j). Depending on the form of the linear predictors η j , we obtain different logit models:

• η j = α j + x T δ j . Individual characteristics x are used with J -1 different slopes δ j . This is the classical multinomial logit model.

• η j = α j + ωT j γ where ωj = ω jω J . Choice characteristics ω j are used with common slope δ. This is the conditional logit model introduced by [START_REF] Mcfadden | Conditional logit analysis of qualitative choice analysis[END_REF].

• η j = α j + x T δ j + ωT j γ. Individual and choice characteristics are used with respectively different slopes and common slope. This is a combination of the two previous parametrizations.

Generalisation of multinomial and conditional logit models

All the classical regression models for categorical data [START_REF] Tutz | Regression for Categorical Data[END_REF] share the generic equations [START_REF] Peyhardi | A new specification of generalized linear models[END_REF])

r j (π) = F (η j ) for j = 1, . . . , J -1, where r is a C 1 -diffeomorphism from the simplex ∆ = {π ∈ (0, 1) J-1 | J-1 j=1 π j < 1} (corner of hypercube)
to an open subset of the hypercube (0, 1) J-1 , π is the vector of probabilities (π 1 , . . . , π J-1 ) T and F is a continuous and strictly increasing cdf. Let us remark that the three logit models (defined in Section 1) share the canonical link function, specified by the reference ratio r j (π) = π j π j + π J and the logistic cdf

F (η) = exp(η) 1 + exp(η)
.

Therefore, these three logit models are specified by the (reference, logistic, Z) triplet with different design matrices Z respectively equal to

Z 1 =    1 x T . . . . . . 1 x T    , Z 2 =    1 ωT 1 . . . . . . 1 ωT J-1    , Z 3 =    1 x T ωT 1 . . . . . . . . . 1 x T ωT J-1    .
The reference ratio is mandatory for non-ordered choices whereas the logistic cdf is not [START_REF] Peyhardi | A new specification of generalized linear models[END_REF]. We thus propose a new class of regression models appropriate for qualitative choices defined by (reference, F , Z i ) models (i = 1, 2, 3) where the cdf F can be selected among e.g. the logistic, Gaussian, Laplace, Gumbel, Gompertz, and Student cdfs (with different degrees of freedom ν ∈ R * + ). The heavy tails of Student distributions may markedly improve the model fit and reduce the classification error rates [START_REF] Peyhardi | A new specification of generalized linear models[END_REF]. Parameter estimates stay interpretable since

π j π J = F 1 -F (η j )
is strictly increasing with η j (we have π j /π J = exp(η j ) in the case of the logistic cdf). Finally, this family of reference models for qualitative choices is easily estimated using the standard Fisher's scoring algorithm.

Fisher's scoring algorithm for reference models

Let us remark that the Fisher's scoring algorithm is simplified in the particular case of the reference ratio (compared to the adjacent, cumulative and sequential ratios) which is a part of the canonical link function. Using the chain rule we obtain the score

∂l ∂β = Z T D(y -π),
and the Fisher's information matrix

E ∂ 2 l ∂β T ∂β = -Z T DCov(Y )DZ,
where

D = diag 1≤j≤J-1 f (η j ) F (η j ){1 -F (η j )} ,
and f is the density function. Remarking that f = F (1 -F ) for the logistic distribution, the Fisher's scoring algorithm turns out to be, in this particular case, the algorithm for the canonical link function.

Application to travel mode demand

The dataset, used by [START_REF] Greene | Econometric Analysis[END_REF], contains 210 observations of choice among J = 4 travel modes between Sydney and Melbourn (Australia): air (1), bus (2), train (3), and car (4). The two individual characteristics are the household income x 1 and the number of people travelling x 2 . The three choice characteristics are the terminal time ω 1 j (ω 1 4 = 0 for car), the amount of time spent traveling ω 2 j and the in-vehicle cost ω 3 j . The sample is choices based so as to balance it among the four choices knowing that the true population is dominated by drivers.

The three logit models and other reference models (i.e. F = logistic) were estimated. Best results were obtained with (reference, Student ν=1 , Z i ) models that markedly outperformed logit models; see Table 1 (we have for instance l = -192.89 for logistic versus l = -169.79 for Student with the same parametrization Z 2 ). The (reference, Student ν=1 , Z 2 ) is the best model according to BIC. The proportions between parameters when significant are approximatively conserved comparing logistic and Student models (α 1 /α 2 1.43 for logistic cdf and α 1 /α 2 1.72 for Student cdf for instance). The interesting difference concerns estimate of the slope γ 1 since γ 1 /γ 2 24 for logistic cdf and γ 1 /γ 2 60 for Student cdf. In the Student case, the terminal time has a stronger impact on the travel mode choice. Finally, the selection of F does not a priori change the sign of parameters but may change the proportion between them. Moreover it may reduce the classification error rate and thus increase the precision of predictions.

TABLE 1 .

 1 Parameters estimates for six reference regression models.

			F = logistic		F = Studentν=1
		Z1	Z2	Z3	Z1	Z2	Z3
	α1	0.9435	4.7399	6.0351	0.6438	13.7305	15.2387
	α2	1.978	3.3062	4.5045	1.9446	7.955	7.3668
	α3	2.4938	3.9532	5.5735	2.354	8.6827	9.5013
	δ 1 1	0.003544		0.007481	0.00496		0.02897
	δ 1 2	-0.03033		-0.0209	-0.02581		-0.004521
	δ 1 3	-0.05731		-0.05923	-0.06026		-0.0571
	δ 2 1	-0.6006		-0.9224	-0.4946		-1.0745
	δ 2 2	-0.9404		-0.1478	-1.0836		0.7765
	δ 2 3	-0.3098		0.2163	-0.2489		0.7777
	γ 1		-0.09689	-0.1012		-0.2548	-0.2597
	γ 2		-0.003995 -0.004131		-0.00426 -0.003878
	γ 3		-0.01391	-0.008667		-0.01849	-0.01746
	L	-253.34	-192.89	-172.47	-253.36	-169.79	-159.02
	BIC	554.8	417.86	409.1	554.84	371.66	382.21
	Error	53.33%	26.19 %	27.14 %	54.29 %	22.38%	21.9%


Transient UWB Antenna Near-Field and Far-Field Assessment From Time Domain Planar Near-Field Characterization: Simulation and Measurement Investigations

I. INTRODUCTION

OR Ground Penetrating Radar applications we are interested in measuring the free space transient response of Ultra Wide Band (UWB) antennas composing the radar [START_REF] Daniels | Surface Penetrating Radar[END_REF] [START_REF] De Jongh | Design and analysis of new GPR antenna concepts[END_REF]. The duration of the transient response is of a crucial importance and the designers try to shorten the antenna free space impulse response while enlarging its frequency bandwidth of matched input impedance (S11<-10dB) [START_REF] Turk | Ultra-wideband TEM horn design antenna design for ground-penetrating impulse radar systems[END_REF]. More generally, GPR using UWB non-dispersive antennas provides easily interpretable radargrams [START_REF] Yarovoy | The dielectric wedge antenna[END_REF] [START_REF] Yarovoy | Reconstruction of the Field Radiated by GPR Antenna into Ground[END_REF].

UWB antennas transient characterization can be carried out using NF or FF techniques [START_REF] Hansen | Plane-wave theory of time-domain fields, near-field scanning applications[END_REF] [START_REF] De Jough | Antenna time-domain measurement techniques[END_REF]. Based on the Huygens principle, the NF measurement techniques make use of the tangential components of electric or magnetic field collected over a scan surface in the vicinity of the AUT. Then, these measured data are transformed to calculate the electromagnetic field at different distances [START_REF] Yaghjian | An overview of near-field antenna measurements[END_REF]. The NF measurements can be conducted in a small and low cost Manuscript received February 19, 2015. M. Serhir is with Laboratoire Genie electrique et electronique de Paris (GeePs) (UMR 8507 CentraleSupelec -CNRS -UPSud -UPMC) 11 rue Joliot-Curie -Plateau de Moulon 91192 Gif sur Yvette Cedex -France, (email: mohammed.serhir@centralesupelec.fr).

anechoic chamber and by time-gating the measured data, one can filter out multiple reflections occurring during the measurement [START_REF] Blech | A time domain spherical near-field measurement facility for UWB antennas employing a hardware gating technique[END_REF] [START_REF] Serhir | Development of pulsed antennas measurement facility: Near field antennas measurement in time domain[END_REF]. This allows the data accuracy enhancement especially for low frequency measurement in small anechoic chambers for which the absorption quality is unsatisfactory. The TD techniques permit also the radiation pattern measurement in non-anechoic environment as presented in [START_REF] Jinhwan | Space Radiation Pattern Reconstruction from Non-Anechoic Measurements Using an Impulse Response of the Environment[END_REF][START_REF] Rammal | Far-Field Reconstruction from Transient Near-Field Measurement Using Cylindrical Modal Development[END_REF][START_REF] Andrieu | Transient characterization of a novel ultrawide-band antenna: the scissors antenna[END_REF][START_REF] Levitas | UWB system for time-domain near-field antenna measurement[END_REF]. They are well adapted for studying the parasitic electromagnetic radiation emitted from electronic devices for electromagnetic compatibility purposes [15].

Even the theoretical development of the NFNF and NFFF transformations are detailed in [START_REF] Hansen | Plane-wave theory of time-domain fields, near-field scanning applications[END_REF][16] [START_REF] Hansen | Planar near-field scanning in the time domain .2. Sampling theorems and computation schemes[END_REF], the experimental utilization of TD NF techniques are less popular than the FD techniques. This is due to difficult experimental implementation combined with the drawback related to the computer time requirement, which could make unrealistic the application of the technique to cases of practical interest. The transient NF measurement systems use wideband and small probes with low interaction with the AUT [START_REF] Yarovoy | Ultrawideband sensor for Electromagnetic Field Measurement in Time Domain[END_REF]. Moreover, in NFNF and NFFF transformation schemes the measured Efield and its time derivative are involved in the calculation process. The accuracy of the NFNF and NFFF transformation results depends on the highest point of the chosen interpolation technique to keep working with the minimum number of measured data [START_REF] Bucci | Near-field far-field transformation in time domain from optimal plane-polar samples[END_REF]. In this paper, we use the reconstruction formula to calculate and interpolate the E-field and its time derivative as well for NFNF and NFFF transformation schemes.

Here we aim at validating NFNF and NFFF transformations using the experimental measurement and the electromagnetic simulation software CST MWS [START_REF]CST Microwave studio 2014 user's guide[END_REF]. An UWB Vivaldi antenna designed for Ground Penetrating Radar (GPR) is used for this purpose. The Vivaldi TD NF data are collected at different distances from the antenna. First, the NFFF transformation results are presented in order to outline the errors due to the measurement surface truncation. Second, we evaluate the accuracy of the NFNF calculation using different set of measured NF data collected at different distances from the AUT. Third, we present the experimental validation of the NFNF transformation scheme we have developed based on the Green's function representation. Our algorithms are applied to data obtained from the measurement of the Vivaldi antenna in The excitation signal spectrum (right), the maximum of the signal spectrum is at 1.75 GHz Supelec planar near-field measurement facility. Finally, the experimental validation aims also at comparing the simulation and measurement NF results. This paper is structured as follows. In Section II, the planar NFNF and NFFF transformations formulations are described.

In Section III, we present comparisons of FF and NF transformation results with the actual fields in the time domain using the simulation software CST MWS. Then, the measurement facilities used for our validation purpose are described and the NFNF transformation routine is validated using planar and cylindrical experimental setups in Section IV. Finally, concluding remarks are presented summarizing the relevance of our contribution in Section V. All theoretical formulations are expressed in the S.I. rationalized system with 𝑒 -𝑖𝜔𝑡 time dependence.

II. PLANAR NEAR-FIELD TRANSFORMATION FOR FORWARD PROPAGATION

Let us consider UWB Vivaldi antenna. This is used to validate NFNF and NFFF transformation algorithms. The antenna has been optimized for the frequency band (0.5 GHz -3 GHz) and its S11 is presented in Fig. 1. In this paper we have adopted the NFNF and NFFF transformation algorithms based on the Green's function representation where NF and FF are determined directly in the time-domain. This choice has been made to clearly distinguish our work from the NFNF or NFFF transformations based on modal expansions (plane wave expansion in the frequency domain) [15].

A. Planar Near-Field to Near-Field transformation

Using the NF tangential components measured over the planar surface at the distance xmeas from the AUT the NFNF transformation algorithm allows calculating the transient field at different distances in the half space x ≥ xmeas. The NFNF transformation is expressed as [6, eq. ( 5

E x r t R c cR t E r t R dydz x x x E r t R c R E x r t R c cR t R y z x E r t R c R                                                          Ny Ny   (1)
The calculated field ( , ) E r t in the left side of ( As presented in Fig. 2

  0 1 0 0 0 0 0 1 0 0 0 0 0 0 sin ( , ) ( , ) ( ).
                                                              (2)
This means that in the time window

  0 0 1 t t t t N t     
and provided the NF is sampled respecting the Nyquist criterion / max t    , the time derivative of the E-field is expressed as

  1 0 0 0 0 0 0 () ( , ) ( , ) 1 with ( ) cos ( ) () t N mn mn l E h t r t E r t l t t t t t l t h t h t t t t l t t                                   (3)
In our analysis we use ( 2) and (3) to calculate the E-field and its time derivative between time-samples (at t R c  ) provided the sampling Nyquist criterion is respected. This criterion, which is linked to spatial NF sampling Δy0 = Δz0 = cΔt, is sufficient to accurately perform the NFNF transformation.

B. Planar Near-Field to Far-Field transformation

The time domain NFFF transformation is formulated as [6, eq. ( 5

                               (4)
The antenna radiation pattern at the far-field distance . As expressed in (4), the infinite integration is truncated and we implicitly suppose that the field is weak beyond the measurement surface edges. However, one of the advantages of time-domain near-field measurements is that one can obtain exact early-time far-field pattern values even though the field outside the scanned area is strong.

In the next sections, we present the validation results of the NFNF and NFFF transformation using simulation and measurement results and the E-field time derivative calculated using (3).

III. ELECTROMAGNETIC SIMULATION RESULTS

We have developed a routine expressing (4) to study the NFFF transformation results when using different NF data collected at different measurement planes xmeas .The Vivaldi antenna has been simulated using the transient simulator of the CST MWS to provide the NF data. The AUT is excited using a cosinusoidal voltage modulated by Gaussian signal to cover the AUT frequency band (0.5 GHz -3.5 GHz),

    22 0 4( ) 0 0 cos ( ) tt e t t t e       with 00 2f
  , f0 = 1.75 GHz, τ = 0.8.10 -9 s and the delay t0=1.75 τ. The excitation signal and its spectrum normalized at f0 = 1.75 GHz are presented in Fig.

3.

The antenna TD NF is collected over the plane x = xmeas respecting regularly spaced positions Δy0 and Δz0. Based on [START_REF] Hansen | Plane-wave theory of time-domain fields, near-field scanning applications[END_REF] the bandwidth of the Gaussian signal is approximately equal to 12/2πτ. One can sample without any loss of information the radiated field respecting the Nyquist criterion Δy0 = Δz0 = λmin/2, where λmin = c/fmax, and fmax is the maximum frequency expressed as fmax = f0+12/2πτ = 4.14 GHz.

Using the CST MWS we generate the antenna transient radiated NF over the plane xmeas = 30 cm. This scan surface is delimited by Dzmax = 100 cm and Dymax = 75 cm as defined in Fig. 2 and the NF in the cut planes ymeas = 0 and zmeas = 0 are presented in Fig. 4 for Ey and Ez components. The Vivaldi antenna is linearly polarized and the E-field is predominantly radiated through the Ey component and the Ez component is negligible comparatively to Ey.

Since the S11 (Fig. 1) of the Vivaldi antenna demonstrates a pass-band behavior (0.5 GHz -3 GHz), we are interested in evaluating the accuracy of the NFFF transformation as a function of the NF sampling frequency. To do so, we compare the actual FF obtained from the CST MWS with the NFFF transformation results using 3 different NF sampling frequencies f1 = 3 GHz, f2 = 3.5 GHz, fmax = 4.14 GHz. In Figs. [START_REF] Yarovoy | Reconstruction of the Field Radiated by GPR Antenna into Ground[END_REF] we present the FF comparison at the observation directions (θ = π/2, φ = 0), (θ = π/4, φ = 0) and (θ = π/2, φ = π/4) for the FF components Etheta and Ephi. As seen in Figs. 5 good agreements are noticed between the Ephi (co-polar) component of the actual FF and the ones resulting from NFFF transformations. The use of the sampling frequency f1 = 3 GHz is responsible of the aliasing effects clearly seen for the cross-polarization component Etheta. In addition, the differences between the curves of Fig. 5 are due to measurement surface truncation. For the far-field observation point (θ = π/2, φ = 0) the difference between the actual FF and the calculated one starts at t = 5.11τ that corresponds to the arrival time when the center of the pulse reached the measurement edge Dymax = 75 cm. The error corresponding to the measurement edge Dzmax = 100 cm is visible around t = 6.10τ.

Using the sampling frequency fmax = 4.14 GHz, we compare the NFFF results calculated from NF data collected at different planes xmeas = 20 cm, 30 cm, 40 cm, 50 cm. These calculated FF are compared in Figs. 6 with the actual one at the directions (θ = π/2, φ = π/6) and (θ = π/3, φ = 0). As presented in Figs. 6, good agreements are noticed for 2.5τ ≤ t ≤ ttr, where ttr depends on the measurement distance xmeas and the FF observation point (θ, φ). This effect is known as the truncation error and the FF reliable region is defined by αθ and αφ expressed in Fig. 2 depends on the measurement distance.

Based on the presented simulation results, the developed NFFF transformation routine has been validated using the CST MWS. The NFFF transformation results have shown a satisfactory accuracy while using the Nyquist criterion for NF sampling.

The next comparisons aim at validating the TD NFNF transformation routine. For this, we provide tangential NF data at xmeas = 20 cm for -75 cm ≤ ymeas ≤ 75 cm and -100 cm ≤ zmeas ≤ 100 cm sampled using fmax = 4.14 GHz. Thereafter, we perform the NFNF transformation routine to calculate the field at 50 cm from the AUT. The NFNF transformation results (ENFtoNF) are compared with the actual NF (Eref) at the plane cuts (x = 50 cm, y = 0, -100 cm ≤ z ≤ 100 cm) and (x = 50 cm, -75 cm ≤ y ≤ 75 cm, z = 0) using the error expressed as:

            0 0 max max max max
             . ( 8
)
The error values are presented in Figs. [START_REF] De Jough | Antenna time-domain measurement techniques[END_REF] for the Ey and Ez components. As it is seen, the NFNF transformation errors are very low for the Ey components (co-polar) and reaches 5% at the edges of the measurement surface (y = ±75 cm and z = ±100 cm) for Ez component (cross-polar). These errors are due to NF measurement surface truncation. In Figs. 8, the calculated AUT transient radiated fields Ey and Ez are compared with the actual ones at different observation points A(x = 50 cm, y = 0, z = 0), B(50 cm, 25 cm, 0), C(50 cm, 0, 70 cm), D(50 cm, 50 cm, 0), E(50 cm, 0, 98 cm) and F(50 cm, 72 cm, 0). The calculated NF agrees very well with the actual one even for the observation points E and F situated in the area of the high error values (maximum 5%) of Figs. 7.

To get a deep insight into the origin of the NFNF transformation errors we calculate the field at x = 50 cm from the AUT using NF data collected over the scan planes xmeas = 10 cm and 30 cm. These are compared with the actual NF (the CST MWS) at the plane cuts (x = 50 cm, y = 0, -100 cm ≤ z ≤ 100 cm) and (x = 50 cm, -75 cm ≤ y ≤ 75 cm, z = 0). The resulting error values are presented in Figs. 9 (a) and (b). From Figs. [START_REF] De Jough | Antenna time-domain measurement techniques[END_REF] and Figs .9 we can conclude that the calculated results are completely unreliable outside of a certain spectral region. The erroneous ripples are due to the discontinuity of the measured field at the edge of the truncated surface. Hence, the Fig. 8. The actual field comparison with transient Ey(V/m) and Ez(V/m) resulting from the NFNF transformation. The NF data is measured at xmeas=20cm and transformed to calculate the antenna transient response at the spatial positions A(x = 50cm, y = 0, z = 0), B(50 cm, 25 cm, 0), C(50 cm, 0, 70 cm), D(50 cm, 50 cm, 0), E(50 cm, 0, 98 cm) and F(50 cm, 72 cm, 0). entire calculated pattern is always affected by errors, and it is not possible to define a region where the error is completely zero. However, the concept of the spectral reliable far-field region is usually applied to refer to the region in which the error is not negligible but is low [START_REF] Martini | Reduction of Truncation Errors in Planar Near-Field Aperture Antenna Measurements Using the Gerchberg-Papoulis Algorithm[END_REF].

From Figs. 7 and Figs. 9, the highest error value over the principal component Ey stay stable around 3%. In Figs. 7 and Figs. 9 (a) corresponding to NF data collected at xmeas = 10 cm and 20 cm, the truncation error is higher in the region -50 cm ≤ y ≤ 50 cm and -75 cm ≤ z ≤ 75 cm in comparison with the results of xmeas = 30 cm. The truncation error is directly linked to the level of the discontinuity of the measured field at the edge of the truncated surface. At xmeas = 10 cm and 20 cm this discontinuity is more important than for xmeas = 30 cm.

Up to now, we have studied the accuracy of NFFF and NFNF transformation routines using electromagnetic simulation software (the CST MWS). In the next section we describe the experimental validation of the developed transformation routines.

IV. MEASUREMENT RESULTS

The Vivaldi antenna used in the simulation presented above has been fabricated and Supelec planar NF range to measure its radiated field is presented in Fig. 10. The Vivaldi antenna and the measurement probe are placed inside the anechoic chamber of dimensions 4 m x 5 m x 3 m. We measure the tangential near-field components Ey and Ez while moving the probe respecting a regular spacing Δy = Δz = 3 cm between measurement points over the planar surface (-60 cm ≤ ymeas ≤ 60 cm and -60 cm ≤ zmeas ≤ 60 cm) at xmeas = 10 cm in front of the antenna. For our measurement investigations we have used the wideband probe "En-Probe EFS-105" [22]. This probe is composed of a small dipole (dimensions 6.6 mm x 6.6 mm) connected to the base unit through a fiber optic link and the base unit is connected to the measuring VNA through a short cable. This probe is dedicated to NF measurement covering the frequency band (5 MHz -3 GHz). The probe performances can be extended to 3.5 GHz with slight measurement sensitivity degradation in the frequency band (3 GHz -3.5 GHz).

The antenna NF measurement is carried out in the frequency band (0.3 GHz -3.5 GHz) using a Vector Network Analyzer (Agilent ENA 5071b). Thereafter the measured frequency-domain NF complex data (amplitude and phase) are Fourier transformed to set the TD NF of the Vivaldi antenna when excited by the pulse given in Fig. 2. In Fig. 10, we compare the measured TD NF at the point (xmeas = 10 cm, ymeas = 0, zmeas = 0) with the one obtained from the CST MWS. The Ey component (co-polar) fits well the CST results. However, the Ez component (cross-polar) shows some discrepancies. The NF probe used in this measurement campaign is supposed to behave like a point source. For probe calibration procedure, we have normalized the measured NF magnitude (co-polar) in order to reach the CST MWS maximum NF level. This is performed at the principal direction (xmeas = 10 cm, ymeas = zmeas = 0) and the normalization coefficient is applied to measured NF data (Ey and Ez) over the measurement surface.

In Figs. 11 we compare the measured NF at xmeas = 10 cm with the CST MWS NF data after the calibration step. The agreements are very satisfactory for Ey component meanwhile the cross-polarization (Ez) shows some differences compared with the CST results. This is due to many reasons: in part to the antenna fabrication accuracy, measurement positioning errors (alignment), cross-polarization of the probe is not well defined, probe sensitivity for low field levels….

Using the NFNF routine we transform the NF data measured at plane xmeas = 10 cm to calculate the field at the spatial positions A(r = 2.44 m, θ = π/2, φ = 0), B(r = 2.44 m, θ = π/2, φ = π/6), C(r = 2.44 m, θ = π/3, φ = 0) and D(r = 2.44 m, θ = π/6, φ = 0). The reference antenna radiation patterns used for these comparisons have been measured at Supelec cylindrical NF measurement facility placed in an anechoic chamber of dimensions 6 m x 7 m x 8 m. This measurement facility is presented in Fig. 12 and the E-field is measured using the same probe as the planar NF range (En-Probe EFS-105) at the plane cuts θ = π/2 and φ = 0. For our comparison, we consider the CST NF data (xmeas = 10 cm) which are transformed to reach the spatial positions A, B, C and D. These comparisons are carried out for the Ey components (copolar) and are presented in Figs. [START_REF] Andrieu | Transient characterization of a novel ultrawide-band antenna: the scissors antenna[END_REF]. As it is seen, good agreements are noticed between the measurement results. The NFNF transformation yields satisfactory results when using NF data measured at 10 cm from the AUT and transformed in order to calculate the transient response at the positions A, B, C. The differences observed for the position D are due to the NF measurement surface truncation error. The CST MWS NF collected at 10 cm is also transformed to calculate the field at the positions A, B, C and D. This has been done in order to verify that our calibration procedure stay correct and no additional calibration step is needed to compensate for the probe spatial response (the probe acts actually as a point source).

The NFNF transformation allows calculating Ex, Ey and Ez components in the half space x ≥ xmeas in front of the AUT using the tangential components measured at a given distance xmeas. Indeed, using Vivaldi antenna transient NF measured in the planar NF facility at xmeas = 10 cm from the AUT we perform NFNF transformation to set the transient field over the planar surface at x = 40cm (-60 cm ≤ y ≤60 cm, and -60 cm ≤ z < 60 cm). These are compared with the CST MWS results at the plane cuts (x = 40 cm, y = 0, -60 cm ≤ z ≤ 60 cm) and (x = 40 cm, -60 cm ≤ y ≤ 60 cm, z = 0). As is it seen in Figs. 14, the NFNF transformation results fit well the simulated Ey component (co-polar). The normal component Ex for the plane cut y = 0 agrees well with the simulated results. However, low level fields (cross-pol) present many discrepancies compared with the CST MWS results. Generally, the results of NFNF and NFFF transformations present some difficulties to set accurately the low-level crosspolarization field since the cross-polarized NF component is difficult to measure accurately with the EFS-105 probe (Fig. 10). accuracy for far-field calculation.

-

  AbstractIn this paper we present the experimental validation of Time-Domain (TD) Near-Field to Near-Field (NFNF) and Near-Field to Far-Field (NFFF) transformations for UWB antenna transient characterization. The used computation schemes for near-field to near-or far-field transformations are based on the Green's function representation of the radiated field where NF and FF are directly calculated in the time domain. The first step of the validation process comprises the electromagnetic simulation results dedicated to evaluate the accuracy of NFFF and NFNF transformations. The second step uses near-field measurement data of a Vivaldi antenna to validate the developed computation schemes whereas the measurement and calibration procedures are fully described. The measured NF data using Supelec planar and cylindrical near-field facilities are also compared with the electromagnetic simulation transient results of the Vivaldi antenna. Index Terms-Time domain, antenna measurement, near-field far-field transformation,

FFig. 3 .

 3 Fig. 3. The time domain excitation signal applied to the Vivaldi antenna (left). The excitation signal spectrum (right), the maximum of the signal spectrum is at 1.75 GHz

Fig. 1 .

 1 Fig. 1. The S11 of the Vivaldi antenna used for our transient analysis (optimized for the GPR application). This present the comparison between the S11 of the antenna modeled in the CST MWS and the measured one (fabricated antenna).

Fig. 2 .

 2 Fig. 2. Near-field measurement configuration at the distance xmeas from the AUT. The tangential components Ey and Ez are regularly recorded over the scan surface in the zy plane. The far-field validity area is dependent on the angular domain described by αθ and αφ.

Fig. 4 .

 4 Fig. 4. The Vivaldi antenna Ey(V/m) and Ez(V/m) radiated near-field as a function of the time at xmeas=30 cm for the plane cuts zmeas = 0 and ymeas = 0

  direction (θ, φ) is written as a function of the time derivative of the tangential measured NF meas Et  . The spherical coordinates (θ, φ) are associated with the planar coordinates system (x, y, z),

Fig. 5 .Fig. 6 .

 56 Fig. 5. The transient Etheta(V/m) and Ephi(V/m) far-field radiation pattern comparison considering different sampling frequencies for the spatial directions: (θ = π/2, φ = 0), (θ = π/4, φ = 0) and (θ = π/2, φ = π/4). The NF data are collected at xmeas=30 cm.

Fig. 7 .

 7 Fig. 7. The error values of Ey and Ez at the plane cuts y = 0 and z = 0 resulting from the NFNF transformation. The NF data is measured at xmeas = 20 cm and transformed to calculate the transient field at 50 cm from the AUT.

Fig. 10 .

 10 Fig. 10. (left) The Supelec planar measurement setup used for NF antenna measurement. The Vivaldi antenna is placed at 10 cm from the measuring probe (Enprobe EFS 105). The probe displacements (translation in y and z directions) are controlled by automated process. (right) The Vivaldi antenna transient radiated field (Ey (left) and Ez (right)) in the boresight direction is presented and compared with the CST MWS results after the calibration procedure (normalization).

Fig. 9 -

 9 Fig. 9-a The error values of Ey and Ez at the plane cuts y = 0 and z = 0 resulting from the NFNF transformation: the NF data is measured at xmeas = 10 cm and transformed to calculate the field at 50 cm from the AUT Fig. 9-b The error values of Ey and Ez at the plane cuts y = 0 and z = 0 resulting from the NFNF transformation: the NF data is measured at xmeas= 30 cm and transformed to calculate the field at 50 cm from the AUT

Fig. 11 .Fig. 12 .

 1112 Fig. 11. The comparison between the measured Vivaldi near-field (Ey (V/m) and Ez (v/m)) at the distance xmeas= 10 cm from the Vivaldi antenna with the CST MWS results for the plane cuts ymeas= 0 and zmeas= 0

Fig. 14 .

 14 Fig. 14. The NFNF results at x = 40 cm compared with the CST MWS results at the plane cuts y = 0 and z = 0. NF data (Ey and Ez) have been measured at xmeas = 10 cm from the AUT and transformed to reach the planar surface x = 40 cm. The comparison comprises Ex, Ey and Ez transient field components.

Fig. 13 .

 13 Fig. 13. The NFNF results at the spatial points A(r = 2.44 m, θ = π/2, φ = 0), B(r = 2.44 m, θ = π/2, φ = π/6), C(r = 2.44 m, θ = π/3, φ = 0) and D(r = 2.44 m, θ = π/6, φ = 0) compared with the CST MWS and the directly measured field in the cylindrical facility.

  , the measurement surface is defined for In our analysis, we use the reconstruction formula to interpolate and assess the time derivative of the measured field at t R c  . Let us consider NF data measured over the scan surface at xmeas from the AUT. This NF is a band limited function that can be accurately interpolated using the reconstruction formula. Namely, the field

		D	ymax		y	meas		D	ymax		and			zmax D		meas z		zmax D	. The NF is
	sampled regularly using the Nyquist criterion	/ max t  	,
	00 y z    	/ max  c	with			max	being the maximum angular
	pulsation of the AUT excitation signal, Ny=Dymax/Δy0 and
	Nz=Dzmax/Δz0.
		From (1) the first step in NFNF transformation is the
	interpolation of the measured field meas 0 ( E r	mn	, ) t and its time
	derivative to set their values at the instant t R c 	. The
	accuracy of NFNF transformation depends on the chosen
	interpolation method. 0 ( E r	mn	, ) t at the instant t in the
	time window	00 t t t   		N	t	 1   is calculated using the t
	field samples		(	0	,	0					)

mn E r t l t

 with the cardinal series as follows

[6][17]