Blind equalisation of multilevel PAM data for nonminimum phase channels via second- and fourth-order cumulants

Signal Processing 31 (1993) 313-327 Elsevier

313

Blind equalisation of multilevel PAM data for nonminimum phase channels via second- and fourth-order cumulants Fu-Chun Zheng, Stephen McLaughlin and Bernard Mulgrew Department of Electrical Engineering, University of Edinburgh, UK Received 12 March 1992 Revised 22 June 1992

Abstract. This paper introduces a new blind equalisation algorithm for the pulse amplitude modulation (PAM) data transmitted through nonminimum phase (NMP) channels. The algorithm itself is based on a noncausal AR model of communication channels and the second- and fourth-order cumulants of the received data series, where only the diagonal slices of cumulants are used. The AR parameters are adjusted at each sample by using a successive over-relaxation (SOR) scheme, a variety of the ordinary LMS scheme, but with a faster convergence rate and a greater robustness to the selection of the 'step-size' in iterations. Computer simulations are implemented for both linear time-invariant (LTI) and linear time-variant (LTV) NMP channels, and the results show that the algorithm proposed in this paper has a fast convergence rate and a potential capability to track the LTV NMP channels. Zusammenfassung. Vorgestellt wird ein neuer voraussetzungsfreier ('blinder') Equalizer-Algorithmus fiir die Dateniibertragung mit Hilfe der Puls-Amplitude-Modulation fiber nicht-minimalphasige Kan/ile. Der Algorithmus basiert auf einem nichtkausalen autoregressiven Modell (AR) des Kommunikationskanals und der Kumulierenden 2. und 4. Grades der empfangenen Datenkette, wobei nut die Diagonalen der Kumulierenden herangezogen werden. Die AR-Parameter werden fiir jeden Abtastwert adaptiert; hierbei wird ein sukzessives Uberrelaxationsschema (SOR) eingesetzt, also eine Variante des gew6hnlichen KleinsteQuadrate-Algorithmus, die jedoch schneUer konvergiert und in bezug auf die Auswahl der 'Schrittweite' der Iterationsschritte robuster ist. Das Verfahren wurde sowohl fiir lineare zeitinvariante als auch fiir lineare zeitvariante nicht-minimalphasige Ubertragungskanfile simuliert. Wie die Ergebnisse zeigen, weist der hier vorgestellte Algorithmus giinstige Konvergenzeigenschaften auf und kann grunds/itzlich aucb im Fall der zeitvarianten Ubertragungskan~ile eingesetzt werden. R6sum6. Nous introduisons dans cet article un algorithme nouveau d'6galisation aveugle pour des donnees ~ modulation d'amplitude puls6e (PAM) transmises par des canaux d phase non minimum (NMP). L'algorithme lui-m~me est bas6 sur un module AR non causal des canaux de communication et les cumulants de deuxi+me et quatri+me ordre des s6ries de donn6es regues, les tranches diagonales des cumulants 6tant seules utilis6es. Les param6tres AR sont ajust6s d chaque 6chantillon ~t I'aide d'une technique de sur-relaxation successive (SOR), une variante de la technique LMS ordinaire mais poss6dant un taux de convergence accru et une plus grand robustesse vis-a-vis de la s61ection du facteur de mise d jour dans les iterations. Des simulations sur ordinateur sont implant6es d la fois pour des canaux NMP lin6aires invariant dans le temps (LTI) et des canaux NMP lin6aires variant dans le temps (LTV), et les r6sultats montrent que l'algorithme propos6 dans cet article poss6de un taux de convergence 61ev6 et la capacit+ potentielle de suivre des canaux NMP LTV.

Keywords. Blind equalisation; non-minimum phase; higher order spectra; higher order statistics.

Correspondence to." Dr. Stephen McLaughlin, Department of Electrical Engineering, The University of Edinburgh, The King's Building, Mayfield Road, Edinburgh EH9 3JL, Scotland, UK.

0165-1684/93/$06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved

314

F.-C. Zheng et al. / Equalisation of multilevel PAM data

I. Introduction

Developments in modern communications, e.g., in multiplexing and multipoint network communications [12, 22], have necessitated the need for blind equalisation techniques. The main advantage of this technology is that the training stage, which is necessary in conventional equalisation schemes but not always possible, is removed. In fact, over the past two decades there has been a growing interest in blind equalisation both in theory and in practical applications (see [1, 2, 12-14, 17, 22] and references therein). The problem of blind equalisation for minimum phase (MP) channels has been well understood. In that case it reduces to a classical innovation restoration problem [2, 16] and the well-known prediction scheme is successful [17]. In practice, however, not all channels are MP. As a result, the equalisation of nonminimum phase (NMP) channels has been investigated, and a range of different techniques have been proposed. The most promising of these are based on higher-order cumulant (HOC) analysis (alternatively termed higher-order spectral (HOS) analysis, HOS being the Fourier transformation of HOC) [7, 13-15, 21, 27]. The vital reason why HOC can be used for this purpose is that the phase property of signals is preserved in HOC up to a linear phase shift (see [5, 13-15, 18, 20, 21, 27, 30] and references therein), whereas in autocorrelation functions this is not the case. Many HOC-based algorithms for blind equalisation (deconvolution) and identification of NMP channels have appeared over the past decade. As in the autocorrelation function domain, these algorithms can be divided into two categories: MA-based and AR-based. Considering the linearity of the algorithm and the convenience of equalisation implementation, AR model-based schemes [5, 15, 18] are preferable. Since the considered channels are NMP, the AR model here is noncausal. However, there are several drawbacks for the existing noncausal AR based method [18]: (i) both diagonal and nondiagonal cumulant slices are used. Since the estimation of nondiagonal Signal Processing

