Patent ID: 12218711

FIG.1shows a telecommunication system comprising a set of Nt transmitting antennas T1, . . . , TNtand a set of Nr receiving antennas R1, . . . , RNr.

Each of the transmitting antennas T1, . . . , TNttransmits electromagnetic signals (generally representing data to be transmitted, encoded by symbols) in an associated communication channel C, where these signals are received by the different receiving antennas R1, . . . , RNr.

In the example described herein, the set of transmitting antennas T1, . . . , TNtis a Uniform Linear Array (ULA). The transmission of the electromagnetic signals between the transmitting antennas T1, . . . , TNtand the receiving antennas R1, . . . , RNris here physically modeled by a flat wave propagated (in direct line) along a defined direction between the centroid of the transmitting antennas T1, . . . , TNtand the centroid of the receiving antennas R1, . . . , RNr.

This situation corresponds in particular to the case of the so-called MIMO (“Multiple-Input Multiple-Output”) systems that comprise a plurality of transmitting antennas T1, . . . , TNtand a plurality of receiving antennas R1, . . . , RNr.

In this context of the MIMO systems, the propagation of the signals in the communication channels can be characterized by an angle of departure θtx, an angle of arrival θrx, and a delay of propagation T between the plurality of transmitting antennas T1, . . . , TNtand the plurality of receiving antennas R1, . . . , RNr.

The present invention takes place in a context of multi-path communication channel (for example, here, P paths) between the transmitting antennas T1, . . . , TNtand the receiving antennas R1, . . . , RNr.

The received signals y are written as:
y=XtTh+nwith Xtthe vector of the transmitted signals, θtxthe angle of departure, n a vector characterizing a thermal noise and the notation T corresponding to the matrix transposition operator.

In this case, the communication channel C is modeled by the vector h expressed as:

h=∑p=1Pβp⁢et(tp)where etis the steering vector characteristic of the antenna array, βpthe gain of the electromagnetic signals and tpthe variable defined by tp=cos (θtx).

In an extended version of the MIMO model, the vector h modeling the propagation channel is expressed as:

h=∑p=1Pβp⁢er(rp)⊗et*(tp)⊗ef(τp)with et, erand efthe steering vectors associated with the directions t, r and T, respectively, the variables t=cos (θtx), r=cos (θrx), T the delay of propagation and ⊗ the Kronecker product.

The propagation channels are hence characterized using these three variables. The object of the present invention is to determine them. The following disclosure describes the determination of one of these variables, for example, the angle of departure θtx, by considering a modeling according to only one of the three variables. As an alternative, a two-variable or three-variable modeling can be used.

FIG.2functionally illustrates a channel estimation device1according to an exemplary embodiment of the invention.

This channel estimation device1comprises a control unit2. The control unit2also comprises a detection device5and an estimation device7.

The control unit2comprises a processor20and a memory22.

The detection device5and the estimation device7are formed by a set of functional modules. For example, the detection device5comprises a signal reception module and a signal detection module. The estimation device7comprises the detection device5and an estimation module.

Each of the different modules described is for example implemented by means of computer program instructions adapted to implement the module in question when these instructions are executed by the processor of the control unit2.

The memory of the control unit2is for example adapted to memorize for example pilot signals, here linear signals, used to test the communication channels.

FIG.3is a flowchart showing an example of channel estimation method that can be implemented in the context described hereinabove.

As shown inFIG.3, the channel estimation method starts with step E2, during which the transmitting antennas T1, . . . , TNttransmit electromagnetic signals Xtto the receiving antennas R1, . . . , RNr.

At step E4, the signals y are received by the receiving antennas R1, . . . , RNr. These signals are transmitted in the plurality of communication channels. These received signals y are given by the expression introduced hereinabove.

The channel estimation method then continues with step E6. During this step, the control unit2determines noisy values z. These noisy values z are function of the transmitted signals Xtand the received signals y. They are representative of the transmission of the signals through the communication channels. On the other hand, these noisy values z are calculated without the use of the pilot signals. They are expressed for example according to the formula:
z=XtTywith the notation T corresponding to the matrix transposition operator.

