Husain Mubashir
email: mubashir.rehmani@lip6.fr

Rehmani

Channel Assortment Strategy for Reliable Communication in Multi-Hop Cognitive Radio Networks

I. INTRODUCTION

Channel selection plays a vital role in efficient and reliable data dissemination. In the context of Cognitive Radio Network (CRN) [START_REF] Akyildiz | Next generation/dynamic spectrum access/cognitive radio wireless networks: a survey[END_REF], channel selection is more challenging due to the traffic pattern and channels' occupancy of primary radio (PR) nodes. Moreover, Cognitive Radio (CR) transmissions should not degrade the reception quality of PR nodes and should be immediately interrupted whenever a neighboring PR activity is detected [START_REF] Khalife | Probabilistic path selection in opportunistic cognitive radio networks[END_REF]. Thus, it is essential and, however, extremely challenging for CR nodes to correctly select channels allowing reliable communication.

In contrast with previous approaches which focused on single-hop wireless infrastructure, the main problem we tackle here is the channel assortment to undergo efficient and reliable data dissemination throughout a multi-hop CRN. Basically, we want to increase the transmission coverage of CR nodes, and consequently, the global network dissemination coverage, by selecting good qualified channels for communicating. More specifically, our channel assortment strategy, named 'SURF' (an allusion for riding on channels), empowers CR nodes with the ability to infer, based on information regarding PR occupancy, the less PR-occupied channel to use. On the flip side, although less occupied channel by PR nodes provides higher space sharing to CR nodes, if all CR nodes decide to switch to it, a higher CR competition will be perceived, increasing the probability of collision and message losses. In this way, how to find a good compromise between transmission opportunity in terms of PR occupancy and the number of CR neighbors, constitutes the main goal of this paper.

II. CHANNEL ASSORTMENT STRATEGY

We consider a multi-hop CRN that has a set of: PR nodes, CR nodes, and network links connecting CR nodes. We assume that a set of frequency channels is available in the network and may be used by PR or CR nodes. We consider a slotted-based channel in which each channel is divided into total number of τ t = τ o + τ a time slots, where τ o and τ a are the slots occupied by PR nodes and the available slots, respectively. In essence, in order to be able to communicate in a CRN, CR nodes must be cognitively aware of their surroundings and able to make decisions about "when, where, and how" to communicate. The SURF goal is to support nodes in their decision about "where" to transmit or overhear, i.e. which channel to use. Here, the role of overhearing is for receiving data.

In a multi-hop context, at each hop, CR nodes will select for data dissemination the best channel in the channel classification provided by SURF. In order to limit energy consumption and communication overhead, the SURF assortment is performed in a distributed way and is based only on information locally inferred by CR nodes. In fact, by implementing the same strategy at the sender and the receiver, SURF helps both of them tune to the appropriate channel for undergoing transmissions or reception without the need of any prior information exchange or synchronization.

Consider the set C of total frequency channels. Since, the number of channels changes with time and location due to PR activity, each CR node has its own available channel set among total frequency channels C, which we consider provided by the spectrum sensing block [START_REF] Yucek | A survey of spectrum sensing algorithms for cognitive radio applications[END_REF]. We call this available channel set as Acs, where ∀Acs ∈ C . Acs describes the channels set a CR node has access to. After sensing them and according to the SURF strategy, each CR node classifies the sensed channels according to their availabilities and selects for transmission and/or overhearing the best weighted channel. SURF assigns weights for channels according to the Eq. 1:

∀i ∈ C :

P (i) w = e -P R (i) o × CR (i) o (1) P (i)
w describes the availability level of a channel (i) and is calculated based on the occupancy of PR (i.e. P R We rely on methods in [START_REF] Yucek | A survey of spectrum sensing algorithms for cognitive radio applications[END_REF] and assume that the spectrum opportunity map of PR nodes and CR neighbor (CR

(i) n ) information is available to CR nodes. Note that CR (i)
n is the number of potential CR neighbors competing for the same channel. This map is then used to compute the percentage of the channel i occupied by PR nodes (i.e., the ratio between the number of PR nodes and the total number of time slots τ t ), or P R as . If a central entity could control the channel access for CR nodes, using a channel with higher available space could allow achieving higher number of 1-hop CR receivers. Since no central entity is present, there is no way to prevent collision and message losses when the number of CR nodes competing for the same channel increases. To cope with this and to allow nodes selecting channel having a good compromise between the number of CR neighbor receivers and the number of competing CR transmitters, we use a Tenancy Factor, named β. β is used to compute the CR occupancy CR

Algorithm 1 Algorithm for CR occupancy's computation

if CR (i) n <β then CR (i) o ← τa/(β -

CR (i) n ) else if CR (i) n = β then CR (i) o ← CR (i) as else CR (i) n >β then CR (i) o ← τa/τt.CR (i) n end if
number of CR neighbors on a particular channel for which the communication is still performed with a good probability of success. The goal is then, to maximize the chances of selecting channels that have higher but enough number of CR neighbors (close to β) allowing reliable communication.

Algorithm 1 shows how CR nodes calculates CR occupancy according to β and the number of CR neighbors competing for the channel. We performed extensive simulations in order to determine the appropriate value of β to be used. We ask readers to refer to [START_REF] Rehmani | Toward Reliable Contention-aware Data Dissemination in Multi-hop Cognitive Radio Ad Hoc Networks[END_REF] for more details.

III. PERFORMANCE ANALYSIS Results are generated from an average of 1000 simulations, with 95% of confidence intervals. We consider 30 PR nodes; uniformly distributed among the existing channels and τ t =6 total time slots for each channel. The transmission range is set to R = 250m and the average CR neighbor density d avg set to 20. The number of CR nodes is fixed to N=70 and randomly deployed within a square area of a 2 =707x707m 2 . TTL is used as the maximum number of hops required for a packet to traverse the whole network and is set to 2a R , i.e. T T L = 6. The Acs size is then set to 3 and 8 for 5 and 15 existing channel (Ch), respectively. The appropriate value of β is set to 10 and 18 when using Ch = 5 and Ch = 15 respectively, which provides the best tradeoff among number of receivers and loss ratio (cf. [START_REF] Rehmani | Toward Reliable Contention-aware Data Dissemination in Multi-hop Cognitive Radio Ad Hoc Networks[END_REF]).

We compare SURF with a random strategy (RD) and two variants of selective broadcasting strategy [START_REF] Kondareddy | Selective broadcasting in multi-hop cognitive radio networks[END_REF]: without (SB) and with (CA) centralized authority. CA can be thus used as an upper bound in message dissemination comparison, as it maximizes the number of receptions by performing overhearing over multiple channels, simultaneously. Since, our goal is to efficiently disseminate the data, therefore, we have chosen two performance metrics to compare these four strategies: (i) the average delivery ratio, which is the ratio of packet received by a particular CR node over total packets sent in the network and (ii) the average number of accumulative CR receivers per hop.

Fig. 1(a) compares the number of accumulative CR receivers at each hop of communication (i.e until TTL=0) for the four strategies. SURF allows the message dissemination to 55% of nodes in the network. It can be clearly seen that SURF outperforms RD and SB and compared to CA, only provides a decrease of 25% in performance (1.5 less transmissions). Note that CA have higher number of receivers at the cost of multiple transmissions (2.5 transmissions), while in RD and SURF, each CR node transmits only once.

Fig. 1(b) compares delivery ratio of RD, SB, CA, and SURF, as a function of CRs' ID. SURF increases considerably the delivery ratio compared to RD and SB and only reduces in 20% the performance compared to CA. In particular, for Ch=5 and Ch=15, SURF guarantees the delivery of approximately 60% of messages (single transmission), contrarily to less than 20% for the RD and SB (single and 2.5 transmissions, resp.) and 80% for CA (2.5 transmissions). These results show the good level of network connectivity provided by SURF, suitable for reliable dissemination.

IV. FUTURE WORKS

  (i) o ) and CR (i.e. CR (i) o ) nodes in this channel. If a CR node finds two or more channels having identical values of P w then it selects the one that has lower P R (i) o . If the problem persists, then, a channel is randomly selected among them.

  remaining available percentage of the channel i, i.e. 1 -P R (i) o , gives then the spectrum opportunity available for channel sharing among CR nodes, named CR (i)

Fig. 1 .

 1 (i) o of each channel i. β provides the upper bound in terms of Comparison of random (RD), selective broadcasting (SB), centralized authority (CA) and our strategy (SURF), when channels=5 and channels=15 in a CRN with 70 CR nodes. (a) Hop Count and average number of accumulative receivers. (b) CR Nodes' ID and Average delivery ratio.

We intend in future to investigate SURF performance under dynamic traffic by consideration of data rates and traffic volume generated by CR nodes and optimize data dissemination delay. Moreover, prediction and history can also be accounted in SURF to enhance the performance.

V. ACKNOWLEDGEMENTS The author would like to thank Aline C. Viana, Hicham Khalife, and Serge Fdida for the helpful discussions and suggestions, and for their very careful reading of the manuscript.


Bruno Torrésani
email: bruno.torresani@cmi.univ-mrs.fr

Emilie Villaron
email: emilie.villaron@cmi.univ-mrs.fr

HARMONIC HIDDEN MARKOV MODELS FOR THE STUDY OF EEG SIGNALS

A new approach for modelling multichannel signals via hidden states models in the time-frequency space is described. Multichannel signals are expanded using a local cosine basis, and the (timefrequency labelled) coefficients are modelled as multivariate random variables, whose distribution is governed by a (hidden) Markov chain. Several models are described, together with maximum likelihood estimation algorithms.

The model is applied to electroencephalogram data, and it is shown that variance-covariance matrices labelled by sensor and frequency indices can yield relevant informations on the analyzed signals. This is examplified by a case study on the characterization of alpha waves desynchronization in the context of multiple sclerosis disease.

INTRODUCTION