cumulant slices is generally noisier than that of diagonal slices [3], using nondiagonal slices may result in degradation of algorithm performance; (ii) the form of algorithms are too complicated; and (iii) the convergence rate can be too slow. Based on the motivation of overcoming the above disadvantages, we introduce in this paper a new blind equalisation algorithm in which only autocorrelation functions (second-order cumulants) and the diagonal slice of fourth-order cumulants are involved. As a result, the form of the concerned matrices become very regular and simple, and the computation requirement also decreases. In the situation studied in this paper, the transmitted data are assumed to be of multilevel pulse amplitude modulation (PAM) type: equally distributed over some discrete set {+1, 4-3, +5 ....

}.

In adapting (or updating) the noncausal AR parameters, we adopt a successive overrelaxation (SOR) scheme, which is a special variety of the conventional least-mean-square (LMS) adaptive scheme, but which possesses a faster convergence rate and a greater insensitivity to different values of 'step-size', termed 'relaxation factor' here. The general structure of this paper is as follows. In Section 2, the blind equalisation problem is formulated. Section 3 derives the new algorithm. The comparison of SOR with the LMS scheme is made in Section 4. Section 5 analyses the convergence of our algorithm. Several computer simulations are shown in Section 6. Finally, some conclusions are drawn in Section 7. A concise report on some of the results in this paper can be found in [28]. 2. Problem formulation The channel studied in this paper is as shown in Fig. 1, and it is assumed that the channel is stable and linear. Then, it may be expressed as the following MA model with finite order q: q

y ( k ) =ynf(k) + n ( k ) = ~ b i x ( k - i) + n ( k ) , i=0

(2.1)

315

F.-C Zheng et aL / Equalisation of multilevel PAM data

In the situation considered here, as a result of symmetric distribution (or unskewedness) of the above PAM data, the third-order cumulants of x(k), hence that of y(k), theoretically vanish (the definition of HOC and its properties can be found in [4, 20]). Thus, it is necessary to use the fourthorder cumulants of y(k), which, mathematically, may be expressed as

............................................... x(k)

I

Channel

-)

ib, I

y(k) D

I

Fig. 1. The channel model.

where n(k) is white Gaussian noise, and x(k) is the transmitted sequence in which the random data are independently and uniformly distributed over some discrete finite set: {+1, +3 . . . . . + ( 2 M - 1 ) } , which corresponds to a 2M-level PAM scheme. Also, we assume that the transfer function of the above channel has no zeros on the unit circle. The object of blind equalisation is to restore x(k), given only y(k), by determining an inverse filter {0t} which lets (2.2) hold: r2

2(k)= ~ Oty(k-l)~x(k). /

(2.2)

rl

Here, rl is the order of the noncausal part of the inverse filter, and re that of the causal part. In fact, when rl = 0, (2.2) represents the case of MP channels, where the well known prediction technique can be used to realise the blind equalisation. On the other hand, when r2 = - - l , (2.2) models the maximum phase channels, for which the prediction method can also be employed [2]. However, for mixed phase channels (i.e., rl <0 and r 2 / > 0 ) , n o correct blind equalisation can be achieved by using the prediction scheme. As a result, several other schemes have been proposed to solve the problem of blind equalisation of NMP channels, e.g., see [ 1, 2, 5, 18] for details. It should be noticed that, if and only if M is sufficiently large, the sequence x(k) is close to being sub-Gaussian [1, 2]. In practice, however, M is usually not taken to be too large, and this is a limitation of the algorithms reported in [1, 2]. Fortunately, HOC-based methods, such as the one presented in this paper, do not require any assumption about distribution other than that it be non-Gaussian. This provides a further motivation in addition to those mentioned previously for employing HOC.

C4(m,, me, m3) =ElY(k) i~_, y(k+mi)] - G(m,)C2(m3-m2)

- C2(m2)C2(m3 - m,) - C2(m3)C2(m2- ml),

(2.3)

where, C2(') denotes the autocorrelation function of y(k), viz., C2(m) = E[y(k)y(k + m)].

(2.4)

In this paper, we only use the diagonal slice of C4(ml,m2,m3), where ml=m2=m3, and we denote it as c 4 ( m ) = C 4 ( m l , m 2 , m3)lm,

m2 m3

m'

(2.5)

For the autocorrelation function, the diagonal slice c2(" ) of (72(') reduces to (72(") itselt, i.e.,

c2(m) = C2(m).

(2.6)

Thus, we have c4(m) = E[y(k)y3(k + m)]

- 3E[y2(k)]E[y(k)y(k + m)].

(2.7)

On the other hand, it should be noted that determination of the inverse filter { 0t} is equivalent to modelling of channels by a noncausal AR structure, viz., q 1 H(z)= 5~ biz-' ~ r2 ,=0 ~ O~z ~ •

(2.8)

l=rl

Consequently, the determination of the inverse filter {0/} can also be considered as a process of Vol. 31, No. 3, April 1993

F.-C. Zheng et al. / Equalisation o f multilevel P A M data

316

channel identification (AR model parameter identification). In this paper, a new second- and fourth-order cumulant-based SOR algorithm is used to adapt 0;s and realise the corresponding blind equalisation procedure.

and V = (Oi)(r2--rl + 1) = [--c4(--r2) .....

--c4(--1),

gq, --c4(q + 1 ) .....

3. Algorithm derivation

In this section, the block form of our algorithm is described first, and then the adaptive version of the algorithm. The SOR scheme is derived and discussed. To determine the inverse filter parameters 0;s in (2.2) or (2.8), the following set of second- and fourth-order cumulant-based equations was obtained in [26] and [29] under the conditions defined above:

(3.4b)

U = ( UO")(r2-- rI + 1) x (r2-- rl)

go,...,

--C4(--rl)]T.

(3.4c)

In practice, c2(" ) and c4(' ) will have to be replaced with their estimates ~2(') and G(" ), respectively. But we will not explicitly indicate the 'hats' in the following parts of this paper to aid the conciseness of description because it is easy to distinguish. From (3.3), the following normal equation can be obtained: uTuT`

(3.5)

= uTv.

Since UTU is positive definite, we must have det(UxU) ~0.

r2

(3.6)

~, Ojc4(n+j)=-e4(n) j--rl j~O

for n~[-r2, -1] w [q+ 1, - r d ,

(3.1)

r2

Z Ojf,j=g,

From (3.6) and the concerned linear equation theory, the solution of (3.5), then the least squares solution of (3.3) exists uniquely, viz.,

for nE[0, q],

i//= ( u T u ) - I (U T V).

j= rl j#O

where q

fnj = • [c2(n - i)c4(iWj) - c4(n - i)c2(i+j)] i=0

(3.7)

Equation (3.7) is the block form of our algorithm. The elements of U and V can be recursively estimated at each sample (or time step) by using the following expressions:

(3.2a)

c(2k+')(m)

and q

g,= ~ [c4(n-i)c2(i)-c2(n-i)c4(i)].

1

(3.2b)

i=0

The difference between 0/s in (3.1) and Ot's in (2.2) is that normalisation has been done in (3.1) so that 0jb=0 = I. Equation (3.1) can be rewritten in the following more compact form: U7` = V,

(3.3)

-

k-m+l

k+2

E

x(i)x(i+m)

~=0

1 ---[(k+ k+2

1)c~k)(m)+x(k-m+ 1)x(k + 1)]

1 )+ x ( k - m + 1)x(k+ 1) = c~k)(m) 1 - ~

