# Blind Separation of Cyclostationary Sources Sharing Common Cyclic Frequencies Using Joint Diagonalization Algorithm

Authors: Amine Brahmi; Hicham Ghennioui; Christophe Corbier; François Guillet; M’hammed Lahbabi
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We propose a new method for blind source separation of cyclostationary sources, whose cyclic frequencies are unknown and may share one or more common cyclic frequencies. The suggested method exploits the cyclic correlation function of observation signals to compose a set of matrices which has a particular algebraic structure. The aforesaid matrices are automatically selected by proposing two new criteria. Then, they are jointly diagonalized so as to estimate the mixing matrix and retrieve the source signals as a consequence. The nonunitary joint diagonalization (NU-JD) is ensured by Broyden-Fletcher-Goldfarb-Shanno (BFGS) method which is the most commonly used update strategy for implementing a quasi-Newton technique. The efficiency of the method is illustrated by numerical simulations in digital communications context, which show good performances comparing to other state-of-the-art methods.


---

## Body

## 1. Introduction

The target of blind source separation (BSS) is to retrieve the input signals called sources from their mixtures coming to multiple sensors without any preknowledge about the mixing process. BSS is a major problem of signal processing which has been addressed in the last three decades (see [1] for a review). In literature, many approaches have been developed in order to figure out this issue using statistics of second and fourth order, namely, Second-Order Blind Identification (SOBI) [2] and Joint Approximate Diagonalization of Eigen matrices (JADE) [3]. These approaches have proved to establish some limitations in a wide scope of practical situations where the source signals are nonstationary and very often cyclostationary such as radiocommunications, telemetry, radar applications, and mechanics [4]. In fact, according to Ferreol in [5], the stationarity assumption of source signals performs ineffectively BSS problem.Cyclostationarity is a subclass of nonstationarity which distinguishes stochastic processes whose statistics change periodically with time. Thus, it is necessary to consider cyclostationarity to perform BSS. Many methods have been proposed to blindly achieve the separation for cyclostationary sources. Brahmi et al. have solved the problem of blind identification of FIR MIMO systems driven by cyclostationary inputs whose cyclic frequencies are pairwise distinct using joint block diagonalization based on BFGS method in [6]. Liang et al. in [7] use the information provided from the cyclic frequencies so as to separate source signals. Abed-Meraim et al. [8] address the problem of BSS assuming that the source signals are cyclostationary based on an iterative algorithm using the cyclic correlation function of observation signals. This method is useful when each source signal has only one cyclic frequency and the number of the source signals which share a common cyclic frequency is known. Ghennioui et al. [9] have proposed a new approach combining a nonunitary joint diagonalization algorithm to a general automatic matrices selection procedure for the case of unknown and different cyclic frequencies. Ghaderi et al. present in [10] a method for blind source extraction of cyclostationary sources, whose cyclic frequencies are known and share some common ones. Jafari et al. in [11, 12] propose an adaptive blind source separation algorithm for the separation of convolutive mixtures of cyclostationary signals based on natural gradient algorithm. The aforementioned algorithm requires estimating the cycle frequencies of source signals. Boustany and Antoni [13] propose a method for blind extraction of one cyclostationay signal using a subspace decomposition of the observation signals via their cyclic statistics. Capdessus et al. [14] propose an algorithm for the extraction a signal of interest which is cyclostationary one and its cyclic frequency is a priori known. This algorithm relies on second-order statistics of the observation signals. Rhioui et al. propose in [15] a method for the mixing matrix identification for underdetermined mixtures of cyclostationary signals with different cyclic frequencies. Jallon and Chevreuil [16] have come up with a justification for using the common algorithm for the cyclostationary context despite the fact that it has been originally developed for the stationary one. Pham [17] has proposed a new approach based on joint diagonalization of a set of cyclic spectral density of observation matrices. The two last approaches are addressed in the simplest mixture model (noise-free data). Despite the fact that these algorithms are successful under assumed conditions, they have diverse limitations, since, in front of real situations, the cyclic frequencies in most of cases are unknown and may be shared by source signals.The main purpose of this work is to perform blind separation of instantaneous mixtures of cyclostationary source signals which may share one or more common cyclic frequencies whose preknowledge is not needed. By exploiting the particular structure of cyclic correlation matrices of source signals, we show that the considered problem can be rephrased as a problem of joint diagonalization of matrices that have been automatically selected using a new procedure. Then, the joint diagonalization algorithm based on BFGS method [18] is applied on this set of matrices. The rest of this paper is organized as follows. Section 2 formulates the problem of interest. Section 3 puts forward some theoretical preliminaries related to cyclostationarity and NU-JD. Section 4 describes the proposed method for BSS. In Section 5, the performances of the proposed method are numerically evaluated and compared with other existing methods in digital telecommunications context. Finally, in the Section 6, conclusions are drawn.


## 2. Problem Statement

