A novel kurtosis driven variable step-size adaptive algorithm

864

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

Fig. 3. Norm of coefficient-error vector for the SOV filter (white Gaussian input of variance 1, SNR 20 and 40 dB).

=

to estimate the SOV system (treated in [9] and [11])

yk =

0 0:78xk01 0 1:48xk02 + 1:39xk03 + 0:04xk04 + 0:54xk201 + 3:72xk01 xk02 + 1:86xk01 xk03 0 0:76xk01xk04 0 1:62x2k02 + 0:76xk02xk03 0 0:12xk02xk04 + 1:41x2k03 0 1:52xk03xk04 0 0:13x2k04 + nk : We then have N = 4; M = 5;  = 2M = 10; n = 14. Fig. 3

illustrates the evolution of the performance criterion for two values of the SNR. We note that no numerical stability problem occurred during the simulations.

[6] F. Fnaiech, S. Guillon, and M. Najim, “A fast M-D recursive least square adaptive second-order Volterra filter,” in Proc. IEEE Workshop Nonlinear Signal Image Process., Neos-Marmaras, Greece, June 22–24, 1995, pp. 408–411. [7] A. Houacine, “Regularized fast recursive least squares algorithms for adaptive filtering,” IEEE Trans. Signal Processing, vol. 39, pp. 860–870, Apr. 1991. [8] B. H. Khalaj, A. H. Sayed, and T. Kailath, “A unified derivation of square-root multichannel least-squares filtering algorithms,” in Proc. IEEE ICASSP, 1993, pp. 523–526. [9] J. Lee and V. J. Mathews, “A fast recursive least square adaptive secondorder Volterra filter and its performance analysis,” IEEE Trans. Signal Processing, vol. 41, pp. 1087–1102, Mar. 1993. [10] L. Ljung and T. Kailath, “Efficient change of initial conditions, dual Chandrasekhar equations, and some applications,” IEEE Trans. Automat. Contr., vol. AC-22, pp. 443–447, June 1977. [11] V. J. Mathews, “Adaptive polynomial filters,” IEEE Signal Processing Mag., pp. 10–26, July 1991. [12] M. Morf, G. S. Sidhu, and T. Kailath, “Some new algorithms for recursive estimation in constant, linear systems,” IEEE Trans. Automat. Contr., vol. AC-19, pp. 315–323, Aug. 1974. [13] M. Sayadi, A. Chaari, F. Fnaiech, and M. Najim, “A fast M-D Chandrasekhar algorithm for second-order Volterra adaptive filtering,” in Proc. IEEE ICASSP, Atlanta, GA, May 7–10, 1996, pp. 1339–1342. [14] A. H. Sayed and T. Kailath, “Extended Chandrasekhar recursions,” IEEE Trans. Automat. Contr., vol. 39, pp. 619–620, Mar. 1994. [15] G. L. Sicuranza, “Quadratic filters for signal processing,” Proc. IEEE, vol. 80, pp. 1263–1285, Aug. 1992. [16] D. T. M. Slock and T. Kailath, “Numerically stable fast transversal filters for recursive least squares adaptive filtering,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 39, pp. 92–114, Jan. 1991.

VIII. CONCLUDING REMARKS In this work, we have presented a fast algorithm for the multichannel linear and SOV filtering problem based on fast Chandrasekhar equations. This technique is extended from the single-channel structure and could handle a different number of delay elements per input channel. The method has good numerical properties since, instead of propagating through the recursions, the overall covariance matrix, which sometimes loses its positive definiteness property, propagates a factorized increment of this matrix, which has a reduced rank. The adaptive filter was successfully used in multichannel linear and nonlinear system identification, and its convergence and numerical stability are illustrated by simulations. ACKNOWLEDGMENT The authors wish to acknowledge the anonymous referees for their valuable comments and constructive suggestions. REFERENCES [1] M. G. Bellanger, Adaptive Digital Filter and Signal Analysis. New York: Marcel Dekker, 1987. [2] J. Cioffi and T. Kailath, “Fast RLS transversal filters for adaptive filtering,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP32, Apr. 1984. [3] G. Demoment and R. Reynaud, “Fast RLS adaptive algorithms and Chandrasekhar equations,” in Proc. SPIE Adaptive Signal Process., S. Haykin, Ed., July 1991, vol. P-1565, pp. 357–367. [4] X. C. Du, A. Brella, and R. Longchamp, “Fast algorithm of Chandrasekhar type for ARMA model identification,” Automatica, vol. 24, no. 5, pp. 927–934, 1992. [5] D. D. Falconer and L. Ljung, “Application of fast Kalman estimation to adaptive equalization,” IEEE Trans. Commun., vol. COMM-26, pp. 1436–1446, Oct. 1978.

