PARAMETRIC-PEARSON-BASED INDEPENDENT COMPONENT ANALYSIS FORFREQUENCY-DOMAIN BLIND SPEECH SEPARATION*

Hiroko Kato, Yuichi Nagahara, Shoko Araki, Hiroshi Sawada and Shoji Makino
NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation
2-4 Hikaridai, Seika-cho, Soraku-gun, 619-0237, Kyoto, Japan
phone: + (81) 774 93 5138, fax: + (81) 774 93 1945, email: katohi@cslab.kecl.ntt.co.jp
web: www.kecl.ntt.co.jp* School of Political Science and Economics, Meiji University
1-1, Kanda-Surugadai, Chiyoda-ku, 101-8301, Tokyo, Japan
phone: + (81) 3 3296 2118, fax: + (81) 3 3296 2350, email: nagahara@kisc.meiji.ac.jp

ABSTRACTSeparation performance is improved in frequency-domain
blind source separation (BSS) of speech with independent
component analysis (ICA) by applying a parametric Pearson distribution system. ICA adaptation rules include ascore function determined by approximated source distribution, and better approximation improves separation performance. Previously, conventional hyperbolic tangent(tanh) or generalized Gaussian distribution (GGD) was
uniformly applied to the score function for all frequency bins,
despite the fact that a wideband speech signal has different
distributions at different frequencies. To obtain better scorefunctions, we propose the integration of a parametric Pearson distribution system with ICA learning rules. The scorefunction is estimated by using appropriate Pearson distribution parameters for each frequency bin. We consider threeestimation methods with Pearson distribution parametersand conduct separation experiments with real speech signals convolved with actual room impulse responses. Consequently, the signal-to-interference ratio (SIR) of the proposed methods significantly improve over 3 dB compared toconventional methods.
1.

INTRODUCTION

Frequency-domain blind source separation (BSS) based on
independent component analysis (ICA) [1][2] has been proposed for the BSS of convolutive mixtures of speech signals[3][4][5][6][7]. ICA adaptation rules have a score function
determined by approximated source distribution. In
conventional applications of ICA to speech signals, superGaussian distribution is used as a distribution model, typically with the score function tanh [1][2]. Recently, a GGDbased modeling approach that is more adaptable than tanh
has also been presented [8][9]. In these modeling approaches, the score function is uniformly applied to all frequencies.
However, when transforming the speech signals to the frame
series (i.e., the time sequence in the frequency domain) by
short-time Fourier transform (STFT), it is clear that the distribution of the frame series illustrates different patterns foreach frequency. The shapes seem to have a non-Gaussian
appearance, as shown by fat-tailed and skewed characteristics, that resembles various distribution shapes of the Pearson distribution system [10]. A previous study investigated
the application of the Pearson distribution system to ICA(Pearson-ICA) [11]. This approach showed better separation
performance than such conventional nonlinear functions astanh. Furthermore, a nonparametric ICA approach to estimating the source distribution was proposed, and its separation performance was compared to those of several methods
such as Pearson-ICA, Fast ICA, and Kernel-ICA [12]; however, [11] and [12] were performed in the time-domain andused artificial data. In order to achieve BSS of convolutive
mixtures, which is the focus of this paper, time-domain
processing is inefficient.
Therefore, this article proposes a Pearson distribution system approach to frequency-domain BSS. In this approach,the source distribution for each frequency bin is modeled byan adaptive parametric Pearson distribution, and the scorefunction is formed by the parameters of the distribution to
improve separation performance. So far, practical implementation of distribution parameters that depend on the
moment of the signal has been unwieldy. To overcome suchproblems, [13], [14] and [15] have proposed the solution of
introducing a discrimination method for Pearson distribution
types and transforming formulae between the moment anddistribution parameters. Accordingly, we also employ a new
implementation in the proposed score function’s estimation
procedure.
This paper is organized as follows: Section 2 introduces thebasic framework of BSS and conventional nonlinear functions. Section 3 outlines the parametric Pearson distribution
system. Section 4 describes our proposed BSS methods.
Section 5 shows experimental results, and conclusions are
given in Sec. 6.
2.

BLIND SOURCE SEPARATION OF SPEECH

In this paper, we consider the blind source separation (BSS)of speech signals observed in real environments, i.e., the
BSS of convolutive mixtures of speech. In such environ-ments, N source signals si ( n) are observed with their reverberant components and delays by M sensors. Therefore,
observations are modeled as convolutive mixtures:NP

