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

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 themixingmatrix 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 stateof-the-art methods.


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 2 Mathematical Problems in Engineering 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.

Problem Statement
The BSS problem can be modelled in a simple linear instantaneous mixture of  emitted source signals that are received by  sensors (see Figure 1).The relationship of mixing system is given by where x() ∈ C  is the observations vector, s() ∈ C  is the vector of unknown source signals, b() ∈ C  is the additive noise vector, and A ∈ C × is the unknown mixing matrix.The purpose of BSS is to find an estimate Ã of A and retrieve s() from x() only, as where  2 = −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 , . . .,   ]  , ∀  ∈ 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 Ã# such that The following assumptions are held in this paper: (A1) The mixing matrix A has full column rank ( ≥ ).2).A time process {()} is first-order cyclostationary, if its first-order moment is  periodic:

Some Recalls
{()} is second-order cyclostationary, if its second-order moments are  periodic: 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: with where R   () is the cyclic autocorrelation function with  as the time lag and  = { = /,  ∈ 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 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.

Nonunitary Joint Diagonalization of Matrices.
We consider a set of  ( ≥ 2) real or complex square matrices sized ( × ), M = {Λ 1 , . . ., Λ  } which are decomposed such that where (⋅)  stands for the transpose conjugate operator.The matrices D  sized ( × ) are diagonal ones, A is a ( × ) 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]: where ‖ ⋅ ‖ 2  stands for the Frobenius norm and Off diag{⋅} denotes the zero-diagonal operator (see Appendix).

Suggested Method
4.1.Our Approach.From (1), it can be easily shown that the cyclic correlation matrix of observation signals x() has the following decomposition: where (⋅)  is the conjugate transpose operator and R  s () (resp., R  b ()) is the cyclic correlation matrix of source signals (resp., noise signals).The expression of R  b () can be further developed.Knowing that the noise is stationary and using (7), we have where This implies that Therefore, if  ̸ = 0, then Practically, the matrix R  s () has one of the following structures: (a) If the source signals have pairwise distinct cyclic frequencies, then R  s () is a diagonal matrix with only one nonnull element corresponding to the th source at   (for this kind of structure, using R  x (), a detection procedure has been proposed in [20]).
(b) If the source signals share one or more cyclic frequencies, then R  s () at these shared frequencies is a diagonal matrix with  nonnull elements which is our case of interest.Thus, we propose to diagonalize jointly the set of cyclic correlation matrices of observation signals R  x () at different time lags and cyclic frequencies ( ̸ = 0).The question now being asked is the following: how to compose the matrices set M jd 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 signals s() 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 R  s () to perform BSS.When the cyclic frequencies are unknown, we calculate R  x () 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 R  x () which correspond to the structure (b).In other words, the procedure has to detect matrices R x  () = WR  x ()W  (W is  ×  whitening matrix such that WAA  W  = I  ) which have the same number  of eigenvalues as R  s () since all source signals are present at a given common cyclic frequency.One possible way to blindly access to diagonal terms of R  s () consists in computing the eigenvalues of R x  ().In fact, using the whitening matrix definition in [3], we have where the vector eig(M) contains the eigenvalues of the square matrix M, WA is an unitary matrix, and   , for  = 1, . . .,  are the eigenvalues of R  s ().Therefore, we propose the following new criteria: where  is positive constant.One can note using the invariance property of the Frobenius norm under an unitary transformation that ‖R  x  ()‖ 2  = ‖R  s ()‖ 2  .Ideally C 1 equals 1.However, in practice, the matrices R  s () can never be strictly diagonal.Thus, matrices R x  () should be selected as C 1 ≥ 1 −  with  close to zero.Furthermore, in order to detect the matrices where all source signals are present in main diagonal of R  s () and to avoid the low energy matrices R  x  (), we add to C 1 : where det(⋅) denotes the matrix determinant and  is a small positive constant.Finally, if a given R  x  () satisfies C 1 and C 2 , 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 R  x () at different time lags for each cyclic frequency {  /  = 1, . . ., } then to diagonalize simultaneously the set built The joint diagonalization is ensured by the BFGS method which will be detailed in the following subsection.