where

for O<~m<~q;

7' = (~)~_,,~ = [ O r I , Or,+l . . . . .

c~2k+ O(m) = e~k+ l)(-m) 0--1, O1 . . . . .

for -q~
(3.8a) (3.8b)

Or2--1, 0,2] T,

(3.4a) Signal Processing

k+2

c~k+l)(rn) = 0,

other,

(3.8c)

317

F.-C. Zheng et al. / Equalisation of multilevel PAM data

and c~4*+ I)(m) = M~4k+ ')(m) - 3c~k+ 1)(0)c~k+ l)(m), (3.9)

change with the increase of time k. On the other hand, the ordinary LMS type of scheme can be borrowed and used in our problem. From (3.3), we can define

where M4~k+l~(m) is the diagonal slice of the fourthorder moments at the ( k + l)th time step, and it can be estimated as follows:

-

1

k

k+2 1

k+2

m + 1

Z

,-0

× [U(k) ~'(k) - V(k)]

x(i)x3(i+m)

[(k+ 1)M4~k)(m)+x(k-m+ 1)x3(k+ 1)]

(

k+2 for0<<,m<~q;

T ( k + 1)= tP(k)+2r(k)V(k)

(3.12a)

V(k) = UT(k)[ V(k) - U(k) ~(k)],

(3.12b)

and

1 )qx(k-m+l)x3(k+l)

= M4~(m) 1 - ~

(3 1l)

as the cost function at the kth iteration, where U(k), V(k) and tP(k) respectively denote the estimates of U, V and q' at the kth iteration. Consequently, we can obtain the following update to ~(k):

M~4k+l)(m) --

e2(k) = [U(k) q'(k) - V(k)] T

(3.10a)

where r(k) is generally taken from the range 1

M4~k+ I)(m ) 1

-k+2

0 < r(k) < tr[UV(k) U(k) ] . k+l ~ x(i)x3(i+m)

i~ -m

= M4~k)(m) 1 - ~

1 )q x3(k+m+ l ) x ( k + 1) k+2 for-q~
M4~k+~)(m) = 0,

other.

(3.10b) (3.10c)

Thus, (3.7) also provides, in principle, a solution to the adaptive inverse filter parameters with the above recursive estimation of cumulants The direct calculation of the right-hand side of (3.7) at each sample, however, as in autocorrelation function based adaptive filtering methods, involves a considerable amount of computation In order to make full use of the information contained in the current set of parameters, the two most widely used schemes which have been proposed to solve such kind of problems in autocorrelation function domain are the recursive least-squares (RLS) and least-mean-squares (LMS) algorithms [6]. In our problem, however, the RLS type of scheme is not applicable. This is because the size and structure of both U and V are 'time-invariant', i e , they do not

(3.12c)

However, the convergence rate of the above scheme is generally slow, as in the autocorrelation function domain [5, 6]. In addition, the convergence rate is greatly dependent on the choice of r(k). In order to make best use of the information contained in the current set of parameters and obtain a faster convergence rate (without increasing the computation penalty too much), the successive overrelaxation (SOR) iterative scheme was adopted. The block diagram of our equalisation procedure is illustrated in Fig 2. The SOR scheme was introduced simultaneously by Young and Frankel in 1950 [23, 25]. In fact, it is an accelerated scheme of the Gauss-Seidel iterative s..............................................................................................

:

y(k)i

C2(' ) and ~ c4(') updating

SOR (LMS)

Fig. 2. The adaptive blind equalisation scheme. Vol. ~1, No 3, April 1993

318

F.-C. Zheng et al. / Equalisation o f multilevel P A M data

method for solving the large scale sets of linear equations. Like the Gauss-Seidel method, determination of the inverse of matrix is not required in the SOR method. But the SOR method can achieve a faster convergence rate. Now, we apply SOR scheme to (3.5), and the following iterative algorithm can be obtained. Let UT(k) U(k) = A (k) = e~(k) ~ i j )(r2--rl)×(r2--rl)

For the (k + 1)th iteration (or sample): (1) Filling U(k) and V(k) by using (3.8a-c) and (3.9 10c) ; (2) Calculate r2--rl+l

aii=

~

(u}/~))2, i = 1 , 2 . . . . .

r2-rl;

1--1

(3.15a) (3) For i= 1, 2 . . . . .

r2--rl,

(3.13a) if i = 1, then calculate

and

r 2 -- r I

(3.13b)

