Time-Frequency Energy Distributions Meet

Compressed Sensing Patrick Flandrin, Fellow, IEEE and Pierre Borgnat, Member, IEEE Abstract-In the case of multicomponent signals with amplitude and frequency modulations, the idealized representation which consists of weighted trajectories on the time-frequency (TF) plane, is intrinsically sparse. Recent advances in optimal recovery from sparsity constraints thus suggest to revisit the issue of TF localization by exploiting sparsity, as adapted to the specific context of (quadratic) TF distributions. Based on classical results in TF analysis, it is argued that the relevant information is mostly concentrated in a restricted subset of Fourier coefficients of the Wigner-Ville distribution neighbouring the origin of the ambiguity plane. Using this incomplete information as the primary constraint, the desired distribution follows as the minimum ℓ1-norm solution in the transformed TF domain. Possibilities and limitations of the approach are demonstrated via controlled numerical experiments, its performance is assessed in various configurations and the results are compared with standard techniques. It is shown that improved representations can be obtained, though at a computational cost which is significantly increased.

Index Terms-time-frequency, localization, sparsity.

EDICS Category: SSP-NSSP

I. INTRODUCTION

C HIRPS (i.e., transient amplitude and frequency modu- lated (AM-FM) signals) are ubiquitous in nature and man-made systems. They are commonly encountered in a variety of fields ranging from audio (speech, music) to wave physics (whistling atmospherics, gravitational waves), from engineering (radar, sonar) to biology and medicine (EEG epilectic seizures, EMG uterine contractions), or even in mathematics (Weierstrass and Riemann functions), see, e.g., [START_REF] Flandrin | Time-frequency and chirps[END_REF] and references therein. From the point of view of exploratory data analysis, a fine and non parametric characterization of chirps is often needed, e.g., as a hint for some further modeling. Since chirps are waveforms whose structure explicitly involves timevarying properties with respect to both amplitude and spectral content, they naturally call for time-frequency (TF) analysis tools. Because of the Fourier uncertainty between time and frequency, any TF approach is however faced with difficulties which become all the more challenging as signals are of a shorter duration. Albeit standard in its motivation, this problem has not yet received its ultimate solution, and the purpose of this paper is to revisit it in a new way that combines basics of P. Flandrin and P. Borgnat are with Université de Lyon, Physics Department (UMR 5672 CNRS) of École Normale Supérieure de Lyon, 46 allée d'Italie, 69364 Lyon Cedex 07 France. Tel.: +33 (0)4 72 72 81 60; Fax: +33 (0)4 72 72 80 80; E-mail: {flandrin,pborgnat}@ens-lyon.fr. Part of this paper was first presented at the IEEE Int. Conf. on Acoust. Speech and Signal Proc. ICASSP-08 (Las Vegas, NV) [START_REF] Borgnat | Time-frequency localization from sparsity constraints[END_REF] and at the 8th IMA Int. Conf. on Math. in Signal Proc. (Cirencester, UK) [START_REF] Flandrin | Sparse time-frequency distributions of chirps from a compressed sensing perspective[END_REF].

"classical" TF analysis with recent advances from Compressed Sensing (CS) (for a general overview of CS, the interested reader is referred to [START_REF]Special Section on Compressive Sampling[END_REF], [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF], as well as to the resources webpage [6]).

