多輸入多輸出有限脈衝系統的盲判別與對等化(II)

行政院國家科學委員會專題研究計畫 期中進度報告

多輸入多輸出有限脈衝系統的盲判別與對等化(2/3)

中 華 民 國 94 年 5 月 31 日

Blind Identification and Equalization of MIMO FIR Systems (II)

A Progress Report

NSC-92-2213-E-009-085 NSC-93-2213-E-009-041 Principal investigator: Ching-An Lin

We describe briefly the research work in progress, the preliminary results obtained so far, and the expected final results at the end of the first two years’s research. Topics of research we have undertaken include

• Blind identification of MIMO channel using periodic precoding. • Reverse link power control in CDMA.

• Space-time block coding scheme

1

The basic idea is to induce cyclostationarity at the transmitter and to exploit the linear relation between the covariance of the received data. This is an extension of our pervious result on SISO channel. A preliminary version of the result had been presented at 2005 ISCAS (Kobe Japan). Quite a few people at the conference are interested in our result. One of the feature attract most attention: the ability of the proposed algorithm to identify the channel when the number of outputs is less than the number of inputs.

The complete result has been submitted to IEEE Transaction on Signal Processing (date: 2005/4/19). The manuscripts is attached at the end of this report.

Currently we are working on an algorithm for blind identification of MIMO channel without using periodic precoding. Preliminary simulation results are promising. We expect to finish the investigation and submit a paper at the end of this budget year.

2

In CDMA multi-user environment, the signal transmitted by a mobile will cause interference at the base station for other user. This is commonly called multiple access interference (MAI), which is due to the non-orthogonality of the codes that identify the individual user. The MAI makes power control important in maintaining a desired level of quality of service (QOS).

The current practice is to use a very conservation control strategy in order to guarantee stability at the expense of the achievable performance. In other words, if an more aggressive

strategy is used it may be possible to have faster tracking of the intended signal-to-interference (SIR) level. Of course, in a feedback loop an aggressive control strategy runs the risk of causing instability if the controller is not properly designed. We believe feedback control theory should be used to design power control strategy so as to improve quality of service or to increase the network capacity.

The power control problem can be modelled as a decentralized control problem, with in-terconnected dynamics and distributed control. The inin-terconnected dynamics has diagonal dominance structure: the correct correlations gives the intended strong links and the incor-rect corincor-rections can be modelled as weak links. There is stability theory based on diagonal dominance available. The purpose of our research is to investigate the possibility of using the theory do analyze the performance of power control algorithm. We will try to establish quantitative relation between the interference level and the achievable performance

We expected some preliminary results at the end of the current budget year.

3

This is a topic we just started to investigate. In the MIMO setup (multiple transmit and receive antenna), the coding is important in using the channel diversity while maintain a reasonable transmission rate. We will focus on the topics of identification and equalization in a space-time coded transmission.

We are currently studying the fundamental aspect of space-time coding scheme and we do not expect any significant process at the end of this budget year.

4

Blind Identification of MIMO Channels Using Optimal Periodic Precoding∗

Ching-An Lin and Yi-Sheng Chen

Department of Electrical and Control Engineering National Chiao-Tung University

Hsinchu, Taiwan

Abstract

We propose a method for blind identification of MIMO FIR channels that exploits cyclo-stationarity of the received data induced at the transmitters by periodic precoding. It is shown that, by properly choosing the precoding sequence, the MIMO FIR transfer functions, with Mt inputs and Mr outputs, can be identified up to a unitary matrix ambiguity. The

transfer functions need not be irreducible or column reduced, and there can be more outputs (Mr ≥ Mt) or more inputs (Mr < Mt). The method exploits the linear relation between the

covariance matrix of the received data and the “channel product matrices”. The method is shown to be robust with respect to channel order overestimation. The proposed algorithm requires solving linear equations and computing the nonzero eigenvalues and eigenvectors of a Hermitian positive semidefinite matrix. The performance of the algorithm, and indeed the identifiability, depend on the choice of the precoding sequence. We propose a method for optimal selection of the precoding sequence which takes into account the effect of additive channel noise and numerical error in covariance estimation. Simulation results are used to demonstrate the performance of the algorithm.

Key Words : MIMO channel, blind identification, transmitter induced cyclostationarity, periodic precoding

1

Introduction

Blind identification of SISO frequency selective channels exploiting transmitter induced cyclostationarity of the second-order statistics of the received data is first proposed in [1, 2]. Since then, various schemes have been proposed to induce cyclostationarity at the transmitter and to blindly identify SISO [3]-[8] and MIMO channels [9]-[13]. One way to induce cyclo-stationarity at the transmitter is by periodic precoding, i.e., multiplying the source symbols

with a periodic sequence before transmission [2, 5, 6, 8],[9]-[11]. For SISO channels, blind identification methods based on periodic precoding are shown to be robust with respective to channel order overestimation and impose no restriction on the locations of channel zeros [2, 5, 6, 8].

In the MIMO context, Chevreuil and Loubaton [9] proposes a scheme that multiplies each input by a constant modulus complex exponential precoding sequence to induce conjugate cyclostationarity at the transmitter. The scheme is used to reduce the MIMO channel iden-tification preblem to several SIMO ones, which are then solved by the subspace method [16]. Each SIMO channel is required to be free from common zeros and only real symbols can be used. B¨olcskei et. al. [10] proposes a method that can identify each of scalar channels up to a phase ambiguity using non-constant modulus periodic precoding sequences. The method imposes no restriction on channel zeros and is insensitivity to channel order overestimation. However, no general procedure for the design of the precoding sequences is given. The method is extended to the multicarrier case [11].

In this paper, we propose a method for blind identification of MIMO FIR channels using periodic precoding as a means to induce cyclostationarity. It is shown that, by properly choosing the precoding sequence, the MIMO FIR transfer functions, with Mt inputs and Mr

outputs, can be identified up to a unitary matrix ambiguity. The transfer functions need not be irreducible or column reduced [14, 15], and there can be more outputs (Mr ≥ Mt) or more

inputs (Mr < Mt). The method exploits the linear relation between the covariance matrix

