2. Manuscript

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

On the Relationship of Non-parametric Methods for Coherence Function Estimation Chengshi Zheng, Mingyuan Zhou, and Xiaodong Li Institute of Acoustics, Chinese Academy of Sciences, Beijing, 100190, China

Abstract In this paper, several seemingly disparate non-parametric magnitude squared coherence (MSC) estimation methods, including the Welch’s averaged periodogram, the minimum variance distortionless response (MVDR), and the canonical correlation analysis (CCA) methods, are treated in a unified way, which makes it simpler to understand the methods and their properties. This uncovered relationship also brings out a new class of MSC estimators in terms of non-linear functions of the covariance matrix. Key words: Non-parametric; Magnitude squared coherence (MSC); Minimum variance distortionless response (MVDR); Canonical correlation analysis (CCA). PACS: 43.60.Hj

1

Introduction

The power spectrum density (PSD) estimation has a relatively long history and the researchers have made great efforts to improve the estimation accuracy [1-7]. Yet, the magnitude squared coherence (MSC) estimation has frequently received less attention [8-13]. As a matter of fact, the MSC estimation also has been widely used in various fields [8, 13], such as speech enhancement, speech recognition, signal detection and estimation, and communications. In recent years, the MSC estimation gradually becomes an attractive research topic [10, 11]. The same as the PSD estimation, there are basically two categories of techniques for the MSC estimation. One is the non-parametric approach, including the Welch’s averaged periodogram method (will be referred as Welch-MSC)[9], the minimum variance distortionless response method (will be referred as MVDR-MSC) [10], and the canonical correlation analysis method (will be Preprint submitted to Elsevier

21 April 2008

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

referred as CCA-MSC) [11]. The other is the parametric method which estimates the MSC by an ARMA model [12, 13]. Compared with the parametric method, the non-parametric approach is more popular for its implementation simplicity and moderate performance. The three non-parametric MSC methods have been proposed independently in the literature. It is well-known that the MVDR-MSC method has higher spectral resolution than the Welch-MSC method, but the property of the CCAMSC method is still unclear. In this paper, to make it simpler to understand the three non-parametric MSC methods, they will be treated in a unified way. This uncovered relationship also brings out a new class of MSC estimators in terms of non-linear functions of the covariance matrix. To be mentioned, Pisarenko also proposed the spectral estimators in terms of non-linear functions of the covariance matrix [6], the intrinsic difference between the Pisarenko’s spectral estimators and the new class of MSC estimators is the direct or indirect usages of non-linear functions of the covariance matrix. This paper is organized as follows. In Section 2, we exploit the relationship of the non-parametric MSC methods, and a new class of MSC estimators are also presented in this section. Conclusions are given in Section 3.

2

Theory

2.1 The non-parametric MSC methods

Let x1 (n) and x2 (n) be two stationary complex time series whose MSC, γx21 x2 (wk ) is to be estimated. Then γx21 x2

|Sx1 x2 (wk )|2 . (wk ) = Sx1 x1 (wk ) Sx2 x2 (wk )

(1)

where Sx1 x1 (wk ) and Sx2 x2 (wk ) are the respective spectra of x1 (n) and x2 (n), Sx1 x2 (wk ) is the cross-spectrum between x1 (n) and x2 (n), and wk = 2πk/K, k = 0, 1, · · · , K − 1. To rewrite (1) in the matrix way, defining the K Fourier vectors 

1 fk = √ 1 e−jwk · · · e−jwk (L−1) L 2

T

,

(2)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where L is the window length parameter. The Welch-MSC is estimated as W elch 2 γ x1 x2

(wk ) =

 2  H  fk Rx1 x2 fk 

(fkH Rx1 x1 fk ) (fkH Rx2 x2 fk )

,

(3)

where superscript H denotes Hermitian conjugate of a vector or a matrix, and Rxi xj are the covariance matrices of xi (n) and xj (n), which are defined as 



Rxi xj = E xi (n)xH j (n) ,

i, j = 1, 2,



(4)

T

xi (n) xi (n − 1) · · · xi (n − L + 1) , and i, j = 1, 2. The MVDR-MSC [10] and the CCA-MSC [11] are calculated as

where xi (n) =

 2  H −1  f fk Rx1 x1 Rx1 x2 R−1 k x2 x2  M V DR 2  , γx1 x2 (wk ) = 

fkH R−1 x1 x1 fk

CCA 2 γ x1 x2

(5)

fkH R−1 x2 x2 fk

 

2 

−1/2 (wk ) = fkH R−1/2 x1 x1 Rx1 x2 Rx2 x2 fk  .

(6)

2.2 The relationship of the non-parametric MSC methods Considering fkH fk = 1 and Rxi xi , i = 1, 2 are Hermitian symmetric matrices. (6) becomes CCA 2 γ x1 x2



(wk ) =

−1/2 Rx x fk 1 1 −1/2 H f Rx x fk k 1 1