It is worth stressing the fact that the approach reported here differs in several respects from related attempts. First, as compared to Gabor-based approaches (see, e.g., [START_REF] Pfander | Sparsity in time-frequency representations[END_REF]) which assume sparsity in linear representations, we make here the choice of working with quadratic energy distributions as TF representations, the rationale being to favor structured components in the form of localized chirps along TF trajectories rather than independent atoms. Second, whereas most "classical" CS-based methods involve random measurements, a deterministic approach will be followed here for the selection of a small subset of observed values used as constraints. For sake of incoherence between the targeted representation and its measurements, such a selection will be done in the Fourier domain in the spirit of the emblematic examples of phantom recontruction from Fourier slices discussed in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] or of "Compressed Sensing MRI" [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [START_REF] Lustig | Compressed Sensing MRI[END_REF], thus contrasting with other possibilities of deterministic sensing [START_REF] Applebaum | Chirp sensing codes: Deterministic compressed sensing measurements for fast recovery[END_REF]. A major difference exists however with such Fourier-based CS approaches since, as it will be made explicit in Sect. IV-B2, the addressed problem will be shown to differ from a classical reconstruction by moving to the construction of a TF distribution which would not be attainable otherwise, regardless of the number of measurements in the Fourier domain. Finally, whereas specific improvements to Fourier-like TF methods could be obtained by means of model-based approaches, we will adopt here a model-free perspective, the only assumption being that the representation which of interest for is intrinsically sparse in the TF domain since the energy of a chirp is expected to be essentially distributed along a 1D curve of the 2D plane.

More precisely, the paper is organized as follows. In Section II, the issue of chirp localization is addressed and classical TF solutions are briefly recalled. Considering that localized TF distributions of chirps are sparse in the TF plane, Section III gives the rationale for revisiting the problem from a CS perspective. This is detailed further in Section IV where the effectiveness and performance of the approach are discussed, in particular with respect to the selection of Fourier samples for which a rule-of-thumb criterion is proposed. Whereas this paper is not focused on the algorithmic part of the method, some computational issues are addressed, and limitations as well as possible extensions are finally discussed in Section V.

II. TIME-FREQUENCY LOCALIZATION

A. AM-FM signals as time-frequency trajectories

If we consider a signal made of the superimposition of a finite number of AM-FM components:

x(t) = K k=1 a k (t) e iϕ k (t) ,
it is natural to attach to it an idealized TF distribution (TFD) ρ x (t, f ) which essentially distributes the total energy along TF trajectories according to:

ρ x (t, f ) = K k=1 a 2 k (t) δ (f -φk (t)/2π) . (1)
In such a picture, each component is characterized at each time instant by essentially one instantaneous frequency (which, in a first approximation, can be identified to the phase derivative), weighted by the corresponding instantaneous power.

Except for very special cases, there is no general methodology to automatically get a distribution as in [START_REF] Borgnat | Time-frequency localization from sparsity constraints[END_REF]. In the case of a single component (K = 1), it is well-known [START_REF] Flandrin | Time-Frequency/Time-Scale Analysis[END_REF] that a perfect localization can be attained for pure FM signals with a linear modulation (a 1 (t) = 1 and φ1 (t) = 2π(f 0 + α t)) by using the Wigner-Ville Distribution (WVD):

W x (t, f ) = +∞ -∞ x t + τ 2 x * t - τ 2 e -i2πf τ dτ. (2)
Although this property can be extended to some forms of nonlinear FMs (e.g., Bertrands' distributions for powerlaws [START_REF] Flandrin | Geometry of affine time-frequency distributions[END_REF]), it is generally at the expense of a substantially increased complexity in the definition (and the computation) of the distributions, with furthermore the limitation of being adapted to some specific type of FM only and to not extend to multicomponent situations. For this last point, the wellknown drawback of energy distributions is to obey a quadratic superposition principle which creates cross-terms in between any two components of a signal, and thus significantly reduces the readability of Wigner-type distributions [START_REF] Hlawatsch | The interference structure of Wigner distribution and related time-frequency representations[END_REF].

B. Classical techniques of TF localization

The aforementioned difficulties have led to many developments during the last 20 years but, unfortunately, since both localization and creation of cross-terms result from the very same mechanism [START_REF] Flandrin | Cross-terms and localization in time-frequency energy distributions[END_REF], it turns out that trying to impose simultaneously localization and cross-terms reduction is faced with a trade-off that can be viewed as a form of time-frequency uncertainty principle. The simplest way to understand where this trade-off comes from (and how to manage it) is to interpret the WVD in its 2D Fourier transform domain.

By definition, the WVD admits a 2D Fourier transform which is referred to as the ambiguity function (AF) and reads

A x (ξ, τ ) = +∞ -∞ x t + τ 2 x * t - τ 2 e i2πξt dt.
If we introduce the TF shift operator T ξ,τ which acts on signals x(t) ∈ L 2 (R) as

(T ξ,τ x) (t) := x(t -τ ) e -i2πξ(t-τ /2) ,
we readily get that A x (ξ, τ ) = x, T ξ,τ x , i.e., that the AF is nothing but the inner product between the analyzed signal and its TF shifted versions. By construction, the AF can thus be viewed as a TF correlation function and, as such, it exhibits most properties of a correlation function, including hermitian symmetry and the fact that its modulus is maximum at the origin. Moreover, in the case of multicomponent signals, the total AF consists of both auto-components neighbouring the origin of the plane and cross-components mostly located at a TF distance from the origin which directly depends on the TF separation between the individual components and that are the Fourier images of the undesired cross-terms in the TF plane. This observation early prompted [START_REF] Flandrin | Some features of time-frequency representations of multicomponent signals[END_REF] to propose improvements upon the WVD by weighting the AF around the origin of the plane prior to applying an inverse 2D Fourier transform: the more restricted the weighting domain, the more effective the cross-terms suppression but, at the same time, the more decreased the TF localization. This procedure (which defines Cohen's class [START_REF] Flandrin | Time-Frequency/Time-Scale Analysis[END_REF] on geometrical grounds and can be improved by adapting the weighting domain to the signal structure [START_REF] Baraniuk | Signal-dependent time-frequency analysis using a radially Gaussian kernel[END_REF]) proved reasonably effective, but it basically faces a trade-off between cross-terms reduction and localization. This motivated the search for improved techniques and, amongst the other approaches that have been proposed, the one that is referred to as reassignment plays a prominent role. The rationale behind reassignment (see [START_REF] Flandrin | Time-Frequency Reassignment -From Principles to Algorithms[END_REF] for a comprehensive presentation of reassignment techniques) is that it basically exploits the phase information that is usually discarded in simple quadratic distributions such as spectrograms (or scalograms). In this respect, reassigned spectrograms (RSP) proved in particular extremely efficient to approach (1) and will thus be used in Section IV for the sake of comparison.

III. LOCALIZATION FROM SPARSITY CONSTRAINTS

If the analyzed signal is given in discrete-time and supposed to be of dimension N in time, its TFD is of dimension N 2 when computed over N frequency bins. However, assuming that K ≪ N , i.e., that the number of components is much smaller than the dimension of the signal, the targeted TFD which is supposed to satisfy (1) is distributed over the plane in a very sparse way, with only K 1D trajectories where at most KN values are expected to be non-zero. Imposing such a sparsity is therefore a new way of approaching the problem.

A. Principle

Based on the arguments recalled in Section II-B, the principle of the proposed sparsity-based approach to TF localization is thus very simple: it simply consists of selecting a suitable collection of AF samples neighbouring the origin of the plane in a given domain Ω(ξ, τ ) and searching for the sparsest TFD ρ such that its 2D Fourier transform F{ρ} coincides with the original AF over Ω. As has been mentioned previously, masking AF around the origin has the desired effect of reducing cross-terms, but at the expense of degrading localization if a crude Fourier inversion is performed without any constraint. Selecting a few samples only of the AF around 1) is given in the top row, in between the Wigner-Ville Distribution (WVD) and the reassigned spectrogram (RSP) that has been "optimized" as explained in Sect. IV-A. The corresponding total ambiguity function (AF) of size 128 × 128 (defined as the 2D Fourier transform of the WVD) is given on the left of the top row, with its restriction to a domain Ω consisting of 13 × 13 samples neighbouring the origin of the plane on the left of the bottom row. The rest of the bottom row presents the results obtained from optimizations based on the only knowledge of this restriction of the AF, with the exact minimum ℓ 1 -norm solution according to [START_REF] Flandrin | Time-frequency and chirps[END_REF] in between the minimum ℓ 2 -norm solution and the approximate minimum ℓ 1 -norm solution according to (4) with ǫ = 0.05 x 2 . For all diagrams, amplitudes are grey coded logarithmically, with a dynamic range of 18 dB.

the origin is therefore the first ingredient for ensuring crossterms reduction, and it is the further assumption of few chirp components that is instrumental for guaranteeing TF localization. In this respect, the recourse to CS algorithms can be viewed as a form of constrained AF extrapolation.

B. Constraints

Looking for a perfectly localized solution such as (1) would require minimizing the total number of non-zero coefficients, i.e., the ℓ 0 -"norm" of the TFD. While this turns out not to be practicable from a computational viewpoint, a series of recent works (see, e.g., [START_REF] Donoho | Compressed sensing[END_REF], [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF], and [6] for a comprehensive list of references) have shown that a near-optimal solution can be attained at a sensibly more affordable cost by minimizing the ℓ 1 -norm, reducing the problem to the solving of a linear program. It is this technique which is proposed to be followed here, the desired localized TFD ρ x (t, f ) being therefore the solution of the constrained minimization problem:

ρ x = arg min ρ ρ 1 ; F{ρ} -A x = 0| (ξ,τ )∈Ω . (3)
The primary constraint which is given by (3) imposes a strict equality over Ω in the AF domain. This however can be relaxed [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF] according to

ρ x = arg min ρ ρ 1 ; F{ρ} -A x 2 ≤ ǫ| (ξ,τ )∈Ω , (4)
where ǫ is a user-specified bound. Both possibilities (3) and (4) will be considered in the following, with the corresponding ℓ 1 solutions referred to as, respectively, ℓ 1 -eq (for ℓ 1 with equality) and ℓ 1 -app (for ℓ 1 with approximation).

One can remark that, the present approach being based on optimization, it is possible to easily impose further constraints besides (3) or ( 4). This has not been followed in this work, but it should certainly be interesting to explore further such a possibility (with respect, e.g., to regularity, marginals, . . . ).

IV. CS-BASED TF ANALYSIS IN ACTION

The feasability of the method has been tested on simple, yet informative examples. Since the purpose of this paper is not to contribute in a new way to the solving of optimization problems such as those stated in (3) and (4), the reported examples are based on a straightforward application of existing solutions. For the sake of completeness, basics of the corresponding methods are however briefly recalled in the Appendix. All the computations presented in this section have been made in MATLAB, with the TIME-FREQUENCY TOOLBOX [21] for the TF computations and the ℓ 1 -MAGIC TOOLBOX [START_REF] Candès | ℓ 1 -MAGIC: Recovery of Sparse Signals via Convex Programming[END_REF] for the optimization, but alternative algorithmic procedures have been considered as well, with similar results: this will be detailed further in Section IV-D.

A. A first comparison

Figure 1 illustrates the principle of the proposed approach and compares different TFDs in the case of a 128 points waveform made of the superimposition of a linear and of a sinusoidal FM, both modulated in amplitude with a Gaussian. In this example, the different optimizations have been based on the only knowledge of the 13 × 13 Fourier samples of the WVD neighbouring the origin of the AF plane, i.e., on a subset of about 1% only of the total number of AF coefficients. With so few AF coefficients, it is clear that a crude Fourier inversion (which corresponds to the ℓ 2 -norm solution of the optimization problem) ends up with a very poor TFD. From a qualitative point of view, it then turns out from this figure that the minimum ℓ 1 -norm solutions (either exact, according to (3) or approximate, according to (4)) are much more effective in terms of localization, even as compared to the reassigned spectogram which is known to usually behave best for this kind of signal. It has to be noted that, since a reassigned spectrogram depends on the choice of a window length, experiments have been conducted with various lengths and the reported quantitative results correspond to the best performance that has been obtained. In practice, such an optimization cannot be achieved, but making use of this "optimum" reassigned spectrogram has been chosen so as to get a form of bounds on the achievable performance (such bounds are likely not to be attained in real situations). Comparing the exact and approximate ℓ 1 solutions, it also appears that the former is in some sense "too sparse" as compared to the latter, with energy contributions well localized on the TF trajectories, but in a discontinuous, spiky fashion.

B. Selection of Fourier samples

The result reported in Figure 1 can be thought of as very peculiar in the sense that it heavily relies on the specific choice of a domain, regarding area (why 13 × 13 Fourier samples?) as well as shape (why a fixed square domain?). All the appreciations commented above can then be quantified further in terms of both the achieved performance with respect to the actual model, and the influence of the AF domain from which the optimization is conducted.

1) Heisenberg cardinality: Concerning the cardinality of the domain (first assumed to be square), the result is reported in Figure 2 which displays a localization measure (the Rényi entropy of order 3, which has proved to be a good measure of localization for energy TFDs [START_REF] Baraniuk | Measuring time-frequency information content using the Rényi entropies[END_REF]) together with the ℓ 1distance to the model, both as functions of the relative sparsity measure card(Ω)/N 2 . What is revealed by this figure is that both minimum ℓ 1 -norm solutions are generally better localized (i.e., have a smaller Rényi entropy) than the other considered TFDs, with even an entropy that might be smaller than the model one whose value is in this case 6.37. In particular, the exact ℓ 1 solution has always the minimum entropy but, as evidenced by Figure 1, this is due to an oversparse, discontinuous structure which results in a larger distance to the actual model. The evolution of this distance shows that the best behavior is obtained with the approximate minimum ℓ 1 -norm solution, the minimum being obtained when using only about 1% of the N 2 AF samples. The small value of this sparsity measure is indeed a consequence of the fact that, as it will be justified further in Sect. IV-B2, the required number of measurements in the Fourier domain scales as card(Ω) ≈ N . Additionnally, it is worth noticing that, while the WVD we started with is known to attain negative values and whereas no positivity constraint has been imposed, the minimum ℓ 1norm solutions happen to be almost positive.

A qualitative interpretation of why a cardinality of the order of the signal length leads to the best result can come from an uncertainty argument. Indeed, AFs (as TFDs) are known to obey uncertainty principles which prevent them to be pointwise localized in their respective planes [START_REF] Price | Bounds on the volume and height distribution of the ambiguity function[END_REF]. Auto-terms of a multicomponent signal are thus expected to exist over a domain (neighbouring the origin of the plane, see Section II-B) whose area is at least of the order of a "Heisenberg cell", i.e., roughly N if the AF of a N points signal is computed over a N × N grid. Most of the information in such a cell being coded by the modulus and the phase of the AF, there is no point in enlarging further the domain, since the main effect of such an enlargement would be to possibly include undesired cross-terms influences in the constraints.

2) Oracle: Some further appreciation can be gained from a comparison with an "oracle", based on the assumed knowledge of the ideal distribution. In this case, the linear program ( 4) is run exactly as in Figure 1, while replacing mutatis mutandis the WVD by the TF model [START_REF] Borgnat | Time-frequency localization from sparsity constraints[END_REF]. This is illustrated in Figure 3, where an optimum cardinality (namely |Ω| * ≈ N 2 /100 ≈ N ) is shown to exist, at least qualitatively, for trading off autoterms localization and cross-terms reduction.