UT(k) F(k) = D ( k ) = (dff >)(,~_,o.

2

j=l Then, ~(k) can be updated by Iff~k + l ) :

(3.14a)

Iffl k) "~- 8 ~t~ k),

l= 1, 2 . . . . .

r2--rl + 1,

81#r~k)__ /U r2--~. +1 . (k)A(k) ~(k) Ull z~ I1 , Ull /=1

where -- ].1 ( d ( k )

8gtlk)---Z/5 ,,i -

ao

j=l

aii \

i=l, 2,...

(k) l~c(k+l)__,2-, , ~(k), (k)~ J E uij ~ j ], j=i

if

(3.14b)

, r2-rl,

and p is termed a relaxation factor, which, to guarantee the convergence of the above iterative process, must be chosen from the following range [23, 25]: 0 < p <2.

(3.14c)

When p = 1, SOR reduces to the Gauss-Seidel method. Furthermore, (3.14b) can be rewritten as r2 --r I + 1

P

E

(u P)2

,=1

/=1 i-- 1

x v~k ) - ~

j=l

ubk) ,

r2--rl • (k) ,,.(k)~

z_,

j=i

ulj

W j J,

(Yl4d) where u~(k) and "vi(k) denote the values of u~ and vi after the kth iteration, respectively. In practical programming, further arrangements can be made to reduce the computational load. For example, we can adopt the following algorithmic form: Signal Processing

(3.15b)

(3.15c) (3.15d)

otherwise, calculate A ( k ) _ .4(k) _ R , ( k ) , ( k ) li --'<--ll,i--I °tYi--lUl, i - - l ,

l=1,2 .....

r2-rl+l,

-- [-/ r2--~ + l • (k)A(k)

8gt}k)

a~ik)

,=l

~,}k + 1> = ~,}k) + 8

t~liZ'lli'

~ff).

(3.15e) (3.15f) (3.15g)

Notice that, in step (3) above, the information of Zal.iIA(k) is fully utilised to calculate A~ ). As a result, the computation amount is greatly reduced. Several discussions for the above algorithm are in order: (i) Initial estimate of T(k): In principle, T(0) can be chosen arbitrarily. But habitually, it is set to be 0. (ii) Since the obtained inverse filter is noncausal, the restored data series has a delay o f - r l samples (or time steps), as shown in Fig. 2, viz., r2

2'(k+rl)= ~ Ojy(k+rl-j). j=rl

(3.16)

319

F.-C Zheng et al. / Equalisation of multilevel PAM data

(iii) As a result ofnormalisation of 0/s (or Vi's), the output data of the inverse filter {0j} has an amplitude gain in comparison with the true transmitted data at each time step; thus, an automatic gain control (AGC) mechanism [13, 14, 17] is needed, as in Fig. 2. From (3.16), the power of 2'(k) at the kth time step can be expressed as

(7~,(k) = E[(2'(k)) 21 r2

r2

= Z Z OiO~E[y(k-i)y(k-j)] i-- rl j

rl

r2

r2

= Z

Z

OiOjc2(i--j).

(3.17)

i=rl j--rl

Notice that 00 = l in (3.16-17). Let the power of the transmitted data be o'x, 2 then the AGC factor can be chosen as dAGC= ±1o'2,/o'~,112.

(3.18)

Since the earlier estimates of c2( ) in adaptive procedure may cause o-2, < 0, the operation of taking the absolute values is added in (3.18) to force the blind equalisation to continue. The sign (+ or - ) here is not identifiable as in many existing methods [13, 14]. (iv) Choosing the relaxation factor it: Young obtained an expression for the optimum relaxation factor [23, 25]: itopt-

j l +./1

2 -

p2(y) ,

r, and this criterion does not always work well. In fact, other criteria have been suggested recently [8, 9]. However, this is not such a major problem for the above SOR algorithm since the selection interval for relaxation factor/1 is clear: 0
(3.19)

where p(Y) is the spectral radius of the Jacobi iterative matrix Y associated with (3.5). Since determination of p(Y) is generally very complicated, (3.23) cannot be used in practical environments. Empirical values can usually be found by trial and error. Fortunately, the disparity in the convergence rates caused by different relaxation factor values in our SOR algorithm is considerably less than resulted from different step-size values in the LMS scheme. In addition, it is well known that (3.12c) only represents an ad hoc basis for selection of step-size

4. Comparison of SOR with LMS In the foregoing sections, the adaptive SOR algorithm has been derived. In this section, the relation between SOR and LMS is illustrated, and the comparison between them is also made in terms of computational complexity. Now, we can rewrite (3.12a c) in their component forms: ~,~k+l) = qtlk)+ 2r(k)Vlk) r~"

v(k) -,

=

r I + 1 • (k)(

Z

/= i

-,, k

(4. la) \

r2-rl

,(k)_ x" j=l

.12'¢7'

)

(4.1b)

Vol. 31, No. 3, April 1993

320

F - C . Z h e n g et aL / Equalisation o f m u l t i l e v e l P A M

and

1 0 < r(k) < ~_~, ~_~, +

E

(4•1C)

E

i =1

1=1

Notice that, when we calculate ~/( k + l ) , the values for ~ + ~ ) . . . . . v,.,~k+ ,,,~,+~1) are already i - 2 ~) and v~available, and the latter are the most recently obtained values and generally better in accuracy than the values for ~0*), ~ ) , . . . . ~i-z (k) and ~k_)l. Thus, from the viewpoint of taking advantage of the most current information, we can change the ordinary form of LMS in (4. la, b) into the following accelerated form (alternatively termed the Gauss-Seidel f o r m ) ' - ~(k)

iltlk + 1)= ilt}k) + 2 r i ( k ) v i

(4.2a)

,

~'~ /=1

Uli

SOR equations (3.15a-g)

LMS equations (4.1a c)

-5 5 439 340 -10 10 1679 1280 -15 10 2599 1975 -15 15 3719 2820 -20 20 6559 4960 (AS number is approximately the same)

