Romain Couillet

Mérouane Debbah

Mathematical Foundations of Cognitive Radios

Recently, much interest has been directed towards software defined radios and embedded intelligence in telecommunication devices. However, no fundamental basis for cognitive radios has ever been proposed. In this paper, we introduce a fundamental vision of cognitive radios from a physical layer viewpoint. Specifically, our motivation in this work is to embed human-like intelligence in mobile wireless devices, following the three century-old work on Bayesian probability theory, the maximum entropy principle and minimal probability update. This allows us to partially answer such questions as "what are the signal detection capabilities of a wireless device?", "when facing a situation in which most parameters are missing, how to react?" and so on. As an introductory example, we will present previous works from the same authors following the cognitive framework, and especially the multi-antenna channel modelling and signal sensing.

Introduction

In 1948, Claude Shannon introduced a mathematical theory of communications [START_REF] Shannon | A mathematical theory of communications[END_REF], allowing two to three generations of research to design increasingly sophisticated telecommunication tools, whose purpose is to constantly increase the achievable transmission rate over various communication channels. One of the key conclusions of Shannon was to show that a linear increase in the transmission bandwidths is expected to provide linear growth in the channel transmission capacity, while linear transmit power increase only provide sub-linear capacity growth. As a consequence, the last decades of research in telecommunications led to a situation in which the available transmission bandwidth became dramatically scarce and can only be acquired by service providers at extraordinarily high prices. Then, in the end of the nineties, the conclusions of Foschini [START_REF] Foschini | On Limits of Wireless Communications in a Fading Environment when Using Multiple Antennas[END_REF] and Telatar [START_REF] Telatar | Capacity of multi-antenna Gaussian channels[END_REF] on their work on multiple antenna (MIMO) systems came as a salvation: when increasing the number of embedded antennas in both transmit and receive devices, a potential linear growth (with the number of antennas) in capacity was expected. Since the exploitation of the space dimension can come virtually at a zero cost compared to the exploitation of the frequency dimension, these stunning results rapidly generated lots of research work in the early years of the twentyfirst century. However, practical applications of multiple antenna systems took a long time to be put in place, when it was clearly realized that the exceptional predicted capacity gain could only come at a very strong signal to noise ratio (SNR) and for low correlated channels; for instance, line of sight components in a transmission almost completely annihilates the gain of multiple antenna systems. However, up to this point in the evolution of wireless devices, the initial result from Shannon was still applicable to the most advanced technologies. After the MIMO delusion, Joseph Mitola [START_REF] Mitola | Cognitive radio: making software radios more personal[END_REF] realized that a new virtual dimension could be exploited to increase the achievable transmission rate: making the radio smarter. The basic insight of Mitola was to observe that most allocated bandwidth is not efficiently used in the sense that, most of the time, large pieces of bandwidth are left unoccupied. Emabling the wireless devices to sense the frequency spectrum in a decentralized manner 1 allows for a potentially high increase of spectral efficiency, which we define here as the actual averaged transmission rate over the theoretical capacity. These large-scope ideas from Mitola recently motivated a wide range of research with common denominator the introduction of intelligence in wireless devices. For instance, Haykin [START_REF] Haykin | Cognitive Radio: Brain-Empowered Wireless Communications[END_REF] introduces the concept of interference temperature, which allows to control the level of interference allowed in a network, i.e. if a given user has a rate constraint largely inferior to the effective channel capacity, the excess unused rate could be used by another device, as long as this device does not request more than the available excess rate. This interference temperature brought the new idea of primary and secondary users in a wireless network: primary users are those subscribers who are charged a high price to communicate with high quality of service, while secondary users pay a lower price to communicate over opportunistic excess rates left unused by the primary users, e.g. [START_REF] Ghasemi | Collaborative spectrum sensing for opportunistic access in fading environments[END_REF][START_REF] Cardoso | Vandermonde frequency division multiplexing for cognitive radio[END_REF]. However, all these ideas, revolutionary as they may seem, only scratch the surface of a larger entity that is the cognitive radio. Indeed, if the cognitive radio is defined, as was supposedly the prior idea of Mitola or even more certainly the basic view of Haykin2 , as a radio in which all entities are capable of cognition, then the limitations in the capabilities of these radios is still unknown and not really explored. Concrete works on smart devices date back to Shannon's time as well. Claude Shannon was already interested in ideas such as a robot capable of playing chess [START_REF] Shannon | Programming a Computer for Playing Chess[END_REF]; he provided an original viewpoint of the cognitive abilities of future computers back in 1953 [START_REF] Shannon | Computers and Automata[END_REF] and even constructed a mind-reading machine, the circuitry of which is depicted in [START_REF] Shannon | A Mind-Reading Machine[END_REF]. In this work, we propose to define a fundamental basis for cognitive radios on a physical layer viewpoint, which enables human-like intelligence in wireless devices. This work comes as a rupture compared to previous telecommunication work, as we will no longer rely on Shannon's work, but will rather extend it. The reasons why we escape from Shannon's framework will be explained and justified in the following sections. The additional mathematical tools needed to extend Shannon's theory of information are the theory of Bayesian probabilities, the maximum entropy principle [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF] and the minimal cross-entropy principle [START_REF] Caticha | Maximum Entropy, fluctuations and priors[END_REF][START_REF] Caticha | Lectures on Probability, Entropy and Statistical Physics[END_REF], among others. The remainder of this paper unfolds as follows: in Section 2., we present the key philosophical ideas which lead from Shannon's classical information theory to Jaynes' more general probability theory. In Section 3., we provide two examples of direct application of Jaynes' maximum entropy principle to the problems of channel modelling and signal sensing. Then in Section 4., we discuss the present advantages and limitations of cognitive radios, and provide our conclusions in Section 5.. Notation: In the following, boldface lower-case symbols represent vectors, capital boldface characters denote matrices (I N is the N × N identity matrix). X i j denotes the (i, j) entry of X. The Hermitian transpose is denoted (•) H . The operators trX and |X| represent the trace and determinant, respectively. The symbol E[•] denotes expectation. The operator vec(•) turns a matrix X into a vector of the concatenated columns of X. Finally, the notation P x (y) denotes the probability density function of the variable x in position x = y.

From Shannon to Jaynes

We will first present a simple example to show the inherent limitations of Shannon's theory of information.

Channel Capacity Revisited

Let us consider the simplest communication scheme, modelled as

y = x + n (1)
for some transmit signal x, additive background noise n and receive signal y. The Shannon capacity C of such a system reads C = sup

p x I(x; y) (2)
with p x the probability distribution of the variable x taken in the set of single-variable probability distributions, and I denotes the mutual information [START_REF] Shannon | A mathematical theory of communications[END_REF]. The equality (2) can only be computed if the distribution of n is known. In practice, n is often taken as Gaussian, both for simplicity reasons and because this is somehow often close to the reality. However, there is no actual way to predict the distribution of the noise before transmitting data, and in reality the expression (2) is impossible to compute. This leads to the conclusion that all capacity computations are in fact only approximations of Equation ( 2). Moreover, it is important to observe that what we call noise is in effect the sum contribution of interfering waves with different properties. If part of this noise can be analyzed by the cognitive device 3 , then the capacity will increase. All these primary observations lead to realize that the channel capacity is largely dependent on the prior information available at the receiver. In particular, two identical receivers, facing the same channel, may have different actual capacities depending on the individual channel state information.

Assuming the noise is known to be Gaussian with zero mean, the receiver is left to estimate the noise variance. In general, only approximative values of the SNR are available. Therefore, the channel capacity might be better seen as a rate vector, with entries indexed by every possible values of the SNR and taking different degrees of probability. These degrees of probability differ for each receiver, making the capacity again information-dependent and user-dependent.

As a matter of fact, what one would call "real capacity", that would correspond to the capacity if the receiver knows exactly the noise variance, does not carry in itself any actual significance: as recalled by Jaynes [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF] pp. 634, the channel capacity is not an intrinsic value of the channel but an intrinsic value of the level of knowledge of the system designer4 .

Limitations of Information Theory

Already in 1963, Leon Brillouin [START_REF] Brillouin | Science and information theory[END_REF] it request to the transmitter in order to increase its performance is also an open question, e.g. should the receiver request more pilot symbols at the risk of a huge waste in spectral efficiency, should the receiver request some deterministic information regarding a given parameter of the channel? All these problems do not have deterministic channeldependent answers but depend on the specific knowledge of the transmitter/receiver pair to which some piece of additional information might or might not be valuable.

To partially answer those questions, we propose in the following to introduce first the notion of degrees of belief, which turns every deterministic measurable entity, e.g. the value of the channel capacity, the value of the SNR or the value of the channel fading, into a random variable with an assigned probability distribution: this probability distribution will translate the confidence of the cognitive devices regarding the estimation of the measurable entity in question. Then, we will introduce the notion of relevance which enables to estimate the relative importance of information.

Finally, we will discuss our general view of the capabilities of a cognitive radio.

The Bayesian Approach

As briefly stated in the previous section, we aim at extending the classical Shannon's information theory to enable cognitive devices with the ability of plausible reasoning. That is, a cognitive radio should not rely on empirical (often erroneous) decisions, but rather should be able to express doubt and to reason honestly when provided with limited knowledge. A first step in this approach is to turn empirical decisions into degrees of belief.

2.3..1 Degrees of Belief and the Maximum Entropy Principle