These noisy values z will serve as a base for determining the value t associated with the angle of departure θtx.

In order to estimate the communication channel, the control unit2initializes an index i at the value 0 (step E8). This index i denotes the current iteration. During this step, the control unit2also initializes the value of a so-called “residue” variable resi(step E8). This variable resiis here initialized from the noisy values determined at step E6: res0=z. Each iteration makes it possible to search for a value taken by this residue resi.

As described hereinabove, the propagation of the signals in the communication channels is characterized by the variable t associated with the angle of departure θtx, the variable r associated with the angle of arrival θrxand the propagation delay T.

The channel estimation method then comprises steps of estimating the values of these variables. In the following, the steps presented describe the determination of the value taken by one of these variables (here the variable t associated with the angle of departure θtx).

The different steps are applied in the same way for determining the values of the other variables characterizing the propagation of the signals in the communication channels.

To characterize the propagation of the signals in the communication channel, a signal must first be detected in the communication system.

The channel estimation method then comprises a detection method Det of a signal in the communication system.FIG.4is a flowchart showing an example of a detection method as it can be implemented within the channel estimation method according to the invention.

As shown inFIG.4, at step E10, the detection method comprises a step of dividing the range of values that can be taken by the value t into a plurality of sub-ranges, for example Ktsub-ranges j.

By definition, the value t associated with the angle of departure θx is between −1 and 1. At step E10, the range [−1, 1] is thus divided into a succession of sub-ranges in such a way as to cover the whole range [−1, 1]. The number Kt is for example predetermined, before execution of the detection method, for example as a function of a desired level of performance.

The following of the detection method then consists in testing the presence of a signal on each of the sub-ranges constituting the range of values that can be taken by the value t associated with the angle of departure θtx.

For that purpose, at step E12, the processor of the detection device determines, for each sub-range j, a so-called “correlator” function f(resi, j). This correlator depends on the noisy values z. More particularly, each correlator associated with the sub-range j in question depends on the current value of the residue resi(which is equal to the noisy values z during the first iteration, where i=0). It also depends on the steering vector etassociated with the direction corresponding to the center tjof the sub-range j.

In practice, this correlator can be interpreted as a spatial filter associated with the sub-range in question, making it possible to filter the signals in order for example to distinguish them from the noise on this sub-range j.

According to an embodiment, the correlator is defined by the following expression:
f(resi,j)=|etT(tj)·resi|2with tjthe center of the sub-range in question, etthe steering vector associated with the direction corresponding to the center tjof the sub-range in question, resithe current value of the residue.

According to another embodiment, on each sub-range, each correlator is defined as a sum of windowed correlator functions associated with the sub-range in question. In this case, it is defined by the following expression:

f⁡(r⁢e⁢si,j)=∑k=1Nt⁢λk|etT(tj)·(rtx,k⊙resi)|2(1)with etthe steering vector associated with the direction t, tjthe center of the sub-range in question, resithe current value of the residue, λkscalars and rtx,kvectors and the notation ⊙ an operator symbolizing the term-by-term product between the different elements of the vectors in question (also called Hadamard product).

The scalars λkand the vectors rtx,kcorrespond respectively to the eigenvalues and eigenvectors associated with the matrix R0defined by the expression:

[R0]k,l=sin⁢c⁡(2⁢π⁢Δt(l-k)⁢a→t⁢xλ)with Δtthe width of the sub-range in question, ∥atx∥ the distance between two transmitting antennas, λ the wavelength of the transmitted signals and the notation sinc corresponding to the sinc function defined by sinc(x)=sin(x)/x.

The eigenvalues λkand the eigenvectors rtx,kdepend only on the width Δtof the sub-range j in question. The eigenvalues λkdecrease towards 0.

The expression
|etT(tj)·(rtx,k⊙resi)|2is herein called “windowed correlator function”.

As shown inFIG.4, the detection method then continues with step E14.