A Novel Kurtosis Driven Variable Step-Size Adaptive Algorithm Dimitrios I. Pazaitis and A. G. Constantinides

Abstract—In this correspondence, a new variable step-size LMS filter is introduced. The time-varying step-size sequence is adjusted, utilizing the kurtosis of the estimation error, therefore reducing performance degradation due to the existence of strong noise. The convergence properties of the algorithm are analyzed, and an adaptive kurtosis estimator that takes into account noise statistics and optimally adapts itself is also presented. Simulation results confirm the algorithm’s improved performance and flexibility.

I. INTRODUCTION Stochastic gradient adaptive filters using the least mean square (LMS) algorithm [12] enjoy great popularity due to their inherent simplicity and robustness. The LMS algorithm has found applications in numerous and diverse fields, such as noise canceling, adaptive line enhancing, telecommunications, geophysics, etc., and its performance has been thoroughly investigated in the literature (e.g., [3], [5], [12]). Manuscript received March 31, 1996; revised May 1, 1998. This work was supported by a scholarship from the State Scholarship Foundation (I.K.Y) of the Hellenic Republic. The associate editor coordinating the review of this paper and approving it for publication was Dr. Akihiko Sugiyama. D. I. Pazaitis is with IMEC-VSDM, Leuven, Belgium (e-mail: [email protected]). A. G. Constantinides is with the Signal Processing Section, Department Electrical and Electronic Engineering, Imperial College, London, U.K, (e-mail: [email protected]). Publisher Item Identifier S 1053-587X(99)01351-3.

1053–587X/99$10.00  1999 IEEE

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

865

independent of xn . The identification (estimation) error at time n is denoted by e(n). All signals are real and zero mean, and n is the discrete time index. The adaptive FIR filter coefficients vector is denoted by n ; n = [hn ; . . . ; hn0N +1 ]T . The vector containing the parameters of the unknown system at time n; opt n is assumed to have the same dimensions as n . The error signal is then defined as

H H

H

H

wn 0 XTn Vn

e(n) = en

V

Fig. 1. System identification setup.

A constant step-size , which is also known as convergence factor, governs the stability of the algorithm, as well as the rate of convergence and the steady-state excess mean square error in relation to the optimal Wiener solution. The convergence time is inversely proportional to , whereas the misadjustment is proportional to N , where N equals the number of the weights of the adaptive filter [12]. When a nonstationary environment is considered, the “lag misadjustment” is proportional to 01 . Consequently, the optimal selection of the step size is dictated by a tradeoff between convergence rate and steady-state excess error and has been frequently addressed in the literature (e.g., [1], [12]). To meet the conflicting requirements of fast convergence, good tracking ability, and low misadjustment, time-varying step-size sequences have been proposed (e.g., [4], [7]–[9]). The LMS+F algorithm [10] may also be regarded as a variable step-size algorithm. The main rationale behind these approaches is to somehow sense the distance from the optimum and correspondingly adapt the stepsize value. Various approaches have been proposed based—among others—on the magnitude of the squared error (e.g., [8], [9]) and the gradient of the error surface (e.g., [4], [7]). It seems to us, however, that the performance of most of the proposed algorithms is susceptible to degradation in the presence of strong noise. Furthermore, very fine tuning of the algorithm parameters is usually required to achieve satisfactory performance. In this correspondence, a new time-varying step-size selection method is proposed that provides noise immunity under Gaussian noise conditions. To achieve this, we resort to the use of higher order statistics and, more precisely, to the fourth-order cumulant, which is known as kurtosis. The reason for using cumulants is twofold. Cumulants are additive in their arguments, i.e., cumulants of sums equal sums of cumulants (hence, the name cumulant), and second, all cumulants are blind to any kind of a Gaussian process, i.e., all cumulants of order greater than two of a Gaussian process equal zero. Thus, by introducing cumulants, we immunize in a way the step-size selection algorithm against Gaussian or Gaussian-like noise. In a real environment, even if the noise is not Gaussian, the existence of a large number of additive interferences “moves” the probability density function of the noise—according to the central limit theorem—toward Gaussianity. Next, we formulate the adaptive system identification problem and introduce our variable step-size LMS algorithm. II. A VARIABLE STEP-SIZE LMS ALGORITHM The structure of the system identification setup is depicted in Fig. 1. The input signal is denoted as xn , and wn is the observation noise. We assume that wn is symmetrically distributed and statistically

