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.

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 of emitted source signals that are received by sensors (see Figure 1). The relationship of mixing system is given bywhere is the observations vector, is the vector of unknown source signals, is the additive noise vector, and is the unknown mixing matrix. The purpose of BSS is to find an estimate of and retrieve from only, aswhere , denotes pseudoinverse matrix, 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), , represents the phase vector (corresponding to phase shift ambiguity in complex domain of the source signals), and is square diagonal matrix containing the elements of the vector . Thus, one should look for a separating matrix such that

The following assumptions are held in this paper:(A1)The mixing matrix has full column rank .(A2)The source signals are zero-mean, cyclostationary, and mutually uncorrelated and may share one or more common cyclic frequencies.(A3)The noise signals 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 Figure 2).

A time process is first-order cyclostationary, if its first-order moment is periodic: is second-order cyclostationary, if its second-order moments are periodic:where is the mathematical expectation and the () stands for complex conjugation. In this case, the correlation function is decomposed into Fourier series as follows:withwhere is the cyclic autocorrelation function with as the time lag and 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 asWe 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 of real or complex square matrices sized , which are decomposed such thatwhere stands for the transpose conjugate operator. The matrices sized () are diagonal ones, is a full column rank matrix, and is its Moore-Penrose pseudoinverse. The NU-JD problem consists of estimating the matrix from only the matrix set . 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 stands for the Frobenius norm and 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 has the following decomposition:where is the conjugate transpose operator and (resp., ) is the cyclic correlation matrix of source signals (resp., noise signals). The expression of can be further developed. Knowing that the noise is stationary and using (7), we havewhereThis implies thatTherefore, if , thenPractically, the matrix has one of the following structures:(a)If the source signals have pairwise distinct cyclic frequencies, then is a diagonal matrix with only one nonnull element corresponding to the th source at (for this kind of structure, using , a detection procedure has been proposed in [20]).(b)If the source signals share one or more cyclic frequencies, then 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 at different time lags and cyclic frequencies (). The question now being asked is the following: how to compose the matrices set 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 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 to perform BSS. When the cyclic frequencies are unknown, we calculate 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 which correspond to the structure (b). In other words, the procedure has to detect matrices ( is whitening matrix such that ) which have the same number of eigenvalues as since all source signals are present at a given common cyclic frequency. One possible way to blindly access to diagonal terms of consists in computing the eigenvalues of . In fact, using the whitening matrix definition in [3], we havewhere the vector contains the eigenvalues of the square matrix , is an unitary matrix, and , for are the eigenvalues of . 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 . Ideally equals 1. However, in practice, the matrices can never be strictly diagonal. Thus, matrices should be selected as with close to zero. Furthermore, in order to detect the matrices where all source signals are present in main diagonal of and to avoid the low energy matrices , we add to :where denotes the matrix determinant and is a small positive constant. Finally, if a given satisfies and , 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 at different time lags for each cyclic frequency then to diagonalize simultaneously the set builtThe 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:where is the th matrix which belongs to the set 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 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 guesses and and a given number of iterations , the following instructions are iterated until converges to the solution.() Look for a search direction by solving () Perform a line search to find the optimal step-size (positive, a small enough number) in the previous direction (found in the first step).() Update() Set() Using the rank-one updates using gradient evaluations, the Hessian matrix is estimated as follows: and denote the inverse and transpose of a matrix, respectively, is the complex absolute gradient matrix of the cost function given in (20) which is defined as (see [21]) where stands for the complex conjugate of the complex matrix and is the partial derivative operator. 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 can be initialized with the identity matrix and 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 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 is a 4th-degree polynomial in whose expression iswhere its coefficients are given bywith denoting the trace operator and . 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.

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

4.4. Summary

The principle proposed is summarized in Algorithm 1.

Set building:
(i) Detect the useful matrices using the selection procedure;
(ii) Compute the set matrices for each cyclic frequency :
;
Initialization:
(i) Given an initial guess and an approximate Hessian ;
(ii) Given the number of iterations ;
(iii) Given a small positive threshold ;
for do
 Compute whose expression is given by Eq. (26);
 Perform a line search to find the optimal step-size ;
 Set ;
 Compute whose expression is given by Eq. (24);
if 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 aswhere are i.i.d and zero-mean random binary sequences, represent the symbol periods, they are, respectively, equal to and , are the normalized carrier frequencies, are the carrier phases, they, respectively, equal and , and is a triangular waveform such thatThe source signals share the cyclic frequencies , , and their multiples. The mixing matrices were randomly generated. We use in the first two simulations and in the last one:The signal-to-noise ratio (SNR) is computed as . We assume that the source cyclic frequencies are unknown and consequently we use the detection procedure described in (17) and (18) with and for all simulations to construct the set of matrices to be joint diagonalized which correspond to the source cyclic correlation matrices with a diagonal structure. In order to evaluate the quality of the estimation, one can measure the Moreau-Amari index presented in [22] defined aswhere is the element of . The closer to zero in linear scale ( in logarithmic scale), the higher the separation accuracy. Regarding the charts, is given in decibel and estimated by averaging 100 independent trials. The proposed method 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 [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.

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.

Appendix

A. Definitions and Properties

Let us consider three square matrices , , and and two rectangular matrices () and () and a square matrix . The following properties introduced in [21] are used in the coming developments:() .() .() .() .() .() .() .() .) ,where denotes Kronecker product, is the vectorization operator which when it is applied on a given square matrix , it concatenates its columns in a column vector such that denotes the “transformation” matrix defined aswhere is the () matrix whose all elements equal 1, is the () identity matrix, and are two matrix operators defined as

B. Coefficients of the 4th-Degree Polynomial

As we highlighted previously, we prefer to hide the dependency upon the iteration for simplicity. Using properties and , the cost function is expressed asFrom properties , , and , we find thatLet us introduce the four following matrices , , , and :

In the same way, we substitute for ; therefore we find that

As a result, (B.4) is nothing but a fourth-degree polynomial in .

Disclosure

This work reflects only the author’s view and the project partnership or the European Commission is not responsible for any use that may be made of the information it contains.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was prepared with the support of the Erasmus Mundus-Al Idrisi II project financed by the Erasmus Mundus Programme of the European Union.