US20040247054A1

US20040247054A1 - Method and device for the fourth-order, blind identification of an under-determined mixture of sources

Info

Publication number: US20040247054A1
Application number: US10/814,808
Authority: US
Inventors: Anne Ferreol; Laurent Albera; Pascal Chevalier
Original assignee: Thales SA
Current assignee: Thales SA
Priority date: 2003-04-01
Filing date: 2004-04-01
Publication date: 2004-12-09
Also published as: EP1465338A3; EP1465338A2; FR2853480B1; ATE408869T1; DE602004016587D1; DK1465338T3; EP1465338B1; US7336734B2; FR2853480A1

Abstract

A method for the fourth-order, blind identification of at least two sources in a system comprising a number of sources P and a number N of reception sensors receiving the observations, said sources having different tri-spectra. The method comprises at least the following steps: a step for the fourth-order whitening of the observations received on the reception sensors in order to orthonormalize the direction vectors of the sources in the matrices of quadricovariance of the observations used; a step for the joint diagonalizing of several whitened matrices of quadricovariance in order to identify the spatial signatures of the sources. Application to a communication network.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention relates especially to a method for the learned or blind identification of a number of sources P that is potentially greater than or equal to the number N of sensors of the reception antenna.

It can be used for example in the context of narrow-band multiple transmission.

It is used for example in a communications network.

It can be applied especially in the field of radio communications, space telecommunications or passive listening to these links in frequencies ranging for example from VLF to EHF.

FIG. 1 is a schematic drawing exemplifying an array of several reception sensors or receivers, each sensor receiving signals from one or more radio communications transmitters from different directions of arrival

Each sensor receives signals from a source with a phase and amplitude that are dependent on the angle of incidence of the source and the position of the sensor. FIG. 2 is a drawing exemplifying the parametrization of the direction of a source. This direction is parametrized by two angles corresponding to the azimuth angle θ and the elevation angle Δ.

2. Description of the Prior Art

The past 15 years or so have seen the development of many techniques for the blind identification of signatures or source direction vectors, assumed to be statistically independent. These techniques have been developed in assuming a number of sources P smaller than or equal to the number of sensors N. These techniques have been described in the references [1][3][7] cited at the end of the description. However, for many practical applications such as HF radio communications, the number of sources from which signals are received by the sensors is increasing especially with the bandwidth of the receivers, and the number of sources P can therefore be greater than the number of sensors N. The mixtures associated with the sources are then said to be under-determined.

A certain number of methods for the blind identification of under-determined mixtures of narrow-band sources for networks have been developed very recently and are described in the references [2] [7-8] and [10]. The methods proposed in the references [2] and [7-8] make use of information contained in the fourth-order (FO) statistics of the signals received at the sensors while the method proposed in the reference [10] make use of the information contained in one of the characteristic functions of the signals received. However, these methods have severe limitations in terms of the prospects of their operational implementation. Indeed, the method described in the reference [2] is very difficult to implement and does not provide for the identification of the sources having the same kurtosis values (standardized fourth-order cumulant). The methods described in the references [7-8] assume that the sources are non-circular. These methods give unreliable results in practice. Finally, the reference method [10] has been developed solely for mixtures of sources with real (non-complex) values.

The object of the present invention relates especially to a new method for the blind identification of an under-determined mixture of narrow-band sources for communications networks. The method can be used especial y to identify up to N ²−N+1 sources from N identical sensors and up to N²sources with N different sensors, in assuming only that the sources have different tri-spectra and non-zero, same-sign kurtosis values (this hypothesis is practically always verified in the context of radio communications).

SUMMARY OF THE INVENTION

The invention relates to a method for the fourth-order, blind identification of at least two sources in a system comprising a number of sources P and a number N of reception sensors receiving the observations, said sources having different tri-spectra, wherein the method comprises at least the following steps:

a step for the fourth-order whitening of the observations received on the reception sensors in order to orthonormalize the direction vectors of the sources in the matrices of quadricovariance of the observations used,

a step for the joint diagonalizing of several whitened matrices of quadricovariance in order to identify the spatial signatures of the sources.

The number of sources P is for example greater than the number of sensors N.

The method can be used in a communication network.

The method according to the invention has especially the following advantages:

it enables the identification of a number of sources greater than the number of sensors:

for identical sensors: N ²−N+1 sources having different tri-spectra and non-zero and same-sign kurtosis values,

for different sensors (arrays with polarization diversity and/or pattern diversity and/or coupling etc) N ²sources having different tri-spectra and non-zero, same-sign kurtosis values,

the method is robust with respect to Gaussian noise, even spatially correlated Gaussian noise,

it enables the goniometry of each source identified, using a wavefront model attached to the signature, with a resolution potentially higher than that of existing methods,

it enables the identification of I (N ²−N+1) cyclostationary sources if the sensors are identical and I×N²cyclostationary sources if the sensors are different: with polarization diversity and/or pattern diversity and/or coupling, where I is the number of cyclical frequencies processed,

using a performance criterion, it enables the quantitative evaluation of the quality of estimation of the direction vector of each source and a quantitative comparison of two methods for the identification of a given source,

using a step for the selection of the cyclical frequencies, it enables the processing of a number of sources greater than the number of sources processed by the basic method.

BRIEF DESCRIPTION OF THE DRAWINGS

Other features and advantages of the invention shall appear more clearly from the following description along with the appended figures, of which: [0026]
FIG. 1 shows an example of a communication network, [0027]
FIG. 2 shows parameters of a source, [0028]
FIG. 3 is a functional diagram of the method according to the invention, [0029]
FIG. 4 shows an example of spatial filtering, [0030]
FIGS. 5 and 6 show examples of variations of the performance criterion as a function of the number of samples observed, comparing the performance of the method with two prior art methods. [0031]
FIGS. [0032] 7 to 9 show three alternative embodiments of the method described in FIG. 3 implementing a selection of the cyclical frequencies.

MORE DETAILED DESCRIPTION