Previously, it is possible to note here that, in the two embodiments described hereinabove, the correlators observe a maximum when the current residue resiis collinear to the steering vector etassociated with the direction corresponding to the center tjof the sub-range j in question.

During this step, the processor of the detection device identifies the sub-range corresponding to the maximum likelihood of the correlator determined at step E12. More precisely, the control unit2determines the maximum value reached by the correlators among the correlators determined for each sub-range j.

Once the maximum value of the correlators determined, the processor of the detection device compares this determined maximum value with a predetermined threshold. This predetermined threshold is function of a level of noise associated with the studied range. It is for example here a Gaussian noise distributed in all the directions.

According to the first embodiment described hereinabove, the processor of the detection device therefore here compares with the predetermined threshold the maximum value among each of the correlator values calculated according to the first embodiment introduced earlier.

As an alternative, according to the second embodiment, the processor of the detection device compares, with the predetermined threshold, the maximum sum of windowed correlator functions determined among the different correlator function sums determined for each sub-range j.

In practice here, during this step, the control unit2identifies the sub-range on which the power of the transmitted signal is the highest and such that this signal cannot be considered as noise.

If, at step E14, the maximum of likelihood of the correlator (either the maximum value among the correlators according to the first embodiment, or the maximum sum among the windowed correlator function sums determined according to the second embodiment) is higher than the predetermined threshold, the processor of the detection device identifies the sub-range containing the searched value t of the angle of departure (step E16). We hence consider here that a signal has been detected in the communication system.

As shown inFIG.3, the channel estimation method then continues with step E20. During this step, the control unit2receives from the detection device the information about the detection of a signal in the communication system.

In the case where a signal has been detected, the channel estimation method then comprises an estimation method Est (described hereinafter and shown inFIG.5) for estimating the value associated with the variable characterizing the propagation of the signals in the communication channel; more particularly here, the estimation method relates to the estimation of the value t associated with the angle of departure θtxrelating to the detected signal.

If, at step E14, the maximum of likelihood of the correlator is lower than the predetermined threshold, it is considered that no signal has been detected. At step E20, the control unit2then receives the information that no signal has been detected by the detection device. The channel estimation method then continues with step E22. This absence of detection here forms a stop condition.

When the stop condition has been obtained, the communication channel is characterized on the basis of the set of signals previously detected at the current index i. The vector h is then defined by the following expression:

h=∑p=1Pβp,i-1⁢et(τp)with βp,ithe gain estimation of the electromagnetic signals associated with the estimation Tpof the value t obtained during the previous iteration of the method (preceding the current index i). We hence have here P=i−1. Each parameter βp,ihas hence been obtained during the previous iterations.

As an alternative, another stop condition can be defined, for example by determining the norm of the residue and by identifying when the latter is lower than a predefined threshold.

FIG.5is a flowchart showing an example of estimation method as it can be implemented within the channel estimation method.

In the case where a signal has been detected in the communication system (and hence a sub-range has been identified), the estimation method Est starts at step E30. During this step, the method of the estimation device determines, from the sub-range identified at step E14of the detection method, a first estimation t0of the value t of the angle of departure θtx.

According to an embodiment, this first estimation t0corresponds to the center tjof the identified sub-range.

Another embodiment is based on the values of the windowed correlators determined on the identified sub-range. More particularly, the first estimation t0of the value t of the angle of departure θtxis based on a test function ftdetermined on each sub-range j and depending of the windowed correlators previously determined:

ft(r⁢e⁢si,j)=∑k=1Ntλk⁢❘"\[LeftBracketingBar]"(rtx,k⊙resi)T·et(tj)|2

The determination of the first estimation t0is then based on a comparison of the windowed correlators associated with the sub-range j identified at step E14, hence corresponding to the signal detected during the detection method Det.

FIG.6shows the evolution of this test function ftas a function of the different values possible for the value t of the angle of departure θtx.

The curve f corresponds to this test function, taking into account all the windowed correlators. The curve a corresponds to the contribution of the first windowed correlator to the function ft. The curve b corresponds to the contribution of the second windowed correlator to the function ft.