In the Bayesian philosophy, contrary to the orthodox probability philosophy, deterministic parameters of a system, e.g. a weight, a height, the channel delay spread, which a cognitive entity needs to evaluate, must be characterized by the degrees of belief attached to all possible values for this parameter. Therefore this gives a clear meaning for instance to the probability that the height of the Eiffel Tower is 50 m. As a consequence, assuming a cognitive telecommunication device is not aware of the intensity σ 2 of the background noise, instead of expressing the achievable transmission rate as the scalar C = log(1 + σ 2 ), which is therefore irrelevant to the communication device, it would be more adequate to consider the "vector" C(x) = log(1 + x), x ≥ 0, attached to a degrees of belief function, i.e. a probability density function, p(x) for each potential noise variance x. Two fundamental questions arise at this point: (i) how to use the vector C(x)?, and (ii) how to compute p(x)? Answering (i) is a matter of decision theory, in the sense that different requirements might come into play to decide on the actual transmission rate to use: if reliability is needed, one will decide to transmit at a rate log(1 + x) such that x 0 p(t)dt is less than a given (small) value, while if performance with low reliability is sought for, then x will take a larger value. This part of the cognitive radio spectrum will not be covered in this contribution. Question (ii), on the contrary, is the point of interest in the present paper. Given the total amount of prior information at the cognitive device, how to assign degrees of belief in a systematic way? The answer to this question partially appears in the work of Shannon [START_REF] Shannon | A mathematical theory of communications[END_REF] but is better explained and developed by Jaynes [START_REF] Jaynes | On the rationale of maximum-entropy methods[END_REF] thanks to the introduction of the maximum entropy principle (MaxEnt) [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. The key idea behind MaxEnt is to find a density function p, which fulfills the constraints imposed by the prior information I while introducing no additional (unwanted) information. In other words, this density function should maximize the ignorance about unknown parameters of the cognitive device, while satisfying the constraints given in I. In Jaynes' terms, this density function is maximally non-committal regarding missing information. This function translating ignorance is proven by Jaynes and more accurately later by Shore and Johnson [START_REF] Shore | Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy[END_REF] to be the entropy function H,

H(p) = -log(p(t))p(t)dt (3)
When the information contained in I is of statistical nature, such as first or second order statistics, the function p which maximizes the entropy while satisfying the constraints in I is unique and can be computed with Lagrangian multipliers. An example will be given in Section 3.1..

2.3..2 Relevance

The problem of relevance of information is a second topic in the establishment of foundations for cognitive radios. If cognitive devices were to act like human beings, they should be able to request additional information when they do not have enough evidence to take decisions. For instance, to obtain a more accurate estimate of the noise variance σ 2 in or- der to have more confidence on the achievable transmission rates, an intelligent device could require the transmitter to stop transmitting so that it can estimate σ 2 . But this would be an expensive waste in spectral efficiency, so it could alternatively request deterministic information on a dedicated channel from the transmitter. How accurate this information must be is then another problem. To be able to decide on what question to ask to the transmitter, the cognitive device needs to be able to judge the relevance of every possible question. This notion of questions, or inquiries, is a philosophical topic upon which little literature and very few concrete results exist. In 1978, Cox [START_REF] Cox | Of inference and inquiry, an essay in inductive logic[END_REF], who is also at the origin of the immense work from Jaynes on Bayesian probability theory [START_REF] Cox | Probability, frequency and reasonable expectation[END_REF], mathematically defined a question as the set of possible answers to this question. Therefore, a question will be relevant if its answers carry valuable information. Assuming the set of questions is seen as an ordered set, with the largest questions being the most relevant (since their answers carry potentially more cogent information), a cognitive device can decide which appropriate request to formulate to the transmitter. The work on relevance and questions is however still in its infancy, but we insist that those are fundamental needs to the cognitive radio field; for instance, interesting contributions are found in the works of Knuth [START_REF] Knuth | Bayesian Inference and Maximum Entropy Methods in Science and Engineering[END_REF], who uses lattice theory to create partial orders of finite sets of questions, which is seen as the dual (in the lattice theory terminology) of the set of answers to those questions.

2.3..3 What is a Cognitive Radio?

In our viewpoint, a cognitive radio must ideally be able to adapt to its environment, by gathering all cogent information about the propagation channel, the transmitted signal etc. while never producing undesirable empirical information. This would therefore relieve the telecommunication field from all ad-hoc methods, based on empirical decisions concerning unknown parameters. This does not mean that a cognitive device is not prone to making errors; however, these potential errors will never originate from erroneous system assumptions, but rather from lack of information, which would generate broad maximum entropy distributions 5 . If more cogent information is provided to a cognitive device, it will integrate it and increase its decision capabilities. In a way, the more signals a cognitive communication device is fed with, the more efficient it is; this would mean for instance that cognitive devices age wisely: the older the cognitive device, the more efficient.

Regarding for instance signal sensing, the first steps of which will be detailed in Section 3., we expect a cognitive device to process the received signals as follows,

1. initialization: integrate all cogent information about the communication channel, the properties of the supposedly received informative signal etc. and compute the degrees of belief associated to all relevant variables.

2. update loop: when the cognitive device is fed with incoming signals, it shall update its degrees of belief regarding all the previous variables and provide the overall probability that the received signal originates from a coherent data source.

3. decision: using some criterion from decision theory, e.g. the evidence for the presence/absence of a coherent data source is more than a given threshold, the cognitive device declares whether data originating from a coherent source have been received.

This protocol does not necessarily provide the most efficient sensing strategy in specific situations (sometimes it might provide a quick response, sometimes traditional algorithms might provide faster responses), but it provides the most honest way to treat the signal sensing problem. It is important to note that no signal detection strategy can be proven superior to any other as long as too much information on the communication environment is missing. If a given algorithm could be proven better than the Bayesian strategy, this would mean this algorithm has an information advantage; honesty would then require that the Bayesian strategy be aware of this additional piece of information. The significant advantage of the Bayesian philosophy and the maximum entropy principle over classical methods is that they do not to take any empirical guess to solve a problem. Therefore, instead of being either luckily very good, or unluckily very bad depending on the accuracy of this "guess", they perform as best as their prior information allows them to.

Also, a cognitive device ought to be capable of requesting information when it faces a situation where it crucially lacks cogent information; for instance, a cognitive mobile phone in a low network coverage situation, should be able to request information (or even help) to the neighboring cell phones which enjoy better coverage. The interest of this request would be measured by its relevance. Adding the possibilities of formulating inquiries might eventually lead to enabling cognitive devices with the ability of discussing, instead of just transmitting and receiving. Bidirectional communications used to be a point of deep interest when it was realized that Shannon's theory of communication is in fact precisely a theory of transmission, in which past transmitted information is assumed uncorrelated with subsequent transmitted information. In 1973, Marko proposed a generalization of Shannon's information theory framework to encompass bidirectional communications [START_REF] Marko | The Bidirectional Communication Theory -A Generalization of Information Theory[END_REF], in the objective to accurately model the social interactions among animals and especially human beings. The lead was then followed by Massey [START_REF] Massey | Causality, Feedback and Directed Information[END_REF] who extends information theory to include feedback in the expression of Shannon's mutual information.

Examples of Application

The most elementary requirement of a cognitive radio lies in its sensing capabilities. When a waveform is received at the cognitive device, it must be capable of deciding whether this waveform originates from a coherent source of information or if this waveform is pure background noise. When little is known by the receiver concerning the surrounding environment, this problem is very intricate and has led to lots of different ad-hoc techniques. Our purpose in the following is to provide a unique way of deciding on the presence of a coherent data source given a specific amount of prior information at the receiver. First, we will discuss channel modelling, which is a necessary step to properly handle the Bayesian signal detection method.

MaxEnt Channel Modelling

3.1..1 Introduction

Channel modelling is an entire field of research in telecommunications, which produces every year lots of new contributions. However, this huge amount of previous work on channel models leads to the following paradoxical conclusion: for a given total information gathered by a cognitive device, there exist many different channel models proposed in the literature. In such a situation, which of those channel models is the cognitive device supposed to trust? In reality, the fundamental difference between all those models lies in the additional hypothesis each of them, explicitly or implicitly, carries; some models might implicitly suggest that channels usually have a short delay spread for a given communication technology, or might suggest that it is very likely to have a strong line of sight component etc. However, if the receiver is not aware of that implicit information, this very information should honestly not be taken into account. What we will provide in the following is a systematic way to model channels, given some cogent information I, which fulfill the constraints imposed by I while being noncommittal regarding unknown parameters. In brief, we will provide the most elementary models compliant with I, without introducing unwanted hypothesis.

3.1..2 Gaussian i.i.d. Channels

Surprisingly enough, we will realize that most of the classical channels in the basic literature fall into the maximum entropy channel modelling methodology. This is the case of Gaussian i.i.d. channels. Indeed, let us assume that the information I known to the cognitive device gathers the following,

1. the transmitter is equipped with n T transmit antennas 2. the receiver is equipped with n R receive antennas

the channel carries an energy E.

The transmission model is

y = ρ n T Hx + n (4)
where x ∈ C n T is the transmitted symbol vector, n ∈ C n R the thermal or interfering noise, ρ the signal to noise ratio (SNR) and H ∈ C n T ×n R the channel we want to model. In mathematical terms, based on the fact that

dH n R ∑ i=1 n T ∑ j=1 |h i j | 2 P H (H) = n T n R E (Finite energy) (5) dP H (H) = 1 (P H (H) is a probability distribution) (6)
what distribution P H6 should the modeler assign to the channel? The modeler would like to derive the most general model complying with those constraints, in other words the one which maximizes our uncertainty while being consistent with the energy constraint. This statement is mathematically expressed by the maximization of the following expression involving Lagrange multipliers with respect to P H

L(P

H ) = -dHP H (H)logP H (H) + γ n R ∑ i=1 n T ∑ j=1 [E -dH|h i j | 2 P H (H)] + β 1 -dHP H (H) (7)
If we derive L(P H ) with respect to P H , we get

dL(P H ) dP H = -1 -logP H (H) -γ n r ∑ i=1 n t ∑ j=1 |h i j | 2 -β = 0 (8)
which yields

P H (H) = e -(β +γ ∑ n R i=1 ∑ n T j=1 |h i j | 2 = e -(β ) n R ∏ i=1 n T ∏ j=1 exp-(γ | h i j | 2 ) = n R ∏ i=1 n T ∏ j=1 P h i j (h i j ) with P h i j (x) = e -(γ|x| 2 + β +1 n R n T ) (9)
One of the most important conclusions of the maximum entropy principle is that, while we have only assumed the knowledge about the variance, this assumption naturally implies independent entries since the joint probability distribution P H simplifies into products of P h i j . Therefore, based on the previous state of knowledge, the only solution to the maximization of the entropy is the Gaussian i.i.d. channel. This does not mean that we have supposed independence of the channel fades in the model, nor does it mean that real channels ought to be i.i.d. if those are known to be of energy E. However, in the generalized L(P P ) expression, there exists no constraint on the dependence of the channel entries and this leads to natural independence as an honest guess on the behaviour of the channel entries. Another surprising result is that the distribution achieved is Gaussian. Once again, Gaussianity is not an assumption but a consequence of the fact that the channel has finite energy.

3.1..3 Other Channel Models

In [START_REF] Guillaud | Maximum Entropy MIMO Wireless Channel Models[END_REF], a more complete survey on MaxEnt channel models is proposed. We will gather in the following the main results.

If the information I at the receiver is the same as previously but the receiver is not aware of the exact value of the channel energy E but knows that it is contained in the interval [0, E max ], then

P H (H) = P H,E (H, E)dE (10) = P H|E (H)P E (E)dE (11)
If P E is assigned a uniform prior on the set [0, E max ], then we obtain 7 [START_REF] Jaynes | On the rationale of maximum-entropy methods[END_REF] Note that the distribution is invariant to unitary transformations, is not Gaussian and moreover the entries are not independent when the modeler has no knowledge on the amount of energy carried by the channel. This point is critical and shows the effect of the lack of information on the exact energy. If the channel covariance matrix Q = E(vec(H)vec(H) H ) is known to the receiver, and therefore is part of the side information I, then, denoting

P H (H) = 1 E max π n R n T ∞ 1 Emax u n R n T -2 e -∑ n R i=1 ∑ n T j=1 |h i j | 2 u du
Q = VΛ Λ ΛV H the spec- tral decomposition of Q, with V = v 1 , ..., v n R n T and Λ Λ Λ = diag(λ 1 , ..., λ n R n T ), P H (H) = 1 ∏ n R n T i=1 πλ i exp n R n T ∑ i=1 | v H i vec(H) | 2 λ i (13)

Signal Detection

Now that channel modelling has been investigated, the multiple antenna signal sensing problem can be completely handled.

3.2..1 Channel State Information

In this problem, the cogent information at the receiver is divided into known parameters, S-i) the receiver has n R antennas.

S-ii) the receiver samples as many as L times the input from the RF interface.