x j (n) = ∑∑ h ji ( p ) si (n − p + 1) ( j = 1,

,M ) ,(1)

i =1 p =1

h ji ( n) : P − taps impulse response from source ito sensor j.
Our goal is to obtain separated signals yk (n) ( k = 1, , N )using only the information provided by observations x j (n) .

In this paper, we handle the case of N = M = 2 (Fig. 1);
however, we can expand our method to any N = M casewithout loss of generality.
This paper employs a frequency-domain approach, that is, ashort-time Fourier transform (STFT) is performed to convert
our problem into a linear instantaneous mixture at each frequency. In the frequency domain, mixtures (1) are modeled
asscore function. If source distributions p are known,score functions (5) are defined as [1][2]:p '(| Y |)
Φ (Y ) = −exp( j∠Y ) ,(5)p(| Y |)
Y : complex number, | ⋅ | : absolute value, ∠Y : argument,

where p '( x) = dp( x) dx . In blind separation, however,source distribution cannot be obtained a priori. Usually, thescore function has to be approximated, and the conventional
tanh is widely used at all frequencies for speech separation
because speech signals have a super-Gaussian distribution[1][2]:

Φ(Y ) = tanh( g | Y |)exp( j∠Y ) .

As mentioned in Sec. 1, conventional GGD has also been
applied to BSS [8]. The GGD-based score function, uniformly applied to all frequencies [9], is represented byN

X j ( f , m) = ∑ H ji ( f ) Si ( f , m) ,Φ (Y ) =| Y |β −1 sign(| Y |) exp( j∠Y ).

(2)

i =1

where f denotes a frequency and m is the frame index. Withmatrices, (2) can be written asX( f , m) = H ( f )S( f , m),
H ( f ) : M × N mixing matrix,
S ( f , m) = [ S1 ( f , m),

T, S N ( f , m)] :

(6)

(7)

β : shape parameter .
As is well known in speech signal processing/engineering,for β = 1 , 2, and 0.5, the GGD becomes, respectively, aLaplacian distribution, whose speech closely follows it, astandard Gaussian distribution, and a Gamma distribution.

STFT of source signal,X( f , m) = [ X 1 ( f , m),

T, X N ( f , m)] :STFT of observed signal.
In a blind scenario, H( f ) and S( f , m) are unknown.
The separation process can be formulated in each frequency f :Y ( f , m) = W ( f ) X( f , m) ,Y ( f , m) = [Y1 ( f , m),(3)

T, YN ( f , m)] :estimated source signal vector,
W ( f ) : N × M separation matrix.

W( f ) is determined so that Y1 ( f , m) , , YN ( f , m) become mutually independent using ICA. After getting theseparated signals (3), we convert the frequency-domain signal Yk ( f , m) into a time-domain signal by using inverse
STFT.
The separation matrix is independently estimated at eachfrequency. An algorithm based on a natural gradient [16] iswidely used for this. The adaptation rule of the i -th iterationis
Wi +1 ( f ) =(4)
Wi ( f ) + η[I − Φ (Y ( f , m))Y H ( f , m) ] i Wi ( f ),Y ( f , m) : average with respect to m,

Figure 1. Frequency-domain speech BSS system ( N = M = 2 )3.

PARAMETRIC PEARSON SYSTEMAPPROACH

The above conventional nonlinear function approach is uniformly applied to all frequencies; however, the actual speech
signal for each frequency has a different distribution. Therefore, Fig. 2 shows the distribution shapes a, b, c, d, e, and f
for bins 3, 5, 30, 125, 250, and 500, respectively. To adapt
these different shapes, we apply the Pearson distribution
system, which is widely used to model such various distributions as Gaussian, Student’s t, gamma, and beta. Distribution parameters are detected by the sample moments asmean, variance, skewness (Skew), and kurtosis (Kurt).
Pearson [10] defined the differential equation related toprobability density function p(| Y |) . If the random variable
is complex value Y , the form is defined by

H : transpose conjugate, η : adaptation stepsize.

Here, Φ(Y ) ( Y ( f , m) is simplified by Y ) indicate the−b0 + b1 | Y |p '(| Y |)exp( j∠Y ) .
=p (| Y |) c0 + c1 | Y | +c2 | Y |2(8)

ing the Pearson system to frequency-domain BSS, our proposed method sets forms (8) and (10) as the score functions.
Pearson IV(κ >1 )
Pearson IV( 0 < κ <1 )Pearson I(0>κ )