Increase 99 399 624 899 1599

and divisions ( M D ) for one iteration is approximately 3(r2 - rj) (r2 - rl + I) + (r2 - rl) ; In (3.15a-g) (programming implementation of (3.14a, b)), the corresponding number of M D becomes 2(r2 - rl + 1) (r2 - rl) + 2(r2 - rl + 1)(r2 - rt - 1)

= 4 ( r 2 - r l + 1)(r2- r l ) - - l.

i- 1

r2--r 1

-0

')- Y

J

j=i

j=l

.o

)• (4•2b)

Comparison of (4•2a, b) with (3•14a) and (3•14d) immediately leads to the conclusion that they are equivalent! Thus, the SOR scheme is actually an accelerated form of the ordinary LMS scheme• Considering the range (0, 2) of p, the above comparison also gives the new range for r(k) in (4.2a) :

O< ri(k)
1

The numbers of additions and subtractions (AS) are very similar to MD. Hence, the SOR scheme approximately requires ( r 2 - rl + 1 ) ( r 2 - r l - 1) = ( r z - r t ) 2 - 1 more M D and the similar number more of AS. Table 1 shows this point further by several concrete examples. Apparently, the increase of computation load is normally not very large. In addition, this can be considered as the penalty to be paid for the faster convergence rate.

5. Consistency and convergence •

(4.2c)

E

/=1

Notice that the range of r~(k)'s in the new form of LMS scheme is greatly augmented• On the other hand, it can be seen from the following discussion that the amount of computation for the SOR becomes slightly larger than that for LMS. Since the computational amounts involved in cumulant estimation are the same for both the SOR and LMS, we only need to compare the computational complexities of (3.15a-g) and (4. la-c). In (4. la-c), the number of multiplications Signal Processing

r2

+ (r=- r, + 1) + ( r 2 - r,)

r2 r I + 1

V} k ) =

Table 1 Computational complexitycomparison between SOR and LMS (MD number) rt

2

data

In this section, the consistency and convergence behaviour of our algorithm are analysed• It is shown that the algorithm is normally consistently convergent• First, let us notice that the assumptions in Section 2 imply that: (1) For the transmitted data series x(k), E[x(k)] = 0, and the diagonal slices of the ith-order cumulants of x ( k ) are limitedly large, where i = 2, 3, 4, 5, 6; and (2) F o r the channels: They are absolutely summable, viz., ~iqo [bil < ~ • Then, according to

F.-C. Zheng et al. / Equalisation of multilevel PAM data

the relevant theorem in [10], we can conclude that the estimators both in (3.8) and in (3.9) are unbiased and consistent. Also, it is assumed that the affect of additive noise is removed. Considering the nonsingularity of UT(k)U(k) and the results in [5, 24], we have lim tP(k)= lim ( U T ( k ) U ( k ) ) - ' ( U T ( k ) V ( k ) ) k~yJ

k~oc,

=q~.

o

asymptotically unbiased, for any finite value of k the estimate may be significantly biased [5]. Now, the remaining problem is: does the SOR scheme converge? and if yes, under what conditions? To answer these two questions, we only need to cite the following theorem. T H E O R E M . For a set of linear equations

(5.1)

Thus, the estimator in (3.7), then that in (3.5), is asymptotically unbiased and consistent. It should be noted, however, that although the estimator is I

I

I

321

(5.2)

Aq'=D,

if A is symmetric and positive definite and the relaxation factor is in the range 0 < p < 2, then the S O R scheme converges.

I

0

I

I

A. ~ ; "

=

1

A

-10-20-

I

.-

-~o-

;

L; ...... LMsL

-20-

~--~"- 3 0 -

~' -30 -

.~-40-

~ -40-_50_60__~

-50-

"%""'~'~ ~

""--'"~"-""

- ...........~'-"-~'"

................................

-60-

I3

-70 0.10 0

2.10

I3

I

4.10

6.10 3

I

I

I

I

2.10 ~

4.10 3

6.10 3

-70

8.10 3

1.10 ~

0.10 c

/c

I

I

1.0

I

I

I

I

I

B. ........ LMS I

.... 1 L---~

0,8"

~0.6-

r~ 0 . 4 -

~0.6-

~

~....,

0.4-

----

SOR

..... .....

0.2__r--- L 0,0 0.10 0

1.10"

Fig. 4(a). The M S E curves of the SOR and L M S algorithms for Example 2.

B. 0.8 ¸

8.10 3

k

Fig. 3(a). T h e M S E curves of the SOR and L M S algorithms for Example 1.

1.0

--I

.... i I

2.10

3 r--

I

4.10 3

6-10

3

I

8.10

0.0

3

1.10 ~

0.10 0

1 3

1 3

2.10

4.10

i

i

6.10 3

8.10 3

1.10~

/c

Fig. 3(b). T h e error rate curves of the S O R and L M S algorithms for Example 1.

Fig. 4(h).

The error

rate curves of the SOR rithms

and LMS

algo-

for Example 2. Vol 31, No. 3, April 1993

F.-C. Zheng et al. / Equalisation of multilevel P A M data

322

The rigorous proof for the above theorem can be found in [25]. Clearly, our algorithm completely satisfies the conditions of the above theorem; consequently, the convergence of our algorithm is guaranteed.

publically available report about its fourth-order cumulant version, the CN method bears at least the three drawbacks mentioned in Section 1 of this paper, deducing from its third-order cumulant version which has been well reported [5, 18]. As to the non-HOC-based blind equalisation techniques, their performance is in general inferior to the HOC-based approaches, although the former is computationally much simpler. This has been well documented (see [13 15, 30] and the references therein). As a consequence, we did not carry out any comparison of our algorithm with other blind equalisation techniques in the following simulations.

6. Simulation examples In this section, several Monte Carlo simulation examples are presented to demonstrate the performance of the algorithm proposed above. In the description below, the signal-to-noise ratio (SNR) is defined as SNR = 10

1

2 2 Ogl0(O'ynJO'n).

(6.1) 6.1. The L T I channels

Here, a2,r denotes the variance of the noise free output series Ynf, and 0-2 that of the additive noise, as shown in Fig. 1. In addition, it should be pointed out that, among the existing HOC based blind equalisation techniques, the only one which is comparable to our algorithm is the Chiang-Nikias (CN) method [5, 18]. Although there has been no detailed and

E X A M P L E 1. The mixed phase channel transfer function is assumed to be H(z) = 1 -2.3z -1 +0.6z -2.

(6.2)

The transmitted data is a four-level PAM series: equally distributed over {+ 1, +3 }, and the SNR is

A ° 1

1.5

1

1

1

1

:

:

:

:

-

1.5

i

l

l

:

:

:

1.5

:

:

:

-

:

:

0

. . . . . . . . .

1.0

. . . . . . . . .

0.5

0.5

. . . . . . . . .

0.5

. . . . . . . . .

0.0-

:

1.5

0,0

7.5

15.0

.

.

.

.

.

v

0.0

:

:

:

:

:

~

:

~

:

0.5,

..

.. . .

.

15.0

:

:

:

1.5

:

1.0

.

:

.

:

.

.

~

.

~

.

.

:

:

0

;

:

7.5

•

15.0

k/103

~

.

~

0.0

i

~

:

:

1.0~

.

. . . . . . . . .

0.5

0.0

:

:

;

:

0.0

.

:

:

7.5

:

;

15.0

0.0 . . . . . . 0.0 7.5

convergence

traces

.

.

.

.

15.0

kilo ~

k/103 The

..

7.5 k/103

1.0,

5(a).

l

1.0+

k/lO 3

Fig.

:

1.0

0,0. 0.0

signal Processing

:

of the inverse

filter coefficients

:

0.0

:

:

;

i

;

n

n

7.5

15.0

k/103 in Example

3 for the SOR

algorithm

with

p = 1.2.

323

F.-C. Zheng et al. ,/Equalisation o f multilevel P A M data I

the approximately fastest converging speed. Clearly, SOR-based algorithm is o f a far faster convergence rate than the LMS-based one. But the MSE curve of the SOR algorithm is not as smooth as that of the LMS algorithm.

B.

-10 -20 ~" -30 -

~

EXAMPLE 2. Let the channel transfer function

-40-

be -50-

H(z) =

1 + 4 . 5 z -1 + 2 z -2.

(6.3)

-60- 70 0.0.100

I

3

I

5.0.10

z.

1.0.10

1.5.10 ~"

Here, we set/1 = 1.5 and r(k) = 0.9/tr[UT(k)U(k)]. The other assumptions are the same as Example 1.

k

Fig. 5(b). The MSE curve of the SOR algorithm for Example 3 with # = 1.2.

I

0 - 1 0 - ,,

A.

I

0

} "-'""" 1~=1"0