For a clear understanding of the object of the invention, the following example is given by way of an illustration that in no way restricts the scope of the invention for a radio communications network in a multiple-transmission, narrow-band context, with sources having different tri-spectra (of cumulants). [0033]
Each sensor of the network, formed by N receivers, receives a mixture of P narrow-band (NB) sources which are assumed to be statistically independent. On this assumption, the vector of the complex envelopes of the signals at output of the sensors is written as follows: [0034] $\begin{matrix} x (t) = \sum_{p = 1}^{P} s_{p} (t) a_{p} + b (t) = A s (t) + b (t) & (1) \end{matrix}$
where s[0035] _p(t) is the signal of the p^thsource as well as the p^thcomponent of the vector s(t), b(t) is the assumed Gaussian noise vector with any covariance, α_pis the signature or the direction vector of the p^thsource and A (N×P) is the matrix of the vectors α_p(direction vectors of the sources).
It is an object of the invention especially to identify the direction vectors α[0036] _pof each of the sources when, especially, the number of sources P is potentially greater than the number of sensors N.
From this identification, it is then possible to apply techniques for the extraction of the sources by the spatial filtering of the observations. The blind extraction is aimed especially at restoring the information signals conveyed by the sources in not exploiting any a priori information (during normal operation) on these sources. [0037]
Fourth-Order Statistics [0038]
The method according to the invention makes use of the fourth-order statistics of the observations corresponding to the time-domain averages Q[0039] _x(τ₁,τ₂,τ₃)=<Q_x(τ₁,τ₂,τ₃)(t)>, on an infinite horizon of observation, of certain matrices of quadricovariance, Q_x(τ₁,τ₂,τ₃)(t), , sized (N²×N²). The elements, Q_x(τ₁,τ₂,τ₃)[i, j, k, l](t) of these matrices are, for example, defined by the relationship:
Q _x(τ₁,τ₂,τ₃)[i,j,k,l][t]=Cum(x _i(t),x _j(t−τ ₁)*,x _k(t−τ ₂)*, x _l(t−τ ₃)) (2)
where * is the conjugate complex symbol, x[0040] _i(t) is the i^thcomponent of the vector x(t), <.> is the operation of time-domain averaging on an infinite horizon of observation and (τ₁,τ₂,τ₃) is a triplet of delays. Assuming that Q_x(τ₁,τ₂,τ₃)[i, j, k, l] is the element [N(i−1)+j, N(k−1)+l] of the matrix Q_x(τ₁,τ₂,τ₃), assuming that the noise is Gaussian and using the expression (1) in the expression (2), the matrix of quadricovariance Q_x(τ₁,τ₂,τ₃) is written as follows:
Q _x(τ₁,τ₂,τ₃)=(A{circle over (×)}A*)Q _s(τ₁,τ₂,τ₃)(A{circle over (×)}A*)^H (3)
where Q[0041] _s(τ₁,τ₂,τ₃) is the averaged matrix of quadricovariance of s(t) with a dimension (P²×P2), A=[α₁. . . α_P], {circle over (×)} is the Kronecker product and ^Hdesignates the transpose and conjugate. Assuming statistically independent sources, the matrix Q_s(τ₁,τ₂,τ₃) is formed by at least P⁴−P zeros and the expression (3) is simplified as follows: $\begin{matrix} Q_{x} (τ_{1}, τ_{2}, τ_{3}) = \sum_{p = 1}^{P} c_{p} (τ_{1}, τ_{2}, τ_{3}) (a_{p} \otimes a_{p}^{*}) {(a_{p} \otimes a_{p}^{*})}^{H} & (4a) \\ = A_{Q} C_{s} (τ_{1}, τ_{2}, τ_{3}) A_{Q}^{H} & (4b) \end{matrix}$
where A[0042] _Qis a matrix sized (N²×P) defined by A_Q=[(α₁{circle over (×)}α₁*), . . . , (α_p{circle over (×)}α_p*)], C_s(τ₁,τ₂,τ₃) is a diagonal matrix sized (P×P) defined by C_s(τ₁,τ₂,τ₃)=diag[c₁(τ₁,τ₂,τ₃), . . . , c_p(τ₁,τ₂,τ₃)] and where c_p(τ₁,τ₂,τ₃) is defined by:
c _p(τ₁,τ₂,τ₃)=<Cum(s _p(t), s _p(t−τ ₁)*, s _p(t−τ ₂)*, s _p(t−τ ₃))> (5)
The expression (4b) has an algebraic structure similar to that of the correlation matrix of the observations used in the algorithm SOBI described in the reference [1]. The notation used here below will be Q[0043] _x=Q_x(0, 0, 0), c_p=c_p(0, 0, 0), C_s=C_s(0, 0, 0) in order to deduce the relationship (4b) therefrom:
Q_x=A_QC_sA_Q ^H (6)
It is assumed here below that the number of sources P is such that P≦N[0044] ², that the matrix A_Qis of full rank, that the averaged cumulants c_p, 1≦p≦P, are non-zero (non-Gaussian sources) and same-sign cumulants and that, for any pair (i, j) of sources, there is at least one triplet of delays (τ₁,τ₂,τ₃) such that |τ₁|+|τ₂|+|τ₃|=0 and
c _i(τ₁,τ₂,τ₃)/|c _i |≠c _j(τ₁,τ₂,τ₃)/|c _j| (7)
Fourth-Order Whitening Step [0045]
The first step of the method according to the invention, called FOBIUM, consists of the orthonormalization, in the matrix of quadricovariance Q[0046] _xof the expression (6), of the columns of the matrix A_Q, considered to be virtual direction vectors of the sources for the array of sensors considered. To this end, the method considers the eigen-element decomposition of the P rank hermitian matrix Q_xgiven by
Q_x=E_xΛ_xE_x ^H (8)
where A[0047] _xis the real longitudinal diagonal, with a dimension (P×P), of the P non-zero eigenvalues of Q_x, and E_xis the matrix sized (N²×P) of the associated, orthonormal eigenvectors. For a full-rank matrix A_Q, it can be shown that there is equivalence between assuming that the kurtosis values of the sources have a same sign ε (ε=±1) and assuming that the eigenvalues of Λ_xalso have a same sign ε. In this context, it is possible to build the following whitening matrix T sized (P×N²):
T=(Λ _x)^−1/2 E _x ^H (9)
The whitening matrix sized (P×N[0048] ²) is defined from the square root of the real diagonal matrix sized (P×P) of the P non-zero eigenvalues of the matrix of quadricovariance and of the transpose of the matrix of the associated eigenvectors with a dimension (P×N²) where (Λ_x)^−1/2is the inverse of the square root of Λ_x. From the expressions (6) and (8), it is deduced that:
εTQ _x T ^H =TA _Q(εC _s)A _Q ^H T ^H =I _P (10)
where I[0049] _Pis the identity matrix with a dimension (P×P) and where εC_s=diag[|c₁|, . . . , |c_p|]. This expression shows that the matrix TA_Q(εC_s)^1/2with a dimension (P×P) is a unitary matrix U. It is deduced from this that:
TA _Q =U(εC _s)^−1/2 (11)
Fourth-Order Identification Step [0050]
From the expressions (4b) and (11), it is deduced that: [0051]
W(τ₁,τ₂,τ₃)=TQ _x(τ₁,τ₂,τ₃)T ^H =U(εC _x)^−1/2 C _s(τ₁,τ₂,τ₃)(εC _s)^−1/2 U ^H (12)
where W(τ[0052] ₁,τ₂,τ₃) is the matrix of quadricovariance whitened at the fourth order by the matrix Q_x. This expression shows that the unitary matrix U diagonalizes the matrices T Q_x(τ₁,τ₂,τ₃) T^Hand that the associated eigenvalues are the diagonal terms of the diagonal matrix (εC_s)^−1/2C_s(τ₁,τ₂,τ₃) (εC_s)^−1/2. For a given triplet of delays (τ₁,τ₂,τ₃), the matrix U is unique, give or take one permutation and one unit diagonal matrix, when the elements of the matrix (εC_s)^−1/2C_s(τ₁,τ₂,τ₃) (εC_s)^−1/2) all are different. If not, the method uses a set of K triplets (τ₁ ^k,τ₂ ^k,τ₃ ^k), 1≦k≦K, defined as follows: for all pairs of sources (i, j), there is at least one triplet (τ₁ ^k,τ₂ ^k,τ₃ ^k), such that the condition of the equation (7) is verified. In these conditions, the unitary matrix U is the only matrix U_solwhich, to the nearest permutation and to the nearest unit diagonal matrix, jointly diagonalizes the K matrices T×Q_x(τ₁ ^k,τ₂ ^k,τ₃ ^k)×T^H. Consequently, the matrix U_sol, which resolves the above problem, is written as a function of the unitary matrix U as follows:
U_sol=UΛΠ (13)
where Λ and Π are respectively the unit diagonal matrix and the permutation matrix referred to here above. From the equations (11) and (13), it is possible to deduce the matrix A[0053] _Qto the nearest unit diagonal matrix and to the nearest permutation, which is expressed by:
T ^# U _sol =[b ₁ . . . b _P ]=E _xΛ_x ^1/2 U _sol =A _Q(εC _s)^1/2Λ (14)
where T[0054] ^# is the pseudo-inverse of the matrix T. Each column, b_l(1≦l≦P), of the matrix T^# U_solcorresponds to one of the vectors μ_q|c_q|^1/2(α_q{circle over (×)}α_q*), 1≦q≦P, where μ_qis a complex scalar value such that |μ_q|=1. Consequently, in converting each column b_lof the matrix T^# U_solinto a matrix B_lwith a dimension (N×N such that B_l[i, j]=b_l((i−1)N+j) (1≦i, j≦N), it is deduced therefrom that:
B _l=μ_q |c _q|^1/2α_qα_q ^Hpour (1≦l, q≦P) (15)
The matrix B[0055] _lis built from the vector b_land depends on a complex scalar value, the square root of the cumulant and the direction vector of the q^thsource and of its conjugate.
In this context, the direction vector α[0056] _qof the q^thsource is associated with the eigenvector of B_lassociated with the highest eigenvalue.