The BSS problem can be modelled in a simple linear instantaneous mixture ofn emitted source signals that are received by m sensors (see Figure 1). The relationship of mixing system is given by(1)xt=Ast+bt,where x(t)∈Cm is the observations vector, s(t)∈Cn is the vector of unknown source signals, b(t)∈Cm is the additive noise vector, and A∈Cm×n is the unknown mixing matrix. The purpose of BSS is to find an estimate A~ of A and retrieve s(t) from x(t) only, as(2)s~t=A~#xt=PΔdiag⁡e-jΦst,where j2=-1, ·# denotes pseudoinverse matrix, P is a permutation matrix (corresponding to an arbitrary order of restitution of the sources), Δ is a diagonal matrix (corresponding to arbitrary scaling for the recovered sources), Φ=ϕ1,…,ϕnT, ∀ϕi∈R represents the phase vector (corresponding to phase shift ambiguity in complex domain of the source signals), and diag⁡(a) is square diagonal matrix containing the elements of the vector a. Thus, one should look for a separating matrix A~# such that(3)A~#A=PΔdiag⁡e-jΦ.Figure 1
The relationship between sources and observations.The following assumptions are held in this paper:(A1)
The mixing matrixA has full column rank (m≥n).(A2)
The source signals are zero-mean, cyclostationary, and mutually uncorrelated and may share one or more common cyclic frequencies.(A3)
The noise signalsb(t) are stationary, zero-mean random signals, and mutually uncorrelated with the source signals.


## 3. Some Recalls

### 3.1. Cyclostationary Signals

Cyclostationary signals are frequent random signals with statistical parameters that alter in periodic manner in time. A coupling of stationary random process with a deterministic periodic one gives rise to a cyclostationary process (see Figure2).Figure 2
Cyclostationary process.A time process{x(t)} is first-order cyclostationary, if its first-order moment is T periodic:(4)mxt=Ext=Ext+T.{x(t)} is second-order cyclostationary, if its second-order moments are T periodic:(5)Rxt,τ=Ext+τ2x∗t-τ2=Rxt,τ+T,where E{·} is the mathematical expectation and the (∗) stands for complex conjugation. In this case, the correlation function is decomposed into Fourier series as follows:(6)Rxt,τ=∑ξRxξτej2πξt,with(7)Rxξτ=limT→∞⁡1T∫-T/2T/2Rxt,τe-j2πξtdt,where Rxξ(τ) is the cyclic autocorrelation function with τ as the time lag and ψ={ξ=q/T,q∈Z} as the cyclic frequencies set. In the frequency domain, cyclostationary signals are characterised by the spectral correlation density (SCD), which is nothing but the Fourier transform of the cyclic autocorrelation function defined as(8)Sxξf=∑τRxξτe-j2πfτ.We note a spectral repetition of cyclostationary signals due to the correlation between their spectral components at a specified frequency from each other which is the cyclic frequency.


### 3.2. Nonunitary Joint Diagonalization of Matrices

We consider a set ofK(K≥2) real or complex square matrices sized (m×m), M={Λ1,…,ΛK} which are decomposed such that(9)Λi=ADiAH,orDi=BΛiBH,where ·H stands for the transpose conjugate operator. The matrices Di sized (n×n) are diagonal ones, A is a (m×m) full column rank matrix, and B is its Moore-Penrose pseudoinverse. The NU-JD problem consists of estimating the matrix A from only the matrix set M. It is worth noting that, in practice, the set of matrices is built by some sample estimated statistics that are corrupted by estimation errors due to noise and finite sample size effects. Thus, they are only approximately jointly diagonalizable. The NU-JD problem could be solved by considering the following cost function introduced in [19]:(10)FjdB=∑i=1KOffdiag⁡BΛiBHF2,where ·F2 stands for the Frobenius norm and Offdiag⁡· denotes the zero-diagonal operator (see Appendix).


## 3.1. Cyclostationary Signals

Cyclostationary signals are frequent random signals with statistical parameters that alter in periodic manner in time. A coupling of stationary random process with a deterministic periodic one gives rise to a cyclostationary process (see Figure2).Figure 2
Cyclostationary process.A time process{x(t)} is first-order cyclostationary, if its first-order moment is T periodic:(4)mxt=Ext=Ext+T.{x(t)} is second-order cyclostationary, if its second-order moments are T periodic:(5)Rxt,τ=Ext+τ2x∗t-τ2=Rxt,τ+T,where E{·} is the mathematical expectation and the (∗) stands for complex conjugation. In this case, the correlation function is decomposed into Fourier series as follows:(6)Rxt,τ=∑ξRxξτej2πξt,with(7)Rxξτ=limT→∞⁡1T∫-T/2T/2Rxt,τe-j2πξtdt,where Rxξ(τ) is the cyclic autocorrelation function with τ as the time lag and ψ={ξ=q/T,q∈Z} as the cyclic frequencies set. In the frequency domain, cyclostationary signals are characterised by the spectral correlation density (SCD), which is nothing but the Fourier transform of the cyclic autocorrelation function defined as(8)Sxξf=∑τRxξτe-j2πfτ.We note a spectral repetition of cyclostationary signals due to the correlation between their spectral components at a specified frequency from each other which is the cyclic frequency.