S-iii) the signal sent by the transmitter has a constant unit mean power. It is quite important to note that this hypothesis is very weak and should be made more accurate for communications schemes that are known only to use either QPSK, 16-QAM, 64-QAM modulations for instance.

S-iv) the MIMO channel has a constant mean power.

We similarly define additional information the receiver may be aware of V-i) the transmitter possesses (and uses) n T antennas.

V-ii) the noise variance σ 2 is known.

7 the assignment of uniform priors on variables defined on a continuous space is a very controversial point of the maximum entropy theory, which is longly discussed in [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF]. Another classically used prior, which solves the problem of invariance to variable change is the so-called Jeffreys' uninformative prior [START_REF] Jeffreys | An Invariant Form for the Prior Probability in Estimation Problems[END_REF].

3.2..2 Signal Model

Given a certain amount of sampled signals, the objective of the signal detection methods is to be able to optimally infer on the following hypothesis:

• H 0 . Only background noise is received.

• H 1 . Informative data added to background noise is received.

Given hypothesis S-iii), the only information on the transmitted signal (under H 1 ) is its unit variance. The maximum entropy principle claims that, under this limited state of knowledge, the transmitted data must be modelled as i.i.d.

Gaussian [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF]. The data vector, at time l ∈ {1, . . . , L}, is denoted

s (l) = (s (l) 1 , . . . , s (l) n T ) T ∈ C n T .
The data vectors are stacked into the receive matrix S = [s (1) , . . . , s (L) ]. If the noise level σ 2 is known, then either under H 0 or H 1 , the background noise must be represented, due to the same maximum entropy argument as before, by a complex standard Gaussian matrix Θ Θ Θ ∈ C n R ×L (i.e. a matrix with i.i.d. standard complex Gaussian entries θ i j ) [START_REF] Tulino | Random Matrix Theory and Wireless Communications[END_REF]. Under H 1 , the channel matrix is denoted H ∈ C n R ×n T with entry h i j being the link between the j th transmitting antenna and the i th receiving antenna. The model for H follows the Max-Ent channel modelling rules. In the present situation, only the constant mean power (or equivalently, the energy) of the channel is known. Therefore H will be modelled as i.i.d. Gaussian, following the reasoning in the previous section. The received data at sampling time l are given by the n T × 1 vector y (l) that we stack, over the L sampling periods, into the matrix Y = [y (1) , . . . , y (L) ] ∈ C n R ×L . This leads for H 0 to the model,

Y = σ Θ Θ Θ (14)
And for H 1 to

Y = [H , σ I N ] S Θ Θ Θ (15)
We also denote by Σ Σ Σ the covariance matrix

Σ Σ Σ = E[YY H ] (16)
= L HH H + σ 2 I n R (17) = U (LΛ Λ Λ) U H (18)
where Λ Λ Λ = diag ν 1 + σ 2 , . . . , ν n R + σ 2 , with {ν i , i ∈ {1, . . . , n R }} the eigenvalues of HH H and U a certain unitary matrix. Our intention is to make a decision on whether, given the received data matrix Y, the probability for H 1 is greater than the probability for H 0 . This problem is usually referred to as hypothesis testing [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF]. The decision criterion is based on the ratio

C(Y) = P H 1 |Y (Y) P H 0 |Y (Y) (19)
which we need to decide is whether lesser or greater than 1.

3.2..3 Results and Experiments

At this point in the derivation, computing C resorts to mere mathematical integration. The details of the calculus are given in [START_REF] Couillet | Bayesian Inference for Multiple Antenna Cognitive Receivers[END_REF]. We only provide here the results. First, assume σ 2 and n T are known, then, denoting x 1 , . . . , x n R the eigenvalues of YY H ,

P Y|H 0 (Y) = 1 (πσ 2 ) n R L e -1 σ 2 trYY H (20) = 1 (πσ 2 ) n R L e -1 σ 2 ∑ n R i=1 x i (21)
and

P Y|H 1 (Y) = Σ Σ Σ P Y|Σ Σ ΣH 1 (Y, Σ Σ Σ)P Σ Σ Σ (Σ Σ Σ)dΣ Σ Σ (22) = U(n R )×(R + ) n R P Y|Σ Σ Σ,H 1 (Y, UΛ Λ ΛU H )P Λ Λ Λ (Λ Λ Λ)dUdΛ Λ Λ (23)
which, after complete derivation, using in particular the Harish-Chandra identity [START_REF] Balantekin | Character expansions, Itzykson-Zuber integrals, and the QCD partition function[END_REF], gives (24) with P(k) the ensemble of permutations of k, sgn(b) the sign of the permutation b,

P Y|H 1 (Y) = α ∑ a⊂[1,n R ] e ∑ n T i=1 xa i σ 2 ∏ a i ∏ j =a 1 ... j =a i (x a i -x j ) × ∑ b∈P(n T ) (-1) sgn(b)+1 n T ∏ l=1 J n R -L-2+b l (n T σ 2 , n T x a i )
J k (x, y) = +∞ x
t k e -t-y t dt [START_REF] Massey | Causality, Feedback and Directed Information[END_REF] and

α = (n R -n T )!n T (2L-n T +1)n T /2 e n T 2 σ 2 - ∑ n R i=1 x i σ 2 n R !π n R L σ 2(n R -n T )(L-n T ) ∏ n T -1 j=1 j! (26)
These expressions are rather complex but show that the Bayesian signal detection, within the state of knowledge I, only depends on the eigenvalues x 1 , . . . , x n R of the Gram matrix YY H of the received data Y. A comparison with the classical power detector, e.g. [START_REF] Urkowitz | Energy detection of unknown deterministic signals[END_REF][START_REF] Kostylev | Energy detection of a signal with Random Amplitude[END_REF][START_REF] Quan | Spatial-spectral joint detection for wideband spectrum sensing in cognitive radio networks[END_REF], which consists in summing all individual powers received on the antenna array is provided in Figure 1. In the latter, n T = 1 and the comparison is made between the difference "correct detection rate minus false alarm rate" computed from Monte Carlo simulations for both Bayesian and classical signal detectors. We observe a slight gain in performance due to the novel Bayesian detector. Especially, for a low false alarm rate (which is often demanded in practice), we observe a large gain in correct detection rate provided by the Bayesian detector. This statement is however only valid for n T = 1. When T is larger, then the channel hardening effect reduces the gain of the Bayesian detector. This is shown in Figure 2 in which n T = 2. Now, if the noise power σ 2 is not perfectly known (this is classically the situation since knowledge of the noise power implies prior identification of the background noise), the probability distribution must be updated by marginalizing over σ 2 , from the lower bound σ 2 -to the upper bound σ 2 + on σ 2 . Therefore,

P Y|I = 1 σ 2 + -σ 2 - σ 2 + σ 2 - P Y|σ 2 ,I (Y, σ 2 )dσ 2 (27)
which is too involved to compute, but can be numerically estimated. An example is provided in Figure 3 in which the intervals [σ 2 -, σ 2 + ] are taken increasingly large. In the latter, correct detection rate against false alarm rate is depicted for different values of σ 2 -and σ 2 + . It is observed that the range of ensured correct detection gets increasingly narrower when [σ 2 -, σ 2 + ] is large. Note that this situation cannot be compared against classical power detection methods which do not provide solutions when σ 2 is not perfectly known.

Discussion

In addition to these first two studies on maximum entropy considerations for cognitive radios, the authors proposed more practical studies on maximum entropy OFDM channel estimation [START_REF] Couillet | A maximum entropy approach to OFDM channel estimation[END_REF], maximum entropy carrier frequency offset estimation [START_REF] Couillet | Information theoretic approach to synchronization: the OFDM carrier frequency offset example[END_REF], minimal update channel estimation [START_REF] Couillet | Minimal Update for OFDM Channel Estimation[END_REF] etc. From all those studies, we draw the following conclusions,

• quite often, classical techniques, in particular in the channel estimation field, are rediscovered using Max-Ent. However, it is important to note that, even if the final formulas are the same in the classical and Bayesian MaxEnt approach, the philosophical conclusions are very different. Usually classical methods derive from empirical parameter settings, which could have been chosen differently, while Bayesian approaches give unique deterministic solutions, which stem from honesty in the treatment of prior information.

• the MaxEnt principle allows one to marginalize over all parameters when those are not perfectly known. As a consequence, while classical solutions are found anew, those methods can usually be extended to cope with the lack of information on some key variables. For instance, in the signal sensing proposed in Section 3. and completed in [START_REF] Couillet | Bayesian Inference for Multiple Antenna Cognitive Receivers[END_REF], the situations where noise variance and number of transmit antennas are not perfectly known can be easily handled, whereas classical methods stumble on these problems and solve them by empirical (possibly largely erroneous) parametrizations.