Many signal classes exhibit specific features that are both time and frequency localized. This motivates the use of time-frequency signal representations, such as short time Fourier or Gabor transforms (see [START_REF] Carmona | Practical Time-Frequency Analysis: continuous wavelet and Gabor transforms, with an implementation in S, volume 9 of Wavelet Analysis and its Applications[END_REF][START_REF] Flandrin | Temps-Fréquence. Traité des Nouvelles Technologies[END_REF], ?] and references therein), or orthonormal bases such as MDCT or Wilson bases. Such elementary waveform systems have more recently been exploited in sparse regression algorithms (matching pursuit and variants, basis pursuit,...) in various contexts. However, real signals often exhibit significant correlations between time-frequency coefficients, which are often poorly accounted for by classical sparse expansion methods. In the case of multichannel signals, such as EEG, MEG,... signals, correlations also exist between channels: information may be somewhat localized in the channel domain, but close channels (sensors) are often significantly correlated.

Sparse regression approaches are often not very good at describing time-correlated or frequency-correlated data, contrary to channel correlations for which dedicated techniques have been developed (see for example [START_REF] Fornasier | Recovery algorithms for vector valued data with joint sparsity constraints[END_REF][START_REF] Kowalski | Décompositions parcimonieuses et persistantes de signaux multicanaux. applications aux signaux MEEG[END_REF]10] and references). Motivated by applications to EEG/MEG signal analysis, we discuss here a stochastic model implementing simultaneously channel and frequency correlated signals, following the lines of [?, 1], with explicit modelling of non-stationarity. Indeed, standard stochastic signal models often rely on (at least implicit) stationarity assumptions. However localized phenomena in EEG signals such as alpha or gamma oscillations cannot be accounted for with such assumptions.

The approach we propose models the distribution of timefrequency coefficients (here MDCT coefficients) in terms of latent variables. The latter control the covariance matrices of Gaussian vectors of fixed-time vectors of multi-channel, multi-frequency, MDCT coefficients. In the framework of application to EEG signals, the latent variables describe some hidden mental state of the subject. This model may be seen as a way to introduce nonstationarity in approaches such as that of [?, 1], where covariance matrices are expanded as linear combinations of tensor products of frequency and channel matrices. Inference for the proposed model can be done using fairly classical algorithms, which yield estimates for the model parameters, together with maximum likelihood estimates for the sequences of latent variables.

We illustrate the model and algorithms with a case study of EEG signal analysis, namely the detection of alpha waves in rest EEG for multiple sclerosis patients and control subjects. We show that the latent variable estimation is precise enough to confirm a biological hypothesis, namely the existence of inter-hemispherical alpha wave desynchronization on multiple sclerosis subjects, which may be a potential target for diagnosis.

HARMONIC HIDDEN MARKOV MODEL

We describe in this section a stochastic model for vector-valued (i.e. multichannel) signals, based on a local Fourier basis expansion with random coefficients, and discuss the corresponding inference problem. Throughout this paper, we denote by B = {u t f , (t, f ) ∈ Λ} a frame of time-frequency atoms, to be used to expand signals. We shall limit ourselves here to finite-dimensional local Fourier bases (also called MDCT bases, see e.g. [?] for a summary), of the form

u t f [n] = 2 L w[n -tL] cos π( f + 1/2) n -tL L , (1)
t = 0, . . . N t -1 (resp. f = 0, . . . N f -1) denoting the time (resp. frequency) index. Here w is a window function, subject to compatibility relations that ensure the family to be an orthonormal basis of the signal space H = R N t N f . Notice that the model is straightforwardly extended to Gabor frames (in which case the estimation algorithms have to be modified though).

The model

