- Email: [email protected]
HOS-based orthogonal subspace algorithm for causal ARMA system identification
Signal Processing 80 (2000) 535}542
HOS-based orthogonal subspace algorithm for causal ARMA system identi"cation Girish Ganesan!,1,*, K.V.S. Hari" !Department of Systems and Control, Uppsala University, P.O. Box 27, S-751 03, Uppsala, Sweden "Department of ECE, Indian Institute of Science, Bangalore 560012, India Received 15 June 1998; received in revised form 8 October 1999
Abstract In this paper a new method, based on subspaces of a cumulant matrix, is proposed for the blind identi"cation of an ARMA system which is driven by white non-Gaussian noise. The relationship between the cumulants and the impulse response is exploited to arrive at a relation between a cumulant matrix and a matrix consisting of impulse responses. This will lead to the formulation of a new algorithm for impulse response recovery and hence estimation of the ARMA parameters. The proposed method is shown to do well at low SNR values and when impulse response recovery is the key issue, it is shown to o!er reduction in computational complexity. ( 2000 Elsevier Science B.V. All rights reserved. Zusammenfassung In dieser Arbeit wird eine neue Methode, die auf UnterraK umen einer Kumulantenmatrix beruht, zur blinden Identi"kation eines von einem wie{en, nicht-Gau{prozess angeregten ARMA Systems vorgeschlagen. Es wird die Beziehung zwischen den Kumulanten und der Impulsantwort ausgenutzt, um eine entspreehende Beziehung zwischen einer Kumulantenmatrix und einer Matrix bestehend aus des impulsantwort herzuleiten. Dies fuK hrt zu der Formulierung eines neuen Algorithmus zur Bestimmung der Impulsantwort und somit zur SchaK tzung der ARMA Systemparameter. Es wird gezeigt, da{ die vorgeschlagene Methode bei niedrigem SNR gut funktioniert. Au{erdem wird der Rechenaufwand reduziert wenn die Bestimmung der Impulsantwort die entscheidene Rolle spielt. ( 2000 Elsevier Science B.V. All rights reserved. Re2 sume2 Une meH thode nouvelle, baseH e sur les sous-espaces d'une matrice de cumulants, est proposeH e dans cet article pour l'identi"cation aveugle d'un syste`me ARMA exciteH par un bruit blanc non gaussien. La relation entre les cumulants et la reH ponse impulsionnelle est exploiteH e pour arriver a` une relation entre une matrice de cumulants et une matrice constitueH e de reH ponses impulsionnelles. Ceci conduit a` la formulation d'un algorithme nouveau pour la reconstruction de la reH ponse impulsionnelle et, de ce fait, l'estimation des parame`tres ARMA. La meH thode proposeH e est montreH e e( tre performance a` de
1 This work was done when the author was with the Department of Electrical Communication Engineering, Indian Institute of Science, Bangalore. * Corresponding author. Tel.: #46-18-471-7840; fax: #46-18-503-611. E-mail addresses: [email protected] (G. Ganesan), [email protected] (K.V.S. Hari) 0165-1684/00/$ - see front matter ( 2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 9 ) 0 0 1 5 0 - 4
536
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
faibles valeurs de SNR et, lorsque la reconstruction de la reH ponse impulsionnelle est l'objectif cleH , on observe une reH duction de la complexiteH en termes de calculs. ( 2000 Elsevier Science B.V. All rights reserved. Keywords: ARMA system identi"cation; Higher-order statistics; Cumulants; Orthogonal subspace methods
1. Introduction Higher-order statistics (HOS) are particularly useful in the identi"cation of non-minimum phase MA and ARMA systems. In parametric estimation methods for the bispectrum, especially in the estimation of bispectra of sinusoids, HOS based AR modeling is used. An excellent introduction to HOS and its applications can be found in [4,5]. In [7], a linear combination of cumulant slices is used to estimate the impulse response of an ARMA process. From the impulse response values, the ARMA parameters are estimated. In this paper, this idea is further developed to arrive at an orthogonal subspace based framework for estimating the impulse response. A similar approach has been proposed for identifying MA systems in [6], and for identifying AR systems in [2].
2. The ARMA impulse response}cumulant relations Consider an ARMA(p, q) process p q x(n)"! + a(i)x(n!i)# + b(k)e(n!k), a(0)"b(0)"1. i/1 k/0
(1)
The innovation process e(n) is assumed to be zero-mean, stationary, white, non-Gaussian with "nite mth order cumulants c . The impulse response, h(n), satis"es the equation m
G
b(n), 0)n)q, p + a(i)h(n!i)" 0 otherwise. i/0
(2)
Let C (i , i ,2, i ) denote the mth order cumulant of the process x(n). Then according to the m,x 1 2 m~1 Brillinger}Rosenblatt summation formula: = m~1 C (i , i ,2, i )"c + < h(n#i ), i ,0, k 0 m,x 1 2 m~1 m n/~= k/0
m*2.
(3)
For a causal system, this can be written as = C (i, i ,2, i )" + h(n#i)h(n; i ,2, i ), m,x 2 m~1 2 m~1 n/0
(4)
where m~1 h(n; i ,2, i )"c h(n) < h(n#i ). 2 m~1 m k k/2
(5)
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
537
Consider a weighted linear combination of the cumulant slices N N N + w (i , i )C (i, i , i )#2, C (i)"w C (i)# + w (i )C (i, i )# + 3 2 3,x 2 4 2 3 4,x 2 3 w 2 2,x i2 /~M i3 /~M i2 /~M where w , w (i ), etc. are the weights. M and N are arbitrary integers which decide the cumulant slices to be 2 3 2 used in the linear combination. It is possible to write C (i) as cross-correlation of the impulse response and w a causal sequence. That is, = C (i)" + h(n#i)g(n), w n/0 N N N g(n)"w c h(n)# + w (i )h(n; i )# + + w (i , i )h(n; i , i )#2. 2 2 3 2 2 4 2 3 2 3 i2 /~M i2 /~M i3 /~M Now consider the vector . c "[C (!P)2C (0)2C (P)]T 2P`1C1 w w w w Then,
C
C
2,x F
(0)
C
(0, j) 3,x F
C (0, j, k) 4,x F
2
(7)
(8)
DC
C (!P) C (!P, j) 2 C (!P, j, k) 2 2,x 3,x 4,x F F F
c " w
(6)
2
C (P) C (P, j) C (P, j, k) 2 2 2,x 3,x 4,x hggggggggggiggggggggggj S 2P`1C=
w
2 F
D
w (j) "Sw, 3 F
(9)
w (j, k) 4 F hij w =C1
where j, k,2 run from !M to N. w is the weight vector. By considering the right-hand side of (6),(9) can be written as
C
0
0
F
F
0 c "Sw" w h(0) F
h(0) h(1)
DCD
h(0)
h(1)
F
F
2
F
F
2
F
2
F
F
2
F
F
F
2
F
2
g(1)
h(P) h(P#1) 2 h(2P) h(2P#1) 2 hgggggggggigggggggggj H
2P`1C=
g(0)
"Hu.
(10)
g(P)
F hij u
=C1
Let P*max(p, q). From the structure of H it is seen that the "rst P#q#1 columns of H are independent. Using (2) for n'q, any column after the P#q#1st one, is dependent on the previous p#1 columns (since P*p) and hence the rank of H is P#q#1. It follows that Rank(S))P#q#1. Assumption 2.1. The "rst P#q#1 values of h(n) are non-zero.
538
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
Theorem 2.1. Let P*max(p, q), M*P!q and N*q. Then if U is the orthogonal subspace2 of S, U is also 0 0 orthogonal to H. Proof. Given in the appendix. Corollary 2.1. The rank of S is P#q#1. Proof. Given in the appendix.
3. Algorithms In this section a new algorithm is proposed, for determining the ARMA parameters, given the noisy observations of the output. The assumption, that consistent estimates of the cumulants are used, is made. 3.1. Orthogonal subspace method (OS) The VF method is presented in [7]. A new algorithm, the OS algorithm, based on Theorem 2.1 is proposed next. Theorem 3.1. Let U be the orthogonal subspace of S and let h "[h(1), h(2),2, h(2P)]T. Let U (i, j) denote the 0 04 0 i, jth element of U . Further dexne 0
C
U (P#2,1) 0 F
2
U (2P#1,1) 0 F
0 F
U (P#2, P!q) 2 U (2P#1, P!q) 0 0 0 A " F 2 2 04 U (2,1) 2 2 0 F U (2, P!q) 0
2
2
0 F
2
0 F
U (2P,1) 0 F
2 U0 (2P, P!q)
U (2P#1,1) 0 F U (2P#1, P!q) 0
D
,
b "[U (P#1,1),2,U (P#1, P!q),U (P,1),2,U (P, P!q),2,U (1,1),2,U (1, P!q)]T. 04 0 0 0 0 0 0
(11)
Then, A h "!b . 04 04 04
(12)
Proof. The theorem follows immediately from Theorem 2.1, by stacking UTh "0, 1)k)P#1. 0 k It can be shown that A is a full rank matrix [2]. 04 2 The columns of U constitute a basis for the Orthogonal Subspace of S. The size of U is 2P#1]m, where m is the dimension of the 0 0 orthogonal subspace of S.
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
539
3.2. Impulse response recovery and ARMA parameter estimation The orthogonal subspace algorithm is as follows: f the matrix S is estimated from the given data. f The orthogonal subspace U is obtained using the Singular Value Decomposition (SVD). 0 f Using Eq. (12), the impulse response can be estimated. If the exact model order is known, then the ARMA parameters can be estimated [1]. On the other hand if all that is known is an upper bound on the model order (p, q), then this algorithm recovers the "rst 2P#1 values of h(n) where P is the upper-bound on max(p, q). If the Vidal}Fonollosa method is used in such a case, then to recover the "rst 2P#1 values of h(n) the size of the cumulant matrix must be doubled. This increases the computational complexity. But when the exact model order is known and the ARMA parameters are to be estimated, then the Vidal}Fonollosa method is computationally less complex than the OS method. The above method proposed in this paper places a restriction that the "rst P#q#1 values of the impulse response are non zero, where P"max(p, q) for an ARMA(p, q) model. While the method gives better results at low signal-to-noise ratio (SNR) values and is computationally e$cient compared to other methods for impulse response recovery, the authors are not aware of any practical situation where such an assumption holds.
4. Numerical studies Simulation studies are presented to study the performance of the methods. For the simulations the driving white noise is exponentially distributed and is zero-mean adjusted with unit variance and the additive noise is white Gaussian. Only third-order cumulants are considered. One thousand Monte Carlo runs are carried out. The performance measures considered are sample mean, sample variance and the number of outliers. The sample mean and sample variance are calculated after removing the outliers. The method for calculating outliers is as explained below. Let a denote the vector consisting of true parameters. Let a( denote the vector consisting of the estimated i parameters at the ith Monte Carlo run. Then the percentage error incurred, (DDa!a( DD/DDaDD)]100, is cali culated. If it is greater than a tolerance value, the run is considered as an outlier. For this study, a tolerance of 20% has been used. The "gure of 20% was arrived at after considering the performance of di!erent algorithms for the `best casesa which did not produce any outliers. 4.1. Examples Example 4.1. An ARMA(2,2) system, B(z)/A(z) given by A(z)"1!0.8z~1#0.65z~2 and B(z)"1!z~2 is considered. The methods considered are the OS and VF methods. For both methods P"5, M"8 and N"2 are used. The number of data points is 15 000. SNRs"0 and 30 dB. The results are given in Table 1. To make clear the concept of outliers, histograms have been plotted for the estimate of a(1)"!0.8 for 30 dB SNR, in Fig. 1. From the "gure it can be seen that in the case of the OS method some estimates are very much away from the true value. They are termed as outliers and have been ignored when the mean and variance were estimated. The number of the ignored estimates is called the number of outliers. For the 0 dB case the same procedure was followed. But in this case both methods produced outliers.
540
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
Fig. 1. Histograms of the estimate of the parameter a(1)"!0.8 for Example 1 with 30 dB SNR.
Table 1 Sample mean and variance (in brackets) of Example 4.1 for the OS and VF methods. No. of samples used is 15 000. M"8, N"2 Method OS
SNR 0 dB 30 dB
VF
0 dB 30 dB
a(1)"!0.8
a(2)"0.65
!0.8003 (0.0041) !0.7992 (0.0005)
0.6626 (0.0097) 0.6603 (0.0011)
!0.7728 (0.0038) !0.7993 (0.0002)
0.6181 (0.0087) 0.6479 (0.0005)
b(1)"0
b(2)"!1.0
No. of outliers
0.0316 (0.0130) 0.0472 (0.0115)
!0.9326 (0.0136) !0.8974 (0.0089)
757
!0.0608 (0.0145) !0.0822 (0.0015)
!0.9675 (0.0151) !1.0870 (0.0006)
857
302
0
Table 2 Sample mean and variance (in brackets) of Example 4.2 for the OS and VF methods. M"8, N"2. No. of data points used is 15 000. SNR is 30 dB Method
h(1)"0.8000 h(2)"!1.0100 h(3)"!1.3280 h(4)"!0.4059 h(5)"0.5385
Method
OS
VF
0.8573 (0.0126) !0.8828 (0.0076) !1.2693 (0.0071) !0.4277 (0.0046) 0.4751 (0.0029)
0.8029 (0.0158) !1.0059 (0.0064) !1.3252 (0.0025) !0.4094 (0.0056) 0.5343 (0.0034)
h(6)"0.6946 h(7)"0.2057 h(8)"!0.2870 h(9)"!0.3633 h(10)"!0.1041
OS
VF
0.6741 (0.0030) 0.2289 (0.0022) !0.2619 (0.0022) !0.3539 (0.0031) !0.1085 (0.0043)
0.6948 (0.0023) 0.2077 (0.0030) !0.2869 (0.0027) !0.3652 (0.0025) !0.1041 (0.0027)
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
541
Table 3 Sample mean and variance (in brackets) of Example 4.2 for the OS and VF methods. M"8, N"2. No. of data points used is 15 000. SNR is 0 dB Method
h(1)"0.8000 h(2)"!1.0100 h(3)"!1.3280 h(4)"!0.4059 h(5)"0.5385
Method
OS
VF
0.8189 (0.0126) !0.9119 (0.0110) !1.2873 (0.0095) !0.3681 (0.0109) 0.5323 (0.0059)
0.7528 (0.0133) !1.0040 (0.0110) !1.2484 (0.0075) !0.3428 (0.0075) 0.4964 (0.0179)
h(6)"0.6946 h(7)"0.2057 h(8)"!0.2870 h(9)"!0.3633 h(10)"!0.1041
OS
VF
0.7165 (0.0097) 0.2232 (0.0084) !0.2820 (0.0113) !0.3853 (0.0111) !0.1267 (0.0183)
0.6921 (0.0131) 0.1818 (0.0172) !0.2755 (0.0167) !0.3495 (0.0064) !0.1180 (0.0118)
Example 4.2. The same system as in the above example is used for estimating the impulse response h(1) to h(10). Assume that the exact model order is not known, but it is known that max(p, q))2. For the VF method P"10 and for the OS method P"5. For both methods, M"8 and N"2. The number of data points used is 15 000. SNRs"0 and 30 dB. The results are given in Tables 2 and 3. For 30 dB The number of outliers are 273 for the VF method and 301 for the OS method. For 0 dB they are 978 and 919, respectively. 5. Discussion It should be noted that since the sample mean and sample variance have been calculated after the removal of outliers, if the number of outliers is quite high, the sample mean and sample variance will not make much sense as the number of runs considered in calculating the statistics becomes quite small. From the tables it is observed that, in the case of parameter estimation: In the noise free case the number of outliers is high for the OS method. But at 0 dB OS method performs better than the VF method. But when impulse response recovery is the key issue, the OS method gives comparable performance with the VF method. But the computational complexity in computing the cumulant matrix is far less for the OS method.
Appendix A Proof of Theorem 2.1. Choose w such that g(n)"c(n)b(q#n!k), 0)k)P#q#1. This can be done by choosing w (i )"a(!i #q!k) and setting all other weights to zero. Then, 3 2 2 N q~k g(n)" + w (i )h(n;i )" + a(!i #q!k)c h(n)h(n#i ). 3 2 2 2 3 2 i2 /~M i2 /~p`q~k Let l"!i #q!k. Then, from (2), 2 p g(n)"c h(n) + a(l)h(n#q!k!l)"c h(n)b(n#q!k). 3 3 l/0
542
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
Now c(n)"c h(n). b(n#q!k)"0 for n'k. Let h denote the kth column of H. The claim now is if U is 3 k 0 orthogonal to h , 0)j(k, where 0(k)P#q#1, and c(k)O0, then U is orthogonal to h . To prove j 0 k this set g(n)"c(n)b(n#q!k). Then Sw"Hu"c(0)b(q!k)h #c(1)b(q#1!k)h #2#c(k!1)b(q!1)h #c(k)b(q)h . (13) 0 1 k~1 k Since UTSw"0 and U is orthogonal to h for 0)j(k, this implies 0 0 j c(k)b(q)UT h "0. (14) 0 k If h(k)O0, then c(k)O0. Hence from Assumption 2.1 it follows that c(k)O0 for 0)k)P#q. Thus UT h "0. Setting k"0, in Eq. (13), it is easily seen that, UTh "0. It follows by induction that U is 0 0 0 0 k orthogonal to the "rst P#q#1 columns of H, which form a basis for the column space of H. This proves Theorem 2.1. In the proof it was assumed that third-order cumulants are used in the formation of the S matrix. This was only done to simplify the proof so as to make it more clear. This can be done for any mth order cumulants by suitably adjusting the non-zero values of w. Proof of Corollary 2.1. Note that the rank of S)P#q#1. Also any vector orthogonal to S is orthogonal to h , 0)k)P#q#1. If the rank of S is less than P#q#1 then there will be more than P!q vectors k orthogonal to S and hence to h . This is not possible since the rank of H is P#q#1. Thus it follows that the k rank of S"P#q#1.
References [1] J.A.R. Fonollosa, J. Vidal, System identi"cation using a linear combination of cumulant slices, IEEE Trans. Signal Processing SP-41 (1993) 2405}2412. [2] G. Ganesan, HOS based orthogonal subspace algorithms for causal IIR system identi"cation, Master's Thesis, Indian Institute of Science, Bangalore, India, June 1998. [4] J.M. Mendel, Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications, Proc. IEEE 79 (1991) 278}305. [5] C.L. Nikias, M. Raghuveer, Bispectrum estimation: a digital signal processing framework, Proc. IEEE 72 (1987) 869}891. [6] L. Srinivas, K.V.S. Hari, FIR system identi"cation based on subspaces of a higher order cumulant matrix, IEEE Trans. Signal Processing SP-44 (1996) 1485}1491. [7] J. Vidal, J.A.R. Fonollosa, Impulse response recovery of linear systems through weighted cumulant slices, IEEE Trans. Signal Processing SP-44 (1996) 2626}2632.
HOS-based orthogonal subspace algorithm for causal ARMA system identi"cation Girish Ganesan!,1,*, K.V.S. Hari" !Department of Systems and Control, Uppsala University, P.O. Box 27, S-751 03, Uppsala, Sweden "Department of ECE, Indian Institute of Science, Bangalore 560012, India Received 15 June 1998; received in revised form 8 October 1999
Abstract In this paper a new method, based on subspaces of a cumulant matrix, is proposed for the blind identi"cation of an ARMA system which is driven by white non-Gaussian noise. The relationship between the cumulants and the impulse response is exploited to arrive at a relation between a cumulant matrix and a matrix consisting of impulse responses. This will lead to the formulation of a new algorithm for impulse response recovery and hence estimation of the ARMA parameters. The proposed method is shown to do well at low SNR values and when impulse response recovery is the key issue, it is shown to o!er reduction in computational complexity. ( 2000 Elsevier Science B.V. All rights reserved. Zusammenfassung In dieser Arbeit wird eine neue Methode, die auf UnterraK umen einer Kumulantenmatrix beruht, zur blinden Identi"kation eines von einem wie{en, nicht-Gau{prozess angeregten ARMA Systems vorgeschlagen. Es wird die Beziehung zwischen den Kumulanten und der Impulsantwort ausgenutzt, um eine entspreehende Beziehung zwischen einer Kumulantenmatrix und einer Matrix bestehend aus des impulsantwort herzuleiten. Dies fuK hrt zu der Formulierung eines neuen Algorithmus zur Bestimmung der Impulsantwort und somit zur SchaK tzung der ARMA Systemparameter. Es wird gezeigt, da{ die vorgeschlagene Methode bei niedrigem SNR gut funktioniert. Au{erdem wird der Rechenaufwand reduziert wenn die Bestimmung der Impulsantwort die entscheidene Rolle spielt. ( 2000 Elsevier Science B.V. All rights reserved. Re2 sume2 Une meH thode nouvelle, baseH e sur les sous-espaces d'une matrice de cumulants, est proposeH e dans cet article pour l'identi"cation aveugle d'un syste`me ARMA exciteH par un bruit blanc non gaussien. La relation entre les cumulants et la reH ponse impulsionnelle est exploiteH e pour arriver a` une relation entre une matrice de cumulants et une matrice constitueH e de reH ponses impulsionnelles. Ceci conduit a` la formulation d'un algorithme nouveau pour la reconstruction de la reH ponse impulsionnelle et, de ce fait, l'estimation des parame`tres ARMA. La meH thode proposeH e est montreH e e( tre performance a` de
1 This work was done when the author was with the Department of Electrical Communication Engineering, Indian Institute of Science, Bangalore. * Corresponding author. Tel.: #46-18-471-7840; fax: #46-18-503-611. E-mail addresses: [email protected] (G. Ganesan), [email protected] (K.V.S. Hari) 0165-1684/00/$ - see front matter ( 2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 9 ) 0 0 1 5 0 - 4
536
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
faibles valeurs de SNR et, lorsque la reconstruction de la reH ponse impulsionnelle est l'objectif cleH , on observe une reH duction de la complexiteH en termes de calculs. ( 2000 Elsevier Science B.V. All rights reserved. Keywords: ARMA system identi"cation; Higher-order statistics; Cumulants; Orthogonal subspace methods
1. Introduction Higher-order statistics (HOS) are particularly useful in the identi"cation of non-minimum phase MA and ARMA systems. In parametric estimation methods for the bispectrum, especially in the estimation of bispectra of sinusoids, HOS based AR modeling is used. An excellent introduction to HOS and its applications can be found in [4,5]. In [7], a linear combination of cumulant slices is used to estimate the impulse response of an ARMA process. From the impulse response values, the ARMA parameters are estimated. In this paper, this idea is further developed to arrive at an orthogonal subspace based framework for estimating the impulse response. A similar approach has been proposed for identifying MA systems in [6], and for identifying AR systems in [2].
2. The ARMA impulse response}cumulant relations Consider an ARMA(p, q) process p q x(n)"! + a(i)x(n!i)# + b(k)e(n!k), a(0)"b(0)"1. i/1 k/0
(1)
The innovation process e(n) is assumed to be zero-mean, stationary, white, non-Gaussian with "nite mth order cumulants c . The impulse response, h(n), satis"es the equation m
G
b(n), 0)n)q, p + a(i)h(n!i)" 0 otherwise. i/0
(2)
Let C (i , i ,2, i ) denote the mth order cumulant of the process x(n). Then according to the m,x 1 2 m~1 Brillinger}Rosenblatt summation formula: = m~1 C (i , i ,2, i )"c + < h(n#i ), i ,0, k 0 m,x 1 2 m~1 m n/~= k/0
m*2.
(3)
For a causal system, this can be written as = C (i, i ,2, i )" + h(n#i)h(n; i ,2, i ), m,x 2 m~1 2 m~1 n/0
(4)
where m~1 h(n; i ,2, i )"c h(n) < h(n#i ). 2 m~1 m k k/2
(5)
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
537
Consider a weighted linear combination of the cumulant slices N N N + w (i , i )C (i, i , i )#2, C (i)"w C (i)# + w (i )C (i, i )# + 3 2 3,x 2 4 2 3 4,x 2 3 w 2 2,x i2 /~M i3 /~M i2 /~M where w , w (i ), etc. are the weights. M and N are arbitrary integers which decide the cumulant slices to be 2 3 2 used in the linear combination. It is possible to write C (i) as cross-correlation of the impulse response and w a causal sequence. That is, = C (i)" + h(n#i)g(n), w n/0 N N N g(n)"w c h(n)# + w (i )h(n; i )# + + w (i , i )h(n; i , i )#2. 2 2 3 2 2 4 2 3 2 3 i2 /~M i2 /~M i3 /~M Now consider the vector . c "[C (!P)2C (0)2C (P)]T 2P`1C1 w w w w Then,
C
C
2,x F
(0)
C
(0, j) 3,x F
C (0, j, k) 4,x F
2
(7)
(8)
DC
C (!P) C (!P, j) 2 C (!P, j, k) 2 2,x 3,x 4,x F F F
c " w
(6)
2
C (P) C (P, j) C (P, j, k) 2 2 2,x 3,x 4,x hggggggggggiggggggggggj S 2P`1C=
w
2 F
D
w (j) "Sw, 3 F
(9)
w (j, k) 4 F hij w =C1
where j, k,2 run from !M to N. w is the weight vector. By considering the right-hand side of (6),(9) can be written as
C
0
0
F
F
0 c "Sw" w h(0) F
h(0) h(1)
DCD
h(0)
h(1)
F
F
2
F
F
2
F
2
F
F
2
F
F
F
2
F
2
g(1)
h(P) h(P#1) 2 h(2P) h(2P#1) 2 hgggggggggigggggggggj H
2P`1C=
g(0)
"Hu.
(10)
g(P)
F hij u
=C1
Let P*max(p, q). From the structure of H it is seen that the "rst P#q#1 columns of H are independent. Using (2) for n'q, any column after the P#q#1st one, is dependent on the previous p#1 columns (since P*p) and hence the rank of H is P#q#1. It follows that Rank(S))P#q#1. Assumption 2.1. The "rst P#q#1 values of h(n) are non-zero.
538
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
Theorem 2.1. Let P*max(p, q), M*P!q and N*q. Then if U is the orthogonal subspace2 of S, U is also 0 0 orthogonal to H. Proof. Given in the appendix. Corollary 2.1. The rank of S is P#q#1. Proof. Given in the appendix.
3. Algorithms In this section a new algorithm is proposed, for determining the ARMA parameters, given the noisy observations of the output. The assumption, that consistent estimates of the cumulants are used, is made. 3.1. Orthogonal subspace method (OS) The VF method is presented in [7]. A new algorithm, the OS algorithm, based on Theorem 2.1 is proposed next. Theorem 3.1. Let U be the orthogonal subspace of S and let h "[h(1), h(2),2, h(2P)]T. Let U (i, j) denote the 0 04 0 i, jth element of U . Further dexne 0
C
U (P#2,1) 0 F
2
U (2P#1,1) 0 F
0 F
U (P#2, P!q) 2 U (2P#1, P!q) 0 0 0 A " F 2 2 04 U (2,1) 2 2 0 F U (2, P!q) 0
2
2
0 F
2
0 F
U (2P,1) 0 F
2 U0 (2P, P!q)
U (2P#1,1) 0 F U (2P#1, P!q) 0
D
,
b "[U (P#1,1),2,U (P#1, P!q),U (P,1),2,U (P, P!q),2,U (1,1),2,U (1, P!q)]T. 04 0 0 0 0 0 0
(11)
Then, A h "!b . 04 04 04
(12)
Proof. The theorem follows immediately from Theorem 2.1, by stacking UTh "0, 1)k)P#1. 0 k It can be shown that A is a full rank matrix [2]. 04 2 The columns of U constitute a basis for the Orthogonal Subspace of S. The size of U is 2P#1]m, where m is the dimension of the 0 0 orthogonal subspace of S.
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
539
3.2. Impulse response recovery and ARMA parameter estimation The orthogonal subspace algorithm is as follows: f the matrix S is estimated from the given data. f The orthogonal subspace U is obtained using the Singular Value Decomposition (SVD). 0 f Using Eq. (12), the impulse response can be estimated. If the exact model order is known, then the ARMA parameters can be estimated [1]. On the other hand if all that is known is an upper bound on the model order (p, q), then this algorithm recovers the "rst 2P#1 values of h(n) where P is the upper-bound on max(p, q). If the Vidal}Fonollosa method is used in such a case, then to recover the "rst 2P#1 values of h(n) the size of the cumulant matrix must be doubled. This increases the computational complexity. But when the exact model order is known and the ARMA parameters are to be estimated, then the Vidal}Fonollosa method is computationally less complex than the OS method. The above method proposed in this paper places a restriction that the "rst P#q#1 values of the impulse response are non zero, where P"max(p, q) for an ARMA(p, q) model. While the method gives better results at low signal-to-noise ratio (SNR) values and is computationally e$cient compared to other methods for impulse response recovery, the authors are not aware of any practical situation where such an assumption holds.
4. Numerical studies Simulation studies are presented to study the performance of the methods. For the simulations the driving white noise is exponentially distributed and is zero-mean adjusted with unit variance and the additive noise is white Gaussian. Only third-order cumulants are considered. One thousand Monte Carlo runs are carried out. The performance measures considered are sample mean, sample variance and the number of outliers. The sample mean and sample variance are calculated after removing the outliers. The method for calculating outliers is as explained below. Let a denote the vector consisting of true parameters. Let a( denote the vector consisting of the estimated i parameters at the ith Monte Carlo run. Then the percentage error incurred, (DDa!a( DD/DDaDD)]100, is cali culated. If it is greater than a tolerance value, the run is considered as an outlier. For this study, a tolerance of 20% has been used. The "gure of 20% was arrived at after considering the performance of di!erent algorithms for the `best casesa which did not produce any outliers. 4.1. Examples Example 4.1. An ARMA(2,2) system, B(z)/A(z) given by A(z)"1!0.8z~1#0.65z~2 and B(z)"1!z~2 is considered. The methods considered are the OS and VF methods. For both methods P"5, M"8 and N"2 are used. The number of data points is 15 000. SNRs"0 and 30 dB. The results are given in Table 1. To make clear the concept of outliers, histograms have been plotted for the estimate of a(1)"!0.8 for 30 dB SNR, in Fig. 1. From the "gure it can be seen that in the case of the OS method some estimates are very much away from the true value. They are termed as outliers and have been ignored when the mean and variance were estimated. The number of the ignored estimates is called the number of outliers. For the 0 dB case the same procedure was followed. But in this case both methods produced outliers.
540
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
Fig. 1. Histograms of the estimate of the parameter a(1)"!0.8 for Example 1 with 30 dB SNR.
Table 1 Sample mean and variance (in brackets) of Example 4.1 for the OS and VF methods. No. of samples used is 15 000. M"8, N"2 Method OS
SNR 0 dB 30 dB
VF
0 dB 30 dB
a(1)"!0.8
a(2)"0.65
!0.8003 (0.0041) !0.7992 (0.0005)
0.6626 (0.0097) 0.6603 (0.0011)
!0.7728 (0.0038) !0.7993 (0.0002)
0.6181 (0.0087) 0.6479 (0.0005)
b(1)"0
b(2)"!1.0
No. of outliers
0.0316 (0.0130) 0.0472 (0.0115)
!0.9326 (0.0136) !0.8974 (0.0089)
757
!0.0608 (0.0145) !0.0822 (0.0015)
!0.9675 (0.0151) !1.0870 (0.0006)
857
302
0
Table 2 Sample mean and variance (in brackets) of Example 4.2 for the OS and VF methods. M"8, N"2. No. of data points used is 15 000. SNR is 30 dB Method
h(1)"0.8000 h(2)"!1.0100 h(3)"!1.3280 h(4)"!0.4059 h(5)"0.5385
Method
OS
VF
0.8573 (0.0126) !0.8828 (0.0076) !1.2693 (0.0071) !0.4277 (0.0046) 0.4751 (0.0029)
0.8029 (0.0158) !1.0059 (0.0064) !1.3252 (0.0025) !0.4094 (0.0056) 0.5343 (0.0034)
h(6)"0.6946 h(7)"0.2057 h(8)"!0.2870 h(9)"!0.3633 h(10)"!0.1041
OS
VF
0.6741 (0.0030) 0.2289 (0.0022) !0.2619 (0.0022) !0.3539 (0.0031) !0.1085 (0.0043)
0.6948 (0.0023) 0.2077 (0.0030) !0.2869 (0.0027) !0.3652 (0.0025) !0.1041 (0.0027)
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
541
Table 3 Sample mean and variance (in brackets) of Example 4.2 for the OS and VF methods. M"8, N"2. No. of data points used is 15 000. SNR is 0 dB Method
h(1)"0.8000 h(2)"!1.0100 h(3)"!1.3280 h(4)"!0.4059 h(5)"0.5385
Method
OS
VF
0.8189 (0.0126) !0.9119 (0.0110) !1.2873 (0.0095) !0.3681 (0.0109) 0.5323 (0.0059)
0.7528 (0.0133) !1.0040 (0.0110) !1.2484 (0.0075) !0.3428 (0.0075) 0.4964 (0.0179)
h(6)"0.6946 h(7)"0.2057 h(8)"!0.2870 h(9)"!0.3633 h(10)"!0.1041
OS
VF
0.7165 (0.0097) 0.2232 (0.0084) !0.2820 (0.0113) !0.3853 (0.0111) !0.1267 (0.0183)
0.6921 (0.0131) 0.1818 (0.0172) !0.2755 (0.0167) !0.3495 (0.0064) !0.1180 (0.0118)
Example 4.2. The same system as in the above example is used for estimating the impulse response h(1) to h(10). Assume that the exact model order is not known, but it is known that max(p, q))2. For the VF method P"10 and for the OS method P"5. For both methods, M"8 and N"2. The number of data points used is 15 000. SNRs"0 and 30 dB. The results are given in Tables 2 and 3. For 30 dB The number of outliers are 273 for the VF method and 301 for the OS method. For 0 dB they are 978 and 919, respectively. 5. Discussion It should be noted that since the sample mean and sample variance have been calculated after the removal of outliers, if the number of outliers is quite high, the sample mean and sample variance will not make much sense as the number of runs considered in calculating the statistics becomes quite small. From the tables it is observed that, in the case of parameter estimation: In the noise free case the number of outliers is high for the OS method. But at 0 dB OS method performs better than the VF method. But when impulse response recovery is the key issue, the OS method gives comparable performance with the VF method. But the computational complexity in computing the cumulant matrix is far less for the OS method.
Appendix A Proof of Theorem 2.1. Choose w such that g(n)"c(n)b(q#n!k), 0)k)P#q#1. This can be done by choosing w (i )"a(!i #q!k) and setting all other weights to zero. Then, 3 2 2 N q~k g(n)" + w (i )h(n;i )" + a(!i #q!k)c h(n)h(n#i ). 3 2 2 2 3 2 i2 /~M i2 /~p`q~k Let l"!i #q!k. Then, from (2), 2 p g(n)"c h(n) + a(l)h(n#q!k!l)"c h(n)b(n#q!k). 3 3 l/0
542
G. Ganesan, K.V.S. Hari / Signal Processing 80 (2000) 535}542
Now c(n)"c h(n). b(n#q!k)"0 for n'k. Let h denote the kth column of H. The claim now is if U is 3 k 0 orthogonal to h , 0)j(k, where 0(k)P#q#1, and c(k)O0, then U is orthogonal to h . To prove j 0 k this set g(n)"c(n)b(n#q!k). Then Sw"Hu"c(0)b(q!k)h #c(1)b(q#1!k)h #2#c(k!1)b(q!1)h #c(k)b(q)h . (13) 0 1 k~1 k Since UTSw"0 and U is orthogonal to h for 0)j(k, this implies 0 0 j c(k)b(q)UT h "0. (14) 0 k If h(k)O0, then c(k)O0. Hence from Assumption 2.1 it follows that c(k)O0 for 0)k)P#q. Thus UT h "0. Setting k"0, in Eq. (13), it is easily seen that, UTh "0. It follows by induction that U is 0 0 0 0 k orthogonal to the "rst P#q#1 columns of H, which form a basis for the column space of H. This proves Theorem 2.1. In the proof it was assumed that third-order cumulants are used in the formation of the S matrix. This was only done to simplify the proof so as to make it more clear. This can be done for any mth order cumulants by suitably adjusting the non-zero values of w. Proof of Corollary 2.1. Note that the rank of S)P#q#1. Also any vector orthogonal to S is orthogonal to h , 0)k)P#q#1. If the rank of S is less than P#q#1 then there will be more than P!q vectors k orthogonal to S and hence to h . This is not possible since the rank of H is P#q#1. Thus it follows that the k rank of S"P#q#1.
References [1] J.A.R. Fonollosa, J. Vidal, System identi"cation using a linear combination of cumulant slices, IEEE Trans. Signal Processing SP-41 (1993) 2405}2412. [2] G. Ganesan, HOS based orthogonal subspace algorithms for causal IIR system identi"cation, Master's Thesis, Indian Institute of Science, Bangalore, India, June 1998. [4] J.M. Mendel, Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications, Proc. IEEE 79 (1991) 278}305. [5] C.L. Nikias, M. Raghuveer, Bispectrum estimation: a digital signal processing framework, Proc. IEEE 72 (1987) 869}891. [6] L. Srinivas, K.V.S. Hari, FIR system identi"cation based on subspaces of a higher order cumulant matrix, IEEE Trans. Signal Processing SP-44 (1996) 1485}1491. [7] J. Vidal, J.A.R. Fonollosa, Impulse response recovery of linear systems through weighted cumulant slices, IEEE Trans. Signal Processing SP-44 (1996) 2626}2632.