On the other hand, MaxEnt calculus and final solutions can turn very rapidly extremely mathematically involved, as exemplified by the final signal sensing formula in Section 3.. This is a major problem, and the subject of most criticism towards Bayesian approaches. A missing part in these Max-Ent approach would be a systematic method which, from the general (very involved) solution, would provide approximate solutions. Quite remarkably, Caticha provides a vision of the maximum entropy principle, or more precisely the minimum cross entropy principle, which might help decide on the optimal approximation taken from a set of possible approximations [START_REF] Tseng | Maximum entropy approach to the theory of simple fluids[END_REF]. These considerations might lead to such systematic approximation methods. Another point of concern in the MaxEnt framework lies in the many integrals that may need to be computed when little is known on the surrounding environment. With the increasing capabilities of modern computers, numerical approximations might help to compute those integrals, but these approximations would only be valid if not so many integrals are considered; two reasons explain this fact: first, the complexity increase due to additional integrals is exponential in the number of integrals and second, small errors in inner integrals tend to lead to large errors when integrated many times (this is often referred to as the curse of dimensionality). As a consequence, while the first MaxEnt results provided by the authors show significant performance increase, many problems remain to be solved for cognitive radios to be fully intelligent, both on fundamental philosophical considerations (many questions raised in the introduction of the present paper are left unanswered) and on practical applications.

Conclusion

In this paper, we introduced the fundamentals of cognitive radios under a physical layer viewpoint. These fundamen-tals are based on the extension of Shannon's information theory to the Bayesian probability theory and the maximum entropy principle, which enables the cognitive devices with plausible (human-like) reasoning. Through the first-step studies of maximum entropy channel modelling and signal sensing, we paved the path to the establishment of strong theoretical grounds to the realm of cognitive radios.

Fig. 1 .Fig. 2 .RomainFig. 3 .

 123 Fig. 1. Detection amplitude comparison in SIMO -M = 1, N = 8, L = 20, SNR = -10 dB

  Let us go back to the channel capacity example above. If the receiver is provided with some additional information concerning the transmission medium, like the typical channel delay spread, the channel Doppler spread, the number of reflections, the presence of buildings in the neighborhood, how much does this affect the channel capacity is an open issue, which cannot be solved within Shannon's framework. And if the receiver experiences a poor decoding rate, what kind of information should

	realized the fundamen-
	tal limitation of Shannon's information theory. In his own
	words, 'The methods of [information] theory can be suc-
	cessfully applied to all technical problems concerning infor-
	mation: coding, telecommunication, mechanical computers,
	etc. In all of these problems we are actually processing in-
	formation or transmitting it from one place to another, and
	the present theory is extremely useful in setting up rules and
	stating exact limits for what can and cannot be done. But we
	are in no position to investigate the process of thought, and
	we cannot, for the moment, introduce into our theory any
	element involving the human value of the information. This
	elimination of the human element is a very serious limita-
	tion, but this is the price we have so far had to pay for being
	able to set up this body of scientific knowledge. The restric-
	tions that we have introduced enable us to give a quanti-
	tative definition of information and to treat information as
	a physically measurable quantity. This definition cannot
	distinguish between information of great importance and a
	piece of news of no great value for the person who receives
	it." -Leon Brillouin, 1963.
	Within the realm of cognitive devices, this situation in which
	information carries relevance, which depends on whom re-
	ceives it, typically arises.

so to limit the needs for control signaling.

remember that the title of his main contribution on cognitive radios[START_REF] Haykin | Cognitive Radio: Brain-Empowered Wireless Communications[END_REF] refers to "brain-empowered" radios

there is no reason why a cognitive device would not be able to infer on what the noise is made of.

the system designer can be seen as a virtual entity sharing the knowledge of both transmitter and receiver.

when little is known on a given parameter, the maximum entropy distribution attached to this parameter will be broad in the sense that no specific value is preferred to any other, while when more information is available on this parameter, the maximum entropy distribution will be very peaky around the exact value of the parameter.

it is important to note that we are concerned with P H|I where I represents the general background knowledge (here the variance) used to formulate the problem. However, for the sake of readability, P H|I will be denoted P H .

Acknowledgment

This work is partially supported by the European Commission in the framework of the FP7 Network of Excellence in Wireless Communications NEWCOM++.


B Fernández-Domínguez

B Pietras

N Patterson

J S Thomas

N Orr

M Chartier

W Catford

N L Achouri

J-C Angélique

N I Ashwood

A Banu

B Bastin

J Brown

R Borcea

S Franchoo

M Freer

L Gaudefroy

B Laurent

M Labiche

R C Lemmon

F Negoita

S Paschalis

E S Paul

M Petri

P Roussel-Chomaz

M Staniou

M J Taylor

L Trache

Structure around the island of inversion with single-neutron knockout reactions at GANIL

The nuclear structure of the 31 Mg nucleus has been studied with the singleneutron knockout reaction. We report on the preliminary results of an experiment performed with the EXOGAM array coupled, for the first time, to the SPEG spectrometer at GANIL. We present a provisional result for the inclusive single-neutron knockout cross section of σ inc = 90(5) mb. Preliminary exclusive cross sections for the measured bound states, including the ground state, are also presented. Finally, preliminary longitudinal momentum distributions for the ground state and first excited state are also shown. These results are compared to Monte Carlo Shell-Model calculations in the sd-pf region.

Introduction

The so-called "island of inversion" region of the nuclear chart, where transitions from normal to intruder ground state configurations across the N = 20 shell gap occur, was first placed around 20 ≤ N ≤ 22 and 10 ≤ Z ≤ 12. However, for the magnesium isotopes (Z = 12), the prediction of the neutron boundaries for this region is still not clear. Recent experimental efforts [1][2][3][4][5] have provided evidence to suggest that 29 Mg and 30 Mg are outside the island of inversion, and that 31 Mg is within the region. G. Neyens et al. [1] and M. Kowalska et al. [2] have determined the ground state spin, magnetic moment, and the low-lying structure of 31 Mg by hyperfine-structure and β -NMR measurements. Comparisons with advanced shell model calculations show that the experimental observables are correctly reproduced when allowing two neutron excitations across the N = 20 shell gap. The ground state wave function is, therefore, interpreted as dominated by a purely 2p-2h intruder configuration. A complementary β decay study on 31 Mg carried out by F. Maréchal et al. [3] finds weak feeding to the first excited states of 31 Al, which are of sd nature, supporting the interpretation of the previous references [1,2].

Single-neutron knockout reactions performed at MSU by R. Terry et al. [4] have allowed a mapping of the ground state configuration of 30 Mg. The intruder admixture of the ground state, measured by occupation numbers of neutrons in the f p single-particle states, was found to be small with a determined value of 0.60(12) neutrons. The structure of the low-lying states of the 30 Mg nucleus has been studied by means of β decay [5], Coulomb excitation [6] and β -delayed neutron decay [7]. The energies of the 0 + 1 and 2 + 1 , as well as the B(E2) values, are consistent with the sd shell model, confirming the 0 ω (or 0p-0h) structure of these states.

The study of the single-neutron knockout reaction of 31 Mg is a very interesting case since it is in the region where dramatic changes in single-particle structure are predicted. Indeed, 31 Mg is the only magnesium isotope between the normal, sd shell dominated configuration of 30 Mg, and the expected onset of the island of inversion at N = 20. The current measurement represents the "missing link" to the understanding of the neutron-deficient boundary in the Mg chain around the island of inversion. The technique of single-neutron knockout reactions is well suited to give quantitative insight into the single-particle components of the nuclear many-body wave function.

Experimental Set-Up

To study the structure of neutron-rich nuclei in the sdpf shell region with knockout reactions, the EX-OGAM array (a high efficiency Ge detector array), supplemented with NaI detectors, was coupled to the SPEG spectrometer for the first time at GANIL (France).

The primary beam of 77.5 MeV/u 36 S, delivered by the CSS1 and CSS2 cyclotrons with a beam power of 1 kW (∼ 10 13 particles per second) was focused on the 12 C SISSI target of 500 mg/cm 2 thickness, rotated by 45 degrees. The secondary beams were selected by the alpha spectrometer with an achromatic degrader at the intermediate focal plane. The central Bρ of the second dipole of the alpha spectrometer was set to 2.80 Tm, while the magnetic rigidity of the first dipole was tuned to optimise the production of the light neutron-rich isotopes, namely 23 O and 25 F. However, due to the relatively large momentum acceptance of the alpha spectrometer, elements from Oxygen to Magnesium were transmitted to the SPEG spectrometer.

The secondary cocktail beam was continuously monitored with a thin scintillator and a drift chamber. A secondary 12 C reaction target of thickness 170 mg/cm 2 was placed at the focal point in between the SPEG analyser and the SPEG dipole. The spectrometer was operated in an achromatic mode at the target position.

Isotopic identification of the ions after the reaction target was achieved with a set of detectors near the SPEG focal plane: two drift chambers, measuring the particle trajectories (X,Y); one ionisation chamber, measuring energy loss (∆E); and a plastic scintillator, measuring the residual energy (E). The nuclear charges of the elements were deduced from a ∆E vs E identification matrix. The isobaric identification was obtained from time-of-flight measurements between the plastic at the end of the focal plane of SPEG and the Radio Frequency from the cyclotron. Finally, to separate the projectile-like reaction products from the incoming beam particles, an additional ∆E vs focal plane position selection was applied.

The gamma rays following the de-excitation of the fragments emitted in an excited state were detected at the reaction target position with a hybrid gamma-ray detection array, consisting of 12 NaI and 8 Ge detectors. The array was arranged in a symmetric configuration: two rings of 4 highly segmented EXOGAM clovers each, were placed around 45 • forward and backward at distances 215 and 134 cm, respectively. In addition, 4 clusters of 3 NaI detectors were placed at 90 • at a 310 cm distance to the target. A view of the geometrical configuration is presented in Fig. 1. Total photopeak efficiency for the full EXOGAM-Ge array was 3.3% with an average energy resolution of 2.7% at 1.3 MeV.

Analysis and Preliminary Results