Figure 2. Histograms a, b, c, d, e and f of the frame series for several
frequency bins f =3, 5, 30, 125, 250, and 500. STFT frame size: 512,
Sampling rate: 8 kHz

Note that form (8) has the same shape as score function (5)of the ICA. That is, if the coefficients of (8) can be estimatedby an appropriate method through the observed data at eachfrequency, we can obtain the score function to approximate
the source distribution at each frequency. For estimation ofseveral kinds of distribution families in a Pearson system,[13] introduced parameter κ using Skew and Kurt obtainedby the data:κ=(Skew) (Kurt + 3)
. (9)4(2 × Kurt − 3 × (Skew) 2 − 6)(4 × Kurt − 3 × (Skew) 2 )2

2

For 0 > κ , 0 < κ < 1 , 1 < κ , the types are discriminated to I,
IV, and VI, respectively. In Fig. 2, panels a and b represent
Types IV and VI. Panels c to f represent the shapes of Type I.
In our preliminary consideration of the STFT series of realspeech data, we calculated κ values for each frequency bin(Fig. 3). Figure 3 shows the distribution of the frame series asclassified by Types I, VI, and IV. Type I was widely detectedin the mid- and high-frequency bins; however, the height and
tail of the distribution were different for each bin, as shownin panels c, d, e, and f of Fig. 2. The distribution figures areeach different J-shaped distributions.
Then the score functions of Types I, IV and VI were applied,
as described below in terms of distribution parameters:
I : Φ (| Y |) =IV : Φ (| Y |) =VI : Φ (| Y |) =−( p + q + 2) | Y | + ( p + q − 2) a + ( p − 1)b| Y | −(2a + b) | Y | + a ( a + b)2

2b | Y | −2bµ − 2bτδ| Y | −2 µ | Y | − µ + τ222exp( j∠Y ),(c + 1) | Y | −(c + 1) a − ( β − 1)α| Y | −(2a − α ) | Y | + a ( a − α )2exp( j∠Y ),exp( j∠Y ).

(10)
Parameters p, q, b, a, µ ,τ , δ , c, β , and α can be calculated by
the moments of the frame series (see [13][14]). When apply-

Figure 3. κ values for frequency bin. The horizontal axis indicatesfrequency bin within 512 STFT frame size and the vertical axis
indicates κ value calculated by (9), which was applied to each series obtained by STFT.

4.

PROPOSED METHODS

To estimate the score function mentioned above, we propose
the following three methods:
Method 1: Minimization of cross-correlation
Here, we use the score function (8). To estimate the Pearson
parameters {b0 ( f ), b1 ( f ), c0 ( f ), c1 ( f ), c2 ( f )} , we select theparameters that minimize the sum of the absolute values ofthe off-diagonal components of [I − Φ (Y)Y H ] in (4); thatis,∑ | φ (Y ( f , m)Y ( f , m) | (i ≠ j ) .
i

i

*
j

(11)

These off-diagonal components represent the higher-ordercross-correlation of the outputs. If output signals are wellseparated, they become mutually independent, and the valueof (11) becomes 0.On the other hand, when the separation isincomplete, the absolute value of off-diagonal components
is far from zero. Therefore, we can use the off-diagonalcomponents as measures of separation performance. According to this measure, we search for the Pearson system parameters {b0 ( f ), b1 ( f ), c0 ( f ), c1 ( f ), c2 ( f )} that minimize(11) in an arbitrary range with a grid search. Using theseparameters, we determine the score function of ICA corresponding to (8) and estimate the separating matrix with (4).
In this paper, Pearson parameters are determined by a gridsearch.
Method 2: Estimation of appropriate Pearson distribution
typeThis method directly determines the appropriate Pearson
type for each bin in the learning process of the separation
matrix. Score function (10) is used in this method. The
specific calculation procedure for each frequency proceeds

specific calculation procedure for each frequency proceeds
as follows:ˆ ( f ) using (6) and set ini1) Estimate separation matrix W
ˆ (f);
tial value W ( f ) = W
02) Calculate κ (see (9)) by the skewness and kurtosis of
the absolute value of Y( f , m) obtained with (4);3) Following κ , the appropriate Pearson distribution typeis specified and the parameters of the score function defined in (10) are calculated by the moments of the STFT
frame series according to [13] [14];4) Renew W ( f ) by (4); and
5) Iterate procedures 2) to 4) until there is a convergenceof (4).
Compared with Method 1, the computational burden is significantly reduced because it is not necessary to perform gridsearch.
Method 3: Combining Methods 1 and 2
In Fig. 3, the κ values for low-frequency bins are unstable.
The distribution shape of the frame series at high frequency
is illustrated by the various J-shape patterns shown as TypeI. In a preliminary investigation, we found that the individual histograms, corresponding to the frequencies of the estimated parameters b0 ( f ) , b1 ( f ) , c0 ( f ) , c1 ( f ) , and c2 ( f ) ,have similar tendencies for all speaker combinations at highfrequency. On the other hand, at low frequency, the distribution types depend on each speaker. From this fact, we propose another method, Method 3, that combines Methods 1
and 2. Here, the score function (10) is applied at low
frequencies, while score function (8) is applied at high frequencies. However, score function (8) is not estimated by agrid search; instead, the pre-estimated mean value for eachparameter is used for (8), thus significantly reducing thetime needed for calculation. Concretely, the procedure is1) Calculate mean values {b0 ( f ), b1 ( f ), c0 ( f ), c1 ( f ), c2 ( f )}of parameters estimated by applying Method 1 to arbitrary data combinations;2) Define f 0 as the boundary point;3) Apply Method 2 to low-frequency f ≤ f 0 according to
the appropriate Pearson type for each frequency bin;and4) To high-frequency f > f 0 , input averaged parameters{b0 ( f ), b1 ( f ), c0 ( f ), c1 ( f ), c2 ( f )} for each bin directly

