0 EURASIP Journal on

Advances in Signal Processing

a SpringerOpen Journal

RESEARCH Open Access

Channel estimation for MIMO multi-relay systems using a tensor approach

Xi Han1^, André LF de Almeida2^ and Zhen Yang3t

Abstract

In this paper, we address the channel estimation problem for multiple-input multiple-output (MIMO) multi-relay systems exploiting measurements collected at the destination only. Assuming that the source, relays, and destination are multiple-antenna devices and considering a three-hop amplify-and-forward (AF)-based training scheme, new channel estimation algorithms capitalizing on a tensor modeling of the end-to-end communication channel are proposed. Our approach provides the destination with the instantaneous knowledge of all the channel matrices involved in the communication. Instead of using separate estimations for each matrix, we are interested in a joint estimation approach. Two receiver algorithms are formulated to solve the joint channel estimation problem. The first one is an iterative method based on a trilinear alternating least squares (TALS) algorithm, while the second one is a closed-form solution based on a Kronecker least squares (KRLS) factorization. A useful lower-bound on the channel training length is derived from an identifiability study. We also show the proposed tensor-based approach is applicable to two-way MIMO relaying systems. Simulation results corroborate the effectiveness of the proposed estimators and provide a comparison with existing methods in terms of channel estimation accuracy and bit error rate (BER).

Keywords: Channel estimation; Cooperative MIMO; Relaying; Tensor decomposition

1 Introduction

Cooperative communications have been considered as a promising concept to improve the link performance in modern wireless communication systems due to spatial diversity gains, enhanced coverage, and increased capacity [1-4]. In this context, relaying has been commonly accepted as a key technique to improve system performance by overcoming channel impairments, such as fading, shadowing, and path loss, in wireless fading channel environments [4-6]. By resorting to relay-assisted cooperation, multiple wireless links between mobile stations and base stations are established to create a virtual multiple-input multiple-output (MIMO) system [7]. In the simplest relay processing strategy, the relay stations amplify and forward the received data towards the base station. In this work, we adopt amplify-and-forward (AF) relaying due to its simplicity of implementation [5]. This strategy is preferable when fixed relay stations have

Correspondence: jinlucky333@gmail.com +Equal contributors

11nformation Network Center, Beijing University of Posts and Telecommunications, Beijing, China

Full list of author information is available at the end of the article

a limited computation capacity as opposed to the base station.

The overall link reliability of cooperative diversity schemes strongly depends on the accuracy of channel state information (CSI) associated with the multiple hops involved in the overall communication. Moreover, the use of common precoding techniques at the source and/or relays generally requires instantaneous CSI knowledge of the different channels to optimize transmission [8,9]. In practice, however, the CSI is unknown and has to be estimated with the aid of training sequences [10,11]. For two-hop relaying systems, the associated channel matrices can be estimated in separate LS estimation stages that operate sequentially at the destination [10]. When the communication involves additional hops, such a sequential LS estimation approach still applies by using additional transmission phases. The main problem is that channel estimation errors accumulate across the consecutive stages. In [11], a closed-form solution was proposed for the joint estimation of the channel matrices in a two-hop MIMO relaying system, avoiding error propagation.

Springer

© 2014 Han et a l.; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.

A few recent works have developed efficient receiver algorithms based on tensor analysis for channel estimation and/or symbol detection in cooperative systems [12-16]. In [12], a training sequence-based channel estimation algorithm is proposed for two-way relaying systems with multiple antennas at the relays. Recently [14], a channel estimation algorithm based on parallel factor (PARAFAC) model [17,18] was developed for two-hop MIMO relay systems. The approach allows estimation of the channel matrices associated with both hops by resorting to training sequences. Other few recent works have developed tensor-based receivers for one-way two-hop cooperative systems [13,15,16]. In particular, the approach of Ximenes et al. [16] assumes a Khatri-Rao space-time (KRST) coding [19] at the source node, and a semi-blind receiver is proposed by assuming the existence of a direct link between the source and the destination.

The approach of Roemer and Haardt and Rong et al. [12,14] allows a joint estimation of the channel matrices by resorting to training sequences. With the idea of avoiding the use of training sequences at the users' and relays' transmissions, the work [13] proposed a blind receiver for uplink multiuser cooperative diversity systems based on a PARAFAC model for the received signal. However, [13] is limited to a clustered relaying scenario, where relays belonging to the same cluster have the same spatial signature. The common feature of all these works is on the assumption of only two hops (source-to-relays and relays-to-destination). To further extend the coverage area and combat channel impairments such as path-loss and shadowing, it may be advantageous to introduce an additional hop along with an extra communication phase by means of three-hop relaying [5]. We highlight that the interest of the proposed work is on the joint channel estimation problem (i.e., joint channel and symbol estimation is not addressed here). The joint channel estimation problem was addressed in [12] for a two-way relaying system and in [14] for a one-way two-hop system. From a tensor modeling viewpoint, the common feature of both works is on the use of the PARAFAC model. Herein, we focus on a one-way three-hop multi-relay system, while resorting to a PARATUCK2 model to derive the proposed algorithms.

In this work, novel channel estimators are proposed for MIMO multi-relay systems. Assuming that the source, relays, and destination are multiple-antenna devices and considering a three-hop AF-based training scheme, new channel estimation algorithms capitalizing on a multilinear structure of the end-to-end communication channel are proposed. The proposed approach is based on a PARATUCK2 tensor model [20] of the data collected at the destination only, which allows the channel matrices to be jointly estimated at the destination. Two receiver algorithms are formulated to solve the channel estimation

problem. The first one is an iterative channel estimation method based on a trilinear alternating least squares (TALS) algorithm derived from a PARATUCK2 tensor model of the received data, while the second one is a closed-form solution based on a Kronecker least squares (KRLS) factorization. The proposed approach provides an extension of the idea recently proposed in [14] to a more general scenario with two-tier relaying using MIMO AF relays. Identifiability of the channel matrices is also examined in this work, and a useful lower-bound on the channel training length is derived. In contrast to conventional pilot-assisted LS channel estimation, where the channel matrices are estimated separately in consecutive stages, our proposed algorithms make a more efficient use of cooperative diversity by providing a joint estimation of all the channel matrices. As will be clear later, such a joint channel estimation is possible due to the use of the tensor approach to model the end-to-end system.