This figure thus makes it possible to evaluate the contribution of each windowed correlator in the test function ft(resi, j). Studying each of these contributions of each windowed correlator thus makes it possible to evaluate the first estimation t0of the value t of the angle of departure θtx. In particular, the phase difference associated with each windowed correlator makes it possible to locate the portion of the sub-range that contains the searched value t (for example, on the left or the right of the sub-range center). The amplitude of each windowed correlator makes it possible to identify the position of the searched value t in the sub-range in question.

In practice, this first estimation t0of the value t associated with the angle of departure θtxcorresponds to an approximate estimation of this value.

At this step of the estimation method, we have an approximate estimation of the value t associated with the angle of departure θtx(we also know, thanks to the implementation of the method of detection of the sub-range in which is located this value t).

The following steps of the estimation method therefore have for object to refine this first estimation t0.

A conventional solution consists in using optimizing methods such as the Newton-Raphson method or the gradient descent. However, these methods are based on the use of functions having properties of convexity.

Here, the studied function (previously called correlator) defined by
f(resi,t)=|etT(t)·resi|2has a behavior similar to a Dirac function near its maximum (FIG.7). Conventional optimization methods are not suitable for processing this type of function. It is hence necessary in a first time to modify the correlator so as to then be able to apply thereto the conventional optimization methods.

The estimation method thus continues with step E32during which the processor of the estimation device determines a modified function fmod_i. This modified function fmod_iis determined on the basis of a scalar product between the vector associated with the noisy values and the steering vector. More precisely, here, the vector associated with the noisy values z here corresponds to the vector associated with the current value of the residue resi.

The modified function fmod_ihas properties of convexity, that is to say that this modified function fmod_iis either convex, or concave. These properties will hence allow implementing the conventional optimization methods.

Moreover, this modified function fmod_ireaches a maximum value for the value t corresponding to the maximum of the initial correlator (that is to say the non-modified shape).

In practice here, the determination of this modified function fmod_iamounts to determine the result of a convolution operation between the correlator |etT(t)·resi|2and a so-called “convolution kernel” function fn. This convolution kernel fnhere has for example the shape of a section of parabola (FIG.8).

The modified function fmod_i, to which will be applied the optimization method, here depends on a sum of windowed correlators (as introduced during the detection method):

fmo⁢d⁢_⁢i(resi,t)=∑k=1Ntλk′⁢❘"\[LeftBracketingBar]"etT(t)·(rtx,k′⊙resi)❘"\[RightBracketingBar]"2with λkscalars and rtx,kvectors depending of the chosen convolution kernel fn.

FIG.9represents the variation of this modified function fmod_ias a function of the different possible values of the variable t. The modified function fmod_ihere has the shape of a parabola centered on the value t associated with the searched angle of departure.

A conventional optimization method is hence applied to this modified function fmod_i. The Newton-Raphson method is here applied, decomposed into the following steps E34to E44. The object of this method is to determine the position of the maximum of the modified function fmod_i.

As shown inFIG.6, the estimation method Est continues with step E34of initializing an index l to the value 0. This index l denotes the current run of this optimization method. During this step, the control unit2also initializes the value of a variable tl. Here, the first estimation t0determined at step E30is used as the initialization value.

At step E36, the processor of the estimation device determines, for the current run, the values of the first and second derivatives of the modified function fmod_ias the value of the current variable tl. In other words, using the conventional notations, the processor of the estimation device determines the values fmod_i′ (z, tl) and fmod_i″(z, tl).

Then, at step E38, the processor of the estimation device determines the value of the variable tl+1defined by the following expression:

tl+1=tl-fmo⁢d′(z,tl)fm⁢od″(z,tl)

The method thus continues with step E40during which the processor of the estimation device evaluates if the determined value tl+1corresponds to the maximum of the modified function fmod_i. For that purpose, the processor calculates the quantity [tl+1−tl]. If this quantity is higher than a predetermined value ε(|tl+1−tl|>ε), it cannot be considered that the convergence towards the maximum is reached.