into (8).
For choosing the best f 0 in advance, we compared SIR values when using f 0 between 0 and 200 bins and selected
the f 0 that provided the highest SIR.
5.

EXPERIMENTAL RESULTS

5.1 Experimental conditions
We conducted separation experiments with real speech signals and measured room impulse responses. The speech data
were convolved with impulse responses measured in an actual room (Fig. 4) whose reverberation time was 130 ms.

Figure 4. Room layout used for experiments

As original speech, we used Japanese sentences of 3 seconds
spoken by male and female speakers. We made observation
signals with (1) and investigated four combinations of speakers. The STFT frame size was 512, and the frame shift was
256 at a sampling rate of 8 kHz. To solve the permutation
problem of frequency-domain ICA, we employed the direction of arrival and correlation approach [7], and to solve thescaling problem we used the minimum distortion principle[17]. For numerical analysis, we arranged four data sets:
female and female (f & f), two types of female and male (f& m), and male and male (m & m). As an evaluation measure, we used the signal-to-interference ratio (SIR) as a separation performance measure:∑ n yii2 (n) ,(12)
SIR i = 10log10∑ k ≠i (∑ n yik (n))2yik (n) = ∑ l uik (l ) sk (n − l ), uik (l ) = ∑ j =1 ∑ t =0 wij (t )h jk (l − t ),M

L −1

wij (t ) : inverse discrete Fourier transform of
ij component of separation matrix W( f ),L : STFT frame size.
5.2
Results
Results using Methods 1, 2, and 3, conventional tanh, and
GGD-based modeling methods with the above four types ofdata sets are summarized in Table 1. The SIR for GGD wasthe value obtained for the best SIR for β ∈ [0.5,1.0] . For
conventional nonlinear functions, the GGD-based modelingmethod was slightly better than tanh. For the proposed Pearson approach, we obtained maximum improvement in separation performance of around 3 dB better than conventional
tanh and around 2 dB better than conventional GGD.
Method 3, combining Methods 1 and 2, showed better separation performances than that by using only Method 2, eventhough the mean parameters used in Method 3 were firstestimated by using only two data combinations. The mean
parameters were applied to all data combinations, including
the data that were not used to estimate the mean parameters.
In fact, the SIR values sometimes became unstable in thecalculation applying only Method 2. Therefore, Method 3’s
results suggest that using the mean parameter values preestimated by Method 1 at high frequencies did not lead toinstability and that at low frequencies the parameters estimated with data moments worked well.
Furthermore, we considered the computational time needed
to perform these methods. We obtained the results using
Matlab® profile report and summarize them in Table 2. The
CPU clock speed was 594 MHz. Methods applying conventional nonlinear functions to the score function were fasterthan Methods 1 and 2. In particular, Method 1 took about
two hours because the calculation algorithm included a gridsearch to obtain optimized parameters. By reducing thisoptimization procedure, Method 3 could work at a reasonable computation speed, thus improving performance.
Table 1: SIR (dB) values for conventional nonlinear functions andthree proposed methods

