Estimating Correlation Coeﬃcient Between Two
Complex Signals Without Phase Observation
Shigeki Miyabe1(B) , Notubaka Ono2 , and Shoji Makino1
1

University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan{miyabe,maki}@tara.tsukuba.ac.jp
2
National Institute of Informatics (NII)/The Graduate University for Advanced
Studies (SOKENDAI), 2-1-2 Hitotsubashi, Chiyoda-ku, Tokyo 101-8430, Japan
onono@nii.ac.jp

Abstract. In this paper, we propose a method to estimate a correlationcoeﬃcient of two correlated complex signals on the condition that only the
amplitudes are observed and the phases are missing. Our proposed method
is based on a maximum likelihood estimation. We assume that the original complex random variables are generated from a zero-mean bivariate
complex normal distribution. The likelihood of the correlation coeﬃcientis formulated as a bivariate Rayleigh distribution by marginalization overthe phases. Although the maximum likelihood estimator has no analytical form, an expectation-maximization (EM) algorithm can be formulatedby treating the phases as hidden variables. We evaluate the accuracy ofthe estimation using artiﬁcial signal, and demonstrate the estimation of
narrow-band correlation of a two-channel audio signal.
Keywords: Correlation · Complex signal · Maximum likelihood · EM
algorithm

1

Introduction

Correlation of complex sequences plays an important role in array signal processing and almost all of its application [1]. The estimation of correlation is easily
obtained by a simple product sum of the signal sequences. However, the correlation estimation by the product sum cannot be used when the phase observationis unreliable or unavailable. For example, in asynchronous distributed acoustic
sensing systems which gather signals observed by multiple independent recording devices, the biases of the sampling frequencies of the individual devices cause
drift of the phases [2,3]. The eﬀect of the drift on the amplitude is not serious,
but the phase is strongly aﬀected [4]. Also, when we analyze the output of thenonlinear signal processing in the amplitude domain, such as nonnegative matrix
factorization (NMF) [5], the phases are often missing. Although many phase estimation methods are studied [6], the accurate estimation of the correlation cannot
be guaranteed.
The goal of this paper is to estimate a correlation coeﬃcient of two complex
signal channels from the observation of amplitude without information of thec Springer International Publishing Switzerland 2015

E. Vincent et al. (Eds.): LVA/ICA 2015, LNCS 9237, pp. 421–428, 2015.
DOI: 10.1007/978-3-319-22482-4 49

422

S. Miyabe et al.

phase. Since the correlation coeﬃcient is expressed by a nonnegative number,it cannot be used for estimation of phase diﬀerence, and the usage is somewhat limited. Still, the correlation coeﬃcient gives important information, andthe correlation estimation from amplitude is useful for speciﬁc purposes. For
example, the estimation can be used to reduce the computational cost of the
maximum likelihood compensation of drift [3], which requires large computational power and memory. As discussed in [3], only the limited frequency bins
with high correlation contributes much to the observation, and estimation of
the correlation only from the amplitude is informative for the eﬃcient analysis discarding the useless frequency bins. Also, the correlation estimation from
amplitude is expected to be useful for evaluation of channel capacity [7], for theestimation of SNR to optimize of the coeﬃcients of the signal enhancement such
as SS or Wiener ﬁlter, and so on.
To estimate a correlation coeﬃcient between two complex random variables
without phase observation, we propose a maximum likelihood estimation assuming that the original complex variables are generated from a zero-mean bivariate
complex normal distribution. We show that the likelihood is given by the marginalization of the arguments, appearing as a bivariate Rayleigh distribution [7],whose parameter estimation algorithm has not been derived to the best of ourknowledge. We derive an expectation-maximization (EM) algorithm to maximize
the likelihood of the bivariate Rayleigh distribution, and obtain the maximum
likelihood estimator of the correlation coeﬃcient.

2Statement of Problem