The compared evolutions of the ℓ 1 solutions derived from the WVD and the model are clearly similar in a first regime, whereas they diverge after the domain cardinality has attained some critical value, the higher localization expected for autoterms being hampered by the off-spring of cross-terms that appear as spikes in the representation due to their oscillatory nature. If no such spiky terms would be present, a useful criterion for choosing the cardinality of Ω could be the entropy: the smaller the entropy, the more localized the TFD. In order to prevent possible values of entropy due the spiky structure stemming from undesired additional cross-terms, the sole entropy has to be penalized in some way. This has been investigated in [START_REF] Flandrin | Sparse time-frequency distributions of chirps from a compressed sensing perspective[END_REF] (but not reported here) by using total variation as such a penalty function, ending up with the identification of an optimum cardinality |Ω| * that is consistent with the results displayed in Figures 2 and3.

At this point, it might be worth stressing the fact that the approach described here differs from a reconstruction problem (as those considered, e.g., in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [START_REF] Lustig | Compressed Sensing MRI[END_REF]) in the sense that there would be no point in exactly recovering the whole WVD from a limited set of measurements (in fact, such a perfect reconstruction is of course obtained when the entire ambiguity plane is chosen as the measurement domain, see the far right solutions of Figure 3). The situation is much more that of a construction problem in which is created some idealized object which does not exist per se prior optimization. From the oracle perspective, a suitable selection of AF samples should thus explicitly correspond to a subset of the Fourier transform of the ideal localized TFD, and no more of the WVD. Based on fundamental CS results for partial Fourier coefficients (see [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] and the review in [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]), the cardinality of the domain Ω should be therefore card(Ω) = O KN log(N 2 ) for the recovery of K AM-FM trajectories of N points each in a TF domain of size N 2 . In practice, the typical value obtained from this formula seems to be higher that the previous findings reported in Figure 2 and those based on the oracle. This is first related to the construction nature of our problem, when the CS results are given for a reconstruction problem. Second, a lower limit for a ℓ 1 minimization problem to find the sparsest solution is given by a phase transition behaviour, as discussed in [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. In that case, the cardinality of the domain Ω could be as low as card(Ω) = O (KN log(N/K)) in this problem. With the numerical values used in Figure 1, this leads to a magnitude of order 10KN (up to a constant of O( 1)), that is of the same order of magnitude as the Heisenberg cardinality discussed before.

3) Domain geometry: Up to now, emphasis has been put on the area of a domain Ω assumed to be square. This is of course an unnecessary restriction that has to be questioned. To this end, Figure 4 displays ℓ 1 distances to the model, as in Figure 2, with a domain Ω whose area is kept fixed and equal to its "optimum" value 169, but with a rectangular or elliptic shape with a varying aspect ratio. The latter being computed as "r = Doppler semi-axis/delay semi-axis", it turns out that the best performance is obtained when the selected ambiguity domain is slightly elongated along the delay axis, with almost uniformly improved results within the range 1/3 r 1 (i.e., equivalently, -0.4 log 10 (r) 0). While such numerical values are likely not to be universal and might depend on the analyzed signal, the interpretation is believed to be general in the sense that aspect ratios r < 1 tend to favor TF coherence in the time direction, a feature that is characteristic of most chirps which are naturally considered as "frequencies that vary with time". Would the signal TF structure be organized differently as, e.g., in dispersive systems where frequencydependent epochs are the most relevant signatures, aspect ratios greater than one would be better. When some anisotropy is expected in a TFD, the recipe is that this should be taken into account as a prior when choosing Ω.

A further remark can be done in the same direction, which points to a basic difference between the approaches based on linear decompositions and the quadratic one proposed here. If sparsity was assumed to exist with respect to some dictionary of TF atoms, their identification (via matching or basis pursuit) would end up with a collection of elementary components (typically, Gabor logons) that could be located along the TF trajectories of interest but as independent components rather than as a globally coherent structure. From a quadratic TF perspective [START_REF] Flandrin | Cross-terms and localization in time-frequency energy distributions[END_REF], and in the spirit of the "cross-terms deleted representations" [START_REF] Qian | Wigner distribution decomposition and crossterms deleted representation[END_REF], this corresponds to a situation where the inner interference terms that guarantee the coherently organized TF structure of chirps are just ignored. In this respect, a flexible improvement can come from the taking into account of such terms that, by nature, have short (TF) distance interactions [START_REF] Qian | Decomposition of the Wigner distribution and time-frequency distribution series[END_REF]. This is also what is done in the present approach since, by nature, inner interference terms mix up with auto-terms in the vicinity of the origin of the ambiguity plane, without any clear-cut separation. They are thus automatically taken into account when restricting Ω to an Heisenberg cell, with a representation all the more satisfactory in terms of chirp continuity as the geometry of the domain is adapted to the TF correlation structure of the analyzed signal.

This last remark prompts to look for even more adaptivity in the domain geometry, beyond convex domains whose main axes would necessarily coincide with those of the ambiguity plane. At first sight, it might be expected that the use of adapted kernels (as proposed, e.g., in [START_REF] Baraniuk | Signal-dependent time-frequency analysis using a radially Gaussian kernel[END_REF]) would prove useful. This however seems not to be the case, which does not necessarily come as a surprise. Indeed, besides an assumed sparsity of the solution, the other ingredient for a successful CS-based approach is that of its incoherence with the measurements on which it is based. Operating in the Fourier domain clearly Fig. 3. Qualitative selection of the optimal cardinality of Ω -Top row: TF solutions obtained from the WVD as in Figure 1, for increasing values of the relative cardinality card(Ω)/N 2 . The evolution shows that too small a cardinality does not permit to gain enough information about auto-terms whereas a too large one leads to take into account more and more undesired cross-terms, converging eventually to the actual WVD when the support of Ω identifies to the whole AF plane. Bottom row: companion "oracle" solutions obtained the same way, but based on the ideal distribution in place of the WVD. Comparison of both evolutions evidences that an optimum cardinality |Ω| * exists for trading off auto-terms localization and cross-terms reduction (about 1% of the AF support in the present case). For all diagrams, amplitudes are grey coded logarithmically, with a dynamic range of 18 dB. goes this way, but adapting the selected subset to the signal structure goes the opposite.

C. More examples

Another example is shown in Figure 5 about the application of the method to a synthetic signal with crossing chirps, a case that is known to be difficult to recover from quadratic TF distributions (due to the interferences, [START_REF] Hlawatsch | The interference structure of Wigner distribution and related time-frequency representations[END_REF]). The optimum magnitude for the Heisenberg cardinality is, for this example with two linear chirps with Gaussian modulated amplitudes, card(Ω) ≃ 0.005 N 2 , in agreement with the previous discussions. The distribution resulting of ℓ 1 -app is satisfactorily close to the model one, see the quantitative results given in the following Simulation examples of this kind could of course be multiplied, and we will give only one more (already considered as a test-signal for reassignment methods in [START_REF] Flandrin | Time-Frequency Reassignment -From Principles to Algorithms[END_REF]). It consists of a signal that is similar to the one used in the first part of this paper, while having a structure that is more complex (and difficult to deal with) in three respects:

1) it has three components, with one of them embedded in the cross-terms resulting from the interference between the two other ones; 2) the newly introduced component is a Gabor logon, which a priori does not enter directly the AM-FM model considered so far (no TF "trajectory" in the plane); 3) some noise is added (with a SNR fixed to 10 dB). As for the previous example, results are plotted in Figure 6, with a performance summarized in the following 7 is concerned with some real data, namely the classical benchmark signal of a bat echolocation call 1 of effective length 400, zero-padded to N = 512 data points. In this case too, the approximate minimum ℓ 1 -norm solution (4) with card(Ω) = 23 × 23 ≈ N compares very favorably with a reassigned spectrogram in terms of localization, with even some smoother regularity along TF trajectories.

D. Computational cost

Given the above reported findings, the new discussed approach is no doubt attractive in terms of its ability to give sharply localized TFDs in the case of transient, short duration AM-FM multicomponent signals. There is however a price to pay for this performance, which is a quite heavy computational 1 The authors wish to thank Curtis Condon, Ken White, and Al Feng of the Beckman Institute of the University of Illinois for the bat data and for permission to use it in this paper. cost. To illustrate this point, one can report the average computation times associated to the simulations of Figure 2, all computations having been performed with MATLAB R2008a in similar conditions (Intel Core 2 Duo 3.06 GHz on Mac OS 10.5). The result is that, whereas such an average time is 0.007 s for a WVD and 0.23 s for a reassigned spectrogram, it goes up to 50 s for the exact ℓ 1 solution and even to 152 s for the approximate one, using the implementation of the ℓ 1 -MAGIC TOOLBOX [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. It thus turns out that computing a CS-based TFD may differ from a classical approach by several orders of magnitude.

There is an increasing literature about alternative algorithms that could be used to solve programs having the form of (4) and this could lead to improvement in that respect. For instance, using the GPSR algorithm proposed in [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF] to solve the quadratic program associated to (4), the computational cost goes down to 43 s (using the GPSR-BB with continuation) and, for a suitable choice of the regularization parameter (see Appendix), the results are identical in terms of ℓ 1 distance to the model and Rényi entropy of the resulting sparse TF distribution. Other approaches (see a comprehensive list on the CS resource page [6]) based on greedy strategies (for instance refinements on Matching Pursuit, e.g., [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF], [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]), or on different strategies for convex relaxation such as Iterative Thresholding (e.g., [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems[END_REF], [START_REF] Daubechies | Accelerated Projected Gradient Method for Linear Inverse Problems with Sparsity Constraints[END_REF], [START_REF] Blumensath | Iterative Thresholding for Sparse Approximations[END_REF] and contributions improving those), could lead to more efficient recovery procedures in the current context, though we do not pursue here this question further on. Our discussion about the obtained non-negativity of the minimum ℓ 1 -norm solutions can be complemented by noting that using some algorithms would make it possible to impose a condition on non-negativity of the solution, as proposed, for instance, in the GPSR method [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF]. The issue of the choice of the better algorithm for the study proposed here is only sketched and will be part of further work.

V. CONCLUSION AND PERSPECTIVES

A new approach has been proposed for getting sharply localized TFDs in the case of multicomponent AM-FM signals by making use of the assumed sparsity of their energy distribution in the TF plane. Once the principle of the method has been given, there is clearly plenty of room for more thorough investigations and further developments.

As has been briefly mentioned in Section III-B, one particular interest of the present approach which is based on optimization, is that further constraints can be imposed besides (3) or (4). One can think, e.g., of the marginalization properties attached to unit cross-sections in the AF plane [START_REF] Flandrin | Time-Frequency/Time-Scale Analysis[END_REF]. This however is not necessarily relevant in the context of a sharp TF localization since, in the case of multicomponent signals, highly oscillatory behaviours along TF trajectories will be favored. A more interesting variation is to favor regular time evolutions, what can be achieved explicitly by imposing specific AF cross-sections or, implicitly, by using as a starting point a time smoothed WVD.

Such extensions would also suggest to conduct further experiments in comparison with distributions within Cohen's class (with fixed or adapted kernels): a comprehensive treatment of this point is clearly beyond the scope of this paper, but this is a perspective that is certainly worth investigating in future developments. In a different direction, one can mention that CS-based approaches have been very recently proposed for the purpose of nonstationary spectral estimation in the case of stochastic processes [START_REF] Jung | Compressive spectral estimation for nonstationary random processes[END_REF], [START_REF] Jung | Compressive nonstationary spectral estimation using parsimonious random sampling of the ambiguity function[END_REF], in place of the timefrequency localization problem that has been considered here for deterministic chirp signals. While their objectives are somehow different, it turns out that the CS-based rationale of the two approaches is quite similar, and it would be interesting to merge them in a common perspective.

Whereas the heavy computational load can be considered as a severe drawback, it has been mentioned that improving upon this limitation is a current topic of very active research, with advances in this direction expected to be most benefitial for the present method. Moreover, it is worth stressing that the proposed optimization-based approach is not primarily intended to long duration signals. It is rather relevant for transient, highly nonstationary signals and, as such, it is expected to provide some nice addition to the existing toolbox of non-parametric TF methods aimed at exploratory data analysis.

APPENDIX

There has been a considerable interest in the solving of optimization problems such as the ones in eqs. ( 3) and ( 4) in the past years, thanks to the emergence of the "Compressed Sensing" (also called "Compressive Sampling", and referred to as CS) [START_REF]Special Section on Compressive Sampling[END_REF], [START_REF] Donoho | Compressed sensing[END_REF], [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. The purpose of this Appendix is to provide the reader with an overview of this domain and of some methods used to solve these optimization problems, and to refer to the articles proposing them.

A. Review on Compressed Sensing and optimization problems with ℓ 1 norm

The minimization problem ℓ 1 -eq of eq. ( 3) is known as "Basis Pursuit" since it corresponds to [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]:

min ρ ρ 1 ; subject to O{ρ} = b, (5)
where ρ ∈ R n , b ∈ R k and O is an operator associated to a k × n matrix representation. The ℓ 1 -app problem of eq. ( 4) is a relaxed version of ℓ 1 -eq, where the equality constraint is replaced by a quadratic constraint:

min ρ ρ 1 ; subject to O{ρ} -b 2 ≤ ǫ. (6)
Both problems are related in the sense that eq. ( 5) has the form (6) with ǫ = 0.

Compressed Sensing (CS) deals with these problems when the number k of linear measurements available is much smaller than the number n of unknowns [START_REF]Special Section on Compressive Sampling[END_REF], and many methods have been proposed to solve both of them. Note that, in our context, O being a Fourier transform on the limited set Ω in the ambiguity domain, the number of measurements k = card (Ω) is indeed much smaller than n = N 2 , the number of points of the TFD; the situation is thus equivalent to a CS one.

B. A short lecture guide for optimization with ℓ 1 -norm

A review of all the algorithmic procedures that have been proposed so far is outside the scope of this Appendix, see for instance discussions in [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF], [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF], [START_REF] Candès | ℓ 1 -MAGIC: Recovery of Sparse Signals via Convex Programming[END_REF] or the PhD thesis [START_REF] Tsaig | Sparse solution of underdetermined linear systems: Algorithms and applications[END_REF] for a comprehensive discussion on modern algorithms. Here, we will discuss how the two optimization problems stated above can be solved for our work, using classical techniques from the CS field.

When b is real, the problem ℓ 1 -eq can be recast and solved as a Linear Program, as shown in the seminal work on Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], and solved using interior point methods. In our context, the operator O is the Fourier transform and the constraint b corresponds to a restriction to some observations in the Fourier domain. This situation was considered in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], proving that sparse vectors can then be recovered exactly through ℓ 1 minimization (see also [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [START_REF] Lustig | Compressed Sensing MRI[END_REF] for use of CS principles in the Fourier domain). When separating on different coordinates in b the mean and the real and imaginary parts of the Fourier transform, b turns out to be real and the primal-dual interior point method proposed in l1eq_pd.m of [START_REF] Candès | ℓ 1 -MAGIC: Recovery of Sparse Signals via Convex Programming[END_REF] can be used to solve [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF], once recast as a constrained linear program.

As argued in the article, the solutions of ( 5) being often too sparse to be considered as a good TFD, solving the problem (6) is more relevant. It is a convex optimization problem, an instance of a second order cone program as considered in [START_REF] Donoho | Stable Recovery of Sparse Overcomplete Representations in the Presence of Noise[END_REF], [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF]. Many approaches are now standard to solve it, including the log-barrier algorithm implemented in [START_REF] Candès | ℓ 1 -MAGIC: Recovery of Sparse Signals via Convex Programming[END_REF], interior-point algorithms as proposed in [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] or Iteratively-Reweighted Least Squares (see Appendix B of [START_REF] Donoho | Stable Recovery of Sparse Overcomplete Representations in the Presence of Noise[END_REF]). Also, using results from convex optimization, it is shown that the problem ℓ 1 -app is related to the LASSO problem (minimization of the error with a ℓ 1 -norm constraint) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and to the Basis Pursuit Denoising (BPDN) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]:

min ρ λ ρ 1 + 1 2 O{ρ} -b 2 2 . (7)
Namely, solutions of (6) are either 0 or a minimizer of (7) for some λ ≥ 0 [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF]. As a consequence, any method that can find a solution to the BPDN problem of eq. ( 7) can be adapted for solving our original ℓ 1 -app problem. For instance, homotopy methods originally proposed for the LASSO problem were considered in [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF], [START_REF] Tsaig | Sparse solution of underdetermined linear systems: Algorithms and applications[END_REF]. Various modern instances of (Orthogonal) Matching Pursuit were proposed also to solve BPDN, e.g. [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF], [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], as well as iterative Shrinkage and Thresholding methods, e.g. [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems[END_REF], [START_REF] Daubechies | Accelerated Projected Gradient Method for Linear Inverse Problems with Sparsity Constraints[END_REF], [START_REF] Blumensath | Iterative Thresholding for Sparse Approximations[END_REF].

Given the tremendous amount of works done on the algorithmic part of BP and BPDN problems, it was not possible to compare the behaviours of all the algorithms in our context. We limited ourselves to two approaches for the solving of ℓ 1 -app with Fourier domain measurements:

• The l1qc_logbarrier.m software of the ℓ 1 -MAGIC toolbox [START_REF] Candès | ℓ 1 -MAGIC: Recovery of Sparse Signals via Convex Programming[END_REF] that looks for a solution of ℓ 1 -app as a second order cone program, solved with a log-barrier algorithm. • The GPSR_BB.m software of the GPSR toolbox [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF] as an instance of an efficient gradient projection algorithm that solves BPDN (expressed as a bound-constrained quadratic program) for a sequence of values of λ. This allows a reduction of the computational time to obtain solutions of ℓ 1 -app, at the expense of not knowing a priori how to select λ. When properly selecting λ, the results obtained for the TFD built from solving ℓ 1 -app, as formulated in eq. ( 4), were similar. Such a proper choice can be made by following the guidelines reported in [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF], i.e., λ = α F -1 {A x } ∞ , with α ≈ 0.1. In practice, it happened that a smaller value α ≈ 0.02 gave the best results for the considered examples. The question of finding in an automatic, data-driven way, the best regularization parameter and, more generally, the best algorithm for the specific context of TF analysis (especially in terms of trade-off between computational cost and localization), is left open to future works.

Fig. 1 .

 1 Fig. 1. Example with a 128 points 2-component signal -Different TFDs are displayed in the case of a 128 points signal whose TF model (1) is given in the top row, in between the Wigner-Ville Distribution (WVD) and the reassigned spectrogram (RSP) that has been "optimized" as explained in Sect. IV-A. The corresponding total ambiguity function (AF) of size 128 × 128 (defined as the 2D Fourier transform of the WVD) is given on the left of the top row, with its restriction to a domain Ω consisting of 13 × 13 samples neighbouring the origin of the plane on the left of the bottom row. The rest of the bottom row presents the results obtained from optimizations based on the only knowledge of this restriction of the AF, with the exact minimum ℓ 1 -norm solution according to[START_REF] Flandrin | Time-frequency and chirps[END_REF] in between the minimum ℓ 2 -norm solution and the approximate minimum ℓ 1 -norm solution according to (4) with ǫ = 0.05 x 2 . For all diagrams, amplitudes are grey coded logarithmically, with a dynamic range of 18 dB.

Fig. 2 .

 2 Fig.2. Rényi entropy and ℓ 1 -distance -Using the same signal as in Figure1, the localization properties and proximity from the model of the different TFDs are quantified in terms of the Rényi entropy of order 3 (left, with the model entropy in thick black line) and ℓ 1 -distance (right), as a function of the relative number of AF samples card(Ω)/N 2 used in the optimizations.

Fig. 4 .

 4 Fig.4. Influence of the anisotropy of Ω -Given a domain Ω (rectangular or elliptic) whose area is kept fixed (and equal to 169), the ℓ 1 distances of the resulting TFDs to the actual model are plotted as a function of the aspect ratio "Doppler semi-axis/delay semi-axis" of Ω in the ambiguity plane.

Fig. 5 .

 5 Fig. 5. Crossing chirps example -This Figure compares, for a signal composed of the 2 linear crossing chirps represented in the model (top left), its Wigner-Ville Distribution (top right), its reassigned spectrogram (bottom left) and the approximate minimum ℓ 1 -norm solution (4) with ǫ = 0.05 x 2 (bottom right). The number of data points is N = 128 and the optimization is based on the knowledge of the 9×9 Fourier samples of the WVD neighbouring the origin of the AF plane, i.e, card(Ω) ≃ 0.005 N 2 . For all diagrams, amplitudes are grey coded logarithmically, with a dynamic range of 18 dB.

Fig. 6 .

 6 Fig. 6. Noisy 3-components example -This Figure compares, for a signal composed of 2 chirps and 1 Gabor logon represented in the model (top left), its Wigner-Ville Distribution (top right), its reassigned spectrogram (bottom left) and the approximate minimum ℓ 1 -norm solution (4) with ǫ = 0.05 x 2 (bottom right), when the observation is the signal embedded in white Gaussian noise with SNR = 10 dB. The number of data points is N = 256 and the optimization is based on the knowledge of the 14 × 14 Fourier samples of the WVD neighbouring the origin of the AF plane, i.e, card(Ω) ≃ 0.0035 N 2 . For all diagrams, amplitudes are grey coded logarithmically, with a dynamic range of 18 dB.

Fig. 7 .

 7 Fig. 7. Real data example -This Figure compares a reassigned spectrogram (left column) and the approximate minimum ℓ 1 -norm solution (4) with ǫ = 0.05 x 2 (right column) in the case of a bat echolocation call of effective length N = 400, the optimization being based on the knowledge of the 23 × 23 Fourier samples of the WVD neighbouring the origin of the AF plane. The bottom row displays enlarged versions of the distributions within the boxes in the top row. For all diagrams, amplitudes are grey coded logarithmically, with a dynamic range of 18 dB.

Table :

 :

		model RSP ℓ 1 -app
	ℓ 1 -distance to model	0	1.016 1.164
	Rényi entropy	7.435 7.769 7.853

Table :

 :

		model RSP ℓ 1 -app
	ℓ 1 -distance to model	0	1.599 1.554
	Rényi entropy	8.065 9.298 9.406
	As a complement to the synthetic examples of Figures 1,
	5 and 6, Figure


Azadeh Moghtaderi

Patrick Flandrin

Pierre Borgnat

TIME-VARYING SPECTRUM ESTIMATION OF UNIFORMLY MODULATED PROCESSES BY MEANS OF SURROGATE DATA AND EMPIRICAL MODE DECOMPOSITION

Keywords: Bandlimited signals, spectral analysis, stochastic processes, time-frequency analysis

We propose a new estimate of the time-varying spectra of uniformly modulated processes. The estimate is based on a resampling scheme which incorporates empirical mode decompositions and surrogate data techniques. The performance of the method is studied via simulations.

INTRODUCTION

Estimation of the time-varying spectra (TVS) of nonstationary processes remains a challenging problem. This is because (i) many different types of nonstationarity exist and (ii) in contrast with the stationary ergodic case, time averaging cannot be used as a substitute for ensemble averaging to reduce fluctuations in estimates. Although no "universal solution" is expected to exist for all types of nonstationarity, specific approaches can be developed for restricted classes.

In this paper, we focus on the class of uniformly modulated processes (UMPs). Theoretically, UMPs have a simple mathematical form, since they can be decomposed into a stationary process multiplied by a modulating function. More pragmatically, UMPs have been used successfully as models of various real-world processes, e.g., in seismology (see [START_REF] Herbst | Periodogram analysis and variance fluctuation[END_REF], [START_REF] Fujita | On a model of earthquake ground motions for response analysis and some examples of analysis through experiment[END_REF] and [START_REF] Dargahi-Noubary | Special ratio discriminant and information theory[END_REF]). We propose here a method for estimating the TVS of a UMP X X X. The method combines an estimate of the TVS obtained by BCMOTIFS estimator introduced in [START_REF] Moghtaderi | Multitaper Methods for Time-Frequency Spectrum Estimation and Unaliasing of Harmonic Frequencies[END_REF] with surrogate data techniques, and proceeds in three steps. The first step estimates the modulating function of X X X via an empirical mode decomposition. The second step determines the "stationary component" of X X X; surrogate data techniques are then used to create additional "virtual" realizations of X X X. These realizations are then used to compute a more accurate estimate of the TVS of X X X which is based on the arithmetic averaging of BCMOTIFS estimates for each realization.

In Section 2, we briefly review the UMPs, empirical mode decompositions, and surrogate data techniques. In Sections 3 and 4, we describe in detail the estimation method sketched in the previous paragraph. Finally, in Section 5, we examine the performance of the method on simulated data.

PRELIMINARIES

Uniformly modulated processes

Let X X X = {X t } t∈Z be a discrete-time, R-valued, zero-mean, nonstationary process. We say that X X X is a uniformly modulated process (UMP) if there exists a zeromean stationary process Y Y Y = {Y t } t∈Z with spectrum S Y Y Y and a sequence {C t } t∈Z of positive real numbers, such that X t = C t Y t for each t ∈ Z. We refer to C t as the modulating function of X X X, and Y Y Y is the stationary component of X X X. We assume that the modulating function is "slowly-varying" in comparison with Y Y Y More precisely, the modulating function is a bandlimited baseband process while Y Y Y is a broadband process without any low-frequency oscillations.

The time-varying spectrum (TVS) [START_REF] Mélard | Contributions to evolutionary spectral theory[END_REF]

of a UMP X X X is de- fined by T X X X (t, f ) C 2 t S Y Y Y (f ) for (t, f ) ∈ Z × [-1/2, 1/2].

Surrogate data techniques

Surrogate data techniques comprise a type of resampling technique. Given a single realization of a stochastic process Y Y Y , surrogate data techniques can produce an arbitrary number of "virtual" realizations of Y Y Y , called surrogates, with similar statistical properties. Originally, surrogates were used to test for nonlinearity [START_REF] Theiler | Testing for nonlinearity in time series: the method of surrogate data[END_REF]; more recently, they were used in [START_REF] Xiao | Testing stationarity with timefrequency surrogates[END_REF][START_REF] Xiao | Testing stationarity with surrogates-a one-class SVM approach[END_REF] to test for nonstationarity. In the latter, the rationale is that, given the same global empirical spectrum, a nonstationary process differs from a stationary one by some structure in time which carries over to the "phase" of the spectrum. Randomizing the phase and keeping the magnitude unchanged leads therefore to a "stationarized" process, while many other realizations can be obtained due to randomization of the phase. In a nonstationary context, this allows the construction of a statistical reference corresponding to the null hypothesis of stationarity. However if the process under study is stationary-as it will be assumed in the following-surrogates can be viewed as virtual realizations.

Empirical mode decompositions

Let S t be an arbitrary signal. Briefly, the empirical mode decomposition (EMD) [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] is a model-free and fully data-driven technique which decomposes S t , comprised of superimposed oscillations into its zero-mean oscillatory components. This is achieved by a "fine-to-coarse" recursive scheme: The fastest local oscillations (identified through neighbouring local extrema) are subtracted from the signal, yielding a residual signal to which the same procedure can be applied. Extracted in this way, each of the components is referred to as an intrinsic mode function (IMF) of S t . The recursion stops when the residual signal has no more oscillation. Denoting IMFs by

M (i)
t for 1 ≤ i ≤ I max , we write

S t = Imax i=1 M (i) t + ρ t , (1)
where ρ t is the residual signal. For MATLAB code and further details concerning implementation of the EMD, see [START_REF] Rilling | On empirical mode decomposition and its algorithms[END_REF].

TIME-VARYING SPECTRUM ESTIMATION BY MEANS OF SURROGATE DATA

Let X X X be a UMP with modulating function C t and stationary component Y Y Y . In this section, we describe a general technique to estimate T X X X , based on the use of surrogate data. Let X = X 0 , X 1 , . . . , X N -1 be a realization of X X X, and let C t be a nonzero estimate of C t for each t. Set Ỹt = X t / C t for each t. Provided that each C t estimates C t accurately (see Section 4 for details), Ỹ = Ỹ0 , Ỹ1 , . . . , ỸN-1 can be regarded as a realization of Y Y Y . We use Ỹ to obtain J surrogate data sets

Ỹ j = Ỹ j 0 , Ỹ j 1 , . . . , Ỹ j N -1 , where 1 ≤ j ≤ J. Define X j = C 0 Ỹ j 0 , C 1 Ỹ j 1 , . . . , C N -1 Ỹ j N -1
, where again 1 ≤ j ≤ J. We call X j the jth nonstationary surrogate data obtained from X and regard this data as elements of the ensemble of X X X. For each 1 ≤ j ≤ J, let T j X X X be an estimate of T X X X based solely on X j . We take the arithmetic average of the T j X X X as an estimate of T X X X , which we call the averaged surrogate-based estimate:

T av X X X (t, f ) 1 J J j=1 T j X X X (t, f ), (t, f ) ∈ Z × [-1/2, 1/2].
In this paper, we choose each T j X X X to be the BCMOTIFS estimator [START_REF] Moghtaderi | Multitaper Methods for Time-Frequency Spectrum Estimation and Unaliasing of Harmonic Frequencies[END_REF][START_REF] Moghtaderi | Evolutionary spectrum estimation for uniformly modulated processes with improved frequency resolution[END_REF]. The rationale behind this choice is as follows: For general processes, theoretical and simulation results in [START_REF] Moghtaderi | Multitaper Methods for Time-Frequency Spectrum Estimation and Unaliasing of Harmonic Frequencies[END_REF] indicate BCMOTIFS has lower bias and variance than other estimators of the TVS, especially near the boundaries of the time-frequency region. In particular, for UMPs, BCMOTIFS is approximately unbiased [START_REF] Moghtaderi | Multitaper Methods for Time-Frequency Spectrum Estimation and Unaliasing of Harmonic Frequencies[END_REF].

MODULATING FUNCTION ESTIMATION

Let X X X and X be as in Section 3, and assume that X t = 0 for each t. To proceed with the technique proposed in Section 3, we must estimate each C t . In this section, we describe how C t and X can be "decoupled," leading to an estimate of C t .

Define

X * t = log |X t |, so that X * t = log C t + log |Y t | for each t. Write X * = log |X 0 |, log |X 1 |, . . . , log |X N -1 |.
We can compute the EMD of X * to obtain its IMFs M (i) t , and write X * t as in Eq. ( 1) when S t is replaced by X * t . Note that the applicability of EMD method in general, does not require particular forms of oscillations and therefore appropriate for the log-transformed signal. Loosely speaking, the assumption made in Section 2.1 states that the oscillations of log C t are much slower than those in log |Y t |. As a result the oscillations of log C t are accurately described by a set of high-order IMFs such that for 1 ≤ i * ≤ I max we have

log C t ≈ Imax i=i * M (i) t + ρ t . (2)
Following above, estimating C t reduces to determining i * .

We now describe three approaches to determining i * . The energy of the ith IMF is defined by

E (i) = N -1 t=0 |M (i) t | 2
for each i. For typical broadband signals, E (i) is decreasing in i [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF]. Identifying the smallest index i + 1 such that E (i+1) > E (i) can therefore provide information about the value of i * . This is called the energy approach.

Denote the number of zero crossings of the ith IMF by Z (i) . It can be observed that, in the absence of low-frequency oscillations, broadband signals give rise to a family of IMFs satisfying R (i+1) ≈ 2, where R (i+1) = Z (i) /Z (i+1) whenever it is defined. Furthermore, we assume that the random variables R (i+1) are approximately equally distributed. We call R (i+1) the ith ratio of zero crossing numbers (ith RZCN). Hence the smallest index i + 1 for which R (i+1) is "significantly different" from 2 can provide information about the value of i * . The question is how to quantify "significantly different." To answer this question, we construct 13 processes (all approximately broadband without low-frequency oscillations), including (i) nine fractional Gaussian noise processes with Hurst exponents H = 0.1, 0.2, . . . , 0.9, (ii) two AR(2) processes, and (iii) nonstationary processes, the first being AR(2) with time-dependent coefficients and the second being frequency-modulated. For each process, we create 5000 realizations of length N = 1000, and compute the EMDs and RZCN of each. We then compute the empirical distribution of

R = [ R 1 R 2 • • • R 5000 ] where denoting the ith RZCN of the bth realization by R (i+1) b , 1 ≤ i ≤ I b max -1, we have R b = [R (2) b R (3) b • • • R (I b max ) b
] for 1 ≤ b ≤ 5000. Fig. 1 shows the empirical distribution of R computed for each process and shown in different colors. The result of our simulations encourages the idea that, regardless of the type of broadband process, the distribution of R remains unchanged. As a result of this, we may then propose a common threshold test for determining when R (i+1) is significantly different from 2. For each process, we compute the empirical distribution of R (i+1) over all the realizations and for each i = 1, 2 . . . , I max -1 respectively where here we assume I max = 10. In this case, if, as an example, for a particular realization, we have I max = 8, we assume R (8) = R (9) = R (10) and if I max = 12, we discard R (11) . We then calculate 5% and 95% significance level of each distribution, and refer to them as ith right and left thresholds respectively. Fig. 2 shows the right and left thresholds for each process (different colors). Any RZCN which is outside of the two standard deviation of the averaged ith left and right thresholds within 13 process is significantly different from 2. Finally, as mentioned earlier, the smallest index i+1 where R (i+1) is significantly different from 2 determines i * . The problem is that since the selection of the left and right thresholds are entirely based on an empirical result, it is always possible that the smallest i + 1 is a false detection and not the correct i * . This approach is called the ratio approach.

The energy and ratio approaches can be combined in order to reduce the number of false detects as follows: For each 1 ≤ i ≤ I max -1, we compute each index i + 1 such that E (i+1) > E (i) . We also evaluate every index i + 1 where R (i+1) is significantly different from 2. We then choose i * to be the smallest common index in both approaches. The value of i * is then plugged into Eq. ( 2) in order to evaluate log C t . The estimated modulating function is C t = exp(log C t ). This combined approach is called the energy-ratio approach.

NUMERICAL EXAMPLES

Let X X X (k) = {X (k) t } t∈Z , 1 ≤ k ≤ 5 be UMPs with modulat- ing functions C (1) t = e - (t-500) 2 2(200) 2 , C (2) t
= 1 + t/T C (3) t = 2 + sin(2πpt), C (4) t
= 1.5 + cos(2πqt) C (5) t = 2 -e - (t-500) 2 2(200) 2 ,
where T = 200, p = 0.002 and q = 0.001. The stationary components of the X X X (k) are as follows: Y Y Y (1) , Y Y Y (2) , and Y Y Y (3) are fractional Gaussian noise with Hurst parameters H = 0.2, H = 0.5, and H = 0.8, respectively, and Y Y Y (4) and Y Y Y (5) are the AR(2) processes

Y (4) t = 0.2Y (4)
t-1 + 0.5Y

t-2 + ζ t Y (5) t = 0.8Y (5) t-1 -0.4Y (5) t-2 + ǫ t . (4)
Here, {ζ t } t∈Z and {ǫ t } t∈Z are independent white noise processes with variance 10 4 .

For each k, we create 5000 realizations of length N = 1000 of X X X (k) and compute the EMD for each realization. We then apply the energy, ratio, and energy-ratio approaches to each realization and each UMP to evaluate i * , denoted i * ,a (k), where a = 1, 2, 3 indicate energy, ratio and energyratio approaches respectively. Using our theoretical knowledge of the modulating function for each UMP, we pick i * which gives the minimum L 2 -distance between log C t from Eq. ( 2) and theoretical log modulating function. We denote i * obtained from minimum L 2 -distance by i † (k). Table 1 shows the number of times in 5000 where we have obtained i * ,a (k

) = i † (k), i * ,a (k) = i † (k) + 1, i * ,a (k) = i † (k) -1, or |i * ,a (k) -i † (k)| > 1.
It is clear from Table 1 that the energy-ratio approach outperforms the other two approaches. We now take X X X (4) for further analysis. Let X = X 0 , X 1 , . . . , X N -1 be a realization of X X X (4) and X * the log-transform of X . We compute the EMD of X * to obtain its IMFs and then evaluate E (i) and R (i+1) for each i. Fig. 3 shows the plot of E (i) , 1 ≤ i ≤ I max , (top plot) and R (i+1) , 1 ≤ i ≤ I max -1, (bottom plot) where I max = 9 in this example. Applying the energy approach we determine that there are 4 indexes i + 1 where i = 3, 5, 7, 8 which satisfies E (i+1) > E (i) . On the other hand, applying the ratio approach, we find out that there are three indexes i + 1 for i = 6, 7, 8 where R (i+1) is significantly different from 2. These indexes are marked by red triangles in Fig 3 . Applying the energy-ratio approach, we see that the first common index between the two approaches is i * = 8.

Using i * = 8, we evaluate C t (see Section 4). The estimate C t is then used to create J = 50 surrogate Table 1. Comparison between energy, ratio and energy-ratio approaches using the theoretical log C t for 5 UMPs. i) , i = 1, 2, . . . , 9 (top) and R (i+1) , i = 1, 2, . . . , 8 (bottom). The red triangles mark the indexes which satisfy the conditions in energy (top) and ratio (bottom) approaches. data X j , 1 ≤ j ≤ 50 (see Section 3). Using 50 nonstationary surrogate data, we then estimate T X X X using the averaged surrogate-based estimate T av X X X (see Section 3). Fig. 4 shows the performance of this estimator in comparison with the BC-MOTIFS estimator for X .

CONCLUSION

In this paper, we proposed a new scheme to estimate the timevarying spectra of uniformly modulated processes. The estimate proceeds by computing "virtual" realizations, using surrogate data and empirical mode decompositions. Simulation results suggest that the new estimator performs well, and affords a significant improvement over BCMOTIFS.

Fig. 1 .

 1 Fig. 1. Empirical distribution of R computed for different processes (different colors). Apart from the expected peak at 2, there exist several smaller but visible peaks. They appear in the presence of high order IMFs which have small values of zero crossing numbers.

Fig. 2 .

 2 Fig. 2. The right and left thresholds (marked by stars) computed for different processes (different colors). Solid lines indicate averaged left and right thresholds and dashed lines indicate two standard deviation of the averaged thresholds.

  a(1) = i † (1) -1 583 484 891 |i * ,a(1)i † (1)| > 1 703 1672 199 i * ,a(2) = i † (2) 1866 1420 2697 i * ,a(2) = i † (2) + 1 1181 735 948 i * ,a(2) = i † (2) -1 658 495 949 |i * ,a(2)i † (2)| > 1 1295 2350 406 i * ,a(3) = i † (3) 1919 1339 2710 i * ,a(3) = i † (3) + 1 1386 902 985 i * ,a(3) = i † (3) -1 373 521 910 |i * ,a(3)i † (1)| > 1 1322 2238 395 i * ,a(4) = i † (4) 2636 1893 3316 i * ,a(4) = i † (4) + 1 1100 980 643 i * ,a(4) = i † (4) -1 232 492 740 |i * ,a(4)i † (4)| > 1 1032 1635 301 i * ,a(5) = i † (5) 2615 1487 3598 i * ,a(5) = i † (5) + 1 1056 902 637 i * ,a(5) = i † (5) -1 295 363 529 |i * ,a(5)i † (5)| > 1 1034 2248 236

Fig. 3 .

 3 Fig. 3. E(i) , i = 1, 2, . . . , 9 (top) and R (i+1) , i = 1, 2, . . . , 8 (bottom). The red triangles mark the indexes which satisfy the conditions in energy (top) and ratio (bottom) approaches.

Fig. 4 .

 4 Fig. 4. Theoretical TVS of X X X (4) (top), estimate of the TVS using X (middle), estimate of the TVS using X j , 1 ≤ 1 ≤ 50 (bottom). The parameters used in BCMOTIFS estimator are N W = 4, K = 7 and B = 201.


Fanny Lignon
email: fanny.lignon@univ-lyon1.fr

TEACH THEM ALL I -Introduction

In France, public schools are mixed. But co-education of boys and girls, however, is not easy. Currently, the curricula are based on an official text published in 2006: the common core of knowledge and skills [START_REF] De Robien | Ecole et collège: tout ce que nos enfants doivent savoir. Le socle commun des connaissances et des compétences[END_REF]. It urges teachers to teach, among other things, respect for others (civility, tolerance, rejection of prejudice and stereotypes) and respect for the other sex. The theme is important, but, in the opinion of most teachers, difficult to address in class.

During the twentieth century, using games became a teaching method. However, if teaching through play is now a common and legitimate practice, teaching the game is still in its infancy. Nevertheless, the common core of knowledge and skills has required from teachers that they educate children to understand the media and teach using media and technology.

An idea: bring together two difficulties to solve a problem! Educating gender through video games, and vice versa! I will present my work in two stages:

1) An analysis of gender discourse conveyed by a certain category of video games.