In comparison with conventional (multi-stage) LS channel estimation [10], the proposed tensor-based estimators have two distinguishing features: i) they avoid accumulation of channel estimation errors since all the channel matrices are estimated simultaneously (either iteratively or in closed-form), and ii) they can operate under less restrictive (and more flexible) conditions on the required number of antennas at the relays and/or destination, as will be clear from our identifiability analysis. Our approach also includes the PARAFAC-based channel estimator of [14] as a particular case. We also show that the proposed tensor modeling approach copes with a two-way MIMO multi-relaying communication system, where the TALS and KRLS channel estimators can be applied.

This paper is organized as follows. In section 2, the system model and working assumptions are described. Section 3 formulates the proposed approach. The data model is recast using tensor analysis, and the two channel estimation algorithms (TALS and KRLS) are derived. Identifiability of the channel matrices is also examined in this section. In section 4, we provide an extension of the proposed tensor-based signal model to a two-way MIMO relaying scenario. Numerical results are presented and discussed in section 5, and the conclusions are drawn in section 6.

Notation: Scalars are denoted by lowercase letters (a, b,...), vectors as lowercase boldface letters (a, b,...), matrices as uppercase boldface letters (A, B,...), and tensors as calligraphic letters (A, B,...). AT and A^ stand for transpose and pseudo-inverse of A, respectively. To retrieve the element (i, j) of A, we use a(i, j). The ith row of A e CIxR is denoted as A(,,:) while its rth column is denoted by A(:,r). The operator D, (A) forms a diagonal matrix out of the ith row of A. The Khatri-Rao (columnwise Kronecker) product between A e CIxR and B e

CJxR, i.e., A ❖ B = [A(:,d ® B(:,1),...,Ar ® Bm] e

£ÏJxR

2 System model

We consider a three-hop MIMO AF communication system where the source node transmits information to the destination node with the aid of R1 relays in the first tier and R2 relays in the second tier. As shown in Figure 1, the source and destination nodes are equipped with Ns > 2 and Nd > 2 antennas, respectively, and half-duplex relays are considered. The qth relay of tier 1, which receives data from the source node, is equipped with Iq antennas, q = 1,..., R1, while the pth relay of tier 2, which receives data from tier 1 relays, is equipped with Jp antennas, p = 1,..., R2. The total number of antennas that transmit in second and third phases are denoted by N1 = Ii + •• •+Ir1 and N2 = J1 + • •• + Jr2, respectively.

Some key assumptions are now given: (i) relays are synchronized at the symbol level. More specifically, the timing offset is assumed to be within one symbol period, so that timing information is acquired only through some form of (rough) coarse synchronization; (ii) fading is assumed to be frequency flat, and the data block size is smaller than the channel coherence time so that the channel is considered as time invariant; (iii) the direct links between the source (resp. tier 1 relays) and the destination node are not availablea. This situation is evidenced in the current uplink of IEEE 802.16j.

2.1 Data model

The communication between source and destination is accomplished in three hops. In the first hop, the modulated signal vector us (t) e CNsX1 is transmitted to R1 relays. The received signal at the qth relay of tier 1 can be written as

yq (t) = HSq)Us (t) + vq (t)

the qth relay of tier 1, H

channel between the source and the qth tier 1 relay, and V?) (t) e CIqX1 is an additive noise vector. Noise samples are modeled as independent and identically distributed complex Gaussian random variables with zero mean and unit variance.

In the second hop, the source stops transmission and all the R1 relays of tier 1 amplify their received signals with diagonal AF matrices G(1),..., G(R1) and simultaneously forward the resulting signals to the tier 2 relays. The received signal vector at the pth relay of tier 2 is then given by

yp (t + 1) = £ H

G(q)y(q)(t) + vip) (t + 1) (2)

where H(p,q) e CJp xIq is the MIMO channel linking the R1 tier 1 relays to R2 tier 2 relays, while vi? (t + 1) e CIpx1 denotes the corresponding noise vector. In the third hop, the source and all tier 1 relays are silent, while the tier 2 relays process the received signal vector with the diagonal AF matrices J(1),..., J(R2) and forward their amplified signals to the destination. The received signal vector at the destination is then given by

yrd(t + 2) = £ HfP J^yi? (t + 1) + vrd (t + 2), (3)

where H^ e CNd xJp is the MIMO channel linking the pth tier 1 relay to the destination, and vrd (t + 2) e CNdx1 the corresponding additive noise term. Let us define the multi-relay (block) channel matrices

Hrd , ... , Hrd2

H(1,1)

e CNdxN2,

e CN2xN1, (5)

where y(q) (t) e CIqX1 is the received signal vector at

£iqxNs is the MIMO

H(R2,1) . . . H(R2,R1)

-ti-rr * *rr

HsTr ^ [HSr)T,...,HS^] e C

iNs xN1

and let G = bdiag[G(1),...,G™] e CN1xN1 and J = bdiag [J(1),..., e CN2xN2 be the two diagonal matrices that collect the AF coefficients of the overall multi-relay system. Using these definitions, and using (1) and (2), we can rewrite (3) as follows:

yrd (t + 2) = HrdJHrrGHsrUs (t) + vrd (t + 2),

where Vrd (t + 2) = Vsr (t) + v„ (t + 1) + (t + 2) is the total noise at the destination, which contains the filtered noise contributions from the multiple relays, with Vsr (t) = HrdJHrrGVsr, Vrr (t + 1) = HrdJVrr (t + 1),

Vrr (t + 1) = Vsr (t) = ^

vrr1)T (t + 1),..., vr?2)T (t + 1)

, vSf1)T (t)

Hrd = Y^:

where Yi e CNd xL2 is the received signal matrix at the destination during the first training stage. In the second

stage, S1d is transmitted from all tier 1 relays to the destination via AF processing at the tier 1 relays. Defining Y2 e CNd xL1 as the data received from tier 1 relays at the second training stage, an LS estimate of Hrr can be obtained as