Suppose there are two correlated complex random variables Xi ∈ C, i = 1, 2,
which have means of zero, variances σi2 , correlation ρ and the uniform arguments
Θi as
E [Xi ] = 0,2
E |Xi | = σi2 ,(1)

E [X1 X2∗ ] = σ1 σ2 ρ,
1, − π ≤ θi < π,fΘi (θi ) =
2π

(3)

(2)

(4)
∗

where E [·] is the expectation operator, | · | is the absolute value, {·} is thecomplex conjugate, σi2 is the variance of Xi , Θi = ∠Xi , ∠{·} is the argument,
and fA (a) denotes the probability density of a random variable A whose sample
is denoted as a. Hereafter, we denote random variables and the samples withupper and lower case letters, respectively.
Suppose the complex random variables X1 , X2 are unavailable, but we canobserve their absolute values Y1 , Y2 :Yi = |Xi | .

(5)

Estimating Correlation Coeﬃcient Between Two Complex Signals

423

Our goal is to estimate the correlation coeﬃcient |ρ|, the absolute value ofthe cross correlation ρ between X1 and X2 , under the condition that only theabsolute values are observed. It is obvious that the maximum likelihood estimator of the variance σi2 can be obtained as the average of the square of theobservation:
N


1 
2
2
yi (n) ,(6)
σi2 = E |Xi | ←N n=1
where (n), n = 1, . . . , N denotes the index of the N observations. Note that
we omit the sample index (n) when the explicit declaration is unnecessary. In
contrast to the estimation of variance, the correlation coeﬃcient |ρ| is diﬀerent
from the correlation of the absolute observations:|ρ| =|E [X1 X2∗ ]|E [Y1 Y2 ]
=
.
σ1 σ2
σ1 σ2

(7)

Thus the following mean of the product of the observed absolute samples does
not give a good estimation:
|ρ| ←

N

1
y1 (n) y2 (n).
N σ1 σ2 n=1(8)

Therefore, the estimation of the correlation coeﬃcient with the absolute observation is not trivial.

3

Correlation Estimation Assuming Bivariate Complex
Normal Distribution

3.1

Probabilistic ModelIn this section, we discuss the estimation of the correlation coeﬃcient |ρ| fromthe absolute samples y1 , y2 assuming that the unobserved complex samples x
with the statistics given by (1)–(3) are generated from a zero-mean bivariate
complex normal distribution as
σ22 |x1 |2 +σ12 |x2 |2 −2σ1 σ2 Re[ρ∗ x1 x∗2]
exp −σ12 σ22 (1−|ρ|2 )