The inclusive knockout cross section of σ inc = 90(5) mb was obtained from the identified projectile-like fragments reaching the final SPEG focal plane, relative to the beam current of 31 Mg and the number of atoms in the target. The uncertainty in the cross section of 5% is mainly dominated by the normalisation procedure and software gates. A velocity measured in an event-by-event basis in SPEG (average β = 0.317) was used to correct for the Doppler shift. The final Doppler-corrected and add-back reconstructed energy spectrum obtained from the EXOGAM array in coincidence with 30 Mg knockout fragments is presented in Fig. 2. A chi-squared fitting procedure was used to minimise the difference between the data and a combined simulated Geant4 response of the entire array, together with an exponential function to describe the background. From singles spectra analysis, six gamma-ray transitions have been determined in the 12 C( 31 Mg, 30 Mg+γ)X reaction. Observed transitions at energies of 988(2), 1484(1), 1822(2), 1981(2), 3534(6) keV, are consistent with the results observed in previous measurements [5][6][7]. A new transition at 958(2) keV, which is well resolved in the doublet near 1 MeV, has been identified in the present work. Gamma-gamma coincidences have allowed the construction of a tentative level scheme presented in Fig. 3. The new transition has given rise to a new level in 30 Mg at 2442(2) keV, not previously observed. Exclusive cross sections to the final states in 30 Mg populated in the knockout reaction are summarised in Table 1. The intensity of each transition is calculated by normalising the number of detected gamma-rays, corrected by the effective efficiency, to the number of reaction products. The branch, b, for the direct population of each level was obtained by intensity balance. Only statistical uncertainties are shown. The uncertainties from the gamma-ray efficiency is expected to be lower than 4%. The ground state cross section was obtained by subtracting from the inclusive cross section all of the observed excited-state contributions. 6) 3534(6) 14.9(8) 14.9(8) 13.4 (10) Table 1: Results for the 12 C( 31 Mg, 30 Mg)X reaction. The proposed energy levels, and transition energies as well as the transition intensities I γ , direct population branches b, and the resulting cross sections σ are given.