H H

(1)

opt where is the weight error vector. The n = n0 n symbol ( )T stands for the vector transposition operation n = T [xn ; . . . ; xn0N +1 ] is the input vector, and N is the number of adaptive coefficients. The popular LMS adaptive algorithm is a steepest descent gradient search type algorithm that uses, at each iteration, the instantaneous value of the gradient and seeks to minimize the mean square error E fe2 (n)g. The weight vector is updated using

Hn

+1

=

X

Hn + nXn e(n)

(2)

where n is the time-varying step size. In its original form [12], the step-size assumes a constant value. In our proposed algorithm, the step size is a function of the kurtosis of the error. The following two types of functions are proposed for adjusting the step size of the LMS algorithm

n =  C4e (n)

(3)

and

n = max

1

0 e0jC

nj

( )

(4)

where

C4e (n) = E fe4 (n)g 0 3E 2 fe2 (n)g

(5)

is the kurtosis of the error,  is a positive constant, and j 1 j denotes the absolute value. The constant max in (4) can be chosen as the maximum allowable value of n that ensures convergence. Each of the above functions exhibits certain advantages. Update (4) is selfrestricted to values lower than or equal to max , requiring, therefore, no supervision. On the other hand, update (3) is simpler to use, but to ensure convergence, an upper bound max on n may have to be imposed. It is also noted that the constant max in update (4) not only represents the maximum allowable value but also scales the term in the parenthesis. Thus, in certain environments, values smaller than the stability bound may prove more beneficial. If desired, a lower bound on the step-size sequence may be applied in order to ensure a minimum level of tracking speed. From (1), we see that the kurtosis of the error equals the sum of two terms: the noise kurtosis and the kurtosis of the convolution of the weight error vector with the input signal. Under the Gaussian noise assumption, the first term equals zero, whereas the second one approaches zero as the filter converges and the weight error vector n tends to zero. Thus, the algorithm starts with a large value for n , which decreases as it approaches the optimum, permitting, theoretically, exact convergence to the optimal solution. A question that might arise is whether it would be more beneficial to use the third-order cumulant (skewness) of the error signal. The reason for preferring the fourth-order cumulant (kurtosis) lies in the fact that the third-order cumulant of a symmetric distribution always equals zero. Thus, under certain conditions, the kurtosis of the second term in (1) might be zero, therefore freezing the adaptation process of the filter weights. Furthermore, experimental results will confirm the

V

866

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

improved performance characteristics of the kurtosis-driven step-size compared with step-size updates based on lower than the fourth-order statistics. The use of the exponential function to control n was first proposed by Karni and Zeng in [7], where the norm of the instant gradient, rather than the kurtosis, was used. Their algorithm can be interpreted as a new type of algorithm, which considers a performance function that includes an additional exponential function term in comparison with that of LMS [6].

R

is symmetric, it may be rotated Since the input covariance matrix into a diagonal matrix by a unitary transformation

R = Q3Q01 = Q3QT Q

In this section, the behavior of the proposed variable step-size adaptive algorithm is investigated. Two types of convergence are examined, namely, convergence in the mean and convergence in the mean square sense. The first type of convergence suffices for the mean weight error to equal zero, whereas the second one ensures a steady-state error with finite variance. To facilitate the analysis, a number of assumptions and approximations are introduced. The most common assumption is that the various input vectors come from mutually independent, zero mean, Gaussiandistributed sequences [3], [5], [12]. Although this is hardly true, especially in tapped-delay filters, where consecutive input vectors share N 0 1 entries, it is extensively used in the literature to simplify the analysis and has produced reliable results. The unknown system (optimal weight vector) opt is considered time varying, and a random walk model is adopted to model its evolution

H

=

Hn

opt

+ A(n)

(6)