H

 −1/2 2  R−1/2 f H Rx x fk   x1 x1 k 2 2 Rx1 x2  H −1/2  −1/2 H f f R f Rx x fk  k x1 x1 k  k 2 2  −1/2 −1/2 H

Rx1 x1

Rx x fk 1 1 −1/2 f H Rx x fk k 1 1

Rx x fk 2 2 −1/2 f H Rx x fk k 2 2



Rx2 x2

−1/2 Rx x fk 2 2 −1/2 H f Rx x fk k 2 2



.

(7)

From (7), we deduce the two filters of the CCA-MSC, CCA

gi,k =

R−1/2 xi xi fk −1/2

fkH Rxi xi fk

,

i = 1, 2.

(8)

While the two filters of the MVDR-MSC and Welch-MSC are, M V DR

W elch

gi,k =

R−1 xi xi fk , H −1 fk Rxixi fk

gi,k = fk =

R−0 xi xi fk , H −0 fk Rxi xi fk 3

i = 1, 2.

i = 1, 2.

(9)

(10)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

From (8) to (10), the only difference is the exponential of the covariance matrices. Thus, we introduce a new class of MSC estimators by the following two filters, R−αxi fk α gi,k = H xi−α , i = 1, 2. (11) fk Rxixi fk and, finally, the MSC is estimated as  2  H −α  f fk Rx1 x1 Rx1 x2 R−α k x2 x2  α 2    . γx1 x2 (wk ) =

fkH R1−2α x1 x1 fk

fkH R1−2α x2 x2 fk

(12)

It is clear that (12) becomes the Welch-MSC method as α = 0, it reduces to the CCA-MSC method as α = 0.5 [11], and it turns into the MVDR-MSC method as α = 1 [10]. Therefore, the three MSC estimation methods, which are seemingly disparate, have intrinsic relation in nature. Eqs. (11) and (12) show that the new class of MSC estimators use the nonlinear functions of the covariance matrix indirectly. In the other hand, the Pisarenko’s spectral estimators use the non-linear functions of the covariance matrix directly as  

1/α 

Sxi xi ,α (wk ) = fkH Rαxi xi fk 

, i, j = 1, 2.

(13)

2.3 Properties of the novel MSC estimators α γx21 x2 (wk ) Property 1: 0 ≤ α γx21 x2 (wk ) ≤ 1 Proof: Since the covariance matrices Rx1 x1 and Rx2 x2 are assumed to be positive definite, it is clear that α γx21 x2 (wk ) ≥ 0. Define the vectors, α

fi,k = R(1−2α)/2 fk , i = 1, 2. xi xi

(14)

Substitute (14) into (12), the MSC is now:  2  α  H −1/2 α ( f1,k ) R−1/2  R R ( f ) x x 2,k x 1 2 x x x 1 1 2 2 α 2   . γx1 x2 (wk ) = H H

(α f1,k ) (α f1,k ) (α f2,k ) (α f2,k )

(15)

Using the Cauchy-Schwartz inequality, it is easy to prove α γx21 x2 (wk ) ≤ 1. For a detailed proof, see [10]. Property 2: As the parameter α decreases from 1 to 0, the signal mismatch problem reduces at the expense of a decrease in frequency resolution. 4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Proof: The signal mismatch problem means that the peaks/valleys in the true spectrum of the signal can be canceled if they are located at non-sampled frequencies, which is similar to the well-known problem of direction of arrival (DOA) mismatch in beamforming applications [11, 14, 15]. In the adaptive beamforming processing, the diagonal loading (DL) of the covariance matrix [15, 16, 17] and the notch widening techniques [18, 19] have been proposed to provide more robust performance. To prove the property 2, we will show the relationship of the DL method and the exponential of the covariance matrix (will be referred as E-CM) method by using singular value decomposition (SVD). Assume that the input signal x(n) consists of the desired signal s(n) and the uncorrelated white Gaussian noise (WGN) d(n), which is x(n) = s(n) + d(n).

(16)

The covariance matrix of x(n) is expressed as, Rxx = Rss + Rdd = Rss + λd I.

(17)

where Rss , Rdd are the covariance matrices of the desired signal s(n) and the WGN d(n), respectively. I is the identity matrix, and λd is the power of the WGN . To simplify the analysis, we consider that Rss is rank deficient with rank (Rss ) = P < L. Then the SVD of Rxx becomes Rxx = QT diag {λxi } Q.

(18)

where Q is a unitary matrix, and the eigenvalues λxi are given by ⎧ ⎪ ⎨ λs + λd i

λx i = ⎪ ⎩

f or i = 1, 2, · · · , P f or i = P + 1, · · · , L

λd

.

(19)

where λsi ,for i = 1, 2, · · · , P , are the P eigenvalues of Rss . The Posterior signal-to-noise ratio (PSNR) is defined as P SNRi =

λs i + λd λd

f or i = 1, 2, · · · P.

(20)

For the DL method, the eigenvalues of Rxx + εI can be written as λDL xi =

⎧ ⎪ ⎨ λs + λd + ε i ⎪ ⎩λ +ε d

f or i = 1, 2, · · · , P f or i = P + 1, · · · , L

