Weighted minimum noise subspace method for blind system identification1

Signal Processing 69 (1998) 15—27

Weighted minimum noise subspace method for blind system identification1 Karim Abed-Meraim, Yingbo Hua* Department of Electrical and Electronic Engineering, The University of Melbourne, Parkville, Victoria 3052, Australia Received 4 July 1996; received in revised form 18 May 1998

Abstract Blind identification of single-input multiple-outputs FIR system has attracted much attention during the past few years. This interest is due to the fact that accurate system parameter estimates can be obtained using only the second-order statistics of the system outputs. In this paper, a weighted minimum noise subspace W-MNS estimate of the FIR system parameters is developed. The asymptotic variances of the W-MNS system parameter estimates are derived, and a statistically optimal weighting matrix is proposed. The optimal W-MNS estimates of the unknown system parameters are shown to be significantly more accurate than the existing MNS estimates. ( 1998 Elsevier Science B.V. All rights reserved. Zusammenfassung Die blinde Identifikation von single-input multiple-output FIR Systemen hat wa¨hrend der letzten Jahre breite Beachtung erlangt. Dieses Interesse liegt darin begru¨ndet, da{ exakte Scha¨tzungen von Systemparametern bereits mit Hilfe der statistischen Eigenschaften zweiter Ordnung der Ausgaben des Systems gewonnen werden ko¨nnen. In diesem Artikel wird eine gewichtete, minimale Sto¨runterraum- (weighted minimum noise subspace, W-MNS) Scha¨tzung der FIR Systemparameter entwickelt. Es werden die asymptotischen Varianzen der W-MNS Systemparameter hergeleitet und eine im statistischen Sinne optimale Gewichtungsmatrix vorgeschlagen. Es wird gezeigt, da{ die optimalen W-MNS Scha¨tzungen der unbekannten Systemparameter deutlich genauer sind als die bekannten MNS Scha¨tzungen. ( 1998 Elsevier Science B.V. All rights reserved. Re´sume´ L’identification aveugle de syste`mes FIR entre´e-simple sorties-multiples a attire´ beaucoup d’attention depuis quelques anne´es. Cet inte´reˆt est duˆ au fait que des estime´es pre´cises des parame`tres du syste`mes peuvent eˆtre obtenues a` l’aide des statistiques de deuxie`me ordre seules des sorties. Nous de´veloppons dans cet article une estimation W-MNS de type sous-espace de bruit minimum ponde´re´ des parame`tres d’un syste`me. Les variances asymptotiques des estime´es W-MNS

* Corresponding author. Tel.: 0061 3 9344; fax: 0061 3 9344 6678; e-mail: [email protected]. 1 This work has been supported by the Australian Research Council and the Australian Cooperative Research Center for Sensor Signal and Information Processing. 0165-1684/98/$ — see front matter ( 1998 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 8 ) 0 0 0 8 4 - X

16

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15–27

sont de´rive´es, et une matrice de ponde´ration statistiquement optimale est propose´e. Les estime´es W-MNS optimales des parame`tres inconnus du syste`me sont montre´s eˆtre significativement plus pre´cises que les estime´es MNS existantes. ( 1998 Elsevier Science B.V. All rights reserved. Keywords: System identification; Minimum noise subspace; Weighting matrix; Performance analysis

1. Introduction In many data communication systems such as narrowband TDMA digital cellular radio systems, data signals are often transmitted through unknown channels which may introduce severe linear distortion due to time varying multi-paths. Therefore, reliable communication often requires the identification of the channel impulse response. Such identification can facilitate channel equalization as well as maximum likelihood sequence detection. Traditionally, channel identification is performed using training sequences. However, there are a variety of situations where the available channel input training sequence may be too short or even non-existent for channel identification. Blind channel identification (i.e., determining the channel response based only on the channel outputs, without the use of a training sequence) can play useful roles in these situations. Earlier approaches to blind channel identification exploit the higher order statistics of the output (see [9,18,22] for example). Such approaches require generally a large number of data samples which limits their applicability in many practical settings. The recently proposed methods by Gardner [8] and Tong et al. [20] allow blind channel identification using only second order statistics. Since these early contributions, numerous alternative second order statistics based methods have been proposed, e.g. [4,11,14,15,21]. This paper deals with a particular family of blind estimation techniques, referred to as subspace estimation methods [12,15]. The minimum noise subspace (MNS) method [12] is a computationally fast version of the subspace (SS) method [2,15]. We focus here on the MNS method and introduce a class of MNS estimates which depend on a weighting matrix. The MNS method is a special

member in this class. The asymptotic performance of the weighted minimum noise subspace (W-MNS) method is derived as a function of the weighting matrix. Then, the statistically optimal weighting matrix is derived which minimizes the asymptotic variance of estimation error. The rest of this paper is organized as follows. Following the problem formulation in Section 2, we first review in Section 3 the original MNS method [12]. Then, in Section 4, we introduce a class of W-MNS methods, and derive the explicit expression of the optimal weighting-matrix value in terms of the channel parameters. In Section 5, the practical implementation details of the proposed algorithm are given. The last section is devoted to numerical simulations which lend support to the theoretical findings.

2. Problem statement 2.1. Notation AT A* range (A)M Aj A ij I 0 diag(x) deg( f (z))

Transpose of A Conjugate transpose of A Orthogonal complement subspace to the column range of A Moore—Penrose pseudo-inverse of A (i, j)th component of A Identity matrix of appropriate dimension Zero matrix of appropriate dimension Diagonal matrix formed with the elements of the vector x Lowest degree M of the polynomial M $%& + expansion f (z)" f (k)z~k k/0

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15—27

17

E( ) ) Statistical expectation operator Re (x), Im (x) Real part and imaginary part of x, respectively

Also in Eq. (3), N (n) is the vector of the observation i noise; and s(n) the vector of the input samples, i.e., s(n)"[s(n),2,s(n!N!M#1)]. The q outputs of the system can then be put into

2.2. Data model

Y(n)"T (h)s(n)#N(n), N

Consider a discrete multichannel system with FIR impulse responses Mh ( ) )N driven by a single i input s( ) ). The outputs of the system are given by

where

M y (n)" + h ( j)s(n!j)#n (n), i"1,2,q, (1) i i i j/0 where q*2 is the number of channels, M the maximum degree of the q channel polynomial filters, i.e., M"max deg(h (z)), where h (z)" i i i +Mi h ( j)z~j and Mn (n), i"1,2,qN the additive j/0 i i noise assumed to be temporally and spatially white, i.e., E(n (t)n*(s))"d d p2. i j i,j t,s In the context of communications, such a model can be formulated either by using q physical receivers (antenna array) where each FIR filter models the propagation channel from the transmitter to each individual receiver, or by using q virtual receivers obtained by temporally oversampling the channel outputs at q times the baud rate (see [15,20] for more details). In this paper we study the estimation of the channels Mh ( ) )N from the outputs My ( ) )N assumi i ing the polynomial transfer function h(z)" [h (z),2,h (z)]T satisfies the following condition: 1 q h(z)O0 for each z, (2) or in other words, the scalar polynomials Mh (z), i"1,2,qN have no common zeros. i For later reference, we consider a vector of N successive samples from the ith output of the system: Y (n)"T (h )s(n)#N (n), (3) i N i i where Y (n)"[y (n),2,y (n!N#1)]T; T (h ) is an i i i N i N](N#M) block Sylvester matrix,

C

h (0) i T (h )" F N i 0

2 } 2

h (M) i

2

0

}

F

h (0) i

2

h (M) i

D

. (4)

C D

Y (n) 1 Y (n) Y(n)" 2 , F

(5)

C D

T (h ) N 1 T (h ) T (h)" N 2 N F

Y (n) q

(6)

T (h ) N q

and N(n) is the system noise vector. T (h) is known N as a generalized Sylvester (GS) matrix of order N associated with h(z).

3. The MNS method Let R be the covariance matrix of the qN-dimensional vector Y(n): $%& E(Y(n)Y *(n))"T (h)R T *(h)#p2I , R" qN N s N

(7)

$%& E(s(n)s*(n)) and p2 is the additive noise where R " s power. In the sequel, we assume R '0 (positive s definite). It was shown in [5] that, for N'M, the GS matrix has the full-column rank N#M. Thus, the eigen-decomposition of the covariance matrix R can be written as R"U K U *#p2U U *, s s s n n

(8)

where U "[u ,2,u ]: the signal eigenvectors, s 1 M`N U "[u , ,u ]: the noise eigenvectors, n M`N`1 2 qN K "diag(j ,2,j ): the signal eigenvalues, s 1 M`N with j *2*j 'p2. 1 M`N

18

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15–27

Fig. 1. Illustration of the estimation procedure of the q!1 noise vectors.

It is easy to show that range(U )"range(T (h)) and s N range(U )"range(T (h))M. The range space of U is n N n referred to as the noise subspace. The standard subspace (SS) method [15] determines the transfer function h(z) by solving the orthogonality relation U *T (h)"0 n N

(9)

in a least squares sense. This estimate is uniquely (up to a constant scalar) equal to h under the assumption (2). The SS method requires the knowledge of the whole noise subspace. Hua et al. [12] recently show that for a q-channel FIR system only q!1 properly chosen noise vectors are necessary and sufficient. Each of the q!1 noise vectors can be found by computing the least eigenvector of a covariance matrix corresponding to a distinct pair of channels (see Fig. 1). More precisely, as shown in [12], the MNS estimation method proceeds as follows: f We first select q!1 distinct pairs from the q channel outputs My (n), i"1,2,qN. The q!1 i pairs must span a tree which connects all q channel outputs. The channel outputs are the nodes of the tree as shown in Fig. 2.

Fig. 2. Illustration of a tree connecting q"5 channel outputs.

f For each pair of channel outputs i"(i , i ), we 1 2 compute the covariance matrix 1 T~N`1 R(i)" + Y (n)Y *(n), i i ¹!N#1 n/1

(10)

$%& [Y T(n),YT ¹ being the sample size and Y (n)" i2 i i1 (n)]T.Then, we compute its least dominant eigenvector € . 6 i T ,€T ]T where each subvector has the f Let € "[€ i1 6 i2 i 6 dimension 6 N]1 (i.e., € "[€ (0),2,€ (N!1)]T, 6 ik 6 ik 6 ik qN]1 for k"1,2). Then define ‘zero-padded’

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15—27

vectors

C D

€ (1) i €" F i € (q) i

G

€ 6 i1 , where € (k)" € i 6 i2 0

k"i k"i

1

2 otherwise

19

weighted subspace fitting (WSF) method [23] in the narrow-band DOA estimation context, this estimation procedure can be improved using weighted least-squares criterion. The existence and expression of the optimal weighting matrix are given next.

(11) and form a qN](q!1) matrix V of the q!1 n vectors M€ N, i.e., V "[€ ,2,€ ]. i n 1 q~1 $%& f Estimate the channel parameter vector h" [hT,2,hT]T where h "[h (0),2,h (M)]T by q i i i 1 minimizing the least squares criterion hK "arg min ET *( f )V E2 N n f

(12)

under the constraint E f E"1. The significant computational advantage of the MNS method over the SS method is obvious. In particular, the SS method requires a full eigendecomposition of a qN]qN matrix, but the MNS method computes the single least dominant eigenvector of a 2N]2N matrix in parallel for each of q!1 pairs of channel outputs. In [1,12], it has been established that under the data model assumptions and for N'M, (i) the MNS method yields the unique estimate of h, and (ii) q!1 is the smallest number of vectors from the noise subspace in order for an equation like V *T ( f )"0 to yield n N the unique estimate of h. More precisely, we have the following two lemmas [12]. Lemma 1. If there are d qN]1 vectors M€ for i i"1,2,dN such that M€*T (h)"0 for i"1,2,dN i N where T (h) is the qN](M#N) GS matrix assoN ciated with h(z), then h is (possibly) unique up to a constant scalar only if d*q!1. Lemma 2. ¹he MNS method yields the unique (up to a constant scalar) estimate of the channel responses under the data model assumptions. In the above least-squares procedure Eq. (12), all the available relations contribute to the global objective function with the same weight. Similar to the

4. The W-MNS method 4.1. Preliminary results Let y "[y(1),2, y(¹)] be a sequence of ¹ zeroT mean random variables whose probability density function depends on a parameter vector h . Let 0 S( y ,h) be a parametric statistic vector of the obserT vation y defined as a linear function of the paraT meter vector h, and such that lim S( y ,h )"0 a.s. T 0 T?= We consider the problem of estimating (up to a scalar factor) the parameter vector h using the weighted least squares (WLS) criterion: hK "arg min S( y ,h)*WS( y ,h) W T T ,h,/1

(13)

(14) "arg min h*Q(T)h, W ,h,/1 where W is a positive weighting matrix, and Q(T) is W the quadratic form associated to the WLS criterion. The best choice of such a weighting matrix is given by the following theorem [16]. Theorem 1. Assume that the statistics sequence S( y ,h ) is asymptotically normal (i.e., [Re(S( y ,h ))T, T 0 T 0 Im(S( y ,h ))T]T is asymptotically normal), and let T 0 R be the asymptotic covariance matrix of S(y ,h ): s T 0 $%& lim ¹E[S( y ,h )S( y ,h )*]. R" s T 0 T 0 T?= ¹hen the parameter vector estimate corresponding to the weighted least squares solution Eq. (13) is consistent and asymptotically normal (i.e., [Re(hK )T, Im(hK )T]T is asymptotically normal) with W W

20

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15–27

mean h and asymptotic covariance matrix given by 0 $%& lim ¹E[(hK !h )(hK !h )*] R" h W 0 W 0 T?= "QjQ R Qj, W W sW W

(15)

where Q " lim Q(T), W W T?=

Q R " lim Q(T)R . W sW W sW T?=

(16)

Furthermore, the asymptotic covariance matrix R is h minimal if the weighting matrix is given by ¼ "R~1. s 015 In this case, the minimal parameter covariance matrix is given by R015"Qj . W015 h Remark. Theorem 1 still holds if we replace the optimal weighting matrix R~1 by an asympS totically consistent estimate [23]. Based on this result, we derive in the next section the optimal weighting matrix associated with the least-squares minimization criterion for the blind channel identification problem using the W-MNS method.

4.2. Statistical analysis and optimal weighting Before proceeding, we need the following additional assumptions: A1. The scalar polynomials h (z) satisfy: i deg(h (z))"M, ∀i and (h (z), h (z)) are coprime i i j for each pair iOj. In addition, it is assumed without loss of generality that EhE"1. A2. The window length parameter N is set to N"M#1. A3. The input source signal is an i.i.d. sequence of complex circular random variable with finite fourth order moment, and the additive noise is complex circular Gaussian distributed. In fact, the assumptions A1 and A3 can be relaxed at the price of more effort and notation. In the

sequel, the symbol ª is consistently used to denote estimated statistics, and for notational convenience, the dependence on the window length N is not indicated. Before going further, we need to define more precisely the solution hK that we consider. Note indeed that if hK is a solution of the least squares criterion, so is e*(hK , where / is an arbitrary real number; the objective function is unaffected by a multiplication of the arguments of f by a unitnorm complex scalar. In practice, a specific value of the phase should be imposed by some arbitrary constraint such as e*hK *0 where e is an arbitrary non-zero vector. The distribution of hK of course will depend upon the choice of vector e. However, any meaningful performance criterion will be insensitive to this choice, since this phase is unidentifiable from second order statistics. The most convenient (but admittedly artificial) choice for vector e turns out to be e"h. Implicitly, we define hK to be the determination of hK verifying hK *h*0. Now, let consider the least squares criterion Eq. (12). Using vector notation, Eq. (12) can be rewritten as hK "arg min Evec(T *( f )VK )E2, n ,f,/1

(17)

where ‘‘vec’’ is the matrix vectorization operator which transforms any m]n matrix into mn]1 vector obtained by concatenating the columns of the matrix. We propose here to generalize the least squares approach by introducing in the above criterion a positive definite weighting matrix W. The channel parameter vector is then estimated by minimizing the WLS criterion hK "arg min vec(T *( f )VK )*W vec(T *( f )VK ). W n n ,f,/1

(18)

Using Section 4.1 notations, the random variable sequence y , the parameter h , and the parametric T 0 statistics S( y ,h) correspond to T y "[ y(1),2, y(¹)], where y(t)"[y (t),2,y (t)]T: T 1 q The observations, h "h"[hT,2,hT]T: The channel parameter vecq 1 0 tor,

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15—27

C DC D C D

S( y , f )"vec(T *( f )VK ) T n

C

T *( f )€K F

"

T *( f )€K

1

DC

T *( f )€K 1 1 F

"

q~1

T *( f

)€K q~1 q~1

D

:

The statistics vector, $%& [T( f )T,T( f )T]T (k"(k ,k ) being where T ( f )" k2 1 2 k k1 the kth pair of channels), and ¹ is the sample size. The optimal weighting matrix is given by the inverse of the asymptotic covariance matrix R of s S( y ,h). This matrix R is given by Theorem 2. T s Theorem 2. ºnder the assumptions A1—A3, the asymptotic covariance matrix of the statistics S( y ,h) is given by T

C

R" s

R1,1

2

R1,q~1

F

2

F

Rq~1,1 2 Rq~1,q~1

D

,

(19)

where Rij are the (2M#1)](2M#1) asymptotic correlation matrices of the ith and jth subvectors $%& of S( y ,h) (i.e., T *(h )€ and T *(h )€K , where T *(h )" T i 6i j j j * * [T (h ),T (h )]): i1 i2 Rij" lim ¹E(T *(h )€K €K *T(h )) j i 6i6j T?= "p~4 + (€*Rji(!q)€ )Tj(h )Rij(q)Tj(h )* s Z 6j i j 6i q| 1 n * (€ (u)Sji(u)€ (u))T (u)Sij(u)T *(u) du, " j i j 2pp4 6i s ~n 6

P

(20) p2 being the input source power and s $%& E(Y (n#q)Y*(n)), Rij(q)" i 6j € (u)"[v (u),v (u)]T, i2 6 i1 6i M with v (u)" + € (l )e~+lu, 6 ik 6 ik l/0

21

h (u) h (u) * j1 Sij(u)"p2 i1 s h (u) h (u) i2 j2 di , j di , j 1 2 , #p2 1 1 di , j di , j 2 1 2 2 T (u)"[T (u),T (u)], i2 i i1 M T (u)" + T (l )e~+lu, ik ik l/0 with Tj(h )"[T (0),2, T (M),T (0),2,T (M)]. i i1 i1 i2 i2 Proof. See Appendix A. Discussion. To provide some insights into the actual performance of the method, we now give some comments on the above expression. f The weighting can be seen as an approach to compensate the effect of ill matrix conditioning due to close common zeros of the channel transfer functions h (z), i"1,2,q. By using the optii mal weighting matrix, the estimation procedure becomes quite insensitive to the ill conditioning problem (see Fig. 7 in Section 5). In fact, if a pair of channels have close common zeros, the corresponding Sylvester matrix T(h ) becomes nearly i singular and consequently Tj(h ) becomes ‘large’ i as well as Rij (see Eq. (20).) Therefore, the optimal weighting matrix given by the inverse of R will, qualitatively, provide ‘more weighting’ s (i.e., larger weighting coefficients) to the noise vectors associated with well conditioned Sylvester matrices T(h ) than to those corresponding to i ill conditioned Sylvester matrices. f The performance of the W-MNS approach does not depend on the higher order statistics of the source signals. In fact, this is a property [6] of all subspace based algorithms. f Under assumptions A1 and A2, the noise polynomial vector € (z) is given (up to a scalar factor) by 6i € (z)"[!h (z),h (z)]T. i2 i1 6i Therefore in Eq. (20), we have

C

di , j €*(z)Sji(z)€ (z)"p2€*(z) 1 1 6i 6j 6j di , j 2 1

D

di , j 1 2 € (z). di , j 6 i 2 2 (21)

22

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15–27

As we can see, the asymptotic covariance matrix is proportional to the noise power. This confirms the fact that we can perform the exact estimation with a finite sample size in the noiseless case. f In addition, we can see from Eq. (21) that two noise vectors associated with two pairs of channels not sharing any common channel are asymptotically decorrelated, i.e., Rij"0 if i"(i ,i ) and j"(j ,j ) do not share any com1 2 1 2 mon member. f Under assumption A2, the covariance matrices Rij(q) are given by

approximated by p2 R + + (€*Kji(!q)€ ) J(q). (24) ij p2 Z 6 j 6i s q| In this case, the two-step approach is not needed since the optimal weighting matrix does not depend explicitly on the unknown channel parameters.2 Also, since R given by Eq. (24) is ij Toeplitz, R~1 can be computed in O(q3M2) flops s using fast inversion algorithms, e.g. [7].

5. Simulation results $%& E(Y (n#q)Y *(n)) Rij(q)" i j "p2T(h ) J(q)T *(h )#p2Kij(q), s i j

(22)

$%& $%& E(s(n#q)s*(n)) and p2Kij" where p2J(q)" s $%& [N T(n),N T(n)]T. E([N (n#q)N *(n)]) with N (n)" i2 i i i i1 By replacing Eq. (22) in Eq. (20), we can rewrite the asymptotic covariance matrix R in the S time-domain as follows:

Validation of the asymptotic covariance expressions. We compare here the theoretical expressions with those obtained by Monte-Carlo simulations. We consider a three-variate (q"3) MA model, with coefficients given by (M"2) h (z)"(!0.0218#0.3867i)!(0.1546#0.2317i)z~1 1 #(0.2948!0.1372i)z~2,

Rij"p~4 + (€*Rji(!q)€ )Tj(h )Rij(q)Tj(h )* s Z 6j i j 6i q|

h (z)"!(0.3653#0.1286i)!(0.0048#0.2785i)z~1 2 #(0.2886#0.11497i)z~2,

p2 " + (€*Kji(!q)€ )[p2J(q) p4 Z 6 j 6i s s q| #p2T (h )jK ij(q)Tj(h )*]. i j

We present here some numerical simulations to first validate the expressions obtained for the asymptotic channel parameter covariance matrices and then to give some comparisons between the MNS and SS methods and the W-MNS method.

(23)

The computation of R~1 requires a total of s O(q3M3) flops. f The optimal covariance matrix depends on the unknown parameters we need to estimate. In practice, we must proceed in two steps: (i) estimate the unknown parameters using the original MNS algorithm; (ii) use the estimated parameters to calculate an asymptotically optimal weighting matrix and then perform the parameters estimation using the optimally weighted least squares criterion. However, it is interesting to notice that at high signal-to-noiseratio (SNR), the covariance matrix R can be s

h (z)"(0.2252!0.3151i)#(0.1465!0.2369i)z~1 3 #(0.0576#0.3200i)z~2. (25) The output observation noise is an i.i.d. sequence of zero-mean Gaussian variables; the width of the temporal window is N"3 and the input signal is an i.i.d. sequence of zero-mean, unit-variance QAM4 variables independent from the noise. For each experiment, N "100 independent r Monte-Carlo runs are performed, the number of

2 Note that under assumption A2, both J(q) and Kij(q) are known matrices (see [2]). Their expressions are omitted here for sake of simplicity. Moreover, the sum in Eq. (24) is finite due to the fact that Kij(q)"0 for DqD*M#1.

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15—27

samples ¹ is held constant (¹"1000), and the noise variance is varied between 0 dB and !30 dB. The simulation performance is measured by the mean-square error (MSE) defined as

S

¹ Nr + EhK !hE2. N r 1 The asymptotic value of MSE should be given by the square-root of the trace of the asymptotic covariance matrix, i.e., c(h;p2)"JTrace(Rh). Fig. 3 compares the simulation MSE (dashed line) and the asymptotical value c(h;p2) (solid line) for the MNS algorithm as a function of $%& p2EhE2/(qp2), where p2 is the input source SNR" s s power. Fig. 4 shows a similar comparison for the optimal W-MNS method. These figures demonstrate that there is a close agreement between the theoretical and experimental values. MSE"

Performance comparisons. We now present a comparison of the W-MNS method against the other methods. We consider a 4-variate (M"4) MA model. The first channel is given by the GSM test channel [17] with 6 delayed paths (M"5), the other three channels are generated by assuming a plane propagation model for each path with corresponding

23

electric angles uniformly distributed in [0, p/3] (i.e., h (z)"+M h (t)e+khtz~t with h 3[0, p/3]). t/0 1 k t The output observation noise is an i.i.d. sequence of zero-mean Gaussian variables and the input signal is an i.i.d. sequence of zero-mean, unit-variance QAM4 variables independent from the noise. Fig. 5 compares the performances of the MNS method (solid line) and the W-MNS method (dashed line). It shows the asymptotic MSE c(h; p2) as a function of the SNR in dB. This shows a high performance gain of the W-MNS method in comparison with the MNS method. Fig. 6 compares the performances of the WMNS method (solid line), the SS method (dashed line), and a weighted subspace (WSS) method [3] (dotted line). It should be noted that besides a performance gain over the MNS method, the W-MNS method is less accurate than the SS and WSS method (computationally the MNS method is the fastest). In the last experiment, the SNR is set to 10 dB and the channel coefficients are set to h (z)"(!0.0218#0.3867i)!(0.1546#0.2317i)z~1 1 #(0.2948!0.1372i)z~2, h (z)"!(0.3653#0.1286i)!(0.0048#0.2785i)z~1 2 #(0.2886#0.11497i)z~2,

Fig. 3. MSE of the MNS method versus SNR.

24

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15–27

Fig. 4. MSE of the optimal W-MNS method versus SNR.

Fig. 5. Comparison of the MNS (——) and the W-MNS (- - - -) methods.

h (z)"(!0.0218#0.3867i) 3 !(0.1546#0.2317i)jz~1 #(0.2948!0.1372i)z~2. The parameter j is used to control the distance between the roots of the scalar polynomials h (z) 1

and h (z). It is varied within the range [0,0.9]. Fig. 7 3 shows the asymptotic MSE as a function of the parameter j. It shows that generally the performance gain of the W-MNS method becomes more significant when the channels become more ill conditioned.

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15—27

Fig. 6. Comparison of W-MNS (——), SS (- - -) and WSS (.....) methods.

Fig. 7. Performance of the MNS (——) and W-MNS (- - - -) methods versus the parameter j.

25

26

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15–27

Finally, we note that although the performance gain is expected of the W-MNS method asymptotically, we have to be cautious when interpreting the results. In fact, when the number of samples is small, the MNS method can be more accurate than the W-MNS method (a similar situation happens to the SS method, see [3,13]).

6. Conclusion This paper presents a weighted version of the MNS estimation method for blind identification of multichannel FIR filters. It exploits some well known statistical results concerning the optimal weighting matrix of a least squares minimization problem. The explicit derivations of the optimal weighting matrix are presented. Also the asymptotic performance analysis of the proposed algorithm has been carried out, and the closed-form expression for the asymptotic covariance of the channel estimates has been obtained. Numerical simulations are provided which demonstrate the potential of the proposed method and support our theoretical findings. Finally, we remark that although the MNS method is a much (orders of magnitude) faster version of the SS method (see discussions following Eq. (12),) the W-MNS method is approximately as costly as the SS method. This paper has concluded that the additional computational cost associated with the WMNS method does not trade well for the accuracy improvement in comparison to the SS method.

Appendix A. Proof of Theorem 2

2N, are jointly asymptotically normal with zero mean and the covariance: $%& lim ¹EM(RK (i)!R(i))(RK (j)!R(j))N C (i, j) " cd ab cd ab abcd T?= "+ Rij (q)Rji (!q)#i (i, j), (26) cb abcd ad q|Z i (i, j) abcd $%& + cum(Y (n),Y * (n),Y (n#q), Y * (n#q)) " i,a i,b j,c j,d q|Z (27) "i(R(i)!p2d )(R(j)!p2d ), cd ab ab cd where cum(x ,2,x ) stands for the cumulant of the 1 k random variable (x ,2,x ), d is the Kronecker sym1 k $%& cum(s(n),s(n)*,s(n),s(n)*) is the kurtosis bol, and i" of the random variable s(n). h By definition, €K is the least eigenvector of RK . i i Using standard perturbation results ([19], pp. 293—295), we have €K "€ !C (RK (i)!R(i))€ #o (¹~1@2), P i 6i 6i 6i where C is the pseudo-inverse of the noise free i covariance matrix R(i), i.e., C "(R(i)!p2I)j"p~2Tj(h )*Tj(h ), s i i i and thus, by using the orthogonality between € and 6i range (T(h )), i.e., T(h )*€ "0, we have i i 6i T(h )*€ "!T(h )*C (RK (i)!R(i))€ #o (¹~1@2), i 6i i i P 6i "!Tj(h )(RK (i)!R(i))€ #o (¹~1@2). (28) i P 6i Using Eq. (28) and the definition of the asymptotic correlation matrices Rij we can write Rij" lim ¹E(T *(h )€K €K *T(h )) j i 6i6j T?=

Under the assumption A3, the empirical estimate of the covariance matrix RK given by Eq. (10) is i asymptotically complex Gaussian. We have the following Lemma (see [10], Theorem 14, p. 228).

" lim ¹E(Tj(h )(RK (i)!R(i))€ €* i 6i6j T?= ](RK (j)!R(j))Tj(h )*). (29) i Now, let us consider the (a,b)th entry of Rij. By using Lemma 3 result, we have

Lemma 3. ºnder the assumption A3, the random variables J¹(RK (i)!R(i)), 1)i)q!1, 1)a,b) ab ab

Rij " + lim ¹E(Tj(h ) (RK (i)!R(i)) kl ab i ak kl k,l,m,n T?= ](€ €*) (RK (j) !R(j) )Tj(h )* ) mn i nb 6 i 6 j lm mn

K. Abed-Meraim, Y. Hua / Signal Processing 69 (1998) 15—27

" + Tj(h ) (€ €*) Tj(h )* C i nb klmn i ak 6 i 6 j lm k,l,m,n "+ (€*Rji(!q)€ )[Tj(h )Rij(q)Tj(h )*] i j ab 6j 6i q|Z #i[Tj(h )(R(i)!p2I)€ €*(R(j)!p2I)Tj(h )*] . i i ab 6i6j (30) Using again the orthogonality relation T(h )*€ "0, i 6i we have (R(i)!p2I)€ "p2T(h )T(h )*€ "0. s i i 6i 6i From Eqs. (30) and (31), we finally obtain

(31)

(32) Rij"+ (€*Rji(!q)€ )Tj(h )Rij(q)Tj(h )*, i j 6i 6j Z q| which can be rewritten in the frequency domain as shown by equation Eq. (20). References [1] K. Abed-Meraim, Y. Hua, Blind identification of multiinput multi-output system using minimum noise subspace, IEEE Trans. Signal. Process. 45 (1) (January 1997) 254—258. [2] K. Abed-Meraim, P. Loubaton, E. Moulines, A subspace algorithm for certain blind identification problems, IEEE Trans. Inform. Theory 43 (March 1997) 499—511. [3] K. Abed-Meraim, J.F. Cardoso, A. Gorokhov, P. Loubaton, E. Moulines, On subspace methods for blind identification of single-input multiple-output FIR systems, IEEE Trans. Signal Process. 45 (1) (January 1997) 42—55. [4] K. Abed-Meraim, P. Duhamel, D. Gesbert, P. Loubaton, S. Mayrargue, E. Moulines, D. Slock, Prediction error methods for time-domain blind identification of multichannel FIR filters, in: Proc. ICASSP, Vol. 3, 1995, pp. 1968—1971. [5] R. Bitmead, S. Kung, B.D.O. Anderson, T. Kailath, Greatest common division via generalized Sylvester and Bezout matrices, IEEE Trans. on Automatic Control 23 (6) (1978) 1043—1047. [6] J. Cardoso, E. Moulines, A robustness property of DOA estimators based on covariance, IEEE Trans. Signal. Process. 42 (November 1994) 3285—3287.

27

[7] K.A. Gallivan, S. Thirumalai, V. Vermaut, P. Van Dooren, High performance algorithms for Toeplitz and block Toeplitz matrices, J. of Linear Algebra and Its Applications 241—243 (1996) pp. 343—388. [8] W. Gardner, A new method of channel identification, IEEE Trans. Commun. 39 (August 1991) 813—817. [9] G. Giannakis, J. Mendel, Identification of non-minimum phase systems using higher-order statistics, IEEE Trans. Acoust. Speech Signal Process. 37 (1989) 360—377. [10] E. Hannan, Multiple Time Series, Wiley, New York, 1970. [11] Y. Hua, Fast maximum likelihood for blind identification of multiple FIR channels, in: Proc. of the 28th Asilomar Conference, CA, 1994, pp. 415—419; IEEE Trans. Signal Process. 44 (March 1996) 661—672. [12] Y. Hua, K. Abed-Meraim, M. Wax, Blind system identification using minimum noise subspace, IEEE Trans. Signal Processing (and also in 8th IEEE workshop on SSAP, Corfu, Greece, June, 1996) 45 (3) (March 1997) 770—773. [13] M. Kristensson, B. Ottersten, Statistical analysis of a subspace method for blind channel identification, in: Proc. ICASSP, Atlanta, Vol. 5, May 1996, pp. 2435—2438. [14] H. Liu, G. Xu, L. Tong, A deterministic approach to blind identification of multichannel FIR systems, in: Proc. ICASSP, Vol. 4, 1994, pp. 581—584. [15] E. Moulines, P. Duhamel, J. Cardoso, S. Mayrargue, Subspace methods for the blind identification of multichannel FIR filters, IEEE Trans. Signal Process. 43 (February 1995) 516—525. [16] B. Porat, Digital Processing of Random Signals, PrenticeHall, Englewood Cliffs, NJ, 1994. [17] J. Proakis, Digital Communications, McGraw-Hill, New York, 1989. [18] O. Shalvi, E. Weinstein, New criteria for blind deconvolution of nonminimum phase systems (channels), IEEE Trans. Inform. Theory 36 (2), 1990. [19] G. Stewart, Introduction to Matrix Computations, Academic Press, New York, 1973. [20] L. Tong, G. Xu, T. Kailath, A new approach to blind identification and equalization of multipath channels, in: Proc. of the 25th Asilomar Conference, CA, 1991, pp. 856—860. [21] A.J Van der Veen, S. Talwar, A. Paulraj, Blind identification of FIR channels carrying multiple finite alphabet signals, in: Proc. ICASSP, 1995, pp. 1213—1216. [22] A. Venetsanopoulos, S. Alshebeili, E. Cetin, Cumulant based identification approaches for non-minimum phase FIR systems, IEEE Trans. Signal Process. 41 (April 1993) 1576—1588. [23] M. Viberg, B. Ottersten, T. Kailath, Detection and estimation in sensor arrays using weighted subspace fitting, IEEE Trans. Signal Process. 39 (November 1991) 2436—2449.