where A(n) is a white, zero mean vector disturbance process with 2 a covariance matrix equal to A . We assume that the triplet f n; A(n); wng are statistically independent random processes, and we also introduce the assumption that the step-size sequence n is statistically independent of n ; n and e(n) [8], [9]. In the ideal version of the algorithm, where the actual kurtosis is utilized, this assumption is valid. In practical applications, however, where the error sequence fe(n)g is considered ergodic and a finite length window is used to estimate the kurtosis, this assumption does not hold. Experimental results have shown, nevertheless, that under slow convergence conditions or by increasing the time span of the window, the assumption tends to hold. For the algorithm to converge in the mean sense, the mean weight error vector should converge to the zero vector. A sufficient condition for this is [12]

X H

E fn g <

2

(7)

max

R

where max is the largest eigenvalue of the correlation matrix of the input signal sequence. We proceed now to investigate the evolution of the mean square T , taking expectations, and error. Forming the product n+1 n +1 using the step-size independence assumption leads to

V V

E fVn+1 VnT+1 g T T T = E Vn Vn 0 E Vn Vn E fn gE Xn Xn 0 E fngE Xn XnT E Vn VnT 2 T T T + E n E Xn Xn Vn Vn Xn Xn T 2 T 2 + E fA(n)A (n)g + E n E Xn Xn E wn :

R

Vn0 1 Vn0T 1 0 0T 0 E V0 V0T E fn g3 = E Vn Vn n n 0 E fng3E Vn0 Vn0T 0 0T 0 0T 0 0T + E fA0 (n)A0T (n)g + E n E Xn Xn Vn Vn Xn Xn + E n 3E wn (10) +

+

2 2

2

0 = T n ; n0 = T n , and n0 = T n . where n Bershad and Qu [2] have shown that the weights of the adaptive filter are jointly Gaussian. Here, the components of n are assumed to be independent Gaussian random variables [5]. Under this assumption, the fourth term on the right hand of equation (10) can be written as [5]

V

Q V X

A

Q X

Xn0 Xn0T Vn0 Vn0T Xn0 Xn0T

E

Q A V

3E Vn0 Vn0T 3 0 0T + 3 tr 3E Vn Vn

=2

(11)

where tr[ 1 ] stands for the trace of the bracketed matrix. If we let n be a vector constituted by the diagonal terms of the 0T g and be a column vector of ones of the same matrix E f n0 n length as n , then from (10) and (11), we have

G 3 VV G

I

X

3

Q

H

opt n+1

R

where is the orthonormal modal matrix of , and = diag[1 ; 2 ; . . . n ] is the diagonal matrix of eigenvalues of . Thus, in rotated (primed) coordinates, premultiplying all vectors by T leads to

E

III. CONVERGENCE ANALYSIS OF THE PROPOSED ALGORITHM

(9)

Gn

+1

=

1

I 0 23E fng + 32 E n (2I + 11T ) Gn 2 + E n 3 1w + 31A : 2

2

2

2

(12)

The mean and mean square value of n are given by

E fn g = E  C4e (n) E 2n = E  C4e (n)

= 2

jE fe4 (n)g 0 3E 2 fe2 (n)gj 2 4 2 2 2 =  jE f(e (n)g 0 3E fe (n)gj (13)

and

E fn g = E max 1 0 e0jC (n)j 0jE fe (n)g03E fe (n)gj = max 1 0 e E 2n = E 2max 1 0 e0jC (n)j 2 2 0jEfe (n)g03E fe (n)gj = max 1 0 e

(14) 2

for the first (3) and second (4) type of function, respectively. It is observed that the mean square value of the step size equals the square of its mean value. Using its definition and the independence assumption, it can be shown that

E fe2 (n)g = w2 + tr 3E

Vn0 Vn0T

=

w2 + 1T Gn :

(15)

From the above equation, it is seen that the excess mean squared error at each time instant is given by (8)

1 G

excess mse = T n :

(16)

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

An approximate expression for the fourth-order error moment is derived in [8]

E fe4 (n)g = 3w4 + 6w2 tr 3E Vn0 Vn0T 0 0T 2 + 6 tr E 2 3V0 V0T + 3 tr E 3Vn Vn n n 4 0 0 6 tr E f(I 0 n01 3) gE fVn01 gE Vn0T01 32 : (17)

867

IV. A SELF-ADJUSTED KURTOSIS ESTIMATION The estimation of the kurtosis using (5) limits the usefulness of the algorithm to Gaussian noise conditions. To overcome this problem and extend the applicability of the algorithm, a modified version is proposed. The modified algorithm is more flexible in that it can track down changes in noise statistics and adapt itself accordingly to retain the desirable performance characteristics of fast convergence and low steady-state error, regardless of the type of noise encountered. The proposed error kurtosis estimator has the form

C4e (n) = E en4

Equations (12) to (17) completely describe the mean squared error behavior. Now, by setting

F = I 0 23E fng + 32 E

2n

I 11T )