5

.

(21)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where ε ≥ 0 is the DL level. For the E-CM method, the eigenvalues of Rαxx become = λE−CM xi

⎧ ⎪ ⎨ (λs + λd )α i ⎪ ⎩ (λd )α

f or i = 1, 2, · · · , P f or i = P + 1, · · · , L

.

(22)

Eqs. (21) and (22) show that both the DL and the E-CM methods adjust the PSNR. When the PSNR is much larger than 1, the DL method decreases the PSNR by the DL level ε, and the E-CM method decreases the PSNR by the parameter α if only α is smaller than 1. Thus, the E-CM method may provide more robust performance by properly selecting the parameter α just as the DL method. Furthermore, as the parameter α decreases from 1 to 0, the PSNR larger than 1 will also reduce, thus the robustness will be improved at the expense of a decrease in frequency resolution. We recommend that the parameter α of the E-CM method is selected between 0 and 1.

3

CONCLUSIONS

In this paper, we reveal that the three non-parametric MSC methods, including the Welch-MSC, MVDR-MSC, and CCA-MSC methods, could be treated in a unified way. The uncovered relationship also brings out a new class of MSC estimators in terms of non-linear functions of the covariance matrix, and we further show that the E-CM method is an alternative approach to the DL method in Section 2.3. It seems that most of the PSD estimation methods are suitable to estimate the MSC, e.g. the Welch’s averaged periodogram method, the MVDR method, and the ARMA method. Yet, some of the PSD estimation methods have not been implemented to estimate the MSC, such as the Pisarenko’s spectral estimators [6]. Further work will be concentrated on which kinds of the PSD estimators are suitable for the MSC estimators.

References [1] J.G. Proakis and D.G. Manolakis, Digital Signal Processing: Principles, Algorithms, and Applications, Prentice Hall, Upper Saddle River, NJ, 1996. [2] P. Stoica and R.L. Moses, Introduction to Spectral Analysis, Prentice Hall, New Jersey, 1997. [3] J. Li and P. Stocia, An adaptive filtering approach to spectral estimation and SAR imaging, IEEE Trans. Signal Process. 44(6)(1996) 1469-1484.

6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[4] L. Xu, J. Li and P. Stocia, Complex amplitude estimation in the known steering matrix and generalized wavform case, IEEE Trans. Signal Process. 54(5)(2006) 1716-1726. [5] J. Capon, High-resolution frequency-wavenumber spectrum analysis, Proc. of the IEEE, 57(8)(1969) 1408-1418. [6] V.F. Pisarenko, On the estimation of spectra by means of nonlinear functions of the covariance matrix, Geophy. J. R. Astronom. Soc., 28(1972) 511-531. [7] C. Vaidyanathan and K.M. Buckley, Performance analysis of DOA estimation based on nonlinear functions of covariance matrix, Signal Processing, 50(1996) 5-16. [8] G.C. Carter, Coherence and time delay estimation, Proc. of the IEEE, 72(2)(1987) 236-255. [9] G.C. Carter, C. Knapp, and A.H. Nuttall, Estimation of the magnitude-squared coherence function via overlapped fast fourier transform processing, IEEE Trans. Audio and Electroacostics, 21(4)(1973) 331-344. [10] J. Benesty, J. Chen, and Y. Huang, Estimation of the coherence function with the MVDR approach, in IEEE International Conference on Acoustic, Speech, and Signal Processing, Toulouse, France, May 2006. [11] I. Santamaria and J. Via, Estimation of the magnitude squared coherence spectrum based on reduced-rank canonical coordinates, in IEEE International Conference on Acoustic, Speech, and Signal Processing, Honolulu, Hawaii, USA, April 2007. [12] Y.T. Chan and R.K. Miskowicz, Estimation of coherence and time delay with ARMA models, IEEE Trans. Acoust. Speech Signal Process. 32(2)(1984) 295303. [13] S.V. Narasimhan, G.R. Reddy, E.I. Plotkin and M.N.S. Swamy, Group delay based magnitude square coherence estimation by an ARMA model, Signal Processing, 46(1995) 285-296. [14] H. Cox, Resolving power and sensitivity to mismatch of optimum array processors, Journal of the Acoustical Society of America, 54(3)(1973) 771-785. [15] H. Cox, R.M. Zeskind, and M.M. Owen, Robust adaptive beamforming, IEEE Trans. Acoust., Speech, and Signal Process. 35(10)(1987) 1365-1376. [16] J. Gu, Robust beamforming based on variable loading, Electronics Letters, 41(2)(2005) 55-56. [17] J.R. Guerci, Theory and application of covariance matrix tapers for robust adaptive beamforming, IEEE Trans. Signal Process. 47(4)(1999) 977-985. [18] R.J. Mailloux, Covariance matrix augmentation to produce adaptive array pattern troughs, Electronics Letters, 31(10)(1995) 771-772. [19] M. Zatman, Production of adaptive array troughs by dispersion synthesis, Electronics Letters, 31(25)(1995), 2141-2142.

7