Hrr = (Hrd J^ Y2SHd.

Finally, So is transmitted from the source to the destination via the two tiers of relays. The destination collects the received data in Y3 e CNsxLo. An estimate of Hsr is then found as

vsr" (t), •

Note that, since this work is concerned with channel estimation, the AF matrices G and J cannot be optimized at the transmission (source and relays). Therefore, for simplicity, we have assumed that these matrices are diagonal. The use of non-diagonal AF matrices in the proposed approach is left for a future work. Note also that, once the channels are estimated, the design of full AF matrices can be done, e.g., based on the SVD of the channel matrices, following the idea of [9] or on the mean-square error (MSE) criterion [21]. If simplified AF schemes are used, where only power allocation is done, G and J are diagonal matrices, the coefficients of which can be designed as a function of the mean channel and noise powers [5] or optimized from power allocation strategies, as shown recently in [22].

2.2 Conventional LS estimation method

The simplest approach to estimate the effective channel Heff = HrdJHrrGHsr (including the amplifying factors) is based on training sequences. If separate estimations of the multi-relay channels Hrd, Hrr, and Hsr are required, for instance, to optimize the source precoding matrix and the relays' AF matrices, three separate LS estimation stages should operate sequentially at the destination. The method would work similarly to that of Kong and Hua [10]. Denote So e CNsXLo as the training sequence matrix sent by the source node, while S1d e CNixLi and S2d e CN2XL2 are the training sequence matrices sent by the relays at tiers 1 and 2, respectively. Assume that orthogonal training sequences are used in all stages, which implies training sequences of length L0 > Ns, L1 > N1 and L2 > N2 at the source, tier 1 and tier 2 relays, respectively. In the first stage, S2d is transmitted from all tier 2 relays to the destination. The LS estimate of Hrd is obtained as

Hsr = (HrdJHrrG)^ Y3sH.

This method requires 6 transmission phases to provide the destination with all the channel matrices (1 phase for estimating Hrd, 2 phases for estimating Hrr and 3 phases for estimating Hsr). Note that the channel estimation errors accumulate across the consecutive stages, due to the dependency between successive channel estimates. Moreover, this method requires Nd > N2 > N1 for the uniqueness of the LS estimates of Hrr and Hsr. In the following, we adopt a different path to solve this problem by capitalizing on tensor analysis. The idea is to provide the destination with a joint estimate of all the partial channels Hrd, Hrr, and Hsr by exploiting the tensor structure of the end-to-end signal model. The proposed approach allows channel estimation to be performed under less restrictive conditions on the number Nd of receive antennas at the destination compared with the conventional LS estimator, while avoiding error accumulation.

3 Proposed approach

In order to derive the proposed channel estimators, we first recast the formulation of the system model by resorting to multi-way (tensor) analysis. First, let us divide the overall training period into K time blocks. In every time block, the same training sequence matrix S0 e CNsXLo is transmitted by the source node. In the //th time block, the relays of tiers 1 and 2 use the AF matrices G/ and J/, respectively, // = 1,..., K. Let us define E e CKxN1 and F e CKxN2 as channel training matrices such that D/(E) = G/ and D/(F) = J/, where D/(•) forms a diagonal matrix out of the / th row of its matrix argument. Otherwise stated, the rows of E (resp. F) hold the AF coefficients of the R1 (resp. R2) relays associated with the different time blocks. Then, the signal received at the destination during the /th time block can be written as:

Y/ = HrdD/ (F) HrrD/ (E) HsrSo + V/, (11) // = 1,..., K,

where Vk = HrdDk (F) HrrDk (E) Vsr,k + HrfDk (F) Vrr,k, Vsr,k e CN1 xLo is the noise matrix at the relays during the kth time block, Vrr,k e CN2 xLo is the noise matrix at the second hop relays for the k-th time block, and Vrd,k e CNd xLo is the noise matrix at the destination for the kth time block.

Regarding the structure of the channel training matrices E e CKxN1 and F e CKxN2, unless otherwise stated, their columns are chosen as length-K random sequences following a uniform distribution between [ —1, 1]. These sequences are defined beforehand and known at the destination node. With such a choice, the signals transmitted by the relays across the K time blocks have random phases and are subject to limited power fluctuations. Clearly, this design is not optimal for minimizing the mean square error of the channel estimation. Determining an optimum design for these matrices is a difficult problem and is not pursued in this work. Nevertheless, extensive computer simulations have demonstrated that this choice yields very good results. For convenience, we will come back later to the problem of choosing E and F from a channel iden-tifiability viewpoint. A more elaborated design of these matrices will be then proposed.

Upon reception of the data matrix Yk, k = 1,..., K, an unstructured estimate of the end-to-end channel during the kth time block is first obtained at the destination. Multiplying both sides of (11) with the known training sequence matrix Sj/ yields

H = Yk SH e CNd

= HraD^ (F) HrrDk (E) Hs

k = 1, ••• , K. Let us introduce

H = Hk + Vk SH,

■ + Vk SH,

Here, we show that this tensor model can be exploited to derive novel channel estimators for a cooperative MIMO relaying system. Now, let us define

H[1] = [vec (H1), ••• ,vec (Hk)] e CNdNsxK

where H[i] is a matrix 'unfolding' for the tensor H obtained by stacking column-wise its K slices. Define also

Wk = Dk (F) HrrDk (E) e CN2xN1.

Substituting (14) into (15), and applying property vec (ACB) = (BT <g> A) vec(C), we get

H[1] = (HsTr ® Hrd) [vec (W1), ••• ,vec (WK)]

= (HsTr ® Hrf) diag (vec (Hrr)) (et © FT) (17)

where ET © FT =

E(1,:) ® F(1,:),

, E(7K,:) ® FCK,:)

CN2N1 xK (18)

E(k,:) e C1xN1 (resp. F(k,:) e C1xN2) denote the kth row of E (resp. F), and © is the Khatri-Rao (columnwise Kronecker) product.