E level (keV) E γ (keV) I γ (%) b (%) σ(
The longitudinal momentum distributions of the knockout fragment 30 Mg were measured in the SPEG spectrometer. The eikonal calculations describing the reaction mechanism have been folded with the experimental filter including the beam spread, energy loss and angular straggling, and the resolution of the spectrometer. This observable characterises the relative orbital angular momentum of the neutron removed from the projectile ground state wave function. The preliminary experimental results for removal to the ground state and first excited state of 30 Mg are compared to the theoretical calculations in Fig. 4 and Fig. 5, respectively. Due to the large acceptance of the SPEG spectrometer, acceptance corrections are less than 2-3%. The ground state distribution is in good agreement with ℓ = 0, confirming the spin and parity assignment of 1/2 + for 31 Mg found by G. Neyens et al. [1]. The first excited state presents a much broader distribution, suggesting an assignment of ℓ = 2, consistent with the assignment reported in [5].

Discussion

In Table 2 exclusive cross sections are compared to the single-particle cross sections , σ sp , calculated using the reaction theory described in references [8,9]. Deduced spectroscopic factors are compared to Monte Carlo shell model calculations using the SDPF-M interaction [10] provided by Y. Utsuno et al [START_REF] Utsuno | Private Communication[END_REF].

The results obtained from the SDPF-M calculations show a significant disagreement with the experimental values. The measured spectroscopic factor for the ground state is around 5 times larger than the shell model predictions. On the other hand, the first-excited state shows a 33% lower value than the predicted one. Although not all the experimental uncertainties have been applied, in total they are expected to be lower than 10%-still not significant enough to account for the differences between theory and results. There is still no clear explanation for the disagreement. For the ground state exclusive cross section (hence, deduced spectroscopic factor) to be in closer agreement with the theoretical prediction, a significant fraction of gamma transition strength would have to have gone unobserved. These presumably higher-lying, unobserved transitions would lead to a reduction in the proportion of the inclusive cross section attributed to the direct population of the ground state of 30 Mg. This seems unlikely.

Conclusions

Preliminary results on single-neutron knockout measurements from 31 Mg have been discussed. Longitudinal momentum distributions for the ground state and first excited state show a good agreement with ℓ = 0 and ℓ = 2 values, respectively. Preliminary spectroscopic factors are compared to SDPF-M calculations. The significant discrepancies between the theory predictions and the experimental results could lead to some revision of the picture of 31 Mg.

Fig. 1 :

 1 Fig. 1: Hybrid gamma detection array mounted around the target position in SPEG.

Fig. 2 :

 2 Fig. 2: Germanium γ-ray spectrum of 30 Mg. Fig. 3: Preliminary level scheme of 30 Mg

Fig. 3 :

 3 Fig. 2: Germanium γ-ray spectrum of 30 Mg. Fig. 3: Preliminary level scheme of 30 Mg

Fig. 4 :Fig. 5 :

 45 Fig. 4: Longitudinal momentum distribution for the ground state in 30 Mg

Table 2 :

 2 Theoretical and experimental spectroscopic factors E th c (MeV) J π nlj C 2 S th E exp c (MeV) l σ sp (mb) C 2 S exp

	g.s.	0 + 2s 1/2 0.10	g.s.	0 28.3	0.51(4)
	1.56	2 + 1d 3/2 0.92	1.48	2 13.7	0.62(6)


V Kerzérho

P Cauvet

S Bernard

F Azais

M Comte

M Renovell
email: renovell@lirmm.fr

A Multi-Converter DFT technique for Complex SIP: Concepts and Validation

This paper presents a new technique called "Analog Network of Converters" that allows to test a set of ADCs and DACs embedded in a complex circuit as SiP and SoC. It presents an experimental validation of this new concept that permits to reduce drastically the testing time and requires only a low cost digital ATE.

I. INTRODUCTION

The convergence of both communication and consumer markets increases silicon usage and needs for integrating various analogue or mixed-signal blocks into a single System-in-Package (SiP) or System-on-Chip (SoC): mobile phones, laptops, audio players… Figure 1 gives an example of such a complex mixed-signal system integrated into a single package. The integration of many different functions into a single package offers several clear benefits but, on the other hand, implies very significant test challenges. Analogue blocks are usually tested with respect to their specified parameters. Consequently, the main difficulty comes from the performance requirements on the test instruments. Indeed, analogue testing consists in a long sequence of parameter characterization that is performed using very expensive instruments able to accurately measure analogue signals. In addition to these required expensive instruments, we should note that controllability and observability of deeply embedded analogue blocks are limited and the possibility of external testing may be limited. Also, as signals become faster and systems are operated at higher speeds, external testing becomes more susceptible to noise, crosstalk and probing problems.

To overcome these problems, several authors have proposed different BIST techniques where signals are internally generated and/or analysed [START_REF] Toner | A BIST scheme for an SNR test of a sigmadelta ADC[END_REF][START_REF] Toner | A BIST Technique for a Frequency Response and Intermodulation Distortion Test of a Sigma-Delta ADC[END_REF][START_REF] Ohletz | Hybrid Built In Self Test (HBIST) for Mixed Analog/Digital Integrated Circuits[END_REF][START_REF] Sunter | A simplified polynomial-fitting algorithm for DAC and ADC BIST[END_REF][START_REF] Azais | Towards an ADC BIST scheme using the histogram test technique[END_REF][START_REF] Azais | Implementation of a linear histogram BIST for ADCs[END_REF][START_REF] Arabi | A New Built-In Self Test Approach For Digital-to-Analog and Analog-to-Digital Converters[END_REF]. Another possible and less expensive solution consists in using DFT techniques to internally transform the analogue signals into digital signals that are made controllable and observable from the chip I/Os [START_REF] Ohletz | Hybrid Built In Self Test (HBIST) for Mixed Analog/Digital Integrated Circuits[END_REF][START_REF] Nagi | A Signature Analyzer for Analog and Mixed-Signal Circuits[END_REF][START_REF] Mahoney | DSP-Based Testing of Analog and Mixed-Signal Circuits[END_REF]. As a result, only digital signals are externally handled by an inexpensive "digital" test equipment (Low Cost Tester).

In current systems, it should be mentioned that converters (Analogue to Digital Converters-or ADCs-and Digital to Analogue Converters -or DACs) are among the main components of any mixed-signal chip. Nowadays, many ADCs and DACs may be implemented in a complex SoC or SiP. For instance, Figure 1 proposes a microphotography of the PNX8327 device for set-top box applications (digital audio and video decoder) where 2 ADCs and 7 DACs are embedded in the same SiP.

In this context, this paper proposes an original DFT technique called "Analogue Network of Converters" (ANC) that permits to test the whole set of embedded ADCs and DACs. An extremely small circuitry is added to the original chip allowing to apply a fully digital test approach to the System-in-Package/System-on-Chip. In the remainder of the paper, section 2 gives the fundamental principle of the ANC technique where the test of a set of n DACs and m ADCs is associated to a system of equations where the converter characteristics are the unknowns. In section 3, the proposed ANC technique is validated through simulations and measurements. Finally, section 4 gives some concluding remarks. A crucial parameter, for ADC and DAC testing, is the Integral Non Linearity (INL). In addition, a long list of dynamic parameters is also considered [START_REF] Mahoney | DSP-Based Testing of Analog and Mixed-Signal Circuits[END_REF][START_REF] Ozev | Seamless test of digital components in mixed-signal paths[END_REF][START_REF] Janik | Estimation of A/D Converter Nonlinearities from complex Spectrum[END_REF]. For most of the application domains, two of the key dynamic parameters are:

• Total Harmonic Distortion (THD),

• Spurious Free Dynamic Range (SFDR),

The THD and SFDR are computed [START_REF] Mahoney | DSP-Based Testing of Analog and Mixed-Signal Circuits[END_REF] from the harmonic values appearing in the spectrum of the output signal. And the INL may also be derived from the same harmonic values [START_REF] Ozev | Seamless test of digital components in mixed-signal paths[END_REF].

That means that a very common way to estimate the dynamic parameters of a given converter is to perform a spectral measurement, i.e. to apply a single tone sine wave signal to the converter input and compute the FFT of the output signal. The obtained harmonic values are then used to compute the dynamic parameters.

Considering for instance the test of a single ADC using accurate instruments, it has been demonstrated [START_REF] Ozev | Seamless test of digital components in mixed-signal paths[END_REF] that the output signal can be represented by [START_REF] Toner | A BIST scheme for an SNR test of a sigmadelta ADC[END_REF]. This equation includes an ideal sampled sine wave x(n) and the sum of all the harmonic values introduced by the converter errors.

( ) ( ) ( ) ( ) ∑ ≥ θ + θ + = 0 k 0 n converter k k cos H n x n s (1)
In equation ( 1), n is the sample index, θ 0 the initial phase shift, converter k H the amplitude of the k th harmonic and θ n is the nominal sampling phase

n M P 2 n ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ π = θ (2)
where M is the number of samples and P the number of periods in the record.

Given the above comments, it clearly appears that accurate measurement of the set of harmonic values converter k H of the output signal is a crucial point for any converter testing technique. The above equation may also apply to the test of a single DAC, because the analogue output signal is converted into a digital sample set before the FFT analysis.

B. Analog Network of Converters

Considering a complex system with several ADCs and DACs, the objective of this paper is to measure the harmonic values converter k H of each converter output signal using a fully digital way [START_REF] Kerzérho | A DFT technique to test a complete set of ADCs and DACs embedded in a complex SiP or SoC[END_REF]. To be fully digital from an outside chip perspective, a very simple analogue circuitry is added to the system to: -realize the sum of any combination of DAC outputs, -connect this sum to any combination of ADC inputs.

This DFT technique requests an extremely simple and limited circuitry as illustrated in Figure 2. A simple OPAMP-based analogue adder can be used to implement the proposed DFT. The multiplexer control signal Ii allows to connect the corresponding DACi. In the same way, the multiplexer control signal Oj allows to connect the corresponding ADCj. When n DACs are connected with m ADCs, this is called a configuration C(n,m). Using configuration C(1,1), the spectrum of the output signal can be computed and we can extract the values of the harmonics

( ) 1 , 1 C k H
. But in this case, the output signal includes both the errors of DAC1 and the errors of ADC1. In other words, the spectrum includes the harmonic contribution of DAC1 as well as the harmonic contribution of ADC1. So, due to the linearity of the system, we can write the following equation:

( ) In this equation, we assume that the harmonic amplitudes created by the DAC are negligible with respect to the fundamental amplitude of the signal. Thus, we can consider the signal driving the ADC as a single tone signal. This working hypothesis will be verified in the validation phase described in section 4.

1 ADC k 1 DAC k m k measure k H H H H + = = (3) DAC1 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 ADC1 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 Σ DAC2 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 DAC3 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 DACN 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 ADC2 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 ADC3 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 ADCM 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 0 I 1 I 2 I 3 ..I N O 1 O 2 O 3 … O M DIGITAL DIGITAL ANALOGUE
Thanks to equation ( 3), we obtain a relation between the harmonic contribution of the different converters. Indeed, in equation (3), the left member is known (it corresponds to the amplitude of the k th spectral bin measured at the output of the ADC), while the right member represents the unknowns.

This small example demonstrates the relationship between one configuration and its resulting equation, which leads to the fundamental idea of the ANC DFT technique. By using different configurations C(n,m) we are able to obtain a set of different equations. So, with an adequate set of configurations (i.e system of equations), we expect to be able to fully determine the set of unknowns, i.e. the individual harmonic contribution of each converter. The ANC DFT technique creates a duality between the configurations and the equations allowing the estimation of the harmonic contributions of each converter. The next section explores the space of possible configurations to obtain an appropriate set of equations.

C. Test Strategy with C(1,1) & C(2,1) configuration

The ANC principle consists in using different hardware configurations in terms of converter interconnections. Then, the idea is to find an adequate test setup to discriminate the influence of each converter on the final response. In practice, the only test setup parameters we can easily control are the phase and the amplitude of the digital stimulus. In this section, two configurations, using DAC1, DAC2 and ADC1, are studied in order to discriminate their harmonic contributions.

In summary, the proposed test strategy, described in [START_REF] Kerzérho | A DFT technique to test a complete set of ADCs and DACs embedded in a complex SiP or SoC[END_REF], is composed of five successive tests. Each test consists in an acquisition and a spectral analysis (by Fast Fourier Transform) to evaluate harmonic bins. We obtain a 5equation system for each harmonic bin:

FS k FS k a , m k 1 Hadc 1 Hdac H + = FS k FS k b , m k 1 Hadc 2 Hdac H + = 2 / FS k 2 / FS k c , m k 1 Hadc 2 Hdac H + = (4)
( )

2 / FS k 2 / FS k FS k d , m k 1 Hadc k cos 2 Hdac 1 Hdac H + π + = ( ) ( ) 2 FS k 1 2 / FS k FS k e , m k k cos 1 Hadc k cos 2 Hdac 1 Hdac H ϕ + ϕ + = with ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - π = ϕ 4 1 cos ar 2 1 , ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - π = ϕ 4 1 cos ar 2
The three first tests use a C(1,1) configuration. The two first equations result from a test at full-scale (FS) through DAC1 and DAC2 respectively, then ADC1 (the converters are considered of same analogue full-scale range). The third equation results from a test at half full-scale (FS/2) through DAC2 and ADC1. The two last tests use a C(2,1) configuration. The two last equations result from tests that combine the outputs of the two DACs. The output range of DAC1 is full-scale while the output level of DAC2 is half full-scale. The two DAC input sine-waves are relatively phase shifted, in order to obtain output signal with a relative control phase shifted. After the sum of this relative phase shifted DAC output, we obtain the following particular amplitudes of the input signal of ADC1:

• full-scale for the fourth test

• half full-scale for the fifth test.

This system of independent equations is sufficient to calculate the values of the required harmonic contributions (

FS k FS k FS k 1 Hadc , 2 Hdac , 1 Hdac
). It allows thus a fully independent characterization of the three converters of the C(2,1) configuration in terms of harmonic contributions.

D. Configurations C(n,m)

Thanks to the converter characterization obtained from the C(1,1) and C(2,1) configuration, it is then easy to test every other converter embedded in the same complex chip. The idea is to use one of the three previously characterized converters as a measurement instrument whose non-ideal features are now well known.

III. VALIDATION

Simulation results have provided a first validation of the proposed method [START_REF] Kerzérho | A DFT technique to test a complete set of ADCs and DACs embedded in a complex SiP or SoC[END_REF]. In this section, we show the results of the first experiments, evaluating the effectiveness of the method in the estimation of ADC distortion components. The experimental setup is first introduced, then the experimental protocol is described and finally the results are presented. The performance of the proposed test strategy is discussed in terms of difference found on SFDR and THD results.

A. Experimental setup

In order to experiment the test strategy, we need at least two DACs and one ADC to establish a C(2,1) configuration (Figure 3).

For a faster implementation, we have decided to carry out the first experiments using a 2-channel Arbitrary Wave Generator (AWG) to play the role of the two DACs. Each channel is composed of a 12-bit DAC and an amplifier.

In addition, two low-pass filters are placed between the AWG outputs and the adder in order to reduce the second and third harmonics generated by the AWG. These two harmonics are too high to emulate correctly two 12bit DACs.

A combiner/splitter HP11667B is used to sum the signals.

The stand-alone ADC used for experimentation is a TDA9910 Philips. For comparison, reference measurements have been collected using a standard setup, as shown in Figure 4.

The test signal is generated by an AWG, and the spurious components (noise and harmonics) are removed by a bandpass filter centred on the test frequency.

B. Experimental results

The objective is to compare the values of the harmonic components of the ADC evaluated using the proposed strategy to the ones obtained using a standard test configuration. To this purpose, we have conducted this experiment on 3 different ADC samples, chosen in a THD range from -70 to -63dB. Each data acquisition is achieved 20 times and averaged, in order to reduce the noise influence. The harmonic components are individually evaluated, from the 2 nd to the 10 th order, then the THD and SFDR parameters are calculated. The results presented in Table 1 are experimental results for one ADC. The harmonics 2 to 10 have been measured by the classical method as well as the new method. The last column corresponds to the difference between the two measures.

From the results shown in Table 1, we can conclude that the new method correlates with the reference when the level of the harmonic is higher than -80dB, which corresponds to the level of the highest noise component of this device. To estimate the performances of the new method, three components with different characteristics have been chosen: a component with a high quality THD (THD=-70.4dB), another component with THD at the tolerance limit (THD=-66.4dB) and a component that wouldn't pass the test (THD=-63.2dB, beyond device specifications). Analysing the results of Tables 2 and3, the maximal difference between the two evaluation methods is 1dB for THD and 2.5dB for SFDR. Correlating the results of harmonic measurements represents a big challenge, because of their high sensitivity to the variations of the input sine wave amplitude [START_REF] De Vries | Decreasing the sensitivity of ADC test parameters by means of wobbling[END_REF]. Due to this sensitivity a maximum error of reproducibility of +/-3dB on individual harmonics and SFDR is accepted for a 12-bit resolution. For THD parameter, a +/-1.5dB variation is generally accepted. As a consequence, the results obtained by the new method are very acceptable and promising.

ADC

IV. CONCLUSION This paper has introduced the novel concept of "Analogue Network of Converters" (ANC) to test distortion of embedded converters. Thanks to this approach, we have demonstrated that it is possible to achieve a fully digital test of the set of converters and consequently to use low cost testers. The technique has been validated not only by electrical simulation but also by silicon experiments.

Figure 1 :

 1 Figure 1: SiP for set-top box applications

Figure 2 :

 2 Figure 2: The ANC DFT technique

Figure 3 :

 3 Figure 3: Experimental setup

Figure 4 :

 4 Figure 4: Reference measurement setup

Table 1 :

 1 First harmonics estimation for one ADC Table2and 3 give a comparison between the two methods, concerning THD and SFDR measurements.

	Harmonic	Reference	Method	Measurement
	number	values (dB)	results (dB)	difference (dB)
	2	-69.2	-68.1	-1.1
	3	-67.9	-67.8	0.1
	4	-72.5	-70.8	-1.7
	5	-70.8	-71.8	1.0
	6	-95.7	-86.2	-9.5
	7	-78.2	-79.4	1.2
	8	-102.9	-97.6	-5.2
	9	-78.7	-78.9	0.2
	10	-92.0	-86.7	-5.3

Table 2 :

 2 THD estimation for 3 ADCs

		THD reference	THD	Measurement
	number	(dB)	measured (dB)	difference
				(dB)
	1	-63.2	-63.6	0.4
	2	-66.4	-66.0	-0.4
	3	-70.4	-69.4	-1
	ADC	SFDR	SFDR	Measurements
	number	reference (dB)	measured (dB)	difference
				(dB)
	1	67.8	65.3	2.5
	2	66.9	66.9	0
	3	70.7	72.1	-1.4

Table 3 :

 3 SFDR estimation for 3 ADCs


Carles Fernandez Prades
email: carles.fernandez@cttc.cat

Pau Closas
email: pau.closas@cttc.cat

Jordi Vilà Valls

Carles Fernández-Prades

Jordi Vilà-Valls
email: jordi.vilavalls@gipsa-lab.grenoble-inp.fr

Nonlinear filtering for ultra-tight gnss/ins integration

come   L'archive ouverte pluridisciplinaire

Index Terms-GNSS/INS integration, cubature rules, synchronization, square-root Cubature Kalman filter

I. INTRODUCTION

Accuracy and availability are the main objectives of a navigation system. Since Global Navigation Satellite Systems (GNSS) performance is not sufficient for all applications in geodesy and navigation, it is common to hybridize them with Inertial Navigation Systems (INS), in order to combine the only short-time stability of an inertial sensor and the longtime stability but noise behavior of a GNSS receiver, taking advantage of their complementary nature.

In most GNSS/INS integrations the position and velocity estimates of a GNSS receiver are used as observations in a navigation filter for the estimation of INS errors, making possible a reduction of GNSS-noise, bridging of GNSS outages, and even GNSS measurements are also usable when fewer than four satellites are tracked. This combination, known as loosely coupled integration, provides good performance combined with simplicity of integration. The mathematical problem posed by this integration can be cast into a nonlinear filter, that traditionally has been solved by means of the Extended Kalman Filter (EKF).

In case of high dynamic applications or in jamming environments, information about the receiver dynamics measured by the inertial sensor can be used to enhance the This work has been partially supported by the Spanish Science and Technology Commission: CENIT2007-2002 (TIMI) and TEC2008-02685/TEC (NARRA), and by the European Commission in the framework of the COST Action IC0803 (RFCSET) and the FP7 Network of Excellence in Wireless COMmunications NEWCOM++ (contract n. 216715).

pseudorange/Doppler estimations of the GNSS receiver. This principle is known as tightly coupled integration, and it poses a nonlinear problem. Examples of this approach can be found in [1]- [START_REF] Li | Sigma-point Kalman filtering for tightly coupled GPS/INS integration[END_REF], making use of the UKF, Sigma-Point Kalman Filter (SPKF), and Particle Filtering (PF) methods. In [START_REF] Bernal | Tight GNSS/INS integration as a constrained least-squares problem[END_REF], the approach is cast into a convex, linearly constrained leastsquares problem and solved via dual decomposition.

An extension of this concept is the ultra-tightly coupled integration. In this case the integrating navigation filter is implemented as one element of the receiver tracking loop. An optimal controlling value for the numerical oscillator can be computed by using inertial information in combination with the in-phase (I) and quadrature (Q) signals at the receiver. Examples of this approach can be found in [START_REF] Petovello | Comparison of vector-based software receiver implementations with application to ultra-tight GPS/INS integration[END_REF]- [START_REF] Babu | Ultra-tight GPS/INS/PL integration: a system concept and performance analysis[END_REF], where the errors in position and velocity between the measured and estimated values act as a link between I & Q samples and INS measurements, and lead to slight modifications of the tracking loops.

This paper proposes a new approach to ultra-tight GNSS/INS integration by deriving the evolution of the synchronization parameters from the INS measurements, considering the eventual biased behavior of the latter, and taking advantage of the numerical efficiency of the Squareroot Cubature Kalman Filter (SCKF), a nonlinear Gaussian filter recently developed [START_REF] Arasaratnam | Cubature Kalman Filters[END_REF], that could substitute the structure of conventional tracking loops.

II. SYSTEM MODEL

II-A. GNSS model

Considering the signals coming from M s satellites in a multipath channel, the received signal is

r(t) = Ms-1 i=0 M (i)-1 m=0 α m,i (t)
e j(2πfc(t-τm,i(t))+φm,i)

× s T,i (t -τ m,i (t)) + w(t) (1)
where {•} denotes the real part of the complex-valued quantity in the brackets, M (i) is the number of propagation paths from the i th satellite to the receiver, f c is the carrier frequency, {α m,i (t), τ m,i (t), φ m,i } are the amplitude, delay and carrier phase of the i th satellite in the m th path, w(t) is additive white Gaussian noise, and s T,i (t) is the Direct-Sequence Spread-Spectrum (DS-SS) baseband signal transmitted by the i th satellite. In order to compute the position, the receiver should estimate the pseudoranges (directly related to time delay of the line-of-sight (LOS) paths, τ 0,i (t)) and pseudorange rate (related to the evolution of the carrier phase). Considering the orthogonality between spreading codes, and the usual parallel architecture of a GNSS receiver, we will simplify the signal model to

r(t) = α 0 (t) e j(2πfc(t-τ0(t))+φ0) s T (t -τ 0 (t)) + n(t)
(2) where the subindex 0 denotes the direct path, and n(t) gathers all the unwanted signals (multiuser interference, multipath and termal noise) in a term which is considered Gaussian but not white, thus admitting a degree of correlation of the unwanted terms with the desired signal.

When this signal impinges the antenna, it passes through a Low Noise Amplifier (LNA) and a downconversion stage to an intermediate frequency. This can be modeled as

y(t) = (r(t) * h LN A (t)) e -j(2πf LO t+φ LO ) (3)
where * stands for the convolution operator, h LN A (t) is the impulse response for the RF front-end filter, f LO is the frequency of the local oscillator and φ LO is its carrier phase.

Defining

f IF = f c -f LO , ψ 0 = φ 0 -φ LO , s R (t)
as the filtered version of s T (t), a 0 (t) as the signal amplitude, and η(t) as the filtered and down-converted noise term n(t), we can write

y(t) = a 0 (t) s R (t -τ 0 (t)) e j(2πf IF t-2πfcτ0(t)+ψ0) +η(t) (4
) Sampling this signal at t k = kT s , we obtain the discrete version of the received signal

y k = a 0,k s R (t k -τ 0,k )e j(2πf IF t k -2πfcτ 0,k +ψ0) + η k (5)
where the subindex k refers to t k instants. Defining p (s) and p (u) as the position vector of the satellite and the user receiver, respectively, and neglecting atmospheric effects, the delay can be written as

τ 0 (t) = 1 c p (s) (t Tx ) -p (u) (t) = t -t Tx (6)
being c the speed of light, t Tx is the transmission instant, and • stands for the Euclidean norm. Approximating ( 6) by its second-order Taylor expansion around an arbitrary point t 0 we get

p (s) (t Tx ) -p (u) (t) ≈ p (s) (t 0 -τ 0 (t 0 )) -p (u) (t 0 ) + + d dt p (s) (t 0 -τ 0 (t 0 )) -p (u) (t 0 ) (t -t 0 ) + + 1 2 d 2 dt 2 p (s) (t 0 -τ 0 (t 0 )) -p (u) (t 0 ) (t -t 0 ) 2 (7)
Substituting equations ( 7) and ( 6) into equation ( 5) yields

y k = a 0,k s R (t k -τ 0,k )e j(2π fk t k + ψk ) + η k (8)
where we have made use of the following definitions:

fk =f IF - f c c v r (9) ψk = - 2πf c c p (s) (t 0 -τ 0 (t 0 )) -p (u) (t 0 ) + -v r t 0 + a r (t k -t 0 ) 2 + ψ 0 (10) v r = d dt p (s) (t -τ 0 (t)) -p (u) (t) t=t0 (11)
a r = 1 2 d 2 dt 2 p (s) (t -τ 0 (t)) -p (u) (t) t=t0 (12)
Terms v r and a r are usually referred to as radial velocity and radial acceleration, respectively, while the term f d = fc c v r is commonly known as the Doppler shift.

II-B. INS model

Inertial measurement units (IMUs) typically contain three orthogonal rate-gyroscopes and three orthogonal accelerometers, measuring angular velocity and linear acceleration respectively. Let us assume a strapdown INS system in with the inertial sensors are mounted rigidly onto the receiver device, and therefore outputs quantities measured in a body frame. A typical strapdown navigation algorithm works as follows: the rate-gyroscope measurements ω b (t) = ω bx (t) ω by (t) ω bz (t)

T (that is, the angular velocity) are integrated to obtain the orientation of the tracked device. Rate-gyroscopes are known to suffer from bias, that is, the average output (usually given in o /h) from the gyroscope when it is not undergoing any rotation is not zero. A constant bias error, when integrated, causes an angular error which grows linearly with time. This effect can be modeled by a bias term, in the form

ωb (t) = ω b (t) -δω b (t), (13)
where δω b (t) = δω bx (t) δω by (t) δω bz (t) T .

Since those measurements are given with respect to a body frame and we need them expressed in a global frame, a rotation matrix C g b such that ωg (t) = C g b ωb (t) must be tracked through time. The rate of change of C g b at time t is given by d dt C g b (t). Using the small angle approximation, it comes out that the orientation algorithm must solve the differential equation

d dt C g b (t) = C g b (t)Ω(t) (14)
where Ω(t) is the skew-symmetric form of the angular vector ωb (t):

Ω(t) =   0 -ω bz (t) ωby (t) ωbz (t) 0 -ω bx (t) -ω by (t) ωbx (t) 0   (15)
Taking into account that IMUs provide samples of the angular velocity at a constant rate of period T i , equation ( 14) can be solved as

C g b (t + T i ) = C g b (t)e t+T i t Ω(t)dt (16)
Solving the integral as t+Ti t Ω(t)dt = Ω(t)T i and performing a Taylor expansion of the exponential term yields

C g b (t + T i ) =C g b (t) I - sin (| ωb T i |) | ωb T i | Ω(t)T i + + 1 -cos (| ωb T i |) | ωb T i | (Ω(t)T i ) 2 (17)
To track the position of an INS the acceleration signal a b (t) = a bx (t) a by (t) a bz (t)

T obtained from the accelerometers is projected into the global frame of reference. Again, accelerometers are also biased. A constant bias error, when double integrated, causes an error in position which grows quadratically with time. This is modeled with a bias term

ǎb (t) = a b (t) -δa b (t), (18)
where δa b (t) = δa bx (t) δa by (t) δa bz (t) T . Thus, the acceleration in a global frame can we written as ǎg (t) = C g b (t)ǎ b (t). Acceleration due to gravity is then subtracted and the remaining acceleration is integrated once to obtain velocity, and again to obtain displacement. In discrete form:

v g,k =v g,k-1 + T i (ǎ g,k -g g ) (19) p g,k =p g,k-1 + T i v g,k (20)
We also define i.i.d zero-mean Gaussian terms w ω,k and w a,k modeling gyroscopes and accelerometers noise.

II-C. GNSS/INS coupling

We propose here a new GNSS/INS coupling approach, which consists on the derivation of the evolution of the synchronization parameters from the IMU measures. Taking into account the geometry of the satellite-user link, we can write that

τ k = τ k-1 - 1 c (p g,k -p g,k-1 ) T d k-1 + + 1 c p (s) (t k -τ k ) -p (s) (t k-1 -τ k-1 ) T d k-1 + w τ,k = τ k-1 + β (τ ) k + w τ,k (21)
where d k is the unitary vector pointing from the user receiver to the satellite:

d k = p (s) (t k -τ k ) -p g,k p (s) (t k -τ k ) -p g,k (22)
The term p (s) (t k -τ k ) represents the position from where the satellite emitted the signal received at t k , an information that can be inferred from the navigation message. Since τ k is what we are looking for, we make the approximation

p (s) (t k -τ k ) ≈ p (s) (t k -τ k-1
). Proceeding in the same way, from ( 9) we can derive

fk = fk-1 + f c c v (s) (t k-1 -τ k-1 ) -v g,k-1 T d k-1 + -v (s) (t k -τ k ) -v g,k T d k + w f,k = fk-1 + β (f ) k + w f,k (23)
where v g,k is the velocity given by the IMU, and v (s) (t) is the velocity of the satellite that can be extracted from the navigation message. Again, we will accept

v (s) (t k -τ k ) ≈ v (s) (t k -τ k-1
). The phase evolution can be predicted from [START_REF] Van Der Merwe | Sigma-point Kalman filters for probabilistic inference in dynamic state-space models[END_REF] as

ψk = ψk-1 - 2πf c c p (s) (t k -τ k ) -p g,t k + -p (s) (t k-1 -τ k-1 ) -p g,t k-1 + + v (s) (t k-1 -τ k-1 ) -v g,k-1 T d k-1 t 0 + -v (s) (t k -τ k ) -v g,k T d k t 0 + -a (s) (t k-1 -τ k-1 ) -a g,k-1 T d k-1 (t k-1 -t 0 ) 2 + + a (s) (t k -τ k ) -a g,k T d k (t k -t 0 ) 2 + w ψ,k = ψk-1 + β (ψ) k + w ψ,k (24)
where a g,k is the acceleration measured by the IMU, and a (s) (t) is the acceleration of the satellite, also extracted from the navigation message. Terms w τ,k , w f,k and w ψ,k are white Gaussian i.i.d. noise sequences. We take a state-space approach in order to use an optimal Bayesian filtering solution. In equation ( 8), we have defined the discrete measurement equation. As we want to consider the biased behavior of the IMU, we introduce the bias of the rate-gyroscopes, δω b,k , and accelerometers, δa b,k , into the state, using a random walk evolution. We define the vector state

x k = τ k , fk , ψk , δω T b,k , δa T b,k T , (25)
and its evolution as

x k = Φ k x k-1 + u k + w k ( 26
)
where Φ k = I 9×9 is the identity matrix, the input vector

is u k = β (τ ) k , β (f ) k , β (ψ)
k , 0 1×6 T , and the noise vector is

w k = w τ,k , w f,k , w ψ,k , w T ω,k , w T a,k T .

III. SIGMA-POINT NONLINEAR FILTERING

When considering a nonlinear filtering problem, we can consider different suboptimal solutions. A common solution is the well-known EKF, a solution that is useful for mild nonlinearities but fails otherwise. Alternative solutions have been proposed in the last decade, among them, the family of Sequential Monte Carlo (SMC) methods [START_REF]Sequential Monte Carlo methods in practice[END_REF] and the family of Sigma-Point Filters (SPFs) [START_REF] Van Der Merwe | Sigma-point Kalman filters for probabilistic inference in dynamic state-space models[END_REF].

SMC methods provide a framework to deal with nonlinear, non-Gaussian problems using a stochastic sampling approach to numerically approximate the integrals of the optimal solution. This broad suitability comes at the expense of a high computational load, that makes this solution difficult to embed in digital light processors or real-time applications.

SPFs use a deterministic sampling to approximate the integrals involved in the Bayesian filter solution. The key point is to assume that the measurement and process noise are independent random Gaussian processes. This leads to Gaussian state transition and measurement likelihood densities, which in turn reverts to a Gaussian posterior density [START_REF] Ito | Gaussian filters for nonlinear filtering problems[END_REF]. Recently, a new SPF method, and its squareroot version, has been derived: the Square-root Cubature Kalman Filter (SCKF) [START_REF] Arasaratnam | Cubature Kalman Filters[END_REF]. The SCKF uses cubature rules to approximate the optimal solution, and propagates the square root of the covariance matrix instead of the covariance itself, thus preserving its symmetry and positive-definiteness, and avoiding its inversion in the update step. The SCKF is numerically stable, accurate, and easily extendable to high dimensional problems, mitigating the curse of dimensionality and divergence effects.

We sketch in Algorithm 1 the square-root, derivativefree scheme. In steps 8 and 17, S = Tria (A) denotes a general triangularization algorithm (for instance, the QR decomposition), where A ∈ R p×q , p < q, and S is a lower triangular matrix. For computational reasons, we prefer to keep the square root as a triangular matrix of the dimension p × p. This can be achieved by the thin QR decomposition [12, § 5.2, Theorem 5.2.2], which has a computational complexity of O qp 2 flops. Algorithm 2 provides the way to generate the deterministic samples, using the third-degree spherical-radial rule, to approximate a standard Gaussian weighted integral.

IV. COMPUTER SIMULATIONS

The proposed integration solution was tested by computer simulation. We used a realistic GPS signal generator which considered aspects such as constellation geometry, receiver/satellites relative motion and time-varying Algorithm 1 Square-root, derivative-free nonlinear Kalman Filter Require: y 1:K , u 0:K , x0 , Σ x,0 = S x,0|0 S T x,0|0 , Σ w,0 , Σ η,0 . Initialization:

1: Define sigma-points and weights {ξ i , ω i } i=1,...,L by using Algorithm 2. 2: Set W = diag( √ ω i )

Tracking:

3: for k = 1 to ∞ do 4:
Time update:

5:
Evaluate the sigma points:

x i,k-1|k-1 = S x,k-1|k-1 ξ i + xk-1|k-1 , i = 1, ..., L. 6:
Evaluate the propagated sigma points using (26): xi,k|k-1 = f (x i,k-1|k-1 , u k-1 ).

7:

Estimate the predicted state:

xk|k-1 = L i=1 ω i xi,k|k-1 . 8:
Estimate the square-root factor of the predicted error covariance: S x,k|k-1 = Tria Xk|k-1 S Σ w,k-1 , where:

S Σ w,k-1 is a square-root factor of Σ w,k-1 such that Σ w,k-1 = S Σ w,k-1 S T Σ w,k-1 ,

and

Xk|k-1 = x1,k|k-1 -xk|k-1 • • • xL,k|k-1 -xk|k-1 W.

9:

Measurement update:

10:

Evaluate the sigma points:

x i,k|k-1 = S x,k|k-1 ξ i + xk|k-1 , i = 1, ..., L.

11:

Evaluate the propagated sigma points using (8): ỹi,k|k-1 = h(x i,k|k-1 , u k ).

12:

Estimate the predicted measurement: ŷk|k-1 = L i=1 ω i ỹi,k|k-1 .

13:

Estimate the square-root of the innovation covariance matrix: S y,k|k-1 = Tria Y k|k-1 S Σ η,k , where: S Σ η,k denotes a square-root factor of Σ η,k such that Σ η,k = S Σ η,k S T Σ n,k , and

Y k|k-1 = ỹ1,k|k-1 -ŷk|k-1 • • • ỹL,k|k-1 -ŷk|k-1 W.

14:

Estimate the cross-covariance matrix Σ xy,k|k-1 = X k|k-1 Y T k|k-1 , where:

X k|k-1 = x 1,k|k-1 -xk|k-1 • • • x L,k|k-1 -xk|k-1 W.

15:

Estimate the Kalman gain K k = Σ xy,k|k-1 /S T y,k|k-1 /S y,k|k-1 .

16:

Estimate the updated state xk|k = xk|k-1 + K k y k -ŷk|k-1 . 17:

Estimate the square-root factor of the corresponding error covariance: typical MEMS-based IMU) over 1 second. The nominal Carrier-to-Noise density ratio was set to C/N 0 = 46 dB-Hz, although we took into account that the instantaneous value of such parameter will depend on the elevation angle and on satellite's antenna radiation pattern. On the other hand, the considered GPS receiver had a pre-correlation filter of 2 MHz, a sampling frequency of 5.7143 MHz, an intermediate frequency of 4.308 MHz, and the tracked signal corresponded to that of the C/A code transmitted in the L1 link. With this setup, the receiver was able to track the synchronization parameters of the satellite, with an initial ambiguity given by the acquisition process, i.e., a standard deviation of half the sampling period. Figure 1 shows the Root Mean Squared Error (RMSE) results obtained with the SCKF in estimating the pseudorange. In Figure 1 we also plotted the Cramér-Rao Bound (CRB), which is the minimum variance that any unbiased estimator can achieve.

