➠➡
A SPATIO-TEMPORAL FASTICA ALGORITHM FOR SEPARATING CONVOLUTIVEMIXTURES
Scott C. Douglas

Hiroshi Sawada and Shoji Makino

Department of Electrical Engineering
Southern Methodist University
Dallas, Texas 75275 USA

NTT Communication Science Laboratories
NTT Corporation
Kyoto 619-0327 Japan

ABSTRACTThis paper presents a spatio-temporal extension of the well-known
fastICA algorithm of Hyvärinen and Oja that is applicable to bothconvolutive blind source separation and multichannel blind deconvolution tasks. Our time-domain algorithm combines multichannel spatio-temporal prewhitening via multi-stage least-squares linear prediction with a fixed-point iteration involving a new adaptive
technique for imposing paraunitary constraints on the multichannel separation filter. Our technique also allows for efficient reconstruction of individual signals as observed in the sensor measurements for single-input, multiple-output (SIMO) BSS tasks. Analysis and simulations verify the utility of the proposed methods.
1. INTRODUCTION
The fastICA algorithm of Hyvärinen and Oja [1] is one of the mostwell-known and popular procedures for both independent component analysis (ICA) and blind source separation. For an m-element
linear non-Gaussian signal mixture, the procedure consists of asignal prewhitening stage followed by a set of m fixed-point iterative procedures that extract independent components using a
non-Gaussianity signal measure. Coefficient vector orthogonality
is used to guarantee uniqueness of the extracted components. The
algorithm enjoys a number of useful properties, including fast convergence, guaranteed global convergence for certain mixing conditions and contrasts, and robust behavior when noise is present.
Recently, several researchers have explored spatio-temporalextensions of spatial-only ICA and BSS procedures to attack problems in which the sources are mixed both in space and in time.
These convolutive signal mixtures appear in the design of wideband antenna arrays for wireless communications and in the design of microphone arrays for audio source separation and localization, among other applications. The simplest of these extensions treat the separation task in the (discrete) Fourier domain andapply existing spatial-only complex-valued ICA and BSS methods
within each frequency bin. These methods suffer from a numberof problems, most notably permutation, amplitude, and scaling inconsistencies across different frequency bins in the separated outputs. Although progress is being made in addressing these issues
through post-processing [2], such extensions require significant effort beyond the separation step to reconstruct the sources.
A more elegant solution is to develop convolutive BSS algorithms using a time-domain separation criterion. An example ofthis approach is the information-theoretic natural gradient convolutive BSS and multichannel blind deconvolution algorithm developed in [3]. While this procedure can be successful, the source
distributions must be approximately known, and the number of

0-7803-8874-7/05/$20.00 ©2005 IEEE

sources must be exactly known as they are simultaneously extracted. These issues make such information theoretic approaches
less practical.
In this paper, we present a spatio-temporal extension of theaforementioned fastICA algorithm to both convolutive BSS andmultichannel blind deconvolution tasks. The proposed time-domain
algorithm combines multichannel whitening via multi-stage leastsquares linear prediction with a fixed-point iteration involving a
new adaptive technique for imposing paraunitary constraints on themultichannel separation filter. A unique feature of our approachis its ability to easily and individually reconstruct the sources as
they appear in the observed signal mixtures for the single-input,
multiple-output (SIMO) BSS separation task. Analytical results
regarding the adaptive paraunitary constraint procedure are provided, and simulations indicating the usefulness of the proposed
approach for convolutive BSS are given.
2. PROBLEM FORMULATION AND EXISTING WORK
We first describe the SIMO BSS task [4]. Let si (k), 1 ≤ i ≤ mdenote m spatially-independent source signals, such that si (k) is
independent of sj (l) for i = j. These sources are measured in an
n-dimensional signal mixture with n ≥ m as
xj (k)=

νj (k) +∞m X
Xajip si (k − p)(1)i=1 p=0

for 1 ≤ j ≤ n, where {aijp } are the coefficients of the multichannel mixing system and νj (k) is uncorrelated Gaussian sensor
noise. The goal is to extract estimates of the sources as they appear
in the signal mixtures, ideally given by the mn signal set
xij (k)=∞X

ajip si (k − p)

(2)

p=0

for 1 ≤ i ≤ m and 1 ≤ j ≤ n. In practice, each xbij (k) is
estimated from a set of separated source signals yi (k) asxbij (k)=M
X