-20C.

-10-

-

-

'"

I

I---- ~,=1.5 I

'',.

/~=0,8 /.L=I,8

L::-- ,,_-,.9 I

~" -30 -

-20-

.~-40-

~., - 3 0 -

~

I

-50-

'.,,, t-,. "'-';,.

-40-60-50,

,

-70

j

I 2.10 3

0.100

I 4.10 B

-60-

-70

0.0-100

I 6.10 3

I 8.10 3

I*I0 L

k

I

5.0.10

I

3

1.0.I0

t.

1.5.I0 z"

Fig. 6(a). The MSE curves of the SOR algorithm for Example 4 with different values of p.

/e

Fig. 5(c). The MSE curves of the SOR algorithm for Example 3 with/1 =0.8 and/~ =1.8.

I

1.0

0.8

taken to be 40 dB. The above SOR-based algorithm is applied to the received data series y(k), where we set r l = - 1 0 , r2 = 10 and / t = 1.3. The mean-square-error (MSE) curve is drawn in Fig. 3(a), and the error rate (ER) curve in Fig. 3(b). For the convenience of comparison, the simulation of the LMS-based version of our algorithm (i.e., (3.12a-c)) is also presented, and the corresponding results are also illustrated in Fig. 3(a) and Fig. 3(b), where we similarly set rl = - 1 0 and r2 = 10, but r(k) = 0.5/tr[ UT(k) U(k)]. Notice that/~ (k) has been chosen to make the LMS algorithm achieve

Qa

I ' ~,=1.0 I I---- ~=1.5 I

. . . . . -i

r--I

I

L2-

o.6-

I

i ii ~ 0.4-

--i ..... I_. . . . .

0.20.0

0.100

.........

i

,

,

2.103

4.103

6.103

8.103

1.10"

k

Fig. 6(b). The error rate curves of the SOR algorithm for Example 4 with different values of/1. Vol. 31, N o . 3, A p r i l 1993

F.-C. Zheng et al. / Equalisation o f multilevel P A M data

324

In the similar way to Example 1, the obtained results are illustrated in Fig. 4(a) and Fig. 4(b), respectively. This example demonstrates further the performance of the SOR algorithm much more.

H ( z ) = 1 + 3z -1 + 2z -2 + z -3.

(from [5]) H ( z ) = - 0 . 7 + 1.63z -t - 0.9z -2.

(6.5)

Let r ~ = - 1 0 , r2=10 and S N R = 3 0 d B . The transmitted data is the same as in Example 3: a two-level PAM series. The SOR algorithm is implemented using different relaxation factor values: /~ =0.5, 1.0, 1.5 and 1.9. The corresponding MSE and ER curves are shown in Fig. 6(a) and Fig. 6(b), respectively. The robustness of our algorithm to/~ can be seen again.

3. The channel transfer function is

EXAMPLE

4. The channel transfer function is

EXAMPLE

