Blind identification of bilinear systems

Abstract-This paper is concerned with the blind identification of a class of bilinear systems excited by non-Gaussian higher order white noise. The matrix of coefficients of mixed input-output terms of the bilinear system model is assumed to be triangular in this work. Under the additional assumption that the system output is corrupted by Gaussian measurement noise, we derive an exact parameter estimation procedure based on the output cumulants of orders up to four. Results of the simulation experiments presented in the paper demonstrate the validity and usefulness of our approach.

484 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51,NO. 2, FEBRUARY 2003 Blind Identification of Bilinear Systems Nicholas Kalouptsidis, Senior Member, IEEE, Panos Koukoulas, Member, IEEE, and V. John Mathews, Fellow, IEEE Abstract-This paper is concerned with the blind identification of a class of bilinear systems excited by non-Gaussian higher order white noise. The matrix of coefficients of mixed input-output terms of the bilinear system model is assumed to be triangular in this work. Under the additional assumption that the system output is corrupted by Gaussian measurement noise, we derive an exact pa­rameter estimation procedure based on the output cumulants of or­ders up to four. Results of the simulation experiments presented in the paper demonstrate the validity and usefulness of our approach. Index Terms-Bilinear systems, blind identification, high-order statistics, nonlinear system identification. I. INTRODUCTION IDENTIFICATION of nonlinear systems is of primary importance in today's applications since many signals of interest are generated by nonlinear sources or are processed by nonlinear systems. There are several situations in which the inherent nonlinearities and distortions cannot be tolerated at a given level of performance, and hence, nonlinear processing techniques need to be employed. Such important examples include nonlinear echo cancellation, predistortion of nonlinear channels, equalization of communication channels where distortion is produced due to operation of amplifiers near to saturation region, linearization of loudspeaker nonlinearities, enhancement of noisy images, edge extraction, distortions in magnetic recording systems, motion of moored ships in ocean waves, control of industrial processes, physiological models, nuclear fission, and others [1], [7], [12], [25], [29], [32]. Conventional identification is concerned with the determina­tion of an unknown system on the basis of input-output informa­tion in an uncertain environment. A given excitation drives the unknown system and the resulting response is measured. On the other hand, blind identification is concerned with the determi­nation of an unknown system on the basis of output information only. In this latter case, information about the input that gener­ates the measured output is limited. For instance, it may be a priori known or assumed that the input is white noise. The need for tractable computational methods requires that the class of nonlinear models is properly restricted. Polynomial systems form a popular class of nonlinear systems [32]. Under Manuscript received January 14, 2002; revised September 12, 2002. The associate editor coordinating the review of this paper and approving it for publication was Prof. Derong Li. N. Kalouptsidis is with the Department of Informatics and Telecommunica­tions, Division of Communications and Signal Processing, University of Athens, Athens, Greece (e-mail: [e-mail]). P. Koukoulas is with the Electronics Division, Department of Radars and Automations, Hellenic Civil Aviation Authority, Ministry of Transport and Communications, Athens, Greece (e-mail: [e-mail]). V. J. Mathews is with the Department of Electrical and Computer En­gineering, University of Utah, Salt Lake City, Utah 84112 USA (e-mail: [e-mail]). Digital Object Identifier 10.1109/TSP.2002.806971 relatively mild conditions, such systems are known to possess the universal approximation capability [4], [13]. This class of systems is defined by input-output relationships of the form where u(n) and y(n) represent the input and output signals, re­spectively, and ./'(■••) is a polynomial in N + M + 1 variables. Polynomial systems can be broadly classified into recursive and nonrecursive systems. Nonrecursive polynomial systems are ob­tained from (1) if / depends only on u{n - i), / = 0.1........A'. In this case, (1) takes the form of a truncated Volterra series expansion [36]: y(n) = h0 + ^2 hi{k\)u(n - ki) + E E h2(k1, k2)u(n - ki)u(n - k2) + ■•■+ E E E hp(h,k2,...,kP) x u(n - ki)u{n - k2) ■ ■ ■ u(n - kp) (2) where hm{k\,..., km) represents the mth-order Volterra kernel of the system, and maxjTVi.... ,NP} represents the memory of the system. Conventional identification of a truncated Volterra series aims at estimating the Volterra kernels from either knowledge of the relevant statistics of the input and output signals or measurement of the input and output signal. The mean squared error (MSE) formulation and the least squares error (LSE) formulation en­able the computation of the Volterra kernels via a linear system of equations. Algorithms for the estimation of the parameters of Volterra models based on input-output data have been exten­sively studied in the past [2], [10], [13], [14], [15], [17], [27], [31]-[33], [39], [40]. Most of these methods view the resulting linear regression as a multichannel setup. The Volterra param­eters are then obtained by linear multichannel parameter esti­mation algorithms in batch or in adaptive form. Cumulants and polyspectra are employed in [19]-[22] to estimate symmetric Volterra kernels. These works derive closed-form solutions for the estimates when the input is a stationary, Gaussian, zero mean stochastic process or a linear process. In a similar manner, the identification of Volterra systems of second and third order for general stochastic inputs is treated in [14], [20], [21], and [33]. In general, for blind identification, the output statistics depend nonlinearly on the kernels even when the system is linear. Prob­ably because of the complexity associated with such problems, 1053-587X/03$17.00 © 2003 IEEE KALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 485 little is known about the blind identification of general Volterra systems [9], [24], [35]. Nonrecursive polynomial systems such as the truncated Volterra series expansion encounters serious limitations in practical applications due to the large number of coefficients that need to be estimated. Recursive polynomial models, just like linear IIR filters, can accurately represent many nonlinear systems with greater efficiency than truncated Volterra series representation. A special class of recursive nonlinear models is the class of bilinear systems. The input-output relationship of a bilinear system is given by t K" ' t ' Kb + b(i)u(n - i) where a(i), b(i), and c(i,j) represent the system coefficients, and the set \K„, Kb, Kcy, Kcu } corresponds to the order of the system. Several practical systems have been modeled by bilinear systems [2], [32]. The input-output bilinear representation of (3) is not equivalent to the original state-space bilinear model yk = Cxk + Duk + vk (4) where v and w denote the measurement and process noise. If (3) is transformed into state-space, it involves polynomial non­linearities between state variables. Conventional identification of bilinear state-space models has been studied in [5] and [6] using subspace identification methods. The more general class of state-affine systems, which provide finite dimensional re­alization of Volterra systems with separable kernels, has been treated in [8] using cumulants. Conventional identification methods for input-output bi­linear systems fall into equation error and output error methods. Equation error algorithms are straightforward to develop, and the mean square estimation error surface has a unique minimum. However, this unique minimum is, in general, biased. Output error algorithms are capable of estimating the coefficients without bias. Such enhanced performance is, however, determined by error surfaces that are nonlinear functions of the coefficient values. Consequently, they may contain local minima, and the estimation algorithms may not necessarily converge to the global minimum of their error surfaces. The parameter estimation for both types of methods can be carried out by the LMS algorithm, the extended least squares algorithms, or their variants [11], [32]. A different approach using cross-cumulant information is pursued in [23] and [41]. This approach divides the identification problem into successive solutions of triangular linear systems of equations by considering appropriate slices of the cross-cumulant sequences for each subproblem. Blind identification of bilinear systems has attracted limited attention so far [30], [37]. In these works, closed-form expressions that relate measurable statistics of the output signal to the unknown parameters are derived for a very restricted class of nonlinear system models and for Gaussian inputs. Consequently, the most common approach to estimating the parameters of the model is to resort to some form of numerical search algorithm that operates in an iterative manner [38]. In this paper, we consider the problem of blind identifica­tion of an input-output bilinear system where the matrix of co­efficients of mixed terms is lower triangular. A new algorithm for the identification of bilinear system parameters is presented. The algorithm employs five stages and utilizes output cumulants up to order 4. The derivations are based on the application of the Leonov-Shiryaev theorem [28] to the output cumulants. The rest of the paper is organized as follows. Section II con­tains a formal statement of the blind identification problem. The structure of the solution is described in Section III. The algo­rithm is described in Section IV. The details of the derivation of the identification structure are provided in Appendices A-C. A simulation example that verifies the accuracy of the derivations and demonstrates the quality of the estimates is given in Sec­tion V. Finally, Section VI contains our concluding remarks. II. Problem Statement We consider bilinear systems of the form Kb + b(i)u(n - i) where y(n) is the output of the system, u(n) the input, and v(n) the measurement noise. The input signal u(n) cannot be ac­cessed for measurement. The first term in (5) is characterized by the parameter vector a of size Ka Extending standard terminology, we will refer to this term as the linear AR part of the bilinear model. Similarly, we will call the second term in (5), which is produced by the parameter vector b of size Kb + 1 the linear MA part of the bilinear system. Finally, the third term is called mixed part and is accountable for the nonlinear be­havior of the system. The parameters of the mixed part are de­scribed by a lower triangular matrix C with entries c(i, j) and size KCy x Kcu. The objective of this paper is to estimate the system param­eters a, b, and C using output information only. To make the analysis tractable, we make the following assumptions. 1) The measurement noise v(n) is a zero mean Gaussian white process and is independent of the input signai u(n).486 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51, NO. 2, FEBRUARY 2003 Fig. 1. Block diagram representation of the proposed algorithm. 2) The input signal is a non-Gaussian white process with zero mean value. This means that the cumulants of the input exist and are given by cum[u(n), u{n-h),..., u{n-lk-i)] =7k^(h, • • •, h-1) (8) where <$(Zi,..., h-i) is the (k - 1)-dimensional unit sample signal, and 7& denotes the signal intensity of order k. For technical reasons that will become clear during the derivations, we assume that 74 ^ 373/72. 3) The parameter vectors a and b, as well as the coefficient matrix C, are such that y(n) is a stationary process. Suf­ficient conditions for the stationarity of bilinear processes are derived in [3], [18], and [26]. 4) In order to find a unique solution and overcome the inherent scaling ambiguity of blind identification, we assume that 6(0) = 1. We further assume that c(Kcu, Kcu) 7^ 0 and that Kcu > K&. Finally, we also assume that the system orders Ka, Kcu, and Kcy are known. Based on the above assumptions, a closed-form solution is developed for the estimation of the parameters a, b, and C, using cumulants of y(n) up to order 4. The main components of the method are presented next. III. Organization of the Blind Estimation Algorithm A block diagram representation of the blind estimation algo­rithm is provided in Fig. 1. It is formed by the cascade of several components. The first component is the "HOS estimator." It is fed with ?/(n), which is the measurable output of the system we seek to identify, and estimates cumulants up to order 4. Effi­cient procedures for estimating cumulants both in terms of sta­tistical and computational performance in the time as well as in the frequency domain have been extensively covered in the literature and will not be repeated here [34]. The coefficients of the linear AR part can be directly estimated from the cumulants. This function is performed by the box termed "AR estimator" in Fig. 1. The calculation of the coefficients of the linear MA part, the mixed part, and the statistics of the input requires knowl­edge of the AR coefficients and combinations of the cumulants derived using the estimated AR coefficient values. The combi­nations are generated by two "auxiliary" filters denoted by D and S in Fig. 1. The outputs of the two auxiliary filters and of the AR estimator are fed into the "Mixed-MA" estimator to evaluate remaining parameters. The mixed-MA estimator contains an initialization module and the main module. The initialization module computes the Kcuth (last) column of C and b(Kcu), which is the last entry of the MA part, both scaled by the input variance. This module also calculates additional relationships between input cumulant intensities that are needed by the main module. The main module recursively computes the remaining columns of the mixed part together with the MA part. The functionality of every component is detailed in the next section. IV. Algorithm Description Our method utilizes suitably chosen slices of the output cumulants to estimate the system parameters. We will derive several relationships between these output cumulants and the unknown parameters using a list of properties presented in Appendix A. We first define the following input-output cross- cumulant sequences ^i(^i) =cum[?/(n),w(n-mi)] (9) # 2( ^ 1, ra2) = cum[2/(77), y (n -mi ) , w ( n -m 2)] ( 10) KALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 487 <74(7711, m2,7773,777,4) = cum[?/(n), 7/(77 - 777,1), y(n - n7,2) The following relationships between the output cumulants and the system parameters are derived in Appendix B. Ka c^(k) =cum[y(n),y(n-l1)] = 'y^a(i)c^\li-i) v=0 3=1 + c(hj)92{i - h:j - h): /l>0. (13) ,7=1 *=3 43) (ii, /2) =cum[y (n), y (n - ), y(n - i2)] IU + Yb®92<h-h,i-h) ?=o j=i *=3 Kcu Key + E E C^1-' ‘j)c<y'> ~h) j=l ^ j=l *=3 " ^ ' ' ' ' /1 > 0, l2> 0. 0, y(n-k), V(n-12), y(n-/3)] = ^a(i)c(4)(/1-7, l2-i,h-i) Kb + - - i~k) «=0 + E E c(*'^)42)(l2-i)92(k~k,j~h) j=l *=j + EEc(i'i)cf(|i_®'!3-^i(i-i2) j=i Zi>0, /2>0, Z3>0. (15) The above equations are considerably simplified if the lags are properly restricted. To this end the following proposition is useful. The proof is given in Appendix B. Proposition 1: The following relations hold: 92 { l,7722,r i.niQ.niQ.r )=(), 777,1 <0 (16) ) =0, mi > rrii2, m2<0 (17) ) =0, 7771 , m2 > ms rn3<0 and (18) ) =0, 777,1 > 77?4, 777,2 > *774, ?™3 > ™4 77)4 <0 (19) Next, we describe in detail each component of Fig. 1. A. AR Estimator The AR parameters a(i) are determined with the aid of the following proposition. Proposition 2: Let Li,L2,L-^ > Kcu. Then IQ i=l IQ (14) r(3) 44)( i=i ' ^ '' : (20) and (21) (22) i=i The proof is a direct application of Proposition 1 and (13)-(15). Proposition 2 states that the higher order output cumulant se­quences behave in a manner that is identical to the covariance function of an autoregressive signal for sufficiently large values of the lag I. This property enables the computation of the a(i) parameters via one of the above relations and a linear system Toeplitz solver such as a variant of the Levinson algorithm [13]. The simplest implementation for the AR parameter estimator re­lies on second-order statistics and (20). Collecting R successive output autocovariance lags in the range Kcu < < Kcu + R, we obtain an overdetermined system of linear equations in the unknown parameters. / (2) C, + i)\ C, v 4 2)( 488 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51,NO. 2, FEBRUARY 2003 C „(2)/ -'V \ „,(2) / 0/ (2), t + R-2) ... 42)( + 2--ft'a ) , + R-Ka)) (23) Ka Likewise we define D3(h,l2) = -ya{i)cf\h-i,l2-i) Kcu S3(m,Kcu) = a(l-m)D-i(l,Kcu) / cf2) V V (2)/ y \ ,(2)( -y v ,(2)( -y \ ,(2)( -a v „(2) / + 1 - KCy) ^ + 2 - KCy) c u ~\~ R ~ Ke y) / (30) from which we can easily solve for {a(l), a(2), • • • i general, we obtain more accurate estimation performance by choosing R, > Ka. B. Mixed-MA Estimator The estimation of the rest of the unknown parameters explic­itly utilizes the estimated values of the AR coefficients. Two sets of auxiliary variables are first calculated using the auxiliary filters marked D and S in Fig. 1. The first set of variables are defined as Given measurements of the cumulant values in the above range, we can create an overdetermined set of linear equations in the unknown parameters 72c(Kcu, Kcu). • • •. 72c(Kcy, Kcu) by choosing R, > Kcy - Kcu + 1 equations. Recall that c(i,j) = 0 if i < j, and thus, only the parameters explicitly involved in (30) are nonzero and need to be estimated. Solving for 72c(-, Kr„) is straightforward. b) Estimation of 73/72-' The following relationship is es­tablished in Appendix C: (24) (25) 72 x a (3 )f i - i, L-i)\/D3(Kcu, (31) where a(0) = -1 for all three definitions. In a similar manner, we define three new sets of auxiliary variables as linear combi­nations of D2, D3 and T) \ as follows: where L is such that D3(KCU,L) / 0. The right-hand side of (31) involves quantities available from previous steps. c) Calculation of 72&(-^c«), 74/72 and 75/72-' The rest of the calculations in the initialization module are performed in a similar way. The following equations are derived in Appendix C and can be used directly to estimate the three remaining quantities. cu) - D2 , Kcy -. - , , (32) (27) (28) (33) and 75 1 Kcv L>Kcu, m = Kcu,Kcu-1... (29) We are now ready to estimate the remaining parameters. The mixed-MA estimator contains an initialization module and a main module. These modules are described separately next. 1) Initialization Module: The initialization module esti- a) Computation of ~f2c{-,Kcu): The following system of linear equations is derived in Appendix C: Kcv Kcv (34) 2) Main Module: The main module involves [(Kcu)/(2)] recurrent steps, where [(.)] denotes the largest integer smaller than or equal to (.). During the mth step, we estimate the mth and (Kcu - m)th columns of C, b(m) and b(Kcu - m). The input variance 72 is also estimated in the main module. At theKALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 489 recursions, all the unknown parameters are end of [Kc. estimated. Foreach m = 1,2,..., [(A', „ )/(2)], the computations in this module are organized into three stages. In what follows, we out­line the steps involved in each of the stages. All the derivations are given in Appendix C. Stage 1: The first stage utilizes a linear system of equations to determine the following set of parameters during the mill step: 1) the mth column of C of length Kcy - m + 1, denoted by i))T; 2) the first m entries of the (Kcu - to) th column of C scaled by the input variance. We denote these parameters as 72ckcu_TO = (72c(Kcu-m, Kcu - m),... t 1: Kcu to)) . (35) 3) a linear combination of the remaining terms of the Kcu - to column c'2k m with the last column of C, cKru of length - Kcu + 1 and given by KCy Key x ^2 c(®> fi ^2 j ~Kcu+rfi x |(i -Kcu + m - n, m - n, s + m - n) + c^\Kcu-i) j - Kcu + rn) + c(j -Kcu+ rn, j -Kcu+ rn) - 72 X c(n,j - Kcu + - n. s + m - n). (39) The matrix F has dimensions R x Kcy - m + 1. The (-s 1. .s2) element of F, with 1 < < R and I < .s2 < Kcy - rn + 1, is given by the expression 4) The scalar quantity m ^ )y2 T'2 ( ' ) (37) The computation of the parameters in items 1-4 is performed by solving an overdetermined system of linear equations of the form Key + 72 y^c(n,Kcu) Kcv (F G P Q) | 72>- = E (38) where the vector of the unknown parameters contains 2Kcy - Kcu + 3 elements, and the right-hand-side vector E has length R > 2Kcy - Kcu + 3. The s-th element of this vector is Kr,, ' + ^2 a(i - Kcu + rn)y2 x ^ ^ Kcu)cy ^ (i ti, Kcu s ri) : ^ KCV ; ; " -7* J2 Key x c(nJ ~ Kcu + rn) + 73 ^2c(n,Kcu)c(jf')(n-Kcu + l-s2,s1 - s2 + 1) ^LCy * . . . , Key + Vis "^2 c(n, Kcu)c(y\Kcu + Si - n) + 13^2 C(n-Kcfi ' Kcy + 2/72 "^2 cin^c^c^iKcu-n.Kcu + Si-n). (40) x af (to - n, s + to - n) - j2 x, The matrix G has size R x rn, and its («i, s2)th entry is 81,82 = clf\m + 1 - "S2j - S2 + TO + 1). (41)490 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51, NO. 2, FEBRUARY 2003 The matrix P has T> x Kcy - K, „ + 1 entries given by Psi,s2 = ~ s'2- S1 _ s2 + !)• (42) Finally, the vector Q has R entries given by I<r Kc. Kc <ls r{rn) = b(rn) - + 72c(m, m). 72 74 I \ . 73_ I s - c(m,m)H-----y > c[n,rn) ^2 ^ 72 73 b(j - Kcu + m) H-----c(j - Kcu + m, j - Kcu + m) . 72 x c(n,j - Kcu + m)c^p(i - Kcu + m - n,m - n) Kc ! = 72 Y c(*! Kcu)c(y \Kcu + .Si - i). (43) -73 X] c^! Kcu) Y c(n, m)42) (i - Kcu + m - n) K Kcy + Y a(* - Kcu + rn)72 ^ c(n, Kcu)c^'1 (i - n) Stage 2: The second stage determines b(rri) and 72. Recall from (37) that (44) Suppose rri = rri* is the first integer for which c(m*, rri*) ^ 0. Then, for every m < m*, (44) gives b{m) = r(771)72/73. Thus, all b{m) for 777 < m* are determined. Next, suppose m = rri*. It is shown in Appendix C that b{rn) and 72 also satisfy a linear equation of the form Kcu ( Kc S:s(Kcu - rri, Kcu) - 72 ^ b(J) + yYc^1'^ cy X + m,) c^(m - n) - 72 x |1: /\ : • // £ -i.K, + 72 373c(Kcu, Kcu)c(m, rn) + c(m, 777)72 (45) We note that all quantities in the above equation except b(rri) and 72 are either measurable statistics of the output signal or parameters that have been estimated in previous steps. Conse­quently, we can solve the system of two linear equations (44) and (45) to estimate b(m) and 72. The determinant of the ma­trix associated with the two equations is c(rn,m)c(Kcu,Kcu) ^74 - - J . (46) Since c(m, m) / 0 for m = m*, the assumptions stated in Sec­tion II ensure that the determinant does not vanish. Therefore, b(rn*) and 72 are uniquely determined. The remaining parame­ters 6(777) for m > m* are readily computed from (44). Stage 3: This stage completes the estimation of the (Kcu - to)th column of C. It also estimates b(Kcu - m). Recall that we estimated the product of 72 and the first m entries of the (Kcn - rri) th column of C in step 2 of Stage 1. Since we computed 72 in Stage 2 of the recursion, it is now straightforward to estimate the first rn elements of the (Kcu - m)th column of C. Similarly, we note that we estimated dm = 72c|-cu_m + b(r (47) in step 3 of Stage 1. Since all variables except in have been estimated at this time, we can solve for the entries of the (Kcu - m)th column that were not computed earlier from the above equation. The only other parameter that is estimated in the mill recursion is b(Kcu - rn). This parameter is estimated from the auxiliary sequence S2 and the expression (- m) = S2(Kcu - m) - yy2 ^ c(i, Kcu - m.) KcKALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 491 X c(j,j) - Kcu + m) + (74 + 7!) x c(j - Kcu + m,j - Kcu + m) + 2/73 Key x y c(n,j-Kcu+m) Key x c(n, j - JsTcu + to) -72 53 53 Key x 53 c(n,j - Kcu+rri) X 42)(* - ifcu + to - n). C = / 0.1 0 0 \ 0 -0.05 0 0 0 0.3 \-0.1 0.05 0.1/ TABLE I True and Estimated Parameters for a PRBS Input Sequence of 16383 Samples (100 Monte Carlo Runs) Parameters True Value Mean Variance o(l) -0.1 -0.0976 4.58 1(T7 o(2) 0.02 0.0225 4.41 10"7 6(1) -0.4 -0.3801 6.57 10"4 c(l,l) 0.1 0.0999 1.67 10-4 c( 2,1) 0 -0.0032 2.14 10"5 c(3,l) 0 -0.0015 8.88 10"7 c(4,l) -0.1 -0.1063 4.22 10"5 c( 2,2) -0.05 -0.0267 6.92 10"5 c( 3,2) 0 0.0121 4.61 10"5 c(4,2) 0.05 0.0594 3.97 10"5 c(3,3) 0.3 0.2741 4.08 10"4 c(4,3) 0.1 0.0911 4.47 10"5 72 3.24 3.4946 4.78 10"2 (48) This completes the set of calculations necessary to perform the blind estimation of the bilinear system parameters. V. Simulation Results In this section, we present the results of a simulation experi­ment illustrating the performance of the algorithm. The method is applied to a bilinear system of the form 2 1 y(n) = a(i)y(n - i) + b(i)u(n - i) 34 where Ka = 2, = 1, Kcu = 3, and Kcy = 4, with a = [-0.1 0.02], b = [1 - 0.4], and The input sequence u(n) is a pseudorandom binary sequence (PRBS) generated by a linear feedback shift register. The char­acteristic polynomial of the register is a primitive polynomial. To reduce the realization dependency, the parameter estimates were averaged over 100 Monte Carlo runs. For each experiment, a new PRBS input is generated of length 214 - 1. The mean and the variance of the estimated parameters against the true ones are shown in Table I. VI. Concluding Remarks This paper dealt with the blind identification of bilinear systems with measurements corrupted by Gaussian noise. The excitation is non-Gaussian white noise. The parameters are determined via a sequence of linear systems involving cumulant slices of orders less than four. Simulations validating the proposed method were supplied. One issue that should be pointed out with regard to this work is the need for good experiment design conditions. This, in our case, translates to inputs with good white characteristics in higher order cumulants. Recent work in this direction utilizing dual BCH sequences, Gold sequences, and sequences generated by modulo 2 addition of maximal length sequences of relatively prime periods is reported in [16]. Appendix A Basic Properties of Input Output Cumulants In this appendix, some basic properties of input-output cu­mulants are derived. They are heavily used in the derivation of the blind estimation algorithm. Property 1: Let zi(n- 1), Z2(n - 1),..., Zk(n- 1) be func­tions of u{n - i) and y(n - i) for i > 1. Recall that u{n) and y(n) are the input and output signals, respectively, of the bilinear system and that u(n) is a higher order white sequence with zero mean value. Then cum[u(n),zi(n - l),Z2(n - 1),... ,zk{r = 0. (49) Proof: The conclusion follows immediately from the fact that u(n) is white.492 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51, NO. 2, FEBRUARY 2003 Property 2: Let zi(n - 1). z2(n - 1)...., zk{n - 1) be as defined in Property 1. Then cum[«(n), u(n),..., u(n), z\(n - 1) z2(n - Izk{n - 1)] = 0. (50) Proof: If a random variable X is independent of the random variables Y'i. Y> ■■ ■ ■ ■ V'/,, then [34] (51) Since u(n) is independent of z\(n-l),.... zk(n-1), the result follows. Property 3: For the same set of definitions of the signals as in Property 1 cum['(/(rt), u(n), z\(n - 1), z^ri - 1),..., Zk{n - 1)] = 0. (52) Proof: Let . Ka ' ' ' ' Kb ' ' ' ' (53) so that the bilinear expression becomes y(n) = u{n) + zk+i(n - 1) + v(n). (54) Substituting (54) for y(n) in the cumulant expression gives cum[u(n), u(n), z\(n- 1), .22(n - 1),..., zk{ri - 1)] +cum[i;(?7,), u(n), zi(n - 1), z2(n - 1),.... zk(ti - l)\. (55) All terms on the right-hand side are zero due to Properties 1 and 2 and independence assumptions of v(n) and u(n). This completes the proof. Multilinearity of cumulants [34] leads to the following straightforward generalization of the above result: Property 4: Let y(n) and vi n ) be the output and input of the bilinear system in (5). Then Proof: Suppose first that N = P = 1, and let y(n) = u(n) + zk+i(n - 1) + v(n), as in (54). Substituting for y(n) and employing Property 1, we get cum[y(n) ,u(n)] = cum[«(n), u{n)] + cum^+i (n - 1), u{n)\ Multilinearity of cumulants proves the generalization given as Property 4. Properties 1-4 in combination with the Leonov- Shiryaev theorem [28] form the main tools for the computation of the output statistics. Appendix B Derivation of Output Cumulant expressions We establish (13)-(15) and Proposition 1 in this Appendix. Using (5) and the multilinearity of cumulants, we obtain 42Hk) =cum[y(n),y(n - k)] Ka = ^ a(*)cum[y(n - i), y(n - /1)] Kb + ^ &(*)cum[«(n - i),y(n - Zi)] y c(i,j)cum[y(n - i)u(n - j), y(n - Zi)] + cum['u(ri), y(n - Zi)]. (59) The last term on the right-hand-side of (59) is zero since v(n) and y(n - /|) are independent variables. Application of the Leonov-Shiryaev theorem [28] and the zero mean assumption of the input signal to the third term on the right-hand side of (59) give cum[y(n - i)u(n - j),y(n - Zi)] = cum[t/(n - h) y{n - i),u(n - j)] + ycum[y(n - h), u{n - j)]. (60) Substituting the above result into (59) and making use of the cross-cumulant definitions in Section IV result in (13). A similar approach is used to derive (14). Thus c^\h,l2) =cum[y(n),y(n - h),y(n - l2) (57) where y(n) appears N times, and u(n) appears P times in the above expression. = a(i)cum[y(n - *). y(n - li),y(n - l2)\ Kh + ^6(i)cum[u(n - i),y(n - h),y(n - l2)] + E E c(i, j)cum[y(n - i) x u(n - j),y(n - k), y(n - l2)} + cum[u(n),y(n - h),y(n - l2)]. (61)KALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 493 Applying the Leonov-Shiryaev theorem to the third term on the right-hand side gives cum[r/(n - i)u(n - j), y(n - 1i), y(n - I2)] =cum[j/(n - 11), y(n - i),y(n - l2),u(n - j)] + cum[y(w - i),y(n - l2)}cwa[y(n - l\),u(n - j)\ + ycam[y(n - h),y(n - l2),u(n - j)}. (62) Substituting the above result and the cross-cumulant definitions in (61), along with the use of the independence property of cu- mulants, we obtain (14). Equation (15) is derived in a similar manner. Next, we turn to Proposition 1. First, we present a lemma describing the recursive structure of the cross-cumulant se­quences. The proof of the proposition is a direct consequence of the lemma. Lemma 1: The following recursions hold. Ka ' . ■ E 0 +2/72 E E c(*, mi > 0 , 0 : mi< 0 ■ E K')s( 0 Kcu +74 EC(.j;:j)b(mi-jy>(r m 1 > m,2 > 0 0; m 1 > 0, m2 = 0 s 0; 777,1 > m 2 <() 7-2-j) Ka ' . ' ■ . ■ kcu „(3) + E^ ( ) ( ' * ' xS(m3-i) Kru K +73 E E c(m)42)( (65) and (63) +ra E c(i,m3)#(mi-m2) m2 > m3, m3 > 0 0; mi >0, rn2>0. m3 = 0 „ 0; mi > 777,3, m2 >777,3, 777,3 < 0 E a(*).94(mi-m2-*, m3-i, m,4-i) I<c + E c!('<,m4)[7244) X (mi -i, rri2 - i, mj-i) ■*, 777.3 - i) +73 c(3)^ m3 -* J (66) (64) (2)^ *) ( ) m2 > Wi4, m3 > 777.4, m4 > 0 0; mi>0, m2>0, m.3>0, 0; mi > iri4, m2 > TO4, , m3 >m4, m4 < 0 . Proof: We start by proving (63). Property 1 ensures that .91(777-1) = 0 for 777,1 < 0. The initialization stage .91 (0) = 72 follows directly from Property 4. Substituting the bilinear equation (5) into the definition of c/i gives .91(^1) =cum[r/(?7,), 7/(71 - m-i)] = 53 a(*)CUm[?Xn - ®): W'(n - ml)494 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51, NO. 2, FEBRUARY 2003 Kb + &(*)cum[?x(r?, - *). u(n - mi)] c(i, j)cum[y(n - i)u(n - j),u(n - mi)] + cum[v(n),u(n - mi)] ; mi > 0. (67) Since v(n) and u(n - rrii) are independent, the last term on the right-hand side of the above expression is zero. Using the Leonov-Shiryaev theorem on the third term of the right-hand side yields cum[j/(n - i)u{n - j),u(n - mi)] = cum[j/(n - i) (68) Among the terms of the form cum[j/(n-*), u{n-j), u(n-m!)], with 1 < j < Kcu and j < i < Kcy, only the term cum [//( ;/ - mi), u(n-m,i), u(n - mi)] is nonzero due to Property 1. When i = j = rrii, Property 4 implies that cum[(/(ri - mi),u(n - mi),u(n - mi)] = 73. Applying these results along with the fact that 11( 71) is white leads to (63). Next, we prove (64). Properties 1 and 3 ensure that .92(^1, m2) = 0 for m,i > m2 < 0. The initialization .92(0,0) = 73 follows directly from Property 4. Substituting the bilinear equation (5) into the definition of g2 gives 9 2( I =cum[y(n),y(n - rni),u(n - 'i Kb 73 Yj ~ i)S(m2 - i) kcu kc c.(i, j)cvm[y(n - i), u ( n-j ) , y ( n - mi ) , u(n - m 2)] Kcu Kcy +EEc('«') x cum[y(n - i),y(n - mi)]cum[w(n - j),u(n - m2)] +ziEc(*^') x cum[y(n - i), u(n - rn2)]cum[u(n - j), y(n - mi)] + V EE c(i,j)cum[u(n-j), y(n - m 1), u(n - m2)]. (71) Arguing as before, we find that cum[t/(n - i),u(n - j),y(n - mi),u(n - m2)] = 0 for j ^ m2 and any i. Moreover, for j = m2, the only nonzero term occurs when i = j = mi = m2. Hence, the first term becomes I<CU • (72) In a similar manner, we can show that the second term is 72 (73) a(i)cvm[y(n - i),y(n - rrii).u(n - rri2)] ^=1 Kb + 6(i)cum[u(n -*), y(n -mi), u(n - m2)] + 'Yj c(i,j)cum[y(n-i)u(n-j), 3=1 y(n - mi). u(n - m,2)] mi>m2> 0. (69) The last term on the right-hand side of the above equation is zero. Let us focus now on the second term in the right-hand side of the same equation. For any m2 > 0, mi > m2 and i < m2, cum[«(n - -i), y(n - mi),u(n - m2)\ = 0 because of Property 2. If i > m2, the above cumulant is again zero due to Property 1. Hence, the second term becomes (70) and is nonzero only when 0 < m2 < Kb and mi = m2 = i. Using the Leonov-Shiryaev theorem, the third term can be expanded as and that the third term is Recall that for the calculations of interest here, i > j, and m 1 > m2. Since .91 (mi) = 0 for m, 1 < 0, the only nonzero term results when i = j = rnt = m2. Hence, the previous term becomes 72c(m2,m2)£(mi - m2). (75) Reasoning as before, the last term of (71) becomes yi'i Y c(hj)S(rrii - j)S(m2 - j). (76) Substituting the above results in (69) and making use of the cross-cumulant definitions, we obtain (64). Similar arguments are employed to establish (65) and (66). APPENDIX C MIXED-MA ESTIMATOR Initialization Module We begin the derivations by considering the relationship for cf\hM) for h = Kcu and l2 = L > Kcu. Using the simpli­fications possible through Lemma 1, we get r‘(3) I KcvKALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 495 We observe from the definition of T)-> (/ |. U) in (25) that (78) ^ Kc y + 272 5] c(i,Kcu)c^(Kcu-i). (84) Substituting this result in (77) gives K cv Finally, (34) follows from (15) with h = l2 = h = Kcu. Then Successive evaluation of (78) for R > Kcy - Kcu + 1 values of L leads to (30). To derive (31), we evaluate (15) for li = Kcu, l> = Kcu, and h = L> Kcu. Applying Lemma 1 to the various terms in (15) for these choices of the parameters, we get I<a Kcy . ■ . . +373 jr c(i,Kcu)c^(Kcu-i). (85) Kc + 7.3 53 c(i,Kcu)c£\L - *). (80) Substituting the definitions for Di(Kcu, Kcu, L) from (26) and the expression for D3(KCU, L) from (78) in (80) results in (81) This completes the derivation of the initialization module. Main Module We first substitute (24) in (13) and apply Lemma 1 to get Kcv, ( Kcy £>2(o = 5>c? - o U(j) + //53r,;'--/: We then multiply both sides of (86) by a,(I - to) and add the result over I in the range m < Z < Kcu. This operation gives The expression for 73/72 in (31) follows from this result. Next, we consider (13) for l\ = Kcu. Again, using Lemma 1, we get . ■ . Ka ' . . ' . , Kcu Kcu = - 53 - m) 53gid ~\b(j)+v 53 ■?)) 0 < m < Kcu. (87) Substituting for D> (K, „) from (24) in the above equation and rearranging the terms gives If we change the order of summations in the first term of the right-hand side, we have Kr Kcu A„ Kcv * ^ (83) which is identical to (32). In a similar manner, (33) follows from (14) for the case when h = h = Kcu. Then l=r Kcu I<CV J Key j-m \ i-j / V l-m / Kcu / Key \ / j-m \ =53 w+y-£c(M) _53a(n)si(j-m~n) n=0 (88) 496 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51, NO. 2, FEBRUARY 2003 Substituting (63) for // i on the right-hand side of the above equa­tion transforms this term to i=j " : ^2 c(l,j - rn) (89) Likewise, by changing the order of summations in the second term of the right-hand side of (87), we find that b=m j=l v=j j=m ^ fc=m ' ' ^ ^2 X! c(®> fi a('fiv2(®_m_rt> J _m_n) j • (90) Applying (64) to the above equation, its right-hand side becomes ' ' ' Kcv ' +yi3 ^2 c(l'j-rn) l-j - m +72 Y2 ci<^fi Yl c{l,j (91) If we substitute (89) and (91) into (87), we get Kc ! J2 C(1^-'1 h=j-m Kc i - rn) h=j-m KCU Kcy Kcy +12^2^2 c{i,fi j=m r=j+1 k=j^m Using a similar approach, we next derive expressions for S3(rn, Kcu) and S^rri, Kcu, L). We only outline the derivation of S3(m,Kcu) below. Application of Lemma 1 to (14) for l2 = Kcu gives Ka 43) (l.Kcu)=J2 a(*)43) (I - i, Kcu-i) Kh ' , ' iz ' (} ( / 1 i=j (93) Substituting for Dll. Kcu) from (25) into the above equation results in A'c„ ^ ^ /' . (94) Multiplying both sides of (94) with a(l - rn) and adding the results over I in the range rn < I < K, „ gives S3(m,Kcu) = K b=tn Kcu Kcu Kc 0 < rn < K cu. (92) h=m i-j Kcu Kcu K^y . . ^ . . h=m j=l i=j KALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 497 I<c Kc 0 < m < Kcu. (95) Kc Applying (64) to the above expression results in Kc„ / Kcv \ Kc l=Kc K cv +2/74 EZ c(l-,Kcu-m) l=Kc X Substitution of (97)-(99) into (95) gives S3(m,,Kcu) = As was done before, we change the order of the summations in the first term of the right-hand side for (95), and we obtain Kcu I<cy i=i , f j-m - EZ a(n)92(Kcu - m-n,j-m-n) . (96) 72 Z y'i ' l-i J Ylc^l^~m^c<y)^Kcu~m~1^ , Kcy , i=Kcu 722^ I tAD+y2_!',;'••/•' j=m y fcg EZ ' ' ' ' '( } ' ' Kc. ' i=Kcu ' +1/73 EZ c(l'Kcu-rn) Kc (97) In a similar manner, we can show that the second term of the right-hand side of (95) is equal to 72 EE E c(i,j)c(l,j-m)c(jf'>(i-m-LKcu-m-l) j=m i=j h=j-m j=m m +73 E] E <ii,Kcu)c{LKcu-m)cf]{i-'m-l) i=Kcu l=Kcu-rn l=Kcu-m j=m i=j Kcu K c y K c y + 72EZEZ EZ c(i-j)c(^j~m)cyi'>(i~rn~^Kcu-rn-l) j=m i=j l=j-rn ' Kcu Key ' j=m fc=/-m +73 EZ EZ c(i;Kcu)c(l,Kcu-m)c^\i-m-l) i^=Kcu h=Kcu-m +2/74 A' b=Kcu-m I< Kcy fc=ra i=K (100) (98) Similarly, the third term of the right-hand-side of (95) becomes j=m j=j (99) A similar analysis leads to the following expression for AVu / 72 EZ &(i)+^EZc(*'^ 7'=m+l \ h=j-m Key ‘ ^ i=Kcu ^498 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 51,NO. 2, FEBRUARY 2003 Kc Kcu K c y K c y +72 J2J2Y1 c^(i-m,-l, Kcu - m - l, L - m - l) ^ K cu K cy ^ ^ ^ ' j=j7?+l k=j-m +73 E E c(i, Kcu)c(l, Kcu - m) i=Kcu h=Kcu-m cS^(i - rn - l, L- m - l) l=Kcu-rn Kcu Kcy Kcy +72 E Y1 Y1 c{i,j)c{l,j-rn) j=ra-|-l i=j fc=j-m +-«££E j=m. i=j I-j-m Kcu 0 + E c^-Kc Key ' -H773 E c(hKcu-m) i~m) b=Kcu-m n i=j K Kcy -ya{l-m)l2 E c(i,Kcu)cf\l-i,L-i) k=m, i=Kcu 0 < m < Kcu. (101) The substitution rn -- K, „ -m in (101), (100), and (92) leads to (38), (45), and (48). References [1] S. Benedetto and E. Biglieri, "Nonlinear equalization of digital satellite channels," IEEE J. Select. Areas Commun., vol. SAC-1, pp. 57-62, Jan. 1983. [2] S. A. Billings, "Identification of nonlinear systems-a survey," Proc. Inst. Elect. Eng. D, vol. 127, no. 6, pp. 272-285, Nov. 1980. [3] T. Bose and M. Q. Chen, "BIBO stability of the discrete bilinear sys­tems," Digital Signal Process. Rev. J., vol. 5, no. 3, 1995. [4] S. Boyd, "Volterra Series: Engineering Fundamentals," Ph.D. disserta­tion, Univ. California, Berkeley, 1985. [5] W. Favoreel, B. De Moor, and P. Van Overschee, "Subspace identifica­tion of bilinear systems subject to white inputs," IEEE Trans. Automat. Contr., vol. 44, pp. 1157-1165, June 1999. [6 [7 [8 [9 [10 [11 [12 [13 [14 [15 [16 [17 [18 [19 [20 [21 [22 [23 [24 [25 [26 [27 [28 [29 [30 [31 [32 W. Favoreel, "Subspace methods for identification and control of linear and bilinear systems," Ph.D. dissertation, Dept. Elect. Eng., Katholieke Univ. Leuven, Leuven, Belgium, 1999. W. A. Frank, "An efficient approximation to the quadratic Volterra filter and its application in real-time loudspeaker linearization," Signal Process., vol. 45, no. 1, pp. 97-114, July 1995. G. Gatt and N. Kalouptsidis, "Identification of discrete-time state affine state space models using cumulants," Automatica, vol. 38, pp. 1663-1681, Oct. 2002. G. B. Giannakis and E. Serpedin, "Linear multichannel blind equalizers of nonlinear FIR Volterra channels," IEEE Trans. Signal Processing, vol. 45, pp. 67-81, Jan. 1997. G. O. Glentis, P. Koukoulas, and N. Kalouptsidis, "Efficient algorithms for Volterra system identification," IEEE Trans. Signal Processing, vol. 47, pp. 3042-3057, Nov. 1999. S. Han, J. Kim, and K. Sung, "Extended generalized total least squares method for the identification of bilinear systems," IEEE Trans. Signal Processing, vol. 44, pp. 1015-1018, Apr. 1996. R. Hermann, "Volterra modeling of digital magnetic saturation recording channels," IEEE Trans. Magn., vol. 26, no. 5, pp. 2125-2127, Sept. 1990. N. Kalouptsidis, Signal Processing Systems: Theory and Design. New York: Wiley, 1997. K. I. Kim and E. J. Powers, "A digital method of modeling quadratically nonlinear systems with a general random input," IEEE Trans. Acoust., Speech, Signal Processing, vol. 36, pp. 1758-1769, Nov. 1988. T. Koh and E. J. Powers, "Second-order Volterra filtering and its appli­cation to nonlinear system identification," IEEE Trans. Acoust. Speech, Signal Processing, vol. ASSP-33, pp. 1445-1455, Dec. 1985. N. Kolokotronis, G. Gatt, and N. Kalouptsidis, "On the generation of sequences simulating higher order white noise for a system identifica­tion," IEEE Trans. Signal Processing, submitted for publication. M. J. Korenberg and L. D. Paarmann, "Orthogonal approaches to time- series analysis and system identification," IEEE Signal Processing Mag., vol. 8, pp. 29-43, July 1991. S. Kotsios and N. Kalouptsidis, "BIBO stability criteria for a certain class of discrete nonlinear systems," Int. J. Contr., vol. 58, no. 3, pp. 707-730, Sept. 1993. P. Koukoulas and N. Kalouptsidis, "Nonlinear system identification using Gaussian inputs," IEEE Trans. Signal Processing, vol. 43, pp. 1831-1841, Aug. 1995. -----, "Third order Volterra system identification," in Proc. IEEE Int. Conf. Acoust. Speech, Signal Process., Munich, Germany, Apr. 1997, pp. 2405-2408. -----, "Second-order Volterra system identification," IEEE Trans. Signal Processing, vol. 48, pp. 3574-3577, Dec. 2000. P. Koukoulas, V. Tsoulkas, G. Gatt, and N. Kalouptsidis, "Identifica­tion of discrete-time Volterra series using higher order statistics," in Preprints first Nonlinear Contr. Network Pedagogical School, Athens, Greece, Sept. 1999, pp. 149-201. P. Koukoulas, V. Tsoulkas, and N. Kalouptsidis, "A cumulant based al­gorithm for the identification of input output quadratic systems," Auto- matica, vol. 38, pp. 391-407, Mar. 2002. P. Koukoulas and N. Kalouptsidis, "Blind identification of second order Hammerstein series," Signal Process., vol. 83, pp. 213-234, Jan. 2003. G. Lazzarin, S. Pupolin, and A. Sarti, "Nonlinearity compensation in digital radio systems," IEEE Trans. Commun., vol. 42, pp. 988-999, Feb./Apr. 1994. J. Lee and V. J. Mathews, "A stability condition for certain bilinear sys­tems," IEEE Trans. Signal Processing, vol. 42, pp. 1871-1873, July 1994. -----, "A fast recursive least squares adaptive second-order Volterra filter and its performance analysis," IEEE Trans. Signal Processing, vol. 41, pp. 1087-1102, Mar. 1993. V. P. Leonov and A. N. Shiryaev, "On a method of calculation of semi­invariants," Theory Prob. Appl., vol. 4, pp. 319-329, 1959. P. Z. Marmarelis and V. Z. Marmarelis, Analysis of Physiological Sys­tems- The White Noise Approach. New York: Plenum, 1978. V. J. Mathews and T. K. Moon, "Parameter estimation for a bilinear time series model," in Proc. IEEE Int. Conf. Acoust. Speech Signal Process., Toronto, ON, Canada, May 1991, pp. 3513-3516. V. J. Mathews, "Adaptive polynomial filters," IEEE Signal Processing Mag., vol. 8, pp. 10-26, July 1991. V. J. Mathews and G. L. Sicuranza, Polynomial Signal Pro­cessing. New York: Wiley, 2000.KALOUPTSIDIS et al.: BLIND IDENTIFICATION OF BILINEAR SYSTEMS 499 [33] S. W. Nam and E. J. Powers, "Application of higher order spectral anal­ysis to cubically nonlinear system identification," IEEE Trans. Signal Processing, vol. 42, pp. 1746-1765, July 1994. [34] C. L. Nikias and A. P. Petropulu, Higher-Order Spectral Analysis: A Nonlinear Signal Processing Framework. Englewood Cliffs, NJ: Pren- tice-Hall, 1993. [35] S. Prakriya and D. Hatzinakos, "Blind identification of LTI-ZMNL-LTI nonlinear channel models," IEEE Trans. Signal Processing, vol. 43, pp. 3007-3013, Dec. 1995. [36] M. Schetzen, The Volterra and Wiener Theories of Nonlinear Sys­tems. New York: Wiley, 1980. [37] S. A. Sesay and T. Subba Rao, "Difference equations for higher order moments and cumulants for the bilinear time series model BP(p,0,p,1)," J. Time Series Anal., vol. 12, no. 2, pp. 159-177, Feb. 1991. [38] T. Subba Rao, "On the theory of bilinear time series models," J. Roy. Soc. Ser. B, vol. 43, pp. 244-255, 1981. [39] M. A. Syed and V. J. Mathews, "Lattice algorithms for recursive least squares adaptive second order Volterra filtering," IEEE Trans. Circuits Syst. II, vol. 41, pp. 202-214, Mar. 1994. [40] , "QR-decomposition based algorithms for adaptive Volterra fil­tering," IEEE Trans. Circuits Syst. I, vol. 40, pp. 372-382, June 1993. [41] V. Tsoulkas, P. Koukoulas, and N. Kalouptsidis, "Identification of input output bilinear systems using cumulants," IEEE Trans. Signal Processing, vol. 49, pp. 2753-2761, Nov. 2001. Nicholas Kalouptsidis (SM"85) was born in Athens, Greece, on September 13, 1951. He received the B.Sc. degree in mathematics (with highest honors) from the University of Athens in 1973 and the M.S. and Ph.D. degrees in systems science and mathematics from Washington University, St. Louis, MO, in 1975 and 1976, respectively. He has held visiting positions at Washington Uni­versity; Politecnico di Torino, Torino, Italy; North­eastern University, Boston, MA; and CNET Lannion, Paris, France. He has been an Associate Professor and Professor with the Depart­ment of Physics, University of Athens. In the Fall of 1998, he was a Clyde Chair Professor with the School of Engineering, University of Utah, Salt Lake City. He is currently a Professor with the Department of Informatics and Telecommuni­cations, University of Athens. He is the author of the textbook Signal Processing Systems: Theory and Design (New York: Wiley 1997) and coeditor, with S. Theodoridis, of the book Adaptive System Identification and Signal Processing Algorithms (Englewood Cliffs, NJ: Prentice-Hall, 1993). His research interests are in system theory and signal processing. Panos Koukoulas (M"02) was born in Athens, Greece, on November 17, 1965. He received the B.Sc. degree in physics in 1989, the M.S. degree in electronics and communications in 1993, and the Ph.D. degree in informatics in 1997, all from the University of Athens. From 1997 to 2000, he was with the Access Network & Wireless Communications Department, Intracom S.A., Athens, where he worked on aspects of digital modem design. He is currently with the Radar and Automation Department of the Hellenic Civil Aviation Authority, Athens. His research interests include nonlinear system identification, higher order statistics, and digital signal processing. V. John Mathews (F'02) received the B. E. (Hons.) degree in electronics and communication engineering from the University of Madras, Madras, India in 1980 and the M.S. and Ph.D. degrees in electrical and com­puter engineering from the University of Iowa, Iowa City, in 1981 and 1984, respectively. At the University of Iowa, he was a Teaching/Re­search Fellow from 1980 to 1984 and a Visiting As­sistant Professor with the Department of Electrical and Computer Engineering from 1984 to 1985. He joined at the University of Utah, Salt Lake City, in 1985, where he is now Professor and Chairman of the Department of Electrical and Computer Engineering. His research interests are in adaptive filtering, non­linear filtering, image compression, and application of signal processing tech­niques in communication systems and biomedical engineering. He is the author of the book Polynomial Signal Processing (New York: Wiley), which he co-au­thored with Prof. G. L. Sicuranza, University of Trieste, Trieste, Italy. Prof. Mathews has served as a member of the Signal Processing Theory and Methods Technical Committee, the Education Committee, and the Conference Board of the IEEE Signal Processing Society. He is currently the Vice Presi­dent- for Finance of the IEEE Signal Processing Society He is a past asso­ciate editor of the IEEE Transactions on Signal Processing and the IEEE Signal Processing Letters. He was the General Chairman of IEEE Inter­national Conference on Acoustics, Speech, and Signal Processing (ICASSP) in 2001.