## 3.2. Nonunitary Joint Diagonalization of Matrices

We consider a set ofK(K≥2) real or complex square matrices sized (m×m), M={Λ1,…,ΛK} which are decomposed such that(9)Λi=ADiAH,orDi=BΛiBH,where ·H stands for the transpose conjugate operator. The matrices Di sized (n×n) are diagonal ones, A is a (m×m) full column rank matrix, and B is its Moore-Penrose pseudoinverse. The NU-JD problem consists of estimating the matrix A from only the matrix set M. It is worth noting that, in practice, the set of matrices is built by some sample estimated statistics that are corrupted by estimation errors due to noise and finite sample size effects. Thus, they are only approximately jointly diagonalizable. The NU-JD problem could be solved by considering the following cost function introduced in [19]:(10)FjdB=∑i=1KOffdiag⁡BΛiBHF2,where ·F2 stands for the Frobenius norm and Offdiag⁡· denotes the zero-diagonal operator (see Appendix).


## 4. Suggested Method

### 4.1. Our Approach

From (1), it can be easily shown that the cyclic correlation matrix of observation signals x(t) has the following decomposition:(11)Rxξτ=ARsξτAH+Rbξτ,where ·H is the conjugate transpose operator and Rsξ(τ) (resp., Rbξ(τ)) is the cyclic correlation matrix of source signals (resp., noise signals). The expression of Rbξ(τ) can be further developed. Knowing that the noise is stationary and using (7), we have(12)Rbξτ=Rbτδξ,where(13)δξ=limT→∞⁡1T∫-T/2T/2e-j2πξtdt.This implies that(14)δξ=1ifξ=00else.Therefore, if ξ≠0, then(15)Rxξτ=ARsξτAH.Practically, the matrix Rsξ(τ) has one of the following structures:(a)
If the source signals have pairwise distinct cyclic frequencies, thenRsξ(τ) is a diagonal matrix with only one nonnull element corresponding to the ith source at ξi (for this kind of structure, using Rxξ(τ), a detection procedure has been proposed in [20]).(b)
If the source signals share one or more cyclic frequencies, thenRsξ(τ) at these shared frequencies is a diagonal matrix with n nonnull elements which is our case of interest.Thus, we propose to diagonalize jointly the set of cyclic correlation matrices of observation signals Rxξ(τ) at different time lags and cyclic frequencies (ξ≠0). The question now being asked is the following: how to compose the matrices set Mjd to be joint diagonalized?


### 4.2. New Criteria for Set Construction

It has to be noted that the preknowledge of the cyclic frequencies of source signalss(t) is not required even if such a thing facilitates the resolution of BSS issue. As a matter of fact, we could take advantage of the particular algebraic structure of the matrix Rsξ(τ) to perform BSS. When the cyclic frequencies are unknown, we calculate Rxξ(τ) for an enough number of frequency bins in order to ensure sweeping almost all the cyclic frequencies of source signals; then, we use a new detection procedure to select the matrices Rxξ(τ) which correspond to the structure (b). In other words, the procedure has to detect matrices Rxbξ(τ)=WRxξ(τ)WH (W is n×m whitening matrix such that WAAHWH=In) which have the same number n of eigenvalues as Rsξ(τ) since all source signals are present at a given common cyclic frequency. One possible way to blindly access to diagonal terms of Rsξ(τ) consists in computing the eigenvalues of Rxbξ(τ). In fact, using the whitening matrix definition in [3], we have(16)eig⁡Rxbξτ=eig⁡WARsξτAHWH=eig⁡Rsξτ=θ1,…,θnT,where the vector eig⁡(M) contains the eigenvalues of the square matrix M, WA is an unitary matrix, and θi, for i=1,…,n are the eigenvalues of Rsξ(τ). Therefore, we propose the following new criteria:(17)C1=∑i=1nθi2RxbξτF2≥1-ε,where ε is positive constant. One can note using the invariance property of the Frobenius norm under an unitary transformation that Rxbξ(τ)F2=Rsξ(τ)F2. Ideally C1 equals 1. However, in practice, the matrices Rsξ(τ) can never be strictly diagonal. Thus, matrices Rxbξ(τ) should be selected as C1≥1-ε with ε close to zero. Furthermore, in order to detect the matrices where all source signals are present in main diagonal of Rsξ(τ) and to avoid the low energy matrices Rxbξ(τ), we add to C1:(18)C2=det⁡Rxbξτ=det⁡Rsξτ≥η,where det⁡(·) denotes the matrix determinant and η is a small positive constant. Finally, if a given Rxbξ(τ) satisfies C1 and C2, then it is retained. Once the matrices set is built, it is directly joint diagonalized. The set size of the selected matrices is related to the choice of the threshold values of ε and η and affects directly the separation quality. However, when the cyclic frequencies are a priori known, the operation is much more easier; it is reduced to computing Rxξ(τ) at different time lags for each cyclic frequency {ξi/i=1,…,n} then to diagonalize simultaneously the set built(19)Mjd=Rxξiτ1,…,Rxξiτmax.The joint diagonalization is ensured by the BFGS method which will be detailed in the following subsection.