(2 +

(18)

F

the convergence of (12) depends on matrix . Its convergence is guaranteed (exponential stability) if and only if the eigenvalues of lie within the unit circle. The eigenvalues are given by the solution in det[ 0  ] = 0. Proceeding with the algebraic manipulations and adopting the approach in [3] and [8], it can be shown that the necessary and sufficient condition for convergence is

F

F

Under this condition,

 3 tr1[R] :

Using the fact that becomes

2

R] :

(21)

In theory and in stationary environments, the algorithm reduces its step size as it approaches the optimum set of coefficients, permitting the exact identification of the unknown system. It can be shown, in fact, that by using the above equations, the step size n tends to zero, thus eliminating the excess mean squared error. In practice, however, the real kurtosis of the error signal is unknown and has to be estimated. This is usually done by assuming ergodicity of the error sequence and replacing ensemble averaging by time averaging over a finite window. In reality, therefore, the estimated kurtosis is a random variable and its steady-state variance is a function of the window length. Consequently, the step size sequence n , as a function of the absolute value of the kurtosis, reaches a steady-state value that is different from zero. Its distance from zero decreases asymptotically, as the window length increases to infinity. Nevertheless, the condition for convergence given by (19) still holds. Substituting the values for n given by (4) and expanding the square results in

max

@ C4e (n) @n 4 = n +  sgn E en

1+E

e0jC (n)j 2 0 2E e0jC (n)j 2 1 0 E e0jC (n)j

 3 tr[1R] : (22)

It can be readily seen that the right-hand term is always lower than or equal to unity. Thus, by constraining max to satisfy

max 

1

R]