of the received data and the “channel product matrices”. The method is shown to be robust with respect to channel order overestimation. The proposed algorithm requires solving linear equations and computing the nonzero eigenvalues and eigenvectors of a Hermitian positive semidefinite matrix. The performance of the algorithm, and indeed the identifiability, depend on the choice of the precoding sequence. We propose a method for optimal selection of the precoding sequence which takes into account the effect of additive channel noise and numerical error in covariance estimation. Simulation results are used to demonstrate the performance of the algorithm. The paper generalizes the results for the SISO case discussed in [8].

The paper is organized as follows. Section 2 is problem statement and formulation. In Section 3, we derive the identification method and propose the blind identification algorithm. In Section 4, we discuss optimal selection of the precoding sequence. Simulation results are given in Section 5. Section 6 concludes the paper.

Notations used in this paper are quite standard: Bold uppercase is used for matrices, and bold lowercase is used for vectors. AT _{represents transpose of the matrix A, and A}∗_represents

conjugate transpose of the matrix A. A ⊗ B is the Kronecker product of A and B. 0M ×N is

2

Problem Statement and Formulation

Figure 2.1: MIMO Channel Model

We consider the linear MIMO baseband model of a communication channel with Mt

trans-mitters and Mr receivers shown in Figure 2.1, where each source symbol sequence is multiplied

by an N -periodic sequence, p(n), before transmission. The transmitted signal is

wj(n) = p(n)sj(n), j = 1, 2, · · · , Mt (2.1)

where p(n + N ) = p(n), ∀ n. The discrete time model describing the relation between the transmitted signal wj(n) and the received signal xi(n) has the form of an MIMO FIR filter

with additive noise:

xi(n) = Mt X j=1 Lij X l=0 hij(l)wj(n − l) + vi(n), i = 1, 2, ..., Mr (2.2)

where hij(0), hij(1), · · · , hij(Lij), are the impulse responses of the channel between the jth

transmitter and the ith receiver, and vi(n) is the channel noise seen at the input of the ith

receiver. The equations (2.1) and (2.2) can be written more compactly as

w(n) = p(n)s(n), x(n) =

L

X

l=0

H(l)w(n − l) + v(n) (2.3)

where w(n), s(n) ∈ CMt_{, and x(n), v(n) ∈ C}Mr _{are vector signals formed by stacking the}

respective scalar signals together, e.g., x(n) = [x1(n) x2(n) · · · xMr(n)]

T_{. The ijth element}

of H(l) ∈ CMr×Mt _{is h}

ij(l), and L = maxi,j{Lij} is the order of the MIMO channel. Thus

H(L) 6= 0Mr×Mt. The following assumptions are made throughout.

(A1) s(n) and v(n) are white with zero mean vector sequences, and s(n) and v(n) are temporally and spatially uncorrelated. More precisely, E[s(k)s(j)∗

] = δ(k − j)IMt ∈

RMt×Mt, E[v(k)v(j)∗] = δ(k − j)σ2

vIMr ∈ R

Mr×Mr, E[s(k)v(j)∗] = 0

Mt×Mr, ∀ k, j, where

(A2) An upper bound ˆL of the channel order L is known and the period N > ˆL + 1. (A3) rank([H(0)T _H(1)T _{· · · H(L)}T_]T_)=M

t.

Due to periodic precoding, the input-output relation between the source s(n) and the received signal x(n), described by (2.3), is periodically time-varying. In order to obtain a time-invariant representation, we consider input-output relation between block input and block output of size N [17]. Define block signal ¯x(n) = [x(N n)T_{, x(N n+1)}T_{, · · · , x(Nn+N−1)}T_]T _∈

CMtN_{, and let ¯v(n), ¯w(n), ¯s(n) be similarly defined. Since p(n) is periodic, ¯w(n) = G¯s(n) for}

all n, where G = diag[p(0)IMt, p(1)IMt, · · · , p(N − 1)IMt] ∈ R

MtN ×MtN _{is a diagonal matrix.}

In terms of block signals, (2.3) can be written as ¯

x(n) = H0w(n) + H¯ 1w(n − 1) + ¯¯ v(n) = H0G¯s(n) + H1G¯s(n − 1) + ¯v(n) (2.4)

where H0 is an MrN × MtN block lower-triangular Toeplitz matrix with [H(0)T H(1)T · · ·

H(L)T ₀T

Mr×Mt · · · 0

T Mr×Mt]

T _{∈ C}MrN ×Mt as its first block column (i.e., the first M

t columns),

and H1 is an MrN × MtN block upper-triangular Toeplitz matrix with [0Mr×Mt · · · 0Mr×Mt

H(L) H(L − 1) · · · H(1)] ∈ CMr×MtN _{as its first block row (i.e., the first M}

r rows).

The problem we study in this paper is blind identification of the MIMO channel matrix H = [H(0)T _H(1)T

· · · H(L)T_]T _{using second-order statistics of the received data. The proposed method exploits}

the cyclostationarity induced by periodic precoding at the transmitters. The performance of the proposed identification algorithm (Section 3.3) depends critically on the choice of the pre-coding sequence. We discuss the optimal selection of the sequence that yields the best perfor-mance. We define the following operations that will be used in the derivation of the main result. First, for any m × m matrix A = [ak,l]0≤k,l≤m−1, define Γj(A) = [a0,j a1,j+1 · · · am−1−j,m−1]T

for 0 ≤ j ≤ m − 1, i.e., Γj(A) is the vector formed from the jth super-diagonal of A. Second,

for any Mrn × Mrn matrix B = [Bk,l]0≤k,l≤n−1, where Bk,l is a block matrix of dimension

Mr× Mr, define Υj(B) = [BT_0,j BT_1,j+1 · · · BT_{n−1−j,n−1}]T for 0 ≤ j ≤ n − 1, i.e., Υj(B) is the

matrix formed from the jth block super-diagonal of B.

3

Channel Identification

We study channel identification in this section. In Section 3.1, we derive the proposed method assuming the channel order is known and the noise is absent. We show that by appropriately selecting the periodic precoding sequence, any MIMO channel satisfying (A3) is identifiable up to an Mt × Mt unitary matrix ambiguity. In Section 3.2, we show that