### 4.3. BFGS Based NU-JD Algorithm

We use the coming cost function to figure out the NU-JD problem:(20)FjdB=∑i=1KOffdiag⁡BRxξiτiBHF2,where Rxξi(τi) is the ith matrix which belongs to the set Mjd to be joint diagonalized. We propose an approach based on a BFGS method with an exact computation of the optimal step-size. It estimates the joint diagonalizer matrix B∈Cn×m by minimizing problem given in (20). The BFGS method requires the gradient and the Hessian of the cost function to be computed at each iteration.


#### 4.3.1. Algorithm Principle

From initial guessesB(0) and He(0) and a given number of iterations niter, the following instructions are iterated until B(k) converges to the solution.(S1) Look for a search direction d(k-1) by solving (21)dk-1=-Hek-1-1∇aFjdBk-1.(S2) Perform a line search to find the optimal step-size α (positive, a small enough number) in the previous direction (found in the first step).(S3) Update(22)Bk=Bk-1+αdk-1.(S4) Set(23)sk-1=αdk-1;yk-1=∇aFjdBk-∇aFjdBk-1.(S5) Using the rank-one updates using gradient evaluations, the Hessian matrix is estimated as follows:(24)Hek=Hek-1+yk-1yk-1Tyk-1Tsk-1-Hek-1sk-1sk-1THek-1sk-1THek-1sk-1,k=1,…,niter.·-1 and ·T denote the inverse and transpose of a matrix, respectively, ∇aFjd(B) is the complex absolute gradient matrix of the cost function given in (20) which is defined as (see [21]) (25)∇aFjdB=2∂FjdB∂B∗,where B∗ stands for the complex conjugate of the complex matrix B and ∂/∂B∗ is the partial derivative operator. ∇aFjd(B) was calculated earlier in [19] and found to be equal to (26)∇aFjdB=2∑i=1KOffdiag⁡BRxξiτiBHBRxξiτiH+2∑i=1KOffdiag⁡BRxξiτiBHHBRxξiτi,where K is the number of selected matrices to be block joint diagonalized. It is worth mentioning that He(0) can be initialized with the identity matrix and B(0) has to be full-rank matrix chosen differently from the zero matrix as it is a trivial solution of (20).


#### 4.3.2. Seek the Optimal Step-Size

It is taken for granted that finding a good step-sizeα in search direction is critical issue for decreasing the total number of iterations to reach convergence. The enhanced line search consists of minimization Fjd(Bk) with respect to α. For simplicity, we prefer to hide the dependency upon the iteration k. It is a matter of standard algebraic manipulations to show that Fjd(B(k)) is a 4th-degree polynomial in α whose expression is(27)FjdB=aα4+bα3+cα2+dα+e,where its coefficients are given by(28)a=∑i=1KtrM3Offdiag⁡M3H,b=-∑i=1KtrM0Offdiag⁡M3H-trM3M2,c=∑i=1KtrM1HOffdiag⁡M3H+M3Offdiag⁡M1+trM0M2,d=-∑i=1KtrM1HM2-trM0Offdiag⁡M1,e=∑i=1KtrM1HOffdiag⁡M1,M0=RxξiτiHΓHB+RxξiτiHBHΓ,M1=BRxξiτiBH,M2=Offdiag⁡ΓRxξiτiBH+Offdiag⁡BRxξiτiΓH,M3=ΓRxξiτiHΓH,with tr· denoting the trace operator and Γ=He-1∇aFjd(B). The optimal α can be found by polynomial rooting of the derivative third-order polynomial, namely, by solving with respect to α, (29)∂FjdB∂α=0⟺4aα3+3bα2+2cα+d=0,to which there is three roots, the true minimum can be found by substituting each root back into the polynomial given in (27) and selecting the solution that gives the littlest value.


#### 4.3.3. Algorithmic Complexity

