- Email: [email protected]
Cumulant-based deconvolutionand identification: several new families of linear equations
Signal Processing 30 (1993) 199 219 Elsevier
199
Cumulant-based deconvolution and identification" Several new families of linear equations Fu-Chun Zheng, Stephen McLaughlin and Bernard Mulgrew Department of Electrical Engineering, University of Edinburgh, Edinburgh, UK Received 13 June 1991 Revised 2 March 1992
Abstract. This paper presents several new families of cumulant-based linear equations with respect to the inverse filter coefficients for deconvolution (equalisation) and identification of nonminimum phase systems. Based on noncausal autoregressive (AR) modeling of the output signals and three theorems, these equations are derived for the cases of 2nd-, 3rd and 4thorder cumulants, respectively, and can be expressed as identical or similar forms. The algorithms constructed from these equations are simpler in form, but can offer more accurate results than the existing methods. Since the inverse filter coefficients are simply the solution of a set of linear equations, their uniqueness can normally be guaranteed. Simulations are presented for the cases of skewed series, unskewed continuous series and unskewed discrete series. The results of these simulations confirm the feasibility and efficiency of the algorithms. Zusammenfassung. In der vorliegenden Arbeit werden einige auf Kumulanten basierende lineare Gleichungen angegeben in Hinblick auf die Berechnung der Koeffizienten inverser Filter zur Entfaltung (Entzerrung) und Identifikation von nicht minimalphasigen Systemen. Ausgehend yon einem nichtkausalen autoregressiven (AR) Modell des Ausgangssignals und drei Theoremen werden diese Gleichungen fiir Kumulanten zweiter, dritter und vierter Ordnung hergeleitet; sie k6nnen in identischer oder ~hnlicher Form ausgedriickt werden. Die von diesen Gleichungen abgeleiteten Algorithmen sind einfacher in der Form, fiihren aber zu genaueren Resultaten als vorhandene Methoden. Da die Koeffizienten des inversen Systems die L6sung eines linearen Gleichungssystems darstellen, kann die Eindeutigkeit im allgemeinen garantiert werden. Es werden Simulationsresultate fiir asymmetrisch verteilte Folgen, symmetrisch kontinuierlich und symmetrisch diskret verteilte Folgen wiedergegeben. Die Resultate dieser Simulationen bestfitigen die Zuverl/issigkeit und die Effizienz der Algorithmen. R6sum6. Cet article pr6sente plusieurs nouvelles families d'6quations lin~aires bas6es sur les cumulants satisfaites par les coefficients du filtre inverse pour la d6convolution (6galisation) et l'identification de syst6mes ~ phase non minimale. Ces 6quations, bas+es sur une mod+lisation autor6gressive (AR) des signaux de sortie non causale et trois th6or6mes, sont d6rivbes pour les cas de cumulants d'ordre deux, trois et quatre, respectivement, et peuvent ~tre exprim6es avec des formes identiques ou semblables. Les algorithmes construits li partir de ces 6quations ont une forme plus simple mais peuvent fournir des r6sultats plus precis que les m6thodes existantes. Du fait que les coeflficients du filtre inverse constituent simplement la solution d'un ensemble d'6quations lin6aires, leur unicit6 peut normalement ~tre garantie. Des simulations sont pr6sent6es pour les cas de s6ries asym6triques, de s6ries sym6triques continues, et de s~ries discr6tes sym6triques. Les r6sultats de ces simulations confirment la faisabilit6 et l'efficacit6 des algorithmes. Keywords. Deconvolution; identification; higher-order cumulants.
Correspondence to: Dr. Stephen McLaughlin, Department of Electrical Engineering, The University of Edinburgh, King's Building, Mayfield Road, Edinburgh EH9 3JL, UK. 0165-1684/93/$05.00 © 1993 Elsevier Science Publishers B.V. All rights reserved
200
F.-C Zheng et al. / Cumulant-based deconvolution and identification
1. Introduction
Over the last decade, increasingly more attention has been paid to deconvolution and identification of nonminimum phase systems. A powerful tool for this has emerged: higher-order spectral (HOS) analysis [4, 9, 13, 14, 16]. Higher-order spectra are mathematically defined as the Fourier transforms of higher-order cumulants (HOC). One of the main characteristics of HOC, which is vital to the deconvolution and identification of the nonminimum phase systems, is that the phase property of the signals is preserved. It is well known that the autocorrelation function (i.e., 2nd-order cumulants) is completely blind to the phase property of signals, and thus cannot be used to determine the phase characteristics of nonminimum systerrLs.As in ordinary spectral analysis, higher-order spectral analysis techniques may be divided into two categories: nonparametric (i.e., conventional or Fourierbased) [3, 16, 20] and parametric [4, 7, 9, 12 14, 16, 19, 21]. Parametric estimators can generally achieve higher resolution and lower variance, but require determination of model order [16, 20]. Although there is no need for order selection in conventional estimators, it is difficult to obtain high resolution and low variance using these types of estimators, especially for the case of relatively short data record length. As to model order selection required in parametric approaches, several HOC-based algorithms have been suggested so far [4, 8-10, 13, 15, 26], and in most cases the system order can be determined. Thus, parametric techniques are generally perferred to nonparametric ones. In this paper, an alternative set of parametric deconvolution and identification algorithms is presented. All these algorithms are derived directly from three theorems relating the inverse coefficients to the diagonal slices of cumulants of different orders. The algorithms requirements for the input are only that it be an independent and identically distributed (IID) non-Gaussian random series. The system can be nonminimum or minimum phase, but must have no zeros on the unit circle; 3rd- and Signal Processing
4th-order cumulants are employed respectively in our algorithms. In addition, although our algorithms are mainly oriented to FIR (MA) systems, IIR (ARMA) systems can also be considered. The current situation of parametric identification and deconvolution techniques for nonminimum phase systems has resulted in two basic approaches: MA model based methods and noncausal AR model based methods. One of the advantages of MA based approaches is that MA coefficients can be directly determined from closedform formula in most cases [7, 9, 19]. Unfortunately, the practical results of such closed-form formula are normally too inaccurate for the following two reasons. Firstly, these formula are oversensitive to the estimation error of the cumulants, and secondly, estimates of higher-order cumulants are generally of a relatively higher variance than that of autocorrelation functions [ 1]. An alternative to this is to reduce estimation of MA parameters to a nonlinear minimisation problem as in [6, 9, 12, 21, 26], however, this is a nontrivial problem. There are several existing schemes, such as good initial estimate finding methods [12, 21, 26] and overparametrisation methods [6, 9], but they do not always work well. Simulated annealing algorithms [ 11, 17] can achieve the global minimum and have recently interested many researchers in the optimisation community. However, they do not show much potential in the deconvolution and identification domain, as a result of their slow rate of convergence. As will be shown later, the problem ofmultimodality does not arise for our algorithms. For noncausal AR model based algorithms [4, 14, 16], the approach considered in this paper, the problem is no longer multimodal since a set of overdetermined linear equations results. In addition, since the identification results of this type of algorithms are simply the coefficients of the corresponding inverse filters, it is very convenient to implement deconvolution. In [3, 14], both diagonal and nondiagonal cumulant slices are used. The latter however has two main drawbacks: firstly, the estimation of nondiagonal cumulant slices is generally more noisy than that of diagonal slices
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
[ 1], thus using nondiagonal slices may cause degradation to the performance of the algorithms, and secondly, the forms of the resulting matrices are very complicated, especially when 4th-order cumulants are used. In order to overcome these two points, in our algorithms, only diagonal cumulant slices and autocorrelation functions (2nd-order cumulants) are employed. As a result, both the algorithms involved with the 4th-order cumulants and their counterparts using the 3rd-order cumulants take an identical simple form. Thus, deconvolution of an unskewed series for the case of N M P systems, which has not received as much attention as the case of skewed series, is realised in this paper as easily as that of a skewed series. These features comprise an important aspect of this paper's contribution. On the other hand, as is well known, the Giannakis-Mendel (GM) formula proposed in [9] represented significant progress, and holds an important position, in cumulant based identification theory. Unfortunately, the existence of the nonlinear terms in G M formula forms an obvious flaw: e.g., contradictory results may be generated when the overparameterisation scheme is implemented. In this paper, as a result of the introduction of three new theorems, we have cancelled, rather than only concealed (as in overparameterisation), the nonlinear terms in the G M formula by means of an AR model. As a consequence, this paper offers a feasible solution to the nonlinearity problem existing with the G M formula, and in the meantime opens a new efficient application channel for the G M formula. The general structure of this paper is as follows. In Section 2, the definitions of deconvolution and necessary assumptions are presented. In Section 3, we first derive three theorems regarding the inverse filter coefficients, and then develop the new deconvolution algorithms, each of them comprising a set of linear equations with respect to the inverse filter coefficients. The cases of 2nd-, 3rd- and 4th-order cumulants are studied, respectively. Some remarks are also given in this section. In Section 4, some computer simulation examples are presented. A
201
comparison between our algorithms and the method proposed by Nikias and Chiang in [4, 14] is also performed in Section 4.1. Finally, some conclusions are drawn in Section 5. In addition, some details of the derivation are contained within Appendix A. A concise report of some of the results of this paper can be found in [24].
2. Formulation of blind deconvolution problem The problem of blind deconvolution can be described as follows. Let x(k) be an ergodic series, which is the output of a stable linear time-invariant (LTI) mixed phase MA system S={beli = 0, 1. . . . . q}. Consequently, we must have q
x ( k ) = Z b,w(k-i),
(2.1)
i=0
where q denotes the order of system S, and w(k) is the independent and identically distributed (IID) driving series. Then, the following procedure is termed blind deconvolution: Find an inverse filter S -1 = { Oil i = r l , . . . , r2} which makes r2
w(k) ~ ~ Oix(k- i),
(2.2)
i= rl
given the output series x(k) only, where r~ and r2 is the properly selected order of the noncausal part and causal part, respectively. But we will write ' ~ ' as ' = ' in the following derivation since the distinction is clear. Let us notice that rl = 0 when S is minimum phase and rz = - 1 when S is maximum phase. Thus minimum phase and maximum phase are two special cases of the above model. It is also assumed that system S has no zeros on the unit circle. In addition, it is well known that x(k) is not deconvolvable if S is nonminimum phase and w(k) is a Gaussian IID series [2, 16]. Hence, we always assume in this paper that w(k) is a non-Gaussian IID series with E[w(k)] = ~'~ = 0,
(2.3a)
E[w2(k)] = ~'2,
(2.3b)
E[w3(k)] = ~'3
(2.3c) Vol. 30. No. 2. January 1993
202
F.-C Zheng et aL / Cumulant-based deconvolution and identification
and
2nd-order cumulants, be expressed as
E[w4(k)] - 3E2[w2(k)] = ~4,
(2.3d)
where ~/4~ 0 as a result of the assumption that w(k) is non-Gaussian. In (2.3c), 7'3 can be zero for nonskewed input signal, and in this case we may use the fourth-order cumulants. It should be noted that the assumption that the time series is non-Gaussian is often realistic in practical situations, as indicated in [22, 23].
3. Algorithm derivation
1 Zr2 Ojcz(i+J) = Ibm, ~'2j=r, [0,
C2(mO = E [ x ( k ) x ( k + m0],
(3.1a)
C3(ml'm2)=EIx(k)~x(k+mi)]'i=l
(3.1b)
other.
(3.3)
PROOF. See Appendix A. T H E O R E M 2. I f system S is LTI, causal and stable, and has no zeros on the unit circle, and its input w(k) is zero-mean, skewed and liD, then the relation between S and its inverse filter S -~ can, in terms o f 3rd-order eumulants, be expressed as 1
For the output series x ( k ) defined in last section, its nth-order cumulants can be denoted as Cn(ml,m2 . . . . . m,_~). From the definition of cumulants [16], we have
if ie [0, q],
~ Ojc3(i+j)=IbZi,
~3 j=rl
~0,
i f i e [ 0 , q], other.
(3.4)
PROOF. See Appendix A.
Notice that only a limited number of equations are nontrivial in (3.3) and (3.4). Specifically, for an MA(q) system, c,(m) = 0
for m e ( - ~ ,
4,ml - C2(ml)C2(m3-m2) - C2(m2) C2(m3 - ml) - C2(m3)C2(m2-ml). (3.1c)
Specially, we denote the C(ml, m2 . . . . . mn-l) as
diagonal
slice
of
c,(m) = Cn(m,, m2, . . . , m,-,)l,~, =m2. . . . . . . ,=m.
(3.2) The algorithms proposed are based mainly on the following theorems. T H E O R E M 1. I f system S is LTI, causal and stable, and has no zeros on the unit circle, and its input w(k) is zero-mean and liD, then the relation between S and its inverse filter S 1 can, in terms o f Signal Processing
-q-
1] w [q+ 1, +oo),
(3.5)
where n = 2, 3, 4, thus it can be seen that the equations in (3.3) and (3.4) are nontrivial only when i+rj<<.q and i+r2>~-q, viz., -q-r2<<.i<~q-rl [1, 3, 12]. Thus, there are r z - r l + 2 q + 1 nontrivial equations in (3.3) and (3.4), respectively. Since there are only r2 - rx + 1 coefficients 0is to be determined, many sets of equations can be constructed with respect to the 0is. As a result of the his being unknown now, we can adopt either of the following two schemes to eliminate them. S C H E M E 1. Only adopting the zero equations i.e., the ones with right-hand sides equal to zero in (3.3) and (3.4).
The number of such equations 2 ( r 2 - r l +q) is still much greater than the number of unknowns r2 - rl + 1. Hence, we need only take the following
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
2 ( r 2 - r l - q) equations in order to reduce unnecessary computation: r2
O]c2(i+j) =0, j = rl r2
Ojc3(i+J) =0, j = rl
i= -rz .....
- 1, q + 1. . . . .
-rl.
(3.6)
SCHEME 2. Eliminating bis in (3.3) and (3.4). In this scheme, both zero and non-zero equations in (3.3) and (3.4) are adopted, but bis are eliminated by using the formula proposed by Giannakis and Mendel in [9]. In [9], Giannakis and Mendel derived the following equations (i.e., eq. (24) in [9], but it is slightly modified here) :
Without loss of generality, we can always let 00 = 1. Viz., the inverse filter coefficients have been normalised by 00 here. Then, (3.6) can be written as
AO=B,
c2( r2+rO c3(--r2+rl)
q
c2(-r2-1)
cz(-r2+ 1)
c3(-r2-1)
c3(-r2+1)
c2(0) c3(0)
i=0
A=
c3(--l+r 0
c2(q+l+rO c3(q+l + r 0
c2(0) c3(01
cz(q) ca(q)
c2(0) c3(0)
c2(-1+r2)
cz(q+2)
c2(q+l+r2) c3(q+l+r2)
c2(-rt-l)
c2(--r,+l)
c3(-r,-1)
c3(--r,+l)
j=rl
, 0rf
i=0 r2
q
= ~ Oj ~ c3(n-i)c2(i+j), .j=rl
f2(g2--Fl) c3(r2-r,)
(3.8a 0--1, 01,..
q
E Oj • c2(n--i)c3(i+j)
c3(-1+r2)
c3(q+2)
(3.10)
Applying the results of Theorem 1 and Theorem 2 to (3.10), we can obtain r2
c2(-2) c~(-2)
O=[Orl, 0rl+l .....
i=0
n = --q . . . . . 2q.
i c2(--l+r0
q
?/3 ~ c2(n-i)b 2= 9~2 ~ c3(n-i)bi,
(3.7)
where
203
nE[-q, 2q].
i=0
(3.11/ As before, 00 can still be assumed to be 1. Then, (3.11) can be rewritten as
(3.8b)
r2
2 f.jOj=g., ne[-q, 2q],
and
(3.12)
j=rl j~O
B=-[c2(-r2), ¢3(-r2)
. . . . .
c2(-1), where
c3(-1), c2(q+l), c3(q+l) . . . . , c2(--rl), c3(--rl)] T.
Notice that (3.7) is still overdetermined since rl and r2 can always be taken to be sufficiently large to make 2(r2 - rl - q) ~>r2 - rl, viz., r2 - rl/> 2q. Thus, taking the least-squares solution of (3.7), the inverse filter coefficients can be expressed as 0 = ( A T A ) - 1( A T B ) .
q
(3.8c)
fnj = ~ [c2(n- i)c3(i+j) - c3(n - i)c2(i+j)] i=0
and q
g,= ~ [c3(n-i)cz(i)-c2(n-i)c3(i)]. i=0
( 3 . 9)
Since AVA is positive definite, we must have det(ATA) > 0 [5]. Hence, O can be uniquely determined from (3.9).
Clearly, many sets of equations with respect to 0is can be formed from (3.12) and the zero equations in (3.3) and (3.4), but the following two are the simplest. Vol. 30, No. 2, January 1993
204
F.-C
Z h e n g et al. / C u m u l a n t - b a s e d
Firstly, combining (3.3) with (3.12), we have
d e c o n v o l u t i o n a n d identification
where
r2
Ojc2(n + j ) = -c2(n)
E
U = YL(I)I/-3
(3.19a)
V= YR(I)I,=3.
(3.19b)
and
j=r I j#O
for ne[-r2, -1] w [q+ 1, -rl],
(3.13)
Then, we must have
r2
E
Ojf,j = g,
for n ~ [0, q].
0 = (uTu)-I(uTv).
j=rl j4=O
Let "ct(-r2+rl)
i el(l+&)
"''
ct(--r2-1)
"-"
: c/(-2)
cd0)
c1(r2+l) : cl(0)
c~(-1+ r2)
(3.20)
If the input series w(k) is unskewed, viz., T3 = 0, the 4th-order cumulants must be used in order to deconvolve and identify the nonminimum phase systems. In this case, we employ the following theorem.
fO,r2 YL(/) = fq.r, "'" el(q+ 1 +r~) • . •
cd0)
'''
fq. , cl(q)
cff-rl-1)
fq., cz(q+ 2) ct(-rt+l)
fq,r cl(q + 1 + r2) c,(r2 - r,)
(3.14a)
T H E O R E M 3. I f system S is LTI, causal and stable, and has no zeros on the unit circle, and its input w(k ) is zero-mean, non-Gaussian and liD, then the relation between S and its inverse filter S -1 can, in terms of 4th-order cumulants, be expressed as L
~ Y4j=n
and YR(I) = [-c,(-r2) . . . . . - c t ( - 1 ) , go . . . . .
Ojc4(i+j)=I b~' (0,
ifi6[0, q],
(3.21)
other.
gq,
PROOF. See Appendix A. - c , ( q + 1) . . . . . - c t ( - r l ) ] v.
(3.14b)
Then (3.13) can be more compactly written as PO=Q,
(3.15)
where P = YL(I)I/-2
From this theorem and foregoing discussion, the following descriptions are straightforward. Firstly, combining Theorem 3 with Theorem 1, and along the line of above Scheme 1, the linear equation systems with the same form as (3.7) can be obtained'
(3.16a)
and
AO=B,
(3.22)
where, however, Q= YR(1)]l=2.
-c2(-rz+&)
(3.16b)
c4(-r2 + r,)
Thus, we can determine O = (pTp)-l ( p r o ) .
(3.17)
Secondly, combining (3.4) and (3.12), we can obtain the following equation in a similar way: UO = V, Signal Processing
(3.18)
A=
''
c2(-r2-1)
• " c4(-r2
cz( r 2 + l )
1) c 4 ( - r 2 + 1 )
c2(0) c4(0)
c2(--1 +rt)
C2(--2)
C2(0)
c4(-1 +r,)
c4(--2)
£'4(0)
c4(-- 1 + r2)
c2(q+ 1 + r 0
c=(q)
c2(q + 2)
c2(q+ 1 +r2)
c4(q+ 1 + r 0
c4(q)
c4(q + 2)
c4(q+ 1 +r=)
c2(-r,+l) c4(-&+l)
c2(r2 - r,) c.(r:-r,)
c2(0) c,(0)
• " c2(-r~-l) " c4(-rt-1)
C2(--I +r2)
(3.23a)
205
F.-C. Zheng et al. Cumulant-baseddeconvolution and identification and
and B = - [c2(-r2),
c4(-r2)
.....
c2(-
1 ),
Q = ZR(/)I,=2.
c4(-1), c2(q+ 1), c4(q+ 1) . . . .
Then, corresponding to (3.18), we can also obtain
c2(-rl), c4(-rl)] T.
(3.23b)
Secondly, considering Theorem 3 and the Giannakis-Mendel formula (eq. (29a) in [9]) for the case of 4th-order cumulants (with a slight modification here), q
(3.27b)
UO= V,
(3.28)
where U=
ZL(l
)Ii = 4
(3.29a)
and
q
~4 E c2(n--i) b3= Y2 ~ c4(n-i)bi, i=0
V = ZR(I)It=4.
i=0
n = -q .....
2q,
(3.24)
the following two equations - (3.26) and (3.28), with the same form as (3.15) and (3.18), respectively, can be derived. Now let ct(-r2+rO
c/(--r2--1)
cl(--r2+l)
"'"
ct(0)
el(-1 +r0
c1(--2)
ct(0)
•""
ct(-- 1 + r2)
J~,r I
fo. ,
fo.,
"'"
fo,~2
i
:
Zc(l)= fq,r I ct(q+ 1 +rO
fo,, ct(q)
i ct(O)
c¢(-r~-l)
f~,' ct(q + 2)
:
"'"
fq,'2
• • • ct(q + 1 + r2)
: cl(-rt+l)
: ---
ct(rz-r 0
(3.25a) and ZR(I) = [--c/(--r2) . . . . .
- c , ( q + 1) . . . . .
--Ct(--1), g o . . . . .
gq,
- c t ( - r 0 ] T,
(3.25b)
where f,j and gn are, respectively, given by q
fnj = ~ [c2(n - i)c4(i + j ) - c4(n - i)c2(i + j ) ] i=0
and q
g, = ~', [c4(n- i)c2(i) - c2(n - i)c4(i)]. i=0
Corresponding to (3.15), we have P O = Q,
(3.26)
where
P = Zc( l)ll- 2
(3.27a)
(3.29b)
An interesting point to note is that the results of Theorems 1, 2 and 3 can actually be expressed as the following general form:
1 ~ Ojcz(i+j)=Ibl-', Ylj=~, [0,
ifi~[0, q], other.
(3.30)
Here we have assumed Yt ¢ 0, l = 2, 3, 4. A logical question to ask here is: can the above formula (3.3.30) be generalised to the case of any lth-order cumulants (i.e., l = 2, 3, 4, 5, 6 , . . . ) ? Our computer simulations (i.e., empirically or numerically) have given an affirmative answer. However, a mathematical proof for the case of l/> 5 still remains to be found (the proofs for the cases of l = 2, 3 and 4 are presented in Appendix A). On the other hand, it should be pointed out that the above question has only theoretical significance since there has been no report on an application where the 5th- or still higher-order cumulants have been used. Notice that the equations involved with the 4thorder cumulants are of exactly the same form as their counterparts using the 3rd-order cumulants. As a result, the 4th-order cumulant based algorithms can be implemented almost as easily as the 3rd-order cumulant based ones. Furthermore, it can be found that all the families of linear equations we have obtained above are of a similar elegant form, which makes the programming of the algorithms very regular. Regarding other specific aspects of our algorithms, we have the following remarks. Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
206
REMARK 1. Since all the above six families of equations (i.e., (3.7), (3.15), (3.18), (3.22), (3.26) and (3.28)) are linear, the inverse filter coefficients {Ojlj=rl . . . . . 0 . . . . . r2} can normally be determined uniquely. Notice that AVA, pVp and UTU are positive definite, and therefore nonsingular. In the above equations, c/(" ) will be replaced in a practical implementation by its estimate ~/(. ) (as in (4.2) and (4.3) in Section 4). It can be proved that the estimator ~t(" ) is both unbiased and consistent (provided the noise effects are removed) under the following two conditions [9]: (i) the input w(k) is stationary and zero-mean, and its first six cumulants are absolutely summable; and (ii) the system is stable. Obviously, these conditions are very loose and thus can usually be satisfied. Considering that the 0is are measurable functions of ~z(. ) in the above sets of linear equations, the 0~s themselves must be consistent. This is because the measurable function of consistent estimator(s) is also consistent (for details, see [9] and the references therein).
REMARK 2. All the above algorithms are derived from the MA (or FIR) systems. In the case of A R M A (including AR, similarly hereinafter), there are in principle two schemes to apply the above formulae: (i) Since any stable LTI system can be approximated by a truncated MA model, we can use the truncated impulse response his of A R M A systems as the parameters of the approximating MA models, an idea used in [7, 26]. The advantage of this method is that the order of the AR part of A R M A model is not needed. But on the other hand, this method requires a relatively higher accuracy of cumulant estimation. (ii) As is well known, for an ARMA(p, q) process P
q
x(k)+ ~ a i x ( k - i ) = ~ biw(k-i), i~l
p>~q,
i~O
(3.31) the AR coefficients ais can be determined from the equations below: Signal Processing
p
c3(-n)+ ~ aic3(-n+i)=O,
n>q,
(3.32)
i=l
where the AR order p can be selected by using other algorithms (for example, Hankel matrix based methods) (see [8, 9, 18] for details). Then, the residual process p
y(k) = x(k) + ~ aix(k- i),
(3.33)
i--1
where ~i denotes the estimate of ai, reduces to an MA process: q
y(k) = Z blw(k-i),
(3.34)
i-O
to which the MA based algorithms can be applied.
REMARK 3. Although the algorithms proposed above are derived from the viewpoint of deconvolution, they may conveniently be employed for system identification. As soon as the 0is are determined, the MA parameters bis (for MA systems), or the truncated impulse response his (for A R M A systems), can be directly calculated up to a gain constant by using one of the following two schemes. TH1 scheme. Directly employ Theorem 1 (eq. (3.3)). This scheme is very suitable for the case of low additive noise. When the level of additive noise is very high, the obtained bis or his may be rather poor estimates since the estimation of c2(0) can be severely degraded. INV scheme. Invert the obtained S - l = {0~}. This scheme can generate satisfactory results no matter whether the level of the additive noise is low or high, as will be demonstrated later in the simulation section. The reason why Theorem 2 or 3 is not adopted here is that the estimation of higher-order cumulants (especially the 4th-order cumulants) normally has a relatively higher variance, thus direct employment of Theorem 2 or 3 could deteriorate the estimation accuracy ofbis or his. In addition, Theorem 3 cannot give the signs of bis (or his).
F.-C Zheng et al. / Cumulant-based deconvolution and identification
Obviously, when our algorithms are used to identify the MA parameters of systems, the problem of multimodality no longer exists, and the uniqueness of identification is guaranteed. To this extent, the paper provides a feasible solution to the nonlinearity problem which exists with GiannakisMendel (GM) formula [9]: in [9], the nonlinear terms are only concealed by using overparameterisation, but in our derivation above, they are removed by introducing the three theorems and adopting an AR model. In the meantime, our algorithms obviously reflect a fresh channel of application of the GM formula to identification and deconvolution. 4. Since both the 3rd-order cumulants and autocorrelation functions are used, the above algorithms can deal with both non-minimum and minimum phase systems for deconvolution and system identification. The only constraint is that the system to be deconvolved may not have any zeros on the unit circle. This can be considered as the penalty paid for the linearity of the equations. REMARK
R E M A R K 5. Selection of the inverse filter order rl
and r2. The algorithms proposed above are insensitive to the choice of rl and r2, as long as they are taken to be large enough. However, the above algorithms only involve the general case of mixed phase systems. For the special cases of minimum and maximum phase systems, several minor changes must be made to ensure that the above algorithms are numerically more reliable. (i) For minimum phase systems rl should be set to be zero. Although this is not necessary theoretically, rl <0 does sometimes cause large errors, especially in the presence of a high level of additive noise. (ii) Similarly, for maximum phase systems, r2 should be set to be -1. Then, 0_~, instead of 00, should be assumed to be 1. The corresponding changes also need to be made for all families of linear equations ((3.7), (3.15), (3.18), (3.22), (3.26) and (3.28)).
REMARK
207
6. Determination of MA order q. This
problem depends on the phase property of the considered system. (i) For minimum phase systems, the above algorithms are robust to the choice of q value. The reason for this is that the matrix structures in all above sets of equations become irrelevant with q value when systems are minimum phase. (ii) For nonminimum phase systems, the above algorithms are generally sensitive to the selection of MA order q, since the value of q in this case does directly affect the actual contents of the matrices in all sets of equations. Thus, reliable order determination is required. In fact, several higher-order cumulant based methods for order selection have emerged in recent years [4, 8-10, 13, 15, 26]. But in this paper we do not discuss this problem, but assume that the system order q is known. We will also demonstrate these two points by using simulation results in Section 4. R E M A R K 7. A theoretical analysis of the performances of above six families of linear equations ((3.7), (3.15), (3.18), (3.22), (3.26) and (3.28)) is not available at this time. Heuristically, for the case of 3rd-order cumulants, (3.15) and (3.18) are more stable than (3.7), especially when the method (i) mentioned previously in Remark 2 is used for ARMA systems. However, if it is noted that (3.18) involves more 3rd-order cumulant estimation, which is robust to additive Gaussian noise, it is recommended that (3.18) be used as long as a sufficiently long output series can be obtained to give an acceptable 3rd-order cumulant estimation; certainly, (3.15) can also be tried if the SNR is not too low. If SNR is too low, (3.15) can give less accurate results since it involves more autocorrelation terms. Similar comments can also be made for the case of 4th-order cumulants. In a subsequent section of this paper where we consider simulations, we will empirically confirm these points. R E M A R K 8. All the above algorithms can be eas-
ily made recursive. Our work, which is reported in [25, 27], has shown that the convergence rate of Vol. 30, No. 2, January 1993
208
F-C. Zheng et al. / Cumulant-based deconoolution and identification
the recursive versions of the algorithms is very satisfactory.
4. Simulation examples In this section, we apply the algorithms described above to MA (FIR) system deconvolution. Both nonminimum and minimum phase systems are considered. To study the case of 3rd-order cumulants, an exponentially distributed series is used as the skewed input series w(k). For the case of 4th-order cumulants, a uniformly distributed continuous series and an equally distributed discrete series are used as the unskewed input series w(k). Additionally, as in the above algorithm derivation, w(k) is assumed to be IID in all examples below. The output series x(k) is generated by convolving w(k) with the system impulse response. The additive white noise n(k) to x(k) is assumed to be Gaussian, and the signal-to-noise ratio (SNR) in this paper is defined as SNR = 10 l o g , 0 ( ~ / . \E[n (k)]]
(4.1)
In the simulations below, the lth order cumulants c;(m) in the foregoing equations is substituted with their corresponding estimates d;(m). In order to obtain a more accurate cumulant estimation, the well known segment-average scheme [9, 16] was adopted. For convenience, we briefly describe this scheme as follows. First, segment K output samples into L records of N samples each, viz., K = NL. Then, estimate the cumulants ~i)(m) for the ith record by .,(/) - - 1 n~ ct ( m ) - ~ ~ x(i)(n)[xU)(n + m)];- 1,
(4.2)
~ n=nl
where x(°(n) denotes the samples in the ith record. In addition, here, i = 1. . . . . L, n;= max(0, - m ) and nu = m i n ( N - 1, N - m - 1). Finally, average ~°(m) over all records:
et(m) Signal Processing
=
I L ;=I ~°(m).
(4.3)
In this section, we always take K = N x L = 300 x 30, and perform 20 Monte-Carlo (MC) runs for each example. It should be noted that the above segment-average procedure is not essential for the ergodic series. The more detailed analysis regarding these points is contained in [25, 27] where a recursive version of the algorithms is derived. The 'True' curves in the following diagrams are determined by directly inverting the corresponding system transfer function H(z), i.e., directly expanding 1/H(z). Also, in the following simulations, only the MA systems with known order are considered. In the case of A R M A systems, we can first use the methods mentioned in Remark 2 to transform A R M A problems to MA problems. Then, the application of the MA identification and deconvolution algorithms to the obtained MA models becomes straightforward.
4.1. Skewed w(k): the case of 3rd-order cumulants Here, w(k) is taken to be exponentially distributed (y2 = 1, Y3= 2).
E X A M P L E 1. Mixed phase MA system: the model is x'(k) = w(k) - 2.4w(k - 1) + 0.8w(k- 2) + n(k),
(4.4)
with two zeros: zl = 2 and z2 = 0.4. Set r 1 = - 1 5 and r2=15. For the noise levels which produce SNR = 50 dB and SNR = 10 dB, the first MC run is then implemented: (3.7), (3.15) and (3.18) are solved, respectively. The obtained inverse filter coefficients 0is are shown in Figs. (la), (lb) and (lc), respectively. It can be seen that, for this MA system, (i) the results of three equations are of approximately identical accuracy at the lower noise level (SNR = 50 dB), and
209
IE-C. Zheng et al. / Cumulant-based deconvolution and identification 3.0
(a)
I
I
I
I
3.0
2.5-
2.5-
2.0-
2.0-
1.5-
1.5-
1.0-
1.0-
0.5-
0.5-
0.0
0.0
-0.5
I -10
-15
]5
-
I 0
j
I 5
I 10
-0.5 15
Fig. la. Results of (3.7) in Example 1: mixed phase MA(2) system. I
3.0
I
I
2.52.01.51.00.50.0 i
-15
-10
-
0
j
I
I
I
= -10
-
15
0
j
I 5
I 10
15
Fig. lc. Results o f (3.18) in E x a m p l e 1 : m i x e d phase M A ( 2 ) system.
I
(b)
-0.5
-15
I
(c)
i
i
5
10
15
Fig. lb. Results of (3.15) in Example 1: mixed phase MA(2) system.
(ii) the result of (3.15) is of a relatively higher error in comparison with (3.7) and (3.18) for the higher noise level (SNR = 10 dB). The reason for the latter is that, as mentioned in Remark 7, (3.15) involves relatively more autocorrelation (2nd-order cumulant) terms, among which c2(0) is subject to the effect of the additive Gaussian noise. Nineteen other MC runs are similarly accomplished. Only the general identification results (mean+variance, 20 MC runs, similarly hereinafter) are listed in Table 1. Both the two MA identification schemes (TH1 and INV) mentioned in
Remark 3 are used here. Also, for comparison, the g a i n a ~ u s t m e n t has been done for every MC run: be = bi/bo, where b;s are results given by TH1 and INV schemes. For the sake of space, we did not include b0 in Table 1 since 60 = 1.00000 ± 0.00000 (=b0). The same operations will also be used in all the following examples. From Table 1, it is clear that the TH1 scheme can generate very accurate results only when SNR is not too poor, but the INV scheme can satisfactorily secure the accuracy of the results not only in the high SNR case but also in the low SNR case. The latter actually also indicates that the output, 0ts, of our deconvolution algorithms is very robust. As will be seen, similar observation can be obtained in all of the following examples. In addition, the sensitivity of the algorithms to the order mismatch is tested. Figure (ld)illustrates the results of (3.7), (3.15) and (3.18) when the system order is set to be 3 and SNR = 50 dB. It can be seen that the algorithms are sensitive to the selected system order in the case of N M P systems, as explained in Remark 6. Finally, for comparison, the method suggested by Nikias and Chiang in [4, 14] (we refer to this method .as Nikias-Chiang (N-C) method in this paper) is implemented using the same data. For this example, the number of equations in the N - C Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
210 Table 1
Identification results (mean :t: variance) in Example 1
SNR
Identification scheme
50dB
THI
INV
10dB
THI
INV
Equations
/~]
b2
(3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18 )
-2.405864-0.11309 -2.389334-0.12451 -2.398474-0.12663 -2.441144-0.11344 -2.45467 4- 0.14289 -2.440324-0.12528 -4.327594-0.80167 -4.11011+0.62873 -4.28709 4- 0.79086 -2.168844-0.10448 -2.09096 4- 0.13401 -2.09915 4- 0.12463 -2.40000
0.82076:t:0.05539 0.820404-0.06225 0.815354-0.05617 0.820194-0.04714 0.82767 4- 0.06580 0.81559:t:0.04452 0.409414-0.21913 0.284824-0.30937 0.23018 4- 0.37202 0.773874-0.06390 0.68703 4- 0.07342 0.74242 4- 0.06348 0.80000
True bi
,:3.0
J
I
I
I
I
/
--
i
........ Eq(3.7) .... Eq(3.15)
1~
2.0-
~
1.00.5-
I
I
2.0-
...... Eq(a.IR)
1.51.00.50.0
0.0 i -10
15 -
i 0
j
I 5
i 10
-0.5 -15
15
Fig. ld. Sensitivity o f the algorithms to order mismatch: mixed phase M A ( 2 ) system, algorithm o r d e r = 3.
method is ( r 2 - r O ( q + l ) = 9 0 . Figure l(e) is the inverse filter {0t} generated by the N - C method in the first MC run (SNR = 50 dB and 10 dB). Nineteen other MC runs of this method are then performed. The INV scheme is similarly used in all these MC runs, and the obtained MA parameters (corresponding to Table 1) are listed in Table la. Clearly, the results o f our algorithms are more accurate than those of the N - C method, while our algorithmic forms are much simpler. E X A M P L E 2. Minimum phase MA system: the model is Signal Processing
J
2.5-
True
2
1.5-
-15
I
(e)
2.5-
-0.5
1
3.0
(d)
-10
-
J5
i
0
j
i
5
i
10
15
Fig. le. Results of Nikias-Chiang method in Example 1 : mixed phase MA(2) system.
x'(k) = w(k) + 0 . 3 w ( k - 1) - 0 . 4 w ( k - 2) + n(k),
(4.5)
with two zeros: z] = - 0 . 8 and z2=0.5. Table la M A parameters (mean±variance) obtained from the N - C method in Example 1 SNR
/~1
62
50 dB 10 dB True bi
-3.07428 :~ 0.22263 -3.13369 4- 0.26429 -2.40000
1.32381 4- 0.16828 1.37363 4- 0.20516 0.80000
211
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
First, set r~ = 0 and r 2 = 15. Adding the noise which makes S N R = 10 dB and S N R = 1 dB this time, the first MC run is realised. F r o m (3.7), (3.15) and (3.18), the corresponding inverse filters are determined, which are drawn in Figs. (2a), (2b) and (2c), respectively. It is clear that the result of (3.18) is better than that of (3.15) and (3.7) when the additive noise is very strong (1 dB, here), as noted in Remark 7. Similarly, nineteen other MC runs are implemented, and only the corresponding
I
ii
I
I
I
I
I
I
I
I
I
(c)
0.0-
0
(a)
I
~- 0.5--
--0.5
1.5
I
I
I
I
I
2
4
6
8 j
I
I
I
10
12
14
Fig. 2c. Results of (3.18) in Example 2: minimum phase MA(2) system.
10
~- 0.5-
0.0-
-0.5 0
I
I
I
2
4
6
I
8 j
I
t
I
10
12
14
Fig. 2a. Results of (3.7) in Example2: minimum phase MA(2) system•
~- 0.5-
identification results are listed in Table 2. In terms of accuracy, they are comparable with the results listed in Table 1 although the SNR is much lower in this example. To this extent, the algorithms offer a stronger robustness to the additive noise for the MP systems than for the N M P ones. Next, let us carry out two more tests. First, set rl = - 1 5 and r2 = 15. Figure (2d) demonstrates the obtained inverse filters (SNR = 1 dB). As pointed out in Remark 5, serious errors appear in the anticausal parts of the inverse filters. But the result of (3.18) is still more satisfactory than that of (3.7) and (3.15). In fact, during the simulation of other examples, not demonstrated here, we observed much more serious errors. Thus, for minimum systems, it would be better to take rl = 0. Then, the system order is taken to be 3, and the obtained inverse filters from the three algorithms are illustrated in Fig. (2e), where r~ = 0, r2 = 15 and SNR = 10 dB. Clearly, the algorithms are robust to the order mismatch in the case of MP systems, as previously mentioned in Remark 6.
0.0-
-0.5
0
I
I
1
I
2
4
6
8 j
I
I
I
10
12
14
Fig. 2b. Results o f (3.15) in Example 2: minimum phase MA(2) system.
From the above two examples, it can be seen that all three equations, (3.7), (3.15) and (3.18), are robust to the presence of additive noise, all producing satisfactory performances, though there are slight differences in accuracy. Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
212
Table 2 Identification results (mean 4-variance in Example 2 SNR
Identification scheme
10dB
TH1
Equations (3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18)
INV
1 dB
6~
THI
INV
-0.25984-4-0.01515 0.23703+0.01519 0.26651 +0.01572 0.29561 ±0.01864 0.24355±0.01564 0.30816+0.02161 0.12033-4-0.01275 0.09951 ±0.01292 0.15129+0.01455 0.185964-0.01553 0.090824-0.01496 0.31768±0.03964 0.30000
True bi
I
1.5
I
I
I
b2
I
I
1.5
-........ .... ........
1.0-
I
I
I
L%
-0.5
-lo
-'5
1
I Eq(3.15) [q(~.le)
1.0-
/ L
0.5-
0.0-
I
o
]
i
i
5
10
-0.5
15
Fig. 2d. Algorithm comparison: r, is set to be -15, minimum phase MA(2) system, SNR = 1 dB.
4.2. Unskewed and continuous w(k)." the case o f the 4th-order cumulants
o
I
I
I
I
2
4
6
8
j
3. M i x e d p h a s e d M A
s y s t e m : the
I
i
I
10
12
14
Fig. 2e. Insensitivity of the algorithms to order mismatch: minimum phase MA(2) system, algorithm order = 3.
x ' ( k ) = w ( k ) - 2 . 4 w ( k - 1) + 0 . 8 w ( k - 2) + n(k),
T h e d i s t r i b u t i o n o f w ( k ) is t a k e n to be u n i f o r m over the interval [ - 2 , 2] (72 = 4 / 3 , 73 = 0 a n d 74 = - 3 2 / 1 5 ) . F o r c o m p a r i s o n , the same system transfer functions as in Section 4.1 are c o n s i d e r e d in the following examples.
Signal Processing
I
EqO.7)
~-
0.0
EXAMPLE m o d e l is
I
-Troe
True Cq(3.7) Eq(3.15) Eq(3.1e)
0.5-
-15
I
(e)
(d)
~-
-0.35709+0.00916 -0.34979±0.00904 -0.35989±0.00930 -0.38188±0.00976 -0.38219±0.01009 -0.384474-0.00983 -0.194204-0.00584 -0.19016±0.00562 -0.20296±0.00620 -0.28536±0.01466 -0.24499±0.01016 -0.36954±0.02220 -0.40000
(4.6)
with two zeros: 2" 1 = 2 a n d Z 2 = 0.4. Set r~ = - 1 5 a n d r2 = 15. N o i s e which results in S N R = 50 dB a n d S N R = 15 d B is a d d e d to the o u t p u t series, respectively. F i r s t the M C run is implem e n t e d : solving (3.22), we o b t a i n the inverse filters s h o w n in Fig. (3a); next, b y solving (3.26), the inverse filters s h o w n in Fig. (3b) is o b t a i n e d ;
F.-C. Zheng et al. / Cumulant-based deconvolution and identification 3.0
I
(a)
I
I
I
I
3.0 2.5-
2.0-
2.0-
1.5-
1.5-
1.0-
1.0-
0.5-
0.5-
0.0
0.0
I -15
15
-10
Fig. 3a. Results
-
I 0
of (3.22)
j
in Example system.
-0.5
I
I
5
10
3: mixed
15
phase M A ( 2 )
I
I
I
(c)
2.5-
-0.5
213
/\ t
-15
-10
t5 -
i 0
Fig. 3c. R e s u l t s o f (3.28) in E x a m p l e system.
j
i
t
5
10
15
3: m i x e d p h a s e MA(2)
EXAMPLE 4. Minimum phase M A system: the model is I
3.0
I
I
I
I
x'(k) = w(k ) + 0 . 3 w ( k - 1)
(b) 2.5-
- 0 . 4 w ( k - 2) + n(k),
1 r
2.01.51.00.50.0-
-0.5
i -15
-10
15 -
i 0
Fig. 3b. R e s u l t s o f (3.26) in E x a m p l e system.
j
i
i
5
10
15
3 : mixed phase MA(2)
finally, (3.28) is solved, and the determined inverse filters are demonstrated in Fig. (3c). Also, nineteen other MC runs are realised, and the corresponding identification results are given in Table 3. Notice that, in Table 3, the corresponding variance becomes higher than in Table 1. This is due to the fact that estimation of the 4th-order cumulants normally has higher variance than that of the 3rdorder cumulants. Similar explanation also applies to the following Examples 4 and 5.
(4.7)
with two zeros: Zl = - 0 . 8 and z2 = 0.5. Set rl = 0 and r2 = 15. Noise which results in SNR = 50 dB and SNR = 15 dB is added, respectively. The determined inverse filters in the first MC run by (3.22), (3.26) and (3.28) are shown in Figs. (4a), (4b) and (4c), respectively. Similarly, the general identification results for the twenty MC runs are listed in Table 4. Clearly, this table verifies our observation in Example 2 again: the algorithms are of a stronger robustness to the additive noise for the MP systems than for the N M P ones. It can be seen that all the results obtained by the algorithms in the above two examples are satisfactory, but the results from (3.28) are always better than the ones from the other two equations. As to the corresponding errors, their major sources include additive noise, cumulant estimation errors and numerical errors. 4.3. Unskewed and discrete w(k)." the case of the
4th-order cumulants w(k) is taken to be an equally distributed discrete series over set {4-1, +3} (72 = 5, ~ 4 = 26). The additive noise level is set to be 20 dB. Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
214
Table 3 Identification results (mean ± variance) in Example 3
SNR
Identification scheme
50dB
TH1
INV
15 dB
THI
INV
Equations
b,
/~2
(3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28)
-2.39696±0.26939 -2.365874-0.26746 -2.37549:t:0.27624 -2.471554-0.25671 -2.476325:0.26435 -2.46264±0.25440 -3.00481 5:0.58958 -2.883345:0.52864 -2.93120 ± 0.55689 -2.32376 ± 0.28544 -2.311804-0.31493 -2.31759 4-0.29774 -2.40000
0.82085-4-0.14106 0.813585:0.14183 0.81180±0.14116 0.81589±0.15489 0.80997± 0.16752 0.80254+0.16107 0.724265:0.22767 0.695024-0.24923 0.70533+ 0.23539 0.80739+ 0.17903 0.76410+0.20260 0.78962+ 0.18820 0.80000
True bi
10
I
1.5
1.5 (a) I
I
I
I
I ~
I
I
I
I
I
C
(b)
J
-1.0-
~- 0.5-
"~ 0.5-
0.0-
0.0-
-0.5 0
I
I
I
I
2
4
6
B j
I
I
I
10
12
14
Fig. 4a. Resultsof(3.22) in Example4:minimum phase MA(2) system.
EXAMPLE
5. M i x e d p h a s e d M A
system" the
m o d e l is x ' ( k ) = w ( k ) - 2.4w(k - 1)
+ 0 . 8 w ( k - 2) + n ( k ) ,
(4.8)
with t w o zeros: z, = 2 a n d z2 = 0.4. O b v i o u s l y , this is a typical channel equal±sat±on p r o b l e m . Set r, = - 1 5 a n d r2 = 15. In the presence o f 20 dB additive noise, first M C r u n is c a r r i e d out, a n d the results are as follows. U s i n g (3.22): result Signal P r o c e s s i n g
--0.5
True 50dB 15dB
........ ........
I
I
I
I.
2
4
6
8 j
i
I
I
10
12
14
Fig. 4b. Results of (3.26) in Example 4: minimum phase MA(2) system.
s h o w n in Fig. (5a). U s i n g (3.26): result shown in Fig. (5b). U s i n g (3.28): result shown in Fig. (5c). Then, nineteen similar M C runs are a c c o m p l i s h e d , a n d the c o r r e s p o n d i n g general identification results are listed in T a b l e 5. These five e x a m p l e s s h o w t h a t o u r a l g o r i t h m s can efficiently deal with the p r o b l e m o f d e c o n v o l u tion a n d identification o f b o t h N M P a n d M P systems, while the a l g o r i t h m i c f o r m s are v e r y simple a n d regular. A s to the identification o f M A p a r a m eters, the I N V scheme is very secure n o m a t t e r the noise level is high o r low.
215
F.-C. Zheng et al. / Cumulant-based deconvolution and identification 1.5 ]
I
I
I
I
I
I
,3.0
I
(a)
I
I
I
I
2.51.01
2.0-
~
1.51.0
oolV -0.51 0
0.5 0.0'
I
i
I
I
2
4
6
8
j
i
T
i
10
12
14
-0.5
Fig. 4c. Results of(3.28) in Example4': minimum phase MA(2) system.
15
i
-15
-10
-
i
0
j
i
i
5
10
15
Fig. 5a, Results of (3.22) in Example 5: mixed phase MA(2) system.
Table 4 Identification results (mean 4-variance) in Example 4 SNR
Identification scheme
50 dB
TH1 INV
15dB
THI INV
True bi
Equations (3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28)
6,
62
-0.295604-0.01701 0.29571 +0.01850 0.29571 4-0.01697 0.300634-0.02037 0.30299-4-0.02157 0.29970+0.02161 0.281654-0.01716 0.27771 4-0.01862 0.28433+0.01731 0.294924-0.02049 0.288604-0.02118 0.29824 4-0.02273 0.30000
5. Conclusions Through the above discussion, six families of linear equations have been derived according to three theorems, and they are all of a similar simple form. As a result of the linearity of these equations, the uniqueness of solutions can normally be guaranteed. Thus, the risk of local minima, which exists with some other methods, is removed. In comparison with the existing linear (i.e., noncausal A R based) approaches, the algorithms
-0.404214-0.01607 -0.40408i 0.01639 -0.404414-0.01603 -0.396384-0.01670 -0.402534-0.01610 -0.393685:0.01814 -0.38793~:0.01652 -0.38586+0.01693 -0.389554-0.01652 -0.39428+0.01781 -0.394714-0.01628 -0.396784-0.01984 -0.40000
proposed in this paper possess the following two features: (a) Only the diagonal cumulant slices are used, which in general makes the algorithms more accurate, and also makes the forms of the concerned matrices simple and elegant; (b) 4th-order cumulants are employed in the last three algorithms exactly as the substitute of 3rdorder cumulants in the first three algorithms, which enables the 4th-order cumulant based algorithms (for deconvolution of an unskewed series) to be Vol. 30, No. 2, January 1993
216
F.-C. Zheng et al. / Cumulant-based deconvolution and identification I
3.0
I
I
I
implemented as easily as their 3rd-order cumulants based counterparts (for deconvolution of a skewed series).
I
(b) 2.52,01.5-
Appendix A
i 1.0-
In this appendix, the symbol system used in the text will be followed.
0.50.0 -0.5
I
15
-lO
-15
-
; o
j
A. 1. Proof of Theorem 1
t
i
5
10
15
From the conditions in Theorem 1 and (2.1), we have Fig. 5b. Results of (3.26) in Example 5: mixed phase MA(2) system. 3.0
I
I
I
I
• .,x(i) = E [w(k)x(k + i)]
I
= E [w(k).=~0 b.w(k+i-n)]
(c)
2.5-
(A1.1) (A1.2)
q
= ~ b.E[w(k)w(k+i-n)] 2.0-
(A1.3)
n--0 q
= ~ b.~'26(n- i)
1.5-
(A1.4)
n=0 1.0-
= 0.5-
I~'2bi, if ie[0, q], (0,
0.0
(A1.5)
other.
On the other hand, from (2.2), we have
-0.5
I -15
-I0
15 -
i 0
j
, 5
10
15
• wx(i) = E [ w ( k ) x ( k +
Fig. 5c. Results of (3.28) in Example 5: mixed phase MA(2) sytem.
= E [ j ~ , Ojx(k-j)x(k+i) 1
Table 5 Identification results (mean + variance) in Example 5
SNR
Identification scheme
20dB
TH1
INV
True bi Signal Processing
i)]
Equations
bl
b2
(3.22) (3.26) (3.28) (3.22) (3.26) (3.28)
-2.52986-4-0.40044 -2.48994-4-0.38909 -2.50339-4-0.39954 -2.395934-0.20816 -2.396274-0.22162 -2.38865 4- 0.20661 -2.40000
0.77463+0.11616 0.76475+0.12562 0.765274-0.11811 0.78276+0.11999 0.76706+0.13470 0.77008 4- 0.12863 0.80000
(A1.1) (A1.6)
217
F.-C Zheng et al. / Cumulant-baseddeconvolutionand identification A.3. Proof of Theorem 3
I"2
= ~ 0iEIx(k-j)x(k+i)]
(A1.7)
j=rl r2
= E Ojc2(i+j).
From the conditions in Theorem 3 and (2.1), we have
(A1.8)
j--rl
. . . . (i) = E[w(k)x3(k + i)]
Equating (A1.5) and (A1.8), Theorem 1 holds.
[]
= E w(k)
b~w(k + i - n n
I
A.2. Proof of Theorem 2
[
q b.w(k + i - n
E = w(k
[
(A3.2)
0 q
q
~ ~ bmbnblw(k+i-m)
m=0 n=0 /=0
×w(k+i-n)w(k+i-l)]
(A2.1) (A2.2)
q
q
q
= ~
~
~ bmbnbl
rn=0n=0
n
E:
q
=E w(k) ~
From the conditions in Theorem 2 and (2.1), we have Cbwxx(i)= E[w(k)x2(k + i)]
(A3.1)
(A3.3)
/=0
× E[w(k)w(k + i - m)
q'
×w(k+i-n)w(k+i-l)].
2 Zb b.
(A3.4)
m = O n=O
For the input series w(k), from (3.1c), we have × w(k + i - n)w(k + i - m)] q
=~"
(A2.3) E[w(k) iO, w(k+mi)]
q
Ebmbn
m = 0 n=O
= C~(ml, m2, m3)
x E[w(k)w(k + i - m)w(k + i - n)] (A2.4) q
= ~
+ )'22[~(ml)~(m3-m2)
q
~ bmbnT3~(m-i, n - i )
(A2.5)
+ ~(me)g(m3 - m l ) + 3(m3)g(m2-ml)],
m-0 n=0
~T3b~, if is[0, q], ~0, other.
(A3.5) (A2.6) where C~(ml, m2, m3) is the 4th-order cumulant of w(k). According to [3, 9],
On the other hand, from (2.2), we have
C~(ml,
• wxx(i) = E[w(k)x2(k + i)]
m2,
m3) = ~/4~(ml ,
m2, m 3 ) .
(A3.6)
(A2.1) Then, from (A3.4), we derive
=E[j~Ojx(k-j)x2(k+i)l
(A2.7)
r2
= ~ OjE[x(k-j)xZ(k+i)]
(A2.8)
j=rl r2
= ~ Ojc3(i+j).
. . . . (i) q
= E
q
q
E E b.,b.b,
m=0 n=0 /=0
(A2.9)
j=rl
Equating (A2.6) and (A2.9), Theorem 2 holds.
x [T46(i-m, i - n , i - l )
+ T~{g ( i - m)g(n - l) + g ( i - n)fi(m - 1) []
+ ~ ( i - l)~(m - n) } ]
(A3.7)
Vol. 30, No. 2, January 1993
218
F-C Zheng et al. / Cumulant-based deconvolution and identification 2
/ if i t [0, q],
9,
2
btbi+ E bmbi÷ E bZ,b, ,
)/463÷)'2
m=O
n=O
other (A3.8)
= y4b3+372b, 2 b], n=0 q
t [0,
i f i e I 0 , q], other. (A3.9)
O n the o t h e r h a n d , f r o m (2.2) a n d ( Y l c ) , we have
. . . . (i) = E[w(k)x3(k + i)] =E[
(A3.1)
Ojx(k-j)x3(k+i)]
~
(A3.10)
t-j=rl r2
= Z OjE[x(k-j)x3(k+i)]
(A3.11)
j = rl r2
= ~
Oj[c4(iWj) ÷ 3 c 2 ( i ÷ j ) c 2 ( 0 ) ] (A3.12)
j=rl r2
= ~. Ojc4(i+j)
(A3.13)
j rl r2 +3c2(01 ~ Ojc2(i+j). j--r1
N o t i c i n g T h e o r e m 1 (3.3) a n d q
c2(0)=72 ~ b],
(A3.14)
n=0
we c a n o b t a i n I r2
q
~" O j c 4 ( i + j ) + 3 y 2 b i ~ bZ. n=0
~j=r I
.... (i)=1
rifiE[0'q]'
[ E Ojc4(i+j),
other.
~j=rl
(A3.15) E q u a t i n g (A3.9) a n d (A3.15), we m u s t have r2 Ojc4(i+j) = I r4b3, j=,~ (0, which is j u s t T h e o r e m 3. SignalProcessing
if i t [ 0 , q], other, []
References
(A3.16)
[1] S. Bellini et al., "Blind deconvolution: Polyspectrum or Bussgang techniques?", in: E. Biglieri et al., eds., Digital Communications, Elsevier, Amsterdam, 1986, pp. 251 263. [2] A. Benveniste et al., "Robust identification of a nonminimum phase system: Blind adjustment of a linear equaliser in data communications", IEEE Trans. Automat. Control, Vol. AC-25, No. 3, June 1980, pp. 385 399. [3] D.R. Brillinger and M. Rosenblatt, "Computation and interpretation of kth-order spectra", in: B. Harris, ed., Spectral Analysis of Time Series, Wiley, New York, 1987, pp. 189 232. [4] H.-H. Chiang and C.L. Nikias, "Adaptive deconvolution and identification of nonminimum phase FIR systems based on cumulants", IEEE Trans. Automat. Control, Vol. AC-35, January 1990, pp. 36-47. [5] A.R. Collar and A. Simpson, Matrices and Engineering Dynamics, Ellis Horwood, Chichester, UK, 1987, Chapter 1. [6] B. Friedlander and B. Porat, "Asymptoticallyoptimal estimation of MA and ARMA parameters of non-Gaussian processes from higher-order moments", IEEE Trans. Automat. Control, Vol. AC-35, January 1990, pp. 27 35. [7] G.B. Giannakis, "Cumulants: A powerful tool in signal processing", Proc. IEEE, Vol. 75, No. 9, September 1987, pp. 1333-1334. [8] G. Giannakis and J.M. Mendel, "ARMA order determination via higher-order statistics", in: C.L. Martin, ed., Mathematical Theory of Networks and Systems, NorthHolland, Amsterdam. [9] G.B. Giannakis and J.M. Mendel, "Identification of nonminimum phase system using higher order statistics", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP37, March 1989, pp. 360 377. [10] G.B. Giannakis and J.M. Mendel, "Cumulant-based order determination of non-Gaussian ARMA models", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, August 1990, pp. 1411 1423. [11] S. Kirkpatrick, C.D. Gelalt and M.P. Vecchi Jr., "Optimization by simulated annealing", Science, Vol. 220, No. 4598, 13 May 1983, pp. 671-680. [12] K.S. Lii and M. Rosenblatt, "Deconvolution and estimation of transfer function phase and coefficients for nonGaussian linear process", Ann. Statist., Vol. 10, 1982, pp. 1195 1208. [13] C.L. Nikias, "ARMA bispectrum approach to nonminimum phase system identification", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-36, April 1988, pp. 513-524. [14] C.L. Nikias and H.-H. Chiang, "Higher-order spectrum estimation via noncausal autoregressive modeling and deconvolution", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-36, December 1988, pp. 1911-1913. [15] C.L. Nikias and R. Pan, "ARMA modeling of fourthorder cumulants and phase estimation", Circuit Systems Signal Processing, Vol. 7, No. 3, 1988, pp. 291-325.
F.-C. Zheng et al. / Cumulant-based deconvolution and identification [16] C.L. Nikias et al., "Bispectrum estimation: A digital signal processing framework", IEEE Proc., Vol. 75, No. 7, July 1987, pp. 869-891. [17] W.H. Press et al., Numerical Recipes - The Art of Scientific Computing, Cambridge Univ. Press, Cambridge, 1986, Chapter 10. [18] M.R. Raghuveer and C.L. Nikias, "Bispectrum estimation: A parametric approach", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-33, 1985, pp. 12131230. [19] A. Swami and J.M. Mendel, "Closed-form recursive estimation of MA coefficients using autocorrelations and third-order cumulants", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-37, No. 11, November 1989, pp. 1794-1795. [20] A.M. Tekalp and A.T. Erdem, "Higher-order spectrum factorization in one and two dimensions with applications in signal processing and nonminimum phase system identification", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-37, October 1. ~'~ op. 1537-1549. [21] K.S. Tugnait, "Identification of linear stochastic systems via second- and fourth-order cumulant matching", IEEE Trans. Inform. Theory, Vol. IT-33, May 1987, pp. 393407. [22] A.T. Walden and J.W.J. Hosken, "An investigation of the spectral properties of primary reflection coefficients", Geophysical Prospecting, Vol. 33, 1985, pp. 400-435.
219
[23] A.T. Walden and J.W.J. Hosken, "The nature of the nonGaussianity of primary reflection coefficients and its significance for deconvolution", Geophysical Prospecting, Vol. 34, 1986, pp. 1038-1066. [24] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind deconvolution algorithms based on 3rd- and 4th-order cumulants", Proc. IEEE Internat. Conf. Acoust. Speech Signal Process. (ICASSP), Toronto, Ontario, Canada, 14-17 May 1991, pp. 1573-1576. [25] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "A 2nd- and 4th-order cumulant based blind equalisation algorithm for nonminimum phase channels", Proc. IEEE & EURASIP lnternat. Signal Processing Workshop on Higher Order Statistics, Chamrousse, France, 10-12 July 1991, pp. 79 82. [26] F.-C Zheng, S. McLaughlin and B. Mulgrew, "Blind equalization of nonminimum phase channels: Higherorder cumulant based algorithm", IEEE Trans. Signal Process., to appear. [27] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind equalisation of multilevel PAM data for nonminimum phase channels via 2nd- and 4th-order cumulants", Signal Processing, Vol. 31, No. 3, April 1993, to appear.
Vol. 30, No. 2, January 1993
199
Cumulant-based deconvolution and identification" Several new families of linear equations Fu-Chun Zheng, Stephen McLaughlin and Bernard Mulgrew Department of Electrical Engineering, University of Edinburgh, Edinburgh, UK Received 13 June 1991 Revised 2 March 1992
Abstract. This paper presents several new families of cumulant-based linear equations with respect to the inverse filter coefficients for deconvolution (equalisation) and identification of nonminimum phase systems. Based on noncausal autoregressive (AR) modeling of the output signals and three theorems, these equations are derived for the cases of 2nd-, 3rd and 4thorder cumulants, respectively, and can be expressed as identical or similar forms. The algorithms constructed from these equations are simpler in form, but can offer more accurate results than the existing methods. Since the inverse filter coefficients are simply the solution of a set of linear equations, their uniqueness can normally be guaranteed. Simulations are presented for the cases of skewed series, unskewed continuous series and unskewed discrete series. The results of these simulations confirm the feasibility and efficiency of the algorithms. Zusammenfassung. In der vorliegenden Arbeit werden einige auf Kumulanten basierende lineare Gleichungen angegeben in Hinblick auf die Berechnung der Koeffizienten inverser Filter zur Entfaltung (Entzerrung) und Identifikation von nicht minimalphasigen Systemen. Ausgehend yon einem nichtkausalen autoregressiven (AR) Modell des Ausgangssignals und drei Theoremen werden diese Gleichungen fiir Kumulanten zweiter, dritter und vierter Ordnung hergeleitet; sie k6nnen in identischer oder ~hnlicher Form ausgedriickt werden. Die von diesen Gleichungen abgeleiteten Algorithmen sind einfacher in der Form, fiihren aber zu genaueren Resultaten als vorhandene Methoden. Da die Koeffizienten des inversen Systems die L6sung eines linearen Gleichungssystems darstellen, kann die Eindeutigkeit im allgemeinen garantiert werden. Es werden Simulationsresultate fiir asymmetrisch verteilte Folgen, symmetrisch kontinuierlich und symmetrisch diskret verteilte Folgen wiedergegeben. Die Resultate dieser Simulationen bestfitigen die Zuverl/issigkeit und die Effizienz der Algorithmen. R6sum6. Cet article pr6sente plusieurs nouvelles families d'6quations lin~aires bas6es sur les cumulants satisfaites par les coefficients du filtre inverse pour la d6convolution (6galisation) et l'identification de syst6mes ~ phase non minimale. Ces 6quations, bas+es sur une mod+lisation autor6gressive (AR) des signaux de sortie non causale et trois th6or6mes, sont d6rivbes pour les cas de cumulants d'ordre deux, trois et quatre, respectivement, et peuvent ~tre exprim6es avec des formes identiques ou semblables. Les algorithmes construits li partir de ces 6quations ont une forme plus simple mais peuvent fournir des r6sultats plus precis que les m6thodes existantes. Du fait que les coeflficients du filtre inverse constituent simplement la solution d'un ensemble d'6quations lin6aires, leur unicit6 peut normalement ~tre garantie. Des simulations sont pr6sent6es pour les cas de s6ries asym6triques, de s6ries sym6triques continues, et de s~ries discr6tes sym6triques. Les r6sultats de ces simulations confirment la faisabilit6 et l'efficacit6 des algorithmes. Keywords. Deconvolution; identification; higher-order cumulants.
Correspondence to: Dr. Stephen McLaughlin, Department of Electrical Engineering, The University of Edinburgh, King's Building, Mayfield Road, Edinburgh EH9 3JL, UK. 0165-1684/93/$05.00 © 1993 Elsevier Science Publishers B.V. All rights reserved
200
F.-C Zheng et al. / Cumulant-based deconvolution and identification
1. Introduction
Over the last decade, increasingly more attention has been paid to deconvolution and identification of nonminimum phase systems. A powerful tool for this has emerged: higher-order spectral (HOS) analysis [4, 9, 13, 14, 16]. Higher-order spectra are mathematically defined as the Fourier transforms of higher-order cumulants (HOC). One of the main characteristics of HOC, which is vital to the deconvolution and identification of the nonminimum phase systems, is that the phase property of the signals is preserved. It is well known that the autocorrelation function (i.e., 2nd-order cumulants) is completely blind to the phase property of signals, and thus cannot be used to determine the phase characteristics of nonminimum systerrLs.As in ordinary spectral analysis, higher-order spectral analysis techniques may be divided into two categories: nonparametric (i.e., conventional or Fourierbased) [3, 16, 20] and parametric [4, 7, 9, 12 14, 16, 19, 21]. Parametric estimators can generally achieve higher resolution and lower variance, but require determination of model order [16, 20]. Although there is no need for order selection in conventional estimators, it is difficult to obtain high resolution and low variance using these types of estimators, especially for the case of relatively short data record length. As to model order selection required in parametric approaches, several HOC-based algorithms have been suggested so far [4, 8-10, 13, 15, 26], and in most cases the system order can be determined. Thus, parametric techniques are generally perferred to nonparametric ones. In this paper, an alternative set of parametric deconvolution and identification algorithms is presented. All these algorithms are derived directly from three theorems relating the inverse coefficients to the diagonal slices of cumulants of different orders. The algorithms requirements for the input are only that it be an independent and identically distributed (IID) non-Gaussian random series. The system can be nonminimum or minimum phase, but must have no zeros on the unit circle; 3rd- and Signal Processing
4th-order cumulants are employed respectively in our algorithms. In addition, although our algorithms are mainly oriented to FIR (MA) systems, IIR (ARMA) systems can also be considered. The current situation of parametric identification and deconvolution techniques for nonminimum phase systems has resulted in two basic approaches: MA model based methods and noncausal AR model based methods. One of the advantages of MA based approaches is that MA coefficients can be directly determined from closedform formula in most cases [7, 9, 19]. Unfortunately, the practical results of such closed-form formula are normally too inaccurate for the following two reasons. Firstly, these formula are oversensitive to the estimation error of the cumulants, and secondly, estimates of higher-order cumulants are generally of a relatively higher variance than that of autocorrelation functions [ 1]. An alternative to this is to reduce estimation of MA parameters to a nonlinear minimisation problem as in [6, 9, 12, 21, 26], however, this is a nontrivial problem. There are several existing schemes, such as good initial estimate finding methods [12, 21, 26] and overparametrisation methods [6, 9], but they do not always work well. Simulated annealing algorithms [ 11, 17] can achieve the global minimum and have recently interested many researchers in the optimisation community. However, they do not show much potential in the deconvolution and identification domain, as a result of their slow rate of convergence. As will be shown later, the problem ofmultimodality does not arise for our algorithms. For noncausal AR model based algorithms [4, 14, 16], the approach considered in this paper, the problem is no longer multimodal since a set of overdetermined linear equations results. In addition, since the identification results of this type of algorithms are simply the coefficients of the corresponding inverse filters, it is very convenient to implement deconvolution. In [3, 14], both diagonal and nondiagonal cumulant slices are used. The latter however has two main drawbacks: firstly, the estimation of nondiagonal cumulant slices is generally more noisy than that of diagonal slices
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
[ 1], thus using nondiagonal slices may cause degradation to the performance of the algorithms, and secondly, the forms of the resulting matrices are very complicated, especially when 4th-order cumulants are used. In order to overcome these two points, in our algorithms, only diagonal cumulant slices and autocorrelation functions (2nd-order cumulants) are employed. As a result, both the algorithms involved with the 4th-order cumulants and their counterparts using the 3rd-order cumulants take an identical simple form. Thus, deconvolution of an unskewed series for the case of N M P systems, which has not received as much attention as the case of skewed series, is realised in this paper as easily as that of a skewed series. These features comprise an important aspect of this paper's contribution. On the other hand, as is well known, the Giannakis-Mendel (GM) formula proposed in [9] represented significant progress, and holds an important position, in cumulant based identification theory. Unfortunately, the existence of the nonlinear terms in G M formula forms an obvious flaw: e.g., contradictory results may be generated when the overparameterisation scheme is implemented. In this paper, as a result of the introduction of three new theorems, we have cancelled, rather than only concealed (as in overparameterisation), the nonlinear terms in the G M formula by means of an AR model. As a consequence, this paper offers a feasible solution to the nonlinearity problem existing with the G M formula, and in the meantime opens a new efficient application channel for the G M formula. The general structure of this paper is as follows. In Section 2, the definitions of deconvolution and necessary assumptions are presented. In Section 3, we first derive three theorems regarding the inverse filter coefficients, and then develop the new deconvolution algorithms, each of them comprising a set of linear equations with respect to the inverse filter coefficients. The cases of 2nd-, 3rd- and 4th-order cumulants are studied, respectively. Some remarks are also given in this section. In Section 4, some computer simulation examples are presented. A
201
comparison between our algorithms and the method proposed by Nikias and Chiang in [4, 14] is also performed in Section 4.1. Finally, some conclusions are drawn in Section 5. In addition, some details of the derivation are contained within Appendix A. A concise report of some of the results of this paper can be found in [24].
2. Formulation of blind deconvolution problem The problem of blind deconvolution can be described as follows. Let x(k) be an ergodic series, which is the output of a stable linear time-invariant (LTI) mixed phase MA system S={beli = 0, 1. . . . . q}. Consequently, we must have q
x ( k ) = Z b,w(k-i),
(2.1)
i=0
where q denotes the order of system S, and w(k) is the independent and identically distributed (IID) driving series. Then, the following procedure is termed blind deconvolution: Find an inverse filter S -1 = { Oil i = r l , . . . , r2} which makes r2
w(k) ~ ~ Oix(k- i),
(2.2)
i= rl
given the output series x(k) only, where r~ and r2 is the properly selected order of the noncausal part and causal part, respectively. But we will write ' ~ ' as ' = ' in the following derivation since the distinction is clear. Let us notice that rl = 0 when S is minimum phase and rz = - 1 when S is maximum phase. Thus minimum phase and maximum phase are two special cases of the above model. It is also assumed that system S has no zeros on the unit circle. In addition, it is well known that x(k) is not deconvolvable if S is nonminimum phase and w(k) is a Gaussian IID series [2, 16]. Hence, we always assume in this paper that w(k) is a non-Gaussian IID series with E[w(k)] = ~'~ = 0,
(2.3a)
E[w2(k)] = ~'2,
(2.3b)
E[w3(k)] = ~'3
(2.3c) Vol. 30. No. 2. January 1993
202
F.-C Zheng et aL / Cumulant-based deconvolution and identification
and
2nd-order cumulants, be expressed as
E[w4(k)] - 3E2[w2(k)] = ~4,
(2.3d)
where ~/4~ 0 as a result of the assumption that w(k) is non-Gaussian. In (2.3c), 7'3 can be zero for nonskewed input signal, and in this case we may use the fourth-order cumulants. It should be noted that the assumption that the time series is non-Gaussian is often realistic in practical situations, as indicated in [22, 23].
3. Algorithm derivation
1 Zr2 Ojcz(i+J) = Ibm, ~'2j=r, [0,
C2(mO = E [ x ( k ) x ( k + m0],
(3.1a)
C3(ml'm2)=EIx(k)~x(k+mi)]'i=l
(3.1b)
other.
(3.3)
PROOF. See Appendix A. T H E O R E M 2. I f system S is LTI, causal and stable, and has no zeros on the unit circle, and its input w(k) is zero-mean, skewed and liD, then the relation between S and its inverse filter S -~ can, in terms o f 3rd-order eumulants, be expressed as 1
For the output series x ( k ) defined in last section, its nth-order cumulants can be denoted as Cn(ml,m2 . . . . . m,_~). From the definition of cumulants [16], we have
if ie [0, q],
~ Ojc3(i+j)=IbZi,
~3 j=rl
~0,
i f i e [ 0 , q], other.
(3.4)
PROOF. See Appendix A.
Notice that only a limited number of equations are nontrivial in (3.3) and (3.4). Specifically, for an MA(q) system, c,(m) = 0
for m e ( - ~ ,
4,ml - C2(ml)C2(m3-m2) - C2(m2) C2(m3 - ml) - C2(m3)C2(m2-ml). (3.1c)
Specially, we denote the C(ml, m2 . . . . . mn-l) as
diagonal
slice
of
c,(m) = Cn(m,, m2, . . . , m,-,)l,~, =m2. . . . . . . ,=m.
(3.2) The algorithms proposed are based mainly on the following theorems. T H E O R E M 1. I f system S is LTI, causal and stable, and has no zeros on the unit circle, and its input w(k) is zero-mean and liD, then the relation between S and its inverse filter S 1 can, in terms o f Signal Processing
-q-
1] w [q+ 1, +oo),
(3.5)
where n = 2, 3, 4, thus it can be seen that the equations in (3.3) and (3.4) are nontrivial only when i+rj<<.q and i+r2>~-q, viz., -q-r2<<.i<~q-rl [1, 3, 12]. Thus, there are r z - r l + 2 q + 1 nontrivial equations in (3.3) and (3.4), respectively. Since there are only r2 - rx + 1 coefficients 0is to be determined, many sets of equations can be constructed with respect to the 0is. As a result of the his being unknown now, we can adopt either of the following two schemes to eliminate them. S C H E M E 1. Only adopting the zero equations i.e., the ones with right-hand sides equal to zero in (3.3) and (3.4).
The number of such equations 2 ( r 2 - r l +q) is still much greater than the number of unknowns r2 - rl + 1. Hence, we need only take the following
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
2 ( r 2 - r l - q) equations in order to reduce unnecessary computation: r2
O]c2(i+j) =0, j = rl r2
Ojc3(i+J) =0, j = rl
i= -rz .....
- 1, q + 1. . . . .
-rl.
(3.6)
SCHEME 2. Eliminating bis in (3.3) and (3.4). In this scheme, both zero and non-zero equations in (3.3) and (3.4) are adopted, but bis are eliminated by using the formula proposed by Giannakis and Mendel in [9]. In [9], Giannakis and Mendel derived the following equations (i.e., eq. (24) in [9], but it is slightly modified here) :
Without loss of generality, we can always let 00 = 1. Viz., the inverse filter coefficients have been normalised by 00 here. Then, (3.6) can be written as
AO=B,
c2( r2+rO c3(--r2+rl)
q
c2(-r2-1)
cz(-r2+ 1)
c3(-r2-1)
c3(-r2+1)
c2(0) c3(0)
i=0
A=
c3(--l+r 0
c2(q+l+rO c3(q+l + r 0
c2(0) c3(01
cz(q) ca(q)
c2(0) c3(0)
c2(-1+r2)
cz(q+2)
c2(q+l+r2) c3(q+l+r2)
c2(-rt-l)
c2(--r,+l)
c3(-r,-1)
c3(--r,+l)
j=rl
, 0rf
i=0 r2
q
= ~ Oj ~ c3(n-i)c2(i+j), .j=rl
f2(g2--Fl) c3(r2-r,)
(3.8a 0--1, 01,..
q
E Oj • c2(n--i)c3(i+j)
c3(-1+r2)
c3(q+2)
(3.10)
Applying the results of Theorem 1 and Theorem 2 to (3.10), we can obtain r2
c2(-2) c~(-2)
O=[Orl, 0rl+l .....
i=0
n = --q . . . . . 2q.
i c2(--l+r0
q
?/3 ~ c2(n-i)b 2= 9~2 ~ c3(n-i)bi,
(3.7)
where
203
nE[-q, 2q].
i=0
(3.11/ As before, 00 can still be assumed to be 1. Then, (3.11) can be rewritten as
(3.8b)
r2
2 f.jOj=g., ne[-q, 2q],
and
(3.12)
j=rl j~O
B=-[c2(-r2), ¢3(-r2)
. . . . .
c2(-1), where
c3(-1), c2(q+l), c3(q+l) . . . . , c2(--rl), c3(--rl)] T.
Notice that (3.7) is still overdetermined since rl and r2 can always be taken to be sufficiently large to make 2(r2 - rl - q) ~>r2 - rl, viz., r2 - rl/> 2q. Thus, taking the least-squares solution of (3.7), the inverse filter coefficients can be expressed as 0 = ( A T A ) - 1( A T B ) .
q
(3.8c)
fnj = ~ [c2(n- i)c3(i+j) - c3(n - i)c2(i+j)] i=0
and q
g,= ~ [c3(n-i)cz(i)-c2(n-i)c3(i)]. i=0
( 3 . 9)
Since AVA is positive definite, we must have det(ATA) > 0 [5]. Hence, O can be uniquely determined from (3.9).
Clearly, many sets of equations with respect to 0is can be formed from (3.12) and the zero equations in (3.3) and (3.4), but the following two are the simplest. Vol. 30, No. 2, January 1993
204
F.-C
Z h e n g et al. / C u m u l a n t - b a s e d
Firstly, combining (3.3) with (3.12), we have
d e c o n v o l u t i o n a n d identification
where
r2
Ojc2(n + j ) = -c2(n)
E
U = YL(I)I/-3
(3.19a)
V= YR(I)I,=3.
(3.19b)
and
j=r I j#O
for ne[-r2, -1] w [q+ 1, -rl],
(3.13)
Then, we must have
r2
E
Ojf,j = g,
for n ~ [0, q].
0 = (uTu)-I(uTv).
j=rl j4=O
Let "ct(-r2+rl)
i el(l+&)
"''
ct(--r2-1)
"-"
: c/(-2)
cd0)
c1(r2+l) : cl(0)
c~(-1+ r2)
(3.20)
If the input series w(k) is unskewed, viz., T3 = 0, the 4th-order cumulants must be used in order to deconvolve and identify the nonminimum phase systems. In this case, we employ the following theorem.
fO,r2 YL(/) = fq.r, "'" el(q+ 1 +r~) • . •
cd0)
'''
fq. , cl(q)
cff-rl-1)
fq., cz(q+ 2) ct(-rt+l)
fq,r cl(q + 1 + r2) c,(r2 - r,)
(3.14a)
T H E O R E M 3. I f system S is LTI, causal and stable, and has no zeros on the unit circle, and its input w(k ) is zero-mean, non-Gaussian and liD, then the relation between S and its inverse filter S -1 can, in terms of 4th-order cumulants, be expressed as L
~ Y4j=n
and YR(I) = [-c,(-r2) . . . . . - c t ( - 1 ) , go . . . . .
Ojc4(i+j)=I b~' (0,
ifi6[0, q],
(3.21)
other.
gq,
PROOF. See Appendix A. - c , ( q + 1) . . . . . - c t ( - r l ) ] v.
(3.14b)
Then (3.13) can be more compactly written as PO=Q,
(3.15)
where P = YL(I)I/-2
From this theorem and foregoing discussion, the following descriptions are straightforward. Firstly, combining Theorem 3 with Theorem 1, and along the line of above Scheme 1, the linear equation systems with the same form as (3.7) can be obtained'
(3.16a)
and
AO=B,
(3.22)
where, however, Q= YR(1)]l=2.
-c2(-rz+&)
(3.16b)
c4(-r2 + r,)
Thus, we can determine O = (pTp)-l ( p r o ) .
(3.17)
Secondly, combining (3.4) and (3.12), we can obtain the following equation in a similar way: UO = V, Signal Processing
(3.18)
A=
''
c2(-r2-1)
• " c4(-r2
cz( r 2 + l )
1) c 4 ( - r 2 + 1 )
c2(0) c4(0)
c2(--1 +rt)
C2(--2)
C2(0)
c4(-1 +r,)
c4(--2)
£'4(0)
c4(-- 1 + r2)
c2(q+ 1 + r 0
c=(q)
c2(q + 2)
c2(q+ 1 +r2)
c4(q+ 1 + r 0
c4(q)
c4(q + 2)
c4(q+ 1 +r=)
c2(-r,+l) c4(-&+l)
c2(r2 - r,) c.(r:-r,)
c2(0) c,(0)
• " c2(-r~-l) " c4(-rt-1)
C2(--I +r2)
(3.23a)
205
F.-C. Zheng et al. Cumulant-baseddeconvolution and identification and
and B = - [c2(-r2),
c4(-r2)
.....
c2(-
1 ),
Q = ZR(/)I,=2.
c4(-1), c2(q+ 1), c4(q+ 1) . . . .
Then, corresponding to (3.18), we can also obtain
c2(-rl), c4(-rl)] T.
(3.23b)
Secondly, considering Theorem 3 and the Giannakis-Mendel formula (eq. (29a) in [9]) for the case of 4th-order cumulants (with a slight modification here), q
(3.27b)
UO= V,
(3.28)
where U=
ZL(l
)Ii = 4
(3.29a)
and
q
~4 E c2(n--i) b3= Y2 ~ c4(n-i)bi, i=0
V = ZR(I)It=4.
i=0
n = -q .....
2q,
(3.24)
the following two equations - (3.26) and (3.28), with the same form as (3.15) and (3.18), respectively, can be derived. Now let ct(-r2+rO
c/(--r2--1)
cl(--r2+l)
"'"
ct(0)
el(-1 +r0
c1(--2)
ct(0)
•""
ct(-- 1 + r2)
J~,r I
fo. ,
fo.,
"'"
fo,~2
i
:
Zc(l)= fq,r I ct(q+ 1 +rO
fo,, ct(q)
i ct(O)
c¢(-r~-l)
f~,' ct(q + 2)
:
"'"
fq,'2
• • • ct(q + 1 + r2)
: cl(-rt+l)
: ---
ct(rz-r 0
(3.25a) and ZR(I) = [--c/(--r2) . . . . .
- c , ( q + 1) . . . . .
--Ct(--1), g o . . . . .
gq,
- c t ( - r 0 ] T,
(3.25b)
where f,j and gn are, respectively, given by q
fnj = ~ [c2(n - i)c4(i + j ) - c4(n - i)c2(i + j ) ] i=0
and q
g, = ~', [c4(n- i)c2(i) - c2(n - i)c4(i)]. i=0
Corresponding to (3.15), we have P O = Q,
(3.26)
where
P = Zc( l)ll- 2
(3.27a)
(3.29b)
An interesting point to note is that the results of Theorems 1, 2 and 3 can actually be expressed as the following general form:
1 ~ Ojcz(i+j)=Ibl-', Ylj=~, [0,
ifi~[0, q], other.
(3.30)
Here we have assumed Yt ¢ 0, l = 2, 3, 4. A logical question to ask here is: can the above formula (3.3.30) be generalised to the case of any lth-order cumulants (i.e., l = 2, 3, 4, 5, 6 , . . . ) ? Our computer simulations (i.e., empirically or numerically) have given an affirmative answer. However, a mathematical proof for the case of l/> 5 still remains to be found (the proofs for the cases of l = 2, 3 and 4 are presented in Appendix A). On the other hand, it should be pointed out that the above question has only theoretical significance since there has been no report on an application where the 5th- or still higher-order cumulants have been used. Notice that the equations involved with the 4thorder cumulants are of exactly the same form as their counterparts using the 3rd-order cumulants. As a result, the 4th-order cumulant based algorithms can be implemented almost as easily as the 3rd-order cumulant based ones. Furthermore, it can be found that all the families of linear equations we have obtained above are of a similar elegant form, which makes the programming of the algorithms very regular. Regarding other specific aspects of our algorithms, we have the following remarks. Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
206
REMARK 1. Since all the above six families of equations (i.e., (3.7), (3.15), (3.18), (3.22), (3.26) and (3.28)) are linear, the inverse filter coefficients {Ojlj=rl . . . . . 0 . . . . . r2} can normally be determined uniquely. Notice that AVA, pVp and UTU are positive definite, and therefore nonsingular. In the above equations, c/(" ) will be replaced in a practical implementation by its estimate ~/(. ) (as in (4.2) and (4.3) in Section 4). It can be proved that the estimator ~t(" ) is both unbiased and consistent (provided the noise effects are removed) under the following two conditions [9]: (i) the input w(k) is stationary and zero-mean, and its first six cumulants are absolutely summable; and (ii) the system is stable. Obviously, these conditions are very loose and thus can usually be satisfied. Considering that the 0is are measurable functions of ~z(. ) in the above sets of linear equations, the 0~s themselves must be consistent. This is because the measurable function of consistent estimator(s) is also consistent (for details, see [9] and the references therein).
REMARK 2. All the above algorithms are derived from the MA (or FIR) systems. In the case of A R M A (including AR, similarly hereinafter), there are in principle two schemes to apply the above formulae: (i) Since any stable LTI system can be approximated by a truncated MA model, we can use the truncated impulse response his of A R M A systems as the parameters of the approximating MA models, an idea used in [7, 26]. The advantage of this method is that the order of the AR part of A R M A model is not needed. But on the other hand, this method requires a relatively higher accuracy of cumulant estimation. (ii) As is well known, for an ARMA(p, q) process P
q
x(k)+ ~ a i x ( k - i ) = ~ biw(k-i), i~l
p>~q,
i~O
(3.31) the AR coefficients ais can be determined from the equations below: Signal Processing
p
c3(-n)+ ~ aic3(-n+i)=O,
n>q,
(3.32)
i=l
where the AR order p can be selected by using other algorithms (for example, Hankel matrix based methods) (see [8, 9, 18] for details). Then, the residual process p
y(k) = x(k) + ~ aix(k- i),
(3.33)
i--1
where ~i denotes the estimate of ai, reduces to an MA process: q
y(k) = Z blw(k-i),
(3.34)
i-O
to which the MA based algorithms can be applied.
REMARK 3. Although the algorithms proposed above are derived from the viewpoint of deconvolution, they may conveniently be employed for system identification. As soon as the 0is are determined, the MA parameters bis (for MA systems), or the truncated impulse response his (for A R M A systems), can be directly calculated up to a gain constant by using one of the following two schemes. TH1 scheme. Directly employ Theorem 1 (eq. (3.3)). This scheme is very suitable for the case of low additive noise. When the level of additive noise is very high, the obtained bis or his may be rather poor estimates since the estimation of c2(0) can be severely degraded. INV scheme. Invert the obtained S - l = {0~}. This scheme can generate satisfactory results no matter whether the level of the additive noise is low or high, as will be demonstrated later in the simulation section. The reason why Theorem 2 or 3 is not adopted here is that the estimation of higher-order cumulants (especially the 4th-order cumulants) normally has a relatively higher variance, thus direct employment of Theorem 2 or 3 could deteriorate the estimation accuracy ofbis or his. In addition, Theorem 3 cannot give the signs of bis (or his).
F.-C Zheng et al. / Cumulant-based deconvolution and identification
Obviously, when our algorithms are used to identify the MA parameters of systems, the problem of multimodality no longer exists, and the uniqueness of identification is guaranteed. To this extent, the paper provides a feasible solution to the nonlinearity problem which exists with GiannakisMendel (GM) formula [9]: in [9], the nonlinear terms are only concealed by using overparameterisation, but in our derivation above, they are removed by introducing the three theorems and adopting an AR model. In the meantime, our algorithms obviously reflect a fresh channel of application of the GM formula to identification and deconvolution. 4. Since both the 3rd-order cumulants and autocorrelation functions are used, the above algorithms can deal with both non-minimum and minimum phase systems for deconvolution and system identification. The only constraint is that the system to be deconvolved may not have any zeros on the unit circle. This can be considered as the penalty paid for the linearity of the equations. REMARK
R E M A R K 5. Selection of the inverse filter order rl
and r2. The algorithms proposed above are insensitive to the choice of rl and r2, as long as they are taken to be large enough. However, the above algorithms only involve the general case of mixed phase systems. For the special cases of minimum and maximum phase systems, several minor changes must be made to ensure that the above algorithms are numerically more reliable. (i) For minimum phase systems rl should be set to be zero. Although this is not necessary theoretically, rl <0 does sometimes cause large errors, especially in the presence of a high level of additive noise. (ii) Similarly, for maximum phase systems, r2 should be set to be -1. Then, 0_~, instead of 00, should be assumed to be 1. The corresponding changes also need to be made for all families of linear equations ((3.7), (3.15), (3.18), (3.22), (3.26) and (3.28)).
REMARK
207
6. Determination of MA order q. This
problem depends on the phase property of the considered system. (i) For minimum phase systems, the above algorithms are robust to the choice of q value. The reason for this is that the matrix structures in all above sets of equations become irrelevant with q value when systems are minimum phase. (ii) For nonminimum phase systems, the above algorithms are generally sensitive to the selection of MA order q, since the value of q in this case does directly affect the actual contents of the matrices in all sets of equations. Thus, reliable order determination is required. In fact, several higher-order cumulant based methods for order selection have emerged in recent years [4, 8-10, 13, 15, 26]. But in this paper we do not discuss this problem, but assume that the system order q is known. We will also demonstrate these two points by using simulation results in Section 4. R E M A R K 7. A theoretical analysis of the performances of above six families of linear equations ((3.7), (3.15), (3.18), (3.22), (3.26) and (3.28)) is not available at this time. Heuristically, for the case of 3rd-order cumulants, (3.15) and (3.18) are more stable than (3.7), especially when the method (i) mentioned previously in Remark 2 is used for ARMA systems. However, if it is noted that (3.18) involves more 3rd-order cumulant estimation, which is robust to additive Gaussian noise, it is recommended that (3.18) be used as long as a sufficiently long output series can be obtained to give an acceptable 3rd-order cumulant estimation; certainly, (3.15) can also be tried if the SNR is not too low. If SNR is too low, (3.15) can give less accurate results since it involves more autocorrelation terms. Similar comments can also be made for the case of 4th-order cumulants. In a subsequent section of this paper where we consider simulations, we will empirically confirm these points. R E M A R K 8. All the above algorithms can be eas-
ily made recursive. Our work, which is reported in [25, 27], has shown that the convergence rate of Vol. 30, No. 2, January 1993
208
F-C. Zheng et al. / Cumulant-based deconoolution and identification
the recursive versions of the algorithms is very satisfactory.
4. Simulation examples In this section, we apply the algorithms described above to MA (FIR) system deconvolution. Both nonminimum and minimum phase systems are considered. To study the case of 3rd-order cumulants, an exponentially distributed series is used as the skewed input series w(k). For the case of 4th-order cumulants, a uniformly distributed continuous series and an equally distributed discrete series are used as the unskewed input series w(k). Additionally, as in the above algorithm derivation, w(k) is assumed to be IID in all examples below. The output series x(k) is generated by convolving w(k) with the system impulse response. The additive white noise n(k) to x(k) is assumed to be Gaussian, and the signal-to-noise ratio (SNR) in this paper is defined as SNR = 10 l o g , 0 ( ~ / . \E[n (k)]]
(4.1)
In the simulations below, the lth order cumulants c;(m) in the foregoing equations is substituted with their corresponding estimates d;(m). In order to obtain a more accurate cumulant estimation, the well known segment-average scheme [9, 16] was adopted. For convenience, we briefly describe this scheme as follows. First, segment K output samples into L records of N samples each, viz., K = NL. Then, estimate the cumulants ~i)(m) for the ith record by .,(/) - - 1 n~ ct ( m ) - ~ ~ x(i)(n)[xU)(n + m)];- 1,
(4.2)
~ n=nl
where x(°(n) denotes the samples in the ith record. In addition, here, i = 1. . . . . L, n;= max(0, - m ) and nu = m i n ( N - 1, N - m - 1). Finally, average ~°(m) over all records:
et(m) Signal Processing
=
I L ;=I ~°(m).
(4.3)
In this section, we always take K = N x L = 300 x 30, and perform 20 Monte-Carlo (MC) runs for each example. It should be noted that the above segment-average procedure is not essential for the ergodic series. The more detailed analysis regarding these points is contained in [25, 27] where a recursive version of the algorithms is derived. The 'True' curves in the following diagrams are determined by directly inverting the corresponding system transfer function H(z), i.e., directly expanding 1/H(z). Also, in the following simulations, only the MA systems with known order are considered. In the case of A R M A systems, we can first use the methods mentioned in Remark 2 to transform A R M A problems to MA problems. Then, the application of the MA identification and deconvolution algorithms to the obtained MA models becomes straightforward.
4.1. Skewed w(k): the case of 3rd-order cumulants Here, w(k) is taken to be exponentially distributed (y2 = 1, Y3= 2).
E X A M P L E 1. Mixed phase MA system: the model is x'(k) = w(k) - 2.4w(k - 1) + 0.8w(k- 2) + n(k),
(4.4)
with two zeros: zl = 2 and z2 = 0.4. Set r 1 = - 1 5 and r2=15. For the noise levels which produce SNR = 50 dB and SNR = 10 dB, the first MC run is then implemented: (3.7), (3.15) and (3.18) are solved, respectively. The obtained inverse filter coefficients 0is are shown in Figs. (la), (lb) and (lc), respectively. It can be seen that, for this MA system, (i) the results of three equations are of approximately identical accuracy at the lower noise level (SNR = 50 dB), and
209
IE-C. Zheng et al. / Cumulant-based deconvolution and identification 3.0
(a)
I
I
I
I
3.0
2.5-
2.5-
2.0-
2.0-
1.5-
1.5-
1.0-
1.0-
0.5-
0.5-
0.0
0.0
-0.5
I -10
-15
]5
-
I 0
j
I 5
I 10
-0.5 15
Fig. la. Results of (3.7) in Example 1: mixed phase MA(2) system. I
3.0
I
I
2.52.01.51.00.50.0 i
-15
-10
-
0
j
I
I
I
= -10
-
15
0
j
I 5
I 10
15
Fig. lc. Results o f (3.18) in E x a m p l e 1 : m i x e d phase M A ( 2 ) system.
I
(b)
-0.5
-15
I
(c)
i
i
5
10
15
Fig. lb. Results of (3.15) in Example 1: mixed phase MA(2) system.
(ii) the result of (3.15) is of a relatively higher error in comparison with (3.7) and (3.18) for the higher noise level (SNR = 10 dB). The reason for the latter is that, as mentioned in Remark 7, (3.15) involves relatively more autocorrelation (2nd-order cumulant) terms, among which c2(0) is subject to the effect of the additive Gaussian noise. Nineteen other MC runs are similarly accomplished. Only the general identification results (mean+variance, 20 MC runs, similarly hereinafter) are listed in Table 1. Both the two MA identification schemes (TH1 and INV) mentioned in
Remark 3 are used here. Also, for comparison, the g a i n a ~ u s t m e n t has been done for every MC run: be = bi/bo, where b;s are results given by TH1 and INV schemes. For the sake of space, we did not include b0 in Table 1 since 60 = 1.00000 ± 0.00000 (=b0). The same operations will also be used in all the following examples. From Table 1, it is clear that the TH1 scheme can generate very accurate results only when SNR is not too poor, but the INV scheme can satisfactorily secure the accuracy of the results not only in the high SNR case but also in the low SNR case. The latter actually also indicates that the output, 0ts, of our deconvolution algorithms is very robust. As will be seen, similar observation can be obtained in all of the following examples. In addition, the sensitivity of the algorithms to the order mismatch is tested. Figure (ld)illustrates the results of (3.7), (3.15) and (3.18) when the system order is set to be 3 and SNR = 50 dB. It can be seen that the algorithms are sensitive to the selected system order in the case of N M P systems, as explained in Remark 6. Finally, for comparison, the method suggested by Nikias and Chiang in [4, 14] (we refer to this method .as Nikias-Chiang (N-C) method in this paper) is implemented using the same data. For this example, the number of equations in the N - C Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
210 Table 1
Identification results (mean :t: variance) in Example 1
SNR
Identification scheme
50dB
THI
INV
10dB
THI
INV
Equations
/~]
b2
(3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18 )
-2.405864-0.11309 -2.389334-0.12451 -2.398474-0.12663 -2.441144-0.11344 -2.45467 4- 0.14289 -2.440324-0.12528 -4.327594-0.80167 -4.11011+0.62873 -4.28709 4- 0.79086 -2.168844-0.10448 -2.09096 4- 0.13401 -2.09915 4- 0.12463 -2.40000
0.82076:t:0.05539 0.820404-0.06225 0.815354-0.05617 0.820194-0.04714 0.82767 4- 0.06580 0.81559:t:0.04452 0.409414-0.21913 0.284824-0.30937 0.23018 4- 0.37202 0.773874-0.06390 0.68703 4- 0.07342 0.74242 4- 0.06348 0.80000
True bi
,:3.0
J
I
I
I
I
/
--
i
........ Eq(3.7) .... Eq(3.15)
1~
2.0-
~
1.00.5-
I
I
2.0-
...... Eq(a.IR)
1.51.00.50.0
0.0 i -10
15 -
i 0
j
I 5
i 10
-0.5 -15
15
Fig. ld. Sensitivity o f the algorithms to order mismatch: mixed phase M A ( 2 ) system, algorithm o r d e r = 3.
method is ( r 2 - r O ( q + l ) = 9 0 . Figure l(e) is the inverse filter {0t} generated by the N - C method in the first MC run (SNR = 50 dB and 10 dB). Nineteen other MC runs of this method are then performed. The INV scheme is similarly used in all these MC runs, and the obtained MA parameters (corresponding to Table 1) are listed in Table la. Clearly, the results o f our algorithms are more accurate than those of the N - C method, while our algorithmic forms are much simpler. E X A M P L E 2. Minimum phase MA system: the model is Signal Processing
J
2.5-
True
2
1.5-
-15
I
(e)
2.5-
-0.5
1
3.0
(d)
-10
-
J5
i
0
j
i
5
i
10
15
Fig. le. Results of Nikias-Chiang method in Example 1 : mixed phase MA(2) system.
x'(k) = w(k) + 0 . 3 w ( k - 1) - 0 . 4 w ( k - 2) + n(k),
(4.5)
with two zeros: z] = - 0 . 8 and z2=0.5. Table la M A parameters (mean±variance) obtained from the N - C method in Example 1 SNR
/~1
62
50 dB 10 dB True bi
-3.07428 :~ 0.22263 -3.13369 4- 0.26429 -2.40000
1.32381 4- 0.16828 1.37363 4- 0.20516 0.80000
211
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
First, set r~ = 0 and r 2 = 15. Adding the noise which makes S N R = 10 dB and S N R = 1 dB this time, the first MC run is realised. F r o m (3.7), (3.15) and (3.18), the corresponding inverse filters are determined, which are drawn in Figs. (2a), (2b) and (2c), respectively. It is clear that the result of (3.18) is better than that of (3.15) and (3.7) when the additive noise is very strong (1 dB, here), as noted in Remark 7. Similarly, nineteen other MC runs are implemented, and only the corresponding
I
ii
I
I
I
I
I
I
I
I
I
(c)
0.0-
0
(a)
I
~- 0.5--
--0.5
1.5
I
I
I
I
I
2
4
6
8 j
I
I
I
10
12
14
Fig. 2c. Results of (3.18) in Example 2: minimum phase MA(2) system.
10
~- 0.5-
0.0-
-0.5 0
I
I
I
2
4
6
I
8 j
I
t
I
10
12
14
Fig. 2a. Results of (3.7) in Example2: minimum phase MA(2) system•
~- 0.5-
identification results are listed in Table 2. In terms of accuracy, they are comparable with the results listed in Table 1 although the SNR is much lower in this example. To this extent, the algorithms offer a stronger robustness to the additive noise for the MP systems than for the N M P ones. Next, let us carry out two more tests. First, set rl = - 1 5 and r2 = 15. Figure (2d) demonstrates the obtained inverse filters (SNR = 1 dB). As pointed out in Remark 5, serious errors appear in the anticausal parts of the inverse filters. But the result of (3.18) is still more satisfactory than that of (3.7) and (3.15). In fact, during the simulation of other examples, not demonstrated here, we observed much more serious errors. Thus, for minimum systems, it would be better to take rl = 0. Then, the system order is taken to be 3, and the obtained inverse filters from the three algorithms are illustrated in Fig. (2e), where r~ = 0, r2 = 15 and SNR = 10 dB. Clearly, the algorithms are robust to the order mismatch in the case of MP systems, as previously mentioned in Remark 6.
0.0-
-0.5
0
I
I
1
I
2
4
6
8 j
I
I
I
10
12
14
Fig. 2b. Results o f (3.15) in Example 2: minimum phase MA(2) system.
From the above two examples, it can be seen that all three equations, (3.7), (3.15) and (3.18), are robust to the presence of additive noise, all producing satisfactory performances, though there are slight differences in accuracy. Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
212
Table 2 Identification results (mean 4-variance in Example 2 SNR
Identification scheme
10dB
TH1
Equations (3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18) (3.7) (3.15) (3.18)
INV
1 dB
6~
THI
INV
-0.25984-4-0.01515 0.23703+0.01519 0.26651 +0.01572 0.29561 ±0.01864 0.24355±0.01564 0.30816+0.02161 0.12033-4-0.01275 0.09951 ±0.01292 0.15129+0.01455 0.185964-0.01553 0.090824-0.01496 0.31768±0.03964 0.30000
True bi
I
1.5
I
I
I
b2
I
I
1.5
-........ .... ........
1.0-
I
I
I
L%
-0.5
-lo
-'5
1
I Eq(3.15) [q(~.le)
1.0-
/ L
0.5-
0.0-
I
o
]
i
i
5
10
-0.5
15
Fig. 2d. Algorithm comparison: r, is set to be -15, minimum phase MA(2) system, SNR = 1 dB.
4.2. Unskewed and continuous w(k)." the case o f the 4th-order cumulants
o
I
I
I
I
2
4
6
8
j
3. M i x e d p h a s e d M A
s y s t e m : the
I
i
I
10
12
14
Fig. 2e. Insensitivity of the algorithms to order mismatch: minimum phase MA(2) system, algorithm order = 3.
x ' ( k ) = w ( k ) - 2 . 4 w ( k - 1) + 0 . 8 w ( k - 2) + n(k),
T h e d i s t r i b u t i o n o f w ( k ) is t a k e n to be u n i f o r m over the interval [ - 2 , 2] (72 = 4 / 3 , 73 = 0 a n d 74 = - 3 2 / 1 5 ) . F o r c o m p a r i s o n , the same system transfer functions as in Section 4.1 are c o n s i d e r e d in the following examples.
Signal Processing
I
EqO.7)
~-
0.0
EXAMPLE m o d e l is
I
-Troe
True Cq(3.7) Eq(3.15) Eq(3.1e)
0.5-
-15
I
(e)
(d)
~-
-0.35709+0.00916 -0.34979±0.00904 -0.35989±0.00930 -0.38188±0.00976 -0.38219±0.01009 -0.384474-0.00983 -0.194204-0.00584 -0.19016±0.00562 -0.20296±0.00620 -0.28536±0.01466 -0.24499±0.01016 -0.36954±0.02220 -0.40000
(4.6)
with two zeros: 2" 1 = 2 a n d Z 2 = 0.4. Set r~ = - 1 5 a n d r2 = 15. N o i s e which results in S N R = 50 dB a n d S N R = 15 d B is a d d e d to the o u t p u t series, respectively. F i r s t the M C run is implem e n t e d : solving (3.22), we o b t a i n the inverse filters s h o w n in Fig. (3a); next, b y solving (3.26), the inverse filters s h o w n in Fig. (3b) is o b t a i n e d ;
F.-C. Zheng et al. / Cumulant-based deconvolution and identification 3.0
I
(a)
I
I
I
I
3.0 2.5-
2.0-
2.0-
1.5-
1.5-
1.0-
1.0-
0.5-
0.5-
0.0
0.0
I -15
15
-10
Fig. 3a. Results
-
I 0
of (3.22)
j
in Example system.
-0.5
I
I
5
10
3: mixed
15
phase M A ( 2 )
I
I
I
(c)
2.5-
-0.5
213
/\ t
-15
-10
t5 -
i 0
Fig. 3c. R e s u l t s o f (3.28) in E x a m p l e system.
j
i
t
5
10
15
3: m i x e d p h a s e MA(2)
EXAMPLE 4. Minimum phase M A system: the model is I
3.0
I
I
I
I
x'(k) = w(k ) + 0 . 3 w ( k - 1)
(b) 2.5-
- 0 . 4 w ( k - 2) + n(k),
1 r
2.01.51.00.50.0-
-0.5
i -15
-10
15 -
i 0
Fig. 3b. R e s u l t s o f (3.26) in E x a m p l e system.
j
i
i
5
10
15
3 : mixed phase MA(2)
finally, (3.28) is solved, and the determined inverse filters are demonstrated in Fig. (3c). Also, nineteen other MC runs are realised, and the corresponding identification results are given in Table 3. Notice that, in Table 3, the corresponding variance becomes higher than in Table 1. This is due to the fact that estimation of the 4th-order cumulants normally has higher variance than that of the 3rdorder cumulants. Similar explanation also applies to the following Examples 4 and 5.
(4.7)
with two zeros: Zl = - 0 . 8 and z2 = 0.5. Set rl = 0 and r2 = 15. Noise which results in SNR = 50 dB and SNR = 15 dB is added, respectively. The determined inverse filters in the first MC run by (3.22), (3.26) and (3.28) are shown in Figs. (4a), (4b) and (4c), respectively. Similarly, the general identification results for the twenty MC runs are listed in Table 4. Clearly, this table verifies our observation in Example 2 again: the algorithms are of a stronger robustness to the additive noise for the MP systems than for the N M P ones. It can be seen that all the results obtained by the algorithms in the above two examples are satisfactory, but the results from (3.28) are always better than the ones from the other two equations. As to the corresponding errors, their major sources include additive noise, cumulant estimation errors and numerical errors. 4.3. Unskewed and discrete w(k)." the case of the
4th-order cumulants w(k) is taken to be an equally distributed discrete series over set {4-1, +3} (72 = 5, ~ 4 = 26). The additive noise level is set to be 20 dB. Vol. 30, No. 2, January 1993
F.-C. Zheng et al. / Cumulant-based deconvolution and identification
214
Table 3 Identification results (mean ± variance) in Example 3
SNR
Identification scheme
50dB
TH1
INV
15 dB
THI
INV
Equations
b,
/~2
(3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28)
-2.39696±0.26939 -2.365874-0.26746 -2.37549:t:0.27624 -2.471554-0.25671 -2.476325:0.26435 -2.46264±0.25440 -3.00481 5:0.58958 -2.883345:0.52864 -2.93120 ± 0.55689 -2.32376 ± 0.28544 -2.311804-0.31493 -2.31759 4-0.29774 -2.40000
0.82085-4-0.14106 0.813585:0.14183 0.81180±0.14116 0.81589±0.15489 0.80997± 0.16752 0.80254+0.16107 0.724265:0.22767 0.695024-0.24923 0.70533+ 0.23539 0.80739+ 0.17903 0.76410+0.20260 0.78962+ 0.18820 0.80000
True bi
10
I
1.5
1.5 (a) I
I
I
I
I ~
I
I
I
I
I
C
(b)
J
-1.0-
~- 0.5-
"~ 0.5-
0.0-
0.0-
-0.5 0
I
I
I
I
2
4
6
B j
I
I
I
10
12
14
Fig. 4a. Resultsof(3.22) in Example4:minimum phase MA(2) system.
EXAMPLE
5. M i x e d p h a s e d M A
system" the
m o d e l is x ' ( k ) = w ( k ) - 2.4w(k - 1)
+ 0 . 8 w ( k - 2) + n ( k ) ,
(4.8)
with t w o zeros: z, = 2 a n d z2 = 0.4. O b v i o u s l y , this is a typical channel equal±sat±on p r o b l e m . Set r, = - 1 5 a n d r2 = 15. In the presence o f 20 dB additive noise, first M C r u n is c a r r i e d out, a n d the results are as follows. U s i n g (3.22): result Signal P r o c e s s i n g
--0.5
True 50dB 15dB
........ ........
I
I
I
I.
2
4
6
8 j
i
I
I
10
12
14
Fig. 4b. Results of (3.26) in Example 4: minimum phase MA(2) system.
s h o w n in Fig. (5a). U s i n g (3.26): result shown in Fig. (5b). U s i n g (3.28): result shown in Fig. (5c). Then, nineteen similar M C runs are a c c o m p l i s h e d , a n d the c o r r e s p o n d i n g general identification results are listed in T a b l e 5. These five e x a m p l e s s h o w t h a t o u r a l g o r i t h m s can efficiently deal with the p r o b l e m o f d e c o n v o l u tion a n d identification o f b o t h N M P a n d M P systems, while the a l g o r i t h m i c f o r m s are v e r y simple a n d regular. A s to the identification o f M A p a r a m eters, the I N V scheme is very secure n o m a t t e r the noise level is high o r low.
215
F.-C. Zheng et al. / Cumulant-based deconvolution and identification 1.5 ]
I
I
I
I
I
I
,3.0
I
(a)
I
I
I
I
2.51.01
2.0-
~
1.51.0
oolV -0.51 0
0.5 0.0'
I
i
I
I
2
4
6
8
j
i
T
i
10
12
14
-0.5
Fig. 4c. Results of(3.28) in Example4': minimum phase MA(2) system.
15
i
-15
-10
-
i
0
j
i
i
5
10
15
Fig. 5a, Results of (3.22) in Example 5: mixed phase MA(2) system.
Table 4 Identification results (mean 4-variance) in Example 4 SNR
Identification scheme
50 dB
TH1 INV
15dB
THI INV
True bi
Equations (3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28) (3.22) (3.26) (3.28)
6,
62
-0.295604-0.01701 0.29571 +0.01850 0.29571 4-0.01697 0.300634-0.02037 0.30299-4-0.02157 0.29970+0.02161 0.281654-0.01716 0.27771 4-0.01862 0.28433+0.01731 0.294924-0.02049 0.288604-0.02118 0.29824 4-0.02273 0.30000
5. Conclusions Through the above discussion, six families of linear equations have been derived according to three theorems, and they are all of a similar simple form. As a result of the linearity of these equations, the uniqueness of solutions can normally be guaranteed. Thus, the risk of local minima, which exists with some other methods, is removed. In comparison with the existing linear (i.e., noncausal A R based) approaches, the algorithms
-0.404214-0.01607 -0.40408i 0.01639 -0.404414-0.01603 -0.396384-0.01670 -0.402534-0.01610 -0.393685:0.01814 -0.38793~:0.01652 -0.38586+0.01693 -0.389554-0.01652 -0.39428+0.01781 -0.394714-0.01628 -0.396784-0.01984 -0.40000
proposed in this paper possess the following two features: (a) Only the diagonal cumulant slices are used, which in general makes the algorithms more accurate, and also makes the forms of the concerned matrices simple and elegant; (b) 4th-order cumulants are employed in the last three algorithms exactly as the substitute of 3rdorder cumulants in the first three algorithms, which enables the 4th-order cumulant based algorithms (for deconvolution of an unskewed series) to be Vol. 30, No. 2, January 1993
216
F.-C. Zheng et al. / Cumulant-based deconvolution and identification I
3.0
I
I
I
implemented as easily as their 3rd-order cumulants based counterparts (for deconvolution of a skewed series).
I
(b) 2.52,01.5-
Appendix A
i 1.0-
In this appendix, the symbol system used in the text will be followed.
0.50.0 -0.5
I
15
-lO
-15
-
; o
j
A. 1. Proof of Theorem 1
t
i
5
10
15
From the conditions in Theorem 1 and (2.1), we have Fig. 5b. Results of (3.26) in Example 5: mixed phase MA(2) system. 3.0
I
I
I
I
• .,x(i) = E [w(k)x(k + i)]
I
= E [w(k).=~0 b.w(k+i-n)]
(c)
2.5-
(A1.1) (A1.2)
q
= ~ b.E[w(k)w(k+i-n)] 2.0-
(A1.3)
n--0 q
= ~ b.~'26(n- i)
1.5-
(A1.4)
n=0 1.0-
= 0.5-
I~'2bi, if ie[0, q], (0,
0.0
(A1.5)
other.
On the other hand, from (2.2), we have
-0.5
I -15
-I0
15 -
i 0
j
, 5
10
15
• wx(i) = E [ w ( k ) x ( k +
Fig. 5c. Results of (3.28) in Example 5: mixed phase MA(2) sytem.
= E [ j ~ , Ojx(k-j)x(k+i) 1
Table 5 Identification results (mean + variance) in Example 5
SNR
Identification scheme
20dB
TH1
INV
True bi Signal Processing
i)]
Equations
bl
b2
(3.22) (3.26) (3.28) (3.22) (3.26) (3.28)
-2.52986-4-0.40044 -2.48994-4-0.38909 -2.50339-4-0.39954 -2.395934-0.20816 -2.396274-0.22162 -2.38865 4- 0.20661 -2.40000
0.77463+0.11616 0.76475+0.12562 0.765274-0.11811 0.78276+0.11999 0.76706+0.13470 0.77008 4- 0.12863 0.80000
(A1.1) (A1.6)
217
F.-C Zheng et al. / Cumulant-baseddeconvolutionand identification A.3. Proof of Theorem 3
I"2
= ~ 0iEIx(k-j)x(k+i)]
(A1.7)
j=rl r2
= E Ojc2(i+j).
From the conditions in Theorem 3 and (2.1), we have
(A1.8)
j--rl
. . . . (i) = E[w(k)x3(k + i)]
Equating (A1.5) and (A1.8), Theorem 1 holds.
[]
= E w(k)
b~w(k + i - n n
I
A.2. Proof of Theorem 2
[
q b.w(k + i - n
E = w(k
[
(A3.2)
0 q
q
~ ~ bmbnblw(k+i-m)
m=0 n=0 /=0
×w(k+i-n)w(k+i-l)]
(A2.1) (A2.2)
q
q
q
= ~
~
~ bmbnbl
rn=0n=0
n
E:
q
=E w(k) ~
From the conditions in Theorem 2 and (2.1), we have Cbwxx(i)= E[w(k)x2(k + i)]
(A3.1)
(A3.3)
/=0
× E[w(k)w(k + i - m)
q'
×w(k+i-n)w(k+i-l)].
2 Zb b.
(A3.4)
m = O n=O
For the input series w(k), from (3.1c), we have × w(k + i - n)w(k + i - m)] q
=~"
(A2.3) E[w(k) iO, w(k+mi)]
q
Ebmbn
m = 0 n=O
= C~(ml, m2, m3)
x E[w(k)w(k + i - m)w(k + i - n)] (A2.4) q
= ~
+ )'22[~(ml)~(m3-m2)
q
~ bmbnT3~(m-i, n - i )
(A2.5)
+ ~(me)g(m3 - m l ) + 3(m3)g(m2-ml)],
m-0 n=0
~T3b~, if is[0, q], ~0, other.
(A3.5) (A2.6) where C~(ml, m2, m3) is the 4th-order cumulant of w(k). According to [3, 9],
On the other hand, from (2.2), we have
C~(ml,
• wxx(i) = E[w(k)x2(k + i)]
m2,
m3) = ~/4~(ml ,
m2, m 3 ) .
(A3.6)
(A2.1) Then, from (A3.4), we derive
=E[j~Ojx(k-j)x2(k+i)l
(A2.7)
r2
= ~ OjE[x(k-j)xZ(k+i)]
(A2.8)
j=rl r2
= ~ Ojc3(i+j).
. . . . (i) q
= E
q
q
E E b.,b.b,
m=0 n=0 /=0
(A2.9)
j=rl
Equating (A2.6) and (A2.9), Theorem 2 holds.
x [T46(i-m, i - n , i - l )
+ T~{g ( i - m)g(n - l) + g ( i - n)fi(m - 1) []
+ ~ ( i - l)~(m - n) } ]
(A3.7)
Vol. 30, No. 2, January 1993
218
F-C Zheng et al. / Cumulant-based deconvolution and identification 2
/ if i t [0, q],
9,
2
btbi+ E bmbi÷ E bZ,b, ,
)/463÷)'2
m=O
n=O
other (A3.8)
= y4b3+372b, 2 b], n=0 q
t [0,
i f i e I 0 , q], other. (A3.9)
O n the o t h e r h a n d , f r o m (2.2) a n d ( Y l c ) , we have
. . . . (i) = E[w(k)x3(k + i)] =E[
(A3.1)
Ojx(k-j)x3(k+i)]
~
(A3.10)
t-j=rl r2
= Z OjE[x(k-j)x3(k+i)]
(A3.11)
j = rl r2
= ~
Oj[c4(iWj) ÷ 3 c 2 ( i ÷ j ) c 2 ( 0 ) ] (A3.12)
j=rl r2
= ~. Ojc4(i+j)
(A3.13)
j rl r2 +3c2(01 ~ Ojc2(i+j). j--r1
N o t i c i n g T h e o r e m 1 (3.3) a n d q
c2(0)=72 ~ b],
(A3.14)
n=0
we c a n o b t a i n I r2
q
~" O j c 4 ( i + j ) + 3 y 2 b i ~ bZ. n=0
~j=r I
.... (i)=1
rifiE[0'q]'
[ E Ojc4(i+j),
other.
~j=rl
(A3.15) E q u a t i n g (A3.9) a n d (A3.15), we m u s t have r2 Ojc4(i+j) = I r4b3, j=,~ (0, which is j u s t T h e o r e m 3. SignalProcessing
if i t [ 0 , q], other, []
References
(A3.16)
[1] S. Bellini et al., "Blind deconvolution: Polyspectrum or Bussgang techniques?", in: E. Biglieri et al., eds., Digital Communications, Elsevier, Amsterdam, 1986, pp. 251 263. [2] A. Benveniste et al., "Robust identification of a nonminimum phase system: Blind adjustment of a linear equaliser in data communications", IEEE Trans. Automat. Control, Vol. AC-25, No. 3, June 1980, pp. 385 399. [3] D.R. Brillinger and M. Rosenblatt, "Computation and interpretation of kth-order spectra", in: B. Harris, ed., Spectral Analysis of Time Series, Wiley, New York, 1987, pp. 189 232. [4] H.-H. Chiang and C.L. Nikias, "Adaptive deconvolution and identification of nonminimum phase FIR systems based on cumulants", IEEE Trans. Automat. Control, Vol. AC-35, January 1990, pp. 36-47. [5] A.R. Collar and A. Simpson, Matrices and Engineering Dynamics, Ellis Horwood, Chichester, UK, 1987, Chapter 1. [6] B. Friedlander and B. Porat, "Asymptoticallyoptimal estimation of MA and ARMA parameters of non-Gaussian processes from higher-order moments", IEEE Trans. Automat. Control, Vol. AC-35, January 1990, pp. 27 35. [7] G.B. Giannakis, "Cumulants: A powerful tool in signal processing", Proc. IEEE, Vol. 75, No. 9, September 1987, pp. 1333-1334. [8] G. Giannakis and J.M. Mendel, "ARMA order determination via higher-order statistics", in: C.L. Martin, ed., Mathematical Theory of Networks and Systems, NorthHolland, Amsterdam. [9] G.B. Giannakis and J.M. Mendel, "Identification of nonminimum phase system using higher order statistics", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP37, March 1989, pp. 360 377. [10] G.B. Giannakis and J.M. Mendel, "Cumulant-based order determination of non-Gaussian ARMA models", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, August 1990, pp. 1411 1423. [11] S. Kirkpatrick, C.D. Gelalt and M.P. Vecchi Jr., "Optimization by simulated annealing", Science, Vol. 220, No. 4598, 13 May 1983, pp. 671-680. [12] K.S. Lii and M. Rosenblatt, "Deconvolution and estimation of transfer function phase and coefficients for nonGaussian linear process", Ann. Statist., Vol. 10, 1982, pp. 1195 1208. [13] C.L. Nikias, "ARMA bispectrum approach to nonminimum phase system identification", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-36, April 1988, pp. 513-524. [14] C.L. Nikias and H.-H. Chiang, "Higher-order spectrum estimation via noncausal autoregressive modeling and deconvolution", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-36, December 1988, pp. 1911-1913. [15] C.L. Nikias and R. Pan, "ARMA modeling of fourthorder cumulants and phase estimation", Circuit Systems Signal Processing, Vol. 7, No. 3, 1988, pp. 291-325.
F.-C. Zheng et al. / Cumulant-based deconvolution and identification [16] C.L. Nikias et al., "Bispectrum estimation: A digital signal processing framework", IEEE Proc., Vol. 75, No. 7, July 1987, pp. 869-891. [17] W.H. Press et al., Numerical Recipes - The Art of Scientific Computing, Cambridge Univ. Press, Cambridge, 1986, Chapter 10. [18] M.R. Raghuveer and C.L. Nikias, "Bispectrum estimation: A parametric approach", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-33, 1985, pp. 12131230. [19] A. Swami and J.M. Mendel, "Closed-form recursive estimation of MA coefficients using autocorrelations and third-order cumulants", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-37, No. 11, November 1989, pp. 1794-1795. [20] A.M. Tekalp and A.T. Erdem, "Higher-order spectrum factorization in one and two dimensions with applications in signal processing and nonminimum phase system identification", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-37, October 1. ~'~ op. 1537-1549. [21] K.S. Tugnait, "Identification of linear stochastic systems via second- and fourth-order cumulant matching", IEEE Trans. Inform. Theory, Vol. IT-33, May 1987, pp. 393407. [22] A.T. Walden and J.W.J. Hosken, "An investigation of the spectral properties of primary reflection coefficients", Geophysical Prospecting, Vol. 33, 1985, pp. 400-435.
219
[23] A.T. Walden and J.W.J. Hosken, "The nature of the nonGaussianity of primary reflection coefficients and its significance for deconvolution", Geophysical Prospecting, Vol. 34, 1986, pp. 1038-1066. [24] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind deconvolution algorithms based on 3rd- and 4th-order cumulants", Proc. IEEE Internat. Conf. Acoust. Speech Signal Process. (ICASSP), Toronto, Ontario, Canada, 14-17 May 1991, pp. 1573-1576. [25] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "A 2nd- and 4th-order cumulant based blind equalisation algorithm for nonminimum phase channels", Proc. IEEE & EURASIP lnternat. Signal Processing Workshop on Higher Order Statistics, Chamrousse, France, 10-12 July 1991, pp. 79 82. [26] F.-C Zheng, S. McLaughlin and B. Mulgrew, "Blind equalization of nonminimum phase channels: Higherorder cumulant based algorithm", IEEE Trans. Signal Process., to appear. [27] F.-C. Zheng, S. McLaughlin and B. Mulgrew, "Blind equalisation of multilevel PAM data for nonminimum phase channels via 2nd- and 4th-order cumulants", Signal Processing, Vol. 31, No. 3, April 1993, to appear.
Vol. 30, No. 2, January 1993