Consider a time-frequency basis B = {u t f }, for example a MDCT basis, and N c signals (representing N c sensors's observations) of the form

x c = ∑ t, f a c t f u t f (2)
The (vector) synthesis coefficients a c t f are modeled using hidden Markov models (see [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF] for a review), as dependent Gaussian random vectors, whose distribution is governed by a time dependent latent state X t ∈ {1, . . . N s }. More precisely, set

A t = {(a c t f , f = 0, . . . N f -1), c = 1, . . . N c }.
Then we assume that • Conditional to X = {x t , t = 0, . . . N t -1} the vectors A s are mutually independent zero mean Gaussian vectors, with density

ψ A t (a|X t = k) = N (a; 0, Σ k )
• The latent states (which depend on the time index only, and are common to all channels) form a Markov chain, characterized by its transition matrix π. In addition, for avoiding working with large dimensional problems, we shall be particularly concerned with models in which the covariance matrix takes the form of a Kronecker product

Σ = Σ (c) ⊗ Σ ( f ) , Σ (c) (resp. Σ ( f )
) representing channels (resp. frequencies) covariances. More precisely, we shall assume that

E a c f .a c ′ f ′ = Σ (c) cc ′ Σ ( f ) f f ′ ,
where the time index has been dropped for the sake of simplicity. This is very much in the spirit of [?, 1]. Notice that such a Kronecker product form is not unique, i.e. given Σ k , Σ

k and Σ ( f ) k are defined up to a multiplicative factor. This fact has to be taken care of in the estimation procedure.

Inference

EM algorithm

Since we are using here an orthonormal basis, the inference problem can be solved using classical tools, after numerical computation of the coefficients a c t f = x c , u t f using the freely available matlab toolbox LTFAT [START_REF] Soendergaard | LTFAT, the linear time-frequency analysis toolbox[END_REF].

EM is an efficient iterative procedure to compute the Maximum Likelihood (ML) estimate in the presence of hidden data. In ML estimation, we wish to estimate the model parameters for which the observed data are the most likely. For the sake of completeness we describe here the main points of the algorithm. Each iteration EM involves two processes: The E-step (Expectation), and the M-step (Maximization). In the first one, the log-likelihood of the observed data is estimated given the current estimate of the model parameters, using the forward and backward variables defined here under. In the M-step, the likelihood function maximized leading to the so-called Baum-Welch re-estimation formulas. Convergence is ensured since the algorithm is guaranteed to increase the likelihood at each iteration.

Expectation: For t ∈ 0, .., N t -1, we define the following quantities : the forward variables α t = (α k t ) k=1,..,N s are the normalized distribution of the latent state X t conditional to the observed coefficients

(A r ) r=0,...t : α k t = P(X t = k/(A r ) r=0,...t ) × L t
where L t is the likelihood of the observations until time t.

The backward variable β k t is the likelihood of the observations (A t+1 , A t+2 , ..., A N t -1 ) conditional to X t = k for k = 1, . . . N s and t = 0, . . . N t -1.

The forward variables may be obtained recursively for t = 1, . . . N t -1 while the backward ones may be achieved by recursion too, operating on decreasing indices t = N t -1 down to 0. Thanks to these recursive equations, the complexity of the computation is linear in Given the forward and backward variables, the distribution of the transition (X t , X t+1 ) for t = 0, . . . N tconditional to the observations up to final time N t -1 reads

P(X t = k, X t+1 = l/(A t ) t=0:N t -1 ) = 1 L N t -1 α k t π k,l ψ(A t+1 )β l t+1 (3)
where

L N t -1 = N s ∑ k=1 α k t β k t
for every time t under consideration. As a consequence, the distribution of any hidden state X t , 0 ≤ t ≤ N t -2 conditional to the observations up to final time N t -1 satisfies

P(X t = k/(A r ) r=0:N t -1 ) = 1 L N t -1 α k t β k t (4)
In order to avoid underflow potentially caused by very small values of the probabilities, we used normalized versions of α and β defined by Rabiner in [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF], to which we refer for more details.

Maximization: the maximum likelihood estimates for the model parameters read:

ν k = α k 0 β k 0 L N t -1 (5) π k,l = π k,l ∑ N t -2 t=0 α k t × ψ l (A t+1 ) × β l t+1 ∑ N t -2 t=0 α k t β k t (6) Σ k = ∑ N t -1 t=0 α k t β k t A t t A t ∑ N t -1 t=0 α k t β k t (7)

Estimation for tensor product covariance matrices.

In order to control the curse of dimensionality when all sensors are kept (from 32 to 64 sensors are usually used in EEG), we also consider models in which covariance matrices take the form of Kröneker products of frequency and channel matrices.

Σ k = Σ (c) k ⊗ Σ ( f ) k , k = 1, .., N s
This approach can be seen as a refinement of the general decomposition of the covariance into a sum of Kronecker products introduced by Bijma and De Munck in [?, 1]. The above approaches can be in order to account for that particular structure. Expectation: the E-step remains unchanged but the computation of the density simplifies, owing to the fact that

Σ -1 k = (Σ (c) k ) -1 ⊗ (Σ ( f ) k ) -1 det(Σ k ) = det(Σ (c) k ) N f × det(Σ ( f ) k ) N C
Maximization: the re-estimation of covariance matrices in M-step is achieved using an alternating minimization procedure : consider

A c t = (a c t, f , f = 0, . . . N f -1) and A f t = (a c t, f , c = 1, . . . N c
) and iterate the following steps

• Estimation of Σ (c) given Σ ( f ) k : let Σ ( f ) k = ( t B k B k ) -1 denote a Cholesky decomposition of Σ ( f ) k -1 , define M k t (c, c ′ ) = B k A c t , B k A c ′ t and set Σ (c) k = 1 N f ∑ N t -1 t=0 P(X t = k)M k t ∑ N t -1 t=0 P(X t = k) (8)
• Normalization: set

Σ (c) k = Σ (c) k / Σ (c) k F , (9)
• F denoting the Frobenius norm.

• Estimation of Σ ( f ) given Σ (c) k : Using the Cholesky decomposi- tion Σ (c) k = ( t D k D k ) -1 , define P k t ( f , f ′ ) = D k A f t , D k A f ′ t
and set

Σ ( f ) k = 1 N c ∑ N t -1 t=0 P(X t = k)P k t ∑ N t -1 t=0 P(X t = k) (10)
These estimators are obtained by alternate optimization of the loglikelihood with respect to

Σ (c) k and Σ ( f )
k respectively. The normalization step described here above enables us to raise an indetermination since

Σ ( f ) k and Σ ( f )
k are defined up to a constant. Strictly speaking, this corresponds to a GEM algorithm for which the convergence to local minimum is still ensured by the fact that each step increases the likelihood of the observations with respect to the model.

Remark 1

As an alternative, the full covariance matrix Σ can also be estimated classically, the channel and frequency covariances being estimated afterwards by factorization, i.e. by minimizing a mean square error

Σ k -Σ (c) k ⊗ Σ ( f ) k 2 F under the constraint Σ (c) k F = 1.

Estimation of the hidden states sequence

Giving the estimated MAP parameters, we use the Viterbi dynamic programming algorithm to find the most likely sequence of hidden states

(X MAP 0 , ..., X MAP N t -1 ) = arg max l=(l 0 ,..l Nt -1 )
P(X = l|A t , t = 0, .., N t -1) (11)

We used the algorithm detailed in [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]. This procedure is also linear in t but requires numerical evaluation of probability densities, i.e. inversion of covariance matrices. Therefore, when N f and N c are large, using tensor product decompositions yields more stability in the numerical evaluation.

Numerical simulations

To validate the algorithms, we simulated signals (N c = 8 channels of length 100000 samples each, corresponding to 40 sec of signal) following the above model. For the sake of simplicity, and because of the nature of the real EEG signals we focus on in section 3, we limite ourselves to N s = 2 hidden states. MDCT transform was tuned so as to generate N f = 10 frequency bins for each channel. The frequency covariance matrices were set to the estimated values for the EEG case study presented below (see Fig. 4 in section 3.2) and artificial channel covariance matrices were generated so as to exhibit topographically distinct sources for the two states (see Fig. 1). A reference transition probability p 11 = P(X t = 1|X t-1 = 1) = 0.85 was chosen (corresponding to a typical value for alpha waves duration, i.e. 700ms), simulations were run for various values of p 00 = P(X t = 0|X t-1 = 0) in the interval [0.75; 0.95] (100 runs per value of p 00 ), and relative errors Σ-Σ F Σ F were computed. The results, displayed in Fig. 2, clearly show that the Kronecker product based estimate outperforms the classical one in terms of accuracy (for comparable running time). As expected, when p 00 grows, the average length of zero-state regions increases, and the accuracy of the estimate of Σ 0 (resp. Σ 1 ) grows (resp. decays). The hidden state identification error rates are very comparable for the two methods and always remain at a very low level (less than 0.5%).

APPLICATION TO EEG SIGNALS

The above model, inspired by [START_REF] Molla | An hybrid audio scheme using hidden Markov models of waveforms[END_REF], was motivated by applications to the modeling of EEG signals, and we describe now a case study, devoted to detection of alpha waves in rest EEG signal 1 . Alpha waves are short duration time-localized oscillations (with frequency around 10 Hz) which appear in EEG signals in some specific situations, and can naturally be accounted for by timefrequency representations. They are also topographically localized in specific sensors situated in posterior regions of head. The alpha wave occurrence may be considered a non-stationary effect, i.e. a departure from a stationary background signal. This therefore motivates the use of hidden Markov models as described above.

Problem statement and data

aim at automatic detection of alpha waves, i.e. segmentation of the signal into time epochs where alpha waves are present/absent. Such a task could be achieved by band pass filtering of signals, followed by appropriate thresholding, or time-frequency transform thresholding. In both cases this would however raise the question of the determination of the threshold. The HHMM described above provides a more systematic approach.

The main goal of the current study is to use the model in an unsupervised manner, in order to test its ability to blindly detect alpha waves. The case study below provides a positive answer. In other situations though (i.e. with more complex signals), supervised training may be necessary in order to obtain accurate detection.

For the problem under consideraton, we use the model with two instances for the hidden states: alpha (k = 1)/non-alpha (k = 0) parameters are the initial hidden state distribution ν where ν k = P(X 0 = k)

• the stochastic transition matrix of the Markov Chain Π where

Π k,l = P(X t+1 = l|X t = k) • the covariance matrices Σ ( f ) 0 , Σ ( f ) 1 , Σ (c) 0 , Σ (c)
1 . Once the estimation of these parameters is achieved using the EM algorithm described here above, we determine the most likely hidden states sequence. This sequence is the main ingredient in the case study below since it detects the bursts in the time-frequency images which corresponds to alpha waves. Indeed, it is clearly seen in Fig. 3 that values equal to one correspond to regions in the timefrequency domain with significantly larger MDCT coefficients in the alpha band. The approach was applied to EEG data originating from the CODYSEP dataset, designed to study the impact of sclerosis in inter-hemispherical transfer. The dataset consists in 31 ill subjects (hereafter termed SEP) and 20 controls (TEM); 17 channels EEG signals were collected at a 256 Hz sampling rate.

A subset of 14 ill patients and 16 controls was selected, namely those exhibiting sufficiently similar time-frequency contents, in particular in the 8-12 Hz range. For the sake of precision in the time domain, the hop size N f was set to N f = 8 (i.e. a time resolution of 125 ms), resulting in a moderate frequency resolution (each MDCT basis function having a bandwidth of approximately 4Hz). This choice is certainly questionable, and the selection procedure should probably turned into a more sustematic one. Our choice was motivated by the desire of focusing on th same phenomenon (i.e. wiht the same time-frequency patterns) for all retained subjects.

Estimation results

A subset of sensors was also selected that correspond to relevant scalp locations for observing alpha waves. Corresponding MDCT transforms were computed, and the model parameters (for a twostate model: alpha and non alpha states) were estimated using the algorithms above. Examples of covariance matrices for the two states are displayed in Fig. 4 and5. The frequency covariance matrices shown in Fig. 4 are almost which indicates decorrelation of the frequency bands. The main difference between the two matrices appears for the frequency band 8-12Hz (element (3,3)), which is significantly bigger in 'alpha' state. This state therefore succeds to capture alpha waves. The channel covariance matrices in Fig 5 are also significantly different for some channels: in particular, channels O1 and O2 (matrix [START_REF] Flandrin | Temps-Fréquence. Traité des Nouvelles Technologies[END_REF][START_REF] Flandrin | Temps-Fréquence. Traité des Nouvelles Technologies[END_REF] and (11,11)) are overactivated in 'alpha' state. Again this is coherent, since those sensors correspond to regions where alpha wave signals are most visible. This phenomenon is even more obvious on the graphical representation of topographies (Fig. 6),actually representing the diagonals of matrices of Fig 5 on which the relevant sensors O1 and are clearly more "energetic" in 'alpha' state.

After completion of parameter estimation, a maximum likelihood estimate for the sequence of latent states is obtained via the Viterbi algorithm.

A study of alpha wave desynchronization

Let us start by quoting [START_REF] Cover | Multiple sclerosis patients show a highly significant decrease in alpha band interhemispheric synchronization measured using MEG[END_REF]: "in multiple sclerosis, both axonal damage and demyelination occur. Therefore it can be expected that the connectivity between the different regions in the brains of MS patients will be impaired compared to healthy controls." Losses in synchronization between EEG signals in left and right hemispheres have been reported in the literature, where coherence is generally used to detect the phenomenon. The HHMM method described here is an alternative natural candidate for testing and evaluating it.

To test the desynchronization assumption, hidden states sequences were estimated for all subjects in the dataset, separately for left sensors and rignt sensors. The latter sequences take the form of sequences of zeros (non-alpha state) and ones (alpha state). For each subject, the Hamming distance between the left sensors X (L) t and right sensors X

(R) t was computed D = ∑ t X (L) t -X (R) t
and the empirical distributions of the so-obtained distances were compared using a non parametric test. The figure 7 shows the boxplots of the two families (SEP and TEM) of distances, and seems to indicate a significant difference: the SEP data exhibit a larger leftright assymetry, in accordance with the above mentioned desynchronization hypothesis. To assess statistically the significance of the result, we performed a Mann-Whitney test. The corresponding P-value was found to equal P = 0.0384, which confirms quantitatively the hypothesis of two distinct distributions.

CONCLUSIONS AND PERSPECTIVES

We have described in this paper a new stochastic signal model, based on hidden Markov modelling of vectors of MDCT coefficients of multichannel signals. We have described corresponding estimation algorithms, and an adaptation to the case where covariance matrices take the form of Kronecker products of smaller matrices. Finally, we have tested the algorithms on a case study involving EEG signals.

While the algorithms described here are not new, the model itself and its adaptation to EEG type data is original. Our numerical results on simulations show that the estimators are good enough to handle data with realistic characteristics (length, number of channels,...).

In addition, the case study clearly shows the relevance of the model for real data processing, since it yields a biologically relevant result. Further work will concern the validation of this approach to different types of EEG detection and characterization problems. Applications to BCI (brain computer interfaces) problems are currently under study. Such applications will probably require online versions of the algorithms.

From a more fundamental point of view, it is to be noted that MDCT bases are known for having important limitations in terms of time-frequency resolution. It is therefore desirable to turn to better quality time-frequency expansions, such as Gabor frames. However, the loss of orthogonality in the waveform system makes the estimation more complex, as it will not be possible any more to use the expansion as an orthogonal transform and develop standard algorithms on the expansion coefficients. Dedicated algorithms have to be designed, which is work in progress.

Figure 1 :

 1 Figure 1: Channel covariance matrices for state 'non-alpha' (left) and 'alpha' (right)

Figure 2 :

 2 Figure 2: Evolution of the estimation relative errors as a function of the transition probability p 11

Figure 3 :

 3 Figure 3: MDCT coefficients for a fixed sensor (top) and hidden states sequence estimated (from all sensors) using Viterbi algorithm (bottom)

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Frequency covariance matrices for state 'non-alpha' (left) and 'alpha' (right)

Figure 7 :

 7 Figure 7: Boxplots of the Hamming distances between left and right hidden states for SEP (left) and TEM (right)

Let us stress however that a similar approach can be adapted to much more general situations, in which signals exhibit sufficient time-frequency localized components.

Acknowledgements

This work was partly supported by CNRS, through the NeuroInformatique program's project Codysep, which allowed the recording of the data used in this work, and by the ANR project CO-ADAPT (ANR-09-DEFIS-002). We wish to thank the anonymous reviewers for their very relevant comments, which helped us improving the paper.


Sylvain Nascimbène
email: sylvain.nascimbene@lkb.ens.fr

Nir Navon

Frédéric Chevy

Christophe Salomon

The equation of state of ultracold Bose and Fermi gases: a few examples

, that reveals a significant difference in the high-temperature regime. The normal phase, at low temperature, is well described by a Landau Fermi liquid model, and we observe a clear thermodynamic signature of the superfluid transition. In a second part we apply the same procedure to Bose gases. From a single image of a quasi ideal Bose gas we determine the equation of state from the classical to the condensed regime. Finally the method is applied to a Bose gas in a 3D optical lattice in the Mott insulator regime. Our equation of state directly reveals the Mott insulator behavior and is suited to investigate finite-temperature effects.

Introduction

Ultracold gases are a privileged tool for the simulation in the laboratory of model Hamiltonians relevant in the fields of condensed matter, astrophysics or nuclear physics [START_REF] Bloch | Many-body physics with ultracold gases[END_REF]. As an example, thanks to the short-range character of interactions, ultracold Fermi mixtures prepared around a Feshbach resonance mimic the behavior of neutron matter in the outer crust of neutron stars [START_REF] Bertsch | Many-Body Challenge Problem[END_REF][START_REF] Gezerlis | Strongly paired fermions: Cold atoms and neutron matter[END_REF]. For cold atoms, the density inhomogeneity induced by the trapping potential has long made the connection between the Hamiltonian of a homogeneous system and an ultracold gas indirect. Early experimental thermodynamic studies have provided global quantities averaged over the whole trapped gas, such as total energy and entropy [START_REF] Stewart | Potential Energy of a 40 K Fermi Gas in the BCS-BEC Crossover[END_REF][START_REF] Luo | Measurement of the entropy and critical temperature of a strongly interacting Fermi gas[END_REF], collective mode frequencies [START_REF] Altmeyer | Precision measurements of collective oscillations in the BEC-BCS crossover[END_REF], or radii of the different phases that may be observed in an imbalanced Fermi gas [START_REF] Partridge | Pairing and phase separation in a polarized Fermi gas[END_REF][START_REF] Zwierlein | Direct observation of the superfluid phase transition in ultracold Fermi gases[END_REF][START_REF] Nascimbène | Collective Oscillations of an Imbalanced Fermi Gas: Axial Compression Modes and Polaron Effective Mass[END_REF]. Reconstructing the equation of state of the homogeneous gas then requires to deconvolve the effect of the trapping potential, a delicate procedure that has not been done so far. However, the gas can often be considered as locally homogeneous (local density approximation (LDA)), and careful analysis of in situ density profiles can directly provide the equation of state of a homogeneous gas [START_REF] Shin | Determination of the equation of state of a polarized fermi gas at unitarity[END_REF][START_REF] Ho | Obtaining the phase diagram and thermodynamic quantities of bulk systems from the densities of trapped gases[END_REF][START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Navon | The Equation of State of a Low-Temperature Fermi Gas with Tunable Interactions[END_REF]. In the case of two-dimensional gases, in situ images taken along the direction of tight confinement obviously give access to the surface density [START_REF] Hadzibabic | Berezinskii-Kosterlitz-Thouless crossover in a trapped atomic gas[END_REF][START_REF] Cladé | Observation of a 2D Bose Gas: From Thermal to Quasicondensate to Superfluid[END_REF][START_REF] Gemelke | In situ observation incompressible Mottinsulating domains in ultracold atomic gases[END_REF][START_REF] Bakr | A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice[END_REF] and thus to the equation of state [START_REF] Rath | The equilibrium state of a trapped two-dimensional Bose gas[END_REF]. For threedimensional gases, imaging leads to an unavoidable integration along the line of sight. As a consequence, inferring local quantities is not straightforward. Local density profiles can be computed from a cloud image using an inverse Abel transform for radially symmetric traps [START_REF] Shin | Phase diagram of a two-component Fermi gas with resonant interactions[END_REF]. A more powerful method was suggested in [START_REF] Ho | Obtaining the phase diagram and thermodynamic quantities of bulk systems from the densities of trapped gases[END_REF] and implemented in [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Navon | The Equation of State of a Low-Temperature Fermi Gas with Tunable Interactions[END_REF]: as explained below, for a harmonically trapped gas the local pressure is simply proportional to the integrated in situ absorption profile. Using this method, the low-temperature superfluid equation of state for balanced and imbalanced Fermi gases have been studied as a function of interaction strength [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Navon | The Equation of State of a Low-Temperature Fermi Gas with Tunable Interactions[END_REF]. In this paper we describe in more detail the procedure used to determine the equation of state of a spin-unpolarized Fermi gas in the unitary limit [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF]. We compare our data with recent results from the Tokyo group [START_REF] Horikoshi | Measurement of Universal Thermodynamic Functions for a Unitary Fermi Gas[END_REF], and reveal a significant discrepancy in the high-temperature regime. In a second part we apply the method to ultracold Bose gases. From an in situ image of 7 Li, we obtain the equation of state of a weakly-interacting Bose gas. Finally, analyzing experimental profiles of a Bose gas in a deep optical lattice [START_REF] Fölling | Formation of spatial shell structure in the superfluid to Mott insulator transition[END_REF], we observe clear thermodynamic signatures of the Mott insulator phases.

Measurement of the local pressure inside a trapped gas

In the grand-canonical ensemble, all thermodynamic quantities of a macroscopic system can be derived from the equation of state P = f (µ, T ) relating the pressure P to the chemical potential µ and temperature T . P can be straigthforwardly deduced from integrated in situ images.

Consider first a single-species ultracold gas, held in a cylindrically symmetric harmonic trap whose frequencies are labeled ω x = ω y ≡ ω r in the transverse direction, and ω z in the axial direction. Provided that the local density approximation is satisfied, the gas pressure along the z axis is given by [START_REF] Ho | Obtaining the phase diagram and thermodynamic quantities of bulk systems from the densities of trapped gases[END_REF]:

P (µ z , T ) = mω 2 r 2π n(z), (1)
where n(z) = dx dy n(x, y, z) is the doubly-integrated density profile, µ z = µ 0 -1 2 mω 2 z z 2 is the local chemical potential on the z axis, µ 0 is the global chemical potential. n(z) is obtained from an in situ image taken along y, by integrating the optical density along x (see Fig. 1). As described below, if one independently determines the temperature T and chemical potential µ 0 , then each pixel row of the absorption image at a given position z, provides an experimental data point for the grand-canonical equation of state P (µ z , T ) of the homogeneous gas. The large amount of data obtained from several images allows one to perform an efficient averaging, leading to a low-noise equation of state.

This formula is also valid in the case of a two-component Fermi gas with equal spin populations if n(z) is the total integrated density. The method can be generalized to multicomponent Bose and Fermi gases, as first demonstrated on spin-imbalanced Fermi gases in [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Navon | The Equation of State of a Low-Temperature Fermi Gas with Tunable Interactions[END_REF].

Thermodynamics of a Fermi gas with resonant interactions

In this section we describe the procedure used in [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF] to determine the grand-canonical equation of state of a homogeneous and unpolarized Fermi gas with resonant interactions (a = ∞). We also compare our data with recent measurements from the Tokyo group [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Horikoshi | Measurement of Universal Thermodynamic Functions for a Unitary Fermi Gas[END_REF]. We then study the physical content of the equation of state at low temperature.

Grand-canonical equation of state

In the grand-canonical ensemble, the equation of state of a spin-unpolarized Fermi gas in the unitary limit, can be written as

P (µ, T ) = P (0) (µ, T )h T (ζ), (2)
where P (0) (µ, T ) is the pressure of a non-interacting two-component Fermi gas and ζ = exp(µ/k B T ) is the inverse fugacity. Since P (0) (µ, T ) is known, the function h T (ζ) completely determines the equation of state P (µ, T ). Let us now describe the procedure used to measure it. The pressure profile of the trapped gas along the z axis is directly obtained from its in situ image using equation [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF]. One still has to know the value of the temperature T and global chemical potential µ 0 in order to infer h T (ζ). We use a small amount of 7 Li atoms, at thermal equilibrium with the 6 Li component, as a thermometer.

We then extract µ 0 from the pressure profile, by comparison in the cloud's wings with a reference equation of state. For high-temperature clouds (k B T > µ 0 ), we choose µ 0 so that the wings of the pressure profile match the second-order virial expansion [START_REF] Ho | High temperature expansion applied to fermions near Feshbach resonance[END_REF] (see Fig. 2a):

P (µ, T ) = 2k B T λ 3 dB (T ) e µ/k B T + 4 3 √ 2 e 2µ/k B T + . . . . (3)
For colder clouds, the signal-to-noise ratio is not good enough, in the region where (3) is valid, to extract µ 0 using the same procedure. We thus rather use the equation of state determined from all images previously treated as a reference, since it is accurate on a wider parameter range than (3) (see Fig. 2b). We then iterate this procedure at lower and lower temperatures, eventually below the superfluid transition. By gathering the data from all images and statistical averaging, we obtain a low-noise equation of state in the range 0.02 < ζ < 5 (see Fig. 3a).

Canonical equation of state

In [START_REF] Horikoshi | Measurement of Universal Thermodynamic Functions for a Unitary Fermi Gas[END_REF] a canonical equation of state E(n, T ) expressing the energy E as a function of density and temperature was measured using fits of absorption images taken after a

/2k B T λ -3 dB versus -µ/k B T = V (z)/k B T -µ 0 /k B T (black points).
A wrong choice of µ 0 corresponds in this representation to a translation of the data in abscissa. We adjust µ 0 so that the wings of the pressure profile match a reference equation of state (in red). (a) For high-temperature clouds, we use the second-order virial expansion (3) (b) For a lower temperature pressure profile, we minimize its distance with the averaged equation of state deduced from higher temperature images (in red) in the overlap region.

short time-of-flight. In situ density profiles were deduced by assuming a hydrodynamic expansion. The temperature was extracted from the cloud's total potential energy at unitarity, using the experimental calibration made in [START_REF] Luo | Measurement of the entropy and critical temperature of a strongly interacting Fermi gas[END_REF]. In Fig. 3b the data from [START_REF] Horikoshi | Measurement of Universal Thermodynamic Functions for a Unitary Fermi Gas[END_REF] is plotted as E(n, T )/E (0) (n, T ) as a function of θ = T /T F , where n is the total atom density, T F is the Fermi temperature, and E (0) (n, T ) is the energy of a non-interacting Fermi mixture. The comparison between the two equations of state requires to express our data in the canonical ensemble. The density n = ∂P/∂µ| T is calculated by taking a discrete derivative, and we obtain the black points in Fig. 3b. While the two sets of data are in satisfactory agreement in the low-temperature regime T /T F < 0.4, they clearly differ in the high-temperature regime. The disagreement of the data from [START_REF] Horikoshi | Measurement of Universal Thermodynamic Functions for a Unitary Fermi Gas[END_REF] with the secondand third-order virial expansions calculated in [START_REF] Ho | High temperature expansion applied to fermions near Feshbach resonance[END_REF][START_REF] Liu | Virial expansion for a strongly correlated Fermi gas[END_REF] indicates a systematic error in this regime. This is possibly due to a breakdown of hydrodynamics during the time-of-flight as expected at high temperature.

Ideal Gas Ideal Gas

Fermi liquid behavior in the normal phase

Above the superfluid transition and in the low-temperature regime 0.05 < ζ < 0.5, our data is well modelled by a Fermi liquid equation of state

P FL (µ, T ) = 2 15π 2 2m 2 3/2 µ 5/2 ξ -3/2 n + 5π 2 8 ξ -1/2 n m * m k B T µ 2 , (4)
where ξ n = 0.51(1) and m * = 1.12(3)m respectively characterize the compressibility of the normal phase extrapolated to zero temperature and the effective mass of the lowlying excitations. The agreement with ( 4) is better than 5% in a large parameter range 0.33 µ < k B T < 2 µ. Our value of ξ n is in agreement with the variational Fixed-Node Monte-Carlo calculations ξ n = 0.54 in [START_REF] Carlson | Superfluid Fermi gases with large scattering length[END_REF], ξ n = 0.56 in [START_REF] Lobo | Normal state of a polarized Fermi gas at unitarity[END_REF], and with the Quantum Monte-Carlo calculation ξ n = 0.52 in [START_REF] Bulgac | Quantum Monte Carlo simulations of the BCS-BEC crossover at finite temperature[END_REF]. It is surprising that the quasiparticle mass m * is quite close to the free fermion mass, despite the strongly-interacting regime. Note also that this mass is close to the effective mass m * = 1.20 m of a single spin-down atom immersed into a Fermi sea of spin-up particles (the Fermi polaron) [START_REF] Chevy | Universal phase diagram of a strongly interacting fermi gas with unbalanced spin populations[END_REF][START_REF] Lobo | Normal state of a polarized Fermi gas at unitarity[END_REF][START_REF] Combescot | Normal state of highly polarized Fermi gases: simple many-body approaches[END_REF][START_REF] Prokof | Fermi-polaron problem: Diagrammatic monte carlo method for divergent sign-alternating series[END_REF][START_REF] Combescot | Normal state of highly polarized Fermi gases: full many-body treatment[END_REF][START_REF] Shin | Determination of the equation of state of a polarized fermi gas at unitarity[END_REF][START_REF] Nascimbène | Collective Oscillations of an Imbalanced Fermi Gas: Axial Compression Modes and Polaron Effective Mass[END_REF][START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF].

Superfluid transition

The deviation of the experimental data from (4) for ζ < 0.05 signals the superfluid phase transition. This transition belongs to the U(1) universality class, and the critical region is expected to be wide [START_REF] Taylor | Critical behavior in trapped strongly interacting Fermi gases[END_REF] in the unitary limit. Assuming that our low-temperature data belongs to the critical region, we fit our data with a function

P (µ, T ) = P FL (µ, T ) + A(ζ c -ζ) 2-α H(ζ c -ζ), (5)
where H is the Heaviside function and α ≃ -0.013 is the specific heat critical exponent, measured with a very good accuracy on liquid 4 He [START_REF] Lipa | Very high-resolution heat-capacity measurements near the lambda point of helium[END_REF]. We obtain the position of the superfluid transition ζ c = 0.05, or k B T c /µ = 0.33, in agreement with the value k B T c /µ = 0.32(3) extracted in [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF] using a simpler fit function. We thus confirm more rigorously our previous determination of the superfluid transition. In the appendix we discuss the validity of local density approximation around the superfluid transition. Under our current experimental conditions, the deviation from LDA is very small.

Thermodynamics of a weakly-interacting Bose gas

In this section we apply equation (1) to the case of trapped Bose gases. First we test the method by determining the equation of state of a weakly-interacting Bose gas [START_REF] Caracanhas | Finite temperature correction to the Thomas-Fermi approximation[END_REF][START_REF] Romero-Rochín | Equation of State of an Interacting Bose Gas Confined by a Harmonic Trap: The Role of the Harmonic Pressure[END_REF].

We use an in situ absorption image of a 7 Li gas taken from [START_REF] Schreck | Quasipure Bose-Einstein condensate immersed in a Fermi sea[END_REF] (see Fig. 4a). 7 Li atoms are polarized in the internal state |F = 1, m F = -1 , and held in an Ioffe-Pritchard magnetic trap with ω r /2π = 4970 Hz and ω z /2π = 83 Hz, in a bias field B 0 ≃ 2 G. Thermometry is provided by a gas of 6 Li atoms, prepared in |F = 1 2 , m F = -1 2 , and in thermal equilibrium with the 7 Li cloud.

Determination of the equation of state

The equation of state of a weakly-interacting Bose gas can be expressed, in the grandcanonical ensemble, as:

P (µ, T ) = k B T λ 3 dB (T ) g(ζ),
where ζ = e -µ/k B T is the inverse fugacity and λ dB (T ) = 2π 2 /mk B T is the thermal de Broglie wavelength. The pressure profile is calculated using [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF]. We aim here at measuring g(ζ). We obtain the global chemical potential value µ 0 = 0.10 k B T by fitting the 7 Li profile in the non-condensed region |z| > 50 µm with a Bose function:

P (µ z , T ) = k B T λ 3 dB (T ) g 5/2 (ζ z ), ζ z = e -µ 0 /k B T exp mω 2 z z 2 2k B T , g 5/2 (z) = ∞ k=1 z -k k 5/2 .
Combining the measurement of the pressure profile, of the cloud's temperature T and global chemical potential µ 0 , we obtain the thermodynamic function g(ζ) plotted in Fig. 4b.

Analysis of the equation of state

In the region ζ > 1 the data agrees with the Bose function g(ζ) = g 5/2 (ζ) expected for a weakly-interacting Bose gas. The departure from the thermodynamic function of a classical gas g(ζ) = ζ -1 , and especially the fact that g(ζ) > 1 above the condensation threshold, is a thermodynamic signature of a bosonic bunching effect, as observed in [START_REF] Yasuda | Observation of two-atom correlation of an ultracold neon atomic beam[END_REF][START_REF] Fölling | Spatial quantum noise interferometry in expanding ultracold atom clouds[END_REF][START_REF] Schellekens | Hanbury Brown Twiss Effect for Ultracold Quantum Gases[END_REF]. The sudden and fast increase of our data for ζ 1 indicates the Bose-Einstein condensation threshold. In the local density approximation framework, the chemical potential of a weakly-interacting Bose-Einstein condensate reads:

µ = 4π 2 a 77 m 7 n,
where m 7 is the 7 Li atom mass and a 77 is the scattering length describing s-wave interactions between 7 Li atoms. We neglect here thermal excitations in the condensed region.

g(ζ) = g 5/2 (ζ c ) + λ dB (T ) 4 a 77 (log 2 ζ -log 2 ζ c ). (6)
Fitting our data with the function g(ζ) given by ( 6 The uncertainties take into account the fit uncertainty and the uncertainty related to the temperature determination. The condensation threshold is in agreement with the value ζ c = 1 expected for an ideal Bose gas, the mean-field correction being on the order of 1% [START_REF] Giorgini | Condensate fraction and critical temperature of a trapped interacting Bose gas[END_REF][START_REF] Giorgini | Thermodynamics of a trapped Bose-condensed gas[END_REF]. Our measurement of the scattering length is in agreement with the most recent calculations a 77 = 7(1) a 0 [41].

Extending this type of measurement to larger interaction strengths on Bose gases prepared close to a Feshbach resonance would reveal more complex beyond-mean-field phenomena, provided thermal equilibrium is reached for strong enough interactions.

Mott-insulator behavior of a Bose gas in a deep optical lattice

Here we extend our grand-canonical analysis to the case of a 87 Rb gas in an optical lattice in the Mott insulator regime. By comparing experimental data with advanced Monte Carlo techniques, it has been shown that in many circumstances the local density approximation is satisfied in such a system [START_REF] Trotzky | Suppression of the critical temperature for superfluidity near the Mott transition: validating a quantum simulator[END_REF]. We analyze the integrated density profiles of the Munich group, Fig. 2 of [START_REF] Fölling | Formation of spatial shell structure in the superfluid to Mott insulator transition[END_REF].

Realization of the Bose-Hubbard model with ultracold gases

Atoms are held a trap consisting of the sum of a harmonic potential V h (x, y, z) and a periodic potential V 0 sin 2 (kx) + sin 2 (ky) + sin 2 (kz) , created by three orthogonal standing waves of red-detuned laser light at the wavelength λ = 2π/k = 843 nm. The atoms occupy the lowest Bloch band and realize the Bose-Hubbard model [START_REF] Jaksch | Cold bosonic atoms in optical lattices[END_REF]:

Ĥ = -J i,j â † i âj + U 2 i (â † i âi -1)â † i âi , (7)
with a local chemical potential µ(r) = µ 0 -V h (r). The index i refers to a potential well at position r i , J is the tunneling amplitude between nearest neighbors, and U is the on-site interaction, U and J being a function of the lattice depth [START_REF] Bloch | Many-body physics with ultracold gases[END_REF]. The slow variation of V h (r) compared with the lattice period λ/2 justifies the use of local density approximation.

We consider here the case of a large lattice depth V 0 = 22E r , for which J ≃ 0.003 U ∼ 0, and assume that the temperature is much smaller than U. In this regime the gas is expected to form a Mott insulator: in the interval µ ∈ [(p -1)U, pU], where p is an integer, the atom number per site remains equal to p, and the density is equal to n = p(2/λ) 3 . Integrating Gibbs-Duhem relation between 0 and µ, we obtain that the pressure P is a piecewise linear function of µ:

P (µ, T = 0) = 2 λ 3 µ - p -1 2
U p where (p -1)U < µ < pU.

Determination of the equation of state

We use a series of three images from [START_REF] Fölling | Formation of spatial shell structure in the superfluid to Mott insulator transition[END_REF], labeled a, b and c, with different atom numbers N a = 1.0 × 10 5 , N b = 2.0 × 10 5 and N c = 3.5 × 10 5 (see Fig. 5a). The integrated profiles n(z) are not obtained using in situ absorption imaging but rather using a tomographic technique, providing a ∼ 1 µm resolution. The pressure profile is then obtained using equation [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF].

Each image i = a, b, c plotted as P as a function of -1 2 mω 2 z z 2 provides the equation of state P (µ) translated by the unknown global chemical potential µ 0 i . By imposing that all images correspond to the same equation of state (in the overlapping µ/U region), we deduce the chemical potential differences between the different images µ 0 bµ 0 a = 0.56 U and µ 0 cµ 0 b = 0.61 U (see Fig. 5b). Gathering the data from all images, we thus obtain a single equation of state, translated by µ 0 a which is still unknown. We fit this data with a function translated by µ 0 a from the following function, capturing the Mott insulator physics: P U(λ/2) -3 = 0 for µ < 0 =

n 1 µ U for 0 < µ < δµ 1 = n 1 δµ 1 U + n 2 µ -δµ 1 U for δµ 1 < µ < δµ 1 + δµ 2 = n 1 δµ 1 U + n 2 δµ 2 U + n 3 µ -δµ 1 -δµ 2 U for δµ 1 + δµ 2 < µ,
with µ 0 a , δµ 1 , δµ 2 , n 1 , n 2 and n 3 as free parameters. The value µ 0 a = 1.51 U yielded by the fit thus corresponds to the condition P → 0 when µ → 0. Once it is determined, we obtain the equation of state of the Bose-Hubbard model in the Mott regime, plotted in Fig. 6.

n 0 0 n 0 0 n 1 0.9 1 n 1 0.9 1 n 2 2.0 1 n 2 2.0 1 n 3 3.1 1 n 3 3.1 1 3 2 1 0 1 2 3

Observation of a Mott-insulator behavior

After fitting the value of µ 0 a , the other parameters resulting from the fit exhibit the characteristic features of incompressible Mott phases. The occupation number in the first Mott region is n 1 = 0.9(1) atom per site and the size is δµ 1 = 0.9(1)U. The second Mott region occupation number is n 2 = 2.0(1) and its size is δµ 2 = 1.1(1)U. Finally, the third Mott region occupation number is n 3 = 3.1 [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF]. These values agree with the theoretical values n i = i and δµ i = U, in the T = 0 and J = 0 limits.

Estimation of finite temperature effects

The equation of state deduced from the experimental data is also suited for investigating finite-temperature effects. Since sites are decoupled in the regime J ≪ U, k B T considered in this study, the finite-temperature equation of state is easily calculated from the thermodynamics of a single site [START_REF] Gerbier | Boson Mott insulators at finite temperatures[END_REF][START_REF] Capogrosso-Sansone | Phase diagram and thermodynamics of the three-dimensional Bose-Hubbard model[END_REF]:

P (µ, T ) = k B T (λ/2) 3 log ∞ p=0 exp - Up(p -1)/2 -µp k B T . (8)
Fitting now the experimental data with (8) and T and µ 0 a as free parameters, we deduce: k B T = 0.09 +0.04 -0.09 U. This value is in agreement with a direct fit of the density profiles and number statistics measurements [46]. This temperature is significantly smaller than the temperature k B T * ≃ 0.2 U at which the Mott insulator is expected to melt [START_REF] Gerbier | Boson Mott insulators at finite temperatures[END_REF]. Second, this temperature should be considered as an upper limit because of its uncertainty on the lowtemperature side. Indeed, the finite resolution of the images tends to smear out the sharp structure associated with Mott insulator boundaries, leading to an overestimation of the actual temperature. To overcome this limit, the spin-gradient thermometry proposed in [START_REF] Weld | Spin Gradient Thermometry for Ultracold Atoms in Optical Lattices[END_REF] could be employed.

Summary and concluding remarks

To summarize, we have shown on various examples of Fermi and Bose gas systems how in situ absorption images can provide the grand-canonical equation of state of the homogeneous gas. This equation of state is obtained up to a global shift in chemical potential and we have given several examples for its determination. The method relies on the local density approximation, which is satisfied in many situations, but notable exceptions exist such as the case of the ideal Bose gas. The equation of state given by this procedure allows direct comparison with many-body theories. While we have illustrated here this method on a single-component Bose gas and a spin-balanced Fermi gas, it can easily be generalized to multi-component gases. For instance the phase diagram and the superfluid equation of state of spin-imbalanced Fermi gases have been obtained in [START_REF] Nascimbène | Exploring the thermodynamics of a universal Fermi gas[END_REF][START_REF] Navon | The Equation of State of a Low-Temperature Fermi Gas with Tunable Interactions[END_REF]. We expect this method to be very useful in the investigation of Bose-Bose, Bose-Fermi and Fermi-Fermi mixtures. Finally the equation of state of a Bose gas close to a Feshbach resonance may reveal thermodynamic signatures of beyond-mean-field behavior in Bose-Einstein condensates [START_REF] Papp | Bragg Spectroscopy of a Strongly Interacting 85 Rb Bose-Einstein Condensate[END_REF].

Figure 1 .

 1 Figure1. Scheme of the local pressure measurement: the absorption of a probe beam propagating along the y direction provides a 2D image on the CCD camera. Integration of this image along x provides the doubly-integrated density profile n(z) and, using equation (1), the pressure profile along the z axis.

Figure 2 .

 2 Figure 2. Determination of µ 0 : we plot the data from an in situ image as P/2k B T λ -3

Figure 3 .

 3 Figure 3. (a) Grand-canonical equation of state of a two-component Fermi gas with resonant interactions from [1] (black dots). (b) Canonical equation of state from the Tokyo group[START_REF] Horikoshi | Measurement of Universal Thermodynamic Functions for a Unitary Fermi Gas[END_REF] (open circles) and from the ENS group (black dots). The dashed black line is the ideal gas equation of state, the dot-dashed (solid) black line is the second-(third-) order virial expansion, the solid green line is the Fermi liquid equation (4) and the solid blue line is the fit function[START_REF] Gezerlis | Strongly paired fermions: Cold atoms and neutron matter[END_REF] in the superfluid phase. The superfluid transition occurs at ζ = 0.05.

Figure 4 .

 4 Figure 4. (a) Integrated density profiles n(z) for the 7 Li component (blacks dots) and the 6 Li component (open circles). The solid line is a fit of the 6 Li component with a finite-temperature Thomas-Fermi profile, yielding T = 1.6(1) µK. (b) Thermodynamic function g(ζ) determined from the 7 Li profile. The solid line is a fit of the data with a Bose function in the non-condensed region and a mean-field equation of state in the condensed region (see text). The dashed line is the equation of state of a classical gas g(ζ) = ζ -1 . The difference between the dashed and solid lines around ζ = 1 is a consequence of Bose statistics.

  Gibbs-Duhem relation at fixed temperature dP = ndµ between the condensation threshold ζ c and ζ < ζ c , and imposing the continuity at ζ = ζ c , we obtain the equation of state in the condensed phase:

  ) for ζ < ζ c and with g 5/2 (ζ) for ζ > ζ c , we obtain ζ c = 1.0(1) and a 77 = 8(4)a 0 = 0.4(2) nm.

Figure 5 .

 5 Figure 5. (a) Integrated density profiles n(z) corresponding to image a (open squares) b (black dots) and c (crosses), from [21]. (b) Determination of the global chemical potential difference µ 0 c -µ 0 b by superposing the equations of states given by each image.

3 Figure 6 .

 36 Figure 6. Equation of state of a Bose gas in an optical lattice, in the Mott insulator regime. The solid line is a fit with a piecewise linear function capturing the Mott insulator behavior. The slope dP/dµ provides the density in each of the Mott zone, n 1 = 0.9(1), n 2 = 2.0(1), n 3 = 3.1(1).

Acknowledgments

We are grateful to Fabrice Gerbier and Kenneth Guenter for stimulating discussions. We acknowledge support from ERC (Ferlodim), ESF (Euroquam Fermix), ANR FABIOLA, Région Ile de France (IFRAF), and Institut Universitaire de France.

Appendix: validity of local density approximation

Let us now discuss the validity of local density approximation around the superfluid transition in our experiment. Along the z axis, the correlation length ξ diverges around the transition point z = z c according to

, where ν = 0.67 is the correlation length critical exponent, directly measured in [START_REF] Donner | Critical behavior of a trapped interacting Bose gas[END_REF], and in agreement with ν = (2α)/3. Local density approximation is expected to become inaccurate in the region z cδz < z < z c + δz where δz is given by [START_REF] Taylor | Critical behavior in trapped strongly interacting Fermi gases[END_REF][START_REF] Pollet | Criticality in Trapped Atomic Systems[END_REF]:

z c is on the order of the cloud size along z, and is much larger that k -1 F which is on the order of the inter-particle distance. Given the parameters of our experiments, (k F z c ) -1/(1+ν) ∼ 1% and the size δz where local density approximation is invalid is very small. Given the noise of our data (a few %), the deviation from local density approximation is thus negligible. Investigating the critical behavior at the superfluid transition, such as measuring the critical exponent α, would be an interesting development for this method, as proposed in [START_REF] Pollet | Criticality in Trapped Atomic Systems[END_REF].


N Pous

F Azaïs

L Latorre

On the use of standard digital ATE for the analysis of RF signals

Keywords: test, digital ATE, analog/RF signals, zero-crossing detection, post-processing algorithm, coherent under-sampling I

In this paper, we investigate the use of standard digital ATE for the analysis of FM-modulated RF signals. The key idea is to use the 1-bit digitizer of a digital test channel in order to convert the frequency information contained in a FM-modulated signal into a timing information contained in a digital bit stream; a post-processing algorithm based on the concept of zero-crossing detection is then employed to retrieve this information. Coherent under-sampling is exploited to extend the capabilities of test equipment with a limited sampling frequency for the analysis of high-frequency signals. The proposed approach is evaluated on two different case studies related to LTE and GSM communication standards. Both simulation and hardware experiments are presented to demonstrate the viability of the technique.

INTRODUCTION

More and more, with the development of wireless and multimedia applications, analog and RF functions become essential elements of electronic systems. In this context, it is clear that the test cost of such functions must be kept as low as possible to be competitive in the market. However in this particular case, the test is a time-consuming procedure that requires costly and dedicated test equipment.

In this work, our objective is to reduce the cost of the required test equipment. The fundamental idea is to complement standard low-cost test equipment with signal processing techniques to enable the analysis of analog/RF signals. More precisely, the idea is to use the comparator of a standard digital ATE channel as a 1-bit digitizer that converts the frequency and/or amplitude information contained in the analog/RF signal into a timing information contained in a digital signal; post-processing algorithms can then be developed to extract the analog/RF signal characteristics from the digital bit stream.

An essential motivation of this approach is that the current trend for actual circuits is to operate at lower supply voltages but with ever higher frequencies. This means that modern CMOS chips undergo a degraded voltage resolution whereas they benefit from an improved time resolution. Regarding test aspects, this trend can be translated into the fact that it is more and more difficult to perform accurate voltage measurements while it becomes easier to obtain precise timing measurements.

As a consequence, our idea is to take advantage of this trend by converting the analysis of analog/RF signals from the voltage domain to the time domain [START_REF] Rivoir | Analog to digital signal conversion method and apparatus[END_REF].

Moreover, another observation is that RF devices are often tested twice, at the wafer-level and again at the package-level. However implementing RF tests at wafer-level is extremely costly due to probing issues and the inability to perform multisite testing. It is therefore of great interest to develop test solutions applicable on low-cost test equipment that (i) provide wafer-level test coverage and (ii) permit multi-site testing. Consequently, another important objective is the ability of performing some levels of RF testing using standard digital ATE channels together with the possibility of performing multi-site testing.

Referring to literature, our approach is quite innovative. Indeed, two main strategies have been investigated up-to now to diminish the test cost of analog/RF circuits. A first strategy is to include BIST features in the circuit in order to reduce the requirements of the external test equipment. During the past ten years, a significant research effort has been carried out in this direction for analog circuits [START_REF] Khaled | Frequencybased BIST for analog circuit testing[END_REF][START_REF] Hong | On-chip analog response extraction with 1-bit Ȉ-ǻ modulators[END_REF] and more recently for RF circuits [START_REF] Negreiros | Low cost analogue testing of RF signal paths[END_REF][START_REF] Staszewski | RF Built-in Self Test of a Wireless Transmitter[END_REF]. Another strategy is to reduce the test time, and a number of works have been realized to replace the classical specification-based tests by shorter alternative tests for both analog and RF circuits [START_REF] Variyam | Enhancing test effectiveness for analog circuits using synthesized measurements[END_REF][START_REF] Haider | Low-cost alternate EVM test for wireless receiver systems[END_REF]. In this work, we adopt a slightly different approach that consists in reducing the requirements of the external test equipment but without inserting any test features within the circuit itself.

In this paper, we focus on the analysis of FM-modulated signals using the frequency reconstruction technique introduced in [START_REF] Pous | Exploiting zerocrossing for the analysis of FM modulated analog/RF signals using digital ATE[END_REF]. In order to extend the application range of the method, we introduce the use of coherent under-sampling and we evaluate the performances of the technique on two different case studies related to LTE and GSM communication standards.

The paper is organized as follows. Section II presents the basic principle of using standard digital ATE channels for frequency estimation together with the associated postprocessing algorithm for frequency demodulation. Section III introduces the use of coherent under-sampling to enable the analysis of high-speed signals using test resources with limited operating frequency. Sections IV and V discuss the performances of the technique on the two different case studies. Finally section VI concludes the paper.

978-

A. Basic idea

Our objective is to investigate the possibility of using digital test resources, cheaper than the RF resources, to analyze analog/RF signals. Based on the concept of level-crossing detection, the idea is to complement a standard digital test channel with a signal processing algorithm so that it can operate as an analog/RF receiver. Note that many works can be found in the literature regarding the use of level-crossing in many different application domains such as image and speech processing, wireless communications… However only a limited number of works deals with the use of such a technique for analog/RF test issues. Zero-crossing is used in [START_REF] Brosa | Digital signature proposal for mixed-signal circuits[END_REF] to generate a digital signature associated with Lissajou-based test. Level-crossing is used in [START_REF] Su | Dynamic analog testing via ATE digital test channels[END_REF] to evaluate offset and amplitude of a triangular analog waveform using digital resources. As illustrated in figure 1, the comparator of the digital test channel is used as a 1-bit digitizer and the resulting bit stream is stored in the ATE memory. Depending on the reference level of the comparator, the frequency and/or amplitude information contained in the analog/RF signal is converted into a timing information contained in the digital bit stream. It is then the role of the post-processing algorithm to retrieve the analog/RF signal characteristics. In this paper, we focus on frequency reconstruction, so we set the reference level of the comparator to zero. Similar solutions may be developed for amplitude reconstruction, using a comparator reference level different from zero.

B. Frequency determination

To illustrate the proposed approach, let us first consider an ideal analog sine-wave:

) 2 cos( ) ( t f A t s ⋅ π = (1)
Assume that this sine-wave is sampled by a 1-bit comparator with a zero-reference level at a sampling rate f sample . The resulting signal is a digital signal that switches from logical "0" to "1". As illustrated in figure 2, a Time Stamp (TS) can be associated with each transition, which corresponds to a zero-crossing of the analog sine-wave. Then, for each pair of rising/falling transitions, the signal frequency can be computed with:

(

) 2
where Δ TS corresponds to the time delay between a pair of successive rising and falling transitions. In the ideal case, level-crossings of the analog sine-wave result in unique transitions at the output of the comparator. Time Stamps can therefore be directly associated to the rising/falling transitions. However in practice, the analog signal is not a perfect one but a noisy signal. The resulting digital signal at the output comparator therefore exhibits multiple transitions at the vicinity of zero-crossings, as illustrated in figure 3. An algorithm is consequently required to filter these multiple transitions and associate a single Time Stamp to each zero-crossing. In this objective, a simple and robust algorithm that allows the precise estimation of Time Stamps associated with level-crossings for a noisy signal has been developed. As depicted in figure 3, the algorithm relies on a running average, i.e. counting the ratio between the numbers of "0" and "1" on a given number of samples N avg . A Time Stamp TS R is associated to a rising transition when this ratio reaches 50% with a positive slope; similarly, a Time Stamp TS F is associated to a falling transition when the running average reaches 50% with a negative slope.

C. Frequency demodulation

Let us now consider the case of an FM-modulated signal expressed by: ( ) ( ) where A c and f c are the amplitude and the frequency of the carrier signal, f m is the frequency of the modulating signal (further called "the message"), and β is the modulation depth given by:

t f t f A t S m c c FM ⋅ π β + ⋅ π = 2 sin 2 cos ) ( (3)
m d m f f A = β (4)
and f d corresponds to the maximum frequency deviation.

In order to demodulate this signal, we can use the zerocrossing Time Stamp determination defined in the previous section. More precisely, the principle consists in estimating the deviation of the signal frequency from the carrier frequency, for each period of the FM-modulated signal.

From a practical point of view, the procedure is as follows:

• For each TS R,i / TS F,i pair, i.e. for each period of the carrier signal, we compute the time t max,i . at which the FM-modulated signal reaches a maximum:

2 , , max, i F i R i TS TS t + = (5)
• For each pair of two consecutive maximums, we reconstruct the message with:

1 • © § - - ⋅ β = + c i i m m i m f t t f A t s max, 1 max, 1 ) ( ˆ at 2 1 max, max, + + = i i i t t t ( 6
)
Note that because the reconstruction is based on the processing of the sampled comparator output, the resulting reconstructed signal is a discrete signal both in time and amplitude. The timing discretization depends on the characteristics of the signal to be analyzed, i.e. the ratio between the modulating and carrier signal frequencies. The amplitude discretization not only depends on the signal characteristics but also on the sampling frequency. The quantization step can be expressed by deriving equation ( 6):

β = sample m c m f f f A Q 2 2 (7)
From equation [START_REF] Haider | Low-cost alternate EVM test for wireless receiver systems[END_REF], it is clear that for a given signal, it is possible to reduce the quantization step and so to increase the number of quantization levels, by increasing the sampling frequency.

III. ANALYSIS OF HIGH-FREQUENCY SIGNALS

The accuracy of the proposed technique obviously depends on the accuracy on the frequency estimation, which is clearly related to the ATE sampling frequency f sample . It actually depends on the ratio between the ATE sampling frequency and the frequency of the analog signal to be analyzed. So the higher the sampling frequency is, the higher the number of captured samples per carrier period and the lower the frequency estimation error. As a consequence, this may limit the range of signals that can be analyzed since the sampling frequency is limited by the ATE capabilities. As an example, standard digital ATE can typically operate up to few GHz. Targeting accuracy better than 1% for frequency estimation, 200 samples per carrier period are required [START_REF] Pous | Exploiting zerocrossing for the analysis of FM modulated analog/RF signals using digital ATE[END_REF], and signals that can be analyzed are therefore limited to few tens of MHz.

In order to extend the application range of the method, the idea is to use a coherent under-sampling strategy. Indeed, under-sampling is a commonly-used strategy to address the capture of high-frequency repetitive signals below the Nyquist rate. It basically consists in taking only few samples within a signal period and repeating this operation several times while changing the sampling phase. The resulting multiple samples can then be reordered to obtain a satisfactory sampling of the original signal. Besides, coherent sampling is a technique that guarantee that the maximum amount of information about a particular waveform exists in the sample set (i.e. there are no duplicate samples). Implementing coherent sampling therefore permits to obtain a complete, periodic waveform representation in the sample set with faster acquisition time and less computation than with non-coherent sampling. Consequently, we introduce in this section the implementation of a coherent under-sampling approach with the ATE to address the capture of high-frequency analog waveforms. where N and M are whole positive integer values. Moreover, N and M must be chosen as co-prime integers to guarantee that samples will differ from one message signal period to another.

Referring to the example of figure 4, the coherency constraints can be satisfied with the sampling period set to T sample =3μs and the number of repetitions set to M=3. After reordering, this setup permits to collect 100 uniformly timedistributed samples on one message period, whereas the direct sampling of the signal without repetition produces only 33 samples on one message period. Figure 5 shows the proposed test setup for FMdemodulation using a digital ATE and a coherent undersampling approach. The signal to be analyzed is generated with an RF source (Agilent N9310A) connected to a PS3600 digital channel of the Verigy 93K ATE. In order to implement coherent sampling, synchronization is required between the signal and sample clock frequencies. The 10MHz output of the tester is therefore used to synchronize the RF source. The postprocessing algorithm that includes samples reordering, running average computation and time stamps determination, instantaneous frequency estimation, and finally signal reconstruction is implemented in MATLAB®.

IV. "LTE-BANDWIDTH" CASE STUDY

To validate the proposed approach, we first consider a case study based on the LTE standard. The purpose is not to handle the complete communication standard, but to use the frequencies involved in such applications as a typical RF signal. Consequently referring to one mode of the LTE standard, we consider a FM-modulated signal with a carrier frequency of 2.5GHz, a modulation frequency of 40MHz and a modulation depth β=3. Moreover to evaluate the performances of the technique, we want to analyze the spectrum of the reconstructed message. In order to have an acceptable resolution after the FFT, the signal to be analyzed is made of 32 message periods, which corresponds to 2000 carrier periods; the message period T m considered in the experiments is therefore 800ns. Note that the reconstruction algorithm results in a discrete analog signal with non-uniform time distribution according to equation [START_REF] Variyam | Enhancing test effectiveness for analog circuits using synthesized measurements[END_REF]. This non-uniformity should be considered for the analysis of the signal in the frequency domain. Indeed, applying a classical Discrete Fourier Transform (DFT) algorithm on the reconstructed signal with non-uniform samples results in a spectrum that presents both leakage and harmonics. It is therefore necessary to align data on a uniform time scale [START_REF] Rivoir | Converting non-equidistant signals into equidistant signals[END_REF] or to use an Interpolated-DFT algorithm to eliminate these effects and obtain the correct spectrum of the signal. First experiments have been performed in simulation. Two different time resolutions of 10ps and 1ps have been considered, which corresponds to about 40 and 400 samples per period carrier respectively. Taking into account that the lowest value of our ATE sampling period is 2.5ns, coherent under-sampling must be used to achieve such resolutions. The under-sampling parameters that satisfy the tester constraints and the coherent requirements are given in Table I for both targeted resolutions. Also, note that a clock jitter, modeled as a Gaussian random time-shift of the sampling event with 3σ =1ps, has been included in the simulation in order to consider realistic practical conditions. resolution, with a significant reduction of the quantization in the message amplitude when improving the time resolution from 10ps to 1ps.

Hardware measurements have been performed to further support these results. However due to the limited capabilities of our RF source, hardware measurements are performed with a downscaling of 1,000 with respect to simulation, i.e. the RF source is set to output an FM-modulated signal with 2.5MHz carrier frequency, 40kHz modulation frequency and a modulation depth of 3rad. To maintain the same number of samples as in the simulation experiments, the downscaling of 1,000 is also applied on the ATE sampling frequency, which results in a time resolution of 10ns and 1ns, for 80,000 and 800,000 captured samples respectively. Note that even if the direct acquisition of the samples would be possible for the 10ns resolution, coherent under-sampling is implemented for both targeted resolutions. Results are summarized in figure 7. A good agreement is observed between these experimental results and simulation results. The second case study addresses the demodulation of a FM signal with parameters close to those found in the GSM telecommunication standard. Note that although a GMSK modulation scheme is implemented in GSM in order to control the signal spectral spreading, our study only focuses on the determination of carried symbols based on frequency demodulation.

The typical GSM signal considered here is based on a sinusoidal carrier of frequency f c =914MHz and a frequency deviation of f d =±67.7kHz corresponding to "0" and "1" symbols. Note that this deviation is extremely small with regard to the carrier frequency. In order to discriminate symbols, the precision in the measurement of the carrier instantaneous frequency must therefore be higher than 0.015%. In other words, the targeted time resolution after sampling and reordering should be as low as 0.08ps, which corresponds to about 13,500 samples collected into one single carrier period. Obviously, this is a very exigent case and it is interesting to evaluate the performance of the proposed technique under such demanding conditions.

Additionally it is worth noting that the symbol rate in GSM is about 270kHz; each symbol is therefore emitted during 3.7μs, which corresponds to 3,385 carrier periods. Considering a time resolution of 0.08ps, the amount of data to collect within one single symbol then reaches 46Mbits, which is way beyond ATE capabilities. In this study, first attempts to demodulate a GSM signal is addressed by setting the total amount of captured data to 16Mbits, which corresponds to our ATE memory capability. As a result, the capture window is lower than the duration of one symbol, and it is assumed that there is no change in the carrier frequency during the capture process. Different time resolutions are investigated, i.e. 0.1ps, 0.05ps and 0.01ps, which corresponds to 10,000, 20,000 and 100,000 samples per carrier period respectively. Obviously, coherent under-sampling is implemented to achieve such resolutions. The ATE sampling period and the number of repetitions are adjusted according to the targeted time resolution (ATE sampling frequency is adjusted around its minimum value of 2.5ns). Table II summarizes the under-sampling parameters used both in simulation and for experimental validation. Simulations have been carried out using two signal frequencies of f c0 =913.9323MHz and f c1 =914.0677MHz corresponding to "0" and "1" symbols respectively. Again, a Gaussian jitter with 3σ=1ps is considered for the sampling clock. The post-processing algorithm is applied on the captured samples to estimate, for each carrier period, the instantaneous frequency.

Results are summarized in figure 8, which gives the distribution of the instantaneous frequency calculated for each carrier period, for the three different time resolutions. It can be observed that the instantaneous frequency exhibits a Gaussian distribution well-centered on the expected frequencies for "0" and "1" symbols, for the three different cases. On the other hand, the standard deviation strongly depends on the time resolution: the finer the time resolution, the lower the standard deviation. In particular, the standard deviation reduces down to about 30kHz when using a 0.01ps resolution, which permits to discriminate symbols based on a single period measurement. More generally, these results demonstrate the potentialities of the proposed approach to handle GSM-type signals. However it should be highlighted that the use of a very fine time resolution may not be the more efficient option in terms of acquisition time and subsequent post-processing. An interesting alternative may be to use a degraded time resolution together with an averaging approach. For instance with a 0.1ps resolution, it is possible to achieve a standard deviation comparable to the one obtained with 0.01ps resolution by averaging the estimated frequency over 4 carrier periods only. Finally, hardware measurements have been performed to further validate these results. A simple 1GHz signal send to two distinct tester channels on the same PS3600 board and the post-processing algorithm has been applied to estimate the frequency at each signal period. Results are summarized in figure 9. Frequency distribution is again very well-centered on the expected 1GHz frequency and the standard deviation reduces when the time resolution improves. However compared to simulation results, the standard deviation is about three times higher on channel 1 and more than 5 times higher on channel 2. This difference may come from the jitter on the capture event that can have slightly different values from one channel to another.

VI. CONCLUSION

In this paper, we have validated a technique that permits to reconstruct the message contained in a FM-modulated signal using standard digital ATE channels, converting the analog/RF signal into a bit stream whose transition time information represents the analog/RF signal characteristics. We have also established that coherent under-sampling can efficiently complement the technique to process signals at frequencies much higher than the ATE acquisition rate. Validations have been performed through simulations and experimental measurements, and the viability of the technique has been illustrated on two case studies with very different characteristics.

Figure 1 .

 1 Figure 1. Primary idea

Figure 2 .

 2 Figure 2. Frequency determination from zero-crossing

Figure 3 .

 3 Figure 3. Frequency determination from running average

Figure 4 .

 4 Figure 4. Coherent under-sampling Figure 4 illustrates the principle of coherent under-sampling on a simple case of FM-modulated signal with f c =100kHz and f m =10kHz. To ensure coherency, the sampling frequency f sample , the message frequency f m , the number of samples N and the number of repetitions of the message signal M must fulfill the following relationship:

Figure 5 .

 5 Figure 5. FM-demodulation setup using digital ATE

Figure 6 .

 6 Figure 6. Simulation results for the "LTE" case studyResults are summarized in figure6that gives both the timedomain and the spectral-domain representation of the reconstructed signals for the two resolutions of 10ps and 1ps. It clearly appears that in both cases, the post-processing algorithm permits to retrieve the message. As expected, the quality of the reconstructed signal increases with time

Figure 7 .

 7 Figure 7. Hardware validations for the "LTE" case study

Figure 8 .

 8 Figure 8. Simulation results for the "GSM" case study

Figure 9 .

 9 Figure 9. Hardware validations for the "GSM" case study

1-4244-5833-2/10/$26.00 ©2010 IEEE

	II. PRINCIPLE

TABLE I .

 I UNDER-SAMPLING PARAMETERS FOR "LTE" CASE STUDY

	Final resolution	Number of	ATE sampling	Number of samples
	after reordering	repetitions M	period Tsample	N
	10 ps	601	6.010 ns	80,000
	1 ps	6007	6.007 ns	800,000

TABLE II .

 II UNDER-SAMPLING PARAMETERS FOR "GSM" CASE STUDY

	Final resolution after reordering	Number of repetitions	ATE sampling period	Number of samples / carrier period	Number of carrier periods
	0.1 ps	~ 22,800	~ 2.5 ns	10,000	1,600
	0.05 ps	~ 45,700	~ 2.5 ns	20,000	800
	0.01 ps	~ 228,500	~ 2.5 ns	100,000	160

ACKNOWLEDGMENTS

The authors would like to acknowledge B. Pradarelli, R. Lorival, P. Coscat and L. Ollivier for technical support during hardware validation. Hardware experiments have been performed using the equipment of the CNFM Test Resource Center in Montpellier sponsored by Verigy® and Region Languedoc-Roussillon.