gjip yi (k − p)

(3)

p=0

where M is a filter length parameter and {gjip } must be estimatedor calculated from the separation system, the extracted signals{yi (k)}, and/or the original input signal mixtures {xj (k)}.
Previous work has identified two strategies for finding {bxij (k)}after the separated sources {yi (k)} have been found. The firstV - 165ICASSP 2005

➡

➡
PREWHITENING STAGE

SEPARATION STAGE

STfastICA

kth Prewhitening Stage

1st Prewhitening Stage

original

m separated andconvolutive(potentially)P (z)1mixtures

Q (z)
1

• • •

P (z)

Q (z)

k

W(z)

k

deconvolvedof m sources

sources

SEPARATION(copy)(copy)(copy)(copy)RECONSTRUCTION(copy)"1st" source

set all but-1

as observed in

P

1

(z)

-1

-1

Q

1

(z)

• • •

P

k

T

-1

(z)

Q

k

(z)

-1

W (z

)

1st output

original mixturesto zeroFIRST SIGNAL RECONSTRUCTION STAGE
•
•
•

•
•
•

•
•
•

•
•
•

•
•
•

(copy)(copy)(copy)(copy)(copy)"mth" source

set all but-1

as observed in

P

1

(z)

-1

-1

Q

1

(z)

• • •

P

k

T

-1

(z)

Q

k

(z)

-1

W (z

)

mth output

original mixturesto zeromTH SIGNAL RECONSTRUCTION STAGE

Fig. 1. Block diagram of the combined separation and signal reconstruction system.
strategy uses traditional linear estimation to calculate the {gijp }coefficients, where the {xij (k)} are the desired signals and the{yj (k)} are the reference signals [5]. This approach is complicated if the {yi (k)} are not uncorrelated in time, however, as
it involves m disjoint n(M + 1)-dimensional estimation tasks.
Moreover, it requires signal averaging between the {yj (k)} andthe {xi (k)}. The second calculates the inverse of the separation
system for the {gijp } [6]. This procedure is challenging due to thedifficulty of calculating a multichannel system inverse that does
not exploit a specific system structure. These procedures generally require filter lengths M that are longer than that of either the
separation system or the original mixing channel.
3. THE PROPOSED METHOD
We propose a different strategy. Consider Fig. 1, which shows asignal processing architecture containing a prewhitening stage, aseparation stage, and a signal reconstruction stage. Each of theseprocessing stages is now described.
The goal of the prewhitening stage is to decorrelate the original signal mixtures in both space and time. We propose a multistep prewhitening structure using pairs of multichannel linear systems with transfer function matrices given by
2
3
P11i (z) P12i (z) · · · P1ni (z)
.
6
7
..
6
7
0
P22i (z)(P )7 (4)Pi (z) = Di 6
6
7
.
.
..
..
4
5
0
···
0
Pnni (z)
3
2
0
···
0
Q11i (z)
..
7
6
7
6
.
Q22i (z)(Q)6 Q21i (z)7 (5)Qi (z) = Di 6
7
.
.
..
..
5
4
0
Qn1i (z)···Qnni (z)where the causal {Pi (z)} and {Qi (z)} are multichannel FIR filters of length K, the causal filters {Pjji (z)} and {Qjji (z)} have(P )(Q)
unity zero-lag coefficient values, and Di and Di are (n × n)

diagonal scaling matrices. The coefficients for the jth row of the{Pi (z)} transfer function matrices are calculated by solving aleast-squares multichannel forward linear prediction task, e.g. by
minimizing the output power of the jth output signal. The diago(P )nal entries of Di are then calculated so that the scaled forward
error residuals have unity variances. These scaled error residuals
are used as inputs to the {Qi (z)} multichannel system, in which
the jth row of the {Qi (z)} transfer function matrices are calculated by solving a second least-squares multichannel forward lin(Q)ear prediction task. The diagonal entries of Di are subsequently
calculated so that these error residuals have unity variances. Notethat this proposed method is a block-based procedure.
Several stages of this processing strategy are usually required
because the estimation of Pi (z) and Qi (z) is performed in a disjoint and sequential fashion. The exact number of stages k can
be made adaptive, with a stopping criterion that depends on how
much Pk (z) and Qk (z) differ from identity.
The above prewhitening strategy has an important advantage:
Both Pi (z) and Qi (z) can be easily inverted without calculating
any new filter coefficients, using the linear system equivalent ofbacksubstitution. Thus, so long as the linear system within theseparation stage can be easily inverted, creating the inverse of theentire prewhitening-separation system is straightforward.
The goal of the second stage is to perform separation of the
prewhitened signal mixtures based on non-Gaussianity. Here, we
propose a novel extension of the fastICA algorithm that imposes
structure on the separation system coefficients. Define the prewhitened input signal vector asv(k)vj (k)==[v1T (k) v2T (k) · · · vnT (k)]T[vj (k) vj (k − 1) · · · vj (k − L + 1)]