3 tr[

convergence can be always guaranteed.

0 n E en E en 2

2

2

2

(25)

with  being a scaling factor. A more simplified yet accurate update is given by

n+1 = n +  sgn E en4

0 n E en : 2

2

(26)

V. SIMULATION RESULTS (20)

E fn2 g = E 2 fn g = 2n , the above condition 3 tr[

(24)

The fact that the updates (25) and (26) are self-driven, requiring, therefore, no a priori knowledge of the noise statistics, renders them a very useful tool in unknown or time-varying signal processing environments.

Gn converges to [3]

E fn g = n 

2

n+1 = n +

(19)

G1 = (I 0 F)01 1 A2 1 + E 2n 31w2 2 T 01 = 2E fn gI 0 3E n (2I + 11 ) 2 A2 1 + E 2n 31w2 :

2

where n is a time-varying sequence aiming to trace the changes in the noise characteristics and adapt itself accordingly. Adjustment of the parameter n is carried out according to

I

E 2n 2E fn g

0 n E en

(23)

In this section, we present and analyze the results obtained from simulations. The algorithm is applied in a system identification setup, where the system to be identified is considered nonstationary. Its coefficients assume the initial values

H

opt

= [1:0; 0:268; 0:083;

00:501; 0:941; 0:011;

00:819; 00:376; 0:701; 0:159]

0

and, after that, experience random disturbances. The null vector ( ) is chosen as the initial vector 0 , and all the results are obtained by averaging over an ensemble of 200 runs. Both noise and input sequences are assumed to be zero mean i.i.d. Gaussian sequences, with variances equal to 0.1 and unity, respectively. The kurtosis of the error C4e (n) is estimated by using a moving finite length rectangular window of size M = 50. Alternatively, the second- and the fourth-order error moments can be estimated by using

H

E en2i

=

E en2i01

+ (1

0 )eni; i = 1; 2: 2

(27)

The window length and the forgetting factor control the memory of the system and should therefore be selected according to the application needs. The selection of the parameter  is dictated by a tradeoff between convergence speed and steady-state misadjustment for both types of functions. In nonstationary environments, where abrupt changes are frequently encountered, it is advisable to choose a high value for  to enable the algorithm to quickly track down these changes. In stationary environments, however, small values of  are proven to be more beneficial. As a rule of thumb, in the exponential case, a value in the range of [1] and [10] results in fast response to system changes and low steady-state error. For noise distributions with long tails (e.g., Laplacian), it is advisable to choose a small value in this interval, whereas for noise distributions with short tails (e.g., uniform), larger values are suggested. For similar performance characteristics, the value of  in the linear update has to be chosen approximately an

868

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

Fig. 2.

Misadjustment comparison between the fixed step-size LMS and the proposed variable step size algorithms.

order of magnitude smaller than the value in exponential one. A detailed parameter investigation can be found in [11]. Fig. 2 depicts the behavior of the proposed variable step-size algorithm in comparison with the fixed step-size LMS algorithm for both types of functions, i.e., the linear and the exponential one. The parameters corresponding to the plotted curves are max = 0:07; LMS = 0:048;  = 5:0 for the exponential and  = 0:3 for the linear function. The above parameters were chosen so that the algorithms exhibit similar convergence rates. The sudden changes that happen in the system are due to zero mean uniformly distributed 2 disturbances with variance equal to A = 0:1 (at random the first time but fixed for the following runs) time instants. The signal-to-noise ratio (SNR) was chosen equal to 10 dB, and the system mismatch T (E f n n g) was taken as the performance measure. The performance of the adaptive step-size algorithms is clearly superior to that of the fixed step-size LMS, which is outperformed by more than 10 dB. This is due to the time-varying nature of n , which allows the algorithm to increase its step size and, therefore, its convergence speed whenever it is far from the optimum. It is also observed that the two types of variable step-size algorithms exhibit almost identical behavior, rendering choice a matter of convenience. To gain better insight in the behavior of the algorithm, the evolution of the step-size sequence, n is depicted in Fig. 3 for the exponential update. The step-size sequence corresponding to the linear update behaves similarly. The theoretically expected step-size evolution is also depicted in the same plot, and it is observed the simulation results are in good agreement with the theoretic ones. To obtain the theoretical results, the independence assumption together with the fact that n = E fn g and, thus, E 2 fn g = 2n , were used. In practice, however, neither of the above holds. The exact values of the kurtosis are unknown and have to be estimated by using a moving finite length window. Consequently, the estimator produces a random variable. Using the absolute value of this estimator to constrain n to positive values results in a steady-state value that is different from zero, whereas the theoretic curves continue to converge. We next proceed to compare the performance of the proposed algorithm with that of previously developed methods. In Fig. 4, the proposed algorithm is compared with two other variable step-size

VV

algorithms, which were developed by Harris et al. [4] and Mathews and Xie [9]. For the particular simulation, the parameters for Harris’s algorithm were chosen as  = 2:0; m0 = m1 = 2; max = 0:07; min = 1e08 , which are also the values recommended in [4]. For the algorithm proposed by Mathews and Xie, the parameters were chosen to match the tracking speed of the algorithms, that is,  = 0:0008; max = 0:1; 0 = 0:06. For our proposed algorithm, the exponential update was used with  = 5:0 and max = 0:07. The disturbances 2 of the system are uniformly distributed with variance A = 0:1. The improved performance of the proposed algorithm is clearly observed. Moreover, the proposed algorithm retains its desirable characteristics in higher SNR environments. This is shown in Fig. 5, where the performance comparison is carried out for an SNR value of 20 dB, and the algorithms’ parameters were chosen to result in similar convergence rates. In Fig. 6, the performance of the kurtosis-driven step-size algorithm compares favorably with the performance of algorithms using other, lower order statistics in the exponential term. The plotted algorithms use the mean absolute, the mean square, and the thirdorder moment (cumulant) of the error sequence, and the results provided correspond to an SNR equal to 10 dB. The algorithms’ parameters were selected to match the convergence rate. Thus,  = 5;  = 1;  = 2:3; and  = 4:0 for the absolute value, the second-, the third-, and the fourth-order cumulants, respectively. In all previous experiments, the noise sequence was assumed to be Gaussian distributed. We now consider the algorithm’s behavior under non-Gaussian or time-varying noise conditions. We confine ourselves to the exponential-type step-size update (4), but the observations and conclusions extend themselves to the linear-type update as well. In Fig. 7, the probability distribution of the noise sequence changes from Gaussian to Laplacian and then through uniform to binary and, finally, to a sinusoidal type. The sinusoidal-type noise is a sinusoid with a random phase. The kurtosis was estimated using (24) and (26) with  = 0:01, and the noise variance was kept equal to 0.1 in all cases. The parameter  was given by the empirical rule

=

10 1 3

n

;

3

n

;

if

n < 3

otherwise

(28)

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

869

Fig. 3. Evolution of the step-size sequence of the proposed adaptive step size algorithm and comparison with the theoretically expected curve.

Fig. 4. Performance comparison between various variable step-size algorithms (SNR 10 dB).

which was found to combine simplicity and satisfactory performance. For reasons of comparison, the curves corresponding to Mathews’ and Harris’ algorithms are also presented, and the improved robustness of the proposed algorithm is easily confirmed. To gain better insight in the performance of the algorithm, the evolution of the n sequence is depicted in Fig. 8. For the same purpose, the steady-state values of the

adaptive sequence n are compared with the theoretically expected ones, which correspond to pure noise signals, in Table I. The steadystate values of n were obtained by averaging over the last 5000 samples of each distribution interval. The simulation results confirm the ability of the kurtosis estimator to trace changes in the noise statistics and adapt itself optimally. They also confirm the domination

870

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

Fig. 5.

Performance comparison between various variable step-size algorithms (SNR

= 20 dB).

Fig. 6.

Performance comparison between various variable step size algorithms (SNR

= 10 dB).

of the noise signal in the steady-state error values. It is therefore justified to use the error signal in the kurtosis estimation process instead of the unavailable pure noise signal. The bias effect of the approximately Gaussian-distributed second error component in (1) is negligible in most cases. This fact supports the use of a constant parameter , provided that the noise statistics are stationary and known in advance. However, in unknown or time-varying environments, the

introduction of the adaptive kurtosis estimator significantly improves the performance characteristics of the algorithm at the expense of a small increase in the computational complexity. VI. CONCLUSION A new variable step-size LMS adaptive algorithm was proposed in this contribution. The algorithm differs from previous techniques

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

871

Fig. 7. Performance curves of variable step-size algorithms for time varying noise statistics.

Fig. 8. n sequence evolution of the algorithm shown in Fig. 7.

involving time-varying step-size sequences in that it utilizes the kurtosis of the estimation error signal to adjust its convergence factor, exhibiting, therefore, reduced sensitivity to Gaussian-type noise signals. Compared with the fixed step-size LMS, the proposed scheme exhibits a small increase in complexity. In the linear update case, this is equal to ten extra multiplications, five of which involve a

constant and are therefore easy to implement, and three additions per algorithmic iteration. In the exponential case, one extra addition and multiplication (with a constant), as well as the implementation of the exponential function, are required. The mean squared error behavior and stability conditions were provided, and simulation results have confirmed the superior tracking ability and lower misadjustment of the proposed algorithm com-

872

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999

TABLE I STEADY-STATE VALUES OF THE n SEQUENCE

An LMS Adaptive Second-Order Volterra Filter with a Zeroth-Order Term: Steady-State Performance Analysis in a Time-Varying Environment Mounir Sayadi, Farhat Fnaiech, and Mohamed Najim

pared with the fixed step-size LMS algorithm and other previously developed variable step-size algorithms. To extend the applicability of the algorithm to environments characterized by non-Gaussian or time-varying noise statistics, a timevarying kurtosis estimator was introduced, which traces the noise statistics and optimally adjusts itself. The flexibility of the algorithm and the ability of the adaptive kurtosis estimator to “tune” itself according to the noise signal statistics was also verified, rendering the proposed technique a very useful tool in adaptive signal processing applications.

Abstract— This correspondence studies the steady-state performance of the least mean square (LMS) adaptive second-order Volterra filter (SOVF) with a zeroth-order term for Gaussian inputs. The mean-squareerror (MSE) criterion is evaluated first. Then, SOV LMS algorithm-based updating equations are derived. Next, the steady-state performance of the recursions is analyzed for a random walk model for the unknown system parameters, and the steady-state excess MSE is evaluated. Finally, the theoretical performance predictions are shown to be in good agreement with simulation results, especially for small step sizes. Index Terms—Nonlinear adaptive filtering, performance analysis, random walk, steady-state excess mean square, steepest descent method, time-varying environment, Volterra filter.

I. INTRODUCTION

REFERENCES [1] N. J. Bershad, “On the optimum gain parameter in LMS adaptation,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-35, pp. 1065–1068, July 1987. [2] N. J. Bershad and L. Z. Qu, “On the probability density function of the LMS adaptive filter weights,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 37, pp. 43–56, Jan. 1989. [3] A. Feuer and E. Weinstein, “Convergence analysis of LMS filters with uncorrelated Gaussian data,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-33, pp. 222–229, Feb. 1985. [4] R. W. Harris, D. M. Chabries, and F. A. Bishop, “A variable step (VS) adaptive filter algorithm,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-34, pp. 309–316, Apr. 1986. [5] L. L. Horowitz and K. D. Senne, “Performance advantage of complex LMS for controlling narrow-band adaptive arrays,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, pp. 722–735, June 1981. [6] B.-E. Jun and D.-J. Park, “Novel steepest descent adaptive filter derived from new performance function with additional exponential term,” Signal Process., vol. 36, pp. 189–199, 1994. [7] S. Karni and G. Zeng, “A new convergence factor for adaptive filters,” IEEE Trans. Circuits Syst., vol. 36, pp. 1011–1012, July 1989. [8] R. H. Kwong and E. W. Johnston, “A variable step size LMS algorithm,” IEEE Trans. Signal Processing, vol. 40, pp. 1633–1642, July 1992. [9] V. J. Mathews and Z. Xie, “A stochastic gradient adaptive filter with gradient adaptive step size,” IEEE Trans. Signal Processing, vol. 41, pp. 2075–2087, June 1993. [10] D. I. Pazaitis and A. G. Constantinides, “LMS F algorithm,” IEE Electron. Lett., vol. 31, no. 17, pp. 1423–1424, Aug. 1995. [11] D. I. Pazaitis, “Performance improvement in adaptive signal processing algorithms,” Ph.D. dissertation, Imperial College, London, U.K., Oct. 1996. [12] B. Widrow and S. D. Stearns, Adaptive Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1985.

+

Due to implementation simplicity, LMS is the most widely used adaptive linear filter algorithm for both stationary and nonstationary environments. The LMS convergence and tracking properties have been studied when used for identifying nonstationary linear filters [4], [9]. The MSE was shown to involve a tradeoff between the gradient error (increasing with ) and the lag error (decreasing with ). For a nonlinear adaptive filter structure and specifically for SOVF [7], few papers deal with the convergence analysis of the LMS type adaptive algorithms. In [3], the authors studied the LMS SOVF convergence for the stationary case and evaluated the steady-state excess MSE. The convergence analysis for the LMS second-order Volterra adaptive filter has not been done in the nonstationary context until recently. Moreover, except in [1], zeroth-order term adaptation has not been considered in either stationary or nonstationary cases. In [1], the authors also showed that the zeroth-order term degrades the stability limit and slightly increases the steady-state MSE. However, no quantitative evaluation of the MSE degradation was made. Generally, the previous works assume no loss when the zero-order term is discarded in the second-order Volterra model. The desired filter output is assumed zero-mean in [3] and [6] for simplicity. The zeroth-order term is assumed known and not used in the updating process. Without the zeroth-order term, the desired output mean is needed. The adaptive zeroth-order term is useful when the desired output mean is unknown or is time varying. This correspondence’s main objective is to evaluate the steadystate excess MSE of the SOVF in the nonstationary context. This analysis studies zeroth-order term adaptation for the second-order Volterra model. Our contribution consists of evaluating the steadystate MSE corresponding to both the filter quadratic part and the zeroth-order term. In this correspondence, the time index is denoted in parentheses, the transpose of a matrix is denoted by the superscript t, vectors and matrices are indicated by underlined and capital Manuscript received January 31, 1997; revised March 5, 1998. The associate editor coordinating the review of this paper and approving it for publication was Dr. Mahmood R. Azaimi-Sadjadi. M. Sayadi and F. Fnaiech are with Ecole Sup´erieure des Sciences et Techniques de Tunis, Tunis, Tunisia (e-mail: [email protected]). M. Najim is with Equipe Signal Image, ENSERB, Talence, France (e-mail: [email protected]). Publisher Item Identifier S 1053-587X(99)01352-5.

1053–587X/99$10.00  1999 IEEE