U.S. patent application number 10/813673 was filed with the patent office on 2004-12-23 for method for the higher-order blind identification of mixtures of sources.
Invention is credited to Albera, Laurent, Chevalier, Pascal, Comon, Pierre, Ferreol, Anne.
Application Number | 20040260522 10/813673 |
Document ID | / |
Family ID | 32982145 |
Filed Date | 2004-12-23 |
United States Patent
Application |
20040260522 |
Kind Code |
A1 |
Albera, Laurent ; et
al. |
December 23, 2004 |
Method for the higher-order blind identification of mixtures of
sources
Abstract
A method for the blind identification of sources within a system
having P sources and N receivers comprises at least one step for
the identification of the matrix of the direction vectors of the
sources from the information proper to the direction vectors
a.sub.p of the sources contained redundantly in the m=2q order
circular statistics of the vector of the observations received by
the N receivers. Application to a communications network.
Inventors: |
Albera, Laurent; (Boulogne
Billancourt, FR) ; Ferreol, Anne; (Colombes, FR)
; Chevalier, Pascal; (Courbevoie, FR) ; Comon,
Pierre; (Peymeinade, FR) |
Correspondence
Address: |
LOWE HAUPTMAN GILMAN & BERNER, LLP
1700 DIAGNOSTIC ROAD, SUITE 300
ALEXANDRIA
VA
22314
US
|
Family ID: |
32982145 |
Appl. No.: |
10/813673 |
Filed: |
March 31, 2004 |
Current U.S.
Class: |
702/189 |
Current CPC
Class: |
G06K 9/6243 20130101;
G06F 17/16 20130101; G06F 17/12 20130101; G06F 17/18 20130101 |
Class at
Publication: |
702/189 |
International
Class: |
C12P 021/06 |
Foreign Application Data
Date |
Code |
Application Number |
Apr 1, 2003 |
FR |
03 04041 |
Claims
1. A method of blind identification of sources within a system
including P sources and N receivers, comprising the steps of
identifying the matrix of direction vectors of the sources from the
information proper to the direction vectors a.sub.p of the sources
contained redundantly in the m=2q order circular statistics of the
vector of the observations received by the N receivers.
2. The method according to claim 1, wherein m=2q order circular
statistics are expressed according to a full-rank diagonal matrix
of the autocumulants of the sources and a matrix representing the
juxtaposition of the direction vectors of the sources as
follows:C.sub.m,x=A.sub.q .zeta..sub.m,s A.sub.q.sup.Hwhere 56 m ,
s = diag ( [ C 1 , 1 , , 1 , s 1 , 1 , , 1 , , C P , P , , P , s P
, P , , P ] ) is the full-rank diagonal matrix of the m=2q order
autocumulants 57 C p , p , , p , s p , p , , p des sources, sized
(P.times.P), and where 58 A q = [ a 1 ( q - 1 ) a 1 * a p ( q - 1 )
a p * ] ,sized (N.sup.q.times.P) and assumed to be of full rank,
represents the juxtaposition of the P column vectors 59 [ a p ( q -
1 ) a p * ] .
3. The method according to claim 1, further comprising the
following steps: a): the building, from the different observation
vectors x(t), of an estimate .sub.m,x of the matrix of statistics
C.sub.m,x of the observations, b): decomposing a singular value of
the matrix .sub.m,x, and deducing therefrom of an estimate P;
{circumflex over ( )} of the number of sources P and a square root
60 C ^ m , x 1 / 2 of .sub.m,x, in taking 61 C ^ m , x 1 / 2 = E s
L s 1 / 2 where .vertline...vertline. designates the absolute value
operator, where L.sub.s and E.sub.s are respectively the diagonal
matrix of the P; {circumflex over ( )} greatest real eigenvalues
(in terms of absolute value) of .sub.m,x and the matrix of the
associated orthonormal eigenvectors; c): extracting, from the
matrix 62 C ^ m , x 1 / 2 = [ 1 T , , N T ] T ,of the N matrix
blocks .GAMMA..sub.n: each block .GAMMA..sub.n sized
(N.sup.(q-1).times.P) being constituted by the N.sup.(q-1)
successive rows of 63 C ^ m , x 1 / 2 starting from the
"N.sup.(q-1)(n-1)+1"th row; d): building of the N(N-1) matrices
.THETA..sub.n1,n2 defined, for all
1.ltoreq.n.sub.1.noteq.n.sub.2.ltoreq.N, by
.THETA..sub.n1,n2=.GAMMA..sub- .n1.sup.# .GAMMA..sub.n2 where #
designates the pseudo-inversion operator; e): determining of the
matrix V.sub.sol, resolving the problem of the joint
diagonalization of the N(N-1) matrices .THETA..sub.n1,n2; f): for
each of the P columns b.sub.p of A; {circumflex over ( )}.sub.q,
the extraction of the K=N.sup.(q-2)vectors b.sub.p(k) stacked
beneath one another in the vector b.sub.p=[b.sub.p(1).sup.T,
b.sub.p(2).sup.T, . . . , b.sub.p(K).sup.T].sup.T; g): converting
said column vectors b.sub.p(k) sized (N.sup.2.times.1) into a
matrix B.sub.p(k) sized (N.times.N); h): joint singular value
decomposition or joint diagonalization of the K=N.sup.(q-2)
matrices B.sub.p(k) in retaining therefrom, as an estimate of the
column vectors of A, of the eigenvector common to the K matrices
B.sub.p(k) associated with the highest eigenvalue (in terms of
modulus); i): repetition of the steps f) to h) for each of the P
columns of A; {circumflex over ( )}.sub.q for the estimation,
without any particular order and plus or minus a phase, of the P
direction vectors a.sub.p and therefore the estimation, plus or
minus a unitary trivial matrix, of the mixture matrix A.
4. The method according to claim 1, wherein the number of sensors N
is greater than or equal to the number of sources P and comprising
a step of extraction of the sources, consisting of the application
to the observations x(t) of a filter built by means of the estimate
A; {circumflex over ( )} of A.
5. The method according to claim 2, wherein C.sub.m,x is equal to
the matrix of quadricovariance Qx and wherein m=4.
6. The method according to claim 2, wherein C.sub.m,x is equal to
the matrix of hexacovariance Hx and wherein m=6.
7. The method according to claim 1, comprising a step for the
evaluation of the quality of the identification of the associated
direction vector in using a criterion such that:D(A,
)=(.alpha..sub.1, .alpha..sub.2, . . . , .alpha..sub.P)where 64 p =
min 1 i P [ d ( a p , a ^ i ) ] [ [ ( 17 ) ] ] and where d(u,v) is
the pseudo-distance between the vectors u and v, such that: 65 d (
u , v ) = 1 - u H v 2 ( u H u ) ( v H v ) [ [ ( 18 ) ] ]
8. The use of the method according to claim 1, for use in a
communications network.
9. A use of the method according to claim 1 for goniometry using
identified direction vectors.
10. The method according to claim 2, wherein the number of sensors
N is greater than or equal to the number of sources P and wherein
the method comprises a step of extraction of the sources,
consisting of the application to the observations x(t) of a filter
built by means of the estimate A; {circumflex over ( )} of A.
11. The method according to claim 3, wherein the number of sensors
N is greater than or equal to the number of sources P and wherein
the method comprises a step of extraction of the sources,
consisting of the application to the observations x(t) of a filter
built by means of the estimate A; {circumflex over ( )} to A .
12. The method according to claim 3, wherein C.sub.m,x is equal to
the matrix of quadricovariance Qx and wherein m=4.
13. The method according to claim 4, wherein C.sub.m,x is equal to
the matrix of quadricovariance Qx and wherein m=4.
14. The method according to claim 3, wherein C.sub.m,x is equal to
the matrix of hexacovariance Hx and wherein m=6.
15. The method according to claim 4, wherein C.sub.m,x is equal to
the matrix of hexacovariance Hx and wherein m=6.
16. The method according to claim 2, comprising a step for the
evaluation of the quality of the identification of the associated
direction vector in using a criterion such that:D(A,
)=(.alpha..sub.1, .alpha..sub.2, . . . , .alpha..sub.P)where 66 p =
min 1 i P [ d ( a p , a ^ i ) ] and where d(u,v) is the
pseudo-distance between the vectors u and v, such that: 67 d ( u ,
v ) = 1 - u H v 2 ( u H u ) ( v H v )
17. The method according to claim 3, comprising a step for the
evaluation of the quality of the identification of the associated
direction vector in using a criterion such that:D(A,
)=(.alpha..sub.1, .alpha..sub.2, . . . , .alpha..sub.P)where 68 p =
min 1 i P [ d ( a p , a ^ i ) ] and where d(u,v) is the
pseudo-distance between the vectors u and v, such that: 69 d ( u ,
v ) = 1 - u H v 2 ( u H u ) ( v H v )
18. The method according to claim 4, comprising a step for the
evaluation of the quality of the identification of the associated
direction vector in using a criterion such that:D(A,
)=(.alpha..sub.1, .alpha..sub.2, . . . , .alpha..sub.P)where 70 p =
min 1 i P [ d ( a p , a ^ i ) ] and where d(u,v) is the
pseudo-distance between the vectors u and v, such that: 71 d ( u ,
v ) = 1 - u H v 2 ( u H u ) ( v H v )
19. The method according to claim 5, comprising a step for the
evaluation of the quality of the identification of the associated
direction vector in using a criterion such that:D(A,
)=(.alpha..sub.1, .alpha..sub.2, . . . , .alpha..sub.P)where 72 p =
min 1 i P [ d ( a p , a ^ i ) ] and where d(u,v) is the
pseudo-distance between the vectors u and v, such that: 73 d ( u ,
v ) = 1 - u H v 2 ( u H u ) ( v H v )
20. The method according to claim 6, comprising a step for the
evaluation of the quality of the identification of the associated
direction vector in using a criterion such that:D(A,
)=(.alpha..sub.1, .alpha..sub.2, . . . , .alpha..sub.P)where 74 p =
min 1 i P [ d ( a p , a ^ i ) ] and where d(u,v) is the
pseudo-distance between the vectors u and v, such that: 75 d ( u ,
v ) = 1 - u H v 2 ( u H u ) ( v H v )
Description
BACKGROUND OF THE INVENTION
[0001] 1. Field of the Invention
[0002] The invention relates especially to a method of fourth-order
and higher order self-learned (or blind) separation of sources from
reception with N receivers or sensors (N.gtoreq.2), this method
exploiting no a priori information on sources or wavefronts, and
the sources being P cyclostationary (deterministic or stochastic,
analog or digital sources with linear or non-linear modulation) and
statistically independent sources,
[0003] It can be applied for example in the field of radio
communications, space telecommunications or passive listening to
these links in frequencies ranging for example from VLF to EHF.
[0004] It can also be applied in fields such as astronomy,
biomedicine, radar, speech processing, etc.
[0005] 2. Description of the Prior Art
[0006] The blind separation of sources and, more particularly,
independent component analysis (ICA) is currently arousing much
interest. Indeed, it can be used in many applications such as
telecommunications, speech processing or again biomedicine.
[0007] For example, in antenna processing, if signals sent from a
certain number of sources are received at an array of receivers and
if, for each source, the timing spread of the channels associated
with a different receivers is negligible as compared with the
timing symbol time, then an instantaneous mixture of the signals
sent from the sources is observed on said receivers.
[0008] The blind separation of sources is aimed especially at
restoring the sources assumed to be statistically independent, and
this is done solely on the basis of the observations received by
the receivers.
[0009] Depending on the application, it is possible to retrieve
only the instantaneous mixture, namely the direction vectors of the
sources. This is the case for example with goniometry where said
mixture carries all the information needed for the angular
localization of the sources by itself: the term used then is "blind
identification of mixtures".
[0010] For other applications such as transmission, it is necessary
to retrieve the signals sent from the sources: the expression used
then is separation or again blind or self-learned extraction of
sources.
[0011] Certain prior art techniques seek to carry out a
second-order decorrelation (as can be seen in factor analysis with
principal component analysis (PCA).
[0012] ICA, for its part, seeks to reduce the statistical
dependence of the signals also at the higher orders. Consequently,
ICA enables the blind identification of the instantaneous mixture
and thereby the extraction of the signals sent from the sources,
not more than one of which is assumed to be Gaussian. At present,
this is possible only in complying with certain assumptions: the
noisy mixture of the sources must be linear and furthermore
over-determined (the number of sources P must be smaller than or
equal to the number of receivers N).
[0013] While Comon was the first to introduce the ICA concept and
propose a solution, COM2 in the reference [1] (the different
references are brought together at the end of the description)
maximized a contrast based on fourth-order cumulants, Cardoso and
Souloumiac [2], for their part developed a matrix approach, better
known as JADE, and thus created the joint diagonalization
algorithm.
[0014] Some years later, Hyvarinen et al. presented the FastICA
method, initially for real signals [3], and then in complex cases
[4]. This method introduces a contrast-optimizing algorithm called
the fixed-point algorithm.
[0015] Comon has proposed a simple solution, COM1 [5], to contrast
optimization presented in [6].
[0016] Although these methods perform very well under the
assumptions stated here above, they may nevertheless be greatly
disturbed by the presence of unknown noise, whether Gaussian or
not, that is spatially correlated and inherent in certain
applications such as HF radio communications.
[0017] Furthermore, as stated further above, the above methods are
designed only to process over-determined mixtures of sources. Now
in practice, for example in radio communications, it is not rare to
have reception from more sources than receivers, especially if the
reception bandwidth is great. We then have what are called
under-determined mixtures (P>N).
[0018] Several algorithms have been developed already in order to
process mixtures of this type. Some of them tackle the difficult
problem of the extraction of sources [7-8] when the mixture is no
longer linearly inverted, while others deal with the indication of
the mixture matrix [7] [9-12.
[0019] The methods proposed in [9-11] exploit only fourth-order
statistics while the method presented in [12] relies on the use of
the characteristic second function of the observations, in other
words on the use of non-zero cumulants of any order. As for the
method used in [7], it relies on the conditional maximization of
probability, in this case that of data conditional on the mixture
matrix.
[0020] While these methods perform well, they have drawbacks in the
operational context.
[0021] Thus, the method [9] is difficult to implement and does not
ensure the identification of the direction vectors of sources of
the same kurtosis. The methods [10] and [11] for their part cannot
be used to identify the direction vectors of circular sources. The
method [10], called S3C2, for its part confines the user in a
configuration of three sources and two receivers, ruling out any
other scenario. The method [7] authorizes the identification of
four speech signals with only two receivers. However the samples
observed must be temporally independent and each source must have a
sparse density of probability. Finally, the method [12] is
applicable only in the case of real sources, which is highly
restrictive especially in digital communications. Furthermore, the
algorithm depends greatly on the number of sources, and there is
nothing today to prove that a poor estimation of this parameter
will not impair the performance of the method.
SUMMARY OF THE INVENTION
[0022] The present invention offers a novel approach relying
especially on the exploitation of the totality or practically the
totality of the information proper to the direction vectors a.sub.p
of the sources, contained redundantly in the matrix representing
m=2q order circular statistics of the vector of the complex
envelopes of the signals at output of the receivers.
[0023] The invention relates to a method for the blind
identification of sources within a system comprising P sources and
N receivers. It comprises at least one step for the identification
of the matrix of the direction vectors of the sources from the
information proper to the direction vectors a.sub.p of the sources
contained redundantly in the m=2q order circular statistics of the
vector of the observations received by the N receivers.
[0024] The m=2q order circular statistics are expressed for example
according to a full-rank diagonal matrix of the autocumulants of
the sources and a matrix representing the juxtaposition of the
direction vectors of the sources as follows:
C.sub.m,x=A.sub.q.zeta..sub.m,sA.sub.q.sup.H (11)
[0025] where 1 m , s = diag ( [ C 1 , 1 , , 1 , s 1 , 1 , , 1 , , C
P , P , , P , s P , P , , P ] )
[0026] is the full-rank diagonal matrix of the m=2q order
autocumulants 2 C p , p , , p , s p , p , , p
[0027] des sources, sized (P.times.P), and where
A.sub.q=[a.sub.1.sup.{cir- cle over (X)}(q-1){circle over
(X)}a.sub.1* . . . a.sub.p.sup.{circle over (X)}(q-1){circle over
(X)}a.sub.p*], sized (N.sub.q.times.P) and assumed to be of full
rank, represents the juxtaposition of the P column vectors
[a.sub.p.sup.{circle over (X)}(q-1){circle over (X)}a.sub.p*].
[0028] The method of the invention is used in a communications
network and/or for goniometry using identified direction
vectors.
[0029] The invention has especially the following advantages:
[0030] It enables the blind identification of instantaneous
mixtures, both over-determined (where the number of sources is
smaller than or equal to the number of receivers) and
under-determined (where the number of sources is greater the number
of receivers) as well as the blind extraction of the sources in the
over-determined case;
[0031] At the m=2q order, which is an even-parity value where
q.gtoreq.2, the procedure called BIOME (Blind Identification of
Over and underdetermined Mixtures of sources) can process up to
P=N.sup.(q-1) sources using the array with N different receivers,
once the m order autocumulants of the sources have the same
sign;
[0032] An application of the fourth order method known as ICAR
(Independent Component Analysis using Redundancies in the
quadricovariance), enables the blind identification of
over-determined (P.ltoreq.N) instantaneous mixtures of sources and
their blind extraction, in a manner that proves to be robust in the
presence of a spatially correlated unknown Gaussian noise, once the
sources have same-sign kurtosis (fourth-order standardized
autocumulants);
[0033] an application of the BIOME method at the sixth-order level,
called BIRTH (Blind Identification of mixtures of sources using
Redundancies in the daTa Hexacovariance matrix), enables the blind
identification of instantaneous mixtures, both over-determined
(P.ltoreq.N) and under-determined (P>N), of sources, as well as
the blind extraction of the sources in the over-determined case.
The BIRTH method has the capacity to process up to N.sup.2 sources
from an array with N different receivers, once the sixth-order
autocumulants of the sources have a same sign.
BRIEF DESCRIPTION OF THE DRAWINGS
[0034] Other characteristics of advantages of the method according
to the invention shall appear more clearly from the following
description from a non-restrictive example of an embodiment and the
appended figures, of which:
[0035] FIG. 1 is a drawing exemplifying an m-order implementation
of the method,
[0036] FIGS. 2, 3 and 4 show results of simulation of a
fourth-order implementation of the method,
[0037] FIGS. 5, 6 and 8 show results of simulation of a sixth-order
implementation of the method.
[0038] The following examples are given for the identification
and/or extraction of sources in an array comprising an array
antenna comprising N receivers. It is assumed that this antenna
receives a noisy mixture of signals from P statistically
independent sources for example in narrow band (NB).
[0039] On the basis of these assumptions, the vector x(t) of the
complex envelopes of the signals at output of the receivers is
written, at the instant t 3 x ( t ) = p = 1 P s p ( t ) a p + v ( t
) = A s ( t ) + v ( t ) ( 1 )
[0040] where v(t) is the noise vector, assumed to be centered,
Gaussian, spatially white and unknown,
[0041] s.sub.p(t) and a.sub.p correspond respectively to the
complex, narrow-band, cyclostationary and cycloergodic envelope BE
(with a possible residue of a carrier as the case may be) and to
the direction vector of the source p,
[0042] s(t) is the vector whose components are the signals
s.sub.p(t) and A is the matrix (N.times.P) whose columns are the
vectors a.sub.p.
[0043] Furthermore, the method generally described here below for
the m=2q (q.gtoreq.2) order uses the following assumptions,
numbered H.sub.1-4:
[0044] H1: At any instant t, the sources with complex values
s.sub.p(t) are cyclostationary, cycloergodic and mutually
decorrelated at the m order;
[0045] H2: At any instant t, the components v.sub.n(t) of the noise
are stationary, ergodic, Gaussian and circular;
[0046] H3: A any instant t, s(t) and v(t) are statistically
independent;
[0047] H4: the m order autocumulants of the sources are not zero
and have the same sign.
[0048] With the above assumptions, for a given even-parity order
m=2q, the problem of the blind identification of instantaneous
mixtures of sources consists in finding the matrix A through the
exploitation of certain m order statistics of the observations.
This matrix A is found to the nearest trivial matrix (a trivial
matrix has the form .LAMBDA. .PI. where .LAMBDA. is an invertible
matrix and .PI. is a permutation matrix).
[0049] The blind separation (extraction) of sources consists
especially in determining the separator that is linear and
invariant in time (LIT), W, with a dimension (N.times.P), the
output vector of which has a dimension (P.times.1),
y(t)=W.sup.Hx(t) (2)
[0050] corresponds, plus or minus a trivial matrix, to the best
estimate, (t), of the vector s(t), where the symbol .sup.H
signifies a conjugate transpose.
[0051] The separator W is defined to the nearest trivial matrix
inasmuch as neither the value of the output power values nor the
order in which they are stored changes the quality of restitution
of the sources.
[0052] Before any explanation of the steps of the method according
to the invention, a few reminders are given on the statistics of
the observations.
[0053] Statistics of the Observations
[0054] The method according to the invention uses especially 2, 4,
6, . . . , m even-parity circular statistics of the
observations.
[0055] According to the prior art described in the reference [15],
the expression of the m order cumulants as a function of the
moments of an order lower than m can be simplified.
[0056] Let 4 G d , , e , x f , , g
[0057] be a scalar quantity with a complex value depending in the q
lower indices d, . . . , e and the q higher indices f, . . . , g
having values in {1, 2, . . . , N}. The quantity 5 G d , , e , x f
, , g
[0058] then verifies the following three symmetries:
[0059] any permutation between the lower indices of 6 G d , , e , x
f , , g
[0060] does not modify the value of 7 G d , , e , x f , , g :
[0061] for example, for q=2, 8 G e , d , x f , g = G d , e , x f ,
g ,
[0062] any permutation between the higher indices of 9 G d , , e ,
x f , , g
[0063] does not modify the value of 10 G d , , e , x f , , g :
[0064] for example, for q=2, 11 G d , e , x g , f = G d , e , x f ,
g ,
[0065] permutating all the higher indices with all the lower
indices of 12 G d , , e , x f , , g
[0066] has the effect of conjugating the value of 13 G d , , e , x
f , , g : G f , , g , x d , , e = ( G d , , e , x f , , g ) * .
[0067] Furthermore, the following notation is adopted: the quantity
14 [ r ] G d , , e , x f , , g G h , , i , x j , , k G l , , m , x
n , , o
[0068] designates the linear combination of the r possible and
distinct products (modulo the three symmetries described here
above) of the type 15 G d , , e , x f , , g G h , , i , x j , , k G
l , , m , x n , , o ,
[0069] weighted by the value 1. Each of the r products is built
from the product 16 G d , , e , x f , , g G h , , i , x j , , k G l
, , m , x n , , o
[0070] in using and combining the following two rules of
permutation:
[0071] lower index of one of the terms of the product 17 G d , , e
, x f , , g G h , , i , x j , , k G l , , m , x n , , o
[0072] permutates with a lower index of another term of the same
product 18 G d , , e , x f , , g G h , , i , x j , , k G l , , m ,
x n , , o
[0073] to give another (distinct) product: for example, for
(q.sub.1, q.sub.2)=(2, 2), 19 G d , e , x f , g G h , i , x j ,
k
[0074] gives as other distinct products (modulo the three
symmetries described here above) 20 G h , e , x f , g G d , i , x j
, k , G d , i , x f , g G h , e , x j , k and G h , i , x f , g G d
, e , x j , k ,
[0075] a higher index of one of the terms of the product 21 G d , ,
e , x f , , g G h , , i , x j , , k G l , , m , x n , , o
[0076] permutates with a higher index of another term of the same
product 22 G d , , e , x f , , g G h , , i , x j , , k G l , , m ,
x n , , o
[0077] to give another (distinct) product: for example, for
(q.sub.1, q.sub.2)=(2, 2), 23 G d , e , x f , g G h , i , x j ,
k
[0078] gives the following as other distinct products (modulo the
three symmetries described here above) 24 G d , e , x j , g G h , i
, x f , k , G d , e , x f , k G h , i , x j , g and G d , e , x j ,
k G h , i , x f , g ,
[0079] in order to obtain the totality of the r possible and
distinct products (modulo the three symmetries described here
above).
[0080] The following example, where (q.sub.1, q.sub.2)=(2, 1) et
r=9, illustrates this notation: 25 [ 9 ] G d , e , x g , h G f , x
i = G d , e , x g , h G f , x i + G d , f , x g , h G e , x i + G f
, e , x g , h G d , x i + G d , e , x g , i G f , x h + G d , e , x
i , h G f , x g + G d , f , x g , i G e , x h + G d , f , x h , i G
e , x g + G f , e , x g , i G d , x h + G f , e , x i , h G d , x g
( 3 )
[0081] m Order Statistics
[0082] In the case of potentially non-centered, stationary or
cyclostationary sources, the m order circular statistics of the
vector x(t), given by (1), can be written: 26 C i 1 , i 2 , i q , x
i q + 1 , i q + 2 , i m = C i 1 , i 2 , i q , x i q + 1 , i q + 2 ,
i m ( t ) c ( 4 )
[0083] where the symbol 27 f ( t ) c = lim T .infin. ( 1 / T ) - T
/ 2 T / 2 f ( t ) t
[0084] corresponds to the operation of time-domain averaging, in
continuous time, of f(t) on an infinite horizon of observation. 28
C i 1 , i 2 , i q , x i q + 1 , i q + 2 , i m ( t ) = Cum ( x i 1 (
t ) , x i 2 ( t ) , , x i q ( t ) , x i q + 1 ( t ) * , x i q + 2 (
t ) * , , x i m ( t ) * ) ( 5 )
[0085] where q terms are conjugate and q terms are
non-conjugate.
[0086] The m order statistics described by the expression (5) are
said to be circular because the m order cumulant Cum{x.sub.d(t),
x.sub.e(t), . . . , x.sub.f(t), x.sub.g(t)*, x.sub.h(t)*, . . . ,
x.sub.i(t)*} is computed from as many conjugate terms (x.sub.g(t)*,
x.sub.h(t)*, . . . , x.sub.i(t)*) as non-conjugate terms
(x.sub.d(t), x.sub.e(t), . . . , x.sub.f(t)).
[0087] The m order circular cumulants may be expressed as a
function of the lower-than-m order moments as follows.
[0088] Let 29 M i 1 , i 2 , i q , x i q + 1 , i q + 2 , i m ( t )
and C i 1 , i 2 , i q , x i q + 1 , i q + 2 , i m ( t )
[0089] be the m order moments and circular cumulants associated
with the observation vector x(t), defined by: 30 M i 1 , i 2 , i q
, x i q + 1 , i q + 2 , i m ( t ) = E [ x i 1 ( t ) , x i 2 ( t ) ,
, x i q ( t ) , x i q + 1 ( t ) * , x i q + 2 ( t ) * , , x i m ( t
) * ] C i 1 , i 2 , i q , x i q + 1 , i q + 2 , i m ( t ) = k = 1 m
( - 1 ) k - 1 ( k - 1 ) ! g = 1 G k S g k Part g , x k M [ S g k ]
( t )
[0090]
Part.sub.g,x.sup.k={S.sub.g,x.sup.k(1).orgate.S.sub.g,x.sup.k(1).or-
gate. . . . .orgate.S.sub.g,x.sup.k(k)} designates the "g"th
partition, among "G.sub.k" possible partitions, of
<<k>> subsets 31 S g , x k = i r , i s , , i t , x i u
, i v , , i w ,
[0091] where 1.ltoreq.r.noteq.s.noteq. . . . .noteq.t.ltoreq.q et
q+1.ltoreq.u.noteq.v.noteq. . . . .noteq.w.ltoreq.m, such that 32
Part g , x k = i 1 , i 2 , i q , x i q + 1 , i q + 2 , i m .
[0092] As for the union operator, written as .orgate., it verifies
33 i t , i u , , i v , x i w , i y , , i z i j , i h , , i o , x i
p , i r , , i s = i j , i h , , i o , i t , i u , , i v , x i p , i
r , , i s , i w , i y , , i z .
[0093] In practical terms, two great classes of estimators may be
used to estimate the above statistics: in the case of stationary
sources, it is possible to use a non-skewed and consistent
empirical estimator for potentially centered, ergodic, stationary
sources whereas, in the case of potentially non-centered
cyclostationary sources, it is necessary to use what is called an
exhaustive, non-skewed and consistent estimator for potentially
non-centered, cycloergodic and cyclostationary sources. This
exhaustive estimator is, for example, determined according to an
approach described in the references [13-14].
[0094] Arrangement and Storage of m Order Statistics
[0095] As shown here above, the statistics 34 C d , e , , f , x g ,
h , , k ,
[0096] for 1.ltoreq.d, e, . . . , f, g, h, . . . k.ltoreq.N are
functions with m=2q inputs (m being an even-parity value), it is
then possible to arrange them in a matrix (N.sup.q.times.N.sup.q)
that will be named C.sub.m,x.
[0097] Explicitly, the quantity 35 C d , e , , f , x g , h , ,
k
[0098] is located at the ith row and at the jth column of the
matrix C.sub.m,x in writing i=N[ . . . N[N(d-1)+e-1]+ . . . ]+k t
j=N[ . . . N[N(g-1)+h-1]+ . . . ]+f.
[0099] Arrangement and Fourth-Order Statistics
[0100] In the case of centered stationary sources, the fourth-order
circular statistics of the vector x(t), given by (1), are written
as follows: 36 C d , e , x f , g = Cum { x d ( t ) , x e ( t ) , x
f ( t ) * , x g ( t ) * } = M d , e , x f , g - M d , e , x M x f ,
g - [ 2 ] M d , x f M e , x g , ( 4 )
[0101] In practical terms, these statistics may be estimated by
using a non-skewed and consistent empirical estimator for centered,
ergodic, stationary sources.
[0102] In the case of potentially non-centered cyclostationary
sources, the fourth-order circular statistics of the vector x(t) to
be taken into account are written as follows: 37 C d , e , x f , g
= Cum { x d ( t ) , x e ( t ) , x f ( t ) * , x g ( t ) * } C = M d
, e , x f , g ( t ) C - [ 4 ] M d , x ( t ) M e , x f , g ( t ) C -
M d , e , x ( t ) M x f , g ( t ) C - [ 2 ] M d , x f ( t ) M e , x
g ( t ) C + 2 M d , x ( t ) M e , x ( t ) M x f , g ( t ) C + 2 M x
f ( t ) M x g ( t ) M d , e , x ( t ) C + 2 [ 4 ] M d , x ( t ) M x
f ( t ) M e , x g ( t ) C - 6 M d , x ( t ) M e , x ( t ) M x f ( t
) M x g ( t ) C , ( 5 )
[0103] In practical terms, these fourth-order statistics may be
estimated by using the estimator known as the exhaustive,
non-skewed and consistent estimator for cyclostationary,
cycloergodic and potentially non-centered sources. This exhaustive
estimator is described in 13-14].
[0104] As shown here above, the statistics 38 C d , e , x f , g
,
[0105] for 1.ltoreq.d, e, f, g.ltoreq.N are four-input functions.
It is possible then to arrange them in a (N.sup.2.times.N.sup.2)
matrix that will be called a quadricovariance matrix Q.sub.x given
for example in the reference [13-14].
[0106] Arrangement and Sixth-Order Statistics
[0107] In the case of centered stationary sources, the sixth-order
circular statistics of the vector x(t), given by (1), are written
as follows: 39 C d , e , f , x g , h , i = Cum { x d ( t ) , x e (
t ) , x f ( t ) , x g ( t ) * , x h ( t ) * , x i ( t ) * } = M d ,
e , f , x g , h , i - [ 3 ] M d , e , f , x g M x h , i - [ 9 ] M d
, e , x g , h M f , x i - [ 3 ] M d , e , x M f , x g , h , i + 2 [
9 ] M d , e , x M f , x g M x h , i + 2 [ 6 ] M d , x g M e , x h M
f , x i ( 6 )
[0108] In practical terms, these statistics may be estimated by
using a non-skewed and consistent empirical estimator for centered,
ergodic, stationary sources.
[0109] In the case of centered cyclostationary sources, the
sixth-order circular statistics of the vector x(t) given by (1),
are written as follows: 40 C d , e , f , x g , h , i = Cum { x d (
t ) , x e ( t ) , x f ( t ) , x g ( t ) * , x h ( t ) * , x i ( t )
* } C = M d , e , f , x g , h , i ( t ) C - [ 3 ] M d , e , f , x g
( t ) M x h , i ( t ) C - [ 9 ] M d , e , x g , h ( t ) M f , x i (
t ) C - [ 3 ] M d , e , x ( t ) M f , x g , h , i ( t ) C + 2 [ 9 ]
M d , e , x ( t ) M f , x g ( t ) M x h , i ( t ) C + 2 [ 6 ] M d ,
x g ( t ) M e , x h ( t ) M f , x i ( t ) C ( 7 )
[0110] In practical terms, these sixth-order statistics may be
estimated by using an estimator called an exhaustive, non-skewed
and consistent estimator for cyclostationary, cycloergodic and
centered sources.
[0111] As shown here above, the statistics 41 C d , e , f , x g , h
, k
[0112] for 1.ltoreq.d, e, f, g, h, k.ltoreq.N are six-input
functions. It is possible then to arrange them in a
(N.sup.3.times.N.sup.3) matrix that will be called a hexacovariance
matrix H.sub.x.
[0113] Principle Implemented in the Method According to the
Invention
[0114] The method according to the invention uses especially a
property of multilinearity of the cumulants and the Gaussian nature
of noise which take the form of the following matrix expression
C.sub.m,x=[A.sup.{circle over (X)}(q-1){circle over
(X)}A*]C.sub.m,s[A.sup.{circle over (X)}(q-1){circle over
(X)}A*].sup.H (10)
[0115] where C.sub.m,x and C.sub.m,s are the matrices of the m
order statistics defined earlier, having respective sizes
(N.sup.q.times.N.sup.q) and (P.sup.q.times.P.sup.q), and being
associated with the vectors x(t) and s(t) where A.sup.{circle over
(X)}(q-1) corresponds to an adopted notation defined thus: the
matrix B.sup.{circle over (X)}k designates the matrix B raised to
the power (in the sense of the Kronecker product) k, i.e. in taking
the Kronecker product 42 B k = B B B k times ,
[0116] in writing B.sup.{circle over (X)}0=1.
[0117] The Kronecker product may be recalled here: let A and B be
two matrices respectively sized (L.sub.A.times.C.sub.A) and
(L.sub.B.times.C.sub.B). The Kronecker product D=A{circle over
(X)}B is a matrix sized (L.sub.A L.sub.B.times.C.sub.A C.sub.B)
defined by D=(A.sub.ij B)1.ltoreq.i.ltoreq.LA,
1.ltoreq.j.ltoreq.CA.
[0118] Without departing from the framework of the invention, other
modes of expression associated with other modes of arrangement of
the cumulants may be used:
C.sub.m,x=[A.sup.{circle over (X)}q]C.sub.m,s,t [A.sup.{circle over
(X)}q].sup.H (10a)
[0119] where 1 is chosen arbitrarily between 1 and q and where
C.sub.m,s,t. is the matrix of the m=2q order statistics of s(t)
associated with the index 1 chosen. Each expression conditions the
number of sources potentially identifiable from a given array.
[0120] Here below in the description, the analysis is given in
using the expression of the relationship (10).
[0121] Inasmuch as the sources are independent, the matrix of the m
order statistics associated with the sources, C.sub.m,s, is a
diagonal matrix. However, it turns out to be not a full-rank
matrix. The method according to the invention considers a matrix
determined from a full-rank diagonal matrix of the autocumulants
and from the matrix representing the juxtaposition of the P column
vectors relative to the direction vectors of the sources:
C.sub.m,x=A.sub.q.zeta..sub.m,sA.sub.q.sup.H (11)
[0122] where 43 m , s = diag ( [ C 1 , 1 , , 1 , s 1 , 1 , , 1 , ,
C P , P , , P , s P , P , , P ] )
[0123] is the full-rank diagonal matrix of the m=2q order
autocumulants 44 C p , p , , p , s p , p , , p
[0124] from the P sources, sized (P.times.P), and where
A.sub.q=[a.sub.1.sup.{circle over (X)}(q-1){circle over
(X)}a.sub.1* . . . a.sub.p.sup.{circle over (X)}(q-1){circle over
(X)}a.sub.p*] sized (N.sup.q.times.P) and assumed to be of full
rank, represents the juxtaposition of the P column vectors
[a.sub.p.sup.{circle over (X)}(q-1){circle over (X)}A.sub.p*].
Furthermore, we assume that the matrix
A.sub.q-1=[a.sub.1.sup.{circle over (X)}(q-2){circle over
(X)}a.sub.1* . . . a.sub.p.sup.{circle over (X)}(q-2){circle over
(X)}a.sub.p*], sized (N.sup.(q-1).times.P), is also a full-rank
matrix.
[0125] The method according to the invention enables the
advantageous exploitation and extraction, for example, of the
totality of the information proper to the direction vectors a.sub.p
of the sources, redundantly contained in the matrix of the m=2q
order circular statistics of the observations vector x(t),
C.sub.m,x and more particularly in the matrix A.sub.q.
[0126] The method comprises, for example, the steps described here
below. The samples of the vector x(t) are assumed to have been
observed and the matrix C.sub.m,x is assumed to have been estimated
from these samples.
[0127] Step 1: Singular Value Decomposition of the Matrix
C.sub.m,x
[0128] This step computes the square root C.sub.m,x.sup.1/2 of the
full-rank matrix C.sub.m,x, for example through the eigenvalue
decomposition of the Hermitian matrix C.sub.m,x=E.sub.s L.sub.s
E.sub.s.sup.H where L.sub.s and E.sub.s are respectively the
diagonal matrix of the P greatest (in terms of absolute value) real
eigenvalues of C.sub.m,x and the matrix of the associated
orthonormal eigenvectors. This step shows the relationship existing
between C.sub.m,x.sup.1/2 and A.sub.q: 45 C m , x 1 / 2 = E s L s 1
/ 2 = A q m , s 1 / 2 V H = [ a 1 ( q - 1 ) a 1 * a p ( q - 1 ) a p
* * ] m , s 1 / 2 V H ( 12 )
[0129] where V is a unit matrix, sized (P.times.P), unique for
L.sub.s and E.sub.s as given matrices, and where
.vertline.L.sub.s.Arrow-up bold..sup.1/2, .zeta..sub.sn.sup.1/2 are
square roots respectively of .vertline.L.sub.s.vertline. and
.zeta..sub.m,s (.vertline...vertline. designates the absolute value
operator).
[0130] L.sub.s and E.sub.s are, for example, respectively the
diagonal matrix of the P greatest (in terms of absolute value) real
eigenvalues of C.sub.m,x and the matrix of the associated
orthonormal eigenvectors.
[0131] For a full-rank matrix A.sub.q, it is possible to ascertain
that the hypothesis (H4) is equivalent to assuming that the
diagonal elements of L.sub.s are non-zero and have the same
sign.
[0132] Step 2
[0133] This step consists of the extraction, from the matrix 46 C m
, x 1 / 2 = [ 1 T , , N T ] T ,
[0134] of the N matrix blocks .GAMMA..sub.n: each block
.GAMMA..sub.n sized (N.sup.(q-1).times.P) is constituted by the
N.sup.(q-1) successive rows of 47 C m , x 1 / 2
[0135] starting from the "N.sup.(q-1) (n-1)+1"th row.
[0136] Step 3
[0137] This step entails the building of the N(N-1) matrices
.THETA..sub.n1,n2 defined, for all
1.ltoreq.n.sub.1.noteq.n.sub.2.ltoreq.- N, by
.THETA..sub.n1,n2=.GAMMA..sub.n1.sup.# .GAMMA..sub.n2 where #
designates the pseudo-inversion operator.
[0138] In noting for all 1.ltoreq.n.ltoreq.N,
.PHI..sub.n=diag([A.sub.n1, A.sub.n2, . . . , A.sub.nP]) where
A.sub.ij is the component of A located on the ith row and jth
column, there is equality .GAMMA..sub.n= 48 A ( q - 1 ) n m , s 1 /
2 V H
[0139] for all 1.ltoreq.n.ltoreq.N, and the fact that the matrix V
jointly diagonalizes the N(N-1) matrices 49 n1 , n2 = n1 # n2 = V m
, s - 1 / 2 n1 - 1 n2 m , s 1 / 2 V H ,
[0140] which, it may be recalled, is sized (P.times.P).
[0141] Step 4
[0142] This step consists in determining the matrix V.sub.sol,
resolving the problem of the joint diagonalization of the N(N-1)
matrices .THETA..sub.n1,n2 for example in using a method of
diagonalization described in the reference [2]. The matrix 50 C m ,
x 1 / 2 V sol ,
[0143] where V.sub.sol=V T is a unitary matrix jointly
diagonalizing the matrices .THETA..sub.n1,n2 to the nearest unitary
trivial matrix T, is an estimate of the matrix A.sub.q to the
nearest trivial matrix.
[0144] Different methods, known to those skilled in the art, enable
the extraction from A.sub.q of an estimate A; {circumflex over ( )}
of the mixture matrix A.
[0145] Step 5
[0146] Step 5A
[0147] One procedure consists, for example, in taking the average
of the K=N.sup.(q-1) blocks (.SIGMA..sub.k)* sized (N.times.P) (for
all 1.ltoreq.k.ltoreq.N.sup.(q-1) (the block .SIGMA..sub.k is
constituted by N successive rows of A.sub.q=[.SIGMA..sub.1.sup.T, .
. . , .SIGMA..sub.k.sup.T].sup.T starting from the "N (k-1)+1"th
row), or else in starting only one, for example (.SIGMA..sub.1)*.
This approach enables the estimation, in any order and excepting an
amplitude, of the P direction vectors a.sub.p and therefore the
mixture A matrix to the nearest trivial matrix.
[0148] Step 5B
[0149] Another step consists, for example, for each of the P
columns b.sub.p of 51 A q = [ a 1 ( q - 1 ) a 1 * a p ( q - 1 ) a p
* ]
[0150] in,
[0151] extracting the K=N.sup.(q-2) vectors b.sub.p(k) stacked one
beneath the other such that: 52 b p = [ a p ( q - 1 ) a p * ] = [ b
p ( 1 ) T , b p ( 2 ) T , , b p ( K ) T ] T ( 14 )
[0152] then
[0153] converting said column vectors b.sub.p(k)=(A.sub.ip . . .
A.sub.jp) [a.sub.p{circle over (X)}a.sub.p*] sized
(N.sup.2.times.1) into a matrix B.sub.p(k)=(A.sub.ip . . .
A.sub.jp) [a.sub.p a.sub.p.sup.H] (where 1.ltoreq.i, j.ltoreq.N)
sized (N.times.N) and
[0154] jointly decomposing these K=N.sup.(q-2) matrices into
singular values (singular value decomposition or SVD): the
eigenvector common to the K matrices B.sub.p(k) and associated with
the eigenvalue that is the greatest (in terms of modulus) is
therefore a column vector of the matrix A. It must be noted that
the quantity (A.sub.ip . . . A.sub.jp) is in the present case the
product of (q-2) components of A.
[0155] This step of processing on the P columns b.sub.p of A.sub.q,
enables the estimation, in any order and plus or minus one phase,
of the P direction vectors a.sub.p and therefore, to the nearest
trivial matrix, the mixture matrix A.
[0156] Step 6
[0157] The mixture matrix A representing the direction vector of
the sources contains, by itself, the information needed for the
angular localization of the sources. In this context, from the
estimation of the different columns of A, it is possible to
implement an arbitrary method of goniometry exploiting this
information. Such a method is presented for example in the document
[18].
[0158] Step 7
[0159] To estimate the sources vector s(t) in an over-determined
context (i.e. when P.ltoreq.N), the method applies an LIT type
filter to the observations x(t) explicitly using the estimation of
the mixture matrix A. It is possible, for example, to choose the
FAS filter described in the reference [17], which is optimal in the
presence of decorrelated sources.
[0160] The method comprises, for example, a step 0 which consists
of the building, from the different observation vectors, x(t), of
an estimate .sub.m,x of the matrix of statistics C.sub.m,x of the
observations, according to the method given earlier. In this case,
the steps 1 to 6 of the method are implemented on the estimate
.sub.m,x of the matrix.
[0161] Criterion of Performance
[0162] According to one alternative embodiment, the method
comprises a step using a normal criterion of performance for the
evaluation of the blind identification of mixtures. This criterion
is not global and enables the evaluation of the quality of
identification of each direction vector estimated: it is then
possible to compare two distinct methods of blind identification
with respect to each direction vector, and hence to each source.
This criterion is the extension, to the blind identification of
mixtures, of the criterion based on SINR
(signal-to-interference-plus-- noise ratio) given in the reference
[17] introduced for the blind extraction of sources. It is a
P-uplet described by
D(A, A;{circumflex over ( )})=(.alpha..sub.1, .alpha..sub.2, . . .
.alpha..sub.P) (15)
[0163] where 53 p = min 1 i P [ d ( a p , a ^ i , ) ] ( 16 )
[0164] and where d(u,v) is the pseudo-distance between the vectors
u and v, defined by
d(u,
v)=1-.vertline.<(u,v)>.vertline..sup.2.parallel.u.parallel..sup-
.-2.parallel.v.parallel..sup.-2 (17)
[0165] It may be noted that < . , . > designates the scalar
product defined for two vectors of a same dimension.
[0166] The method described for the m=2q order application can be
applied especially to the fourth-order and sixth-order statistics,
for example according to the examples given here below.
[0167] Application of the Method for the Blind Separation of
Fourth-Order Sources
[0168] An alternative embodiment of the method known as ICAR
(Independent Component Analysis using Redundancies in the
quadricovariance) exploits the m=4 (q=2) order statistics,
corresponding to the matrix of circular quadricovariance of the
Cm,x, written as Q.sub.xj. This method enables the blind
identification of the instantaneous mixture A or the blind
extraction of the sources s(t) when N.gtoreq.P, in other words only
when the mixture is over-determined.
[0169] The model (1) is assumed to be verified along with the
fourth-order hypotheses H.sub.1-4.
[0170] The method called ICAR exploits especially the expression
(11) which, when written for fourth-order statistics, is expressed
as follows:
Q.sub.x=A.sub.2 .zeta..sub.4,sA.sub.2.sup.H (18)
[0171] where 54 4 , s = diag ( [ C 1 , 1 , s 1 , 1 , , C P , P , s
P , P ] )
[0172] is the full-rank matrix of the fourth-order autocumulants
C.sub.p,p,s.sup.p,p of the sources, sized (P.times.P), and where
A.sub.2=[a.sub.1{circle over (X)}a.sub.1* . . . a.sub.p{circle over
(X)}a.sub.p*] sized (N.sup.2.times.P) and assumed to be full-rank,
represents the juxtaposition of the P column vectors [a.sub.p
{circle over (X)}a.sub.p*].
[0173] Furthermore, assuming that the mixture A matrix sized
(N.times.P), is also a full-rank matrix.
[0174] The method performs, for example, the steps 0 to 5 described
in the case of the m=2q order application, in using the following
parameters: Cm,x=Q.sub.x and .zeta..sub.m,s=.zeta..sub.4,s
[0175] In this example of an implementation, the method may also
include a step 0 which consists of: the building, from different
observation vectors x(t), of an estimate {circumflex over (Q)} of
the matrix of quadricovariance Q.sub.x of the observations. The
steps 1 to 6 are then carried out on this estimated value.
[0176] Examples of Results Obtained by Applying the Method to
Fourth-Order Statistics
[0177] FIGS. 2, 3 and 4 show a graph in which the x-axis
corresponds to the number of samples L and the y-axis to the
performance, the results of the simulation presenting the
performance of the method according to the invention in a
fourth-order application presented here above of ICAR (in
implementing the step 5B), ICAR.sub.2 (in implementing the step 5A)
and methods for the blind separation of sources (COM1, COM2, JADE,
FastICA) known to those skilled in the art. The conditions of
simulation are the following:
[0178] It is assumed that signals from P=3 non-filtered sources,
namely one BPSK source and two QPSK sources, are received on a
circular array of N=5 receivers such that R/.lambda.=0.55 (with R
and .lambda. the radius of the array and the wavelength) and such
that the signal-to-noise ratio, SNR, is equal to 20 dB for each
source.
[0179] The noise is Gaussian and spatially non-correlated.
[0180] The three sources are baseband sources and their time symbol
is chosen to be equal to the sampling time.
[0181] The criterion used to obtain the best appreciation of the
results of extraction from the source p for a given method is the
maximum signal-to-interference-plus-noise ratio associated with the
source p, better known as SINRM.sub.p [17]. It may be compared with
the optimum SINRM.sub.p computed by using not the estimated mixture
matrix but, on the contrary, the exact mixture matrix as well as
the exact statistics of the observations. It is this comparison
that is presented in FIGS. 2-4.
[0182] More particularly, FIG. 2 represents the SINRM of the source
1 associated with ICAR, ICAR.sub.2 along with the most efficient
methods currently being used for the blind separation of sources
such as JADE, COM1, COM2 and FastICA.
[0183] FIG. 3 shows the SINRM of the source 2 for the same methods
(ICAR, ICAR.sub.2, JADE, COM1, COM2, FastICA and FIG. 4 shows those
of the source 3.
[0184] In each of the figures, it can be seen that the two methods
ICAR and ICAR.sub.2 perform very well and are slightly more
efficient than the other methods JADE, COM1, COM2 and FastICA. As
for the FastICA algorithm, its best performance relates to the
source 2 for which it converges completely from 550 samples
onwards.
[0185] As for the difference in results between ICAR and ICAR.sub.2
in this configuration it proves to be negligible as compared with
the difference in performance between the ICAR methods and the
JADE, COM1, COM2 and FastICA methods which, however, are very
good.
[0186] Application of the Method to Sixth-Order Statistics
[0187] According to another alternative embodiment, the method uses
sixth-order statistics. This variant known as BIRTH (Blind
Identification of mixtures of sources using Redundancies in the
daTa Hexacovariance matrix) enables the identification from an
array of N receivers, of the direction vectors of at most P=N
sources (in the case of an array with different receivers). It also
enables the extraction of at most P=N sources for which the
direction vector's are explicitly identified.
[0188] The BIRTH method uses certain sixth-order statistics, stored
in the hexacovariance matrix C6,x designated H.sub.x. Thus, this
alternative implementation fully exploits the information proper to
the instantaneous mixture A of the sources, contained in H.sub.x,
especially through an artful writing of H.sub.x relative to the
direction vectors of the sources, this being done by means of the
property of multi-linearity of the cumulants:
H.sub.x=A.sub.3 .zeta..sub.6,s A.sub.3.sup.H (20)
[0189] where 55 6 , s = diag ( [ C 1 , 1 , 1 , s 1 , 1 , 1 , , C P
, P , P , s P , P , P ] )
[0190] is the full-rank matrix of the sixth-order autocumulants
C.sub.p,p,p,s.sup.p,p,p of the sources, sized (P.times.P), and
where A.sub.3=[a.sub.1.sup.{circle over (X)}2{circle over
(X)}a.sub.1* . . . a.sub.p.sup.{circle over (X)}2{circle over
(X)}.sub.p*], sized (N.sup.3.times.P) and assumed to be a full-rank
matrix, represents the juxtaposition of the P column vectors
[a.sub.p.sup.{circle over (X)}2{circle over (X)}a.sub.p*]=.left
brkt-bot.a.sub.p{circle over (X)}a.sub.p{circle over
(X)}a.sub.p*.right brkt-bot.. Furthermore, we assume that the
matrix A.sub.2=[a.sub.1{circle over (X)}a.sub.1* . . .
a.sub.p{circle over (X)}a.sub.p*], sized (N.sup.2.times.P), is also
a full-rank matrix.
[0191] The sixth-order method comprises the steps 1 to 6 and the
step 0 described here above in using the following parameters:
C.sub.m,x=H.sub.x et m=3.
[0192] Simulations
[0193] FIGS. 5, 6, 7 and 8 show a graph in which the x-axis
corresponds to the number of samples L and the y-axis corresponds
to performance, namely the efficient operation of the method
according to the invention in a sixth-order application. The
following are the conditions for obtaining the curves:
[0194] Let it be assumed that P=3 statistically independent
sources, more particularly two QPSK sources and one BPSK source,
all three being unfiltered, are received on a linear array of N=2
receivers such that R/.lambda.=0.55 (with R and .lambda.
respectively being the radius of the array and the wavelength).
[0195] The three sources, assumed to be synchronized, have the same
signal-to-noise ratio, written as SNR and being equal to 20 dB for
each source with a symbol time that is four times the sampling
time.
[0196] The BPSK source is chosen in baseband while the two QPSK
sources have carriers respectively equal to half and one-third of
the sampling frequency.
[0197] The mixture matrix A is chosen so that the column vectors of
the matrix A.sub.3 are linearly independent. The noise for its part
is Gaussian and spatially non-correlated.
[0198] The instantaneous mixture of noisy sources is considered to
be an over-determined mixture because the number of sources is
greater than the number of receivers. The algorithms JADE, COM1,
COM2, S3C2 known to those skilled in the art and the method of the
invention in a sixth-order application are implemented for the
blind identification of the over-determined mixture A.
[0199] The performance criterion .alpha..sub.p, defined by the
equation (18), is computed on 200 operations, and this is done for
each source p (1.ltoreq.p.ltoreq.3): it will thus enable the
comparison of the five methods.
[0200] According to the above assumptions, FIG. 5 shows the
variations in the quantity .alpha..sub.3 resulting from the
algorithms JADE, COM1, COM2, S3C2 and the method according to the
invention, BIRTH, depending on the number of samples. The method
according to the invention, unlike the prior art methods, makes it
possible to identify the directional vector in question.
[0201] FIG. 6 gives a view, in the same context, of the variations
of the triplet D(A, A;{circumflex over ( )})=(.alpha..sub.1,
.alpha..sub.2, .alpha..sub.3), associated with the method according
to the invention in a sixth-order application as a function of the
number of samples. The three coefficients .alpha..sub.p rapidly
decrease to zero as and when the number of samples increases.
[0202] FIG. 5 shows the variations in the quantity .alpha..sub.3
resulting from the prior art methods JADE, COM1, COM2, S3C2 and the
method according to the invention depending, this time, on the
signal-to-noise ratio (SNR) of the source 3. The method BIRTH is
fully successful in identifying the direction vector of the source
3 even for a low value of SNR.
[0203] Finally, let it be assumed that the above P=3 sources are
received on a circular array of N=3 receivers such that
R/.lambda.=0.55.
[0204] FIG. 8 then shows the variations of the quantity
.alpha..sub.3 resulting from the algorithms JADE, COM1, COM2 and
BIRTH has a function of the number of samples: the BIRTH algorithm
works in an over-determined context, namely in the context where
the number of sources is smaller than the number of receivers, and
although sixth order cumulants of the observations must be
estimated, the speed of convergence of BIRTH is the same magnitude
is that of the methods referred to further above.
[0205] Bibliography:
[0206] [1] P. COMON, "Independent Component Analysis, a new
concept?" Signal Processing, Elsevier, vol. 36, No. 3, pp. 287-314,
April 1994.
[0207] [2] J. -F. CARDOSO and A. SOULOUMIAC, "Blind beamforming for
non-gaussian signals," IEE Proceedings-F, vol. 140, No. 6, pp.
362-370, December 1993.
[0208] [3] H. HYVARINEN and E. OJA, "A fast fixed-point algorithm
for independent component analysis," Neural Computation, vol. 9,
No. 7, pp. 1483-1492, 1997.
[0209] [4] E. BINGHAM and H. HYVARINEN "A fast fixed-point
algorithm for independent component analysis of complex valued
signals," Int. J. of Neural Systems, vol. 10, No. 1, pp. 1-8,
2000.
[0210] [5] P. COMON, "From source separation to blind equalization,
contrast-based approaches," in ICISP 01, Int. Conf. on Image and
Signal Processing, Agadir, Morocco, May 3-5 2001, pp. 20-32.
[0211] [6] N. THIRION and E. MOREAU, "New criteria for blind signal
separation," in IEEE Workshop on Statistical Signal and Array
Processing, Pennsylvania, US, August 2000, pp. 344-348.
[0212] [7] T. W. LEE, M. S. LEWICKI, M. GIROLAMI, and T. J.
SEJNOWSKI, "Blind source separation of more sources than mixtures
using overcomplete representations," IEEE Signal Processing
Letters, vol. 6, No. 4, pp. 87-90, April 1999.
[0213] [8] P. COMON and O. GRELLIER, "Non-linear inversion of
underdetermined mixtures," in ICA 99, IEEE Workshop on Indep. Comp.
Anal. and Signal Separation, Aussois, France, Jan. 11-15 1999, pp.
461-465.
[0214] [9] J. F. CARDOSO, "Super-symetric decomposition of the
fourth-order cumulant tensor. Blind identification of more sources
than sensors," in ICASSP 91, Toronto, Canada, May 1991, pp.
3109-3112.
[0215] [10] P. COMON, "Blind channel identification and extraction
of more sources than sensors," in SPIE Conference, San Diego, US,
Jul. 19-24 1998, pp. 2-13.
[0216] [11] L. DeLATHAUWER, P. COMON, and B. DeMOOR, "Ica
algorithms for 3 sources and 2 sensors," in Sixth Sig. Proc.
Workshop on Higher Order Statistics, Caesarea, Israel, Jun. 14-16
1999, pp. 116-117.
[0217] [12] A. TALEB, "An algorithm for the blind identification of
N independent signal with 2 sensors," in ISSPA 01, sixteenth
symposium on signal processing and its applications, Kuala-Lumpur,
Malaysia, Aug. 13-16 2001.
[0218] [13] A. FERREOL, P. CHEVALIER, and L. ALBERA, patent
application entitled "Procede de traitement d'antennes sur des
signaux cyclostationnaires potentiellement non centres," (Method of
antenna processing on potentially non-centered cyclostationary
signals) filed on behalf of THALES, National filing No. 02 05575,
filed on 3rd May 2002.
[0219] [14] A. FERREOL, P. CHEVALIER, and L. ALBERA, "Higher order
blind separation of non zero-mean cyclostationary sources," in
EUSIPCO 02, XI European Signal Processing Conference, Toulouse,
France, Sep. 3-6 2002, pp. 103-106.
[0220] [15] L. ALBERA and P. COMON, "Asymptotic performance of
contrast-based blind source separation algorithms," in SAM 02,
Second IEEE Sensor Array and Multichannel Signal Processing
Workshop, Rosslyn, US, Aug. 4-6 2002.
[0221] [16] S. M. SPOONER and W. A. GARDNER, "The cumulant theory
of cyclostationary time-series, Part. II: Development and
applications," IEEE Transactions on Signal Processing, vol. 42, No.
12, pp. 3409-3429, December 1994.
[0222] [17] P. CHEVALIER, "Optimal separation of independent
narrow-band sources: Concept and Performances" Signal Processing,
Elsevier, vol. 73, pp. 27-47
[0223] [18] P. CHEVALIER, G. BENOIT, A. FERREOL, "DF after Blind
Identification of the source steering vectors: the Blind-Maxcor and
Blind-MUSIC methods", Proc. EUSIPCO, Triestre (Italy), pp
2097-2100, September 1996.
* * * * *