A POLAR-COORDINATE BASED ACTIVATION FUNCTION
FOR FREQUENCY DOMAIN BLIND SOURCE SEPARATION
Hiroshi Sawada

Ryo Mukai

Shoko Araki

Shoji Makino

NTT Communication Science Laboratories
2-4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japanfsawada,ryo,shoko,makig@cslab.kecl.ntt.co.jp
ABSTRACTThis paper presents a new activation function for an ICA
algorithm to process complex-valued signals, which is usedin frequency domain blind source separation. The new activation function is based on the polar coordinates of a complex number, whereas the conventional one is based on the
Cartesian coordinates of a complex number and calculates
the real part and imaginary part separately. The new activation function eliminates an undesirable constraint occurredby the conventional function. In experiments for separating
speech signals in a reverberant environment, we obtained
improved SNRs by using the new activation function.

on the polar coordinates of a complex number, and eliminates the additional constraint discussed above. Experimental results have shown that the new activation function
works well and SNRs (signal-to-noise ratios) are improved
over conventional methods. We discuss several reasons why
improved results were obtained, by looking into what happens in ICA algorithms.
In Section 2, we explain frequency domain BSS in detail. The new activation function is proposed in Section 3.
In Section 4, we show and discuss experimental results. We
conclude this paper in Section 5.
2. FREQUENCY DOMAIN BSS

1. INTRODUCTION
2.1. Problem Formulation
Blind source separation (BSS) is a technique to estimate
original source signals using only sensor observations thatare linear mixtures of the original signals. Independent component analysis (ICA) [1–4] works very well for BSS, if the
mixture is instantaneous (non-convolutive). In a real room
environment, however, sounds are mixed in a convolutive
manner with reverberation, and longer reverberation makes
a BSS problem more difficult. Several methods have been
proposed for convolutive mixtures. One of the major methods is frequency domain BSS [5–14].
In frequency domain BSS, a convolutive mixture problem of the time domain is converted into multiple instantaneous mixture problems of the frequency domain. Then,these instantaneous mixture problems are solved by ICA
in every frequency bin. Since the conversion is performedby using a windowed discrete Fourier transform (DFT), we
have to deal with complex numbers. An activation function
for an ICA algorithm to process complex numbers was proposed [5]. This function, however, imposes an additional
constraint that prevents a learning algorithm from converging. To avoid the additional constraint, another formula
to calculate a gradient of an ICA algorithm was proposed[6, 7], and has been used in [8, 9, 13, 14].
In this paper, we propose a new activation function for
an ICA algorithm to process complex numbers. It is based( ) (1p
Suppose that there are N source signals sp t ;N that are mutually independent, and these signals are observed at M microphones in a real room environment withreverberation. The observed signals cab be written in a convolutive mixture form:)

X( ) = h  ( ) (1   )
N

xq t

qp

sp t ;

q

M ;

p=1

h

where qp represents the impulse response from source p to
microphone q , and  denotes the convolution operator.
The goal of blind source separation is to separate observed signals xq t into N unmixed signals yp t ;p  N that are mutually independent. The separation
has to be done without knowing impulse responses qp nor
original source signals sp t . The unmixing system can becomposed of N  M FIR filters. The unmixed signals are
obtained by())( ) (1

h

()

()=

yp t

w

X w  ( ) (1   )Mpq

xq t ;

p

N ;

q =1

where pq represents the coefficients of the FIR filters. Figure 1 shows a BSS system for the case of N = M = 2.

663

sourcesignals

observed
signals

h11

s1

x1

+

h21

1

unmixedsignals

w11

0.8

y1

+

w21

0.6
0.4

s2

w12

x2

+

h22

mixing system

y2

+

w22

tanh(Y)

0.2

h12

0
−0.2

unmixing system

−0.4
−0.6

Fig. 1. Model of a BSS system

−0.8
−1
−4

−3

−2

−1

(

) =

Xq !; mX()(

e j!t xq t win t

( )

W

i+1

= W i + W i

W

W =  [I h(Y )Y H i] W ;(2)
(Y ) = tanh[  re(Y )] + j  tanh[  im(Y )]; (3)

(Y )

Y

Y

(Y )

W

h(Yp )Yq i = 0 (p 6= q);h(Yp )Yq i = 1 (p = q);Tpq

In an ICA algorithm, an unmixing matrix
improved by the learning rule:

W

where H represents the conjugate transpose of , andreand im
are the real and imaginary parts of ,respectively.
Based on (2),converges to a point that satisfies

Y (!; m) = W (!)X (!; m);where X (!; m) = [X (!; m); : : : ; X (!; m)] . Then,
we can obtain FIR filters w of length T by applying theinverse DFT to all W (! ).
2.3. ICA algorithm

(1)



W( )
M

tanh(Y )is widely used [1, 2] as a nonlinear activation function, where is a scaling parameter to control the nonlinearity of .
In frequency domain BSS, the signals obtained by a DFTandare complex numbers. Thus, the calculations of has been extended for complex numbers [5]:

w

1

4

(Y ) = tanh(  Y )

)

mS ;

where win denotes a window function and S is the shifting
interval of the window. The frame length T of the window
is the same as the number of frequency bins, and also thelength of the FIR filters pq composing the unmixing system.
For each frequency ! , an ICA algorithm is applied to! and frequency
obtain an unmixing N  M matrix
domain N signals!; m , which are estimations of thesource signals at the frequency:

)