S x,k|k = Tria X k|k-1 -K k Y k|k-1 K k S Σ η,k .
The CRB in its standard formulation is valid for ML-like estimators, that is to say when no use of prior information is done. However, this is not our case since we are considering prior data, which is expressed in the state-evolution models proposed in Section II-C. Therefore, we must resort to the Bayesian CRB as a valid benchmark for our algorithm.

From the figure it follows that the obtained pseudorange variances are in accordance to the theoretical bound. In addition, we can see the enhancement given by the use of prior information.

V. CONCLUSIONS

This paper presented a new approach to GNSS/INS integration. The main objective was to establish the pertinent mathematical relationships that allow the ultra-tight coupling of both systems, casting the problem into the Bayesian framework and proposing a solution based on Gaussian nonlinear filtering. Detailed computer simulations confirmed the validity of the approach.

VI. REFERENCES

[1] Y. Yi amd D. A. Grejner-Brzezinska, "Tightly-coupled GPS/INS integration using unscented Kalman filter and particle filter," in Proc. ION

  18: end for delay/Doppler signal parameterizations. In particular, we tested the tracking of the synchronization parameters of one satellite by a static receiver (considering the biases of a Algorithm 2 Generation of Sigma-Points and weights for third-degree spherical-radial cubature rule 1: Set M = dim(x) and L = 2M . 2: Set the cubature points ξ i = √ M I M ×M -I M ×M i , where [•] i=1,...,L indicates the i-th column. 3: Set the cubature weights ω i = 1 2M , i = 1, ..., L.

Fig. 1 .

 1 Fig. 1. RMSE of pseudorange estimates, CRB and Bayesian CRB over time.

  GNSS, Fort Worth, TX, Sept. 2006, pp. 2182-2191. [2] J. Rezaie, B. Moshiri, B. N. Araabi, and A. Asadian, "GPS/INS integration using nonlinear blending filters," in Proc. SICE'07, Kagawa, Japan, Sept. 2007, pp. 1674-1680.