the proposed method is robust with respect to channel order overestimation and we propose an identification algorithm in Section 3.3. The effect of noise and optimal selection of the precoding sequence are discussed in Section 4.

3.1

The Identification Method

We consider the noise free case and assume that the channel order L is known. Equation (2.4) now becomes

¯

x(n) = H0G¯s(n) + H1G¯s(n − 1) (3.1)

With assumption (A1), the covariance matrix of ¯x(n) can be written as R¯x(0) = E[¯x(n)¯x(n)∗] = σs2 H0G2H∗0+ H1G2H∗1

(3.2) Let J ∈ RN ×N _{be the matrix whose first sub-diagonal are all one, i.e., Γ}

1(JT) = [1 1 · · · 1]T ∈

R(N −1)_{, and all remaining entries are zero. The block Toeplitz structures of H0 and H1 allow} us to write H0 =PL_k=0Jk⊗ H(k) and H1 =PL_k=0(JT)N −k⊗ H(k), respectively. Besides, we

define Gp = diag[p(0), p(1), · · · , p(N − 1)] ∈ RN ×N. Hence H0G2H∗0 can be written as

H0G2H∗0 = L X k=0 Jk_{⊗ H(k) G}2p⊗ IMt
XL l=0 Jl_{⊗ H(l)}
∗ = L X k=0 L X l=0 Jk ⊗ H(k)
G2p⊗ IMt
(JT₎l ⊗ H(l)∗
= L X k=0 L X l=0 Jk_G2 p(J T₎l
⊗ (H(k)H(l)∗ ) (3.3)

where we have used the identies (A ⊗ B)∗ _{= A}∗

⊗ B∗

and (A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD) [20, p.190]. Similarly, H1G2H∗1 can be written as

H1GH∗1 = L X k=0 L X l=0 (JT₎N −k_G2 pJ N −l
⊗ (H(k)H(l)∗ ) (3.4)

The following proposition shows that the matrices Jk_G2

p(JT)l and (JT)N −kG2pJN −l have

special structures that allow decomposition of (3.2) into a group of decoupled equations. Roughly speaking, the jth block super-diagonal part of (3.2) involves only the unknown “channel product matrices”, H(k)H(k + j)∗

, k = 0, 1, · · · , L − j. For example, the equa-tions corresponding to the diagonal blocks (j = 0) involve only H(k)H(k)∗

, k = 0, 1, · · · , L. In the proposed identification algorithm, these “channel product matrices” are computed first by solving linear equations, and then the channel impulse response matrices H(k) are com-puted via eigenvalue-eigenvector decomposition.

(a) For l = k + j, where 0 ≤ j ≤ L − k, both Jk_G2

p(JT)l and (JT)N −kG2pJN −l are upper

triangular matrices with only the respective jth upper diagonals nonzero, and Γj JkG2p(J T₎l
= [0 · · · 0_{| {z }} k entries p(0)2 p(1)2 _{· · · p(N − 1 − k − j)}2 | {z } N −k−j entries ]T _(3.5) Γj (JT)N −kG2pJ N −l
= [p(N − k)_| 2 p(N − k + 1)_{z 2 · · · p(N − 1)_}2 k entries 0 · · · 0 | {z } N −k−j entries ]T _(3.6) (b) For l < k, both Γj JkG2p(JT)l
and Γj (JT)N −kG2pJN −l

are lower triangular with zero diagonal matrices.

Proof : See [8].

It follows from (3.5) and (3.6) that Γj JkG2p(J T₎l
_{+ Γ} j (JT)N −kG2pJ N −l
=      [p(N − k)_| 2 · · · p(N − 1)_{{z }}2 k entries p(0)2 _{· · · p(N − 1 − k − j)}2 | {z } N −k−j entries ]T _if j = l − k ≥ 0 0(N −j)×1 if j 6= l − k (3.7) Since Υj JkG2p(JT)l
⊗ H(k)H(l)∗
= Γj JkG2p(JT)l
⊗ H(k)H(l)∗ (3.8) and Υj (JT)N −kG2pJ N −l
⊗ H(k)H(l)∗
= Γj (JT)N −kG2pJ N −l
⊗ H(k)H(l)∗ (3.9) it follows from (3.2)-(3.4) and (3.7)-(3.9) that Υj(R¯x(0)) can be derived as follows.

Υj(R¯x(0)) = Υj H0G2H∗0+ H1G2H∗1
= L X k=0 L X l=0 Υj JkG2p(J T₎l
⊗ (H(k)H(l)∗ )
+ Υj (JT)N −kG2pJ N −l
⊗ (H(k)H(l)∗ )
= L X k=0 L X l=0 {Γj JkG2p(J T )l
+ Γj (JT)N −kG2pJ N −l
} ⊗ H(k)H(l)∗ = L−j X k=0 [p(N − k)2 · · · p(N − 1)2 p(0)2 _{· · · p(N − 1 − k − j)}2]T _{⊗ H(k)H(k + j)}∗ = L−j X k=0 [p(N − k)2IMr · · · p(N − 1) 2_I Mr p(0) 2_I Mr · · · p(N − 1 − k − j) 2_I Mr] T_{H(k)H(k + j)}∗ (3.10) The right hand side of (3.10) is a linear combination of block columns with the channel product

matrices, H(k)H(k + j)∗

, as coefficients. If we define, for 0 ≤ j ≤ L,

then (3.10) can be written in a more compact form (3.12). Υj(R¯x(0)) = MjFj ∀ 0 ≤ j ≤ L (3.12) where Mj ∈ RMr(N −j)×Mr(L−j+1) is defined as Mj =          p(0)2 _{p(N − 1)}2 _{· · · p(N − L + j)}2 p(1)2 _p(0)2 _{· · · p(N − L + j + 1)}2 .. . ... ... ... p(N − 2 − j)2 _{p(N − 3 − j)}2 _{· · · p(N − L − 2)}2 p(N − 1 − j)2 _{p(N − 2 − j)}2 _{· · · p(N − L − 1)}2          ⊗ IMr (3.13)

We note that Mj, 1 ≤ j ≤ L, is obtained from M0 by deleting its last jMr rows and last jMr

columns.

Since N > L + 1, the (L + 1) equations in (3.12) are overdetermined and the equations are also consistent. If Mj is full column rank, then the solution can be obtained as

Fj = (MTjMj)−1MTjΥj(R¯x(0)) (3.14)

If Fj, 0 ≤ j ≤ L, are computed from (3.14), then we have the channel product matrices

H(k)H(l)∗

for 0 ≤ k ≤ l ≤ L. We now consider the computation required to determine the channel impulse response matrices H(0), H(1), · · · , H(L) from Fj.

Let Q be the Hermitian matrix defined by Υj(Q) = Fj for j = 0, 1, · · · , L, and let the

channel matrix H = [H(0)T _H(1)T _{· · · H(L)}T_]T_{. Clearly we have}

Q = HH∗

(3.15) Since rank(H) = Mt by assumption (A3), Q has rank Mt. Since Q is Hermitian and positive

semidefinite, Q has Mt positive eigenvalues, say, λ1, · · · , λMt. We can expand Q as

Q = Mt X j=1 (pλjdj)( p λjdj)∗ (3.16)

where dj is a unit norm eigenvector of Q associated with λj > 0. We can thus choose the

channel matrix to be b H = [pλ1d1 p λ2d2 · · · p λMdMt] ∈ C Mr(L+1)×Mt _(3.17)

We note H can only be identified up to a unitary matrix ambiguity U ∈ CMt×Mt, i.e., bH = HU,

since bH bH∗ _{= HH}∗ _{= Q.}

We note that the matrix Mj, j = 0, 1, · · · , L, is determined by the precoding sequence.

We summarize what we have so far:

(a) If the MIMO channel described by (2.3) satisfies (A1) and (A3) and the channel order L is known, then the channel matrix H can be identified up to a unitary matrix ambiguity. (b) The proposed identification method use the covariance matrix of the received signal Rx¯(0)

as data, and the computations involved are solving linear equations (3.12) and perform-ing eigenvalue-eigenvector decomposition of the Hermitian matrix Q in (3.16).

We note that in the proposed method, the channel condition is assumption (A3), i.e., the channel matrix H is full column rank. Hence the channel needs not be irreducible or column reduced. If Mr ≥ Mt (more outputs), then (A3) is generically satisfied [18, ch.7]. If Mt > Mr

(more inputs), then (A3) is generically satisfied provided (L + 1)Mr ≥ Mt. Besides, in (A1),

for simplicity, we have assumed that the covariance matrix of the source is the identity matrix. The method still applies if the covariance matrix Σ2

s is diagonal. The equation (3.15) becomes

Q = HΣ2 sH

∗ _{and the channel is identifiable up to a unitary matrix ambiguity.}

3.2

Channel Order Overestimation

So far we have assumed that the channel order L is known. If only an upper bound ˆ

L ≥ L is available with N > ˆL + 1, then following the same process given in Section 3.1, the corresponding Mr( ˆL + 1) × Mr( ˆL + 1) matrix Q can be similarly constructed as in (3.15). The

last ( ˆ_{L − L) block columns (i.e., ( ˆL − L)M}r columns) of Q are zero, so are its last ( ˆL − L)

block rows. Hence again, Q is of rank Mt and has Mt positive eigenvalues with the associated

eigenvectors all of the form ˆd = [dT _{0 · · · 0]}T _{∈ C}Mr( ˆL+1) _{where d ∈ C}Mr(L+1)_{. Thus, we can}

determine the actual channel order and impulse response matrices, up to a unitary matrix ambiguity, from the Mt eigenvectors associated with the Mt positive eigenvalues of Q.

3.3

Identification Algorithm

We summarize the proposed method as the following algorithm.

1) Select the precoding sequence p(n) such that each matrix Mj defined in (3.13) is full column

rank.

2) Estimate the autocovariance matrix R¯x(0) via the time average

ˆ Rx¯(0) = 1 K K X i=1 ¯ x(i)¯x(i)∗ (3.18)

where K is the number of data block (i.e., KN is the number of samples for each transmitter). 3) Compute Fj, formed by the channel product matrices, for j = 0, 1, · · · , L, using (3.14).

4) Form the matrix Q as in (3.15), and obtain the channel impulse response (3.17) by com-puting the Mt largest eigenvalues and the associated eigenvectors of Q.

4

Optimal Selection of the Precoding Sequence

In Section 3, we see that in order to identify the channel, the precoding sequence must be selected so that the resulting matrix Mj is full column rank such that Fj can be exactly

solved as (3.14). However, when noise is present, the covariance matrix ˆRx¯(0) contains the

contribution of noise and numerical error is present in the estimation of ˆRx¯(0) by (3.18). This

implies that (3.12) usually has no solution and (3.14) becomes a least squares approximate solution. The choice of Mj will affect error in the computation of Fj since different MTjMj in

(3.14) may have different condition numbers. In this section, we discuss the optimal selection of the precoding sequence, which takes into account the effect of noise and numerical error in estimating ˆR¯x(0), so as to increase the accuracy of Fj and thus reduce the channel estimation

error.

4.1

Optimality Criterion

Now we consider the general case that the noise is present and discuss the design of the precoding sequence p(n). From (2.4) and assumption (A1), the covariance matrix of the received signal is

R¯x(0) = H0G2H∗0+ H1G2H∗1+ σv2IMr⊗ IN (4.1)

From (4.1) and (3.2), we see that noise has only contribution to the diagonal entries of R¯x(0).

Therefore the (L + 1) decoupled groups of equations in (3.12) remain unchanged, except for the j = 0 group, which becomes

Υ0(Rx¯(0)) = Υ0 H0G2H∗0+ H1G2H∗1
+ σv2Υ0(IMr ⊗ IN) = M0F0+ Y (4.2) where Y = σ2 v[IMr IMr · · · IMr] T ∈ RMrN ×Mr _{with unknown σ}2

v. Thus from (3.14), ˆF0, the

least squares approximation of F0, can be written by

ˆ

F0 = (MT0M0)−1MT0 (M0F0+ Y)

| {z }

Υ0(R¯x(0))

= F0+ (MT0M0)−1MT0Y = F0 + Z (4.3)

which is F0plus a perturbation term due to noise. The perturbation term Z is the least squares

solution of the equation M0Z = Y. We note that if every column of Y is orthogonal to every

column of M0, then Z = 0, which implies ˆF0 = F0. But that is impossible since the entries of

M0 are positive and those of Y are nonnegative. Therefore, we seek to appropriately choose

that of M0 as possible. To this end, we first define qki and yi shown below as the columns of M0 and Y, respectively: M0 = " q01 q02 · · · q0Mr | {z } M0(:,1:Mr) q11 q12 · · · q1Mr | {z } M0(:,Mr+1:2Mr) · · · qL1 qL2 · · · qLMr | {z } M0(:,LMr+1:(L+1)Mr) # (4.4) Y = σ_v2[IMr IMr · · · IMr] T _{= [y} 1 y2 · · · yMr] (4.5)

Then, due to the special structure of the block matrix M0 and Y, it is easy to check that qki

is orthogonal to yj, i.e., qTkiyj = 0 for j 6= i, e.g.,

qT 01y2 = [p(0)20 · · · 0 | {z } Mr entries · · · p(N − 1)_{| {z}20 · · · 0_} Mr entries ][0 σ_v2_{0 · · · 0} | {z } Mrentries · · · 0 σ_| v2_{z0 · · · 0_} Mrentries ]T _{= 0} and each qT

kiyi assumes the same value, σv2

PN −1 n=0 p(n)2, for k = 0, 1, · · · , L, i = 1, 2, · · · , Mr, e.g., qT 01y1 = [p(0)20 · · · 0 | {z } Mr entries · · · p(N − 1)_{| {z}20 · · · 0_} Mr entries ][σ2_{v 0 · · · 0} | {z } Mrentries · · · σ_{| {z }}v20 · · · 0 Mrentries ]T _{= σ}2 v N −1_X n=0 p(n)2

Thus we only need to consider the relation between columns of q01 and y1 (the case of k = 0

and i = 1). Define the correlation coefficient γ = q

T 01y1

kq01k2ky1k2

(4.6) Since γ is nonnegative and by Cauchy-Schwarz inequality, 0 ≤ γ ≤ 1. In order to make the perturbation term Z small, we choose q01 so that the correlation coefficient γ is as small

as possible. Based on this point of view, we formulate the optimal selection problem as minimizing γ subject to 1 N N −1 X n=0 |p(n)|2 = 1 (4.7) |p(n)|2 ≥ τ > 0, ∀ 0 ≤ n ≤ N − 1 (4.8) Roughly, constraint (4.7) normalizes the power gain of the precoding sequence of each trans-mitter to 1; constraint (4.8) requires that at each instant, the power gain is no less than τ . Note that the problem of selecting the precoding sequence is identical to the SISO case considered in [8]. Thus the optimal precoding sequence p(n) is a two-level sequence with a single peak in one period [8]. More specifically, for each m, 0 ≤ m ≤ N − 1,

p(n) = ( p

N (1 − τ) + τ , n = m √

is an optimal precoding sequence. Because the precoding sequence is periodic with period N , the single peak can be placed at any one of the N positions which yield the same γ. However the peak location m does significantly affect the numerical condition of the linear equation (3.12) as we discuss next.

4.2

On Selection of m

We now consider the selection of m. We know that different choices of m result in different matrix Mj and affect the numerical computation of Fj, j = 1, 2, · · · , L, in (3.14) and ˆF0 in

(4.3), since different MT

jMj may have different condition number. If the condition number is

large, then the matrix MT

jMj is ill-conditioned and the computation in (3.14) and (4.3) are

sensitive to data error. Let

µ = max

0≤j≤Lκ(M T

jMj) (4.10)

where κ(A) is the condition number of A. Our goal is to choose m so as to minimize the largest condition number of the corresponding matrices MT

jMj, j = 0, 1, · · · , L. Since the

peak appears at one of the N possible positions in the periodic precoding sequence, there are N precoding sequences which may result in N different µ. The following result shows that some choices of m are to be avoided since they result in some Mj being rank deficient and

thus µ = ∞.

Proposition 4.1 : At least one Mj, 0 ≤ j ≤ L, is not full column rank if and only if

N − L + 1 ≤ m ≤ N − 2. Proof : See Appendix A.

Hence if we choose, either 0 ≤ m ≤ N − L or m = N − 1, then each Mj is full column rank

and the channel is identifiable. The following result shows that we can classify the remaining choices into 2 groups that are relevant to the optimal choice of m.

Proposition 4.2 :

(a) Each of the (N − L) choices, m = 0, m = 1, · · · , m = N − L − 1, results in the same µ denoted by µ1.

(b) The two choices m = N − L and m = N − 1 result in the same µ denoted by µ2. Also

µ2 ≥ µ1.

Proof : See Appendix A.

From Proposition 4.2, we know if µ2 > µ1, then we choose case (a); if µ2 = µ1, we

proceed to compare the second largest condition numbers of the set of matrices {MT

jMj}Lj=0

for these two cases and choose the case whose value is smaller. If they are again equal, the same procedure can be done by comparing the third largest condition numbers and so on. Moreover, for 0 ≤ m ≤ N − L − 1 (case (a)), since the condition numbers of MT

jMj are the

(a). Similarly, m = N − 1 can be used to represent case (b). Hence the optimal selection of m reduces to one of two cases: m = 0 or m = N − 1. In other words, the optimal precoding sequence has a peak either at the beginning or at the end.

5

Simulation Results

In this section, we use several examples to demonstrate the performance of the proposed method. The channel normalized root-mean-square error (NRMSE) is defined as

NRMSE = 1 kHkF v u u t1 I I X i=1 k bH(i)_U(i)_{− Hk}2 F (5.1)

where the unitary matrix ambiguity U(i)_{is computed by the least squares method [15], solving}

min

U(i) k bH

(i)_U(i)

− Hk2 (5.2)

I is the number of Monte Carlo runs, and bH(i) _{= [ bH}(i)₍₀₎T _Hb(i)₍₁₎T

· · · bH(i)_(L)T_]T _{is the}

estimate of channel impulse response matrix H. The input source symbols are i.i.d. QPSK signals. The channel noises are white Gaussian. The signal-to-noise ratio (SNR) at the output is defined as SNR = 1 N PN −1 n=0 E[kz(n)k22] E[kv(n)k2 2] (5.3) where z(n) = [z1(n) · · · zMr(n)]

T _{is the signal component of the received signal (see Figure}

1) Simulation 1 – optimal selection of precoding sequences In this simulation, we use the following model

H(z) = " 0.4851 0.3200 −0.3676 0.2182 # | {z } H(0) + " −0.4851 0.9387 0.8823 0.8729 # | {z } H(1) z−1+ " 0.7276 −0.1280 0.2941 −0.4364 # | {z } H(2) z−2 (5.4)

given in [11] to demonstrate the effect of different precoding sequences on the performance of the proposed method. For comparsion, the first sequence is chosen as {0.767 1.07 1.07 1.07}, which satisfies (4.7) and (4.8). The second and third sequences are chosen based on (4.9) for N = 4 and τ = 0.5878 with the two possible peak positions: m = 0 and m = 3. By computation, the corresponding µ for the three cases are 40.0, 4.66 and 22.1, respectively. Thus m = 0 is the optimal selection. Figure 5.1 shows that for SNR=10 dB, there are about 5∼6 dB and 6∼7 dB difference in NRMSE between the optimal one and two others. Figure 5.2 shows the NRMSE versus SNR when the number of samples (for each transmitter) is fixed at 1000. For each sequence, the NRMSE decreases as SNR increases and is roughly constant for SNR ≥ 20 dB. 100 200 300 400 500 600 700 800 900 1000 −24 −22 −20 −18 −16 −14 −12 −10 −8 −6 Number of Samples Channel NRMSE(dB) m=0 (optimal) m=3 non−optimal p(n)

Figure 5.1: Channel NRMSE versus Number of Samples

2) Simulation 2 – robustness to channel order overestimation

In this simulation, we use the same channel model (5.4). For each upper bound ˆ_{L, 0 ≤} ( ˆ_{L − L) ≤ 6, we choose N = ˆ}L + 2, SNR=10 dB, and 1000 samples (for each transmitter)

0 5 10 15 20 25 30 35 40 −24 −22 −20 −18 −16 −14 −12 −10 −8 SNR(dB) Channel NRMSE(dB) m=0 (optimal) m=3 non−optimal p(n)

Figure 5.2: Channel NRMSE versus Output SNR

for simulation. The precoding sequences are chosen as (4.9) with m = 0 and τ = 0.5878. Figure 5.3 shows the NRMSE increases with increasing channel order overestimation. With the channel order fixed at L = 2, the NRMSE increases from -22.5dB to -16dB as the ( ˆ_{L − L)} increases from 0 to 6. The proposed method is quite robust to channel order overestimation since the NRMSE still maintains a low value (about -16dB) when ( ˆ_{L − L) = 6.}

3) Simulation 3 – a 3-input 2-output channel

In this simulation, we use the 3-input 2-output model

H(z) = " 1.6 0.88 0.66 0.8 0.44 0.33 # | {z } H(0) + " −0.44 0.35 0.14 −0.14 0.37 0.23 # | {z } H(1) z−1+ " 0.13 0.01 0.08 0.26 0.02 0.16 # | {z } H(2) z−2 (5.5)

to illustrate the performance of the proposed method for channel with more inputs than outputs. Note that this model is not irreducible [15] because H(0) is not full rank, and it is not column reduced [15] either because H(2) is not full rank. It satisfies the condition (L + 1)Mr ≥ Mt given at the end in Section 3.1 and H is full column rank. The precoding

sequences (N = 4) are given as in Simulation 1 for m = 0 and m = 3. Figures 5.4 shows NRMSE versus the number of data samples with SNR=10 dB. Figure 5.5 shows NRMSE versus SNR with the number of data samples (for each transmitter) fixed at 1000.

0 1 2 3 4 5 6 −23 −22 −21 −20 −19 −18 −17 −16 −15

Overestimated Channel Order

Channel NRMSE(dB)

Figure 5.3: Channel NRMSE versus ( ˆ_{L − L)}

100 200 300 400 500 600 700 800 900 1000 −22 −21 −20 −19 −18 −17 −16 −15 −14 −13 −12 Number of Samples Channel NRMSE(dB) m=0 m=3

0 5 10 15 20 25 30 35 40 −26 −24 −22 −20 −18 −16 −14 −12 SNR(dB) Channel NRMSE(dB) m=0 m=3

Figure 5.5: 3-input 2-output Model: Channel NRMSE versus Output SNR

In this simulation, we use the 2-input 3-output model

H(z) =     1 0 0 1 1 1     | {z } H(0) +     −0.6 −0.5 0 0 −1.2 −1     | {z } H(1) z−1 (5.6)

given in [15]. We compare the performance of the proposed method with the outer-product-decomposition-algorithm (OPDA) method [14, 15]. The SNR is fixed at 10 dB. Note that this model is irreducible because H(0) is full rank, but not column reduced because H(1) is not full rank. The precoding sequences (N = 4) are given as in Simulation 1 for m = 0 and m = 3. As seen in Figure 5.6, with the optimal precoding sequence, the proposed method yields lower NRMSE than that of OPDA method, although with the selection of m = 3 (peak at the end), the resulting NRMSE is about 3dB higher than that of OPDA method.

5) Simulation 5 – comparison with the method in [10]