The method then continues with step E42during which the value tlis actualized by the value tl+1determined at step E38. The index l is also incremented. A new iteration is then implemented and the method restarts at step E36.

On the other hand, if, at step E40, the quantity |tl+1−tl| is lower than the predetermined value ε(|tl+1−tl|<ε), the determined value tl+1can be considered as representing the value of the maximum of the modified function fmod_i(step E44). In other words, this value tl+1corresponds to a finer estimation of the value t associated with the searched angle of departure θtx.

In the example of a parabolic convolution kernel, a single iteration is sufficient to reach the value of the maximum with the Newton-Raphson method.

As an alternative, other optimization methods can be used such as, for example, the gradient descent.

Here, the value associated with the angle of departure, obtained at this step, is denoted Tp.

Once obtained the estimation of the value Tpassociated with the angle of departure, the channel estimation method then continues with step E50(FIG.3).

From the obtained estimation of the value t associated with the direction of departure, the parameter βp,irepresenting the gain estimation of the electromagnetic signals associated with the estimation Tpof the value t obtained by the estimation derived from the current iteration i of the method is determined. More precisely, the determination thereof is based on the calculation of a pseudo-inverse. Indeed, using a matrix notation, the propagation channel h is written: h=Et·b with Etthe matrix containing the steering vectors etand b the matrix associated with the gain estimations βp,i.

The determination of the matrix b associated with the gain estimations βp(and hence the gain estimations βp,ithemselves) then uses a pseudo-inversion according to the following formula: b=(EtH·Et)−1·EtHh.

This parameter βp,iis here recalculated at each iteration of the method for all the values p lower than or equal to the current index i.

The signal detected is then fully characterized and the method is continued by actualizing the value of the residue resipreviously introduced (step E52) to take into account the last signal detected and the estimation associated with the corresponding angle of departure. In other words, the signal, among the signals remaining in the residue, whose power was the highest, is deduced from the noisy values z to obtain a new residue resi+1:

r⁢e⁢si+1=z-∑p=1Pβp,i⁢et(τp)

The new residue resi+1is completely recalculated from the noisy values z at each iteration because the gain estimations βp,iare actualized at each iteration.

The index i is then incremented at step E54and the method then restarts before the detection method Det, as long as the stop condition defined at step E20is not obtained.

These estimated values are used, in a context of demodulation, in order to remove the propagation channel influence in the data.

These estimated values can also be used by the control unit2to configure circuits for processing the electromagnetic signals received by the antennas R1, . . . , RNrof the array of antennas (these processing circuits being included in the control unit2but not shown so as to simplify the disclosure). These estimated values can also be estimated to configure pre-encoders adapted to perform a pre-encoding of the electromagnetic signals to be transmitted by means of the antennas R1, . . . , RNrof the array of antennas (when these antennas also operate in transmission as mentioned hereinabove).

Annex: Demonstration of Formula (1)

The present invention includes a step of detecting a ray whose parameters (here the variable t associated with the direction of departure θtx) are unknown, based on the samples received:
y=XtTh+nwithh=β·et*(t).

Then, the distribution of y knowing the direction of departure t is given by:

p⁡(y|t)=1(2⁢π⁢σn2)Nt⁢e⁢❘"\[LeftBracketingBar]"β❘"\[RightBracketingBar]"2⁢❘"\[LeftBracketingBar]"yH⁢XtT⁢et*(t)❘"\[RightBracketingBar]"2σn2

The direction of departure belongs to the range [−1, 1]. This range is divided into a succession of sub-ranges. The probability of presence of the ray in each sub-range is tested. The probability of presence of the ray in the range

[t¯-Δ⁢t2,t¯+Δ⁢t2]
is given by:

P⁡(t¯-Δ⁢t2<t<t¯+Δ⁢t2)=∫t¯-Δ⁢t2t¯+Δ⁢t2p⁡(y|t)⁢d⁢t

