# An Analytical Multimodulus Algorithm for Blind Demodulation in a Time-Varying MIMO Channel Context

Authors: Steredenn Daumont; Daniel Le Guennec
Journal: International Journal of Digital Multimedia Broadcasting
 (2010)

---

## Abstract

This paper addresses the issue of blind multiple-input multiple-output (MIMO) demodulation
of communication signals, with time-varying channels and in an interception
context. A new adaptive-blind source separation algorithm, which is based on the implementation
of the Multimodulus cost function by analytical methods, is proposed. First a
batch processing analysis is performed; then an adaptive implementation of the
(Analytical Multi-Modulus Algorithm) AMMA and its simplified version named  (Analytical
Simplified Constant Modulus Algorithm) ASCMA is detailed. These algorithms, named
adaptive-AMMA and adaptive-ASCMA, respectively, are compared with the adaptive  (Analytical Constant Modulus Algorithm) ACMA and the MMA (Multi-Modulus Algorithm).
The adaptive-AMMA and adaptive-ASCMA achieve a lower residual intersymbol
interference and bit error rate than those of the adaptive-ACMA and MMA.


---

## Body

## 1. Introduction

The context of this study is the Multiple-Input Multiple-Output (MIMO) interception, and more specifically blind demodulation in reception. The channel is supposed to be not stationary, so the methods used to demodulate the received signals must be fast. Such systems exist, among them are the Blind Source Separation (BSS) implemented with analytical methods.BSS algorithms require few information about emission signals in order to process and recover the transmitted symbols (sources). This type of algorithms estimate an equalization matrix from received signals to obtain the transmitted sources. To find this matrix, different cost function can be minimized, like the Multimodulus (MM) [1, 2], the Constant Modulus (CM) [3, 4], the Simplified Constant Modulus (SCM) [5], and the Multiuser Kurtosis (MUK) [6]. The algorithms in [1–6] use the stochastic gradient to minimize these cost functions, so they are slow to converge and with time-varying environments they are less accurate. So, analytical methods such as the adaptive Analytical Constant Modulus Algorithm (ACMA) [7–10], which converges quickly, should be used with time-varying channels.Another particularity of the BSS algorithms is their ability to restore, in the output, the sources with a certain permutation and a certain rotation or phase. Moreover, when a frequency offset is present, the constellations in output of adaptive-ACMA and ACMA spin. In order to avoid that and track the frequency offset jointly to the blind source separation, an MM or SCM cost functions can be used. To our knowledge, the MM and the SCM cost functions associated with analytical methods have never been tried so far (we have presented these works partly in conference [11]). So, in this paper an adaptive Analytical Multimodulus Algorithm and its simplified version named the adaptive Analytical Simplified Constant Modulus Algorithm are proposed.The performances of the adaptive-AMMA and adaptive-ASCMA are compared with those of the adaptive-ACMA and MMA. With time-selective channels, the adaptive-AMMA and the adaptive-ASCMA allow to obtain a higher signal-to-interference noise ratio (SINR) and a lower bit error rate (BER), than those obtained with the adaptive-ACMA and the MMA. Moreover, both proposed methods perform a source separation and a carrier phase recovery, unlike the adaptive-ACMA.First of all, in this paper, the system model used here is described. Then a batch analysis of the blind source separation algorithm AMMA is performed, and an adaptive form of this algorithm and of its simplified version is presented. Finally, simulation results of the proposed algorithms are compared to those obtained with the adaptive-ACMA and MMA in time-varying envrionments.The Kronecker and Khatri-Rao products have the following properties:vec(abH)=b*⊗a,(A⊗B)(C⊗D)=AC⊗BD,vec(ABC)=(CT⊗A)vec(B),vec(Adiag(b)C)=(CT∘A)b.


## 2. System Model

A MIMO system withNt transmitters and Nr receivers is considered. In this system, the information symbols x(k) are transmitted through the time-varying MIMO Nr×Nt channel H(k) at the time k. Taking into consideration a carrier frequency offset δf, the Nr×1 vector of the received signals y(k) can be expressed as(1)y(k)=H(k)x(k)e-j2πkδfTs+b(k),
where(2)y(k)=(y1(k)⋮yNr(k)),b(k)=(b1(k)⋮bNr(k)),H(k)=(h11(k)h1Nt(k)⋮⋮hNr1(k)hNrNt(k)),x(k)=(x1(k)⋮xNt(k)),
and(i)
y is an (Nr×1) representing the baseband received symbols without any time oversampling,(ii)
x is an (Nt×1) vector representing the transmitted signals (sources). The sources are independent, identically distributed (i.i.d.), and mutually independent with zero mean and unit variance,(iii)
b is an (Nr×1) vector representing the additive Gaussian noise,(iv)
H is an (Nr×Nt) unknown, full column rank, instantaneous linear and time-selective MIMO channel,(v)
1/Ts is the baud rate of the transmitter.


## 3.   The AMMA: A BSS Analytical Method so as to Minimize the MM Cost Function

### 3.1. Hypotheses

The BSS can be performed under the following hypotheses.(i)
Nr≥Nt.(ii)
H has i.i.d. complex components, is unitary unlike the received signals y(k), and must be prewhitened before applying the BSS.(iii)
The noise is additive, white, and Gaussian with a zero mean and a covariance matrixCb=E[bbH]=σb2INr.(iv)
The sources are zero mean discrete-time sequences, with a covariance matrixCx=E[xxH]=INt, and they must be mutually independent at a given time and identically distributed.


### 3.2. BSS's Principle

The BSS objective is to find a matrix(Nr×Nt)W(k) such as(3)z(k)=WH(k)y(k),
where the output vector sequence z(k) contains accurate estimates of the Nt source signals x(k). Each column of W(k) allows to find one source. So, all the columns of W(k) must be orthogonal with one another, in order to obtain independent sources. The matrix channel H must be unitary, otherwise the received signals will have to be prewhitened. Generally H is not unitary, so we will consider a prewithening operation which also reduces the dimension of x(k) from Nr to Nt. An underscore is used here to denote pre-whitened variables. Thus, let y̲(k)=FHy(k) be the prewhitened received signals, where F is an (Nr×Nt) prefilter. Like in [8], the prewhitening algorithm in [12] is used. The initial problem to find W(k) is then replaced by finding a (Nt×Nt) matrix T(k):(4)z(k)=TH(k)y̲(k)=TH(k)FH(k)y(k)=WH(k)y(k)=WH(k)(H(k)x(k)+b(k))=ΞH(k)x(k)+b′(k),
where W=FT and Ξ=HH(k)FT are a global separation matrix. Unfortunately, the calculated matrix T(k) separates the sources, except for a possible permutation of z(k) and an arbitrary phase on each source signal. The different BSS cost functions, described in Section 3.3, are indeed minimized with any permutation and any rotation on sources. Thus after separation, separator outputs present the following form:(5)z(k)=ΦΓx(k),
where Φ is a diagonal matrix and Γ is a permutation matrix which represent the arbitrary phase and the permutation, respectively.In order to findT(k), several cost functions can be used. The next section recalls the well-known (Constant Modulus) CM cost function and (Multimodulus) MM cost function.


### 3.3. Cost Function

#### 3.3.1. CM Cost Function

Many telecommunication signals have a constant modulus (PSK, 4-QAM), normalized to one generally. If the samples of one source are shown on a complex plan, every samples are on an unitary circle, but it is not true if we show a linear combination of several sources. It is this property that cost function CM exploits. Thus, the algorithms using the CM cost function make use of the sources' constant-modulus property:(6)J(T(k))=∑n=1NtE[(|zn(k)|2-R)2]=∑n=1NtE[(|tnH(k)y̲(k)|2-R)2],
where tn represents the nth row of T and R=Ri=Rr=E[|ℜ(x(k))|4]/E[|ℜ(x(k))|2]. This criterion forces BSS outputs to be on a circle with radius R.


#### 3.3.2. MM Cost Function

The cost function is composed of two cost functions; one cost function for real part(ℜ) and one for imaginary part (ℑ) of the equalizer output z(k):(7)J(T(k))=∑n=1Nt(E[(ℜ2(tnH(k)y̲(k))-Rr)2]+E[(ℑ2(tnH(k)y̲(k))-Ri)2]).
The error estimation for the real and imaginary parts is done jointly, and so the cost function uses implicitly the phase of the equalizer output. Thus, the carrier phase recovery is also accomplished with blind equalization, while the equalizer output signal constellation, obtained when using the CM cost function, suffers from an arbitrary phase rotation. So, the MM cost function is most interesting and used later in this paper.


#### 3.3.3. SCM Cost Function

The SCM cost function [5] consists in projection of the equalizer outputs on one dimension (either real or imaginary part). This cost function presents a lower computational complexity compared to the CM and MM cost functions:(8)J(T(k))=∑n=1NtE[(ℜ2(tnH(k)y̲(k))-Rr)2].


### 3.4. Implementation of Cost Function

Several methods exist to implement cost functions. Among them are the stochastic gradient descent and analytical methods.The algorithms using stochastic gradient are adaptive algorithms which are slow to converge. Furthermore, in varying environments, subsequent signals tend to be less accurate because of a maladjustment. These algorithms, using a stochastic gradient in order to minimize the cost functions CM and MM, are named, respectively, (Constant Modulus Algorithm) CMA [3, 4] and (MultiModulus Algorithm) MMA [1].Analytical methods exist in both the adaptive form and batch form. With the adaptive form, the algorithm converges quickly, so it is an interesting option in varying environments.Our choice is to use, for varying environments, an adaptive analytical method so as to minimize the MM cost function. So, a new algorithm named AMMA and which uses this association, is proposed, once a batch analysis has been performed.