SUMMARY OF THE PRINCIPLE OF THE INVENTION

In brief, the different steps of the method according to the invention include at least the following steps: for L vector observations received in the course of the time: x(lTe) (1≦l ≦L), where T[0057] _eis the sampling period.
Estimation [0058]
Step 1: The estimation, through Q;{circumflex over ( )}[0059] _x, of the matrix of quadricovariance Q_x, from the L observations x(lTe), using a non-skewed and asymptotically consistent estimator. Depending on the nature of the sources, the estimator is adapted as follows:
Stationary and centered case: empirical estimator used in the reference [3]. [0060]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0061]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0062]
Whitening [0063]
Step 2: The eigen-element decomposition of the estimated matrix of quadricovariance Q;{circumflex over ( )}[0064] _x, estimating the number of sources P and restricting this eigenvalue decomposition to the P main components: Q;{circumflex over ( )}_x≈E;{circumflex over ( )}_xΛ;{circumflex over ( )}_xE;{circumflex over ( )}_x ^H, where A;{circumflex over ( )}_xis the diagonal matrix containing the P highest modulus eigenvalues and E;{circumflex over ( )}_xis the matrix containing the associated eigenvectors.
Step 3: The building of the whitening matrix: T;{circumflex over ( )}=(Λ;{circumflex over ( )}[0065] _x)^−1/2E;{circumflex over ( )}_x ^H.
Selection of the Triplets [0066]
Step 4: The selection of K triplets of delays (τ[0067] ₁ ^k,τ₂ ^k,τ₃ ^k) where |τ₁ ^k|+|τ₂ ^k|+|τ₃ ^k|≠0.
Estimation [0068]
Step 5: The estimation, through Q;{circumflex over ( )}[0069] _x(τ₁ ^k,τ₂ ^k,τ₃ ^k), of K matrices of quadricovariance Q_x(τ₁ ^k,τ₂ ^k,τ₃ ^k). As in the step 1, this estimation depends especially on the assumptions made on the observations:
Stationary and centered case: empirical estimator used in the reference [3]. [0070]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0071]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0072]
Identification [0073]
Step 6: The computation of the matrices T;{circumflex over ( )} Q;{circumflex over ( )}[0074] _x(τ₁ ^k,τ₂ ^k,τ₃ ^k) T;{circumflex over ( )}^Hand the estimation, by U;{circumflex over ( )}_sol, of the unitary matrix U_solthrough the joint diagonalizing of the K matrices T;{circumflex over ( )} Q;{circumflex over ( )}_x(τ₁ ^k,τ₂ ^k,τ₃ ^k) T;{circumflex over ( )}^H.
Step 7: The computation of T{circumflex over ( )} [0075] ^#U;{circumflex over ( )}_sol=[b;{circumflex over ( )}_l. . . b;{circumflex over ( )}_P] and the building of the matrices B;{circumflex over ( )}_lsized (N×N).
Step 8: The estimation, by α;{circumflex over ( )}[0076] _P, of the signatures α_q(1≦q≦P) of the P sources in carrying out a decomposition into elements on each matrix B;{circumflex over ( )}_l.
Applications [0077]
At the end of the step 8, the method has identified the direction vectors of P non-Gaussian sources having different tri-spectra with same-sign kurtosis values. p<N[0078] ²and P may reach N²−N+1 or N²depending on the type of sensors used.
Using this information, the method may implement a method of goniometry or a spatial filtering of antennas. [0079]
A method of goniometry can be used to determine the direction of arrival of the sources and more precisely the azimuth angle θ[0080] _mfor 1D goniometry and azimuth and elevation angles (θ_m, Δ_m) for 2D goniometry.
FIG. 4 represents a spatial filtering of antennas for spatial filtering structures. It enables especially the optimizing of reception from one or all the sources present by the spatial filtering of the observations. When several sources are of interest to the receiver, we speak of source separation techniques. When no a priori information on the sources is exploited, we speak of blind techniques. [0081]
Verification of the Quality of the Estimates [0082]
According to one alternative embodiment, the method comprises a step of qualitative evaluation, for each source, of the quality of identification of the associated direction vector. [0083]
This new criterion enables the intrinsic comparison of two methods of identification for the restitution of the signature of a particular source. This criterion, for the identification problem, is an extension of the one proposed in [5] for extraction. It is defined by the P-uplet [0084]
D(A, Â)=(α₁, α₂, . . . , α_P) (16)
where [0085] $\begin{matrix} α_{p} = \min_{1 \leq i \leq P} [d (a_{p}, {\hat{a}}_{i})] & (17) \end{matrix}$
and where d(u,v) is the pseudo-distance between the vectors u and v, such that: [0086] $\begin{matrix} d (u, v) = 1 - \frac{{\langle u^{H} v \rangle}^{2}}{(u^{H} u) (v^{H} v)} & (18) \end{matrix}$
In the simulations of FIGS. 5 and 6, there are P=6 statistically independent sources received on a circular array of N=3 sensors having a radius r such that r/λ=0.55 (λ: wavelength). The six sources are non-filtered QPSK sources having a signal-to-noise ratio of 20 dB with a symbol period T=4T[0087] _e, where T_eis the sampling period.
The incidence values of the sources are such that θ[0088] ₁=2.16°, θ₂=25.2°, θ₃=50°, θ₄=272.16°, θ₅=315.36°, θ₆=336.96° and the associated carrier frequencies verify Δf₁T_e=0, Δf₂T_e=1/2, Δf₃T_e=1/3, Δf₄T_e=1/5, Δf₅T_e=1/7 and Δf₆T_e=1/11. The JADE [3], SOBI [1] and FOBIUM method according to the invention are applied and the performance values α_qfor q=1 . . . 6 are evaluated after an averaging operation on 1000 results. For the FOBIUM method, we choose K=4 triplets of delays (τ₁ ^k,τ₂ ^k,τ₃ ^k) where τ₁ ^k=kT_eand τ₂ ^k=τ₃ ^k=0.
In the above assumptions, FIG. 5 shows the variation in α[0089] ₂(performance of the second source) resulting from the JADE (b), SOBI (c) and FOBIUM (a) separators as a function of the number L of samples. The curves show firstly that the JADE and SOBI methods present difficulties in identifying the direction vector of the second source in an under-determined mixture context and that, secondly, that the FOBIUM method performs very well.
FIG. 6 gives a view, in the same context, of the variations of all the α[0090] _p(1≦p≦6) values resulting from the FOBIUM method as a function of L. The curve (index p) is associated with the p^thsource. It can be seen that all the coefficients α_pconverge towards zero and that, asymptomatically, the direction vectors are perfectly identified.
Variants of the Cyclical FOBIUM Method [0091]
FIGS. 7 and 8 show two examples of variants of the method according to the invention known as the cyclical FOBIUM method. [0092]
The idea lies especially in introducing selectivity by the cyclical frequencies into the method presented here above and is aimed especially at the blind identification, with greater processing capacity, of under-determined mixtures of cyclostationary sources. [0093]
The major difference between the [0094] steps 1 to 8 explained here above and this variant is the implementation of a step for the cyclical isolation of the sources by fourth-order discrimination according to their cyclical frequencies. It is thus possible to separately identify the sources associated with a same fourth-order cyclical parameter without being disturbed by the other sources processed separately.
The two variants shown in FIGS. 7 and 8 can be implemented by reiterating the process of cyclical isolation on the “other sources” with other cyclical parameters. The process of cyclical isolation can be applied several times in a third version illustrated in FIG. 9. [0095]
Fourth-Order Cyclical Statistics [0096]
The fourth-order cyclical statistics of the observations or sensor signals used are characterized by the matrices of cyclical quadricovariance Q[0097] _x ^ε(α,τ₁,τ₂,τ₃) with a dimension (N²×N²) where the elements Q_x ^ε(τ₁,τ₁,τ₁)[i, j, k, l], are defined by: $\begin{matrix} Q_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) [i, j, k, l] = < Cum (x_{i} (t), {x_{j} (1 - τ_{1})}^{*}, {x_{k} (t - τ_{2})}^{*}, x_{l} (t - τ_{3})) \exp (- j2 πα t) > Q_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) [i, j, k, l] = < Cum (x_{i} (t), {x_{j} (1 - τ_{1})}^{*}, x_{k} (t - τ_{2}), x_{l} (t - τ_{3})) \exp (- j2 πα t) > Q_{x}^{3} (α, τ_{1}, τ_{2}, τ_{3}) [i, j, k, l] = < Cum (x_{i} (t), x_{j} (1 - τ_{1}), x_{k} (t - τ_{2}), x_{l} (t - τ_{3})) \exp (- j2 πα t) > & (19) \end{matrix}$
It can be seen that Q[0098] _x ^ε(α,τ₁,τ₂,τ₃) is associated with the ε^thfourth-order moment. In stating that Q_x ^ε(α,τ₁,τ₂,τ₃)[i, j, k, l] is the element [N(i−1)+j, N(k−1)+l] of the matrix Q_x ^ε(α,τ₁,τ₂,τ₃) and assuming that the noise is Gaussian, the matrix Q_x ^ε(α,τ₁,τ₂,τ₃), is written as follows, in using (1) and (19): $\begin{matrix} Q_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) = (A \otimes A^{*}) Q_{s}^{1} (α, τ_{1}, τ_{2}, τ_{3}) {(A \otimes A^{*})}^{H} Q_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) = (A \otimes A^{*}) Q_{s}^{2} (α, τ_{1}, τ_{2}, τ_{3}) {(A \otimes A)}^{T} Q_{x}^{3} (α, τ_{1}, τ_{2}, τ_{3}) = (A \otimes A) Q_{s}^{3} (α, τ_{1}, τ_{2}, τ_{3}) {(A \otimes A)}^{T} & (20) \end{matrix}$
Where Q[0099] _s ^ε(α,τ₁,τ₂,τ₃) is a cyclical matrix of quadricovariance s(t) with a dimension (P²×P²), {circle over (×)} is the Kronecker product and ^Tdesignates the transpose. On the assumption of statistically independent sources, the matrix Q_s ^ε(α,τ₁,τ₂,τ₃) is formed by at least P⁴−P zeros and the expression (20) is simplified as follows: $\begin{matrix} Q_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) = \sum_{p = 1}^{P} c_{p}^{1} (α, τ_{1}, τ_{2}, τ_{3}) (a_{p} \otimes a_{p}^{*}) {(a_{p} \otimes a_{p}^{*})}^{H} = A_{Q} C_{s}^{1} (α, τ_{1}, τ_{2}, τ_{3}) A_{Q}^{H} Q_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) = \sum_{p = 1}^{P} c_{p}^{2} (α, τ_{1}, τ_{2}, τ_{3}) (a_{p} \otimes a_{p}^{*}) {(a_{p} \otimes a_{p})}^{T} = A_{Q} C_{s}^{2} (α, τ_{1}, τ_{2}, τ_{3}) B_{Q}^{T} Q_{x}^{3} (α, τ_{1}, τ_{2}, τ_{3}) = \sum_{p = 1}^{P} c_{p}^{3} (α, τ_{1}, τ_{2}, τ_{3}) (a_{p} \otimes a_{p}) {(a_{p} \otimes a_{p})}^{T} = B_{Q} C_{s}^{3} (α, τ_{1}, τ_{2}, τ_{3}) B_{Q}^{T} & (21) \end{matrix}$
where A[0100] _Qand B_Qare matrices with a dimension of (N²×P) defined by A_Q=[(α₁{circle over (×)}α₁*), . . . , (α_p{circle over (×)}α_p*)] and B_Q=[(α₁{circle over (×)}α₁*), . . . , (α_p{circle over (×)}α_p*)], C_s ^ε(α, τ₁,τ₂,τ₃) is a diagonal matrix sized (P×P) defined by C_s ^ε(α, τ₁,τ₂,τ₃)=diag[c₁ ^ε(α, τ₁,τ₂,τ₃), . . . , $\begin{matrix} c_{p}^{1} (α, τ_{1}, τ_{2}, τ_{3}) = < Cum (s_{p} (t), {s_{p} (t - τ_{1})}^{*}, {s_{p} (t - τ_{2})}^{*}, s_{p} (t - τ_{3})) \exp (- j2 π α t) > c_{p}^{2} (α, τ_{1}, τ_{2}, τ_{3}) = < Cum (s_{p} (t), {s_{p} (t - τ_{1})}^{*}, s_{p} (t - τ_{2}), s_{p} (t - τ_{3})) \exp (- j2 π α t) > c_{p}^{3} (α, τ_{1}, τ_{2}, τ_{3}) = < Cum (s_{p} (t), s_{p} (t - τ_{1}), s_{p} (t - τ_{2}), s_{p} (t - τ_{3})) \exp (- j2 π α t) > & (22) \end{matrix}$
It can be seen that the classic quadricovariance of (6) also verifies that Q[0101] _x=Q_x ¹(0, 0, 0, 0), c_p=c_p ¹(0, 0, 0, 0), C_s=C_s ¹(0, 0, 0, 0). Recalling that Q_x[i, j, k, l] is the element [N(i−1)+j, N(k−1)+l] of the matrix Q_x, the following is deduced:
Q_x=A_QC_sA_Q ^H (23)
In stating that Q[0102] _x[i, j, k, l] is the element [N(i−1)+l, N(k−1)+j] of {tilde over (Q)}_xwe obtain the matrix {tilde over (Q)}_xwhich is written as follows:
{tilde over (Q)}_x=B_QC_sB_Q ^H (24)
Whitening Step [0103]
The first step of the cyclical method orthonormalizes the columns of the matrices A[0104] _Qor B_Qcontained in the matrices Q_xor {tilde over (Q)}_xof the expressions (23)(24). The matrices Q_xand {tilde over (Q)}_xare P rank hermitian matrices and verify the following eigen-element decomposition:
Q_x=E_xΛ_xE_x ^Het {tilde over (Q)}_x={tilde over (E)}_x{tilde over (Λ)}_x{tilde over (Q)}_x ^H (25)
Where {tilde over (Λ)}[0105] _xis the diagonal matrix sized (P×P) of the P non-zero values of {tilde over (Q)}_xand {tilde over (E)}_xis the matrix sized (N²×P) of the associated eigenvectors. For a full-rank matrix B_Q, it can be shown that there is an equivalence assuming that that the kurtosis values of the sources have a same sign ε(ε=±1) and assuming that the eigenvalues of {tilde over (Λ)}_xalso have a same sign ε. In this context, the whitening matrix can be built according to {tilde over (T)} with a dimension of (P×N²):
{tilde over (T)}=({tilde over (Λ)}_x)^−1/2 {tilde over (E)} _x ^H (26)
where ({tilde over (Λ)}[0106] _x)^−1/2is the inverse of the square root of {tilde over (Λ)}_x. From the expressions (24) and (25), it is deduced that:
ε{tilde over (TQ)} _x {tilde over (T)} ^H ={tilde over (T)}B _Q(εC _s)B _Q ^H {tilde over (T)} ^H =I _P (27)
This expression shows that the matrix {tilde over (T)}B[0107] _Q(εC_s)^1/2with a dimension (P×P) is a unitary matrix Ũ. It is then deduced from this that:
{tilde over (T)}B _Q =Ũ(εC _s)^−1/2 (28)
It is recalled that the whitening matrix T of Q[0108] _xverifies:
TA _Q =U(εC _s)^−1/2 (29)
Step of Cyclical Isolation [0109]
From the expressions (28)(29) and (21), it is deduced that: [0110] $\begin{matrix} W_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) = T Q_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) T^{H} = {U (ɛ C_{s})}^{- 1 / 2} C_{s}^{1} (α, τ_{1}, τ_{2}, τ_{3}) {(ɛ C_{s})}^{- 1 / 2} U^{H} W_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) = T Q_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) {\tilde{T}}^{T} = {U (ɛ C_{s})}^{- 1 / 2} C_{s}^{2} (α, τ_{1}, τ_{2}, τ_{3}) {(ɛ C_{s})}^{- 1 / 2} {\tilde{U}}^{T} W_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) = \tilde{T} Q_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) {\tilde{T}}^{T} = {\tilde{U} (ɛ C_{s})}^{- 1 / 2} C_{s}^{3} (α, τ_{1}, τ_{2}, τ_{3}) {(ɛ C_{s})}^{- 1 / 2} {\tilde{U}}^{T} & (30) \end{matrix}$
When there are P[0111] ₁sources verifying c_i ^ε(α,τ₁,τ₂,τ₃)≠0 (1≦i≦P₁), the matrix W_x ^ε(α,τ₁,τ₂,τ₃) is a P₁≦P ranking matrix. Thus, the unitary matrices U and Ũ with a dimension P×P may be decomposed into two sub-matrices with a dimension P×P₁and P×(P−P₁) such that:
U=[U ₁ U ₂] and Ũ=[Ũ ₁ Ũ ₂] (31)
where the matrices U[0112] ₁and Ũ₁are sized P×P₁and U₂and Ũ₂are sized P×(P−P₁). The matrices U₁and Ũ₁contain the singular vectors associated with the non-zero singular values of W_x ^ε(α,τ₁,τ₂,τ₃). It is deduced from this that: $\begin{matrix} W_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) = T Q_{x}^{1} (α, τ_{1}, τ_{2}, τ_{3}) T^{H} = {U_{1} (ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{C}}_{s}^{1} (α, τ_{1}, τ_{2}, τ_{3}) {(ɛ {\tilde{C}}_{s})}^{- 1 / 2} U_{1}^{H} W_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) = T Q_{x}^{2} (α, τ_{1}, τ_{2}, τ_{3}) {\tilde{T}}^{T} = {U_{1} (ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{C}}_{s}^{2} (α, τ_{1}, τ_{2}, τ_{3}) {(ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{U}}_{1}^{T} W_{x}^{3} (α, τ_{1}, τ_{2}, τ_{3}) = \tilde{T} Q_{x}^{3} (α, τ_{1}, τ_{2}, τ_{3}) {\tilde{T}}^{T} = {{\tilde{U}}_{1} (ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{C}}_{s}^{3} (α, τ_{1}, τ_{2}, τ_{3}) {(ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{U}}_{1}^{T} & (32) \end{matrix}$
Where the matrix {tilde over (C)}[0113] _s ^ε(α,τ₁,τ₂,τ₃) is a diagonal matrix sized P₁×P₁formed by non-zero, diagonal elements c_i ^ε(α,τ₁,τ₂,τ₃) of the matrix C_s ^ε(α,τ₁,τ₂,τ₃). The matrix {tilde over (C)}_s={tilde over (C)}_s ¹(0,0,0,0) sized P₁×P₁is formed by corresponding elements c_i(1≦i≦P₁). Thus, after a singular value decomposition of W_x ^ε(α,τ₁,τ₂,τ₃), it is possible to determine the matrices T₁and {tilde over (T)}₁from singular values associated with the non-zero singular values and T₂and {tilde over (T)}₂from singular vectors associated with the zero singular values such that:
T₁=U₁Π₁ ^H, T₂=U₂Π₂ ^H, {tilde over (T)}₁=Ũ₁{tilde over (Π)}₁ ^Het {tilde over (T)}₂=Ũ₂{tilde over (Π)}₂ ^H (33)
where the matrices are Π[0114] ₁, Π₂, {tilde over (Π)}₁and {tilde over (Π)}₂are unitary. From W_x ^ε′(α′,τ₁′,τ₂′,τ₃′), it is possible to build a matrix {tilde over (W)}_x ^ε′(α′,τ₁′,τ₂′,τ₃′) depending solely on the sources of cyclical parameters (α,τ₁,τ₂,τ₃,ε) such that c_i ^ε(α,τ₁,τ₂,τ₃)≠0. To do this, the following computation is made: $\begin{matrix} {\tilde{W}}_{x}^{1} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = T_{1}^{H} W_{x}^{1} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) T_{1} = {Π_{1} (ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{C}}_{s}^{1} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{C}}_{s})}^{- 1 / 2} Π_{1}^{H} {\tilde{W}}_{x}^{2} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = T_{1}^{H} W_{x}^{2} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {\tilde{T}}_{1}^{*} = {Π_{1} (ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{C}}_{s}^{2} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{Π}}_{1}^{T} {\tilde{W}}_{x}^{3} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = {\tilde{T}}_{1}^{H} W_{x}^{3} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {\tilde{T}}_{1}^{*} = {{\tilde{Π}}_{1} (ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{C}}_{s}^{3} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{C}}_{s})}^{- 1 / 2} {\tilde{Π}}_{1}^{T} & (34) \end{matrix}$
Similarly, from W[0115] _x ^ε′(α′,τ₁′,τ₂′,τ₃′) it is possible to build a matrix ${\tilde{\tilde{W}}}_{x}^{ɛ^{'}} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'})$
that does not depend on the sources of cyclical parameters (α,τ[0116] ₁,τ₂,τ₃,ε) such as c_i ^ε(α,τ₁,τ₂,τ₃)=0: Other sources. To do this, the following computation is performed: $\begin{matrix} {\tilde{\tilde{W}}}_{x}^{1} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = T_{2}^{H} W_{x}^{1} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) T_{2} = {Π_{2} (ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2 <} {\tilde{\tilde{C}}}_{s}^{1} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} Π_{2}^{H} {\tilde{\tilde{W}}}_{x}^{2} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = T_{2}^{H} W_{x}^{2} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {\tilde{T}}_{2}^{*} = {Π_{2} (ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} {\tilde{\tilde{C}}}_{s}^{2} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} {\tilde{Π}}_{2}^{T} {\tilde{\tilde{W}}}_{x}^{3} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = {\tilde{T}}_{2}^{H} W_{x}^{3} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {\tilde{T}}_{2}^{*} = {{\tilde{Π}}_{2} (ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} {\tilde{\tilde{C}}}_{s}^{3} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} {\tilde{Π}}_{2}^{T} & (35) \end{matrix}$
where the matrix [0117] ${\tilde{\tilde{C}}}_{s}^{ɛ^{'}} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'})$
is a diagonal matrix with dimensions (P−P[0118] ₁)×(P−P₁) formed by diagonal elements c_i ^ε″(α′,τ₁′,τ₂′,τ₃′) such that the corresponding elements c_i ^ε(α,τ₁,τ₂,τ₃) are zero elements. The matrix ${\tilde{\tilde{C}}}_{s}$
sized (P−P[0119] ₁)×(P−P₁) is formed by the corresponding elements ci (1≦i≦P₁).
In particular, in the first version of this variant, it is possible to carry out the cyclical isolation in α′=0 and ε′=1. Writing W(τ[0120] ₁′,τ₂′,τ₃′)=W_x ^ε′(α′,τ₁′,τ₂′,τ₃′), {tilde over (C)}_s(α′,τ₁′,τ₂′,τ₃′)={tilde over (C)}_s ^ε(α′,τ ₁′,τ₂′,τ₃′) and ${\tilde{\tilde{C}}}_{s} (τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = {\tilde{\tilde{C}}}_{s}^{ɛ} (α^{'}, τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}),$
the equations (34) and (35) become: [0121]
{tilde over (W)} _x(τ₁′,τ₂′,τ₃′)=T ₁ ^H W _x(τ₁′,τ₂′,τ₃′)T ₁=Π₁(ε{tilde over (C)} _s)^−1/2 {tilde over (C)} _s(τ₁′,τ₂′,τ₃′)(ε{tilde over (C)} _s)^−1/2Π₁ ^H (36)
[0122] $\begin{matrix} {\tilde{\tilde{W}}}_{x} (τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) = T_{2}^{H} W_{x} (τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) T_{2} = {Π_{2} (ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} {\tilde{\tilde{C}}}_{s} (τ_{1}^{'}, τ_{2}^{'}, τ_{3}^{'}) {(ɛ {\tilde{\tilde{C}}}_{s})}^{- 1 / 2} Π_{2}^{H} & (37) \end{matrix}$
Identification Step [0123]
The equations (34) and (36) show that it is possible to identify the unitary matrices Π[0124] ₁and {tilde over (Π)}₁associated with the sources of cyclic parameters (α,τ₁,τ₂,τ₃) in carrying out the joint SVD (singular value decomposition) of the matrices {tilde over (W)}_x ^ε ^j(α^j,τ₁ ^j,τ₂ ^j,τ₃ ^j) for 1≦j≦K. Thus, to estimate the left unitary matrix, the joint diagonalizing of the matrices is performed:
{tilde over (W)} _x ^ε ^j(α^j,τ₁ ^j,τ₂ ^j,τ₃ ^j){tilde over (W)} _x ^ε ^j(α^j,τ₁ ^j,τ₂ ^j,τ₃ ^j)^Hfor 1≦j≦K (38)
and to estimate the right unitary matrix, the joint diagonalizing of the matrices is performed: [0125]
{tilde over (W)} _x ^ε ^j(α^j,τ₁ ^j,τ₂ ^j,τ₃ ^j){tilde over (W)} _x ^ε ^j(α^j,τ₁ ^j,τ₂ ^j,τ₃ ^j) for 1≦j≦K (39)
To estimate the unitary matrices Π[0126] ₂and {tilde over (Π)}₂associated with the sources not associated with the cyclical parameters (α,τ₁,τ₂,τ₃,ε), the joint SVD of the matrices ${\tilde{\tilde{W}}}_{x}^{ɛ^{j}} (α^{j}, τ_{1}^{j}, τ_{2}^{j}, τ_{3}^{j})$
for 1≦j≦K is performed in jointly diagonalizing the matrices [0127] ${\tilde{\tilde{W}}}_{x}^{ɛ^{j}} (α^{j}, τ_{1}^{j}, τ_{2}^{j}, τ_{3}^{j}) {{\tilde{\tilde{W}}}_{x}^{ɛ^{j}} (α^{j}, τ_{1}^{j}, τ_{2}^{j}, τ_{3}^{j})}^{H}$
and then the matrices [0128] ${{\tilde{\tilde{W}}}_{x}^{ɛ^{j}} (α^{j}, τ_{1}^{j}, τ_{2}^{j}, τ_{3}^{j})}^{H} {\tilde{\tilde{W}}}_{x}^{ɛ^{j}} (α^{j}, τ_{1}^{j}, τ_{2}^{j}, τ_{3}^{j}) .$
Knowing Π[0129] ₁and Π₂, from the equation (33), the unitary matrices U₁, U₂and U are deduced to the nearest permutation matrix, in performing:
U ₁ =T ₁Π₁ , U ₂ =T ₂Π₂ et U=[U ₁ U ₂] (40)
From the equations (9) and (29), it is possible to deduce the matrix A[0130] _Qto the nearest diagonal matrix and permutation matrix, such that:
T ^# U=[b _l . . . b _P ]=E _xΛ_x ^1/2 =A _Q(εC _s)^1/2ΛΠ (41)
where T[0131] ^# is the pseudo-inverse of the matrix T. Each column, b_l(1≦l≦P), of the matrix T^# U is associated with a vector μ_q|c_q|^1/2(α_q{circle over (×)}α_q*), 1≦q≦P, where μ_qis a complex scalar value such that |μ_q|=1. As a consequence, in converting each column b_lof the matrix T^# U into a matrix B. sized (N×N) such that B_l[i, j]=b_l((i−1)N+j) (1≦i, j≦N), it is deduced that:
B _l=μ_q |c| ^1/2α_qα_q ^Hfor (1≦l, q≦P) (42)
In this context, the direction vector α[0132] _qof the q^thsource is associated with the eigenvector of B_lassociated with the greatest eigenvalue.
Recapitulation of the First Version of the Cyclical Procedure [0133]
In short, the steps of the first version of the cyclical method are summarized here below and are applied to L observations x(lTe) (1≦l≦L) of the signals received on the sensors (T[0134] _e: sampling period).
Estimation [0135]
Step-1: The estimation of the matrices Q[0136] _xand {tilde over (Q)}_xfrom the L observations x(lTe). The estimation of these matrices will depend on the following assumptions:
Stationary and centered case: empirical estimator used in the reference [3]. [0137]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0138]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0139]
Whitening [0140]
Step 2: The eigen-element decomposition of the estimates of the matrices Q[0141] _xand {tilde over (Q)}_x. From these operations of decomposition, the estimation of the number of sources P and use of the P main eigenvalues such that: Q_x≈E_xΛ_xE_x ^Het {tilde over (Q)}_x={tilde over (E)}_x{tilde over (Λ)}_x{tilde over (E)}_x ^Hwhere Λ_xand {tilde over (Λ)}_xare diagonal matrices containing the P eigenvalues with the highest modulus and E_xand {tilde over (E)}_xare the matrices containing the associated eigenvectors.
Step 3: The building of the whitening matrices: T=({tilde over (Λ)}[0142] _x)^−1/2E_x ^Hand {tilde over (T)}=({tilde over (Λ)}_x)^−1/2{tilde over (E)}_x ^H
Step 4: The selection of the cyclical parameters (α,τ[0143] ₁,τ₂,τ₃,ε) and the estimation of the matrix Q_x ^ε(α,τ₁,τ₂,τ₃) from the L observations x(lTe). The estimation of this matrix will depend on the following assumptions on the signals:
Stationary and centered case: empirical estimator used in the reference [3]. [0144]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0145]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0146]
Step-5: The computation of a matrix W[0147] _x ^ε(α,τ₁,τ₂,τ₃) of (30) from matrices Q_x ^ε(α,τ₁,τ₂,τ₃), T and {tilde over (T)}. After singular value decomposition W_x ^ε(α,τ₁,τ₂,τ₃), the determining of the unitary matrices T₁and {tilde over (T)}₁associated with the non-zero singular values and T₂and {tilde over (T)}₂associated with the zero singular values.
Selection [0148]
Step-6: The selection of the K triplets of delays (τ[0149] ₁ ^k,τ₂ ^k,τ₃ ^k) where |τ₁ ^k|+|τ₂ ^k|+τ₃ ^k|≠0.
Estimation [0150]
Step-7: The estimation of the K matrices Q[0151] _x(τ₁ ^k,τ₂ ^k,τ₃ ^k) of (2). As in the step-1 this estimation will depend on the assumptions made on the signal such as:
Stationary and centered case: empirical estimator used in the reference [3]. [0152]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0153]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0154]
Identification [0155]
Step-8: The computation of the matrices T[0156] ₁Q_x(τ₁ ^k,τ₂ ^k,τ₃ ^k) T₁ ^Hand the estimation of the unitary matrix U₁(associated with the cyclical parameters (α,τ₁,τ₂,τ₃,ε)) in carrying out the joint diagonalizing of the K matrices: T₁Q_x(α,τ₁ ^k,τ₂ ^k,τ₃ ^k)T₁ ^H.
Step-9: The computation of the matrices T[0157] ₂Q_x(τ₁ ^k,τ₂ ^k,τ₃ ^k) T₂ ^Hand the estimation of the unitary matrix U₂(associated with the other sources) in carrying out the joint diagonalizing of the K matrices T₂Q_x(τ₁ ^k,τ₂ ^k,τ₃ ^k)T₂ ^H.
Step-10: The computation of the unitary matrix U in performing: U=[U[0158] ₁U₂]
Step-11: The computation of T[0159] ^#U=[b_l. . . b_P] and the building of the matrices B_lwith a dimension (N×N) from the columns b_lof T^#U.
Step-12: The estimation of the signatures a, (1≦q≦P) of the P sources in applying a decomposition into elements in each matrix B[0160] _l.
Recapitulation of the Second Version of the Cyclical Procedure [0161]
The steps of the second version of the cyclical FORBIUM method are summarized here below and are applied to L observations x(lTe) (1≦l≦L) of the signals received on the sensors (T[0162] _e: sampling period).
Estimation [0163]
Step-1: The estimation of the matrices Q[0164] _xand {tilde over (Q)}_xfrom the L observations x(lTe). The estimation of these matrices will depend on the following assumptions:
Stationary and centered case: empirical estimator used in the reference [3]. [0165]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0166]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0167]
Step-2: The eigen-element decomposition of the matrices Q[0168] _xand {tilde over (Q)}_x. From these operations of decomposition, the estimation of the number of sources P and the use of the P main eigenvalues such that: Q_x≈E_xΛ_xE_x ^Hand {tilde over (Q)}_x={tilde over (E)}_x{tilde over (Λ)}_x{tilde over (E)}_x ^Hwhere Λ_xand {tilde over (Λ)}_xare diagonal matrices containing the P eigen values with the highest modulus and E_xand {tilde over (E)}_xare the matrices containing the P associated eigen vectors.
Step-3: The building of the whitening matrices: T=(Λ[0169] _x)^−1/2E_x ^Hand {tilde over (T)}=({tilde over (Λ)}_x)^−1/2{tilde over (E)}_x ^H
Step-4: The selection of the cyclical parameters (α, τ[0170] ₁,τ₂,τ₃,ε) and the estimation of the matrix Q_x ^ε(α,τ₁,τ₂,τ₃) from the L observations x((lTe). The estimation of this matrix will depend on the following assumptions on the signals:
Stationary and centered case: empirical estimator used in the reference [3]. [0171]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0172]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0173]
Step-5: The computation of a matrix W[0174] _x ^ε(α,τ₁,τ₂,τ₃) of (30) from the matrices Q_x ^ε(α,τ₁,τ₂,τ₃), T and {tilde over (T)}. After a singular value decomposition of W_x ^ε(α,τ₁,τ₂,τ₃), the determining of the unitary matrices T₁and {tilde over (T)}₁associated with the non-zero significant values and T₂and {tilde over (T)}₂associated with the zero singular values.
Step-6: The selection of K sets of parameters (α[0175] ^k,τ₁ ^k,τ₂ ^k,τ₃ ^k).
Step-7: The estimation of the K matrices A[0176] _x ^εk(α^k,τ₁ ^k,τ₂ ^k,τ₃ ^k) of (19). As in the step-1 this estimation will depend on the assumptions made on the signal such as:
Stationary and centered case: empirical estimator used in the reference [3]. [0177]
Cyclostationary and centered case: estimator implemented in the reference [10]. [0178]
Cyclostationary and non-centered case: estimator implemented in the reference [11]. [0179]
Step-8: The computation of the matrices T[0180] ₁Q_x ^εk(α^k,τ₁ ^k,τ₂ ^k,τ₃ ^k) T₁ ^Hand the estimation of the unitary matrix U₁or Ũ₁(associated with the cyclical parameters (α,τ₁,τ₂,τ₃,ε)) in carrying out the joint SVD of the K matrices: T₁Q_x ^εk(α^k,τ₁ ^k,τ₂ ^k,τ₃ ^k)T₁ ^H.
Step-9: The computation of the matrices T[0181] ₂Q_x ^εk(α^k,τ₁ ^k,τ₂ ^k,τ₃ ^k)T₂ ^Hand the estimation of the unitary matrix U₂or Ũ₂(associated with the cyclical parameters (α,τ₁,τ₂,τ₃,ε)) in carrying out the joint SVD of the K matrices: T₂Q_x ^εk(α^k,τ₁ ^k,τ₂ ^k,τ₃ ^k)T₂ ^H.
Step-10: The computation of the unitary matrix U in performing: U=[U[0182] ₁U₂]
Step-11: The computation of T[0183] ^#U=[b_l. . . b_P] and the building of the matrices B_lwith a dimension (N×N) from the columns b_lof T^#U. Step-12: The estimation of the signatures aq (1≦q≦P) of the P sources in applying a decomposition into elements on each matrix B_l.
Bibliography [0184]
[1] A. Belouchrani, K. Abed—Meraim, J. F. Cardoso, E. Moulines, “A blind source separation technique using second-order statistics”, [0185] IEEE Trans. Sig. Proc., Vol.45, N°2, pp. 434-444, February 1997.
[2] J F. Cardoso, “Super-symmetric decomposition of the fourth order cumulant tensor”, [0186] ICASSP 1991.
[3] J. F. Cardoso, A. Souloumiac, “Blind beam forming for non-gaussian signals”, [0187] IEE Proceedings-F, Vol.140, N°6, pp. 362-370, December 1993.
[4] P. Chevalier, G. Benoit, A. Ferréol <<Direction finding after blind identification of sources steering vectors: The blind-maxcor and blind-MUSIC methods>>, EUSIPCO, Trieste, pp 2097-2100, 1996 [0188]
[5] P. Chevalier, “Optimal separation of independent narrow-band sources: concept and performance”, [0189] Sig. Proc., Elsevier, Vol.73, pp 27-47, 1999.
[6] P. Chevalier, A. Ferréol , “On the virtual array concept for the fourth-order direction finding problem”, [0190] IEEE trans on signal processing, Vol.47, N°9, pp. 2592-2595, September 1999.
[7] P. Comon, “Independent component analysis—a new concept?”, [0191] Sig. Proc., Elsevier, Vol.36, N°3, April 1994.
[8] P. Comon “Blind channel identification and extraction of more sources than sensors”, [0192] SPIE Conf Adv Proc VIII , San Diego, July 1998.
[9] L. De Lathauwer, B. De Moor , J. Vandewalle , “ICA techniques for more sources than sensors”, Proc [0193] IEEE Processing Workshop on Higher Order Statistics, Caesarea, Israel , June 1999.
[10] A. Ferréol, P. Chevalier, “On the behavior of current second and higher order blind source separation methods for cyclostationary sources”, [0194] IEEE Trans. Sig. Proc., Vol.48, N°6, pp. 1712-1725, June 2000.
[11] A. Ferréol, P. Chevalier, L. Albera “Higher order blind separation of non zero-mean cyclostationnary sources”, ([0195] EUSPICO 2002), Toulouse, September 2002.
[12] A. Ferréol, P. Chevalier, L. Albera, “Procédé de traitement d'antennes sur des signaux cyclostationnaires potentiellement non centrés,” (Method of antenna processing on potentially non-centered cyclostationary signals), patent, May 2002. [0196]
[13] A. Ferréol, P. Chevalier “Higher order blind source separation using the cyclostationarity property of the signals”, ICASSP Munich, Vol 5, pp4061-4064, 1997 [0197]
[14] A. Ferréol, P. Chevalier “Higher order blind source separation using the cyclostationarity property of the signals”, ICASSP Munich, Vol 5, pp4061-4064, 1997 [0198]
[10] A. Taleb “An algorithm for the blind identification of N independent signal with 2 sensors”, 16[0199] ^th symposium on signal processing and its applications (ISSPA 2001), Kuala-Lumpur, August 2001.