Generally, the following formula will be calculated

∫t¯-Δ⁢t2t¯+Δ⁢t2p⁡(y|t)⁢dt=∫-∞∞g⁡(t-t_)⁢p⁡(y|t)⁢dt⁢with⁢g⁡(t)={1⁢if⁢t∈[-Δ⁢t2,Δ⁢t2]0⁢else

This probability can be rewritten as:

P⁡(t¯-Δ⁢t2<t<t¯+Δ⁢t2)=1(2⁢π⁢σn2)Nt⁢∫-∞∞g⁡(t-t_)⁢e❘"\[LeftBracketingBar]"β❘"\[RightBracketingBar]"2⁢❘"\[LeftBracketingBar]"yH⁢XtT⁢et*(t)❘"\[RightBracketingBar]"2σn2⁢dt.

A first-order series expansion of the exponential near 0 is performed: ex≈1+x. The probability is now rewritten with introduction of constants K0and K1as:

P⁡(t_-Δ⁢t2<t<t¯+Δ⁢t2)≈K0+K1⁢∫-∞∞g⁡(t-t_)⁢❘"\[LeftBracketingBar]"yH⁢XtT⁢et*(t)❘"\[RightBracketingBar]"2⁢d⁢t

We denote:
z=XtHy

Then the integral can be rewritten as:

∫-∞∞f⁡(t-t_)⁢❘"\[LeftBracketingBar]"yH⁢XtT⁢et*(t)❘"\[RightBracketingBar]"2⁢dt=∫-∞∞g⁡(t-t_)⁢❘"\[LeftBracketingBar]"zH⁢et*(t)❘"\[RightBracketingBar]"2⁢dt=∑k,l=0Nt-1[z]k[z]l⁢∫-∞∞g⁡(t-t_)⁢e2⁢π⁢j⁢(l-k)⁢❘"\[LeftBracketingBar]"a❘"\[RightBracketingBar]"⁢tλ⁢dtwhere |a| represents the distance between two antennas and λ represents the wavelength associated with the frequency of the carrier used. It can be observed that the integral involves the Fourier transform of the kernel function g(t), called

Φg((l-k)⁢❘"\[LeftBracketingBar]"a❘"\[RightBracketingBar]"λ)

The ray presence probability can be written:

P⁡(t_-Δ⁢t2<t<t_+Δ⁢t2)≈K0+K1⁢zH⁢ΔH(t_)⁢R0⁢Δ⁡(t_)⁢z
where
Δ(t)=diag(et(t))
and

[R0]k,l=Φf((l-k)⁢❘"\[LeftBracketingBar]"a❘"\[RightBracketingBar]"λ)

In the case presented here, where

g⁡(t)={1⁢if⁢t∈[-Δ⁢t2,Δ⁢t2]0⁢else
we have:

[R0]k,l=Δ⁢t⁢sin⁢c⁡(2⁢πΔt(l-k)⁢a→txλ)

For each sub-range j, of size Δtand centered on tj, let's:
f(z,j)=zHΔH(tj)R0Δ(tj)zbe the part depending on z of the probability of presence of the ray in the range j.

The matrix R0can be decomposed into eigenvectors and eigenvalues as follows:

R0=∑k=0Nt-1λk⁢rtx,kH⁢rtx,kso that the function fi(z, j) can be written as

f⁡(z,j)=∑k=0Nt-1λk⁢❘"\[LeftBracketingBar]"etT(t¯j)·(rtx,k⊙z)|

During the phase of detection, the ray presence probability Pjfor each range j is calculated:

Pj=P⁡(tj¯-Δ⁢t2<t<tj¯+Δ⁢t2)

The sub-range the most liable to contain the ray is then sought. This operation amounts to retain the sub-range having shown the maximum probability. For that purpose, the complete calculation of Pjis not necessary, it is sufficient to calculate fi(z, j), the variable part of Pj. The value reached by f(z, j) in the sub-range retained is tested with respect to a threshold. The decision of the detection method is given by the test result.