2) Some leads to work from this category of games, and with pupils, on gender.

II -Combat video games: female images, male images

A -A specific corpus Among all the available video games, I chose to work on action games, fighting games. I chose these games because they offer the player to embody human avatars and because their reputation puts them a priori as the antithesis of all learning. A kind of bet! I worked from all the demo CD PS1, PS2 and PS3 published in France between 1995 and 2009 and examined 44 fighting games.

B -My method

1) References

My work is a continuation of that initiated in 2000 by J. Cassell and H. Jenkins [START_REF] Cassell | From Barbie to Mortal Kombat: gender and computer games[END_REF]. It also extends to research of M. Marriott [START_REF] Marriott | Fighting women enter the arena, no holds barred[END_REF] and M. Consalvo [START_REF] Consalvo | Hot dates and fary-tales romances: studying sexuality in video games[END_REF].

2) Quantitive analysis of the number of men and women

Combat games belong to the category of simulation games. With the exception of wresting games, it is to be noticed that they only propose male characters. Those belonging to the category of arcade games always let the player the possibility to choose to embody a man or a woman.The more recent the game is, the greater choice one has. 27 % of the characters directed in the games published for PS1 (1995)(1996)(1997)(1998)(1999)(2000) were women. The number has grown with PS2 (2000-2005) and PS3 (2005)(2006)(2007)(2008)(2009), reaching 33%. This phenomenon can be internally observed within the series themselves.

