Device and method for estimating interference and radiofrequency communication system
11108432 · 2021-08-31
Assignee
Inventors
Cpc classification
H04B17/336
ELECTRICITY
H04B1/1036
ELECTRICITY
H04B1/1027
ELECTRICITY
International classification
H04B1/10
ELECTRICITY
Abstract
A method comprises: Determining a set of all possible configurations of occupation or non-occupation of set of transmission bands, defined as a set of possible vectors satisfying at a time instant a non-overlapping condition of said radiofrequency system, said non-overlapping condition corresponding to the fact that only one interferer, among a set of possible interferers, can be active at a same time on each channel of said set of channels and forming, with contiguous channels, a transmission band, Obtaining measurements of occupation of at least a part of said set of channels, at respective tune instants, Performing probabilities calculations so as to determine, for each transmission band, an estimated activation rate, on the basis of said measurements, said estimated activation rate corresponding to an occupation rate of a transmission band by an interferer within said given observation time window.
Claims
1. A method implemented by a computer for estimating interference on a radiofrequency system using a set of channels (CHi), said interference being caused by interferers of an interfering system using a set of transmission bands (TBi), each of said transmission bands extending on a plurality of contiguous channels of said set of channels, Wherein the method comprises: Determining a set Ω of all possible configurations of occupation or non-occupation of said set of transmission bands, defined as a set of possible vectors Z.sub.k=[Z.sub.1,k, . . . , Z.sub.i,k, . . . Z.sub.I,k] satisfying at a time instant k a non-overlapping condition of said radiofrequency system, said non-overlapping condition corresponding to the fact that only one interferer i, among a set of I possible interferers, can be active at a same time k on each channel of said set of channels and forming, with contiguous channels, a transmission band i, Obtaining measurements X.sub.1, . . . , X.sub.k, . . . , X.sub.K of occupation of at least a part of said set of channels, at respective time instants k: 0<k≤K, where K defines a given observation time window, Performing probabilities calculations so as to determine, for each transmission band, an estimated activation rate i, on the basis of said measurements X.sub.1, . . . , X.sub.k, . . . , X.sub.K, said estimated activation rate i corresponding to an occupation rate of a transmission band i by an interferer within said given observation time window.
2. The method according to claim 1, wherein said radiofrequency system implements a frequency hopping on said channels, and said obtaining measurements X.sub.1, . . . , X.sub.k, . . . , X.sub.K is performed according to a frequency hopping implementation.
3. The method according to claim 1, wherein said probabilities calculations follow an expectation-maximum approach which iteratively approximates a maximum likelihood solution, according to two steps at each iteration t: a) For a fixed estimate τ.sup.(t), define the expectation Q(τ|τ.sup.(t)) such that
Q(τ|τ.sup.(t))=E.sub.Z∈Ω|X,τ.sub.
4. The method according to claim 3, wherein each transmission band extends over J channels and wherein said probabilities calculations comprise steps of 1) computing the probabilities, ∀k, ∀Z.sub.k∈Ω, P(X.sub.k|Z.sub.k), such as
5. The method according to claim 4, wherein an approximation is made that the coefficients Z.sub.i,k are independent according to the interferer index i, and the probabilities calculations of P(Z.sub.k|τ) are simplified into P(Z.sub.k|τ.sup.(t))=Π.sub.i P(Z.sub.i,k|τ.sub.i.sup.(t)) where P(Z.sub.i,k|τ.sub.i.sup.(t))=(1−τ.sub.i.sup.(t)))(1−Z.sub.i,k)=τ.sub.i.sup.(t).Math.Z.sub.i,k and wherein the iterative computations 4) of τ.sup.(t+1) are given by:
6. The method according to claim 4, wherein steps 2) to 4) are performed by: i. Defining a matrix A as follows: ∀Z ∈Ω, j ∈[1,|Ω|] the index of Z in Ω, and ∀i, A(i,j)=Z.sub.i Where |Ω| is the cardinality of Ω, ii. Decomposing A=U.sub.A[Δ.sub.A0.sub.I×|Ω|−I][V.sub.A.sup.TW.sub.A.sup.T].sup.T where V.sub.A is of size I×I and W.sub.A is of size (|Ω|−I)×I iii. Computing H=[V.sub.A.sup.T1.sub.I×1].sup.T and h(τ.sup.(t)=[(Δ.sub.A.sup.−1U.sub.A.sup.Tτ.sup.(t)).sup.T1].sup.T iv. Initializing a vector φ=H.sup.#h(τ.sup.(t)) of size equal to I v. Computing W.sub.H from a decomposition of:
H=U.sub.H[Δ.sub.H0.sub.I+1×|Ω|−I−1][V.sub.H.sup.TW.sub.H.sup.T].sup.T vi. Determining φ on the basis of the calculation of
W.sub.H.sup.T(W.sub.HW.sub.H.sup.T).sup.−1W.sub.H(φ−H.sup.#h(τ.sup.(t)))+H.sup.#h(τ.sup.(t)) vii. Getting P(Z.sub.k|τ.sup.(t))=φ viii. Obtaining θ.sup.(t+1) as:
7. The method according to claim 6, wherein step vi) is implemented by successive iterations for refining y according to two conditions considered alternatively from one iteration to the other: an orthogonal projection of φ on an hyperplane such as:
φ←W.sub.H.sup.T(W.sub.HW.sub.H.sup.T).sup.−1W.sub.H(φ−H.sup.#h(τ))+H.sup.#h(τ) and a projection of φ in an hypercube such as:
φ←max(min(φ,1),0).
8. The method according to claim 1, wherein a number of said contiguous channels forming a transmission band is four, a total number of channels of said set of channels being sixteen.
9. The method according to claim 8, wherein each of said channels extends over 5 MHz, whereas each of said transmission bands extends over 20 MHz with a spread spectrum technology implementation.
10. The method according to claim 1, wherein said radiofrequency system implements an ISM type communication, while the interfering system implements a Wifi type communication.
11. The method according to claim 1, wherein said non-overlapping condition derives from a CSMA/CA multiple access implementation performed by said radiofrequency system, said CSMA/CA multiple access defining communication timeslots, and said measurements X.sub.k being collected at each timeslot k.
12. The method according to claim 1, comprising further a selection for communication of at least one channel among said set of channels, said selected channel being within a transmission band for which said estimated activation rate i is the lowest.
13. Computer program comprising instructions for performing the method according to claim 1, when these instructions are run by a processor.
14. Device for estimating interference on a radiofrequency system using the set of channels (CHi), said interference being caused by said set of possible I interferers of an interfering system using the set of transmission bands (TBi), each of said transmission bands extending on a plurality of contiguous channels of said set of channels, Said device comprising a processing circuit for performing the method according to claim 1.
15. A radiofrequency communication system, comprising the device according to claim 14 for estimating interference susceptible to occur on the set of channels to be used by the radiofrequency communication system.
Description
BRIEF DESCRIPTION OF DRAWINGS
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
DESCRIPTION OF EMBODIMENTS
(10) A problem solved by the invention is the classification of interference (Wifi or not) according to measurements performed by a radio system. In the example disclosed hereafter, the radiofrequency system in stake has the following properties: Frequency hopping on channels of 5 MHz band at a rate around 4 ms Packets of 1.5 ms Within each slot of 4 ms, CSMA/CA or CSMA/CD multiple access (two transmission attempts) Measurements collected at each time slot
(11) The WiFi system properties are Spread spectrum technology on 20 MHz transmission bands WiFi PHY layer packet of 200 μs in general (including acknowledgements ACK) Transmission frame duration of around 100 ms in average (packet at transport layer) CSMA/CA multiple access
(12) It is therefore assumed in the present example of implementation that each current transmission band of the interfering system is to be considered with an observation on bands of 5 MHz among the sixteen channels of the radio system. Each current transmission band extends on a total band of 20 MHz according to the WiFi system properties, which corresponds here to four contiguous channels of the radio system. Of course, the numbers of 16 channels and 4 contiguous channels are examples given here and may admit variants. Also, the transmission bands (also called “W-channels” hereafter) can be overlapping, producing thus a total of 13 transmission bands. By indexing the channels which are 5 MHz wide from 1 to 16, the transmission bands are defined by the aggregations of channels with indexes: [1 2 3 4],[2 3 4 5], . . . ,[12 13 14 15], [13 14 15 16].
(13) The system is able to detect neighbouring interferers and try another connection in case of collision detection. It is worth to add that these interferers are not impacted generally by the current radio system which usually uses directional antennas at the transmitter. The interferers are not deemed to be active according to a frequency hopping scheme but rather to be active on a fixed subset of several channels as shown in
(14) each of the interferers is susceptible to be active on a transmission band formed by a plurality of contiguous channels of the whole set of channels, in particular at respective time instants k, and
(15) only one interferer can be active on each transmission band at a same time k,
(16) and for example when two transmission bands (having a width of 20 MHz) are overlapping on one or several channels (from one to three channels, each having a width of 5 MHz), an additional condition is that only one interferer can be active on each channel at a same time k.
(17) In order to distinguish hereafter the transmission bands of the interfering system (for example WiFi) and the channels of the communicating system (for example ISM), the transmission bands are called “W-channel” hereafter, while the “channels” remain the channels of the communicating system.
(18) An illustration of the coexistence of the radio system packets (uplink UL and downlink DL) and WiFi interference (WINT) is shown in
(19) Hereafter, the following notations are used: k is a time instant, i.e. a hop of the current radio system. More precisely, a time window of K time instants is considered. X.sub.k is an interference observation, over one channel in the present example of embodiment, at time instant k: 0<k≤K. Each X.sub.k can correspond to a power measurement or of a signal to noise ratio, sounded in each channel, or to a simulation according to a given scheme. f.sub.k is the index of the frequency channel sounded at time instant k. Z.sub.i,k is a random variable stating if interference is active on the i-th W-channel (out of I) at a time instant k.
(20) As a result of the CSMA/CA multiple access scheme, two interferers cannot overlap at the same time, which can be mathematically written as:
Z.sub.i,k=1.Math.∀j≠i/i−3≤j≤i+3,Z.sub.i,k=0 (1)
(21) The random variables are independent in time, i.e.,
∀i,∀(k,k′),k≠k′,E[Z.sub.i,kZ.sub.i,k′]=0 (2)
(22) Where E[ ] denotes the mathematical expectation. When this expectation must be performed on a finite set of values, the expectation is replaced and approximated by an arithmetic mean.
(23) It is considered at first a set Ω defined as a set of possible vectors Z.sub.k=[Z.sub.i,k . . . Z.sub.I,k] satisfying at a time instant k the non-overlapping condition.
(24)
(25) For example, when taking a channel index 9 having a D1 dot at the first third of the X-axis (one abscissa representing one possible configuration for a vector Z in the set 0), the W-channel 9 is composed of the channels 9,10,11,12. The D1 dot indicates that the channel 9 is the index used to identify this 20 MHz-wide W-channel. It can be observed that, because two W-channels cannot overlap at the same time, only the W-channels starting at index 13,5,4,3,2,1 can coexist at the same time of the W-channel starting at index 9. According to another example of
(26) It can be observed that 95 allocations with at least one interferer and satisfying the CSMA/CA non overlapping properties between interferers are possible (95 abscissa along the X-axis).
(27) One notation which is used also below also is: τ.sub.i, being the activation rate of the i-th interferer, i.e. τ.sub.i=E.sub.k[Z.sub.i,k]
(28) where E.sub.k[.] denotes the mathematical expectation over the different time instant. When this expectation must be performed on a finite set of values, i.e., when considering a finite time window, the expectation is replaced and approximated by an arithmetic mean.
(29) The problem to solve is the computation of the set τ from the vector of observations X=[X.sub.1, . . . , X.sub.K].
(30) The maximum a posteriori:
(31)
can be converted into a maximum likelihood problem by considering that all the variables τ are equiprobable a priori:
(32)
This involves that before receiving any observation, there is no information to say that one of the activation rates τ.sub.i for the i-th W-channel is higher than another.
Furthermore, the latent variables Z are also equiprobable a priori such that:
(33)
which gives a new optimization problem to solve.
(34) The problem (5) is intractable because of the high dimensionality of the set of possible vectors Z (∈Ω) and τ. Thus, an expectation-maximum approach which iteratively approximates the maximum likelihood solution is used instead, as proposed here after.
(35) That approach comprises two steps, at each iteration t: a) For a fixed estimate τ.sup.(t), define the expectation Q(τ|τ.sup.(t)) such that
Q(τ|τ.sup.(t))=E.sub.Z∈Ω|X,τ.sub.
(36)
(37) The invention proposes then an application of this iterative approach, by implementing computations which are specific to the problem raised above.
(38) Referring to
(39) in step S1, determining the set Ω of possible configurations,
(40) in step S2, initializing τ.sup.(1) with zero values and t at t=1; where t represents an iteration index in the iterative processing.
(41) in step S3, collecting measurements X.sub.1, . . . , X.sub.K,
(42) in step S4, computing the probabilities, ∀k, ∀Z.sub.k ∈Ω, P(X.sub.k|Z.sub.k), such as
P(X.sub.k|Z.sub.k)=p.sub.η(X.sub.k−1+Π.sub.j=0.sup.3(1−Z.sub.f.sub.
where p.sub.η( ) defines the probability relative to the estimation error on the observation X.sub.k.
(43) in step S5, computing the probabilities, ∀k, ∀Z.sub.k ∈Ω, P(Z.sub.k|τ.sub.(t)), with two possible embodiments as described below;
(44) in step S6, computing the values T.sub.Z.sub.
(45)
(46) in step S7, computing τ.sup.(t+1) from the previously computed T.sub.Z.sub.
(47) Steps S5, S6 and S7 are related to a given iteration index t and are repeated iteratively until a stopping condition is met, such as a maximum number of iterations. At the end of each new iteration, the iteration index t is increased.
(48) Two possible embodiments can be implemented for performing the calculation probabilities.
(49) In a first embodiment, the calculation of the probabilities in step S5 can simply be implemented as follows:
P(Z.sub.k|τ.sup.(t))=Π.sub.iP(Z.sub.i,k|τ.sub.i.sup.(t)) (8)
where
P(Z.sub.i,k|τ.sub.i.sup.(t))=(1−τ.sub.i.sup.(t))(1−Z.sub.i,k)+τ.sub.i.sup.(t).Math.Z.sub.i,k. (9)
and in step S7, the increment on t can be such that:
(50)
(51) Where the values T.sub.Z.sub.
(52) In a second embodiment rather based on a matrix computation, the following iterative procedure is applied, as detailed below referring to
(53) i. Define the matrix A as follows:
(54) ∀Z ∈Ω, j ∈[1, |Ω|] the index of Z in Ω and ∀i, A(i,j)=Z.sub.i, and |Ω| is the cardinality of Ω.
(55) (this matrix A visually corresponds then to
(56) ii. Decompose: (obtained with an “SVD algorithm”)
A=U.sub.A[Δ.sub.A0.sub.I×|Ω|−I][V.sub.A.sup.TW.sub.A.sup.T].sup.T
(57) where V.sub.A is of size I×I and W.sub.A is of size (|Ω|−I)×I.
(58) iii. Compute H=[V.sub.A.sup.T1.sub.I×1].sup.T and h(τ.sup.(t))=[(Δ.sub.A.sup.−1U.sub.A.sup.Tτ.sup.(t)).sup.T1].sup.T
(59) iv. Initialize a vector φ=H.sup.#h(τ.sup.(t)) of size equal to I.
(60) v. Compute W.sub.H from the decomposition obtained with an SVD algorithm
H=U.sub.H[Δ.sub.H0.sub.I+1×|Ω|−I−1][V.sub.H.sup.TW.sub.H.sup.T].sup.T
(61) vi. Then, repeat the following two steps until a maximum number of iterations is reached or until no modification of φ is obtained
φ←W.sub.H.sup.T(W.sub.HW.sub.H.sup.T).sup.−1W.sub.H(φ−H.sup.#h(τ.sup.(t)))+H.sup.#h(τ.sup.(t)) A.
φ←max(min(φ,1),0) B.
(62) vii. Get P(Z.sub.k|τ.sup.(t))=φ
(63) viii. And for performing step S7, θ.sup.(t+1) is obtained as
(64)
(65) A detailed description of the embodiments presented above is now given below.
(66) Starting from the expectation maximization procedure presented relatively to main steps a) and b) above, a development of each term according to the specific problem addressed by the invention is proposed hereafter.
(67) The proposed algorithm involves an iterative estimation of the variable τ. In a current iteration t, the variable τ is estimated by τ.sup.(t), and is refined into the output τ.sup.(t+1) of said iteration. In the following, T denotes a mathematical variable and is used in the following for deriving expressions useful for the computation of τ.sup.(t+1). For example, a general expression of P(Z.sub.k|τ) will be expressed, and the computation of P(Z.sub.k|τ.sup.(t)) can be obtained by replacing τ by τ.sup.(t) in the expression of P(Z.sub.k|τ). The random variables X.sub.k are mutually independent in time as well as the random variables Z.sub.k ∈Ω. Thus, it is possible first to decompose equation (6) of general step a) above into:
(68)
(69) Furthermore, P(Z.sub.k|X.sub.k,τ.sup.(t)) can be computed by using the Bayes theorem and the Law of total probability, as follows:
(70)
(71) By using the conditional probabilities, one can state further that
P(X.sub.k,Z.sub.k|τ)=P(X.sub.k|Z.sub.k,τ)P(Z.sub.k|τ) (16)
(72) The optimization problem can be thus stated as follows:
(73) In the considered system, when Z.sub.k ∈Ω is fixed, X.sub.k is not dependent on, which leads to
P(X.sub.k|Z.sub.k,τ)=P(X.sub.k|Z.sub.k) (17)
(74) More precisely, since X.sub.k is only potentially impacted by the interferers indexes [f.sub.k−3, f.sub.k] and their non-concurrent activation at time k, P (X.sub.k|Z.sub.k,τ) can be determined by:
P(X.sub.k|Z.sub.k,τ)=P(X.sub.k|Z.sub.k)=P(X.sub.k|Z.sub.f.sub.
(75) Hereafter, p.sub.η(.) is the notation for the measurement noise probability density function of the model X.sub.k=Y.sub.k+η.sub.k, where:
(76) Y.sub.k=0⇔∀0≤j≤3, Z.sub.f.sub.
(77) Y.sub.k=1 otherwise.
(78) Then:
Y.sub.k=1−Π.sub.j=0.sup.3(1−Z.sub.f.sub.
and in this case:
P(X.sub.k|Z.sub.k)=p.sub.η(X.sub.k−1+Π.sub.j=0.sup.3(1−Z.sub.f.sub.
(79) For example, in most cases, the power observation samples are impacted by a chi-square noise, which involves that p.sub.η(.) is defined as the probability density function, and parametrized according to the inherent thermal noise and number of signal samples used to compute a power observation. Another option is to consider that the noise effect is negligible, in this case p.sub.η(x)=1 if and only if x=0 and has a zero value otherwise. Furthermore, by using equations (12), (15), (16), (23), one can get to:
Q(τ|τ.sup.(t))=Σ.sub.kΣ.sub.Z.sub.
(80) Thus, the optimization equation (7) becomes
(81)
(82) Finally, as presented above with reference to
(83) Computing the probabilities, ∀k, ∀Z.sub.k ∈Ω, P(X.sub.k|Z.sub.k) (step S4)
(84) Computing the probabilities, ∀k, ∀Z.sub.k ∈Ω, P(Z.sub.k|τ.sup.(t)) (step S5)
(85) Computing the values T.sub.Z.sub.
(86)
(87) Computing τ.sup.(t+1) from the solution of equation (19) (step S7).
(88) Steps S5, S6, S7 are repeated iteratively until convergence. It is shown further here that convergence can be reached advantageously. The computation of P(X.sub.k|Z.sub.k), and P(Z.sub.k|τ) is detailed in the following. The computation of P (Z.sub.k|τ.sup.(t)) which is required for T.sub.Z.sub.
(89)
(90) More precisely, the Z.sub.i,k random variables are not independent since they are linked by the CSMA/CA contention protocol: they are mutually independent according to k but not to the interferer index i. However, it is made here the approximation that they are independent according to i, and then the coefficients Z.sub.i,k only depend on τ.sub.i. This approximation is fully justified in the case where few interferers are expected to be present, such as for example a train path in the countryside. In the case where several interferers are expected such as for an urban path, the second embodiment might be preferred.
(91) Thus, in the first embodiment, in step S11, the probability calculation P(Z.sub.k|τ) can be given as:
P(Z.sub.k|τ)=Π.sub.iP(Z.sub.i,k|τ.sub.i) (24)
where, by definition, the activation rate of the probability of having an interference (Z.sub.i,k=1) is τ.sub.i such that:
P(Z.sub.i,k|τ.sub.i)=(1−τ.sub.i)(1−Z.sub.i,k)+Σ.sub.i.Math.Z.sub.i,k (25)
(92) Then, in step S12, the expression of T.sub.Z.sub.
P(Z.sub.i,k|τ.sub.i.sup.(t))=(1−τ.sub.i.sup.(t))(1−Z.sub.i,k)+τ.sub.i.sup.(t).Math.Z.sub.i,k (26)
(93) Finally, in step S13, the expression of τ.sup.(t+1) can be calculated until a convergence condition is reached over t.
(94) More precisely, according to this approximation, the maximization problem (19) is concave. The zero of the derivative of the function to be optimized according to τ.sub.i, allows to get a relationship providing τ.sub.i.sup.(t+1) as
(95)
(96) Since Z.sub.i,k ∈{0,1}, this equation can be simplified into
(97)
which gives
(98)
(99) As a remark, in this case:
Σ.sub.Z.sub.
which leads to
(100)
(101) The calculations according to this first embodiment are simpler and necessitate less computation resources than the second embodiment, provided that the approximation made in step S10 can be justified.
(102) The second embodiment is detailed below with reference now to
τ.sub.i=E[Z.sub.i,k]=Σ.sub.Z.sub.
This property can be used then to evaluate P(Z.sub.k|τ) as detailed below.
(103) Noting hereafter the term φ.sub.j=p(Z.sub.k=A.sub.j|τ), the probability vectors φ and τ are linked by a matrix A having the dimensions I×|Ω| and defined in step S22 as follows:
∀Z∈Ω,j∈[1,|Ω|] is the index of Z in Ω, and ∀i,∀(i,j)=Z.sub.i (33)
A.sub.j denotes the j-th column of A, which is equal to the j-th element of Ω.
(104) Thus, φ is a vector of size |Ω| that characterizes P(Z.sub.k|τ). Furthermore, φ must be a probability vector. Thus, the set of solutions can be written as
(105)
(106) This system (34) is under-determined. It is proposed hereafter two algorithms to find a suitable solution satisfying the constraints.
(107) At first, in step S23, matrix A is decomposed as follows:
A=U.sub.A[Δ.sub.A0.sub.I×|Ω|−I][V.sub.A.sup.TW.sub.A.sup.T].sup.T (35)
(108) where V.sub.A is of size I×I and W.sub.A is of size (|Ω|−I)×I.
(109) Thus, the set of solutions φ satisfying τ=Aφ can be written
φ=A.sup.#τ+W.sub.A.sup.x (36)
where A.sup.# is the Moore-Penrose pseudo-inverse matrix of A, and x is any vector of size |Ω|−I. Thus, the solution φ belongs to a hyperplane having the equation
V.sub.Aφ=V.sub.AA.sup.#τ=Δ.sub.A.sup.−1U.sub.A.sup.Tτ (37)
(110) The constraint Σ.sub.j φ.sub.j=1 also puts φ onto a hyperplane. Thus, the intersection of the two hyperplanes can be written as
Hφ=h(τ) (38)
where H=[V.sub.A.sup.T1.sub.I×1].sup.T and h(τ)=[(Δ.sub.A.sup.−1U.sub.A.sup.Tτ).sup.T1].sup.T.
(111) Then, matrix H can be decomposed also as:
H=U.sub.H[Δ.sub.H0.sub.I+1×|Ω|−I−1][V.sub.H.sup.TW.sub.H.sup.T].sup.T (39)
(112) In step S24, the orthogonal projection of the point y on the hyperplane defined by Hφ=h(τ) is given by
φ=W.sub.H.sup.T(W.sub.HW.sub.H.sup.T).sup.−1W.sub.H(y−H.sup.#h(τ))+H.sup.#h(τ) (40)
(113) Moreover, in step S25, the vector φ must also satisfy the constraint ∀j, 0≤φ.sub.j≤1. Any vector y can be projected in the hypercube defined by the constraint by applying the per dimension max(.,.) and min(.,.) functions as max(min(φ,1),0).
(114) As it is difficult to both satisfy the constraints at the same time, it is proposed to alternate the calculations according to steps S24 and S25 between two operations in order to refine φ, namely: An orthogonal projection of φ on the hyperplane, according to step S24:
φ←W.sub.H.sup.T(W.sub.HW.sub.H.sup.T).sup.−1W.sub.H(φ−H.sup.#h(τ))+H.sup.#h(τ) (41) a projection of φ in the hypercube, according to step S25:
φ←max(min(φ,1),0) (42)
(115) Starting from an initialization of φ given by φ=H.sup.#h(τ), in step S26, the two operations are iterated several times. More precisely, several possible configurations can be tested until both conditions are met (projection on the hyperplane 41, and projection in the hypercube 42). It has been shown that this system converges to one solution actually. The steps S24 and S25 are repeated until a convergence condition is met which is for example given by a maximum number of iterations or a detection that the values φ has not changed in the last two iterations.
(116) The operation described above is used for computing φ by using h(τ.sup.(t)) instead of h(τ). Then, it is possible to extract directly P(Z.sub.k=A.sub.j|τ.sup.(t))=φ.sub.j and compute the T.sub.Z.sub.
(117) The expression of τ.sup.(t+1) can be obtained from the optimization problem (19), by first introducing A such that:
(118)
and then substituting the variable τ.sup.(t+1) by θ.sup.(t+1), and
P(Z.sub.k=Δ.sub.j|τ) with φ.sub.j, so that finally:
(119)
such that Σ.sub.j=1.sup.|Ω| φ.sub.j=1 and ∀j, 0≤φ.sub.j≤1
(120) Then, by using τ.sub.j=1.sup.|Ω| T.sub.Z.sub.
(121)
And finally, τ.sup.(t+1)=Δθ.sup.(t+1) is calculated in step S28.
(122)
(123) Of course, all the algorithms represented by flow charts of