m&m
f&m 1
f&m 2
f&f

tanh
17.71
14.99
15.92
17.27

GGD
17.58
15.89
17.22
17.76

Meth. 1
17.00
15.99
18.61
20.12

Meth. 2
17.54
16.38
17.08
18.62

Meth. 3
18.12
16.42
18.38
19.43

Table 2: Computational time required to perform BSS

conventional nonlinear function
Proposed
Pearson systemapproach6.

Methods
tanh
GGD
Meth. 1
Meth. 2
Meth. 3Time [s]
10.38
18.09
6888.66
56.11
13.07

CONCLUSION

To estimate the frequency-domain separation matrix for ICA,
we proposed a practical parametric Pearson distribution system. This system detects a score function for the source distribution at each frequency. We first confirmed the efficiency of three methods developed to apply the Pearson
system to frequency-domain speech BSS under blind conditions. As the first method, unknown parameters were estimated to minimize the cross-correlation of the separation
matrix. In the second method, we directly calculated thetransform formulae based on κ discrimination. The third
method was a combination of these two methods. The proposed approach significantly improved separation performance compared with conventional tanh and GGD-basedmodeling approaches. Furthermore, the combined method
showed better performance than applying the individual
methods, and its computational speed was reasonable.

REFERENCES[1] S. Haykin, ed. Unsupervised Adaptive Filtering (Volume
I: Blind Source Separation), John Wiley & Sons, 2000.

[2] A. Hyvärinen, J. Karhunen, and E. Oja, Independent
Component Analysis, John Wiley & Sons, 2001.
[3] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Frequency-domain blind source separation,” in Speech Enhancement, eds. J. Benesty, S. Makino and J. Chen, Springer,
2005.
[4] P. Smaragdis, “Blind separation of convolved mixtures in
the frequency domain,” Neurocomputing, vol. 22, pp. 21-34,
1998.
[5] L. Parra and C. Spence, ”Convolutive blind separation of
nonstationary sources,” IEEE Trans. SAP. vol. 8, no. 3, pp.
320-327, May 2000.
[6] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Polar
coordinate based nonlinear function for frequency domain
blind source separation,” IEICE Trans. Fundamentals, vol.
E86-A, no. 3, pp. 590-596, March 2003.
[7] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust
and precise method for solving the permutation problem offrequency-domain blind source separation,” IEEE Trans.
SAP. vol. 12, no. 5, pp. 530-538, Sept. 2004.
[8] K. Kokkinakis and A.K. Nandi, “Multichannel speech
separation using adaptive parametrization of source PDFs,”
eds. C.G. Puntonet and A. Prieto, ICA 2004, LNCS 3195,
Springer-Verlag Berlin Heidelberg, pp. 486-493, 2004.
[9] R. Prasad, H. Saruwatari, and K. Shikano, “Blind separation of speech by fixed-point ICA with source adaptive
negentropy approximation,” IEICE Trans. Fundamentals, vol.
E88-A, no. 7, July 2005.
[10] K. Pearson, “Memoir on skew variation in homogeneous material,” Philos. Trans. Roy. Soc. A, vol. 186, pp. 343414, 1895.
[11] J. Karvanen, J. Eriksson, and V. Koivunen, “Pearson
system based method for blind separation,” Proc. the Second
International Workshop on ICA and BSS, pp. 585-590, 2000.
[12] R. Boscolo, H. Pan, and V. P. Roychowdhury, “Independent component analysis based on nonparametric density
estimation,” IEEE Trans. NN., vol. 15, no. 1, pp. 55-65, Jan.
2004.
[13] Y. Nagahara, “Non-Gaussian filter and smoother basedon the Pearson distribution system,” J. Time Ser. Anal, vol.
24, no. 6, pp. 721-738, 2003.
[14] Y. Nagahara, “The PDF and CF of Pearson type IV distributions and the ML estimation of the parameters,” Stat.
Prob. Letters, vol. 43, pp. 251-264, 1999.
[15] Y. Nagahara, “A method of simulating multivariate nonnormal distributions by the Pearson distribution system andestimation”, Comp. Stat. Data. Anal, vol. 47, issue 1, pp. 129, 2004.
[16] S. Amari, “Natural gradient works efficiently in learning,” Neural Computation, vol. 10, pp. 251-276, 1998.
[17] K. Matsuoka and S. Nakashima, “Minimal distortion
principle for blind source separation,” Proc. ICA2001, pp.
722-727, Dec. 2001.