In this simulation, the performance of the proposed method is compared with that of [10], which also use periodic precoding to induce cyclostationarity at the transmitter. We use the

100 200 300 400 500 600 700 800 900 1000 −26 −24 −22 −20 −18 −16 −14 −12 −10 Number of Samples Channel NRMSE(dB) m=0 (optimal) m=3 OPDA method

Figure 5.6: 2-input 3-output Model: Channel NRMSE versus Output SNR

channel (after normalized)

H(z) = " 0.4082 0.5392 0.6396 0.4264 # | {z } H(0) + " −0.4082 0.5392 −0.4264 −0.6396 # | {z } H(1) z−1+ " 0.8164 0.6470 0.6396 0.6396 # | {z } H(2) z−2 (5.7)

given in [10]. For the proposed method, we use three precoding sequences. Precoding se-quences 1 and 2 are those given in Simulation 1 for m = 0 and m = 3 (N = 4), respectively. Precoding sequence 3 is {1.1 1.1 1.1 1 1 1}, which is the same given in [10]. We use 1200 i.i.d. 4-QAM symbols, 2000 Monte Carlo trials, and the same definitions of SNR and MSE given in [10, 11]. Figure 5.7 shows that the proposed method performs better than the method in [10].

6

Conclusions

We propose a method for blind identification of FIR MIMO channels using periodic pre-coding sequence. Since the cyclostationarity is induced at the transmitter, the identifiability condition imposed on the channel is minimum: it only requires that channel impulse response matrix H = [H(0)T _H(1)T

· · · H(L)T_]T _{is full column rank. The channel transfer matrix is}

