PARAMETER ESTIMATION OF MOVING AVERAGE PROCESSES

USING CUMULANTS AND NONLINEAR OPTIMIZATION

ALGORITHMS

M. Boulouird, M. M. Hassani

Laboratoire d’Electronique et Instrumentation LEI (FSSM/UCAM)

Department of Physics, Faculty of Sciences Semlalia

Avenue Prince My Abdellah, B.P 2390, 40001 Marrakesh, Morocco

G. Favier

Laboratoire d’Informatique, Signaux et Systèmes de Sophia-Antipolis I3S (CNRS/UNSA)

Les algorithmes - Bâtiment Euclide B, 2000 route des Lucioles

B.P 121, Sophia Antipolis Cedex France

Keywords:

MA system identiﬁcation, Higher-Order Statistics, Estimation parameter, Linear algebra solution, Gradient

descent algorithm, Gauss-Newton algorithm, Cumulants.

Abstract:

In this paper nonlinear optimization algorithms, namely the Gradient descent and the Gauss-Newton algo-

rithms, are proposed for blind identiﬁcation of MA models. A relationship between third and fourth order

cumulants of the noisy system output and the MA parameters is exploited to build a set of nonlinear equa-

tions that is solved by means of the two nonlinear optimization algorithms above cited. Simulation results are

presented to compare the performance of the proposed algorithms.

1 INTRODUCTION

Numerous methods have been proposed in the lite-

rature for blind identiﬁcation of MA models using

cumulants. The present paper is concerned with

the linear algebra solutions approach. It consists in

constructing a system of equations obtained from ex-

plicit relations that link third and fourth order cumu-

lants of the noisy output with the MA parameters and

solving this system by the least-squares method ((Al-

shebeili, 1993), (Giannakis, 1989), (Martin, 1996),

(Na, 1995), (Srinivas, 1995), (Stogioglou, 1996), (Tu-

gnait, 1990), (Tugnait, 1991)). In order to take the

redundancy in the unknown parameters vector into

account, (Abderrahim, 2001) proposed a constrained

optimization based solution.

In this paper, we propose another approach to re-

duce this redundancy. It consists in exploiting the non-

linearity existing in the unknown parameters estima-

ted vector. In the literature, the parameters of the vec-

tor to be estimated are regarded as independent, but

actually it isn’t the case. Thus the major contribution

of this paper lies in the estimates of a non redundant

vector of unknown parameters.

The organization of this paper is as follows. The

problem formulation is given in Section 2. In Sec-

tion 3, the resolution with least squares and the non-

linear optimization algorithms, in the event Gradient

descent and Gauss-Newton algorithms, well be deve-

loped. Computer simulation results are given in Sec-

tion 4 to show the effectiveness of the proposed tech-

niques. Finally, the paper is concluded in Section 5.

2 PROBLEM FORMULATION

We consider the discrete, causal, linear time-invariant

process represented on ﬁgure 1, with the following

assumptions :

H.1. The input w(k) is a zero mean, independent

and identically distributed (i.i.d), stationary non-

Gaussian, non measurable real sequence, with

unknown distribution, and :

m,w

(τ

, τ

, . . . , τ

m−1

) = γ

m,w

δ(τ

, τ

, . . . , τ

m−1

)

where :

⋄ C

m,w

(τ

, τ

, . . . , τ

m−1

) is the mth-order cu-

mulant of the input signal of the MA model.

⋄ γ

m,w

6= 0, ∀ m ≥ 2

⋄ γ

m,w

= C

m,w

( 0, 0, . . . , 0

m−1

)

Boulouird M., M. Hassani M. and Favier G. (2005).

PARAMETER ESTIMATION OF MOVING AVERAGE PROCESSES USING CUMULANTS AND NONLINEAR OPTIMIZATION ALGORITHMS.

In Proceedings of the Second International Conference on Informatics in Control, Automation and Robotics - Signal Processing, Systems Modeling and

Control, pages 11-15

DOI: 10.5220/0001182900110015

 SciTePress

⋄ γ

2,w

= σ

= E

w(k)

is the variance of

w(k).

⋄ γ

3,w

= E

w(k)

is the skewness of w (k).

⋄ γ

4,w

= E w(k)

− 3 E w(k)

is the

kurtosis of w(k).

H.2. The additive noise v(k) is assumed to be an

i.i.d Gaussian sequence with unknown variance,

zero-mean, and independent of w (k).

H.3. The non measurable output x(k) is assumed

to be a nonminimum or minimum phase MA

process.

H.4. The order q of the model is assumed to be

known.

w(k)

H(z)

+

x(k)

v(k)

y(k)

