Comparison of Blind Source Separation Methods based on

Iterative Algorithms
J. Rinas, K.D. Kammeyer
Department of Communications Engineering, Universität Bremen, {rinas,kammeyer}@ant.uni-bremen.de

Abstract
In this paper some approaches to detect signal streams of a multi layer transmission system are presented. We will
focus on blind algorithms for the separation of the data stream and improve their performance in an iterative way
in order to gain nearly the same performance as with a known channel matrix. The overall algorithm will remain
blind and does not need any training data.

1 Introduction transmitter receiver

Multi layer transmissions are an up to date technology
to exploit the capacity of a spatial channel. Lately, s0 x0

detector
powerful detection algorithms have been developed. data H x1
data
These algorithms need a suitable estimation of the
transmission channel. In most cases the channel estima-
tion is obtained by including a pilot sequence in the data
stream. However, this will lead to a loss of the payload snT −1 xnR −1
rate that can be transmitted. Therefore we will introduce ν
a new scheme that combines blind source separation
Fig. 1. transmission system
techniques with an efficient detection algorithm in an
iterative way. This will lead to an overall blind detection
of the transmitted symbols. We will study this idea seen in figure 1. This system can be described by
with some variants of the constant modulus algorithm
(CMA) and higher order statistics (HOS) based source x=H·s+ν . (1)
separation approaches. Equation (1) models a transmission of signal
The remainder of the paper is organized as follows: streams s trough a spatial channel H. Let s =
In section 2 we will introduce the transmission system T
[s0 (k), s1 (k), . . . , snT −1 (k)] be the vector of trans-
that was used for all simulations. In section 3 we will mitted symbols at time instant k. We usually use a
present some approaches to achieve a blind separation block processing of the data and use a length of k =
of a multi layer transmission. We will present some 0 . . . L − 1 symbols. For a clear arrangement and short
variants of multiple input multiple output (MIMO) formulas we‘ll drop the index k in the following text.
CMA algorithms and point out connections with clas- The dimensions of H are nR × nT for nR receive
sical blind source separation (BSS) approaches. We and nT transmit antennas.
compare their performance using Monte Carlo simula- ⎡ ⎤
tions. A new approach to combine source separation h0,0 ··· h0,nT −1
⎢ .. .. .. ⎥
techniques with a multi layer detection scheme is H=⎣ . . . ⎦ (2)
presented in section 4. This will lead to an overall blind hnR −1,0 ··· hnR −1,nT −1
symbol detection scheme. Section 5 introduces channel
coding into the presented scheme and illustrates some The entries of H are assumed to be Gaussian dis-
problems in this context. A summary and concluding tributed and statistical independent. They are assumed
remarks can be found in section 6. to be constant for one block of L time instants
in order to model a slow fading channel. ν :=
T
[ν0 (k), ν1 (k), . . . , νnR −1 (k)] is a vector with additive
2 System description Gaussian noise that is assumed to be white.
During this paper we concentrate on a scenario with
In this paper we analyze a multi layer transmission as nT = 4 transmit and nR = 4 receive antennas. We
This work was supported by the German national science founda- transmit L = 200 uncoded and coded QPSK symbols.
tion (DFG) within the project AKOM (project #Ka 841/9-1) These signals have a constant magnitude if the symbol


timing is perfectly known and therefore can be detected cum ei , ei , ej , el , obtained from the extracted signals
by constant modulus approaches. The energy of the e defined in equation (6).
transmitted signals is normalized so that the average n
T −1


received energy of one signal is one. (Therefore we’ll max =
! cum ei , ei , ej , el , 2 (5)
only observe diversity gains and no additional array B
i,j,l=0
gain due to an increasing number of antennas.)
This optimization problem is solved by an eigenvalue
decomposition of Qz and a joint diagonalization of
3 Blind separation approaches for the dominant eigenvectors rearranged as matrices. This
diagonalization leads to the unitary nT × nT matrix B
communication signals and the independent data streams
There are several methods available to separate data- e = BH · z . (6)
streams blindly with more or less knowledge of the
signals. In this paper we concentrate on two families of A similar approach that needs less computational effort
approaches. The classical approaches try to maximize is the SSARS algorithm presented in [4].
the statistical independence of signal streams and make
no use of the discrete alphabet of modulated signals,
whereas the constant modulus approaches considered 3.1.2 fastICA
here use the discrete amplitude level of the QPSK Beside the JADE algorithm the fastICA algorithm is
signals. organized in a different way [5] [6]. The basic idea
of this algorithm is to do a blind source extraction
3.1 Classical blind source separation (BSE) for each component and to prevent the same
The classical blind separation schemes are based on signal from being extracted multiple times. It starts with
the assumption of independent components in the the whitening of the received data in the same way as
received signal streams. This property is fulfilled if presented before in (4).
the data carried within the symbols of the received In order to extract signal number t out of the mixture
streams is random. Established algorithms that lead to a z (4) a randomly initialized extraction vector bt – a col-
separation are the JADE algorithm as a batch algorithm umn vector of the nT × nT matrix B – is generated. In
and the fastICA algorithm that solves the problem order to preserve the signals to remain uncorrelated B
by extracting the independent components step by step. has to be a unitary matrix. Therefore bt is constrained
to form an orthonormal basis using the knowledge of
the vectors b0 . . . bt−1 obtained in former iterations.
3.1.1 JADE To achieve this goal matrix
The joint approximate diagonalization of eigenmatrices Bt = [b0 , b1 , . . . bt−1 ] (7)
(JADE) algorithm [1] is a batch procedure that solves
the separation problem. containing the extraction vectors of the former itera-
The first step of this algorithm is to decorrelate the tions is build. The randomly initialized vector bt is
input streams x. That is the nT × nR matrix W has to orthogonalized to the former detected ones
be calculated to fulfill
bt = bt − Bt BH
t bt (8)
H 1
I = WRxx W (3)
and normalized to a length of one
 H
with Rxx = E xx being the spatial correlation bt = bt / ||bt || . (9)
matrix of the received signals. In the source separation
literature [2] [3] this step is also known as principle In order to determine bt we choose the maximization
component analysis and can be solved using a eigen- of the kurtosis of a single signal as the criterion.
value decomposition of Rxx . The decorrelation can be  
4
computed by a multiplication of the received signal max JfastICA,t (et ) = max E |et |
streams x with the whitening matrix W.
bt bt  4  (10)
= max E bH t z

bt
z = Wx (4)
This can be solved using a fixed point iteration in-
The utilization of second order information is not suffi- cluding the additional constraints (8) and (9) [6]. The
cient to obtain independent signals. Therefore the JADE resulting signal streams e can be extracted by multiply-
algorithm additionally uses 4th order information. It ing the received signal with the matrix of all collected
maximizes some elements of the cumulant matrix Qz = extraction vectors B = BnT −1 .
1 WH denotes the conjugate transpose i.e. the Hermitean of W. e = BH · z (11)
3.2 MIMO CMA off diagonal elements to the general MIMO CMA cost
function (12). This leads to the new cost function
The idea of the constant modulus algorithm (CMA)
is to compare the amplitude of the equalized received nT−1

(i) 2
signal with a reference amplitude [7][8][9]. This com- Jcorr (C) = JCMA (C) + ψk,l (16)
parison expressed in mathematical terms leads to a cost k,l=0; k=l
function that has to be minimized. The general cost
where
function for MIMO CMA is  
 2  (i) (i) (i) ∗

nT −1
2
ψk,l = E ek · el (17)
JCMA (C) = E |et | − 1
t=0  2 

nT −1 H 2 describe the correlation of the extracted signals ek and
(i)
= E c t x − 1 , (i)
t=0 el as elements of a matrix
(12)
⎛ ⎞
where et denotes the equalized component number t 0
(i)
ψ0,1 ...
(i)
ψ0,nT −1
of the output vector e. The component ct of the nR × ⎜ (i) (i) ⎟
⎜ ψ1,0 0 ... ψ1,nT −1 ⎟
nT CMA extraction matrix C = [c0 , c1 , . . . , cnT −1 ] Ψ(i) =⎜ ⎟
corr ⎜ .. .. .. ⎟
spatially filters the signal et . Separation of all signals ⎝ . . . ⎠
can be done by (i) (i)
ψnT −1,0 ψnT −1,1 ... 0 .
e = CH · x . (13) (18)
Note that (18) contains no diagonal elements.
In order to minimize the CMA cost function (12) the
Using this cost function leads to the update equation
steepest descent algorithm is used.
∂   C(i+1) C(i)
C(i+1) = C(i) − μ J C (i) 2
(14)  =  
∂C(i)
CMA  (i) (i) H (i) H
−μ 4E D e − I + Ψcorr ·x·e .
μ is the step factor of the gradient descent. The CMA
(19)
algorithm is initialized by a whitening matrix C0 = W
as described in equation (4) . This starts the CMA
algorithm with separation vectors pointing in directions
of uncorrelated sources of power and improves conver- The computational effort of this approach can be
gence. immense, because the correlation matrix Ψicorr has to
The Matrix C can be calculated using the update be estimated in every iteration.
equation
    H

C(i+1) = C(i) − 4μE D e(i) − I · x · e(i) . 3.2.2 2nd approach: determinant penalty

 (i) (15) Another idea to prevent the CMA algorithm to extract

where the operator  D e
forms the diagonal matrix
the same signal twice can be found in [11]. The
 (i) (i) 2 (i) 2 (i) 2 basic idea of this approach is to minimize the linear
D e = diag e0 , e1 , . . . , enT −1 . The
dependency of the separation vectors contained in C.
main problem with this algorithm is, that it is still The linear dependency is measured here in terms of
possible to resolve the same signal multiple times. This the logarithm of the magnitude of the determinant of
will lead to very bad error rates of the algorithm as can C. This leads to the cost function
be seen in figure 3.

Therefore we will present three approaches to avoid Jdet (C) = JCMA (C) − ln det CH . (20)
this problem. Two of them are already known from
different sources. A new approach presented here has For nearly linear dependent components in C the right
a slight connection to the approach of the fastICA hand side term of (20) will grow to very high values.
algorithm. The corresponding update equation can be calculated
to be
3.2.1 1st approach: correlation penalty C(i+1)
 =C
(i)

 (i)    −1 
One approach to tackle the problem of extracting the (i) H (i) H
−μ 4E D e −I ·x·e − C .
same signal multiple times was presented in [10]. The
idea of this approach is to minimize the correlation (21)
of the extracted signals. This can be realized by cal-
culating the cross correlations of all signals in every This drawback from the computational side of view is
iteration step and adding the squared magnitude of the that this update equation contains a matrix inversion
in every iteration.
2 The superscript (i) denotes the iteration step.
3.2.3 3rd approach: subspace limitation and the permutation problem. In this paper we neglect
this problem and calculate a propper estimation of the
Motivated by the CMA approaches presented above and
permutation matrix P̂ by using all transmitted data.
the idea to successively separate signals we present a
(The quadrant ambiguity problem can be solved by
new scheme to separate constant modulus signals from
e.g. using differential coding, while the permutation
a MIMO system.
problem can be used by including addresses that can
As a first step we decorrelate the received signal
also be used by higher layers.)
streams as described by equation (4). In order to extract
the signal number t from our whitened stream z we 
start with a nR × 1 vector ct that is zero except for P̂ = E s̃ · sH (24)
a one at row t. This vector is orthogonalized to the
former vectors as done by the fastICA algorithm in Using P̂ we get the assignments by successively tak-
equations (8) and (9). ing the maximum absolute values. The corresponding
Then we adjust ct by the steepest descent algorithm phase factors can be used to determine the quadrant
using the CMA cost function (12) trimmed down to ambiguity.
update only ct .
This leads to the update equation
  
(i+1) (i) (i) 2 (i) H
ct = ct − 4μE et − 1 · z · et . −1
10
(22)
After separating one signal we consider the known
separation vector using equations (8) and (9) to select −2
10
the next initialization in order to gain the next signal.
BER

Compared to the update equation of the general

MIMO CMA (15) we need no additional computations
inside the iteration loop. 10
−3
CMA
CMA corr
CMA det
3.3 performance comparison CMA sub
zero forcing
All separation methods presented so far will result in 10
−4
MMSE

separated data streams e. In order to compare their 0 5 10 15 20

Eb /N0 in dB
25 30 35 40

performance in terms of bit error rates (BER) we use

a set up as depicted in figure 2. Fig. 3. BER error rates of all CMA approaches

phase Figure 3 depicts the simulation results of the three

source separation

x0 e0 correct. CMA approaches as well as two reference curves. One

s̃0
blind

reference curve is the linear zero forcing filter and

x1
one is the MMSE solution with ideally known channel
matrix to equalize the spatial mixture. For the CMA
phase approaches we use a step factor of μ = 0.01 and used
xnR −1 enT −1 correct. s̃nT −1 300 iterations (not optimized).
From figure 3 (CMA curve) we can state that
Fig. 2. source separation set up for BER measurements it is absolutely necessary to prevent the CMA from
detecting one signal multiple times. The other 3 CMA
Since the phase of every data stream et has not approaches work well but suffer from an error floor
been taken into account during the separation, every due to the steepest descent algorithm; We reach BER
data stream is rotated by a random phase factor and regions where an application of channel coding should
permutated. be possible. As expected the performance of the al-
et,derot = et · e−(arg E{et })/4
4 gorithms is between the zero forcing and the MMSE
(23)
solution for the equalization of the spatial channel.
Equation (23) points out an approach feasible for QPSK The new subspace limitation approach has the best
signals. The average of the data signal taken to the performance of all CMAs because it includes a very
power of four will remove the data. The result can be hard constraint to avoid same signals.
taken to estimate the phase rotation whereby a discrete Figure 4 shows bit error curves of the classical
phase ambiguity (quadrant ambiguity) will remain. We source separation algorithms. They gain comparable
can take the derotated signals to decide the symbol performance to the CMA with subspace limitation, but
positions s̃t . In order to do the bit demapping of the do not show an error floor, because these algorithms
datastreams we have to solve the quadrant ambiguity don’t use a steepest descent.
 Since this detection loop is initialized in a blind way
and the VBLAST algorithm only decides symbol posi-
−1
10 tions the whole detection scheme remains blind.
Because of quadrant ambiguities of the initially
used symbols s̃0 . . . s̃nT −1 we have to prove that the
10
−2 proposed scheme works. Therefore we define symbol
BER

matrices from the vectors used so far. The transmitted

data streams are arranged in matrix S

10
−3 S = [s(0), s(1), . . . , s(L − 1)] (25)
fastICA and the received and coarsely decided data streams in
JADE
zero forcing
matrix S̃
MMSE
10
−4
0 5 10 15 20 25 30 35 40
S̃ = [s̃(0), s̃(1), . . . , s̃(L − 1)] . (26)
Eb /N0 in dB
Fig. 4. BER error rates of the source separation approaches using S̃ ≈ ΦPS (27)
HOS
We can formulate the quadrant ambiguity of the
decision by equation (27). Whereby P is a random
4 Application of Iteration Tech- permutation matrix and Φ = diag [φ0 , φ1 , . . . φnT −1 ]
is a diagonal matrix modelling the discrete quadrant
niques - uncoded ambiguities. The equal sign is only valid, if there are
The algorithms presented until this point only approxi- no decision errors.
mated linear spatial filters in order to separate the data In order to use the coarse decision results S̃ for
streams. In this section we try to improve the detection channel estimation we calculate
performance by applying a cancellation scheme. This
H̃ = X · S̃+ 3
. (28)
will utilize the finite symbol alphabet that was only
used marginally till now. Assuming no noise in the system will lead to
H̃ = HS · S̃+
initialization iteration  −1 (29)
s̃0 = HS · S̃H S̃S̃H .
s̃1
s̃nT −1 If we introduce equation (27) and assume correct deci-
sions4 we get

−1
H̃ = HS · SH PH ΦH SPΦΦH PH SH . (30)
channel
estimation Because Φ and P are unity matrices (ΦΦH = I and
PPH = I) we can simplify the expression to

−1
H̃(i) H̃ = HS · SH PH ΦH SSH . (31)
(i)
s̃0 If we assume that S contains uncorrelated signal
x0 VBLAST (i) streams of sufficient length L, the terms S · SH will
x1 (MMSE) s̃1
be approximately diagonal matrices. Therefore we can
xnR −1 (i)
s̃nT −1 further simplify to
Fig. 5. improving the channel estimation / detection H̃ = H · PH ΦH . (32)

We will use a system as depicted in figure 5. We start This leads to an estimation of the channel matrix H̃
with coarsely decided symbols s̃0 . . . s̃nT −1 that were in a permutated form where every column contains a
obtained using a blind separation method as depicted quadrant error. If we apply this channel estimation in
(i) (i)
in figure 2. Using this data we produce a first channel the VBLAST algorithm we get symbols s̃0 , . . . s̃nT −1
estimation H̃(0) . This channel estimation will be used (figure 5) with the corresponding discrete phase ambi-
to detect the symbols once more using the VBLAST guities φ0 , . . . , φnT −1 , but this will not influence our
detection algorithm as presented in [12]. We use the further detection and cancellation process, as long as
MMSE variant of the VBLAST algorithm but also other we only want to decide the symbol positions.
multi layer detection schemes are possible e.g. SQRD To summarize: We found an iterative estimation and
[13] [14]. detection scheme that utilizes the finite symbol alphabet
Using the output of the VBLAST detector for improved and remains completely blind.
channel estimation in combination with a new detection 3 The + sign denotes the Moore-Penrose pseudo inverse.
of the data will iteratively lead to better results. 4 The terms in (27) are equal.
4.1 performance of the proposed iterative transmitter receiver
scheme FEC
In order to show the feasibility of our detection ap- s0 x0

detector
proach we present some BER results of the transmis- data H x1
data
sion system presented in figure 1. As an initialization
we exemplary use the output of the JADE algorithm.
The permutation problem was solved in the same way
FEC
as above (24). snT −1 xnR −1
ν
JADE
JADE it=5 Fig. 7. transmission system with channel coding
zero forcing
10
−1 MMSE
VBLAST

code, because a phase rotated QPSK signal will lead to

10
−2 an invalid codeword. Therefore we introduce a further
BER

phase comparison with the transmitted data (after the

hard decision of s̃0 . . . s̃nT −1 ) in order to make the
Viterbi decoding process possible.
−3
10
The decision device inside the iterations using the
VBLAST algorithm can be replaced by a Viterbi de-
tector without any problems.
−4
10
0 5 10 15 20 25 30 35 40
Eb /N0 in dB
JADE
JADE it=5
Fig. 6. BER of JADE with iterative VBLAST detection cJADE
10
−1 cJADE it=5
cVBLAST
Figure 6 depicts the results of our simulations.
Beside the reference curves of the MMSE and zero
forcing detection the raw VBLAST detection with ideal −2
10
known channel matrix was introduced as a further
BER

reference.
Using our new iterative scheme (figure 5) we can
observe a gain of about 10 dB at a bit error rate of 10−3 10
−3

(JADE it=5) compared to the classical source separation

using only the JADE algorithm. For this improvement
we need only 5 iterations of detections and channel −4
10
estimations. We nearly reach the curve of the VBLAST 0 5 10 15 20 25 30 35 40

algorithm with known channel matrix H .5 We have Eb /N0 in dB

to emphasize that the whole detection scheme remains Fig. 8. BER or coded multi layer transmission
blind since no reference data is used to gain the symbol
decisions.
Figure 8 depicts the results of a Monte Carlo sim-
ulation using such a set up. We can find only a small
5 Inclusion of channel coding coding gain between the uncoded (JADE) and the coded
(cJADE) scheme. This is due to the hard decision input
In order to improve the overall robustness against
of the Viterbi decoder. When using VBLAST iterations
noise and to improve the iterative estimation loop we
in order to improve the channel estimation, we can only
introduce channel coding in our system. For channel
see a coding gain at low SNRs. This is a common
coding we use the (5, 7)8 convolutional code with tail
behavior of coded multi layer systems with VBLAST.
bits in each stream as depicted in figure 7.
Since the SNR gets better, the interference cancellation
On the receiving side we replace every decision
process inside the VBLAST can’t be improved by a
device with a Viterbi decoder. But a simple replacement
coder. Based on this first analysis we can state that
will lead to a problem: The input signal of the decision
it is a promising approach to introduce coding in the
device in figure 2 has to cope with a quadrant ambigu-
system. But we have to cope the problem of quadrant
ity of the signal. This can not work with a convolutional
ambiguity of the incoming signal in the input of the
5 We can even decrease the gap to the detection with ideally known decision devices. Therefore further work will be done
channel matrix if we increase the length of the data block L. on tackling this problem.
6 Summary and Conclusions [12] P. W. Wolniansky, G. J. Foschini, G. D. Golden, and R. A.
Valenzuela, “V-BLAST: An Architecture for Realizing Very
High Data Rates Over the Rich-Scattering Wireless Channel,”
Blind source separations techniques are promising in invited paper, Proc. ISSSE-98, 29. September 1998.
connection with MIMO communication systems. Their [13] D. Wübben, R. Böhnke, J. Rinas, V. Kühn, and K. Kammeyer,
“Efficient algorithm for decoding layered space-time codes,”
strength is that they can extract signal streams utilizing IEE Electronic Letters, vol. 37, no. 22, pp. 1348–1350, Novem-
only some statistical properties of the signals. They do ber 2001.
not need any reference data. [14] D. Wübben, J. Rinas, R. Böhnke, V. Kühn, and K. Kammeyer,
“Efficient Algorithm for Detecting Layered Space-Time Codes,”
In this paper we presented in detail some separation in 4th International ITG Conference on Source and Channel
methods that are based on the MIMO constant modulus Coding, Berlin, Germany, January 2002.
property. We presented a new computational efficient [15] A. Stuard and K. Ord, Kendalls Advanced Theory of Statistics
- Volume I: Distribution Theory, 6th ed. New York: Wiley,
CMA approach using a subspace limitation for MIMO 2001.
systems.
We used the blind separation methods to initialize
MIMO detection schemes and to improve channel esti-
mation and data detection in an iterative way. Therefore
we pointed out that these iterations can be done without
reference data, so that the overall symbol detection
algorithm remains blind. We have shown that we gain
about 10 dB compared to the classical separation
approaches and reach nearly the performance of the
VBLAST detection algorithm with known channel ma-
trix.
We presented a first step to combine the proposed
iterative scheme with channel coding and pointed out
directions to further improvements.

References
[1] J.-F. Cardoso and A. Souloumiac, “Blind beamforming for non
Gaussian signals,” IEE Proceedings-F, vol. 140, no. 6, pp. 362–
370, December 1993.
[2] A. Cichocki and S. Amari, Adaptive Blind Signal and Image
Processing. Wiley, 2002.
[3] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component
Analysis. New York: Wiley, 2001.
[4] M. Feng and K. Kammeyer, “Blind Source Separation for
Communication Signals Using Antenna Arrays,” in Proc. IEEE
Int. Conference on Personal Wireless Communications (ICUPC-
98), Florence, Italy, October 1998.
[5] A. Hyvärinen and E. Oja., “A fast fixed-point algorithm for
independent component analysis,” Neural Computation, vol. 9,
no. 7, pp. 1483–1492, 1997.
[6] E. Bingham and A. Hyvärinen, “A fast fixed-point algorithm for
independent component analysis of complex-valued signals,”
Int. J. of Neural Systems, vol. 10, no. 1, pp. 1–8, 2000.
[7] K. D. Kammeyer, Nachrichtenübertragung, 2nd ed. Stuttgart:
Teubner, 1996.
[8] R. Mann Pelz, W. Tobergte, and K. Kammeyer, “Applications
of Constant Modulus Algorithms for Adaptive Equalization
of Time-Varying Multipath Channels,” in International Con-
ference on Digital Signal Processing, Florence, Italy, January
1987, pp. 421–425.
[9] R. Mann Pelz and K. Kammeyer, “A Pole-Zero-Tracking
Constant Modulus Algorithm,” in Proc. IEEE Int. Conf. on
Acoustics, Speech and Signal Proc., Glasgow, Great Britain,
January 1989, pp. 1227–1230.
[10] C. Papadias and A. Paulraj, “A constant modulus algorithm for
multiuser signal separation in presence of delay spread using
antenna arrays,” IEEE Signal Processing Letters, vol. 4, no. 6,
pp. 178–181, June 1997.
[11] C. Mejuto, A. Dapena, and L. Castedo, “Frequency-Domain
Infomax for Blind Separation of Convolutive Mixtures,” in
Second International Workshop on Independent Component
Analysis and Blind Signal Separation, June 2000, pp. 315–320.