3

( )

(! = 0; T1 2; : : : ; T T 1 2)Y(

2

In this formula,  is a step size parameter that has an effecton the speed of convergence, hi denotes the averaging operator, and  is an activation function. The hyperbolic
tangent (Fig. 2)

()

t=0

1

Fig. 2. Hyperbolic tangent

A convolutive mixture in the time domain corresponds to an
instantaneous mixture in the frequency domain. Therefore,
we can apply an ordinary ICA algorithm in the frequency
domain to solve a BSS problem. A T -point windowed DFT
is used to convert time domain signals xq t into frequency
domain time-series signals:
T 1

0

Y

2.2. Framework of Frequency Domain BSS

Y

(4)
(5)

where Yq is the complex conjugate of Yq . The first equation (4) concerns the mutual independence of Yp and Yq .
The second equation (5) makes the average amplitude of Yp
converge to some value near 1 since the range of is fromto 1. Decomposing (5) into real and imaginary parts, we
have

W is gradually1h[re(Yp )]re(Yp ) + [im(Yp )]im(Yp )i = 1h[im(Yp )]re(Yp ) [re(Yp )]im(Yp )i = 0

Ybased on the minimization of the mutual information of[1, 2]. For the calculation of, the natural gradient [3]is widely used:

(6)
(7)

( )
( )

Equation (7) imposes an additional constraint that re Yp
and im Yp should be mutually independent [7]. This constraint is too strong since there may be a case that re Yp
and im Yp are not mutually independent.

( )
( )

W =  [I h(Y )Y T i] W :664

5.73 m1
1.56 m1.15 m3.12 mtanh(abs(Y))

0.8
0.6

loudspeakers(height : 1.35 m)

30

2.15 m

40

4 cm

0.4

microphones 1.15 m(height : 1.35 m)

0.2
0
4

room height : 2.70 m
4

2
2

0
−2
−4

imaginary part of Y

Fig. 4. Layout of a room used to record impulse responses

0

−2
−4

Fig. 3. Hyperbolic tangent
number Y

real part of Y

4. EXPERIMENTS AND DISCUSSIONS

tanh[abs(Y )] for a complex

To show the effectiveness of the new activation function,
we conducted experiments to compare the performance of a
BSS system using the following three combinations:

To avoid constraint (7), another formula was proposed[6, 7] and has been used in [8, 9, 13, 14]:

W =  [diag(h(Y )Y H i) h(Y )Y H i] W : (8)W

Y

According to this formula,converges to a point that satisfies only (4), and the amplitudes of do not change muchduring ICA.

4.1. Conditions for experiments and overall results
Experiments were conducted for speech signals that wereconvolved with impulse responses and then mixed. Figure 4 shows the layout of the room used to record the impulse responses. The numbers of sources and observations
were both 2. The reverberation time of the room could bechanged. We used two sets of impulse responses whose reverberation times were 150 ms and 300 ms. We selected
two speech signals from the ASJ continuous speech corpus.
The lengths of the speech signals were about eight seconds,
and the entire eight seconds of the mixed data were used
for ICA. We used different frame lengths T of a windowed
DFT depending on the reverberation time as shown in Table
1. Since the sampling rate was 8000 Hz, 150 ms and 300
ms correspond to 1200 points and 2400 point, respectively.
To avoid the permutation problem [8, 9] of frequency
domain BSS, we assumed that the first source signal camefrom the left-hand side and the other came from the righthand side. Then, by using the technique of a null beamsuch that the first andformer, we set the initial value of
second rows of the matrix had steep null directivity patterns
Æ, respectively. Amplitude ambiguities
towards Æ andwere solved by the technique proposed in [8].
Table 1 shows results of SNRs for the three combinations. The numbers are the averages of SNRs at two outputs. In the column “Ref”, SNRs measured with the speechsounds used in learning are shown, and SNRs measured
with impulses [14] are shown in “Imp”. We can see that

3. A NEW ACTIVATION FUNCTION
In this section, we propose a new activation function thatsolves the problem of constraint (7) caused by activation
function (3). The new function is based on the polar coordinates of a complex number:(Y ) = tanh[  abs(Y )]  ejangle Y ;(9)where abs(Y ) and angle(Y ) are the absolute values andthe angles of Y , respectively. It consists of two parts: amplitude part tanh[  abs()] and phase part ej angle  . The
amplitude part changes the amplitudes of Y . Figure 3 shows
function values of tanh[  abs(Y )] for a complex number
Y . The phase part concerns the phases of (Y ), and main(

)

()

Y

tains the phases equal to the phases of . We can see thatthe new activation function is a natural extension of ordinary function (1), and produces the same value as (1) for areal number.
If we use this activation function, constraint (7) does
not appear. Let 
angle Yp . Since Yp is a complex
conjugate of Yp ,=() and (2) for W ,Cartesian-I using (3) for () and (2) for W ,Cartesian-diag using (3) for () and (8) for W .

Polar-I using (9) for

W

60

( )

(Yp)Yp = tanh[  abs(Yp)]  ej  abs(Yp)  e j= tanh[  abs(Yp)]  abs(Yp):
Hence, the imaginary part of h(Yp )Yp i in (5) becomes 0.

665

60

= 150

= 300

TR
ms TR
ms
Ref
Imp
Ref
Imp
Polar-I18.3 19.7
12.7 16.3
Cartesian-I
17.9 19.4
12.3 15.6
Cartesian-diag 17.8 18.0
11.9 14.6
T
TR
ms ;
TR
ms
S
,  = 0.1,  = 100, #iteration = 100
Ref: measured with speech sounds
Imp: measured with impulses

1

0.5

Imaginary= 1024 ( = 150 ) 2048 ( = 300 )
= 256

0

−0.5

W

Table 1. SNRs (dB) for different activation functions anddifferent formulas for

−1
−1

−0.5

0

0.5

1

Real

Fig. 5. Values of activation function (3)the “Polar-I” case outperforms the others, and the resultsof “Cartesian-diag” are not as good as others, especially in“Imp”. The reasons of these results are discussed in the next
subsections.

1

0.5

4.2. Comparison between “Polar-I” and “Cartesian-I”

( )

(Y )

Imaginary

In order to see how the two  functions behave, we plot
the values ofin Figs. 5 and 6. These data were obtained at the final learning step in the 200th (1554.7 Hz), whichfrequency bin for the 150 ms cases. We set 
is the parameter to control the nonlinearity of  .
In the case of (3), many samples were put at one of thej,j,j,j . On the otherfour corners:
hand, in the case of (9), the samples were put on the unit
circle ej . We can see that formula (9) represents more of
the information of , especially the phase information, thanformula (3). In addition, we can say that the entropy ofis larger in the case of (9) than in the case of (3). This
fact can be a reason why “Polar-I” outperforms “CartesianI”.
Next, we discuss the additional constraint (7), which is
imposed in the “Cartesian-I” case. Figures 7 and 8 show theHi,absolute values and the imaginary parts of I hrespectively, when using activation function (3). These data
were obtained in the 100th (773.4 Hz) frequency bin for the
150 ms cases. In Fig. 7, we see oscillations that hinder convergence. These oscillations come from the imaginary partsH
i as shown in Fig. 8.
of the diagonals of I h
If we use activation function (9), we can eliminate such
oscillations as discussed in Section 3. Figure 9 shows abH
i when using (9). We cansolute values of I h
see smooth convergence. Clearly, the mutual information
among is well minimized in this case than in the case of
(3).

= 100
( )

1+ 1

1+

−0.5

1

−1
−1

[

[

Y

[

−0.5

0

0.5

1

Real

Y

(Y )

0

Fig. 6. Values of activation function (9)erage amplitude of Yp approach some value near 1, whereas
formula (8) does not change the amplitude much. The initial value of the average amplitude of Yp differs considerably from frequency to frequency, since given signals are
not generally flat in frequency.
afNow, we discuss how the average amplitudes offects the speed of convergence in the case of “Cartesiandiag”. Sincem
m , the amount of change ofin a learning step is(Y )Y ]Y(Y )Y ]Y ( ) = WX ( )Y (m) = WX (m)=  [diag(h(Y )Y i) h(Y )Y i] Y (m).
H

(Y )Y ]YH

X h( ) i ( )Thus, the amount of change of each Yp isYp (m) = 
( )4.3. Comparison between “Polar-I” and “Cartesian-diag”

W

The main difference between the two is in the calculation of. As discussed in Section 2, formula (2) makes the av-

666

6

q =p

Yp Yq Yq m :Here, h Yp Yq i is approximately proportional to the average amplitude of Yq if the mutual information between Yp
and Yq is the same. Therefore, the average amplitude accelerates the convergence speed if it is large, and inversely, it
reduces the speed if it is small.

0.6

0.6
[1,1]
[1,2]
[2,1]
[2,2]

0.5

0.4

Absolute value

Absolute value

0.5

0.3

0.4

0.3

0.2

0.2

0.1

0.1

0
0

20

40

60

80

[1,1]
[1,2]
[2,1]
[2,2]

0
0

100

20

40

Iteration[

Fig. 7. Absolute values of I

60

80

100

Iterationh(Y )Y H i] using (3)[Fig. 9. Absolute values of I

[1,1]
[2,2]

h(Y )Y H i] using (9)

16

0.1

SNR (dB)

Imaginary part

14

0

12

10

Imp, Polar−I
Imp, Cartesian−I
Imp, Cartesian−diagRef, Polar−I
Ref, Cartesian−I
Ref, Cartesian−diag

8
−0.1

6
0

20

40

60

80

0

100

Iteration[

Fig. 8. Imaginary part of Ih(Y )Y H i] using (3)

50

100

150
Iteration

200

250

300

Fig. 10. Convergence speeds

5. CONCLUSIONS

Based on the discussion above, the convergence speed
for each frequency differs considerably in the case of (8),whereas the speed is gradually normalized in the case of (2).
Figure 10 shows convergence speeds for the 300 ms cases.
We can see that more iterations are needed to reach a pointof convergence in the case of “Cartesian-diag” than in thecase of “Polar-I” and “Cartesian-I”. To accelerate the convergence speed in “Cartesian-diag”, we could set a larger
value for step size . However, the larger  might be toolarge to converge smoothly for a frequency bin with large
energy. Therefore, some normalization of step size in eachfrequency bin would be needed for the case of “Cartesiandiag”.

We have proposed a new activation function for a complexvalued ICA algorithm. The function eliminates constraint(7), and enables us to use original gradient formula (2) without oscillations that hinder convergence. Another gradient
formula (8) has been used to avoid the oscillation problem. This formula, however, exhibits irregular convergence
speeds among frequency bins in our experiments. We obtained improved SNRs and good convergence characteristics by using the new activation function with the original
gradient formula.

667

REFERENCES[13] S. Araki, S. Makino, T. Nishikawa, and H. Saruwatari,“Fundamental limitation of frequency domain blindsource separation for convolutive mixture of speech,”in Proc. ICASSP2001, 2001, MULT-P2.3.

[1] A. J. Bell and T. J. Sejnowski, “An informationmaximization approach to blind separation and blind
deconvolution,” Neural Computation, vol. 7, no. 6,
pp. 1129–1159, 1995.

[14] R. Mukai, S. Araki, and S. Makino, “Separationand dereverberation performance of frequency domain
blind source separation for speech in a reverberant environment,” in Proc. Eurospeech2001, Sept. 2001.

[2] T. W. Lee, Independent component analysis - Theory
and applications, Kluwer academic publishers, 1998.
[3] S. Amari, A. Cichocki, and H. H. Yang, “A newlearning algorithm for blind signal separation,” inAdvances in Neural Information Processing Systems.
1996, vol. 8, pp. 757–763, The MIT Press.
[4] S. Haykin, Ed., Unsupervised adaptive filtering, John
Wiley & Sons, 2000.
[5] P. Smaragdis, “Blind separation of convolved mixturesin the frequency domain,” Neurocomputing, vol. 22,
pp. 21–34, 1998.
[6] N. Murata and S. Ikeda, “An on-line algorithm forblind source separation on speech signals,” in Proc.
Int. Symposium on Nonlinear Theory and Its Application (NOLTA ’98), 1998, pp. 923–926.
[7] S. Kurita, H. Saruwatari, S. Kajita, K. Takeda, and
F. Itakura, “Blind signal separation using directivity
pattern,” in Technical Report of Japanese Society for
Artificial Intelligence, Nov. 1999, pp. 21–26.
[8] S. Ikeda and N. Murata, “A method of ICA in time–
frequency domain,” in Proc. ICA ’99, Jan. 1999, pp.
365–370.
[9] S. Kurita, H. Saruwatari, S. Kajita, K. Takeda, and
F. Itakura, “Evaluation of blind signal separation
method using directivity pattern under reverberant
conditions,” in Proc. ICASSP2000, June 2000, pp.
3140–3143.
[10] L. Parra and C. Spence, “Convolutive blind separationof non-stationary sources,” IEEE Trans. Speech Audio
Processing, vol. 8, no. 3, pp. 320–327, May 2000.
[11] M. Z. Ikram and D. R. Morgan, “Exploring permutation inconsistency in blind separation of speech signals
in a reverberant environment,” in Proc. ICASSP2000,
2000, pp. 1041–1044.
[12] F. Asano, S. Ikeda, M. Ogawa, H. Asoh, andN. Kitawaki, “A combined approach of array processing and independent component analysis for blindseparation of acoustic signals,” in Proc. ICASSP2001,
2001, MULT-P2.1.

668