Figure 1: Single-channel system.

The measured noisy MA process y(k) is represen-

ted by the following equations :

x(k) =

i=0

h(i) w(k − i); {h(0) = 1} (1)

y(k) = x(k) + v(k) (2)

For the MA model described by equation (1) with

the assumptions H.1, H.2, H.3, and H.4, the mth and

nth-order cumulants of the MA system output (2) are

linked by the following relation (Abderrahim, 2001) :

max

i=i

min

h(i)

m−s−1

k=1

h(i + τ

)

n,y

(β

, β

, . . . , β

n−s−1

, i+α

, . . . , i+α

) =

n,w

m,w

max

j=j

min

h(j)

n−s−1

k=1

h(j + β

)

m,y

(τ

, τ

, . . . , τ

m−s−1

, j+α

, . . . , j+α

)

(3)

where m > 2, n > 2 and s is an arbitrary integer

number satisfying : 1 ≤ s ≤ min(m, n) − 2,

and











min

= max(0, −τ

, · · · , −τ

m−s−1

)

max

= min(q, q − τ

, · · · , q − τ

m−s−1

)

min

= max(0, −β

, · · · , −β

n−s−1

)

max

= min(q, q − β

, · · · , q − β

n−s−1

)

Setting n = 3, m = 4, and s = 1 in equation (3),

yields

max

i=i

min