### 3.5. The Analytical Multimodulus Algorithm

The following cost function should be minimized:(9)J(T̃)=∑l=1NtE[(ℜ2(zl(k))-Rr)2+(ℑ2(zl(k))-Ri)2].
In order to bring out the real and imaginary parts zl(k), this last term is written as(10)zl(k)=tlHy̲(k),l∈{1,…,Nt},=t̃lTỹ̲(k)+jt̃lTy̲̅(k),
where t̃l=(ℜ(tl)ℑ(tl)),ỹ̲(k)=(ℜ(y̲(k))ℑ(y̲(k))),y̲̅(k)=(ℑ(y̲(k))-ℜ(y̲(k))),y̲̅(k)=Gỹ̲(k) with G=(0NtINt-INt0Nt).The matrixG has the following properties:(11)-G=GT=G-1,G⊗G=(G⊗G)T=(G⊗G)-1.
Using Kronecker product properties and taking into account [9], ℜ2(zn(k)) and ℑ2(zn(k)) can be written as(12)ℜ2(zl(k))=t̃lT(ỹ̲(k)ỹ̲T(k))t̃l=(ỹ̲(k)⊗ỹ̲(k))T(t̃l⊗t̃l),ℑ2(zl(k))=t̃lT(y̲̅(k)y̲̅T(k))t̃l=(y̲̅(k)⊗y̲̅(k))T(t̃l⊗t̃l)=[(G⊗G)(ỹ̲(k)⊗ỹ̲(k))]T(t̃l⊗t̃l).


#### 3.5.1. The Batch Analysis