Applying property vec (Adiag(x)B) = (BT © A x we get from (17) the following expression:

vec(H[1]) = ^1vec (Hrr), (19)

(ET © F^ © (HsTr ® Hrd)

çNdNsK xN1N2

Hk = HrdDk (F) HrrDk (E) Hsr, k = 1,..., K, (14)

is the matrix-of-interest that represents the effective end-to-end channel, Vk is the total noise matrix, and Hk is the noisy observation of Hk. We can assemble the set {H1, ••• , HK} to form a three-way array, or a third-order tensor, H e CNdxNsxK, whose dimensions are Nd (first dimension), Ns (second dimension), and K (third dimension).

Equation (14) corresponds to a PARATUCK2 model of the (noiseless) tensor H [23]. The PARATUCK2 model has first appeared in [20]. A more comprehensive formulation is given in [23], which also details an alternating least squares procedure for estimating its matrix factors.

In addition to the matrix unfolding H[1], it is useful to define two other matrix unfoldings, which collect the information of tensor H. Therefore, let us now define

CNdKxNs, h[3] =

e cNsKxNd.

From (14) and (16), it follows that

H[2] = Hsr =

_HrdWK_

- HsJrWj -

H[3] = Hjd =

-HJWK_

r hJ,.

or, more compactly,

H[2] = (Ik ® Hrd) &2H

H[3] = (Ik ® HsJr) QbH

Û2. =

CN2K XN1, H3 =

g CNiK xN2 (26)

3.1 Identifiability of channel matrices

Identifiability of Hsr, Hrr, and Hrd in the LS sense from H[i], H[2], and H[3] (see Equations (19), (24), and

(EJ © FJ)

(25)), respectively, requires that fi1 =

(Hjr ® Hrd)] e CndnkxNiN2, ZZ[2] = (Ik ® Hrd) ^ e CNdKxNi and Z[3] = (Ik ® HJ) &3 e CN»Kxn2 be full column-rank. These requirements come from the fact that fi1, Z[2], and Z[3] must be left-invertible, from which the following necessary conditions are obtained:

NdNsK > N1N2, NdK > N1, NsK > N2.

From the three inequalities and from the fact that we must have K > 2, the lower bound on the number K of time blocks necessary for identifiability is given by

K > max

N1N2 N1 N2

NdNs Nd Ns

where ¡x~\ is equal to the smallest integer that is greater than or equal to x.

Note that the identifiability of the channel matrices Hsr, Hrr, and Hrd from the unstructured channel tensor H will ensure that the compound channel Hc = HrdHrrHsr e CNdxNs is strictly unique. Note also that conditions NdNsK > N1N2 and NdK > N1 are clearly much less restrictive in terms of the required number Nd of antennas at the destination node, in comparison with

the conventional three-step LS estimator that requires Nd > N2 > N1. Otherwise stated, estimation of the partial channels can be done even in situations where the number of receive antennas is much less than the number of relay antennas (provided that K satisfies condition (28)). This situation may arise in scenarios with denser deployments of relay stations, where the total number of relay antennas exceeds those of source and/or destination antennas. As shown by these inequalities, the possibility of affording fewer receive antennas is compensated by an increase on the number K of training blocks, which represents a trade-off.

Condition (28), although necessary, is not sufficient for identifiability. Since ZM = (Ik ® Hrd) &2 e CNdKxNl and Z[3] = (IK ® HsTr) a3 e CN»KxN2, additionally, must have rank(fi2) = N1 and rank(fi3) = N2, i.e., both H2 e Cn2kxNi and Q3 e CNikxN2 must be full column-rank. Otherwise, Z[2] and Z[3] will be rank-deficient, even if (28) is respected.

Let us assume that the partial channels Hsr, Hrr, and Hrd are full rank matrices, which is a reasonable assumption when the wireless links are assumed to undergo scattering-rich multipath propagation. The following corollaries can then be obtained:

C1 If N1 = N2, identifiability of the partial channels is guaranteed for N1 < Ns and N2 < Nd;

C2 If N1 = 1, identifiability of the partial channels is guaranteed for N2 < Nd and N2 < K;

C3 If N2 = 1, identifiability of the partial channels is guaranteed for N1 < Ns and N1 < K.

Remark: For the first corollary, we can note that if N1 < Ns and N2 < Nd, then H^ ® Hrd is full column-rank, which ensures that fi1 e CN°NsKxN2 is full column-rank due to its Khatri-Rao product structure [24]. Likewise, q2 e Cn2kxN1 and Q3 e CN1kxN2 are also full column-rank in this case, guaranteeing the identifiability of the channel matrices. Regarding the second corollary, it corresponds to a special case of our system model where the first relay tier reduces to a single-antenna relay. In this case, satisfying N2 < Nd and N2 < K ensures that fi1, Z[2], and Z[3] are all full column-rank, so that the three partial channels are identifiable. The same reasoning is valid for the third corollary, which is analogous to the second one.

3.2 Essential uniqjueness

Let {Hsr, Hrr, Hr^ be an alternative set of matrices yielding the same unstructured channel tensor H satisfying the PARATUCK2 model (14). If Hsr, Hrr, and Hra are full rank and the identifiability conditions (27) are satisfied, then Hsr, Hrr, and Hrd are essentially unique. In this case, we have Hsr = AsrHsr, Hrd = Hrd Ard

and Hrr = A^H „A^ where the following relation holds:

(ASrA«) ® (AraA[2^ = IN1N2.

Note that permutation ambiguity does not exist due to the knowledge of the training matrices E and F. The relation (29) can be obtained by replacing the alternative solutions Hsr, Hrd, and Hrr into (14) and then applying some basic manipulations using properties of the Kronecker product. Equation (29) turns into the following relations: ArdA(2 = aIN2 and AgrAr1 = (1/a)IN1, where a is an arbitrary scalar factor. These two relations come from the fact that the Kronecker product between any two diagonal matrices is equal to the identity matrix if and only if these diagonal matrices are (scaled) identity matrices that compensate each other. Consequently, Hsr, Hrr, and Hrd can be recovered in an essentially unique manner up to scaling factors. The scaling ambiguity can be eliminated by normalizing the first column of Hsr or the first row of Hrd to one. Since these ambiguities compensate each other, the compound channel is strictly unique and we have Hc = HrdHrrHsr = HrdHrrHsr = Hc.

3.3 Trilinear alternating least squares algorithm

The TALS algorithm is an iterative estimation method that alternates among the LS estimations of the channel matrices Hsr, Hrr, and Hrd by fitting a PARATUCK2 model from the noisy matrices H[,] = H[,] + V[,], i = 1,2,3. Note that the noise term V[,] is constructed in a way analogous to H[,], i = 1,2,3, following Equations (15) and (21), respectively. The AF training matrices E and F are assumed to be known at the destination and are fixed during the estimation process. From (19), (24), and (25), we respectively obtain the following linear optimization problems:

TALS algorithm (direct estimation of Hsr, Hrr, and Hrd)

1. Set n = 0; _ _ Initialize randomly Hrr (n = 0) and Hrd (n = 0) ; Using (16) and (26), construct Sl2(n = 0) and Q.3(n = 0), respectively;

2. n ^ n + 1;

3. Find an estimate of Hsr using (22), by solving the LS problem

argmin || H[2] — (Ik ® H^) &2H.

■sr ll f

the solution of which at the nth iteration is given by

Hsr (n) = [(Ik ® Hrd (n —1)) H2 (n—1)]tH[2]

4. Find an estimate of Hrd using (23), by solving the LS problem

argmin

H[3] — (Ik ® H

® HsTr) fisHTi

the solution of which at the nth iteration is given by

Hrf (n) =

(Ik ® H^ (n)) H3 (n—1)

5. Find an estimate of hrr = vec (Hrr) using (17), by solving the LS problem

vec (H[1]) — (ET © FT)T © (Hjr ® Hrd)

argmin

rr || F

the solution of which at the nth iteration is given by

Hrr(n) = (et © F^T © (Hr (n) ® Hrd (n)) f x vec (H[1])

6. Rebuild H2 (n) and (n) and repeat Steps 2to5 until convergence.

argmin |vec (H[1^ — H1vec(Hrr)|F,

vec(Hrr)

argmin 1 H[2] — (Ik ® Hrd) &2Hsr|2,

argmin

H[3] — (Ik ® HTr) HsHTd

These LS estimation problems can be solved alternately by estimating one channel matrix at each time, while fixing the other matrices to their values obtained in previous estimation steps. Therefore, each iteration of the algorithm has three estimation steps. The algorithm starts by randomly initializing two out of the three channel matrices and proceeds until convergence. In the following, a summary of the TALS algorithm is provided.

Define e(n) = vec (Hw) — (ET © FT)T © (HTr(n)®

Hrd(n)) hrr(n). The sum of squared residuals (SSR) at the end of the nth iteration is defined as SSR(n) = eH (n)e(n). We declare the convergence of the algorithm when | SSR (n) — SSR (n — 1)| < 10—6, meaning that the model reconstruction error does not significantly change between two successive iterations.

Generally, the ALS algorithm is sensitive to the initialization, and convergence to the global minimum can be slow when all the matrix factors of the model are unknown [25]. However, in our case, we have observed that convergence to the global minimum is always achieved (e.g., within 10 to 30 iterations for medium-to-high SNRs)

due to the knowledge of the AF training matrices E and F.

3.4 Kronecker least squares algorithm

We now derive a closed-form solution to our channel estimation problem by exploiting the mixed Kronecker/ Khatri-Rao factorization structure of the matrix unfolding H[i] defined in (17). Starting from (13), the noisy version of (15) is given by:

H[1] = H[1] + (SH ® In,) V[1],

where V[i] = [vec(Vi),... ,vec(V/C)] e CNdNsxK. Let Z = ET © FT e CNlN2xK denote the combined AF training matrix and assume that ZZH = In1n2. Multiplying both sides of (33) by ZH, we have:

X[i] — H[i]ZH + (SH ® InJ V[i]Z

where X[1] = X[1] have:

(SH ® InJ V[1]ZH. From (17), we

X[1] = (HT ® Hrf) diag (vec (H„)).

x« —

(hsTT(:, «1) ® Hrd(:, «2)) hrr(n2, «1 ) (36)

Defining XXH1,H2 = unvec(xHi,H2) e CNdxNs as a rank-one matrix obtaining by reshaping, we have

XX«1,«2 — hr.

(«2, n1) Hrd (:, «2)Hsr(«1,:)

Consider the singular value decomposition (SVD) of

XH1,H2 = UH1,H2 ^«1,«2 VH1,n2

Hi = 1, ..., N1, «2 = 1, •••, N2. From the rank-one property of XH1,H2, we have:

Hr«2) (:, H2) = U„i,„2 (:, 1), Hi = 1,...,N1, HsHi)(Hi,:) = (V„i,„2 (:,1))r, «2 = 1,..., N2, hrr(H2, Hi) = X„i,„2 (1, 1).

Final estimates of Hrd(:, n2) and Hsr(:, n1) can be obtained by averaging over the Ni and N2 independent estimates, respectively:

1 N1 ^ Hrd(:, «2) — NiJ2 X

Our goal is to directly identify the channel matrices from (35). However, let us first address the deterministic design of the AF training matrices E and F such that ZZH = (ET © FT) (ET © FT)H = INlN2. Assuming K > N1N2, this condition is satisfied by designing Z, for instance, as a discrete Fourier transform (DFT) matrix. Having fixed the structure of Z, we are left with the problem of factorizing this matrix as the Khatri-Rao product between ET and FT. This problem can easily be solved by means of K rank-one matrix factorizations, which admit unique solutions. Note that the /cth column of Z can be written as

Z(:, k) = (E(//,:) ® F(//, :))T e CNlN2x1, k = 1,...,K.

Defining a rank-one matrix Z/ = unvec(Z(:, k)) e CN2xN1, it follows that

Z k = (F(k, :))T E(k,:),

from which E(k,:) and F(k,:) can be determined as the unique right and left singular vectors of Zk, k = 1,..., K. Note that the proposed design, although not optimized to minimize the mean square error of the channel estimation, ensures that the noise characteristics in (33) will not be changed when H[1] is post-multiplied by ZH (i.e., inverse DFT transformation).

Coming back to the channel estimation problem, from (35), let us define x^,n2 e CNNdx1 as the [ (m - 1)N + n2]-th column of %]'ee CN»NdxN1N2, m = 1,...,N1, n2 = 1,..., N2. Note that

Hsr(«1, 0 — N2E H «2—1

T«1)(:, «2),

«1—1

(«1,:),

Hrd — [Hrd(:,1), -.., Hrd (:, N2)],

Hsr — [Xsr(:, 1),..., Hsr(:, N1)]r, ^1,1 (1,1) ••• ^N1,1 (1,1)

.^1,N2 (1,1) ••• ^N1,N2 (1,1)

Note that the columns of the estimated Hsr and Hrd have unit energy while each entry of Hrr concentrates all the energy of the wireless link connecting the source node to the destination node via a given tier 1-tier 2 relay pair. Such an interpretation is useful for designing transmit and receive spatial filters for system optimization as well as for power allocation purposes.

Discussion: The KRLS algorithm involves the computation of N1N2 SVDs to provide rank-one approximations for the matrices X1,1,..., Xn1,n2, of dimensions Nd x Ns, which are constructed from the NlN2 columns X[1]. The distinguishing feature of the KRLS-based estimator is on the closed-form solution to the problem, as opposed to the TALS algorithm that consists of iterative LS estimation steps, which implies a higher computational complexity. However, note that the KRLS algorithm is only applicable under the condition K > N1N2, which is necessary for Z = ET © FT to have orthogonal rows, leading to (35). In contrast, the TALS algorithm can operate under a

«1 ,«2

much lower bound on K, as discussed in Section 3.1. This is clearly a trade-off between both estimators in terms of identifiability conditions and computational complexity. As will be shown in the next section, both estimators provide satisfactory performances, and the choice of the best estimator is rather dependent on the design constraints of the system. For instance, we can say that the TALS estimator is preferable if processing power at the receiver is not too limited, as is often the case with base station reception in outdoor micro- or macro-cells. The KRLS solution would be more likely chosen in indoor scenarios, where channel coherence time is long enough to allow for higher values of K.

4 Extension to two-way MIMO relaying systems

In the previous sections, we have focused on a multi-relay cooperative scheme, where transmission is directed in one direction, i.e., from a specific source to a specific destination via two tiers of multiple relays. In this section, we show that the same modeling approach can be extended to a two-way MIMO relaying scenario, where pilot/data transmission takes place in both directions. In the first phase, two sources simultaneously transmit their data to the multiple relays. Note that, in the two-way case, the relays of each tier receive a superposition of Ns1 + Ns2 signals coming from sources 1 and 2. In the second and third phases, inter-relay communication takes place. More specifically, in phase two, tier 1 relays transmit signals towards tier 2 relays, while tier 1 relays stay silent. In phase three, the opposite happens. Finally, in the fourth communication phase, all the relays transmit to the two sources, and each one of them receives a superposition of N1 + N2 signals.

In the first transmission phase, we assume that training symbol matrices Si e CNs1xL and S2 e CNs2xL are transmitted from sources 1 and 2, respectively. We omit the additive noise terms for convenience of presentation. The signal received at the ith relay tier is given by:

X(i) = И81Г; S1 + HS2r,. S2 = H(i)S, i = 1,2, (48)

where H(i) = [HS1ri HS2ri] e CNiX(N*1 +N*2), and S = [S[ST]t e C(Ns1 +n2)xl. the training sequence Si chosen by source i, is designed to satisfy the following conditions:

(i) SiSf = IN;, i = 1,2,

(ii) S^ = ON1XN2.

A possible construction satisfying these two conditions is based on the normalized DFT matrix of size L x (Ns1 + Ns2), with L > Ns1 + Ns2. This design allows the sources to eliminate the self-interference generated by their own transmission, when receiving the signal back from the relays.

In the second and third phases, where inter-relay communications happen, the signal received at the relays of tier i from the relays of tier j, (i, j) = {(1,2), (2,1)}, can be written as:

Z(i = Hrjri Dk (Ej)Xj

= Hrjri Dk (Ej)Hj)S, (49)

k = 1,...,K, where H,T/ e CNixNj is the MIMO channel

linking the relays of tier j at transmission to the relays of tier i at reception, (i, j) = {(1,2), (2,1)}. Note that channel reciprocity in the inter-relay communications is not a necessary assumption which means that we may have

Hr1r2 = Hr2r1 .

Finally, in the fourth transmission phase, the signals received at sources 1 and 2 are post-multiplied by SH and SH, respectively, to accomplish self-interference elimination, yielding

Yk1) = (HrmDk(F1)Z(k1]) S? + (Hr2s1Dk(F2)Zk2)) S?

= HrmDk (F1) Hr2r1 Dk (E2)Hs2r2

tier 2 — tier 1 relay path

+ Hr2s1 Dk (F2)Hr1r2Dk (E1) Hs2r1

tier 1 -— tier 2 relay path

= H(1,1)Dk(F 1,2)G[1)Dk(-E2,1 )HH(1,2), k = 1,...,k,

Yf = (Hr^Dk(F1)Zk1^ S? + (Hr2s2Dk(F2)Zk2^ S?

= Hr2s2Dk (F2 )Hr1r2Dk (E1)Hs1r1 tier 1 - tier 2 relay path + Hr1s2Dk (F1)Hr2r1Dk (E2)Hs1r2 tier 2 - tier 1 relay path

= H(2,1)Dk(F2,1)G(2)Dk(E 1,2)HH(2,2), k = 1,...,K,

Hii1,1) == [Hr1s1 Hr2s1] e CNS1 x(N1+N2) (52)

Hi(1,2) == [Hsr2r2 Hsr2r1 ]r e C(N1 +n2)xNs2 (53)

Hii2,1) = [Hr2s2 Hr1s2] e CNS2 x(N1+N2) (54)

HH(2,2) = [H^ H^r2]r e C(N1+N2)XNS1 (55)

Grr} == bloclcdiag (Hr2r1 Hr1r2) e C(N1+N2)X(N1+N2) (56)

Gr2) == blockdiag (Hr1r2 Hr2r1 ) e C(N1+N2)X(N1+N2) (57)

Fi,j == [Fi Fj], ^u == [Ei Ej], (i, j) = {(1,2), (2,1)}. (58)

Therefore, we can conclude that the signals received at sources 1 and 2 in the considered two-way MIMO relaying scenario (Equations (50) and (51)) follows a PARATUCK2 model. By analogy with the noiseless part of the one-way

signal model (14), we have the following correspondences between the factor matrices:

(Hrd, Hsr, Hrr) ^ (H(1,1), H(1,2), Grr))

(E, F) ^ (Fu, Ew) (source 1) (59)

(Hrd, HSr, Hrr) ^ (H(2,1), H(2,2), Gr2>)

(E, F) ^ (F2,1, E 1,2) (source 2) (60)

Consequently, the tensor-based channel estimation algorithms proposed in the previous section can be equally applied at each source to estimate the channels H(a), H(i,2) and G®, i = 1,2, from Equations (50) and (51), respectively. If reciprocity is assumed in the two-way relay channels, we have:

Hrisi = Hjri, i = 1,2 (61)

Hris = H?ri, (i, j) = {(1, 2), (2,1)}, (62)

Hrirj = H^, (i, j) = {(1,2), (2,1)}, (63)

which in turn implies H(1,1) = (H(2,2))T = Hs1, H(1,2) = (H(2,i))T = hS2, and Grr} = (g^^T = G. In this particular case, the PARATUCK2 models (50) and (51) become essentially equal, i.e., they depend on the same unknown channel matrices Hs1, Hs2, and G to be estimated. Note, however, that such a reciprocity is not a necessary assumption of our modeling approach, which can be used in the general case of non-symmetrical two-way MIMO relay channels.

5 Numerical results

We now present computer simulation results for assessing the performance of the proposed channel estimator in selected system configurations. The estimator's performance is evaluated in terms of the normalized mean square error (NMSE) of the estimated channel matrices. From the estimated channels, the performance in terms of bit error rate (BER) is calculated by assuming a linear receive filter. The BER and NMSE curves are plotted as a function of the overall signal-to-noise ratio (SNR) at the destination. This SNR is given by the ratio between the powers of the useful signal component and the noise component in Equation (11). For each simulated SNR value, the results represent an average over L = 5,000 Monte Carlo runs. At each run, the channel coefficients are drawn from a circularly symmetric complex-valued Gaussian distribution with zero-mean and unit variance, while the transmitted symbols are drawn from a BPSK sequence. The SNR level at the tier 1 and tier 2 relays are assumed to be 30 dB above the SNR level at the destination.

For purposes of performance evaluation, the scaling ambiguities affecting the estimates of the channel matrices are removed by assuming the first column of Hsr and first row of Hrd contain all one elements, similarly to [11,14]. These scaling ambiguities can be determined as follows. First, we find Asr = D1 (HsTr) and Ard = D1 (Hrd). Then, applying property (AB) ® (CD) = (A ® C)(B ® D) yields (ASr ® Ard) A> ® Ag]) = In1n2, from which we obtain A(1) ® A(2^ = A-1 ® A-1. A solution to this relation is then found as A^^ =

[D1 (HsTr)]_1 and A® = [D1 (Hrd)]_1.

In Figure 2, we depict the NMSE performance for the compound channel of our proposed estimators in comparison with the conventional LS estimator. The parameters are Ns = 2, N1 = 4, N2 = 4, Nd = 6, K = 16, L0 = 30, and the number of transmitted data symbols is N = 1000. We can see that TALS and KRLS have similar performances, which are considerably better than the conventional (three-stage) LS estimator. The worst performance of the LS estimator comes from the error accumulation across successive channel estimation stages, which degrades its overall NMSE performance.

Figure 3 shows the NMSE performance of our proposed estimators in comparison with the two-hop bilinear alternating least squares (BALS) estimator of Rong et al. [14]. This estimator is a special case of the proposed one, where only one tier of relays is used. In this case, model (11) reduces to

Yk = HrdDk (E) HSrSc + Vk, (64)

k = 1,..., K,

and the channel matrices Hsr and Hrd are estimated by means of a BALS algorithm. The parameter setting is the same as that of Figure 2. It can be seen the proposed estimator operates satisfactorily, being able to effectively estimate the three channel matrices. Figure 3 also indicates the proposed estimator performs close to the BALS estimator operating in a two-hop system. A small performance degradation is observed, which is due to the presence of an additional AF transmission phase of our three-hop system, resulting in a higher overall noise contribution at the destination. Note also that the TALS estimator involves three estimation steps while the BALS one has two estimation steps only.

Figure 4 shows the BER performance of a linear zero forcing (ZF) receiver designed from the estimated channel matrices, which are obtained from the TALS, KRLS, or the conventional LS estimators. The ZF receiver operates on data block collected in the received data matrix Y e CKNdxN. The length of the data block is N = 100 symbols, and the remaining system parameters are the

Figure 2 NMSE of estimated compound channel Hc. Proposed estimators (TALS and KRLS) vs. conventional LS estimator.

same as those of the previous experiment. The ZF filter output is given by:

HrdD1 (F) HrrD1 (E) Hs

, HrdDK (F) HrrDK (E) Hs.

This figure shows similar BER performances for TALS and KRLS, which are better than that of the conventional LS algorithm. This result corroborates the effectiveness of

our channel estimators when used with linear receiver for symbol detection. In Figure 5, we evaluate the impact of the number of relay antennas on the BER performance of a linear ZF detector using the proposed TALS channel estimator. The fixed system parameters are Ns = 2, Nd = 6, L0 = 30, and K = 10. It can be seen that the BER performance is considerably improved as the number of relay antennas is increased, corroborating the expected gains of cooperative diversity. Although not plotted in this figure, the BER curves of the KRLS estimator are similar to those obtained with the TALS one.

Figure4 BER performance of linear ZF detector (TALS and KRLS). BER performance of a linear ZF detector at the destination designed from the estimated channel matrices (TALS and KRLS).

Figure 6 depicts the performance of the ZF receiver designed from the perfect CSI for all channel matrices. Two parameter settings are considered, where Nd = 2 and 4, respectively. The other system parameters are fixed to Ns = 2, Ni = N2 = 3, Lo = 6, and K = 9. First, it can be seen that the BER performances are considerably improved as the number of antennas at the destination is increased, owing to the higher spatial diversity, as expected. From these results, we also find that the

TALS and KRLS estimators provide similar results and, more interestingly, their performances are close to that of the perfect CSI case. For instance, for a target BER of 10-1, the SNR gap with respect to the perfect CSI case is less than 2 dB.

6 Conclusions

We have proposed channel estimation algorithms for MIMO AF multi-relay systems. The proposed estimators

Figure 6 BER performance of linear ZF detector versus the number of antennas at the destination. BER performance of a linear ZF detector as a function of the number of antennas at the destination (TALS and KRLS).

are designed to provide the destination (base station) with the instantaneous CSI of all the channels involved in the communication. In contrast to conventional pilot-assisted channel estimation, the proposed algorithms make a more efficient use of cooperative diversity by providing a joint estimation of all the channel matrices thanks to the use of a tensor modeling of the end-to-end system. Such a joint estimation can be accomplished either iteratively (using TALS) or in closed-form (using KRLS). Our numerical results corroborate the effectiveness ofthe proposed algorithms. The TALS estimator has a higher computational complexity than the KRLS one due to its iterative nature. On the other hand, the minimum condition for operation of KRLS (K > N1N2) is more restrictive than the identi-fiability conditions of TALS, which implies more training (i.e., higher number of time blocks) to carry out the joint channel estimation. Both algorithms are suitable to the joint channel estimation problem, and a particular choice is mostly dictated by practical system requirements. We have also provided an extension of the proposed approach to two-way MIMO multi-relay system and verified that such an extension results in the same tensor model as the one-way scenario. Consequently, the proposed algorithms can be applied to one- and two-way multi-relay MIMO schemes.

Endnote

aSince our focus is on the relay channel, direct links are not considered for simplicity. However, the idea proposed in this work can be easily extended to include direct links.

Competing interests

The authors declare that they have no competing interests. Acknowledgements

André L. F. de Almeida is partially supported by CNPq and CAPES. This work was supported by the China's Next Generation Internet Project (CNGI Project) (CNGI-12-03-009) and DNSLAB. This work was also supported by National Natural Science Foundation of China (Grant No. 61 173017).

Author details

1 Information Network Center, Beijing University of Posts and Telecommunications, Beijing, China. 2Department ofTeleinformatics Engineering, Federal University of Ceara, Campus do Pici, B. 725,60455-970 Fortaleza, Brazil.3 School of Computer Science, Beijing University of Posts and Telecommunications, 100876 Beijing, China.

Received: 28 May 2014 Accepted: 24 October 2014 Published: 17 November 2014

References

1. A Sendonaris, E Erkip, B Aazhang, User cooperation diversity - part I: system description. IEEE Trans. Commun. 51 (11), 1927-1938 (2003)

2. A Sendonaris, E Erkip, B Aazhang, User cooperation diversity - part II: implementation aspects and performance analysis. IEEE Trans. Commun. 51(11), 1939-1948 (2003)

3. JN Laneman, DNCTse, GW Wornell, Cooperative diversity in wireless networks: efficient protocols and outage behavior. IEEE Trans. Inform. Theor. 50(12), 3062-3080 (2004)

4. L Cao, J Zhang, N Kanno, in Proc. IEEE GLGBECGM'09. Multi-user cooperative communications with relay-coding for uplink IMT-advanced 4G systems, vol. 1 no. 1 (Honolulu, HI, November 2009), pp. 1-6

5. KJR Liu, AK Sadek, W Su, A Kwasinski, Cooperative Communications and Networking (Cambridge University Press, New York, USA, 2009)

6. R Pabst, BH Walke, DC Schultz, P Herhold, H Yanikomeroglu, S Mukherjee, H Viswanathan, M Lott, W Zirwas, M Dohler, H Aghvami, DD Falconer, GP Fettweis, Relay-based deployment concepts for wireless and mobile broadband radio. IEEE Comm. Mag. 42(9), 80-89 (2004)

7. M Dohler, Y Li, Cooperative Communications: Hardware, Channel and PHY (John Wiley & Sons, West Sussex, United Kingdom, 2010)

8. Y Rong, XTang, Y Hua, A unified framework for optimizing linear nonregenerative multicarrier MIMO relay communication systems. IEEE Trans. Signal Process. 57(12), 4837-4851 (2009)

9. AToding, MRA Khandaker, Y Rong, Joint source and relay optimization for parallel MIMO relay networks. EURASIP J. Adv. Signal Process. 174,1-7 (2012)

10. T Kong, Y Hua, Optimal design of source and relay pilots for mimo relay channel estimation. IEEE Trans. Signal Process. 59(9), 4438-4446 (2011)

11. P Lioliou, M Viberg, M Coldrey, Efficient channel estimation techniques for amplify and forward relaying systems. IEEETrans. Comm. 60(11), 3150-3155 (2012)

12. F Roemer, M Haardt,Tensor-based channel estimation and iterative refinements for two-way relaying with multiple antennas and spatial reuse. IEEETrans. Signal Process. 58(11), 5720-5735 (2010)

13. CAR Fernandes, ALF de Almeida, DB Costa, Unified tensor modeling for blind receivers in multiuser uplink cooperative systems. IEEE Signal Process. Lett. 19(5), 247-250 (2012)

14. Y Rong, MRA Khandaker, YXiang, Channel estimation of dual-hop MIMO relay system via parallel factor analysis. IEEE Trans. Wireless Comm. 11 (6), 2224-2233(2012)

15. ALF de Almeida, CAR Fernandes, D Benevides da Costa, Multiuser detection for uplink ds-cdma amplify-and-forward relaying systems. IEEE Signal Process. Lett. 20(7), 697-700 (2013)

16. LR Ximenes, G Favier, ALF Almeida, YCB Silva, PARAFAC-PARATUCK semi-blind receivers for two-hop cooperative MIMO relay systems. IEEE Trans. Signal Process. 62(14), 3604-3615 (2014)

17. RA Harshman, Foundations of the PARAFAC procedure: model and conditions for an 'explanatory' multi-mode factor analysis. UCLA Working Papers Phonetics 16,1-84 (1970)

18. JD Carroll, J-J Chang, Analysis of individual differences in multidimensional scaling via an N-way generalization of "Eckart-Young" decomposition. Psychometrika 35(3), 283-319 (1970)

19. ND Sidiropoulos, RS Budampati, Khatri-Rao space-time codes. IEEE Trans. Signal Process. 50(10), 2396-2407 (2002)

20. RA Harshman, ME Lundy, Uniqueness proof for a family of models sharing features of Tucker's three-mode factor analysis and PARAFAC/CANDECOMP. Psychometrika 61 (1), 133-154 (1996)

21. BKChalise, YD Zhang, MG Amin, in Proc. ofSPIE'12. Joint Optimization of Source Beamformer and Relay Coefficients Using MSE Criterion, vol. 8404, no. 1, (May 2012)

22. M Mohammadi, M Ardebilipour, Z Mobini, R-A Zadeh, Performance analysis and power allocation for multi-hop multi-branch amplify-and-forward cooperative networks over generalized fading channels. EURASIP J. Wireless Commun. Networking 2013(1), 1-13 (2013)

23. R Bro, Multi-way analysis in the food industry: Models, algorithms and applications. PhD thesis, University of Amsterdam, Amsterdam, 1998

24. ND Sidiropoulos, X Liu, Identifiability results for blind beamforming in incoherent multipath with small delay spread. IEEE Trans. Signal Process. 49(1), 228-236 (2001)

25. A Smilde, R Bro, P Geladi, Multi-way Analysis: Applications in the Chemical Sciences (John Wiley & Sons, West Sussex, England, 2004)

doi:10.1186/1687-6180-2014-163

Cite this article as: Han etal.: Channel estimation for MIMO multi-relay systems using a tensor approach. EURASIP Journal on Advances in Signal Processing 2014 2014:163.

Submit your manuscript to a SpringerOpen journal and benefit from:

7 Convenient online submission 7 Rigorous peer review 7 Immediate publication on acceptance 7 Open access: articles freely available online 7 High visibility within the field 7 Retaining the copyright to your article

Submit your next manuscript at 7 springeropen.com