h(i)h(i + τ

)h(i + τ

3,y

(β

, i + α

) =

3,w

4,w

max

j=j

min

h(j)h(j + β

4,y

(τ

, τ

, j + α

) (4)

where











min

= max(0, −τ

, −τ

)

max

= min(q, q − τ

, q − τ

)

min

= max(0, −β

)

max

= min(q, q − β

)

By setting τ

= τ

= 0 in (4), we get the rela-

tion used in this paper for estimating the parameters

{h(i)}

i=1,2,...,q

of the MA model.

i=0

(i)C

3,y

(β

, i + α

) =

3,w

4,w

max

j=j

min

h(j)h(j + β

4,y

(0, 0, j + α

) (5)

It is important to determine the range of

values of α

and β

so that the cumulants

3,y

(β

, i + α

)}

i=0,··· ,q

, {C

4,y

(0, 0, j +

)}

j=j

min

,··· ,j

max

, and the coefﬁcients {h(j + β

)}

be not all zero for each equation.

By taking account of the property of causality of

the model and the domain of support for third and

fourth order cumulants of an MA(q) process (Mendel,

1991), we obtain :











−q ≤ β

≤ q

−2q ≤ α

≤ q

−2q + β

≤ α

≤ q + β

(6)

Using the symmetry properties of cumulants (Ni-

kias, 1993), the set of values for α

and β

is deﬁned

by :

(

−q ≤ β

≤ 0

−2q ≤ α

≤ q + β

(7)

3 PARAMETER ESTIMATION

3.1 Least-Squares (LS) Solution

Concatenating (5) for all values of α

and β

deﬁned

by (7), we obtain the following system of equations :

Mθ = r (8)

ICINCO 2005 - SIGNAL PROCESSING, SYSTEMS MODELING AND CONTROL

where :

θ = [h(1) · · · h(q) h

(1) h(1)h(2)

· · · h(1)h(q) h

(2) · · · h(2)h(q) h

(3)

· · · h

(q) ǫ

4,3

(1) · · · ǫ

4,3

(q)]

(9)

⋄ ǫ

4,3

= γ

4,w

/γ

3,w

⋄ M is a matrix of dimension

+7q+2

+5q+2

⋄ θ is a vector of dimension

+5q+2

, 1

⋄ r is a vector of dimension

+7q+2

, 1

Assuming that (M

M) is invertible, the unique LS

estimate of θ is

θ = (M

−1

r (10)

Solving (8) provides the parameter estimates

{h(i)}

i=1,··· ,q

as the ﬁrst q components of the

estimated parameter vector

θ in (10).

3.2 Gradient Descent Algorithm

(GDA)

The idea satisfying this paper is to reduce the dimen-

sion of the estimated parameter vector (9). In section

3.1, θ is a vector of

+5q+2

elements. The linear

algebra solutions regard the elements of the vector θ

clariﬁed in relation (9) as independent parameters, but

the dependence of these elements is almost obvious.

To palliate the problem of redundancy in θ, we pro-

pose a new approach based on non-linear optimiza-

tion algorithm. In this part, the parameters vector θ is

a (q + 1) length vector. It has this form :

= [h(1), · · · , h(q), ǫ

4,3

]

(11)

The criterion to be minimized in this case is as fol-

lows :

= kr − φ(θ

The GDA solution has the following form :

i+1

+ λJ

(r − φ(

)) (12)

where :

⋄ r is deﬁned in section 3.1.

⋄ φ is the system of equations obtained by conca-

tenating (5) for all values of α

and β

deﬁned

by (7).

⋄ J is the Jacobian matrix of φ,

J =



∂φ

∂θ



(k,l)

where k = 1, · · · ,

+7q+2

, and l = 1, · · · , q+1.

⋄ λ is the step-size.

The parameter ǫ

4,3

must be estimated since we are

supposed that we don’t know the nature of the distri-

bution of the input signal w(k).

3.3 Gauss-Newton Algorithm (GNA)

This algorithm has this form :

i+1

+ µ(J

−1

(r − φ(

))

(13)

where :

⋄ r, φ, and J are deﬁned in section 3.2.

⋄

has the form of (11).

⋄ µ is the step-size of this algorithm.

4 SIMULATION RESULTS

To demonstrate the effectiveness of the proposed

techniques, let us examine two examples treated in

the literature. In both models the input signal w(k) is

a zero-mean exponentially distributed i.i.d noise se-

quence with γ

2,w

= σ

= 1 and γ

3,w

= 2. We deﬁne

the Signal-to- Noise Ratio as

SN R(dB) = 10 log



E[x

(k)]/E[v

(k)]



For each run, we calculate the Normalized Mean

Square Error (NMSE) deﬁned as

NMSE =

i=1



h(i) −

h(i)



i=1

(i)

where h(i) and

h(i) are respectively the actual and the

estimated impulse responses, respectively. The Error

to Signal Ratio (ESR) in decibels is also used as a

measure of the estimation error. The ESR is deﬁned

ESR(dB) = 10 log

(NMSE)

Example 1 :

y(k) = w(k)−2.3333w(k−1)+0.6667w(k−2)+v(k)

The zeros of the system transfer function H(z) are lo-

cated at 2 and 0.3333. This model has also been used

in (Abderrahim, 2001), (Giannakis, 1989), and (Srini-

vas, 1995).

Additive colored noise is generated as the output of

the following MA(2) model (Abderrahim, 2001) :

v(k) = e(k) + 2.3333e(k − 1) − 0.6667e(k − 2)

where the input sequence e(k) is an i.i.d Gaussian

sequence. We carried out Monte Carlo simulations

with K = 100 different noise sequences, N = 5120

data for each run, and three different values of SNR

(20dB, 10dB, and 0dB). The simulation results are

summarized in Tables 1, 2, and 3.

Example 2 :

y(k) = w(k) + 0.1w (k − 1) − 1.87w(k − 2)

+3.02w(k−3)−1.435w(k−4)+0.49w(k−5)+v(k)

PARAMETER ESTIMATION OF MOVING AVERAGE PROCESSES USING CUMULANTS AND NONLINEAR

OPTIMIZATION ALGORITHMS

Table 1: System identiﬁcation using signal model MA(2),

Mean, Standard deviation, NMSE, and ESR . SNR=20dB.

True Values LS GDA GNA

λ = 8. 10

−5

µ = 0.5

h(1) -2.3333 mean -2.1704 -2.2676 -2.2676

σ 0.0585 0.0149 0.0149

h(2) 0.6667 mean 0.5101 0.6541 0.6541

σ 0.0211 0.0060 0.0060

NMSE mean 0.0093 0.0008 0.0008

σ 0.0043 0.0003 0.0003

ESR -20.3041 -30.9500 -30.9500

Table 2: System identiﬁcation using signal model MA(2),

Mean, Standard Deviation, NMSE, and ESR. SNR=10dB.

True Values LS GDA GNA

λ = 8. 10

−5

µ = 0.5

h(1) -2.3333 mean -2.1013 -2.1158 -2.2624

σ 0.1835 0.4664 0.0511

h(2) 0.6667 mean 0.4987 0.6166 0.6557

σ 0.0632 0.1275 0.0192

NMSE mean 0.0203 0.0478 0.0014

σ 0.0181 0.1490 0.0015

ESR -16.9250 -13.2057 -28.5387

The zeros of the system transfer function H(z) are lo-

cated at −2, 0.7±j0.7 and 0.25±j0.433. This model

has also been used in (Alshebeili, 1993), (Stogioglou,

1996), (Tugnait, 1990), and (Tugnait, 1991).

Additive colored noise is generated as the output of

the following MA(3) model (Abderrahim, 2001) and

(Na, 1995) :

v(k) = e(k)+0.5e(k−1)−0.25e(k−2)+0.5e(k−3)

where the input sequence e(k) is an i.i.d Gaussian se-

quence. In this case N = 10240. The simulation re-

sults are given in tables 4, 5, and 6.

In these simulations, we initialize the nonlinear op-

timization algorithms with the LS solution. The ad-

vantage of this is to avoid the convergence to local

minimum. The good value of the step-size allows also

to avoid this problem of the local minima.

Table 1 shows that the parameters estimation via

the proposed method is much more powerful than LS

solution in term of mean value, standard deviation, or

NMSE, and ESR . The Gradient Descent and Gauss-

Newton algorithms converge to the same values, but

the convergence speed is different. The ﬁgures 2 and

3 illustrate this.

Tables 2 and 5 show that the Gradient Descent al-

gorithm is sensitive to the additive noise. However,

the Gauss-Newton algorithm is more robust to mea-

surement noise.

We can note that for a low SNR (SNR = 0dB),

the estimate of the parameters is poor what is due

to the bad estimate of third and fourth order cumu-

lants. These estimation results could be improved by

processing more data, which allows to get better esti-

mates of the cumulants.

Table 3: System identiﬁcation using signal model MA(2),

Mean, Standard Deviation, NMSE, and ESR. SNR=0dB.

True Values LS GDA GNA

λ = 8. 10

−5

µ = 0.5

h(1) -2.3333 mean -1.3781 -0.8631 -1.0952

σ 0.6033 0.6223 1.0662

h(2) 0.6667 mean 0.3217 0.3006 0.2882

σ 0.1794 0.2092 0.3499

NMSE mean 0.2418 0.4623 0.4964

σ 0.2464 0.2762 0.4890

ESR -6.1661 -3.3511 -3.0421

Table 4: System identiﬁcation using signal model MA(5),

Mean, Standard Deviation, NMSE, and ESR. SNR=20dB.

True Values LS GDA GNA

λ = 6. 10

−9

µ = 0.5

h(1) 0.1 mean 0.0474 0.0779 0.0887

σ 0.0359 0.0384 0.0392

h(2) -1.87 mean -1.4212 -1.2822 -1.7798

σ 0.0845 0.0665 0.1010

h(3) 3.02 mean 2.2489 1.6874 2.8902

σ 0.1258 0.0667 0.1494

h(4) -1.435 mean -1.0779 -1.1768 -1.4206

σ 0.0671 0.0754 0.0876

h(5) 0.49 mean 0.3640 0.3927 0.4536

σ 0.0275 0.0301 0.0320

NMSE mean 0.0651 0.1484 0.0046

σ 0.0221 0.0184 0.0075

ESR -11.8660 -8.2867 -23.3527

By increasing the order q of the model, we use the

third and fourth order cumulants with large lags that

are poorly estimated for a low SNR. The consequence

of this is the bad parameters estimated (table 6) much

poorer than in the case of a model of order 2 (table 3).

These two algorithms are numerically expensive

compared to the Least-Squares algorithm.

5 CONCLUSION

In this paper, a blind identiﬁcation of the MA mo-

dels using Higher-Order Statistics (HOS) is exposed.

The linear algebra solution is compared with the non-

linear optimization algorithms solution. Computer si-

mulation results prove that it is interesting to use this

algorithms in spite of their expensive calculative cost.

REFERENCES

Abderrahim, K., Ben Abdennour, R., Favier, G., Ksouri, M.,

and Msahli, F. (2001). New Results on FIR Identiﬁca-

tion using Cumulants. Journal Européen des Systèmes

Automatisés-JESA, Vol. 35, No. 5, pp. 601-622.

Alshebeili, S. A., Venetsanopoulos, A. N., and Cetin, A.

E. (1993, April). Cumulant Based Identiﬁcation Ap-

proaches for Nonminimum Phase FIR Systems. IEEE

Transactions on Signal Processing, Vol. 41, No. 4, pp.

1576-1588.

ICINCO 2005 - SIGNAL PROCESSING, SYSTEMS MODELING AND CONTROL

Giannakis, G. B., and Mendel, J. M. (1989, March). Iden-

tiﬁcation of Nonminimum Phase Systems Using Hi-

gher Order Statistics. IEEE Transactions on Acous-

tics, Speech, and Signal Processing, Vol. 37, No. 3,

pp. 360-377.

Martin, J. K., and Nandi, A. K. (1996, January). Blind Sys-

tem Identiﬁcation using Second, Third and Fourth Or-

der Cumulants. Journal of the Franklin Institute : En-

ginnering and applied mathematics, Pergamon, Vol.

333(B), No. 1, pp. 1-13.

Mendel, J. M. (1991, March). Tutorial on Higher-Order

Statistics (Spectra) in Signal Processing and System

Theory : Theoretical Results and some Applications.

Proceedings of the IEEE, Vol. 79, No. 3, pp. 278-305.

Na, Y. N., Kim, K. S., Song, I., and Kim, T. (1995, Decem-

ber). Identiﬁcation of Nonminimum Phase FIR Sys-

tems Using the Third- and Fourth-Order Cumulants.

IEEE Transactions on Signal Processing, Vol. 43, No.

8, pp. 2018-2022.

Nikias, C. L., and Mendel, J. M. (1993, July). Signal Pro-

cessing With Higher Order Spectra. IEEE Signal Pro-

cessing Magazine, pp. 10-37.

Srinivas, L., and Hari, K. V. S. (1995, December). FIR Sys-

tem Identiﬁcation Using Higher Order Cumulants - A

Generalized Approach. IEEE Transactions on Signal

Processing, Vol. 43, No. 12, pp. 3061-3065.

Stogioglou, A. G., and McLaughlin, S. (1996, July). MA

Parameter Estimation and Cumulant Enhancement.

IEEE Transactions on Signal Processing, Vol. 44, No.

7, pp. 1704-1718.

Tugnait, J. K. (1990, July). Approaches to FIR System Iden-

tiﬁcation With Noisy Data Using Higher Order Statis-

tics. IEEE Transactions on Signal Processing, Vol. 38,

No. 7, pp. 1307-1317.

Tugnait, J. K. (1991, October). New Results on FIR System

Identiﬁcation Using Higher-Order Statistics. IEEE

Transactions on Signal Processing, Vol. 39, No. 10,

pp. 2216-2221.

Table 5: System identiﬁcation using signal model MA(5),

Mean, Standard Deviation, NMSE, and ESR. SNR=10dB.

True Values LS GDA GNA

λ = 6. 10

−9

µ = 0.5

h(1) 0.1 mean -0.0929 -0.0841 -0.0604

σ 0.1145 0.1193 0.1161

h(2) -1.87 mean -1.0274 -0.9433 -1.2443

σ 0.2466 0.2177 0.3461

h(3) 3.02 mean 1.6917 1.4036 2.1096

σ 0.3516 0.2240 0.5103

h(4) -1.435 mean -0.8187 -0.8680 -0.9787

σ 0.1748 0.1971 0.2740

h(5) 0.49 mean 0.2760 0.2902 0.3151

σ 0.0653 0.0692 0.0887

NMSE mean 0.2122 0.2694 0.1311

σ 0.1105 0.0959 0.1345

ESR -6.7328 -5.6965 -8.8246

Table 6: System identiﬁcation using signal model MA(5),

Mean, Standard Deviation, NMSE, ESR. SNR=0dB.

True Values LS GDA GNA

λ = 6. 10

−9

µ = 0.5

h(1) 0.1 mean -0.2698 -0.2888 -0.3315

σ 0.1083 0.1154 0.1085

h(2) -1.87 mean -0.1297 -0.1187 -0.0863

σ 0.1577 0.1459 0.0811

h(3) 3.02 mean 0.3457 0.3302 0.2627

σ 0.2129 0.2021 0.1582

h(4) -1.435 mean -0.1644 -0.1630 -0.0949

σ 0.1068 0.1061 0.0845

h(5) 0.49 mean 0.0420 0.0427 0.0237

σ 0.0366 0.0359 0.0321

NMSE mean 0.8191 0.8280 0.8732

σ 0.1353 0.1300 0.0936

ESR -0.8667 -0.8199 -0.5887

0 10 20 30 40 50 60

−2.3

−2.28

−2.26

−2.24

−2.22

−2.2

Gradient descent algorithm − h(1)=−2.3333 (True value of h(1))

0 10 20 30 40 50 60

−2.27

−2.26

−2.25

−2.24

−2.23

−2.22

−2.21

−2.2

Gauss−Newton algorithm − h(1)=−2.3333 (True value of h(1))

Estimated parameter h(1)

Figure 2: The convergence of h(1).

0 10 20 30 40 50 60

0.64

0.65

0.66

0.67

0.68

0.69

0.7

Gradient descent algorithm − h(2)=0.6667 (True value of h(2))

0 10 20 30 40 50 60

0.56

0.58

0.6

0.62

0.64

0.66

Gauss−Newton algorithm − h(2)=0.6667 (True value of h(2))

Estimated parameter h(2)

Figure 3: The convergence of h(2).

PARAMETER ESTIMATION OF MOVING AVERAGE PROCESSES USING CUMULANTS AND NONLINEAR

OPTIMIZATION ALGORITHMS