For the gradient matrix whose expression is given in (26), the computational cost is given by Table 1. The computational cost of the optimal step-size is ruled by the computation of the five coefficients of the 4th-degree polynomial; its amount is shown in Table 2. Therefore, the computational cost of BFGS algorithm is given by Table 3. Finally, notice that the global complexity of the BFGS algorithm has to be multiplied by the overall iterations number Ni required to attain the convergence. In practical applications, the computational time necessary to build the set of the K matrices should be counted too.Table 1
The gradient matrix computational cost.
Case Cost per iteration m ≠ n 4 K n m ( m + n ) + 2 K n 2 m ≫ n 4 K n m 2 m = n 8 K n 3 + 2 K n 2Table 2
The computational cost of the optimal step-size.
Case Cost per iteration m ≠ n 24 K m n ( m + n ) + 9 K n 2 ( 1 + n ) m ≫ n 24 K m 2 n m = n 57 K n 3 + 9 K n 2Table 3
The BFGS algorithm computational cost.
Case Cost per iteration m ≠ n 24 K m n ( m + n ) + 9 K n 2 ( 1 + n ) m ≫ n 24 K m 2 n m = n 57 K n 3 + 9 K n 2


### 4.4. Summary

The principle proposed is summarized in Algorithm1.Algorithm 1:Summary of the proposed method.
Set building: (i) Detect the useful matricesRXξ(τ) using the selection procedure; (ii) Compute the set matricesMjd for each cyclic frequency ξi: M j d = { R x ξ i ( τ 1 ) , … , R x ξ i ( τ m a x ) }; Initialization: (i) Given an initial guessB(0) and an approximate Hessian He(0); (ii) Given the number of iterationsniter; (iii) Given a small positive thresholdϵ; for k ← 1 , n i t e r do Compute∇aFjd(B) whose expression is given by Eq. (26); Perform a line search to find the optimal step-sizeα; SetB(k)←B(k-1)-αHe(k-1)-1∇aFjd(B(k-1)); ComputeHe(k) whose expression is given by Eq. (24); if|B(k)-B(k-1)|≤ϵthen break; end if end for


## 4.1. Our Approach

From (1), it can be easily shown that the cyclic correlation matrix of observation signals x(t) has the following decomposition:(11)Rxξτ=ARsξτAH+Rbξτ,where ·H is the conjugate transpose operator and Rsξ(τ) (resp., Rbξ(τ)) is the cyclic correlation matrix of source signals (resp., noise signals). The expression of Rbξ(τ) can be further developed. Knowing that the noise is stationary and using (7), we have(12)Rbξτ=Rbτδξ,where(13)δξ=limT→∞⁡1T∫-T/2T/2e-j2πξtdt.This implies that(14)δξ=1ifξ=00else.Therefore, if ξ≠0, then(15)Rxξτ=ARsξτAH.Practically, the matrix Rsξ(τ) has one of the following structures:(a)
If the source signals have pairwise distinct cyclic frequencies, thenRsξ(τ) is a diagonal matrix with only one nonnull element corresponding to the ith source at ξi (for this kind of structure, using Rxξ(τ), a detection procedure has been proposed in [20]).(b)
If the source signals share one or more cyclic frequencies, thenRsξ(τ) at these shared frequencies is a diagonal matrix with n nonnull elements which is our case of interest.Thus, we propose to diagonalize jointly the set of cyclic correlation matrices of observation signals Rxξ(τ) at different time lags and cyclic frequencies (ξ≠0). The question now being asked is the following: how to compose the matrices set Mjd to be joint diagonalized?


## 4.2. New Criteria for Set Construction

It has to be noted that the preknowledge of the cyclic frequencies of source signalss(t) is not required even if such a thing facilitates the resolution of BSS issue. As a matter of fact, we could take advantage of the particular algebraic structure of the matrix Rsξ(τ) to perform BSS. When the cyclic frequencies are unknown, we calculate Rxξ(τ) for an enough number of frequency bins in order to ensure sweeping almost all the cyclic frequencies of source signals; then, we use a new detection procedure to select the matrices Rxξ(τ) which correspond to the structure (b). In other words, the procedure has to detect matrices Rxbξ(τ)=WRxξ(τ)WH (W is n×m whitening matrix such that WAAHWH=In) which have the same number n of eigenvalues as Rsξ(τ) since all source signals are present at a given common cyclic frequency. One possible way to blindly access to diagonal terms of Rsξ(τ) consists in computing the eigenvalues of Rxbξ(τ). In fact, using the whitening matrix definition in [3], we have(16)eig⁡Rxbξτ=eig⁡WARsξτAHWH=eig⁡Rsξτ=θ1,…,θnT,where the vector eig⁡(M) contains the eigenvalues of the square matrix M, WA is an unitary matrix, and θi, for i=1,…,n are the eigenvalues of Rsξ(τ). Therefore, we propose the following new criteria:(17)C1=∑i=1nθi2RxbξτF2≥1-ε,where ε is positive constant. One can note using the invariance property of the Frobenius norm under an unitary transformation that Rxbξ(τ)F2=Rsξ(τ)F2. Ideally C1 equals 1. However, in practice, the matrices Rsξ(τ) can never be strictly diagonal. Thus, matrices Rxbξ(τ) should be selected as C1≥1-ε with ε close to zero. Furthermore, in order to detect the matrices where all source signals are present in main diagonal of Rsξ(τ) and to avoid the low energy matrices Rxbξ(τ), we add to C1:(18)C2=det⁡Rxbξτ=det⁡Rsξτ≥η,where det⁡(·) denotes the matrix determinant and η is a small positive constant. Finally, if a given Rxbξ(τ) satisfies C1 and C2, then it is retained. Once the matrices set is built, it is directly joint diagonalized. The set size of the selected matrices is related to the choice of the threshold values of ε and η and affects directly the separation quality. However, when the cyclic frequencies are a priori known, the operation is much more easier; it is reduced to computing Rxξ(τ) at different time lags for each cyclic frequency {ξi/i=1,…,n} then to diagonalize simultaneously the set built(19)Mjd=Rxξiτ1,…,Rxξiτmax.The joint diagonalization is ensured by the BFGS method which will be detailed in the following subsection.