BFGS Based NU-JD Algorithm.
We use the coming cost function to figure out the NU-JD problem: where R   x (  ) is the th matrix which belongs to the set M jd 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 ∈ C × 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.

Algorithm Principle.
From initial guesses B (0) and He (0) and a given number of iterations  iter , the following instructions are iterated until B () converges to the solution.
(S 1 ) Look for a search direction d (−1) by solving (S 2 ) Perform a line search to find the optimal step-size  (positive, a small enough number) in the previous direction (found in the first step).
(S 3 ) Update (S 4 ) Set s (−1) = d (−1) ; (S 5 ) Using the rank-one updates using gradient evaluations, the Hessian matrix is estimated as follows: (24) [⋅] −1 and [⋅]  denote the inverse and transpose of a matrix, respectively, ∇  F jd (B) is the complex absolute gradient matrix of the cost function given in (20) which is defined as (see [21]) where B * stands for the complex conjugate of the complex matrix B and /B * is the partial derivative operator.
∇  F jd (B) was calculated earlier in [19] and found to be equal to where  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).

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 F jd (B () ) with respect to .For simplicity, we prefer to hide the dependency upon the iteration .It is a matter of standard algebraic manipulations to show that F jd (B () ) is a 4thdegree polynomial in  whose expression is where its coefficients are given by   with tr[⋅] denoting the trace operator and Γ = He −1 ∇  F jd (B).
The optimal  can be found by polynomial rooting of the derivative third-order polynomial, namely, by solving with respect to , 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.

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 stepsize 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   required to attain the convergence.In practical applications, the computational time necessary to build the set of the  matrices should be counted too.

Summary.
The principle proposed is summarized in Algorithm 1.

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 Set building: (i) Detect the useful matrices R   () using the selection procedure; (ii) Compute the set matrices M jd for each cyclic frequency   : M jd = {R   x ( 1 ), . . ., R   x ( max )}; Initialization: (i) Given an initial guess B (0) and an approximate Hessian He (0) ; (ii) Given the number of iterations  iter ; (iii) Given a small positive threshold ; for  ← 1,  iter do Compute ∇  F jd (B) whose expression is given by Eq. ( 26); Perform a line search to find the optimal step-size ; Set B () ← B (−1) − He −1 (−1) ∇  F jd (B (−1) ); Compute He () whose expression is given by Eq. ( 24 where a =1,2 () are i.i.d and zero-mean random binary sequences,  =1,2 represent the symbol periods, they are, respectively, equal to  1 = 10 and  2 = 8,  =1,2 = 0.2 are the normalized carrier frequencies,  =1,2 are the carrier phases, they, respectively, equal  1 = /6 and  2 = /8, and g() is a triangular waveform such that The source signals share the cyclic frequencies  1 =  2 = 0.2, 5 1 = 10/ 1 = 8/ 2 = 1, and their multiples.The mixing matrices were randomly generated.We use A 1 in the first two simulations and A 2 in the last one: ( The signal-to-noise ratio (SNR) is computed as SNR = −10 log 10 ( 2  ).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 R  s () 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 where  , is the (, )th element of Z = BA.The closer to zero in linear scale (−∞ in logarithmic scale), the higher the separation accuracy.Regarding the charts,  perf (⋅) is given in decibel and estimated by averaging 100 independent trials.The proposed method JD BFGS 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 JD GRAD [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 ( >  or  = ).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.

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.

B. Coefficients of the 4th-Degree Polynomial
As we highlighted previously, we prefer to hide the dependency upon the iteration  for simplicity.Using properties P 2 and P 6 , the cost function F jd (B) is expressed as As a result, (B.4) is nothing but a fourth-degree polynomial in .

Figure 1 :
Figure 1: The relationship between sources and observations.

(
A2) The source signals are zero-mean, cyclostationary, and mutually uncorrelated and may share one or more common cyclic frequencies.(A3) The noise signals b() are stationary, zero-mean random signals, and mutually uncorrelated with the source signals.

Figure 3 :Figure 4 :
Figure 3: Performance index versus SNR (a) and versus number of selected matrices (b) with A 1 .

Table 1 :
The gradient matrix computational cost.

Table 2 :
The computational cost of the optimal step-size.

Table 3 :
The BFGS algorithm computational cost.