3) Analysis of the identity cards of the gaming characters

At first sight a kind of equality seems to be the rule concerning the onomastic treatment of feminine and masculine characters. However, the first names chosen, whether real or imaginary, always point out to the character's gender.

The characters' identification files mention almost in every case the protagonists' age. For women, the figure is systematically between 19 and 22. In some cases, the mention « no data » seems to pertain to some virtual coquetry. In combat games the girls are always young. You will never see a mature woman fighter or an old woman fighter. The average age is comparatively higher for men (26). Some fighters are bald, some are going grey at the temples, others are white-haired. Ageing seems to be a male privilege. Size, weight and measurements of protagonists, boys and girls, are also mentioned. The figures look very much like ideal runway models. These allegedly perfect proportions are perhaps suitable for a model but certainly do not correspond to a woman fighter's physique. Besides, from a sporting viewpoint, only gender and weight are taken into account. What good is it to know the waist, the hip and the bust measurements of one's adversary? One may object that this data is also present in the male characters' files. But why bother about waist, hip and chest measurements for a sportsman! We can suppose that if men are submitted to the same treatment as women, it is above all in order to legitimate the presence of this information concerning women… Last but not least: combat style. If one refers to their files, female characters are attributed styles using the redirection of energy or guile while the male characters are attributed more 'brute force' style. But beyond this interpretation, it is obvious that women practice martial arts which are, with a few exceptions, more obscure than those practiced by male characters. As if girls were unable to win a contest when practicing an acknowledged style. As if they had to turn to different, rare or strange styles.