## 4.3. BFGS Based NU-JD Algorithm

We use the coming cost function to figure out the NU-JD problem:(20)FjdB=∑i=1KOffdiag⁡BRxξiτiBHF2,where Rxξi(τi) is the ith matrix which belongs to the set Mjd to be joint diagonalized. We propose an approach based on a BFGS method with an exact computation of the optimal step-size. It estimates the joint diagonalizer matrix B∈Cn×m by minimizing problem given in (20). The BFGS method requires the gradient and the Hessian of the cost function to be computed at each iteration.


### 4.3.1. Algorithm Principle

From initial guessesB(0) and He(0) and a given number of iterations niter, the following instructions are iterated until B(k) converges to the solution.(S1) Look for a search direction d(k-1) by solving (21)dk-1=-Hek-1-1∇aFjdBk-1.(S2) Perform a line search to find the optimal step-size α (positive, a small enough number) in the previous direction (found in the first step).(S3) Update(22)Bk=Bk-1+αdk-1.(S4) Set(23)sk-1=αdk-1;yk-1=∇aFjdBk-∇aFjdBk-1.(S5) Using the rank-one updates using gradient evaluations, the Hessian matrix is estimated as follows:(24)Hek=Hek-1+yk-1yk-1Tyk-1Tsk-1-Hek-1sk-1sk-1THek-1sk-1THek-1sk-1,k=1,…,niter.·-1 and ·T denote the inverse and transpose of a matrix, respectively, ∇aFjd(B) is the complex absolute gradient matrix of the cost function given in (20) which is defined as (see [21]) (25)∇aFjdB=2∂FjdB∂B∗,where B∗ stands for the complex conjugate of the complex matrix B and ∂/∂B∗ is the partial derivative operator. ∇aFjd(B) was calculated earlier in [19] and found to be equal to (26)∇aFjdB=2∑i=1KOffdiag⁡BRxξiτiBHBRxξiτiH+2∑i=1KOffdiag⁡BRxξiτiBHHBRxξiτi,where K is the number of selected matrices to be block joint diagonalized. It is worth mentioning that He(0) can be initialized with the identity matrix and B(0) has to be full-rank matrix chosen differently from the zero matrix as it is a trivial solution of (20).


### 4.3.2. Seek the Optimal Step-Size

It is taken for granted that finding a good step-sizeα in search direction is critical issue for decreasing the total number of iterations to reach convergence. The enhanced line search consists of minimization Fjd(Bk) with respect to α. For simplicity, we prefer to hide the dependency upon the iteration k. It is a matter of standard algebraic manipulations to show that Fjd(B(k)) is a 4th-degree polynomial in α whose expression is(27)FjdB=aα4+bα3+cα2+dα+e,where its coefficients are given by(28)a=∑i=1KtrM3Offdiag⁡M3H,b=-∑i=1KtrM0Offdiag⁡M3H-trM3M2,c=∑i=1KtrM1HOffdiag⁡M3H+M3Offdiag⁡M1+trM0M2,d=-∑i=1KtrM1HM2-trM0Offdiag⁡M1,e=∑i=1KtrM1HOffdiag⁡M1,M0=RxξiτiHΓHB+RxξiτiHBHΓ,M1=BRxξiτiBH,M2=Offdiag⁡ΓRxξiτiBH+Offdiag⁡BRxξiτiΓH,M3=ΓRxξiτiHΓH,with tr· denoting the trace operator and Γ=He-1∇aFjd(B). The optimal α can be found by polynomial rooting of the derivative third-order polynomial, namely, by solving with respect to α, (29)∂FjdB∂α=0⟺4aα3+3bα2+2cα+d=0,to which there is three roots, the true minimum can be found by substituting each root back into the polynomial given in (27) and selecting the solution that gives the littlest value.


### 4.3.3. Algorithmic Complexity