Let us denoteN a fixed number of received signals, the N rows of (ỹ̲(k)⊗ỹ̲(k))T and (y̲̅(k)⊗y̲̅(k))T are then, respectively, stacked in the N×(2Nt)2 matrix P̃ and P̅.Using the Khatri-Rao product (A∘B=(a1⊗b1a2⊗b2⋯)), P̃ and P̅ can be written as P̃=(Ỹ̲∘Ỹ̲)T and P̅=(Y̲̅∘Y̲̅)T. Let us define the (2Nt)2×1 vectors dl=t̃l⊗t̃l and 1=(1,…,1)T, the cost function is then given by(13)J(T̃)=1N∑l=1Nt[∥P̃dl-R·1∥2+∥P̅dl-R·1∥2]‍
with R=Ri=Rr.The linear system can be rewritten by using a QR decomposition [7, 8].We haveQ, an Ns×Ns unitary matrix, and such as Q1=Ns(10).After doing the QR factorization of(1P̃) and (1P̅) and applying Q to (1P̃) and (1P̅), we obtain(14)Q(1P̃)=Ns(1γ̃T0Õ),Q(1P̅)=Ns(1γ̅T0O̅)
with 1 the column vector composed of Ns ones, 0 a column vector composed of Ns-1 zeros, γ̃T and γ̅T two 1×(2Nt)2 rows vectors and Õ and O̅ two Ns-1×(2Nt)2 matrices, then(15)P̃dl=R1⇔Q(1P̃)(-Rdl)=0⇔{γ̃Tdl=R,Õdl=0,P̅dl=R1⇔Q(1P̅)(-Rdl)=0⇔{γ̅Tdl=R,O̅dl=0.
The cost function becomes(16)JAMMA(dl)=1Ns∑l=1Nt[|γ̃Tdl-R|2+|γ̅Tdl-R|2+∥Õdl∥2+∥O̅dl∥2],
where |γ̃Tdl-R|2 and |γ̅Tdl-R|2 allow to avoid the trivial solution dl=0. By squaring (14), the expressions of vectors γ̃ and γ̅, and of matrices C̃=ÕTÕ and C̅=O̅TO̅ are(17)γ̃=1NsP̃T1=1Ns∑k=1Nsỹ̲(k)⊗ỹ̲(k),C̃=ÕTÕ=1NsP̃TP̃-1NsP̃T1·1Ns1TP̃=1Ns∑k=1Ns(ỹ̲(k)⊗ỹ̲(k))(ỹ̲(k)⊗ỹ̲(k))T-1Ns2(∑k=1Nsỹ̲(k)⊗ỹ̲(k))(∑k=1Nsỹ̲(k)⊗ỹ̲(k))T,γ̅=1NsP̅T1=1Ns∑k=1Nsy̲̅(k)⊗y̲̅(k),C̅=O̅TO̅=1NsP̅TP̅-1NsP̅T1·1Ns1TP̅=1Ns∑k=1Ns(y̲̅(k)⊗y̲̅(k))(y̲̅(k)⊗y̲̅(k))T-1Ns2(∑k=1Nsy̲̅(k)⊗y̲̅(k))(∑ky̲̅(k)⊗y̲̅(k))T.
By using the definition of the estimated covariance matrices of ỹ̲(k) and y̲̅(k), named, respectively, R̂ỹ̲ and R̂y̲̅, and by using the property of the Kronecker product, we obtain the following equalities:(18)γ̃=vec(1Ns∑k=1Nsỹ̲(k)ỹ̲T(k))=vec(R̂ỹ̲),γ̅=vec(1Ns∑k=1Nsy̲̅(k)y̲̅T(k))=vec(R̂y̲̅).
Thus(19)γ̃Tdl=vec(Rỹ̲)T(tl⊗tl)=tlTRỹ̲tl,γ̅Tdl=vec(Ry̲̅)T(tl⊗tl)=tlTRy̲̅tl.
The criterion JAMMA is minimized when γ̃Tdl=R,γ̅Tdl=R,∥Õdl∥2=0, and ∥O̅dl∥2=0. Afterwards, we consider γ̃Tdl=R, γ̅Tdl=R as constraints. Since received signals are pre-whitened and Rỹ̲=Ry̲̅=I thus(20)γ̃Tdl=γ̅Tdl=tlTRỹ̲tl=tlTRy̲̅tl=tlTtl=∥tl∥2=R.
Moreover, ∥dl∥2=(tl⊗tl)T(t̃l⊗tl)=∥tl∥4=R2. Thus, γ̃Tdl=γ̅Tdl=R⇔∥dl∥=R, and to a scaling which is not important the optimization problem becomes(21)tl=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1Nt(∥Õdl∥2+∥O̅dl∥2)=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1Nt(dlTÕTÕdl+dlTO̅TO̅dl)=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1NtdlT(k)(C̃(k)+C̅(k))dl(k).In the appendix, it is shown that vectorsdl,for alll∈{1⋯Nt}, which minimize the function ∑l=1NtdlT(C̃+C̅)dl under the hypotheses dl=t̃l⊗t̃l and ∥d̃l∥=R, are in the image of C̃+C̅ (c.f. Theorem 2 and Lemma 2). Moreover, the minimization of ∑l=1NtdlT(C̃+C̅)dl is equivalent to minimize dlT(C̃+C̅)dl,for alll∈{1⋯2} (c.f. Theorem 1 and Lemma 1).Now, we will see how to create vectorsdl that present a Kronecker structure, in the image of C̃+C̅.


#### 3.5.2. The Batch Method

The vectorsdl which minimize dlT(C̃+C̅)dl are eigenvectors associated with the smallest eigenvalues. However, these vectors have not necessarily a Kronecker structure. Thus, in order to obtain the vectors t̃l in the image and to have a Kronecker struture, we search a Nt×Nt full rank matrix Λ that relates the two bases such so that(22)T̃=argmin∥D-(T̃∘T̃)Λ∥F2
with T̃=(t̃1,…,t̃Nt) and D=(d1,…,dNt) the 2Nt×Nt and (2Nt)2×Nt matrix, respectively. Using the property of the Kronecker product, this minimization problem can be rewritten as(23)∥D-(T̃∘T̃)Λ∥F2=∑n=1Nt∥dn-(T̃∘T̃)λn∥F2=∑n=1Nt∥Dn-T̃ΛnT̃∥F2
with(i)
λn the nth vector of the Λ matrix,(ii)
Λn=diag(λn) is a diagonal matrix,(iii)
Dn=vec(dn).Please note that if the matrixT̃ was square, the minimization of (23) would be equivalent to a joint diagonalization. In order to obtain the T̃ matrix, let us suppose in first phase that it is indeed square. Thus, with the aid of joint diagonalization of Dn matrices, a 2Nt×2NtT̃′ matrix is obtained. This diagonalization can be done by using the algorithm described in [13]. The columns of the matrix T̃′ are dependent; the 2Nt×Nt matrix T̃ is extracted by keeping the columns of T̃′ linearly independent.Now, this algorithm is implemented thanks to an adaptive method. So, we must estimate in an adaptive manner,C̃, C̅, and the eigenvectors of C̃+C̅ associated with the smallest eigenvalues. The methods used are described in the next section. Then the adaptive tracking of the minor subspace of C̃+C̅ is performed thanks to the use of the NOOJA algorithm. Finally the adaptive update of the joint diagonalization is realized.


#### 3.5.3. The Adaptive Method

The functiondlT(C̃+C̅)dl will be minimized, at each time k, but iteratively(24)t̃l(k)=mindl=t̃l⊗t̃l∥dl∥=RdlT(k)(C̃(k)+C̅(k))dl(k),∀l∈{1,…,Nt}.
The matrices C̃ and C̅ were defined by a batch of N samples. In this adaptive method they are converted into an exponential window (λ-scaling). The details are shown in [8].The following method used to estimate iteratively the matrixC̃ and C̅ is a generalization of the Van der Veen method [9](25)C̃(k)=λC̃(k-1)+β(k)c̃(k)c̃T(k)=f(C̃(k-1),c̃(k)),(26)C̅(k)=λC̅(k-1)+β(k)c̅(k)c̅T(k)=f(C̅(k-1),c̅(k)).
Then, C̃(k) and C̅(k) are, respectively, the estimations of the autocorrelation matrices for the signals c̃(k) and c̅(k), with 0<λ<1,β(k)=(α(k-1)/α(k))λ(1-λ),α(k)=λα(k-1)+1-λ and(27)c̃(k)=ỹ̲(k)⊗ỹ̲(k)-p̃(k-1)α(k-1),c̅(k)=y̲̅(k)⊗y̲̅(k)-p̅(k-1)α(k-1),
where(28)p̃(k)=λp̃(k-1)+(1-λ)(ỹ̲(k)⊗ỹ̲(k)),p̅(k)=λp̅(k-1)+(1-λ)(y̲̅(k)⊗y̲̅(k)).
We saw that the solution of this minimization is in the image of C̃+C̅. In order to obtain these vectors in an adaptive manner, a subspace tracking algorithm is used, associated with an adaptive joint diagonalization.Subspace Tracking
Used alone, the (Normalized Orthogonal Oja) NOOJA algorithm of [14] is utilized to track minor subspace spanned by the eigenvectors relating to the smallest eigenvalues. But, when associated with an adaptive joint diagonalization, as described in the next section, the NOOJA algorithm tracks the minor subspace spanned by the eigenvectors in the image of C̃+C̅ relating to the smallest eigenvalues which are different from 0. This algorithm extracts adaptively the minor subspace V relating to the input signal r(k)'s autocorrelation matrix R, by maximization the following cost function:
(29)J(V)=E∥r(k)-VVTr(k)∥2=tr(R)-2tr(VTRV)+tr(VTRVVTV).
In our problem, the minor subspace of C̃+C̅ is searched. So, in order to use the NOOJA algorithm, C̃+C̅ must be expressed as C̃+C̅=E[ccT]. But since the function f (cf. (25) and (26)) is not bilinear, the vector c is unknown. Instead of searching the minor subspace of C̃+C̅, the minor subspace of both C̃ and C̅ will be searched. Since C̅ can not be expressed as
(30)C̅=(G⊗G)C̃(G⊗G)-1.
The matricesC̃ and C̅ are similar. They have therefore the same image. In order to find the minor subspace D=(d1,…,dNt) spanned by the eigenvectors and relating to the smallest eigenvalues different from 0 for both C̃ and C̅, we propose to maximize the following cost function:
(31)J(D)=E∥c̃(k)-DDTc̃(k)∥2+E∥c̅(k)-DDTC̅(k)∥2=tr(C̃+C̅)-2tr(DT(C̃+C̅)D)+tr(DT(C̃+C̅)DDTD).
Since C̃ and C̅ have the same image, the first term of the cost function is sufficient to obtain convergence with a constant channel. Using only the matrix C̃ is equivalent to minimize the cost function Simplified CM [5] in analytical way. This simplification reduces the complexity of the adaptive AMMA and we name it the adaptive Analytical Simplified CMA (ASCMA). Please note that the ASCMA batch analysis is the same with the AMMA batch analysis performed in Section 3.5.1. However, the study is solely about C̃ rather than C̃+C̅.
The maximization ofJ(D) may be achieved by using the gradient-descent technique:
(32)D(k+1)=D(k)-β∇J(k),
where
(33)∇J=(-2(C̃+C̅)+(C̃+C̅(k))DDT+DDT(C̃+C̅))D.
The time instantk has been removed to simplify the equations.
By usingDT(k)D(k)=I, the cost function can be simplified as
(34)J(D)=tr(C̃+C̅)-tr(DT(C̃+C̅)D).
By replacing (32) in (34) and by doing the gradient of J(D(k+1)) with respect to β(k), the optimal stepsize βopt(k) is obtained
(35)βopt(k)=tr(DT(k)(C̃(k)+C̅(k))∇J(k))tr(∇JT(k)(C̃(k)+C̅(k))∇J(k)).
Keeping in mind thatD(k) is orthogonal at each iteration, that is, DT(k)D(k)=I, and replacing C̃(k) and C̅(k) by the estimate c̃(k)c̃T(k) and c̅(k)c̅T(k), respectively, we can write
(36)∇J(k)=p̃(k)z̃T(k)+p̅(k)z̅T(k),
where z̃(k)=DT(k)c̃(k),z̅(k)=DT(k)c̅(k),p̃(k)=z̃(k)-D(k)z̃(k), and p̅(k)=z̅(k)-D(k)z̅(k). The optimal step size becomes
(37)βopt=(∥c̃∥2-∥z̃∥2)∥z̃∥2+(∥c̅∥2-∥z̅∥2)∥z̅∥2(∥c̃∥2-∥z̃∥2)2∥z̃∥2+(∥c̅∥2-∥z̅∥2)2∥z̅∥2.
The instant k has been removed in the βopt(k)'s expression in order to simplify the equation.
The updating ofD(k+1) is performed as
(38)D′(k+1)=D(k)-βopt(k)(p̃(k)z̃T(k)+p̅(k)z̅T(k)),
where the matrix D(k+1) is obtained thanks to a Gram-Schmidt orthogonalization of D′(k+1). The algorithm obtained is described in Algorithm 1.
Once the eigenvectors relating to the smallest eigenvalues of bothC̃(k) and C̅(k) are obtained, the constraint dl=t̃l⊗t̃l must be satisfied. In other words, dl must have a Kronecker structure and must be in the image of C̃(k)+C̅(k).Algorithm 1:Adaptive nullspace of C̃and C̅tracking by using NOOJA.
Initialize  D(0)for  k=1,2,…z̃=DT(k)c̃,z̅=DT(k)c̅,p̃=z̃-D(k)z̃,p̅=z̅-V(k)z̅computeβopt by using(37)D′(k+1)=D(k)-βopt(p̃z̃T+p̅z̅T)d1(k+1)=d1'(k+1)/∥d1'(k+1)∥for  q=2 to Nts=[0⋯0]Tfor  l=1 to q-1s=s+(dlT(k+1)dq'(k+1))dl(k+1)end fordq(k+1)=(dq'(k+1)-s)/∥(dq'(k+1)-s)∥end forend forAdaptive Update of the Joint Diagonalization
The method used to satisfy the constraintdl=t̃l⊗t̃l is done as in [9]. T̃∘T̃ is computed and regarded as the current estimate of the subspace basis (D=T̃∘T̃). Using this basis, the subspace update is performed by using the “Normalized Orthogonal OJA” giving D. The last step of this algorithm is the mapping of the columns dl of D to a Kronecker-product:
(39)dl=t̃l⊗t̃l=vec(t̃lt̃lT),vec-1(dl)=Dl=t̃lt̃lT.
Then a power iteration [15] is applied. This iteration can be written as
(40)t̃l(k+1)=Dlt̃l(k).
At the next iteration, the matrix D in the image, which has a Kronecker structure, is obtained with D=T̃∘T̃. Then, D is transmitted at the NOOJA algorithm, allowing to vectors dl obtained in output NOOJA to be in the nullspace and to have a Kronecker struture.Estimated Symbols
In order to obtain the transmitted symbols, the complex matrixT is rebuilt from T̃(41)tlm=t̃lm+jt̃l(m+Nt),l,m∈1,…,Nt.
In order to prevent the extraction of the same source many times, T is reconditioned after each update. The technique used here is the computation of a singular value decomposition of T: T=∑σjujvjH. Then the singular values of T, which are smaller than 0.5, are replaced by 1
(42)recond(T)=∑σj′ujvjH.
Then, the matrixT is applied on the pre-whitened received signals:
(43)ŝ(k)=TH(k)y̲(k).
The suggested adaptive-AMMA algorithm is described in  Algorithm 2.Algorithm 2:Adaptive Analytical-MMA.
Initialize  F,T̃,p̅=0,p̃=0,α=0for  k=1,2,…Update  F, the prewhitened filter,by usingy(k) and [12] we obtain y̲(k)=FHy(k)Update real vectorsỹ(k)=(ℜ(y(k))ℑ(y(k))),y̅(k)=(-ℑ(y(k))ℜ(y(k))),w̃l=(ℜ(w)ℑ(w))Update vector  c̃ and c̅c̃=ỹ̲(k)⊗ỹ̲(k)-p̃/α,c̅=y̲̅(k)⊗y̲̅(k)-p̅/αp̃=λp̃+(1-λ)(ỹ̲(k)⊗ỹ̲(k)),p̅=λp̅+(1-λ)(y̲̅(k)⊗y̲̅(k))α=λα+1-λD=T̃∘T̃RegardD as a basis of the nullspace of C̃ and C̅,and update it usingc̃ and c̅ (Algorithm 1)for  l=1,…,NtDl=vec-1(dl),t̃l=Dlt̃l,t̃l=t̃l/∥t̃l∥Obtain complex  Nt×Nt  matrix  T(k)for  m=1,…,2Nttlm=t̃lm+jt̃l(m+Nt)end forend forT=recond(T),T̃=(ℜ(T)ℑ(T))TEstimate the sourcesŝ(k)=TH(k)y̲(k)end for


## 3.1. Hypotheses

The BSS can be performed under the following hypotheses.(i)
Nr≥Nt.(ii)
H has i.i.d. complex components, is unitary unlike the received signals y(k), and must be prewhitened before applying the BSS.(iii)
The noise is additive, white, and Gaussian with a zero mean and a covariance matrixCb=E[bbH]=σb2INr.(iv)
The sources are zero mean discrete-time sequences, with a covariance matrixCx=E[xxH]=INt, and they must be mutually independent at a given time and identically distributed.


## 3.2. BSS's Principle

The BSS objective is to find a matrix(Nr×Nt)W(k) such as(3)z(k)=WH(k)y(k),
where the output vector sequence z(k) contains accurate estimates of the Nt source signals x(k). Each column of W(k) allows to find one source. So, all the columns of W(k) must be orthogonal with one another, in order to obtain independent sources. The matrix channel H must be unitary, otherwise the received signals will have to be prewhitened. Generally H is not unitary, so we will consider a prewithening operation which also reduces the dimension of x(k) from Nr to Nt. An underscore is used here to denote pre-whitened variables. Thus, let y̲(k)=FHy(k) be the prewhitened received signals, where F is an (Nr×Nt) prefilter. Like in [8], the prewhitening algorithm in [12] is used. The initial problem to find W(k) is then replaced by finding a (Nt×Nt) matrix T(k):(4)z(k)=TH(k)y̲(k)=TH(k)FH(k)y(k)=WH(k)y(k)=WH(k)(H(k)x(k)+b(k))=ΞH(k)x(k)+b′(k),
where W=FT and Ξ=HH(k)FT are a global separation matrix. Unfortunately, the calculated matrix T(k) separates the sources, except for a possible permutation of z(k) and an arbitrary phase on each source signal. The different BSS cost functions, described in Section 3.3, are indeed minimized with any permutation and any rotation on sources. Thus after separation, separator outputs present the following form:(5)z(k)=ΦΓx(k),
where Φ is a diagonal matrix and Γ is a permutation matrix which represent the arbitrary phase and the permutation, respectively.In order to findT(k), several cost functions can be used. The next section recalls the well-known (Constant Modulus) CM cost function and (Multimodulus) MM cost function.


## 3.3. Cost Function

### 3.3.1. CM Cost Function

Many telecommunication signals have a constant modulus (PSK, 4-QAM), normalized to one generally. If the samples of one source are shown on a complex plan, every samples are on an unitary circle, but it is not true if we show a linear combination of several sources. It is this property that cost function CM exploits. Thus, the algorithms using the CM cost function make use of the sources' constant-modulus property:(6)J(T(k))=∑n=1NtE[(|zn(k)|2-R)2]=∑n=1NtE[(|tnH(k)y̲(k)|2-R)2],
where tn represents the nth row of T and R=Ri=Rr=E[|ℜ(x(k))|4]/E[|ℜ(x(k))|2]. This criterion forces BSS outputs to be on a circle with radius R.


### 3.3.2. MM Cost Function

The cost function is composed of two cost functions; one cost function for real part(ℜ) and one for imaginary part (ℑ) of the equalizer output z(k):(7)J(T(k))=∑n=1Nt(E[(ℜ2(tnH(k)y̲(k))-Rr)2]+E[(ℑ2(tnH(k)y̲(k))-Ri)2]).
The error estimation for the real and imaginary parts is done jointly, and so the cost function uses implicitly the phase of the equalizer output. Thus, the carrier phase recovery is also accomplished with blind equalization, while the equalizer output signal constellation, obtained when using the CM cost function, suffers from an arbitrary phase rotation. So, the MM cost function is most interesting and used later in this paper.


### 3.3.3. SCM Cost Function

The SCM cost function [5] consists in projection of the equalizer outputs on one dimension (either real or imaginary part). This cost function presents a lower computational complexity compared to the CM and MM cost functions:(8)J(T(k))=∑n=1NtE[(ℜ2(tnH(k)y̲(k))-Rr)2].


## 3.3.1. CM Cost Function

Many telecommunication signals have a constant modulus (PSK, 4-QAM), normalized to one generally. If the samples of one source are shown on a complex plan, every samples are on an unitary circle, but it is not true if we show a linear combination of several sources. It is this property that cost function CM exploits. Thus, the algorithms using the CM cost function make use of the sources' constant-modulus property:(6)J(T(k))=∑n=1NtE[(|zn(k)|2-R)2]=∑n=1NtE[(|tnH(k)y̲(k)|2-R)2],
where tn represents the nth row of T and R=Ri=Rr=E[|ℜ(x(k))|4]/E[|ℜ(x(k))|2]. This criterion forces BSS outputs to be on a circle with radius R.


## 3.3.2. MM Cost Function

The cost function is composed of two cost functions; one cost function for real part(ℜ) and one for imaginary part (ℑ) of the equalizer output z(k):(7)J(T(k))=∑n=1Nt(E[(ℜ2(tnH(k)y̲(k))-Rr)2]+E[(ℑ2(tnH(k)y̲(k))-Ri)2]).
The error estimation for the real and imaginary parts is done jointly, and so the cost function uses implicitly the phase of the equalizer output. Thus, the carrier phase recovery is also accomplished with blind equalization, while the equalizer output signal constellation, obtained when using the CM cost function, suffers from an arbitrary phase rotation. So, the MM cost function is most interesting and used later in this paper.


## 3.3.3. SCM Cost Function

The SCM cost function [5] consists in projection of the equalizer outputs on one dimension (either real or imaginary part). This cost function presents a lower computational complexity compared to the CM and MM cost functions:(8)J(T(k))=∑n=1NtE[(ℜ2(tnH(k)y̲(k))-Rr)2].


## 3.4. Implementation of Cost Function

Several methods exist to implement cost functions. Among them are the stochastic gradient descent and analytical methods.The algorithms using stochastic gradient are adaptive algorithms which are slow to converge. Furthermore, in varying environments, subsequent signals tend to be less accurate because of a maladjustment. These algorithms, using a stochastic gradient in order to minimize the cost functions CM and MM, are named, respectively, (Constant Modulus Algorithm) CMA [3, 4] and (MultiModulus Algorithm) MMA [1].Analytical methods exist in both the adaptive form and batch form. With the adaptive form, the algorithm converges quickly, so it is an interesting option in varying environments.Our choice is to use, for varying environments, an adaptive analytical method so as to minimize the MM cost function. So, a new algorithm named AMMA and which uses this association, is proposed, once a batch analysis has been performed.


## 3.5. The Analytical Multimodulus Algorithm

The following cost function should be minimized:(9)J(T̃)=∑l=1NtE[(ℜ2(zl(k))-Rr)2+(ℑ2(zl(k))-Ri)2].
In order to bring out the real and imaginary parts zl(k), this last term is written as(10)zl(k)=tlHy̲(k),l∈{1,…,Nt},=t̃lTỹ̲(k)+jt̃lTy̲̅(k),
where t̃l=(ℜ(tl)ℑ(tl)),ỹ̲(k)=(ℜ(y̲(k))ℑ(y̲(k))),y̲̅(k)=(ℑ(y̲(k))-ℜ(y̲(k))),y̲̅(k)=Gỹ̲(k) with G=(0NtINt-INt0Nt).The matrixG has the following properties:(11)-G=GT=G-1,G⊗G=(G⊗G)T=(G⊗G)-1.
Using Kronecker product properties and taking into account [9], ℜ2(zn(k)) and ℑ2(zn(k)) can be written as(12)ℜ2(zl(k))=t̃lT(ỹ̲(k)ỹ̲T(k))t̃l=(ỹ̲(k)⊗ỹ̲(k))T(t̃l⊗t̃l),ℑ2(zl(k))=t̃lT(y̲̅(k)y̲̅T(k))t̃l=(y̲̅(k)⊗y̲̅(k))T(t̃l⊗t̃l)=[(G⊗G)(ỹ̲(k)⊗ỹ̲(k))]T(t̃l⊗t̃l).


### 3.5.1. The Batch Analysis

Let us denoteN a fixed number of received signals, the N rows of (ỹ̲(k)⊗ỹ̲(k))T and (y̲̅(k)⊗y̲̅(k))T are then, respectively, stacked in the N×(2Nt)2 matrix P̃ and P̅.Using the Khatri-Rao product (A∘B=(a1⊗b1a2⊗b2⋯)), P̃ and P̅ can be written as P̃=(Ỹ̲∘Ỹ̲)T and P̅=(Y̲̅∘Y̲̅)T. Let us define the (2Nt)2×1 vectors dl=t̃l⊗t̃l and 1=(1,…,1)T, the cost function is then given by(13)J(T̃)=1N∑l=1Nt[∥P̃dl-R·1∥2+∥P̅dl-R·1∥2]‍
with R=Ri=Rr.The linear system can be rewritten by using a QR decomposition [7, 8].We haveQ, an Ns×Ns unitary matrix, and such as Q1=Ns(10).After doing the QR factorization of(1P̃) and (1P̅) and applying Q to (1P̃) and (1P̅), we obtain(14)Q(1P̃)=Ns(1γ̃T0Õ),Q(1P̅)=Ns(1γ̅T0O̅)
with 1 the column vector composed of Ns ones, 0 a column vector composed of Ns-1 zeros, γ̃T and γ̅T two 1×(2Nt)2 rows vectors and Õ and O̅ two Ns-1×(2Nt)2 matrices, then(15)P̃dl=R1⇔Q(1P̃)(-Rdl)=0⇔{γ̃Tdl=R,Õdl=0,P̅dl=R1⇔Q(1P̅)(-Rdl)=0⇔{γ̅Tdl=R,O̅dl=0.
The cost function becomes(16)JAMMA(dl)=1Ns∑l=1Nt[|γ̃Tdl-R|2+|γ̅Tdl-R|2+∥Õdl∥2+∥O̅dl∥2],
where |γ̃Tdl-R|2 and |γ̅Tdl-R|2 allow to avoid the trivial solution dl=0. By squaring (14), the expressions of vectors γ̃ and γ̅, and of matrices C̃=ÕTÕ and C̅=O̅TO̅ are(17)γ̃=1NsP̃T1=1Ns∑k=1Nsỹ̲(k)⊗ỹ̲(k),C̃=ÕTÕ=1NsP̃TP̃-1NsP̃T1·1Ns1TP̃=1Ns∑k=1Ns(ỹ̲(k)⊗ỹ̲(k))(ỹ̲(k)⊗ỹ̲(k))T-1Ns2(∑k=1Nsỹ̲(k)⊗ỹ̲(k))(∑k=1Nsỹ̲(k)⊗ỹ̲(k))T,γ̅=1NsP̅T1=1Ns∑k=1Nsy̲̅(k)⊗y̲̅(k),C̅=O̅TO̅=1NsP̅TP̅-1NsP̅T1·1Ns1TP̅=1Ns∑k=1Ns(y̲̅(k)⊗y̲̅(k))(y̲̅(k)⊗y̲̅(k))T-1Ns2(∑k=1Nsy̲̅(k)⊗y̲̅(k))(∑ky̲̅(k)⊗y̲̅(k))T.
By using the definition of the estimated covariance matrices of ỹ̲(k) and y̲̅(k), named, respectively, R̂ỹ̲ and R̂y̲̅, and by using the property of the Kronecker product, we obtain the following equalities:(18)γ̃=vec(1Ns∑k=1Nsỹ̲(k)ỹ̲T(k))=vec(R̂ỹ̲),γ̅=vec(1Ns∑k=1Nsy̲̅(k)y̲̅T(k))=vec(R̂y̲̅).
Thus(19)γ̃Tdl=vec(Rỹ̲)T(tl⊗tl)=tlTRỹ̲tl,γ̅Tdl=vec(Ry̲̅)T(tl⊗tl)=tlTRy̲̅tl.
The criterion JAMMA is minimized when γ̃Tdl=R,γ̅Tdl=R,∥Õdl∥2=0, and ∥O̅dl∥2=0. Afterwards, we consider γ̃Tdl=R, γ̅Tdl=R as constraints. Since received signals are pre-whitened and Rỹ̲=Ry̲̅=I thus(20)γ̃Tdl=γ̅Tdl=tlTRỹ̲tl=tlTRy̲̅tl=tlTtl=∥tl∥2=R.
Moreover, ∥dl∥2=(tl⊗tl)T(t̃l⊗tl)=∥tl∥4=R2. Thus, γ̃Tdl=γ̅Tdl=R⇔∥dl∥=R, and to a scaling which is not important the optimization problem becomes(21)tl=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1Nt(∥Õdl∥2+∥O̅dl∥2)=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1Nt(dlTÕTÕdl+dlTO̅TO̅dl)=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1NtdlT(k)(C̃(k)+C̅(k))dl(k).In the appendix, it is shown that vectorsdl,for alll∈{1⋯Nt}, which minimize the function ∑l=1NtdlT(C̃+C̅)dl under the hypotheses dl=t̃l⊗t̃l and ∥d̃l∥=R, are in the image of C̃+C̅ (c.f. Theorem 2 and Lemma 2). Moreover, the minimization of ∑l=1NtdlT(C̃+C̅)dl is equivalent to minimize dlT(C̃+C̅)dl,for alll∈{1⋯2} (c.f. Theorem 1 and Lemma 1).Now, we will see how to create vectorsdl that present a Kronecker structure, in the image of C̃+C̅.


### 3.5.2. The Batch Method

The vectorsdl which minimize dlT(C̃+C̅)dl are eigenvectors associated with the smallest eigenvalues. However, these vectors have not necessarily a Kronecker structure. Thus, in order to obtain the vectors t̃l in the image and to have a Kronecker struture, we search a Nt×Nt full rank matrix Λ that relates the two bases such so that(22)T̃=argmin∥D-(T̃∘T̃)Λ∥F2
with T̃=(t̃1,…,t̃Nt) and D=(d1,…,dNt) the 2Nt×Nt and (2Nt)2×Nt matrix, respectively. Using the property of the Kronecker product, this minimization problem can be rewritten as(23)∥D-(T̃∘T̃)Λ∥F2=∑n=1Nt∥dn-(T̃∘T̃)λn∥F2=∑n=1Nt∥Dn-T̃ΛnT̃∥F2
with(i)
λn the nth vector of the Λ matrix,(ii)
Λn=diag(λn) is a diagonal matrix,(iii)
Dn=vec(dn).Please note that if the matrixT̃ was square, the minimization of (23) would be equivalent to a joint diagonalization. In order to obtain the T̃ matrix, let us suppose in first phase that it is indeed square. Thus, with the aid of joint diagonalization of Dn matrices, a 2Nt×2NtT̃′ matrix is obtained. This diagonalization can be done by using the algorithm described in [13]. The columns of the matrix T̃′ are dependent; the 2Nt×Nt matrix T̃ is extracted by keeping the columns of T̃′ linearly independent.Now, this algorithm is implemented thanks to an adaptive method. So, we must estimate in an adaptive manner,C̃, C̅, and the eigenvectors of C̃+C̅ associated with the smallest eigenvalues. The methods used are described in the next section. Then the adaptive tracking of the minor subspace of C̃+C̅ is performed thanks to the use of the NOOJA algorithm. Finally the adaptive update of the joint diagonalization is realized.


### 3.5.3. The Adaptive Method

The functiondlT(C̃+C̅)dl will be minimized, at each time k, but iteratively(24)t̃l(k)=mindl=t̃l⊗t̃l∥dl∥=RdlT(k)(C̃(k)+C̅(k))dl(k),∀l∈{1,…,Nt}.
The matrices C̃ and C̅ were defined by a batch of N samples. In this adaptive method they are converted into an exponential window (λ-scaling). The details are shown in [8].The following method used to estimate iteratively the matrixC̃ and C̅ is a generalization of the Van der Veen method [9](25)C̃(k)=λC̃(k-1)+β(k)c̃(k)c̃T(k)=f(C̃(k-1),c̃(k)),(26)C̅(k)=λC̅(k-1)+β(k)c̅(k)c̅T(k)=f(C̅(k-1),c̅(k)).
Then, C̃(k) and C̅(k) are, respectively, the estimations of the autocorrelation matrices for the signals c̃(k) and c̅(k), with 0<λ<1,β(k)=(α(k-1)/α(k))λ(1-λ),α(k)=λα(k-1)+1-λ and(27)c̃(k)=ỹ̲(k)⊗ỹ̲(k)-p̃(k-1)α(k-1),c̅(k)=y̲̅(k)⊗y̲̅(k)-p̅(k-1)α(k-1),
where(28)p̃(k)=λp̃(k-1)+(1-λ)(ỹ̲(k)⊗ỹ̲(k)),p̅(k)=λp̅(k-1)+(1-λ)(y̲̅(k)⊗y̲̅(k)).
We saw that the solution of this minimization is in the image of C̃+C̅. In order to obtain these vectors in an adaptive manner, a subspace tracking algorithm is used, associated with an adaptive joint diagonalization.Subspace Tracking
Used alone, the (Normalized Orthogonal Oja) NOOJA algorithm of [14] is utilized to track minor subspace spanned by the eigenvectors relating to the smallest eigenvalues. But, when associated with an adaptive joint diagonalization, as described in the next section, the NOOJA algorithm tracks the minor subspace spanned by the eigenvectors in the image of C̃+C̅ relating to the smallest eigenvalues which are different from 0. This algorithm extracts adaptively the minor subspace V relating to the input signal r(k)'s autocorrelation matrix R, by maximization the following cost function:
(29)J(V)=E∥r(k)-VVTr(k)∥2=tr(R)-2tr(VTRV)+tr(VTRVVTV).
In our problem, the minor subspace of C̃+C̅ is searched. So, in order to use the NOOJA algorithm, C̃+C̅ must be expressed as C̃+C̅=E[ccT]. But since the function f (cf. (25) and (26)) is not bilinear, the vector c is unknown. Instead of searching the minor subspace of C̃+C̅, the minor subspace of both C̃ and C̅ will be searched. Since C̅ can not be expressed as
(30)C̅=(G⊗G)C̃(G⊗G)-1.
The matricesC̃ and C̅ are similar. They have therefore the same image. In order to find the minor subspace D=(d1,…,dNt) spanned by the eigenvectors and relating to the smallest eigenvalues different from 0 for both C̃ and C̅, we propose to maximize the following cost function:
(31)J(D)=E∥c̃(k)-DDTc̃(k)∥2+E∥c̅(k)-DDTC̅(k)∥2=tr(C̃+C̅)-2tr(DT(C̃+C̅)D)+tr(DT(C̃+C̅)DDTD).
Since C̃ and C̅ have the same image, the first term of the cost function is sufficient to obtain convergence with a constant channel. Using only the matrix C̃ is equivalent to minimize the cost function Simplified CM [5] in analytical way. This simplification reduces the complexity of the adaptive AMMA and we name it the adaptive Analytical Simplified CMA (ASCMA). Please note that the ASCMA batch analysis is the same with the AMMA batch analysis performed in Section 3.5.1. However, the study is solely about C̃ rather than C̃+C̅.
The maximization ofJ(D) may be achieved by using the gradient-descent technique:
(32)D(k+1)=D(k)-β∇J(k),
where
(33)∇J=(-2(C̃+C̅)+(C̃+C̅(k))DDT+DDT(C̃+C̅))D.
The time instantk has been removed to simplify the equations.
By usingDT(k)D(k)=I, the cost function can be simplified as
(34)J(D)=tr(C̃+C̅)-tr(DT(C̃+C̅)D).
By replacing (32) in (34) and by doing the gradient of J(D(k+1)) with respect to β(k), the optimal stepsize βopt(k) is obtained
(35)βopt(k)=tr(DT(k)(C̃(k)+C̅(k))∇J(k))tr(∇JT(k)(C̃(k)+C̅(k))∇J(k)).
Keeping in mind thatD(k) is orthogonal at each iteration, that is, DT(k)D(k)=I, and replacing C̃(k) and C̅(k) by the estimate c̃(k)c̃T(k) and c̅(k)c̅T(k), respectively, we can write
(36)∇J(k)=p̃(k)z̃T(k)+p̅(k)z̅T(k),
where z̃(k)=DT(k)c̃(k),z̅(k)=DT(k)c̅(k),p̃(k)=z̃(k)-D(k)z̃(k), and p̅(k)=z̅(k)-D(k)z̅(k). The optimal step size becomes
(37)βopt=(∥c̃∥2-∥z̃∥2)∥z̃∥2+(∥c̅∥2-∥z̅∥2)∥z̅∥2(∥c̃∥2-∥z̃∥2)2∥z̃∥2+(∥c̅∥2-∥z̅∥2)2∥z̅∥2.
The instant k has been removed in the βopt(k)'s expression in order to simplify the equation.
The updating ofD(k+1) is performed as
(38)D′(k+1)=D(k)-βopt(k)(p̃(k)z̃T(k)+p̅(k)z̅T(k)),
where the matrix D(k+1) is obtained thanks to a Gram-Schmidt orthogonalization of D′(k+1). The algorithm obtained is described in Algorithm 1.
Once the eigenvectors relating to the smallest eigenvalues of bothC̃(k) and C̅(k) are obtained, the constraint dl=t̃l⊗t̃l must be satisfied. In other words, dl must have a Kronecker structure and must be in the image of C̃(k)+C̅(k).Algorithm 1:Adaptive nullspace of C̃and C̅tracking by using NOOJA.
Initialize  D(0)for  k=1,2,…z̃=DT(k)c̃,z̅=DT(k)c̅,p̃=z̃-D(k)z̃,p̅=z̅-V(k)z̅computeβopt by using(37)D′(k+1)=D(k)-βopt(p̃z̃T+p̅z̅T)d1(k+1)=d1'(k+1)/∥d1'(k+1)∥for  q=2 to Nts=[0⋯0]Tfor  l=1 to q-1s=s+(dlT(k+1)dq'(k+1))dl(k+1)end fordq(k+1)=(dq'(k+1)-s)/∥(dq'(k+1)-s)∥end forend forAdaptive Update of the Joint Diagonalization
The method used to satisfy the constraintdl=t̃l⊗t̃l is done as in [9]. T̃∘T̃ is computed and regarded as the current estimate of the subspace basis (D=T̃∘T̃). Using this basis, the subspace update is performed by using the “Normalized Orthogonal OJA” giving D. The last step of this algorithm is the mapping of the columns dl of D to a Kronecker-product:
(39)dl=t̃l⊗t̃l=vec(t̃lt̃lT),vec-1(dl)=Dl=t̃lt̃lT.
Then a power iteration [15] is applied. This iteration can be written as
(40)t̃l(k+1)=Dlt̃l(k).
At the next iteration, the matrix D in the image, which has a Kronecker structure, is obtained with D=T̃∘T̃. Then, D is transmitted at the NOOJA algorithm, allowing to vectors dl obtained in output NOOJA to be in the nullspace and to have a Kronecker struture.Estimated Symbols
In order to obtain the transmitted symbols, the complex matrixT is rebuilt from T̃(41)tlm=t̃lm+jt̃l(m+Nt),l,m∈1,…,Nt.
In order to prevent the extraction of the same source many times, T is reconditioned after each update. The technique used here is the computation of a singular value decomposition of T: T=∑σjujvjH. Then the singular values of T, which are smaller than 0.5, are replaced by 1
(42)recond(T)=∑σj′ujvjH.
Then, the matrixT is applied on the pre-whitened received signals:
(43)ŝ(k)=TH(k)y̲(k).
The suggested adaptive-AMMA algorithm is described in  Algorithm 2.Algorithm 2:Adaptive Analytical-MMA.
Initialize  F,T̃,p̅=0,p̃=0,α=0for  k=1,2,…Update  F, the prewhitened filter,by usingy(k) and [12] we obtain y̲(k)=FHy(k)Update real vectorsỹ(k)=(ℜ(y(k))ℑ(y(k))),y̅(k)=(-ℑ(y(k))ℜ(y(k))),w̃l=(ℜ(w)ℑ(w))Update vector  c̃ and c̅c̃=ỹ̲(k)⊗ỹ̲(k)-p̃/α,c̅=y̲̅(k)⊗y̲̅(k)-p̅/αp̃=λp̃+(1-λ)(ỹ̲(k)⊗ỹ̲(k)),p̅=λp̅+(1-λ)(y̲̅(k)⊗y̲̅(k))α=λα+1-λD=T̃∘T̃RegardD as a basis of the nullspace of C̃ and C̅,and update it usingc̃ and c̅ (Algorithm 1)for  l=1,…,NtDl=vec-1(dl),t̃l=Dlt̃l,t̃l=t̃l/∥t̃l∥Obtain complex  Nt×Nt  matrix  T(k)for  m=1,…,2Nttlm=t̃lm+jt̃l(m+Nt)end forend forT=recond(T),T̃=(ℜ(T)ℑ(T))TEstimate the sourcesŝ(k)=TH(k)y̲(k)end for


## 3.5.1. The Batch Analysis

Let us denoteN a fixed number of received signals, the N rows of (ỹ̲(k)⊗ỹ̲(k))T and (y̲̅(k)⊗y̲̅(k))T are then, respectively, stacked in the N×(2Nt)2 matrix P̃ and P̅.Using the Khatri-Rao product (A∘B=(a1⊗b1a2⊗b2⋯)), P̃ and P̅ can be written as P̃=(Ỹ̲∘Ỹ̲)T and P̅=(Y̲̅∘Y̲̅)T. Let us define the (2Nt)2×1 vectors dl=t̃l⊗t̃l and 1=(1,…,1)T, the cost function is then given by(13)J(T̃)=1N∑l=1Nt[∥P̃dl-R·1∥2+∥P̅dl-R·1∥2]‍
with R=Ri=Rr.The linear system can be rewritten by using a QR decomposition [7, 8].We haveQ, an Ns×Ns unitary matrix, and such as Q1=Ns(10).After doing the QR factorization of(1P̃) and (1P̅) and applying Q to (1P̃) and (1P̅), we obtain(14)Q(1P̃)=Ns(1γ̃T0Õ),Q(1P̅)=Ns(1γ̅T0O̅)
with 1 the column vector composed of Ns ones, 0 a column vector composed of Ns-1 zeros, γ̃T and γ̅T two 1×(2Nt)2 rows vectors and Õ and O̅ two Ns-1×(2Nt)2 matrices, then(15)P̃dl=R1⇔Q(1P̃)(-Rdl)=0⇔{γ̃Tdl=R,Õdl=0,P̅dl=R1⇔Q(1P̅)(-Rdl)=0⇔{γ̅Tdl=R,O̅dl=0.
The cost function becomes(16)JAMMA(dl)=1Ns∑l=1Nt[|γ̃Tdl-R|2+|γ̅Tdl-R|2+∥Õdl∥2+∥O̅dl∥2],
where |γ̃Tdl-R|2 and |γ̅Tdl-R|2 allow to avoid the trivial solution dl=0. By squaring (14), the expressions of vectors γ̃ and γ̅, and of matrices C̃=ÕTÕ and C̅=O̅TO̅ are(17)γ̃=1NsP̃T1=1Ns∑k=1Nsỹ̲(k)⊗ỹ̲(k),C̃=ÕTÕ=1NsP̃TP̃-1NsP̃T1·1Ns1TP̃=1Ns∑k=1Ns(ỹ̲(k)⊗ỹ̲(k))(ỹ̲(k)⊗ỹ̲(k))T-1Ns2(∑k=1Nsỹ̲(k)⊗ỹ̲(k))(∑k=1Nsỹ̲(k)⊗ỹ̲(k))T,γ̅=1NsP̅T1=1Ns∑k=1Nsy̲̅(k)⊗y̲̅(k),C̅=O̅TO̅=1NsP̅TP̅-1NsP̅T1·1Ns1TP̅=1Ns∑k=1Ns(y̲̅(k)⊗y̲̅(k))(y̲̅(k)⊗y̲̅(k))T-1Ns2(∑k=1Nsy̲̅(k)⊗y̲̅(k))(∑ky̲̅(k)⊗y̲̅(k))T.
By using the definition of the estimated covariance matrices of ỹ̲(k) and y̲̅(k), named, respectively, R̂ỹ̲ and R̂y̲̅, and by using the property of the Kronecker product, we obtain the following equalities:(18)γ̃=vec(1Ns∑k=1Nsỹ̲(k)ỹ̲T(k))=vec(R̂ỹ̲),γ̅=vec(1Ns∑k=1Nsy̲̅(k)y̲̅T(k))=vec(R̂y̲̅).
Thus(19)γ̃Tdl=vec(Rỹ̲)T(tl⊗tl)=tlTRỹ̲tl,γ̅Tdl=vec(Ry̲̅)T(tl⊗tl)=tlTRy̲̅tl.
The criterion JAMMA is minimized when γ̃Tdl=R,γ̅Tdl=R,∥Õdl∥2=0, and ∥O̅dl∥2=0. Afterwards, we consider γ̃Tdl=R, γ̅Tdl=R as constraints. Since received signals are pre-whitened and Rỹ̲=Ry̲̅=I thus(20)γ̃Tdl=γ̅Tdl=tlTRỹ̲tl=tlTRy̲̅tl=tlTtl=∥tl∥2=R.
Moreover, ∥dl∥2=(tl⊗tl)T(t̃l⊗tl)=∥tl∥4=R2. Thus, γ̃Tdl=γ̅Tdl=R⇔∥dl∥=R, and to a scaling which is not important the optimization problem becomes(21)tl=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1Nt(∥Õdl∥2+∥O̅dl∥2)=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1Nt(dlTÕTÕdl+dlTO̅TO̅dl)=argmindl(k)=tl(k)⊗tl(k)∥dl(k)∥=R∑l=1NtdlT(k)(C̃(k)+C̅(k))dl(k).In the appendix, it is shown that vectorsdl,for alll∈{1⋯Nt}, which minimize the function ∑l=1NtdlT(C̃+C̅)dl under the hypotheses dl=t̃l⊗t̃l and ∥d̃l∥=R, are in the image of C̃+C̅ (c.f. Theorem 2 and Lemma 2). Moreover, the minimization of ∑l=1NtdlT(C̃+C̅)dl is equivalent to minimize dlT(C̃+C̅)dl,for alll∈{1⋯2} (c.f. Theorem 1 and Lemma 1).Now, we will see how to create vectorsdl that present a Kronecker structure, in the image of C̃+C̅.


## 3.5.2. The Batch Method

The vectorsdl which minimize dlT(C̃+C̅)dl are eigenvectors associated with the smallest eigenvalues. However, these vectors have not necessarily a Kronecker structure. Thus, in order to obtain the vectors t̃l in the image and to have a Kronecker struture, we search a Nt×Nt full rank matrix Λ that relates the two bases such so that(22)T̃=argmin∥D-(T̃∘T̃)Λ∥F2
with T̃=(t̃1,…,t̃Nt) and D=(d1,…,dNt) the 2Nt×Nt and (2Nt)2×Nt matrix, respectively. Using the property of the Kronecker product, this minimization problem can be rewritten as(23)∥D-(T̃∘T̃)Λ∥F2=∑n=1Nt∥dn-(T̃∘T̃)λn∥F2=∑n=1Nt∥Dn-T̃ΛnT̃∥F2
with(i)
λn the nth vector of the Λ matrix,(ii)
Λn=diag(λn) is a diagonal matrix,(iii)
Dn=vec(dn).Please note that if the matrixT̃ was square, the minimization of (23) would be equivalent to a joint diagonalization. In order to obtain the T̃ matrix, let us suppose in first phase that it is indeed square. Thus, with the aid of joint diagonalization of Dn matrices, a 2Nt×2NtT̃′ matrix is obtained. This diagonalization can be done by using the algorithm described in [13]. The columns of the matrix T̃′ are dependent; the 2Nt×Nt matrix T̃ is extracted by keeping the columns of T̃′ linearly independent.Now, this algorithm is implemented thanks to an adaptive method. So, we must estimate in an adaptive manner,C̃, C̅, and the eigenvectors of C̃+C̅ associated with the smallest eigenvalues. The methods used are described in the next section. Then the adaptive tracking of the minor subspace of C̃+C̅ is performed thanks to the use of the NOOJA algorithm. Finally the adaptive update of the joint diagonalization is realized.


## 3.5.3. The Adaptive Method

The functiondlT(C̃+C̅)dl will be minimized, at each time k, but iteratively(24)t̃l(k)=mindl=t̃l⊗t̃l∥dl∥=RdlT(k)(C̃(k)+C̅(k))dl(k),∀l∈{1,…,Nt}.
The matrices C̃ and C̅ were defined by a batch of N samples. In this adaptive method they are converted into an exponential window (λ-scaling). The details are shown in [8].The following method used to estimate iteratively the matrixC̃ and C̅ is a generalization of the Van der Veen method [9](25)C̃(k)=λC̃(k-1)+β(k)c̃(k)c̃T(k)=f(C̃(k-1),c̃(k)),(26)C̅(k)=λC̅(k-1)+β(k)c̅(k)c̅T(k)=f(C̅(k-1),c̅(k)).
Then, C̃(k) and C̅(k) are, respectively, the estimations of the autocorrelation matrices for the signals c̃(k) and c̅(k), with 0<λ<1,β(k)=(α(k-1)/α(k))λ(1-λ),α(k)=λα(k-1)+1-λ and(27)c̃(k)=ỹ̲(k)⊗ỹ̲(k)-p̃(k-1)α(k-1),c̅(k)=y̲̅(k)⊗y̲̅(k)-p̅(k-1)α(k-1),
where(28)p̃(k)=λp̃(k-1)+(1-λ)(ỹ̲(k)⊗ỹ̲(k)),p̅(k)=λp̅(k-1)+(1-λ)(y̲̅(k)⊗y̲̅(k)).
We saw that the solution of this minimization is in the image of C̃+C̅. In order to obtain these vectors in an adaptive manner, a subspace tracking algorithm is used, associated with an adaptive joint diagonalization.Subspace Tracking
Used alone, the (Normalized Orthogonal Oja) NOOJA algorithm of [14] is utilized to track minor subspace spanned by the eigenvectors relating to the smallest eigenvalues. But, when associated with an adaptive joint diagonalization, as described in the next section, the NOOJA algorithm tracks the minor subspace spanned by the eigenvectors in the image of C̃+C̅ relating to the smallest eigenvalues which are different from 0. This algorithm extracts adaptively the minor subspace V relating to the input signal r(k)'s autocorrelation matrix R, by maximization the following cost function:
(29)J(V)=E∥r(k)-VVTr(k)∥2=tr(R)-2tr(VTRV)+tr(VTRVVTV).
In our problem, the minor subspace of C̃+C̅ is searched. So, in order to use the NOOJA algorithm, C̃+C̅ must be expressed as C̃+C̅=E[ccT]. But since the function f (cf. (25) and (26)) is not bilinear, the vector c is unknown. Instead of searching the minor subspace of C̃+C̅, the minor subspace of both C̃ and C̅ will be searched. Since C̅ can not be expressed as
(30)C̅=(G⊗G)C̃(G⊗G)-1.
The matricesC̃ and C̅ are similar. They have therefore the same image. In order to find the minor subspace D=(d1,…,dNt) spanned by the eigenvectors and relating to the smallest eigenvalues different from 0 for both C̃ and C̅, we propose to maximize the following cost function:
(31)J(D)=E∥c̃(k)-DDTc̃(k)∥2+E∥c̅(k)-DDTC̅(k)∥2=tr(C̃+C̅)-2tr(DT(C̃+C̅)D)+tr(DT(C̃+C̅)DDTD).
Since C̃ and C̅ have the same image, the first term of the cost function is sufficient to obtain convergence with a constant channel. Using only the matrix C̃ is equivalent to minimize the cost function Simplified CM [5] in analytical way. This simplification reduces the complexity of the adaptive AMMA and we name it the adaptive Analytical Simplified CMA (ASCMA). Please note that the ASCMA batch analysis is the same with the AMMA batch analysis performed in Section 3.5.1. However, the study is solely about C̃ rather than C̃+C̅.
The maximization ofJ(D) may be achieved by using the gradient-descent technique:
(32)D(k+1)=D(k)-β∇J(k),
where
(33)∇J=(-2(C̃+C̅)+(C̃+C̅(k))DDT+DDT(C̃+C̅))D.
The time instantk has been removed to simplify the equations.
By usingDT(k)D(k)=I, the cost function can be simplified as
(34)J(D)=tr(C̃+C̅)-tr(DT(C̃+C̅)D).
By replacing (32) in (34) and by doing the gradient of J(D(k+1)) with respect to β(k), the optimal stepsize βopt(k) is obtained
(35)βopt(k)=tr(DT(k)(C̃(k)+C̅(k))∇J(k))tr(∇JT(k)(C̃(k)+C̅(k))∇J(k)).
Keeping in mind thatD(k) is orthogonal at each iteration, that is, DT(k)D(k)=I, and replacing C̃(k) and C̅(k) by the estimate c̃(k)c̃T(k) and c̅(k)c̅T(k), respectively, we can write
(36)∇J(k)=p̃(k)z̃T(k)+p̅(k)z̅T(k),
where z̃(k)=DT(k)c̃(k),z̅(k)=DT(k)c̅(k),p̃(k)=z̃(k)-D(k)z̃(k), and p̅(k)=z̅(k)-D(k)z̅(k). The optimal step size becomes
(37)βopt=(∥c̃∥2-∥z̃∥2)∥z̃∥2+(∥c̅∥2-∥z̅∥2)∥z̅∥2(∥c̃∥2-∥z̃∥2)2∥z̃∥2+(∥c̅∥2-∥z̅∥2)2∥z̅∥2.
The instant k has been removed in the βopt(k)'s expression in order to simplify the equation.
The updating ofD(k+1) is performed as
(38)D′(k+1)=D(k)-βopt(k)(p̃(k)z̃T(k)+p̅(k)z̅T(k)),
where the matrix D(k+1) is obtained thanks to a Gram-Schmidt orthogonalization of D′(k+1). The algorithm obtained is described in Algorithm 1.
Once the eigenvectors relating to the smallest eigenvalues of bothC̃(k) and C̅(k) are obtained, the constraint dl=t̃l⊗t̃l must be satisfied. In other words, dl must have a Kronecker structure and must be in the image of C̃(k)+C̅(k).Algorithm 1:Adaptive nullspace of C̃and C̅tracking by using NOOJA.
Initialize  D(0)for  k=1,2,…z̃=DT(k)c̃,z̅=DT(k)c̅,p̃=z̃-D(k)z̃,p̅=z̅-V(k)z̅computeβopt by using(37)D′(k+1)=D(k)-βopt(p̃z̃T+p̅z̅T)d1(k+1)=d1'(k+1)/∥d1'(k+1)∥for  q=2 to Nts=[0⋯0]Tfor  l=1 to q-1s=s+(dlT(k+1)dq'(k+1))dl(k+1)end fordq(k+1)=(dq'(k+1)-s)/∥(dq'(k+1)-s)∥end forend forAdaptive Update of the Joint Diagonalization
The method used to satisfy the constraintdl=t̃l⊗t̃l is done as in [9]. T̃∘T̃ is computed and regarded as the current estimate of the subspace basis (D=T̃∘T̃). Using this basis, the subspace update is performed by using the “Normalized Orthogonal OJA” giving D. The last step of this algorithm is the mapping of the columns dl of D to a Kronecker-product:
(39)dl=t̃l⊗t̃l=vec(t̃lt̃lT),vec-1(dl)=Dl=t̃lt̃lT.
Then a power iteration [15] is applied. This iteration can be written as
(40)t̃l(k+1)=Dlt̃l(k).
At the next iteration, the matrix D in the image, which has a Kronecker structure, is obtained with D=T̃∘T̃. Then, D is transmitted at the NOOJA algorithm, allowing to vectors dl obtained in output NOOJA to be in the nullspace and to have a Kronecker struture.Estimated Symbols
In order to obtain the transmitted symbols, the complex matrixT is rebuilt from T̃(41)tlm=t̃lm+jt̃l(m+Nt),l,m∈1,…,Nt.
In order to prevent the extraction of the same source many times, T is reconditioned after each update. The technique used here is the computation of a singular value decomposition of T: T=∑σjujvjH. Then the singular values of T, which are smaller than 0.5, are replaced by 1
(42)recond(T)=∑σj′ujvjH.
Then, the matrixT is applied on the pre-whitened received signals:
(43)ŝ(k)=TH(k)y̲(k).
The suggested adaptive-AMMA algorithm is described in  Algorithm 2.Algorithm 2:Adaptive Analytical-MMA.
Initialize  F,T̃,p̅=0,p̃=0,α=0for  k=1,2,…Update  F, the prewhitened filter,by usingy(k) and [12] we obtain y̲(k)=FHy(k)Update real vectorsỹ(k)=(ℜ(y(k))ℑ(y(k))),y̅(k)=(-ℑ(y(k))ℜ(y(k))),w̃l=(ℜ(w)ℑ(w))Update vector  c̃ and c̅c̃=ỹ̲(k)⊗ỹ̲(k)-p̃/α,c̅=y̲̅(k)⊗y̲̅(k)-p̅/αp̃=λp̃+(1-λ)(ỹ̲(k)⊗ỹ̲(k)),p̅=λp̅+(1-λ)(y̲̅(k)⊗y̲̅(k))α=λα+1-λD=T̃∘T̃RegardD as a basis of the nullspace of C̃ and C̅,and update it usingc̃ and c̅ (Algorithm 1)for  l=1,…,NtDl=vec-1(dl),t̃l=Dlt̃l,t̃l=t̃l/∥t̃l∥Obtain complex  Nt×Nt  matrix  T(k)for  m=1,…,2Nttlm=t̃lm+jt̃l(m+Nt)end forend forT=recond(T),T̃=(ℜ(T)ℑ(T))TEstimate the sourcesŝ(k)=TH(k)y̲(k)end for


## 4. Simulation Results

The channel used in simulations is a time-selective channel; that is, it is time varying but flat in frequency. The transmitted symbols, resulting from a 4-QAM constellation, are sent through a MIMO Rayleigh fading channel. Two transmitters and four receivers are used (Nt=2 and Nr=4). The performances of proposed adaptive-AMMA are compared with those of the adaptive-ASCMA, the adaptive-ACMA and the MMA, in terms of Bit Error Rate (BER) and Signal to Interference plus Noise Ratio (SINR).In the simulation results of the Figure1, we considered the case of a carrier frequency offset (δfTs=10-3) with SNR=30 dB. The equalizer output constellation of the adaptive-ACMA (Figure 1(e)) spins because of the carrier frequency offset. Normally, the MMA (Figure 1(f)) allows to track the carrier frequency offset, but in this case the carrier phase is too high to obtain a good tracking. Only the proposed adaptive-AMMA and adaptive-ASCMA can remove the frequency offset (Figures 1(c) and 1(d)), without the use of a carrier tracking loop.Constellations with a carrier frequency offset.SNR=30 dB, Nt=2, Nr=4, δfTs=10-3,fdTs=0.001.
(a)
Received signals(b)
Transmitted signals(c)
Adaptive-AMMA(d)
Adaptive-ASCMA(e)
Adaptive-ACMA(f)
MMAFigure2 shows the average BER as a function of SNR (with the term “average” being here a time average), and this for two sources and 1000 runs, each one presenting different randomly varying channels and different noise sequences. The value of the Doppler shift fdTs is 1·1e-3. The simulation here was realized without any carrier frequency offset.Figure 2
BER without carrier frequency offset versus SNR.Nt=2,Nr=4,δfTs=0,fdTs=0.0011.The lowest BER on Figure2 is obtained with the proposed adaptive-AMMA algorithm, followed by the adaptive-ASCMA, then by the adaptive-ACMA and the highest BER is obtained with the MMA.So the proposed adaptive-AMMA algorithm allows to obtain a 4 dB gain in SNR terms compared with that obtained by the adaptive-ACMA with a value of BER equal to10-4.Figure3 illustrates the convergence speed of the adaptive-AMMA. The path gain of the channel and its estimate obtained by the adaptive-AMMA are represented on the same figure. We can notice that estimation of adaptive-AMMA follows the variations of the channel.Figure 3
Gain of pathh11 and h22 without carrier frequency offset versus symbols number. Nt=2,δfTs=0,fdTs=0.0033 and SNR=30 dB.From the results of the Figure4, we can see that the performances of the proposed adaptive-AMMA and the adaptive-ASCMA are very closed in residual SINR terms. Moreover, these two algorithms perform higher residual SINR than the adaptive-ACMA and the MMA.Figure 4
SINR without carrier frequency offset versus symbols number.Nt=2,Nr=4,δfTs=0,fdTs=0.0033 and SNR=15 dB.


## 5. Conclusion

In this paper, we have proposed two new BSS algorithms, named adaptive-AMMA and adaptive-ASCMA, build on the Multimodulus criteria, the Simplified Constant Modulus respectively, and an analytical processing method. The second algorithm is a simplified version of the adaptive-AMMA and it accomplishes performances results closed to the adaptive-AMMA. These algorithms perform better results in BER and SINR terms compared with the adaptive-ACMA and the MMA algorithms. In addition to the blind separation, they can perform a carrier phase recovery simultaneously unlike the adaptive-ACMA.



---

*Source: 307927-2010-02-02.xml*