(6.4)

As in [5], we set the following conditions: r~ = - 1 0 , r2 = 30, SNR = 20 dB, and the transmitted data is a two-level PAM series. Then, the SOR algorithm is used exactly as in Examples 1 and 2. When p is 1.2, the convergence traces of 0j, j = + 1, +2, +3, are drawn in Fig. 5(a), and the corresponding MSE curve in Fig. 5(b). In comparison, the results here are at least as good as that in [5], although the form of the algorithm here is very simple and regular. Another advantage of the algorithm here is that its convergence rate is much less sensitive to/1 than the LMS scheme to r. To see this point, we draw the MSE curves of/.t =0.8 and/1 = 1.8 in Fig. 5(c). Clearly, there is no vital difference between them.

6.2. L T V channels

It should be noted that an exponentially windowed form of the algorithm has been used in these examples. This is an obvious step and the derivation is omitted since such an approach is commonly used. It is also possible to derive a rectangular window form of the algorithm.

EXAMPLE

5. The slowly variant channel with the

transfer function [5]

A. 2.0

*

1.0

*

:

:

:

:

:

:-

4.0

•

0.0

-

:

:

:

;

:

:

~

i

~

:

:

:

4.0

-

0.0.

8.0

:

:

:

:

0.0

:

:

:

:

:

:

:

:

:

"

0.~

0.0

.

.

.

~0

,

:

;

;

;

4.0

:

:

~

"

.

.

.

.

.

conver~n~

.

.

4.0

tra~s

6.0]

.

.

.

.

.

~.50 . . . . .

8.0

.

0.0!

8.0

:

:

:

;

0.0

.

.

.

.

:

:

:

:

0.0

inverse

i

;

;

;

:

4.0

8. 0

.

k/10 3 :

:

:"

0.I

~

t

~

:

:

:

:

:

:

:

:

:

:

:

"

. . . . . . .

,

. . . . .

4.0

-0.1.

8.0

:

;

0.0

filter

coefficients

:

:

~

4.0

k/l~ of the

.

0.0

k/103 The

:

~.15

-0.5 .

7(a).

:

k/10 3

~

Fig.

:

3.0

k/10 3

Signal P r o c e s s i n g

:

2.0

0.0

0.5

:

.

8.0

k/103 in Example

5 ~r

the

SOR

algorithm

with

~ = 1.4.

F.-C. Zheng et al. / Equalisation of multileveI PAM data ,

0

I

,

I

,

I

,

and

B.

-10-

The transmitted data is a two-level PAM series as in Examples 3 and 4. First, set SNR = 30 dB, rl = - 1 0 , r2 = 15 and p = 1.4. Then we apply the algorithm to the received data, where (3.19a-c) and (3.20a c) are used, and the forgetting factor f is taken to be 0.9999. The obtained converging traces of Oj,j=+l, +2, +3, are drawn in Fig. 7(a), and the MSE curve in Fig. 7(b). It is clear that the algorithm tracks the change of the channel very satisfactorily.

~" -30 -40 -

-50 -60-70 20100

4000

6 0 0i 0

800~

k

Fig. 7(b). The MSE curve of the SOR algorithm for Example 5 w i t h p = 1.4.

b22-2,

H ( z ) = bo + blZ -1 +

(6.7c)

bl = 1 + bob2.

-20-

0

325

6. The suddenly changed channel: The transfer function [5] is EXAMPLE