For the gradient matrix whose expression is given in (26), the computational cost is given by Table 1. The computational cost of the optimal step-size is ruled by the computation of the five coefficients of the 4th-degree polynomial; its amount is shown in Table 2. Therefore, the computational cost of BFGS algorithm is given by Table 3. Finally, notice that the global complexity of the BFGS algorithm has to be multiplied by the overall iterations number Ni required to attain the convergence. In practical applications, the computational time necessary to build the set of the K matrices should be counted too.Table 1
The gradient matrix computational cost.
Case Cost per iteration m ≠ n 4 K n m ( m + n ) + 2 K n 2 m ≫ n 4 K n m 2 m = n 8 K n 3 + 2 K n 2Table 2
The computational cost of the optimal step-size.
Case Cost per iteration m ≠ n 24 K m n ( m + n ) + 9 K n 2 ( 1 + n ) m ≫ n 24 K m 2 n m = n 57 K n 3 + 9 K n 2Table 3
The BFGS algorithm computational cost.
Case Cost per iteration m ≠ n 24 K m n ( m + n ) + 9 K n 2 ( 1 + n ) m ≫ n 24 K m 2 n m = n 57 K n 3 + 9 K n 2


## 4.3.1. Algorithm Principle

From initial guessesB(0) and He(0) and a given number of iterations niter, the following instructions are iterated until B(k) converges to the solution.(S1) Look for a search direction d(k-1) by solving (21)dk-1=-Hek-1-1∇aFjdBk-1.(S2) Perform a line search to find the optimal step-size α (positive, a small enough number) in the previous direction (found in the first step).(S3) Update(22)Bk=Bk-1+αdk-1.(S4) Set(23)sk-1=αdk-1;yk-1=∇aFjdBk-∇aFjdBk-1.(S5) Using the rank-one updates using gradient evaluations, the Hessian matrix is estimated as follows:(24)Hek=Hek-1+yk-1yk-1Tyk-1Tsk-1-Hek-1sk-1sk-1THek-1sk-1THek-1sk-1,k=1,…,niter.·-1 and ·T denote the inverse and transpose of a matrix, respectively, ∇aFjd(B) is the complex absolute gradient matrix of the cost function given in (20) which is defined as (see [21]) (25)∇aFjdB=2∂FjdB∂B∗,where B∗ stands for the complex conjugate of the complex matrix B and ∂/∂B∗ is the partial derivative operator. ∇aFjd(B) was calculated earlier in [19] and found to be equal to (26)∇aFjdB=2∑i=1KOffdiag⁡BRxξiτiBHBRxξiτiH+2∑i=1KOffdiag⁡BRxξiτiBHHBRxξiτi,where K is the number of selected matrices to be block joint diagonalized. It is worth mentioning that He(0) can be initialized with the identity matrix and B(0) has to be full-rank matrix chosen differently from the zero matrix as it is a trivial solution of (20).


## 4.3.2. Seek the Optimal Step-Size

It is taken for granted that finding a good step-sizeα in search direction is critical issue for decreasing the total number of iterations to reach convergence. The enhanced line search consists of minimization Fjd(Bk) with respect to α. For simplicity, we prefer to hide the dependency upon the iteration k. It is a matter of standard algebraic manipulations to show that Fjd(B(k)) is a 4th-degree polynomial in α whose expression is(27)FjdB=aα4+bα3+cα2+dα+e,where its coefficients are given by(28)a=∑i=1KtrM3Offdiag⁡M3H,b=-∑i=1KtrM0Offdiag⁡M3H-trM3M2,c=∑i=1KtrM1HOffdiag⁡M3H+M3Offdiag⁡M1+trM0M2,d=-∑i=1KtrM1HM2-trM0Offdiag⁡M1,e=∑i=1KtrM1HOffdiag⁡M1,M0=RxξiτiHΓHB+RxξiτiHBHΓ,M1=BRxξiτiBH,M2=Offdiag⁡ΓRxξiτiBH+Offdiag⁡BRxξiτiΓH,M3=ΓRxξiτiHΓH,with tr· denoting the trace operator and Γ=He-1∇aFjd(B). The optimal α can be found by polynomial rooting of the derivative third-order polynomial, namely, by solving with respect to α, (29)∂FjdB∂α=0⟺4aα3+3bα2+2cα+d=0,to which there is three roots, the true minimum can be found by substituting each root back into the polynomial given in (27) and selecting the solution that gives the littlest value.


## 4.3.3. Algorithmic Complexity