0 5 10 15 20 25 30 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 SNR Channel MSE sequence 1 sequence 2 sequence 3 method in [10]

Figure 5.7: Proposed Method versus the method in [10]

inputs. The method is shown to be robust with respect to channel order overestimation. The performance of the algorithm depends on the precoding sequence which is optimally designed to reduce the effect of noise and error in estimating the covariance matrix of the received data. Simulation results show that the method yields good performance.

Appendix

A

Proof of Proposition 4.1 and 4.2

Preliminary :

For each j, let Nj ∈ R(N −j)×(L−j+1) be similarly defined as (3.13), except that IMr is

re-placed by 1. It can be easily check that there exists permutation matrices Plj ∈ R

Mr(N −j)×Mr(N −j)

and Prj ∈ R

Mr(L−j+1)×Mr(L−j+1)such that P

ljMjPrj = diag[Nj, Nj, · · · , Nj] = Dj ∈ R

Mr(N −j)×Mr(L−j+1)

is a block diagonal matrix with each block of dimension (N −j)×(L−j+1). Since Plj

T = Plj −1 and Prj T = Prj −1 _{[19, p.110], we have M} j = Plj T DjPrj T

. Hence Mj is full column rank if

and only if Nj is full column rank for j = 0, 1, · · · , L.

Also, MT jMj = (PrjD T jPlj)(Plj T_D jPrj T_{) = P} rjD T jDjPrj T _{= P} rjdiag[N T jNj, · · · , NTjNj]Prj T_.