Claims

1. A method of fourth-order, blind identification of two sources in a system including a number of sources P and a number N of reception sensors receiving the observations, the sources having different tri-spectra, comprising the following steps:

a) fourth-order whitening of the observations received on the reception sensors in order to orthonormalize the direction vectors of the sources in the matrices of quadricovariance of the observations used,

b) joint diagonalizing of several whitened matrices of quadricovariance to identify the spatial signatures of the sources.

2. The method according to claim 1, wherein the observations used correspond to the time-domain averaged matrices of quadricovariance defined by:

\begin{matrix} Q_{x} (τ_{1}, τ_{2}, τ_{3}) = \sum_{p = 1}^{P} c_{p} (τ_{1}, τ_{2}, τ_{3}) (a_{p} \otimes a_{p}^{*}) {(a_{p} \otimes a_{p}^{*})}^{H} & (4a) \\ = A_{Q} C_{s} (τ_{1}, τ_{2}, τ_{3}) A_{Q}^{H} & (4b) \end{matrix}

where A_Qis a matrix with a dimension (N²×P) defined by A_Q=[(α₁{circle over (×)}α₁*), . . . , (α_p{circle over (×)}α_p*)], C_s(τ₁,τ₂,τ₃) is a diagonal matrix with a dimension (P×P) defined by C_s(τ₁,τ₂,τ₃)=diag[c_l(τ₁,τ₂,τ₃), . . . , c_p(τ₁,τ₂,τ₃)] and c_p(τ₁,τ₂,τ₃) is defined by:

c _p(τ₁,τ₂,τ₃)=<Cum(s _p(t), s _p(t−τ ₁)*, s _p(t−τ ₂)*, s _p(t−τ ₃))> (5)

3. The method according to claim 2, comprising the following steps:

Step 1: estimating, through Q;{circumflex over ( )}_x, of the matrix Q_x, from the L observations x(IT_e) using a non-skewed and asymptotically consistent estimator.

Step 2: eigen-element decomposition of Q;{circumflex over ( )}_x, the estimation of the number of sources P and the limiting of the eigen-element decomposition to the P main components: Q;{circumflex over ( )}_x≈E;{circumflex over ( )}_xΛ;{circumflex over ( )}_xE;{circumflex over ( )}_x ^H, where Λ;{circumflex over ( )}_xis the diagonal matrix containing the P eigenvalues with the highest modulus and E;{circumflex over ( )}_xis the matrix containing the associated eigenvectors.

Step 3: building of the whitening matrix: T;{circumflex over ( )}=(Λ;{circumflex over ( )}_x)^−1/2E;{circumflex over ( )}_x ^H.

Step 4: selecting K triplets of delays (τ₁ ^k,τ₂ ^k,τ₃ ^k) where |τ₁ ^k|+|τ₂ ^k|+|τ₃ ^k≠0.