For the gradient matrix whose expression is given in (26), the computational cost is given by Table 1. The computational cost of the optimal step-size is ruled by the computation of the five coefficients of the 4th-degree polynomial; its amount is shown in Table 2. Therefore, the computational cost of BFGS algorithm is given by Table 3. Finally, notice that the global complexity of the BFGS algorithm has to be multiplied by the overall iterations number Ni required to attain the convergence. In practical applications, the computational time necessary to build the set of the K matrices should be counted too.Table 1
The gradient matrix computational cost.
Case Cost per iteration m ≠ n 4 K n m ( m + n ) + 2 K n 2 m ≫ n 4 K n m 2 m = n 8 K n 3 + 2 K n 2Table 2
The computational cost of the optimal step-size.
Case Cost per iteration m ≠ n 24 K m n ( m + n ) + 9 K n 2 ( 1 + n ) m ≫ n 24 K m 2 n m = n 57 K n 3 + 9 K n 2Table 3
The BFGS algorithm computational cost.
Case Cost per iteration m ≠ n 24 K m n ( m + n ) + 9 K n 2 ( 1 + n ) m ≫ n 24 K m 2 n m = n 57 K n 3 + 9 K n 2


## 4.4. Summary

The principle proposed is summarized in Algorithm1.Algorithm 1:Summary of the proposed method.
Set building: (i) Detect the useful matricesRXξ(τ) using the selection procedure; (ii) Compute the set matricesMjd for each cyclic frequency ξi: M j d = { R x ξ i ( τ 1 ) , … , R x ξ i ( τ m a x ) }; Initialization: (i) Given an initial guessB(0) and an approximate Hessian He(0); (ii) Given the number of iterationsniter; (iii) Given a small positive thresholdϵ; for k ← 1 , n i t e r do Compute∇aFjd(B) whose expression is given by Eq. (26); Perform a line search to find the optimal step-sizeα; SetB(k)←B(k-1)-αHe(k-1)-1∇aFjd(B(k-1)); ComputeHe(k) whose expression is given by Eq. (24); if|B(k)-B(k-1)|≤ϵthen break; end if end for


## 5. Numerical Simulations

We present simulations to illustrate the effectiveness of the proposed method in the BSS context. We consider two amplitude-modulated source signals defined as(30)sit=∑k∈Zaikgt-kTicos⁡2πfit+ϕi,where ai=1,2(k) are i.i.d and zero-mean random binary sequences, Ti=1,2 represent the symbol periods, they are, respectively, equal to T1=10 and T2=8, fi=1,2=0.2 are the normalized carrier frequencies, ϕi=1,2 are the carrier phases, they, respectively, equal ϕ1=π/6 and ϕ2=π/8, and g(n) is a triangular waveform such that(31)gn=2Tnif0≤n≤T2-2Tn+2ifT2+1≤n≤T-10else.The source signals share the cyclic frequencies f1=f2=0.2, 5f1=10/T1=8/T2=1, and their multiples. The mixing matrices were randomly generated. We use A1 in the first two simulations and A2 in the last one:(32)A1=0.55990.40700.28990.84310.94150.3397,A2=0.84000.99710.70350.1965.The signal-to-noise ratio (SNR) is computed as SNR=-10log10⁡(σb2). We assume that the source cyclic frequencies are unknown and consequently we use the detection procedure described in (17) and (18) with ε=10-3 and η=10-1 for all simulations to construct the set of matrices to be joint diagonalized which correspond to the source cyclic correlation matrices Rsξ(τ) with a diagonal structure. In order to evaluate the quality of the estimation, one can measure the Moreau-Amari index presented in [22] defined as(33)Iperf=1nn-1∑i=1n∑j=1nzi,jF2maxl⁡zi,lF2-1+∑j=1n∑i=1nzi,jF2maxl⁡zl,jF2-1,where zi,j is the i,jth element of Z=BA. The closer to zero in linear scale (-∞ in logarithmic scale), the higher the separation accuracy. Regarding the charts, Iperf(·) is given in decibel and estimated by averaging 100 independent trials. The proposed method JDBFGS is compared to two types of methods. The first one is based on unitary joint diagonalization algorithms: JADE, while the second one is based on the nonunitary joint diagonalization using the optimal gradient descent JDGRAD [9]. Figures 3 and 4 show that the proposed algorithm takes the lead ahead of its competitors in all two cases of the mixing matrix (m>n or m=n). In addition, one can notice in Figure 3 that when the number of selected matrices to be joint diagonalized gets bigger, the performances are better, while the computational cost steps up. The performances clearly increase when the signal-to-noise ratio gets higher.Figure 3
Performance index versusSNR (a) and versus number of selected matrices (b) with A1.
(a) (b)Figure 4
Performance index versusSNR with A2.


## 6. Conclusion

To conclude, we have proposed a new approach dedicated to blindly separating instantaneous mixtures of cyclostationary sources. It operates into three steps: first, we compute the cyclic autocorrelation matrices of observation signals, then, a new detection procedure is used to select specific matrices, and finally a nonunitary joint diagonalization algorithm based on BFGS method is employed to estimate the mixing matrix and recover the sources. One of the main advantages of such an approach is that it applies for source signals which may share one or more common cyclic frequencies. Extensions for further research would be to undertake the blind separation of convolutive mixtures of cyclostationary sources.



---

*Source: 2546838-2017-02-22.xml*