(6)
T

(7)

Furthermore, define
wi

=

T
T T
T
· · · win
]
[wi1
wi2

(8)

wij

=

[wij0 wij1 · · · wij(L−1) ]T(9)as the separation system coefficient vector for the ith system out-V - 166

➡

➡
Table 1: MATLAB implementation of the proposed spatio-temporal fastICA algorithm.

---------------------------------------------------------------+--------------------------------------------function [y,W] = stfica(v,L,numiter);| function [Wi] = orthW(W,n,L,numorth);|
L = L + rem((L+1),2);| i = size(W,2);[n,N] = size(v);| Wi = W(:,i)/norm(W(:,i));W = kron(eye(n),[zeros((L-1)/2,1);1;zeros((L-1)/2,1)]);| for k=1:numorth
V = zeros(n*L,N);|
Wt = zeros(n*L,1);
for i=1:n|
for j=1:i-1
V((i-1)*L+1:i*L,:) = toeplitz([v(i,1);zeros(L-1,1)],v(i,:)); |Wt = Wt + gfun(Wi,W(:,j),n,L);
end
|
end
y = zeros(n,N);|
Wi = 3/2*Wi - 1/2*gfun(Wi,Wi,n,L) - Wt;
for i=1:n| end
Wold = zeros(n*L,1);+--------------------------------------------k = 0;
| function [G,C] = gfun(U,V,n,L);y(i,:) = W(:,i)’*V;|
crit = 1;| Wi = zeros(L,n);Wi(:) = U;
while (crit*(k<numiter))| Wj = zeros(L,n);Wj(:) = V;
Wold = W(:,i);| Ct = zeros((3*L-1)/2,1);
k = k+1;| Z = zeros((L-1)/2,1);
f = y(i,:)’.ˆ3;% OR f = tanh(20*y(i,:)’);| ll = (L+1)/2:(3*L-1)/2;
llr = L:-1:1;
fp = 3*sum(y(i,:).ˆ2); % OR fp = 20*sum(sech(20*y(i,:))); | for i=1:nW(:,i) = V*f - fp*W(:,i);|Ct = Ct + filter(Wi(llr,i),1,[Wj(:,i);Z]);
W(:,i) = orthW(W(:,1:i),n,L,10);| endy(i,:) = W(:,i)’*V;| C = Ct(ll);crit = (abs(abs(W(:,i)’*Wold)-1)>0.0001);| Gt = filter(C(llr),1,[Wj;zeros((L-1)/2,n)]);end| Gt = Gt(ll,:);end| G = Gt(:);---------------------------------------------------------------+---------------------------------------------

put. Then, we compute the ith separated signal sequence asyi (k)=

wiT v(k)(10)

for 1 ≤ k ≤ N , assuming a data record length of N samples.
The coefficient vectors {wi } are updated iteratively using amodified version of the fastICA algorithm as follows:• Step 1: Compute yi (k) in (10) for 1 ≤ k ≤ N .
• Step 2: Update the coefficient vector as
1 X
f (yi (k))v(k) − f  (yi (k))wi ,
N
N

wi

←−

(11)

k=1

Then, (13)–(15) can be expressed using polynomials as» hi(L−1)/2
3
1
W i (z) ←− W i (z)−W Ti (z −1 )W i (z)W i (z)22−(L−1)/2#L−1
i−1 hi(L−1)/2X
W Tj (z −1 )W i (z)W j (z ) (17)+−(L−1)/2hi(L−1)/2
W Tj (z −1 )W i (z)• Step 3: Normalize the length of the coefficient vector as
←−

q

wi

−(L−1)/2(12)

wiT wi•Step 4: While wi is not paraunitary with w1 , w2 , . . ., wi−1 ,
wi

←−

i−1
X
3
1
wi − g(wi , wi ) −g(wi , wj )
2
2
j=1

=

(14)

The first three steps of this procedure are identical in form to
the single-unit fastICA algorithm. To understand Step 4, define the(n × 1)-dimensional ith system vector polynomial as
W i (z)=0

for 1 ≤ j < i.(19)

wi

=

Wi

=

W i (z)|z=ejωˆ˜˛
W 1 (z) W 2 (z) · · · W i−1 (z) ˛z=ejω

(20)
(21)

Then, (17) can be rewritten as

Cij wjk

and the (p, q)th element of Cij is given by
8
m L−1
>
L−1
< XX
wjkl wik(l+p−q) if |p − q| <[Cij ]pq =2 (15)>: k=1 l=00otherwise.

L−1
X

−→

The above constraints are a spatio-temporal extension of the orthonormality constraints imposed on wi in the original fastICA
procedure and imply that the separation system is paraunitary.
As further justification of the iterative procedure for enforcing
paraunitary constraints, let L → ∞, and define

(13)

T
T
T
where g(wi , wj ) = [gij1
gij2· · · gijm
]T withgijk

0

where [·]KJ denotes truncating the polynomials of its argument toorder −J through −K. Extensive simulations of this iterative subprocedure indicate that (17) causeshi(L−1)/2
W Ti (z −1 )W i (z)
−→ 1
(18)

where f (yi ) is the contrast nonlinearity.
wi−(L−1)/2

j=1

T[wi1l wi2l · · · winl ] z−l.

wi,new= 32 wi − 12 ||wi ||2 wi − Wi WH
i wi

(22)

where ·H denotes Hermitian transpose. Define the variables
a0i
a1i
2

=

||wi ||2 − 1

(23)

=

2
||W H
i wi ||(24)wH
i wi .

The condition |a0i | = a1i = 0 imwhere ||w i || =plies that [W 1 (z) · · · W i (z)] form an m-dimensional paraunitary
sequence if [W 1 (z) · · · W i−1 (z)] is already paraunitary. Then,(22) implies that these state variables evolve as

(16)

j=0

V - 167

a0i,new=

a1i,new=

1 2
a0i (a0i − 3) + (a0i − 1)a1i
4
1 2
a a1i .
4 0i

(25)
(26)

➡

➠
0

0

10

4. SIMULATIONS

10
1st Channel
2nd Channel
3rd Channel
4th Channel
5th Channel

−5

10

−5

10

−10

10

−10

−15

10

10

−15

E{a1i}

E{a20i}

10

−20

10

−20

−25

10

−30

10

10

−25

10

−30

10

−35

10

−35

0

2

4
6
number of iterations

10

8

0

2

4
6
number of iterations

82cFig. 2: Evolution of E{ac1i } for (13)–(15).
0i } and E{aThis pair of nonlinear coupled scalar equations can be easily simulated for different initial conditions, and such simulation studies
show that a0i and a1i converge to zero for a wide range of initialvalue pairs. Empirically, we have observed convergence of thissystem if−1 < a0i < 2
a1i < a0i + 1

or

0 < ||wi ||2 < 3
22||WHi wi || < ||w i || ,(27)

which are typically satisfied in practice.
Table 1 provides MATLAB code for implementing the spatiotemporal fastICA algorithm complete with appropriately chosen
stopping criteria, where v is the (n × N ) prewhitened signal matrix, L is the separation system filter length, and numiter is the
maximum number of iterations of the fastICA routine.
The goal of the last set of parallel stages, shown at the bottomof Fig. 1, is to reconstruct the individual sources as they appear
in the original mixtures. The reconstruction of the ith separated
signal involves setting all but the ith signal yi (k) to zero and thenpassing this signal through the inverse of the prewhitening and separation systems. For this calculation, note that
[W(z)Qk (z)Pk (z) · · · Q1 (z)P1 (z)]−1
=

−1
−1
−1
T
−1
P−1
) (28)
1 (z)Q1 (z) · · · Pk (z)Qk (z)W (zdue to the paraunitariness of W(z) as constructed by the separation stage. Because of the triangular structures of the Pi (z) andQi (z) systems, they can be easily inverted.
The method we have described has a number of advantages
over competing approaches:
1. No step size needs to be selected.
2. Knowledge of the source distributions is not needed, so long as
their statistics imply a non-zero contrast value.
3. The number of non-Gaussian sources within the mixture neednot be known a priori.
4. For SIMO BSS, the system inverse used for signal reconstruction is computed directly and exactly.
5. Convergence of the newly-developed spatio-temporal fastICA
procedure appears to be as fast as its spatial-only counterpart; usually fewer than 10 iterations per unit are needed for i.i.d. sources.

We now present numerical evaluations indicating the performanceof the proposed methods. We first investigate the iterative paraunitary constraint scheme in (13)–(15). For these evaluations, wehave chosen n = 10 and L = 51. For each simulation run, wi was
initialized to a 510-element vector containing zero-mean uncorrelated Gaussian noise of variance 10−4 summed with a single nonzero unity-valued “center” tap at position k = L(i−1)+(L+1)/2.
Shown in Fig. 2 for 1 ≤ i ≤ 5 are the average evolutions of2E{acc1i }, computed from the elements of Cij for 1 ≤0i } and E{aj ≤ i, as averaged over 100 different simulation runs. As can be
2seen, convergence to a paraunitary condition given by E{ac0i } ≈ 0
}
≈
0
is
fast,
approachingthemachineprecisionofand E{ac1i
MATLAB in about 10 iterations.
We now demonstrate the ability of the proposed method tosolve the SIMO BSS task. The measurements used for this evaluation were taken from a two-loudspeaker, two-microphone acoustic
laboratory setup similar to that described in [7], in which the room
reverberation time was approximately 130 msec. Mixtures containing a male and a female talker sampled at 8kHz were processedusing the technique in Fig. 1, where k = 4, K = 101, L = 81,
N = 55200, and numiter= 100. After processing, the meansquared errors between the reconstructed sources and the original
mixtures for the two-channel signals were between 10 dB and 11.8
dB below the powers of the original source mixtures for both talkers in the four output channels. In addition to these results, wehave verified that our proposed method successfully deconvolves
convolutive mixtures of i.i.d. signals with binary, uniform, andLaplacian distributions; these results are omitted due to space limitations.
5. CONCLUSIONSIn this paper, we have described a spatio-temporal extension of
the well-known fastICA procedure. Our method employs leastsquares prewhitening methods and a novel iterative scheme formaintaining paraunitary constraints on the separation system. The
procedure enables the reconstruction of the individual sources within the signal mixtures for single-input, multiple-output BSS without explicitly calculating the inverse of the separation system’s
impulse response. A theoretical study of the proposed method’s
capabilities is the subject of current work.
6. REFERENCES[1] A. Hyvärinen and E. Oja, “A fast fixed-point algorithm for independent component analysis,” Neural Computation, vol. 9, no. 7, pp. 1483-1492, Oct. 1997.
[2] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust and precise method
for solving the permutation problem of frequency-domain blind source separation,” IEEE Trans. Speech Audio Processing, vol. 12, pp. 530-538, Sept. 2004.
[3] S. Amari, S.C. Douglas, A. Cichocki, and H.H. Yang, “Multichannel blind deconvolution and equalization using the natural gradient,” Proc. 1st IEEE Workshop Signal Process. Adv. Wireless Comm., Paris, France, pp. 101–104, Apr.
1997.
[4] T. Takatani, T. Nishikawa, H. Saruwatari, and K. Shikano, “SIMO-model-based
independent component analysis for high-fidelity blind separation of acoustic
signals,” Proc. 4th Int. Symp. Indep. Compon. Anal. Signal Separation, Nara,
Japan, pp. 993–998, Apr. 2003.
[5] Y. Inouye and K. Tanabe, “Super-exponential algorithms for multichannel blind
deconvolution,” IEEE Trans. Signal Processing, vol. 48, pp. 881-888, Mar.
2000.
[6] S. Ikeda and N. Murata, “A method of ICA in time-frequency domain,” Proc.
Int. Workshop Indep. Compon. Anal. Blind Signal Separation, Aussois, France,
pp. 365-371, Jan. 1999.
[7] S.C. Douglas, H. Sawada, and S. Makino, “Natural gradient multichannel blind
deconvolution and speech separation using causal FIR filters,” IEEE Trans.
Speech Audio Processing, vol. 13, pp. 92-104, Jan. 2005.

V - 168