4) Audiovisual analysis of body types, clothing, gestures and attitudes. a) Virile men

The combat game heroes' common bodies are athletic ones. All of them are conceived more or less on the same developed pattern: broad-shouldered, with protruding pectoral, abdominal muscles and huge biceps. Clothes worn by the combatants are varied according to shape and fabric. They are set in all tints except pink. Besides traditional combat dress one finds devices meant to protect the body or facilitate its movements. Jackets are very often sleeveless. They are worn open and next to the skin, which allows the combatant to display torso and arms.

But what about the way bodies move? We have observed two types of gestures: sporting, coherent, generally well led gestures and others, improbable or too acrobatic to be realistic. In either case, our heroes are endowed with all the qualities which make good combatants: steady positions, powerfulness, an excellent spring.

In combat games, whether in terms of physique or dress, everything seems to be called up to valorize men's bodies and assert their virility.

b) Mannish women

Men being caricatures of virility, one expects women in the same way to be caricatures of femininity, therefore to wear no element connoting masculinity whatsoever. This is not quite the case. Even though they definitely appear as women, they take on the same objectives as their male counterparts, behave in the same way when confronted with problems, and solve these problems in the same way. They give blows, they receive some. As if a woman were merely a man in a woman's body.

c) Feminine men

Considering what we have just written, it seems uneasy to evoke femininity, even though stereotyped, regarding the video games combatants. Considering the size of their wardrobes, shouldn't we consider their attention to good looks as surprising? Some recent games display a « creation » mode which allows modification of the hero's physical appearance and provides a choice of garments one by one, from shoes to hat and underwear. The range of choice is equivalent for men and women. We must note here that combat games, under cover of fight, allow boys to play with their dolls.