.
(9)
fX1 ,X2 (x1 , x2 ; ρ) =2
π 2 σ12 σ22 1 − |ρ|Then, the joint density of the observation Yi and the unobserved argument
Θ1 , Θ2 is expressed as
σ 2 y 2 +σ 2 y 2 −2σ σ |ρ|y1 y2 cos(θ1 −θ2 −∠ρ)y1 y2 exp − 2 1 1 2 σ21 σ22 1−|ρ|
2
)
1 2(


fY1 ,Y2 ,Θ1 ,Θ2 (y1 , y2 , θ1 , θ2 ; ρ) =,2
π 2 σ12 σ22 1 − |ρ|
(10)

424

S. Miyabe et al.

by applying the polar coordinate conversion to (9). By the marginalization of
the uniform distribution of the arguments Θ1 , Θ2 , the likelihood of the absolute
observation Y1 , Y2 is given as a bivariate Rayleigh distribution, which can befound in many papers, e.g., [7] as a special case of multivariate Nakagami-m
distributions:
fY1 ,Y2 (y1 , y2 ; |ρ|) = π  π−π

−π

fY1 ,Y2 ,Θ1 ,Θ2 (y1 , y2 , θ1 , θ2 ; ρ) dθ1 dθ2⎛⎞⎞⎛
σ22 y12 + σ12 y22y
4y1 y2
2
|ρ|
y
1
2

 I0 ⎝ ⎠ exp ⎝−⎠,=σ12 σ22 1 − |ρ|2
σ1 σ2 1 − |ρ|2σ12 σ22 1 − |ρ|2

(11)
where Iν (·) denotes the modiﬁed Bessel function of the ﬁrst kind with the order ν.
It can be seen that the density of Y1 and Y2 depends on the correlation coeﬃcient|ρ| but not on the argument ∠ρ of the correlation. Therefore, the maximizationof the likelihood gives the estimation of the correlation coeﬃcient |ρ|. However,
the maximum likelihood estimator does not have the analytical form.
3.2

Maximum Likelihood Estimation by EM Algorithm

Here we describe the maximum likelihood estimation of the correlation coeﬃcient|ρ| by the iterative procedure. By treating the observed samples y1 (n) and y2 (n)together with the unobserved arguments θ1 (n) and θ2 (n), we can formulate the
EM algorithm. We treat the arguments Θ1 , Θ2 as hidden variables, and theposterior density of the hidden variables is given by
fY1 ,Y2 ,Θ1 ,Θ2 (y1 , y2 , θ1 , θ2 ; ρ)fY ,Y (y1 , y2 ; ρ)
 1 2
2|ρ|y1 y2 cos(θ1 −θ2 −∠ρ)exp
σ1 σ2 (1−|ρ|2 )


=
.
1 y2
2πI0 σ σ2|ρ|y2
1 2 (1−|ρ| )fΘ1 ,Θ2 |Y1 ,Y2 (θ1 , θ2 |y1 , y2 ; ρ) =(12)

Then, the auxiliary function Q(|ρ|, |ρ̄|) to maximize in each iteration of the EM
algorithm is obtained asN

Q (|ρ| ; |ρ̄|) =n=1

log fY1 ,Y2 ,Θ1 ,Θ2 (y1 (n) , y2 (n) , θ1 (n) , θ2 (n) ; ρ) θ (n),θ (n)|y (n),y (n);|ρ̄|
1



2 |ρ|


y1 (n) y2 (n) λ (n) − N log 1 − |ρ|2 ,σ1 σ2 1 − |ρ|2 n=1g (a)a|b;c =fA|B (a|b; c) g (a) da,

2

1

2

N

∝

D(a)

(13)
(14)

Estimating Correlation Coeﬃcient Between Two Complex Signals

425

where D(a) is the domain of the variable a, and λ(n) is the posterior expectationof cos(θ1 (n) − θ2 (n) − ∠ρ) with the current parameter estimation |ρ̄|, given by
λ (n) = cos (θ1 (n) − θ2 − ∠ρ)θ1 (n),θ2 (n)|y1 (n),y2 (n);|ρ̄|
⎞
⎛
(n)y
(n)
2
|ρ̄|
y
2
1
⎠ ,
= L⎝
2
σ1 σ2 1 − |ρ̄|
L (x) = I1 (x)/I0 (x)

(15)
(16)

Since L(x) is a monotonically increasing function giving L(0) = 0 and limx→∞
L(x) = 1, λ(n) acts as weighting in the update of the estimation of |ρ| in the
M-step:N1|ρ| =y1 (n) y2 (n) λ (n) .
(17)
N σ1 σ2 n=1
By iterating the updates of E- and M-steps given by (15) and (17), respectively,
the estimation of |ρ| converges to a local optimal. Note that σ12 and σ22 areestimated by (6), which can also be derived as the maximization of the auxiliary
function, although the related terms are omitted in (13). Also note that L(x)can be calculated by one-dimensional table lookup.

4

Experimental Results

4.1

Evaluation with Artificial Signal

To evaluate the performance of the proposed method, we conducted a numericalsimulation to estimate the correlation coeﬃcients of artiﬁcial two-channel complex signals generated from pseudorandom numbers. To show the baseline, we
also evaluated the performance of absolute correlation given by (8). In addition,
to show the upper limit with the ideal condition, we evaluated the standard
maximum likelihood estimation under the condition where the original complex
sequences is available, given by
|ρ| ←

N

1
∗
x1 (n) x2 (n) .
N σ1 σ2 n=1(18)

As the evaluation criterion, we calculated the root mean squared errors (RMSEs).
The correlated data are generated by the linear mixture of the two independent pseudorandom numbers with the same variance. We controlled the correlation by manipulating the linear mixture. To evaluate the robustness of the
proposed method against the mismatch of the Gaussian assumption, we also evaluated the artiﬁcial signals generated by super-Gaussian pseudorandom numbers.
The super-Gaussian data are generated from the circular generalized normal distribution [8], whose density of the random variables X1 , X2 is given by   ci 
ci exp − |xξii | 
fXi (xi ) =,(19)
2πξi2 Γ c2i

426

S. Miyabe et al.

where ξi > 0 is the scale parameter, ci > 0 is the shape parameter, and Γ (·)
is a gamma function. The shape parameters are set as c1 = 0.5 and c2 = 0.8,
and scale parameters are adjusted to give the unit variance. With such shapeparameter setting, the sequences have super Gaussian property with long tail,and their kurtosis, a well-used Gaussianity measure, are about 7.4 and 2.3. Notethat the kurtosis is changed after the mixing to give correlation.
We show an example of the estimation in Fig. 1. For both Gaussian andsuper-Gaussian data, the absolute correlation does not give accurate estimation. Although the variance of the proposed method is slightly larger than theestimation with arguments, the estimated correlation coeﬃcients distribute near
the true correlation coeﬃcients and the accuracy of the proposed method ismatch better than the baseline. The estimation of the correlation coeﬃcients ofsuper-Gaussian data tends to underestimate the correlation, but the RMSEs are
similar to the Gaussian case without the model mismatch. Thus it is conﬁrmedthat the proposed method can eﬀectively estimate the correlation coeﬃcients ofthe complex sequences without the observation of arguments.
To examine the eﬀect of the bias and variances of the estimation, we evaluate RMSEs of various numbers of observed samples. The result is shown in
Fig. 2. Although the proposed method has larger variance than the ideal estimation with phase, we can see that the estimation accuracy of the proposed
method improves according to the increase of the number of the samples when
the data is Gaussian. Thus the proposed method can estimate the correlationof Gaussian data eﬀectively with small bias. However, we can see the saturation of the improvement of the accuracy under the condition mismatch with thesuper-Gaussian data.
4.2

Demonstration with Audio Data

As a demonstration of the correlation coeﬃcient estimation in practical signal
processing, we evaluated the estimation of the correlation coeﬃcients of thenarrow band amplitudes of a two-channel audio signal.
We analyzed observation of speech mixture by an array of two microphones.
The data is chosen from the Underdetermined Test dataset of the Signal Separation Evaluation Campaign (SiSEC) [9], a benchmark of speech separation.
Utterances of four female speakers were recorded in a room whose reverberation
time T60 is 250 ms. The spacing of the microphones was 1 m. The recorded data
was 10 s long, sampled with the frequency of 16 kHz. The signal was analyzed
by short-time Fourier transform with the von Hann window of the length 1024
samples and the shift 128 samples. The number of the frames is 1258.
The estimated results are shown in Fig. 3. In contrast to the artiﬁcial signals,
the true correlation coeﬃcients are unknown. Thus it should be noted that thehorizontal axis is not the true correlation coeﬃcients but the estimated correlation coeﬃcients by the ideal estimation in (18) with the original complex signal
given. We can see that the proposed method is much better than the baseline.
Superiority of the proposed method holds for other frame lengths, which stronglyEstimating Correlation Coeﬃcient Between Two Complex Signals

427

Fig. 1. An example of the estimation results of correlation coeﬃcients between two
artiﬁcial signals generated from (a) Gaussian and (b) super-Gaussian pseudorandom
numbers. The number of samples was 100 for each trial. The EM iteration number was
10. The RMSEs are about 0.42, 0.05 and 0.12 for the baseline, the ideal estimation andthe proposed method, respectively for the results in (a), and about 0.05, 0.37 and 0.12
for the results in (b).

Fig. 2. Root mean squared errors for the number of the samples N = 10, 20, 50,
100, 200, 500, 1000, 2000, 5000, 10000. The RMSEs of each condition was calculated from100 trials. The EM iteration number was set to 200.

Fig. 3. Estimation results of narrow band correlation coeﬃcients of a two-channel audio
signal. The EM iteration number was 10. The RMSEs from the ideal estimation areabout 0.28 and 0.20 for the baseline and the proposed method, respectively.

aﬀects the correlation. However, the accuracy is not as good as that of the artiﬁcial signal. Thus, further analysis is required to improve the estimation accuracy
with the realistic signals.

428

5

S. Miyabe et al.

Conclusions

In this paper, we proposed an EM algorithm to obtain the maximum likelihood
estimation of the correlation coeﬃcient of the two correlated complex randomvariables only with the observation of the absolute values. Assuming that theoriginal complex data are generated from a zero-mean bivariate complex normal
distribution, we formulated the likelihood of the correlation coeﬃcient. Although
the maximum likelihood estimator is not analytical, we formulated the EM algorithm by treating the diﬀerence of the arguments as a hidden variable. We evaluated accuracy and the robustness against the model mismatch of the proposed
method by the simulation using artiﬁcial signals. Also, we demonstrated theestimation of the narrow-band correlation coeﬃcients of the two-channel audio
signal. It is conﬁrmed that the proposed method is much more accurate thanthe estimation with the correlation of the absolute values. However, estimationof the real signal was worse than that of the artiﬁcial signal, and the further
analysis and improvement are required. Our future work includes the modiﬁed
maximum likelihood estimation with a non-Gaussian distribution.
Acknowledgment. This work was supported by JSPS KAKENHI Grant Number
23240023.

References
1. Johnson, D., Dudgeon, D.: Array Signal Processing: Concepts and Techniques.
Simon & Schuster, New York (1992)
2. Markovich-Golan, S., Gannot, S., Cohen, I.: Blind sampling rate oﬀset estimationand compensation in wireless acoustic sensor networks with application to beamforming. In: Proceedings of IWAENC, pp. 1–4 (2012)
3. Miyabe, S., Ono, N., Makino, S.: Blind compensation of interchannel sampling frequency mismatch for ad hoc microphone array based on maximum likelihood estimation. Signal Process. 107(2015), 185–196 (2015)
4. Chiba, H., Ono, N., Miyabe, S., Takahashi, Y., Yamada, T., Makino, S.: Amplitudebased speech enhancement with nonnegative matrix factorization for asynchronous
distributed recording. In: Proceedings of IWAENC, pp. 204–208 (2014)
5. Lee, D., Seung, H.: Algorithms for non-negative matrix factorization. In: Proceedings of NIPS, pp. 556–562 (2001)
6. Griﬃn, D., Lim, J.: Signal estimation from modiﬁed short-time Fourier transform.
IEEE Trans. ASSP 32(2), 236–243 (1984)
7. Fraidenraich, G., Lévêque, O., Cioﬃ, J.: On the MIMO channel capacity for the
nakagami-m channel. IEEE Trans. Inf. Theory 54(8), 3752–3757 (2008)
8. Novey, M., Adali, T., Roy, A.: A complex generalized Gaussian distributioncharacterization, generation and estimation. IEEE Trans. Signal Process. 58(3),1427–1433 (2010)
9. Ono, N., Koldovsky, Z., Miyabe, S., Ito, N.: The 2013 Signal Separation Evaluation
Campaign. In: Proceedings of MLSP, pp. 1–6 (2013)