H(z)_S-0.7+~ 1.2lz -1-0.3z -2, - [ - 0 . 5 + l . 2 5 z - l - 0 . 5 z -2,

(6.6)

0~~4000.

where

(6.8)

b0 = -0.7 + 0.000025(k - 1),

(6.7a)

b2 = -0.3 - 0.000025(k- 1)

(6.7b)

The transmitted data is also a two-level PAM series. Again, set S N R = 3 0 d B , r l = - 1 0 , r2=15 and p = 1.1. Equations (3.21a c) and (3.22a-c) are

A 3.0

:

-"

:

:

:

:

:

:

:

"

?

:

5.0

:

:

:

:

:

:

:

:

"

7.0,

:

:

-"

-"

:

.

.

-"

:

:

:

"

2 1.5 ~

2.5

. . . . . . . . . .

0.0

: . : : ; : : : : .

0,0

4.0

0.0

$.0

0.0

k/103 0.6

,

,

:

3.5

4.0

0.0 . . . . .

8.0

.

, .....

0.0

4.0

8.0

k/103

k/103

: : : : : ,

.

"

0.2

:

:

:

:

:

r

=

,

,

41.1

0.0

-0.20'0~ . . . . . . . . . ,

-6.0

0.0

,

,

:

;

:

4.0

k/103

:

:

,

8.0

-0.5 : . . : • : : j , 0.0 4.0 k/103

8.0

0.0

4.0

8.0

k/103

Fig. 8(a). The MSE curve of the SOR algorithm for Example 6 with p = 1.1. VoL 31, No. 3, April 1993

326

I

[

F.-C. Zheng et aL / Equalisation of muhilevel P A M data I

B. -10 -15

-20~ -25-30-

~,

-35-

-40 0

,

2000

,

4000 k

,

6000

8000

Fig. 8(b). The convergence traces of the inverse filter coefficients in Example 6 for the SOR algorithm with/z = 1.1.

used to update the cumulants, where the window length is set to be 2000. The converging traces of O j , j = +1, +2, ±3, are illustrated in Fig. 8(a), and the MSE curve in Fig. 8(b). It can be seen that the change of the channel is closely tracked.

7. Conclusions

A new higher-order cumulant based equalisation algorithm has been described. In the algorithm, only the diagonal slices of cumulants are used. Several conclusions can be drawn: (a) The algorithm is simple in form, and can be easily programmed and implemented. (b) The introduction of the SOR scheme greatly accelerates the convergence rate and improves the adaptive performance of the algorithm. The algorithm also becomes insensitive to the selection of the step-size, viz., the relaxation factor. (c) The algorithm can deal with both LTI and LTV NMP channels. The simulation results confirmed the feasibility and efficiency of the algorithm.

References [1] A. Benveniste and M. Goursat, "Blind equalizers", 1EEE Trans. Comm., Vol. COM-32, August 1984, pp. 871 883.

SignalProcessing

[2] A. Benveniste, M. Goursat and G. Ruget, "Robust identification of a non-minimum phase system: Blind adjustment of a linear equalizer in data communications", IEEE Trans. Automat. Control, Vol. AC-25, June 1980, pp. 385 399. [3] S. Bellini et al., "Blind deconvolution: Polyspectrum or bussgang techniques?", in: E. Biglieri and G. Prati, eds., Digital Communications, Elsevier (North-Holland), Amsterdam, 1986, pp. 251 263. [4] D.R. Brillinger, "An introduction to polyspectra", Ann. Math. Statist., Vol. 36, 1965, pp. 1351 1374. [5] H.-H. Chiang and C.L. Nikias, "Adaptive deconvolution and identification of nonminimum phase FIR systems based on cumulants", IEEE Trans. Automat. Control, Vol. AC-35, January 1990, pp. 36 47. [6] B. Friedlander, "Adaptive algorithms for finite impulse response filters", in: C.F.N. Cowan and P.M. Grant, eds., Adaptive Filters, Prentice Hall, Englewood Cliffs, N J, 1985, Chapter 3. [7] B. Friedlander and B. Porat, "Adaptive IIR algorithms based on higher-order statistics", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-37, April 1989, pp. 485 495. [8] W. Gardner, "Learning characteristics of stochasticgradient-descent algorithms: A general study, analysis, and critique", Signal Processing, Vol. 6, No. 2, April 1984, pp. 113-133. [9] V.A. Gholkar, "Mean square convergence analysis of LMS algorithm", Electron. Lett., Vol. 26, No. 20, September 1990, pp. 1705-1706. [10] G.B. Giannakis and J.M. Mendel, "Identification of nonminimum phase system using higher order statistics", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP37, March 1989, pp. 360-377. [11] G.B. Giannakis and J.M. Mendel, "Cumulant-based order determination of non-Gaussian ARMA models", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, August 1990, pp. 1411-1423. [12] D.N. Godard, "Self-recovering equalization and carrier tracking in two-dimensional data communication systems", IEEE Trans. Comm., Vol. COM-28, November 1980, pp. 1867 1875. [13] D. Hatzinakos and C.L. Nikias, "Polyspectral techniques for blind equalization of multilevel PAM schemes", Proc. IEEE ICC'90, 1990, pp. 1507 1511. [14] D. Hatzinakos and C.L. Nikias, "Blind equalization based on second and fourth order statistics", Proc. 1EEE 1CC'90, 1990, pp. 1512 1516. [15] D. Hatzinakos and C.L. Nikias, "Blind equalisation using a trispectrum-based algorithm", IEEE Trans. Comm., Vol. COM-39, May 1991, pp. 669-682. [16] L. Ljung, "Consistency of the least squares identification method", IEEE Trans. Automat. Control, Vol. AC-22, October 1977, pp. 539 551. [17] O. Macchi and A. Hachicha, "Self-adaptive equalization based on a prediction principle", Proc. IEEE Globecom, 1986, pp. 1641 1645. [18] C.L. Nikias and H.-H. Chiang, "Higher-order spectrum estimation via noncausal autoregressive modeling and

F.-C. Zheng et al. / Equalisation of multilevel P A M data

[19]

[20]

[21]

[22]

[23] [24]

[25]

deconvolution", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-36, December 1988, pp. 1911 1913. C.L. Nikias and R. Pan, "ARMA modeling of fourthorder cumulants and phase estimation", Circuit Syst. Signal Process., Vol. 7, No. 3, 1988, pp. 291 325. C.L. Nikias and M.R. Raghuveer, "Bispectrum estimation : A digital signal processing framework", Proc. IEEE, Vol. 75, July 1987, pp. 869 891. C.L. Nikias and A.N. Venetsanopoulos, "Identification of nonminimum phase communication channels via parametric modelling of third moments", Proc. IEEE ICC'86, 1986, pp. 667 670. Y. Sato, "A method of self-recovering equalization for multilevel amplitude-modulation systems", IEEE Trans. Comm., Vol. COM-23, June 1975, pp. 679 682. R.S. Varga, Matrix lterative Analysis, Prentice Hall, London, 1962. K.Y. Wong and E. Polak, "Identification of linear discrete time systems using the instrumental variable method", IEEE Trans. Automat. Control, Vol. AC-12, December 1967, pp. 707 718. D.M. Young, Iterative Solution o f Large Linear Systems, Academic Press, New York, 1971.

327

[26] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind deconvolution algorithms based on 3rd- and 4th-order cumulants", Proc. IEEE Internat. Conf. Acoust. Speech Signal Process. (ICASSP), Toronto, Ontario, Canada, 14 17 May 1991, pp. 1753 1756. [27] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind equalisation of PAM series via higher-order cumulant fitting", Proc. IEEE Internat. Conf. on Communications (ICC), Denver, Colorado, USA, 23 26 June 1991, pp. 1393 1397. [28] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "A 2nd- and 4th-order cumulant based blind equalisation algorithm for nonminimum phase channels", Proc. IEEE & EURASIP International Signal Processing Workshop on Higher Order Statistics, Chamrousse, France, 10 12 July 1991, pp. 79 82. [29] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Cumulantbased deconvolution and identification: Several new families of linear equations", Signal Processing, Vol. 30, No. 2, January 1993, pp. 199 219. [30] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind equalization of nonminimum phase channels: Higher order cumulant based algorithm", IEEE Trans. Signal Process., to appear.

VoW.31, No. 3, April 1993