Soul Calibur III. (Namco. 2005). Screenshot. d) Feminine women

The virtual combatants' bodies display shapes in relation with their measurements. Hips, breasts and buttocks they do have. Sexual ambiguity is not the matter. All are built according to an identical canon: longlegged and wasp-waisted. Only faces change, somewhat the way Barbie dolls' faces do. Combat games offer thus a unique and simplifying model to contemplate. In this, they are far from being innovating! The women combatants' wardrobe is rather well furnished, from dresses and skirts open to the waist, slacks and shorts peculiarly tight, tank tops and close-fitting bras. Colours are generally showy, shapes often eccentric and the matching frequently daring.

All these garments have one thing in common. They are both sexy and definitely feminine. In combat games, unisex dresses are banned. Shapes mould bodies and valorize them. They even show, locally, nude skin, arms, legs, navel.

These garments are otherwise typically unsuitable. Unlike men, women are not favored with sportswear. They even wear dress rather unfit with the activity they are required to assume. For instance, how adequate is practising martial arts with a long dress? Previously, the female characters were, in fighting games, swift, supple and weak. Today, they are swift, supple and powerful and can fight as equals with their male counterparts. However, everything concurs to bias the player towards contemplating the sportswomen' movements and those movements are, to a great extent, alien to sport. Attitudes, movements and gestures put on stage are often suggestive, if not licentious. You go at once from true sport to sport in chamber. Suppleness and acrobatic prowess then assume another dimension, that of an erotic promise. Indeed, chasing the glimpses has long pertained to the pleasure of these games. Thus, combat games put on stage erotically characterized gestures involved by sport practices. Generous breasts delineate a movement; an open skirt occasionally unveils a glimpse of an upper thigh; panties are glimpsed through somersault. Heroines unwillingly commit unchaste and exciting gestures. Men and women in these games obviously are not in the same boat. Abdominal muscles excepted, a male character's advantages will never be valorized. This is all that Namco's team, according to Zoe Flower, calls "functional beauty in combat" [START_REF] Flower | Getting the girl: the myths, misconceptions, and misdemeanors of female in games[END_REF]!

C -My conclusion

As a matter of fact, men and women I encounter in my study are caricatures. But can we speak of sexism for all that? On the woman side, I answer this question in the affirmative mainly because of the erotic connotations borne solely by the heroines' gestures. Contrary to men, women have been conceived, so it seems, in order to please and only for that purpose. The player is literally invited to manipulate splendid women as sex objects with his joystick. In the case in point we are exactly confronted with discrimination. Nevertheless, to rest on firmer ground, we have to verify if by chance one can find some egalitarian hints... Men and women meet, thus defining a place where they practice together an activity traditionally devoted to men. « Shedding blood » is not a male privilege any more [START_REF] Héritier | Le sang du guerrier et le sang des femmes[END_REF]. Equality manifests itself thus in opposition. We will not reexamine the virtual women' sporting competences, neither will we their evolution. Let us merely note that they proceeded from weakness to vigor, and that they are from now on respectable and fearsome opponents. This purported « equality » however remains somewhat contrived. From weak to wonder woman, it seems that games creators have gone too fast and too far. Maybe combat games aim at some equality, but they only succeed in putting on stage a certain kind of egalitarianism by negating the differences between men and women.

Finally we must resort to an ultimate parameter: the player. Combat games let girls as well as boys adopt one gender or the other, according to their mood. Each game provides then the opportunity to live up to a punctual transvestism, so much so that the caricatural aspect of male and female representations becomes relative. It even becomes possible to explore the idea of equality, all the more so as one of the purposes of such games is to learn how to master each character's special blow to win. To master the choice of gender is not without consequence, and this may be the site of the contestation of the man-woman bipolarity

III -Ideas for the classroom

The fighting games provide material to work on prejudice and gender stereotypes. Why not use them to learn tolerance and respect for the opposite sex? This proposal is a priori quite paradoxical, but not as paradoxical as it may sound. Let us do an experiment. Take a girl, put her in command of a male. Take a boy, put him in command of a female. Let the mixture « brew », then interview! Predict that the debate on gender stereotypes will be richer after this experience than before. For it is one thing to identify and analyze these stereotypes, it is another to live them from within. Many other options are available to teachers:

 Analyze the physical, psychological and social virtual characters set in different types of video games.  Compare characters' evolution and transformation in other media as well (press, television, cinema).  Identify in other cultures, the various elements determining the feminine and the masculine: in other words how the feminine and masculine components emerge and combine, elsewhere.  Create virtual characters that can be either stereotyped, or non-stereotyped, or even put stereotypes into question.  Imagine interactive stories around these characters.  Etc...

IV -Conclusion

Teachers in the XXIth century have to learn to consider action video games and their popularity not as necessary evils but as welcome teaching aids. Go from « Beat them all » to « Teach them all »!

Tekken 4 .

 4 (Namco, 2002). Screenshot.