Let λ(A) denote the spectrum of A [19, p.310], that is, the set of eigenvalues of A. Then λ(MT

Proof of Proposition 4.1 :

If at N − L + 1 ≤ m ≤ N − 2, it can be checked that Nj, j = 2, 3, · · · , L − 1 is not of full

column rank since it has two columns both equal to [τ τ · · · τ]T _{which implies that at least}

one Mj is rank deficient and vice versa.

Proof of Proposition 4.2 :

From the Preliminary, since λ(MT

jMj) = λ(NTjNj), the condition number of MTjMj

is identical to that of NT

jNj, i.e., κ(MTjMj) = κ(NTjNj). Thus we need only compute the

condition number of NT jNj.

Case (a): For m = 0, m = 1, · · · , and m = N − L − 1, we know NT

jNj = a · IL−j+1+ (2b + cj) · [1 · · · 1]T[1 · · · 1] (A.1)

where a = N2_{(1 − τ)}2_{, b = N τ (1 − τ), c}

j = (N − j)τ2. Hence the maximum and minimum

eigenvalues are a+(L−j +1)(2b+cj) and a respectively. Thus the condition number of MTjMj

is 1 + [(L − j + 1)(2b + cj)/a] which is a decreasing function of j. Therefore the corresponding

µ is equal to µ1 = 1 + [(L + 1)(2b + c0)/a].

Case (b): For m = N − L and m = N − 1, we consider the j = 0 case and j 6= 0 case for Nj separately. For j = 0 with m = N − L or m = N − 1, direct multiplication of NT0N0

gives the same matrix as (A.1), and the condition number of MT

0M0 is µ1. For j 6= 0 with m = N − L, direct multiplication of NT jNj yields NT jNj =          a + 2b + cj 2b + cj 2b + cj · · · 2b + cj b + cj 2b + cj a + 2b + cj 2b + cj · · · 2b + cj b + cj ... ... . .. ... ... ... 2b + cj 2b + cj 2b + cj · · · a + 2b + cj b + cj b + cj b + cj b + cj · · · b + cj cj          ∈ R(L−j+1)×(L−j+1) (A.2) The eigenvalues of NT

jNj in ascending order, are αj, a, βj, where a has a multiplicity L−j −1,

and βj = 1₂{(L − j)(2b + cj) + (a + cj) + p [(L − j)(2b + cj) + (a − cj)]2+ 4(L − j)(b + cj)2}, αj = 1₂{(L − j)(2b + cj) + (a + cj) − p [(L − j)(2b + cj) + (a − cj)]2+ 4(L − j)(b + cj)2}. All

of the eigenvalues are positive and real. (A proof is given in Appendix B). It can be similarly shown that for j 6= 0 with m = N − 1, NT

jNj has the same eigenvalues αj, a, βj. Hence for

j = 1, 2, · · · , L, λ(MT

jMj) = {αj, a, βj} and the condition number is

κ(MT jMj) = βj αj = 1 + χ2 j − 4(N − L)b2+ χj q χ2 j − 4(N − L)b2 2(N − L)b2 (A.3)

maximum value is β1/α1. Therefore, combining the two cases (j = 0, j 6= 0), the corresponding

µ is µ2 = max{µ1, β1/α1} ≥ µ1.

B

The Eigenvalues of N

N

for m = N − L

Proof :

Let Aj = NTjNj defined in (A.2), then Aj is positive definite since Nj is full column rank.

It can be checked that the eigenvectors corresponding to (L − j − 1) multiple eigenvalue a are: [1, −1, 0, 0, · · · , 0]T

, [1, 1, −2, 0, · · · , 0]T

, · · · , [1, 1, · · · , 1, −(L−j −1), 0]T_{. The remaining}

eigenvectors are [1, 1, · · · , 1, x]T _{∈ R}L−j+1_{. Hence}

Aj        1 .. . 1 x        =        a + (L − j)(2b + cj) + (b + cj)x .. . a + (L − j)(2b + cj) + (b + cj)x (L − j)(b + cj) + cjx        = λj        1 .. . 1 x        (B.1)

which implies the following two equations

a + (L − j)(2b + cj) + (b + cj)x = λj (B.2)

(L − j)(b + cj) + cjx = λjx (B.3)

Substitute (B.2) into (B.3), we can get an second order equation of x. Solving this equation can lead to two solutions of x. Bring these two x into (B.2) and we can obtain the two eigenvalues βj, αj. In addition, βj ≥ a because of (B.4)

βj = 1 2{(L − j)(2b + cj) + (a + cj) + q [(L − j)(2b + cj) + a − cj]2+ 4(L − j)(b + cj)2} ≥ 1 2{(L − j)(2b + cj) + (a + cj) + q [(L − j)(2b + cj) + a − cj]2} = 1 2{[(L − j)(2b + cj) + (a + cj) + [(L − j)(2b + cj) + a − cj]} = a + (L − j)(2b + cj) ≥ a (B.4)

and αj ≤ a because of the interlacing property [19, p.396].

References

[1] M. K. Tsatsanis and G. B. Giannakis, ”Transmitter induced cyclostationarity for blind channel equalization”, IEEE Trans. Signal Processing, vol. 45, no. 7, pp. 1785-1794, July 1997.

[2] E. Serpedin and G. B. Giannakis, ”Blind channel identification and equalization with modulation induced cyclostationarity”, IEEE Trans. Signal Processing, vol. 46, no. 7, pp. 1930-1944, July 1998.

[3] A. Scaglione, G. B. Giannakis, and S. Barbarossa, ”Redundant filter bank precoders and equalizers Part I: Unification and optimal designs.”, IEEE Trans. Signal Processing, vol. 47, no. 7, pp. 1988-2006, July 1999.

[4] A. Scaglione, G. B. Giannakis, and S. Barbarossa, ”Redundant filter bank precoders and equalizers Part II: Blind channel estimation, synchronization, and direct equalization”, IEEE Trans. Signal Processing, vol. 47, no. 7, pp. 2007-2022, July 1999.

[5] A. Chevreuil, E. Serpedin, P. Loubaton and G. B. Giannakis, ”Blind channel identifica-tion and equalizaidentifica-tion using periodic modulaidentifica-tion precoders: performanve analysis”, IEEE Trans. Signal Processing, vol. 48, no. 6, pp. 1570-1586, June 2000.

[6] E. Serpedin, A. Chevreuil, G. B. Giannakis, and P. Loubaton, ”Blind channel and carrier frequency offset estimation using periodic modulation precoders”, IEEE Trans. Signal Processing, vol. 48, no. 8, pp. 2389-2405, Aug. 2000.

[7] P. Ciblat, A. Chevreuil and P. Loubaton, ”α-repetition/modulation and blind second order identification”, IEEE Trans. Signal Processing, vol. 48, no. 11, pp. 3153-3161, Nov 2000.

[8] C. A. Lin and J. W. Wu, ”Blind identification with periodic modulation: A time-domain approach”, IEEE Trans. Signal Processing, vol. 50, no. 11 pp. 2875-2888, Nov. 2002. [9] A. Chevreuil and P. Loubaton, ”MIMO blind second-Order equalization method and

conjugate cyclostationarity”, IEEE Trans. Signal Processing, vol. 47, no. 2, pp. 572-578, Feb. 1999.

[10] H. B¨olcskei, R. W. Heath, Jr., and A. J. Paulraj, ”Blind channel estimation in spatial multiplexing systems using nonredundant antenna precoding”, Proc. of the 33rd Annual IEEE Asilomar Conf. on Signals, Systems, and Computers, vol. 2, pp. 1127-1132, Oct. 1999.

[11] H. B¨olcskei, R. W. Heath, Jr., and A. J. Paulraj, ”Blind channel identification and equalization in OFDM-based multiantenna systems”, IEEE Trans. Signal Processing, vol. 50, no. 1, pp. 96-109, Jan. 2002.

[12] I. Bradaric, A. P. Petropulu, and K. I. Diamantaras, ”On blind identifiability of FIR-MIMO systems with cyclostationary inputs using second order statistics”, IEEE Trans. Signal Processing, vol. 51, no. 2, pp. 434-440, Feb. 2003.

[13] I. Bradaric, A. P. Petropulu, and K. I. Diamantaras, ”Blind MIMO FIR channel identifi-cation based on second-order spectra correlations”, IEEE Trans. Signal Processing, vol. 51, no. 6, pp. 1668-1674, June 2003.

[14] Z. Ding, ”Matrix outer-product decomposition method for blind multiple channel identi-fication”, IEEE Trans. Signal Processing, vol. 45, no. 12, pp. 3053-3061, Dec. 1997. [15] Z. Ding and L. Qiu, ”Blind MIMO channel identification from second order statistics

using rank deficient channel convolution matrix”, IEEE Trans. Signal Processing, vol. 51, no. 2, pp. 535-544, Feb. 2003.

[16] E. Moulines, P. Duhamel, J-F. Cardoso, and S. Mayrargue, ”Subspace method for the blind identification of multichannel FIR filters”, IEEE Trans. Signal Processing, vol. 43, no. 2, pp. 516-525, Feb. 1995.

[17] C. A. Lin and C. W. King, ”Minimal periodic realization of transfer matrices”, IEEE Trans. Automatic Control, vol. 38, no. 3 pp. 462-466, Mar. 1993.

[18] M. W. Hirsch and S. Smale, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press, 1974.

[19] G. H. Golub and C. F. Van Loan, Matrix Computions, 3rd edition, The Johns Hopkins University Press, 1996.