Step 5: estimating, through Q;{circumflex over ( )}_x(τ₁ ^k,τ₂ ^k,τ₃ ^k), of the K matrices Q_x(τ₁ ^k,τ₂ ^k,τ₃ ^k).

Step 6: computing of the matrices T;{circumflex over ( )} Q;{circumflex over ( )}_x(τ₁ ^k,τ₂ ^k,τ₃ ^k) T;{circumflex over ( )}^Hand the estimation, by U,{circumflex over ( )}_sol, of the unitary matrix U_solby the joint diagonalizing of the K matrices T;{circumflex over ( )} Q;{circumflex over ( )}_x(τ₁ ^k,τ₂ ^k,τ₃ ^k) T;{circumflex over ( )}^H

Step 7: computing T;{circumflex over ( )}^#U;{circumflex over ( )}_sol=[b;{circumflex over ( )}₁. . . b;{circumflex over ( )}_P] and the building of the matrices B;{circumflex over ( )}₁sized (N×N).

Step 8: estimating, through α;{circumflex over ( )}_P, of the signatures aq (1≦q≦P) of the P sources in applying a decomposition into elements on each matrix B;{circumflex over ( )}₁.

4. The method according to claim 1, comprising evaluating quality of the identification of the associated direction vector in using a criterion:

D(A, Â)=(α₁, α₂, . . . , α_P) (16)

where

\begin{matrix} α_{p} = \min_{1 \leq i \leq P} [d (a_{p}, {\hat{a}}_{i})] & (17) \end{matrix}

and where d(u,v) is the pseudo-distance between the vectors u and v, such that:

\begin{matrix} d (u, v) = 1 - \frac{{\langle u^{H} v \rangle}^{2}}{(u^{H} u) (v^{H} v)} & (18) \end{matrix}

5. The method according to claim 1, a fourth-order cyclical after the step a) of fourth-order whitening.

6. The method according to claim 5, wherein the identification step is performed in using fourth-order statistics.

7. The method according to claim 1 wherein the number of sources P is greater than or equal to the number of sensors.

8. The method according to claim 1, comprising goniometry using the identified signature of the sources.

9. The method according to claim 1, comprising spatial filtering after the identified signature of the sources.

10. The use of the method according to claim 1, for use in a communications network.

11. The method according to claim 2, comprising evaluating quality of the identification of the associated direction vector in using a criterion

D (A, \hat{A}) = (α_{1}, α_{2}, \dots, α_{P}) where

α_{p} = \min_{1 \leq i \leq P} [d (a_{p}, {\hat{a}}_{i})]

and where d(u,v) is the pseudo-distance between the vectors u and v, such that:

d (u, v) = 1 - \frac{{\langle u^{H} v \rangle}^{2}}{(u^{H} u) (v^{H} v)}

12. The method according to claim 3, comprising evaluating quality of the identification of the associated direction vector in using a criterion

D (A, \hat{A}) = (α_{1}, α_{2}, \dots, α_{P}) where

α_{p} = \min_{1 \leq i \leq P} [d (a_{p}, {\hat{a}}_{i})]

and where d(u,v) is the pseudo-distance between the vectors u and v, such that:

d (u, v) = 1 - \frac{{\langle u^{H} v \rangle}^{2}}{(u^{H} u) (v^{H} v)}

13. The method according to claim 2, a fourth-order cyclical after the step a) of fourth-order whitening.

14. The method according to claim 2, wherein the identification step is performed in using fourth-order statistics.

15. The method according to claim 2, wherein the number of sources P is greater than or equal to the number of sensors.

16. The method according to claim 2, comprising goniometry using the identified signature of the sources.

17. The method according to claim 3, a fourth-order cyclical after the step a) of fourth-order whitening.

18. The method according to claim 3, wherein the identification step is performed in using fourth-order statistics.

19. The method according to claim 3, wherein the number of sources P is greater than or equal to the number of sensors.

20. The method according to claim 3, comprising goniometry using the identified signature of the sources.