Florent Berthaut
email: florent@hitmuri.netlukedahl@virginia.edu

Luke Dahl

BOEUF : a Unified Framework for Modeling and Designing Digital Orchestras

Keywords: Boeuf, orchestra, collaboration, digital musical instrument, CVE

Orchestras of Digital Musical Instruments (DMIs) enable new musical collaboration possibilities, extending those of acoustic and electric orchestras. However the creation and development of these orchestras remain constrained. In fact, each new musical collaboration system or orchestra piece relies on a fixed number of musicians, a fixed set of instruments (often only one) and a fixed subset of possible modes of collaboration. In this paper, we describe a unified framework that enables the design of Digital Orchestras with potentially different DMIs and an expandable set of collaboration modes. It relies on research done on analysis and classification of traditional and digital orchestras, on research in Collaborative Virtual Environments, and on interviews of musicians and composers. Implementing this framework as a software library will enable advanced musical collaboration modes to be used in any orchestra, even spontaneous ones.

Introduction

Orchestras of Digital Musical Instruments (DMIs) began to appear at the end of the 1970s, with the League of Automatic composer and later The Hub [START_REF] Gresham-Lancaster | The aesthetics and history of the hub: The effects of changing technology on network computer music[END_REF]. With the subsequent generalization of hardware platforms, music software and communication protocols, the number of orchestras has grown. A Digital Orchestra (DO) can be composed of a single multi-user DMI [START_REF] Jordà | The reactable*[END_REF], or of multiple DMIs of the same type, such as in most laptop [START_REF] Wang | Stanford laptop orchestra(slork)[END_REF] or mobile phones [START_REF] Oh | Evolving the mobile phone orchestra[END_REF] orchestras, or of a set of DMIs of different types. In turn, the DMIs may rely on various hardware interaction devices, feedback modalities, software architectures and sound synthesis techniques. DOs enable new modes of collaboration that were not possible with acoustic or electric orchestras. For example, musicians can share audio or control streams, thus processing the output of one instrument in another instrument. They may exchange precomposed building blocks, such as patterns or sounds. Or a musician might even allow their instrument to be controlled by another musician. Finally, collaboration modes present in traditional ensembles, such as non-verbal communication for coordinating musical cohesion or variation [START_REF] Seddon | A comparison of modes of communication between members of a string quartet and a jazz sextet[END_REF], can be augmented in DOs.

However, the design and implementation of Digital Orchestras remains complex, especially when they include diverse musical instruments. In fact, most existing orchestras rely on frameworks specific to a software or hardware platform or even to a single instrument. Even if standard communication protocols such as MIDI or OpenSoundControl are used, each orchestra manages communication between instruments in a specific way, thus prohibiting new instruments from easily joining. Thus, orchestras of mixed instruments are more complex to build, and, as in the case of spontaneous jam sessions, advanced collaboration possibilities may be completely inaccessible. Often these orchestras rely on a limited set of collaboration possibilities, which are reimplemented for each new orchestra or musical piece.

These practices result in a number of impediments to the exploration of musical collaborations using DMIs. Therefore the computer music community would benefit from a framework that is simple enough to be integrated into any instrument, that takes into account existing modes of collaboration, and that can be expanded by adding new collaboration modes. Such a framework would facilitate the creation of orchestras and encourage exploration of new collaboration modes. The contributions presented in this paper are twofold:

1) We present a classification of musical modes of collaboration based on a literature review of research on DOs and Collaborative Virtual Environments (CVE), and on interviews and discussions we conducted with composers and musicians in DOs at the SCRIME at University of Bordeaux, and at CCRMA at Stanford University. This classification, described in section 2.1, allows for the practical analysis of musical collaboration and can be extended to include novel collaboration modes.

2) We then provide a set of components that compose a Digital Orchestra. These components, described in section 2.2, allow for the design of orchestras with access to any mode of collaboration in our classification, as demonstrated in section 2.3.

These form the conceptual basis of the BOEUF framework. Their implementation in the form of a software library, and the integration of this library in various software components will enable advanced musical collaboration in DOs with any set of instruments, even in the context of improvised sessions.

Related work

A number of systems exist for describing orchestras of DMIs and collaborative instruments. Jorda [START_REF] Jordà | Multi-user instruments: models, examples and promises[END_REF] classifies multi-user instruments according to the number of users and whether this can be variable, the flexibility of roles, and whether users can influence others' musical output. Blaine and Fels [START_REF] Blaine | Contexts of collaborative musical experiences[END_REF] describe collaborative interfaces according to their use (focus), properties of the instrument (location, media, physical interface, musical range, level of physicality, and pathway to expert performance), and the structure of interactions between players (scalability, player interaction, directed interaction.) Hattwick and Wanderley [START_REF] Hattwick | A dimension space for evaluating collaborative musical performance systems[END_REF] create a space of collaborative musical systems, with dimensions of texture, equality of roles, centralization of information, the role of physical location, whether interactions are time synchronous, and whether sound production depends on more than one performer. Our framework is intended to span this space, with the exception that ours is intended to enable only synchronous (realtime) interactions. Weinberg [START_REF] Weinberg | Interconnected musical networks: Toward a theoretical framework[END_REF] provides a historical overview of interconnected musical networks. Weinberg's examples, as well as his descriptions of intra-player interdependencies were influential on our development of the BOEUF framework.

Several protocols and software tools have been created to deal with the sharing of musical data for both single instruments and within networked orchestras. For example, Minuit, Jamoma and libMapper all give access to the structure and parameters of networked DMIs, sometimes with features for watching and grabbing parameters. An interesting example is the Digital Orchestra Toolbox [START_REF] Malloch | A network-based framework for collaborative development and performance of digital musical instruments[END_REF] which simplifies the collaborative creation and mapping of DMIs. Most of these tools in turn rely on the OpenSoundControl protocol for network communication. However, while they provide all the generic sharing and mapping features required for networked musical control, these tools do not specifically cover the modes of collaboration used in DOs, and thus fail to provide a common basis for creating orchestras of mixed DMIs.

The BOEUF Framework

The conceptual part of the BOEUF3 framework consists of a set of components that allow for designing any digital orchestra, and enables a set of collaboration modes, which are organised into three categories. Our framework draws on research done on CVEs. A Collaborative Virtual Environment is defined in [START_REF] Snowdon | Collaborative Virtual Environments: Digital Places and Spaces for Interaction[END_REF] as a "computer-based, distributed, virtual space or set of places. In such places, people can meet and interact with others, with agents or with virtual objects". Challenges of CVEs are very similar to those of DOs. As with orchestras, users in CVEs need to perceive each other and to communicate within the virtual environment in order to cooperate for specific tasks. However, unlike research on musical collaboration, CVE research has gone further than merely classifying collaborative applications. Practical implementation-oriented models have been proposed, as have frameworks that aim at helping developers build CVEs. Adapting CVE frameworks and models to the musical domain provides us with insights that inform our own practical collaboration framework.

Modes of collaboration

From the existing work presented above and our interviews, we extract three modes of collaboration between musicians: Cooperation, Communication and Organisation. Figure 1 shows the use of these modes in a selection of orchestras and multi-user instruments and pieces. Sound Bounce [START_REF] Dahl | Sound bounce: Physical metaphors in designing mobile music performance[END_REF] is a piece for the Stanford Mobile Phone Orchestra [START_REF] Oh | Evolving the mobile phone orchestra[END_REF]

Mopho: Sound Bounce X X X X X X Anahata X X X X X Les complementaires X X X X X X Couacs X X X X Reactable X X X X X X X LOLC X X X X X The Hub: X X X Slork: Intellectual Improperty X X X X X X JamiOki X X X X Cobra X X X X X X X
Fig. 1. Modes of collaboration in various orchestras and multi-user instruments electro-acoustic trios at the SCRIME. Couacs [START_REF] Berthaut | First Person Shooters as Collaborative Multiprocess Instruments[END_REF] is a multi-user instrument in the form a 3D first person shooter video game. The Reactable is a multi-user instrument based on a tabletop tangible user interface. LOLC [START_REF] Lee | Collaborative musical improvisation in a laptop ensemble with lolc[END_REF] is an orchestra that relies on live-coding and instant messaging. The Hub is one of the first Digital Orchestras. Slork: Intellectual Improperty is a piece for the Stanford Laptop Orchestra [START_REF] Wang | Stanford laptop orchestra(slork)[END_REF]. JamiOki is a system for playing game pieces, providing instructions to each musician and getting feedback from them. Cobra is a game piece with various collaboration modes.

Cooperation modes describe the coordination of musicians' actions with respect to their instruments. We define three subcategories that correspond to different possible interconnections between musicians and their instruments. These subcategories are inspired by the CVE cooperation framework described in [START_REF] Margery | A general framework for cooperative manipulation in virtual environments[END_REF]. Cooperation modes can therefore be: Independent when multiple musicians do not control the musical output of the same instrument or module. Complementary when two or more musicians can affect the same musical output but at different levels of the audio synthesis graph, i.e. each musician controls a different sound parameter. Concurrent when musicians can affect the same musical output at the same level, i.e. when multiple musicians modify the same musical parameter. Independent modes of cooperation exist whenever two musicians play at the same time. The cooperation in this case consists of the coordination of gestures as each musician performs their own instrument. Complementary modes of cooperation are used in many orchestras. Our interviews with musicians from Les Complementaires and Anahata revealed that they were both using a non-mediated complementary mode of cooperation, by spreading one sound over several musicians, one playing the attack, another the sustained part and another the end of the sound. Similarly, in the Reactable, complementary cooperation occurs when two musicians manipulate different tangibles on the same audio path. Finally, concurrent modes of cooperation are less common as they imply either conflicts or games between musicians. In Couacs [START_REF] Berthaut | First Person Shooters as Collaborative Multiprocess Instruments[END_REF] for example, musicians control avatars in a musical video game and can override the parameters of other musicians by shooting at them. Concurrency handling strategies must be applied, as discussed for example in [START_REF] Fencott | Hey Man, you're invading my Personal Space! Privacy and Awareness in Collaborative Music[END_REF], such as applying the newest action, averaging between actions of different musicians, using a physical model with different weights, or grabbing a parameter for exclusive use.

Communication modes do not directly impact the production of sound, but rather influence the actions of the musicians, in particular those who are involved in one of the cooperation modes. In most orchestras, communication modes are non-mediated and enable musical cohesion or variation [START_REF] Seddon | A comparison of modes of communication between members of a string quartet and a jazz sextet[END_REF]. However DOs often lack visibility [START_REF] Dahl | Wicked problems and design considerations in composing for laptop orchestra[END_REF], making it difficult for performers to see or understand the state of others' activities. This issue is amplified in the case of networked orchestras when musicians are not physically collocated. These difficulties can be addressed by integrating communication modes into the framework. In the field of CVE, similar communication problems have led to the use of concepts such as Embodiment, i.e. using avatars to represent users, CoPresence i.e. the feeling users have of being together in the virtual environment, and Awareness, i.e. the understanding of other users' actions [START_REF] Benford | Managing mutual awareness in collaborative virtual environments[END_REF]. This concept of awareness has also been discussed for the case of DOs in [START_REF] Fencott | Hey Man, you're invading my Personal Space! Privacy and Awareness in Collaborative Music[END_REF]. Modes of communication are divided into three subcategories. Awareness includes all non-intentional communication, such as making musicians' activities visible to each other for the purpose of enabling synchronisation, cohesion or variations. Indications are intentional communicative acts such as demonstrating gestures and intentions or sending commands. Exchanges include transfers between musicians of musical data. An example of awareness in Couacs [START_REF] Berthaut | First Person Shooters as Collaborative Multiprocess Instruments[END_REF] is the use of avatars that provide information on the musical output each player is generating. Awareness is especially important in mixed instruments orchestras, when it is not clear how each musician contributes to the musical output. One example of indications is the system of text messages sent by the conductor to instrumentists in the piece Intellectual Improperty. They are also used in Jamioki, to guide musicians through musical games or improvised pieces. An example of exchange can be found in Sound Bounce where musicians use ball-throwing gestures to pass a sound process from one player to another.

Organisation modes do not have any effect on the music produced but rather impact the communication and cooperation modes. We define three organisation modes: Nomination consists in defining the roles of musicians within the orchestra. Grouping consists in defining a hierarchy of groups of instruments. Selection is the act of choosing a single instrument or a group in the context of cooperation or communication, e.g. selecting which musician to send an indication to. A common example of nomination is the role of a conductor. In most orchestras, this role is fixed and the conductor has a specific interface, such as in Intellectual Improperty. However, roles can also be dynamically changed as described in [START_REF] Jordà | Multi-user instruments: models, examples and promises[END_REF]. For example, a musician from Les Complementaires explained how they use the role of soloist to give priority to one musician at certain moments of a performance. An example of grouping is found in Intellectual Improperty, where the conductor can group musicians. With the selection mode, he then chooses the group he sends a message to (indication).

Orchestra components

In order to design orchestras that enable the collaboration modes described above, our framework includes a generic model of a Digital Orchestra with a set of components : Session, Group, Instrument, Module, Parameter, Output, Meter, and Message. A session represents an instance of a DO. It contains the instruments and the network of possible interactions between instruments. A group is a set of instruments or groups. The parameters common to all instruments in the group can be grabbed and set simultaneously. Similarly messages that are sent to a group are sent to all members. For example a group can be all the musicians with the same instrument in a poly-instrument orchestra. A default group of all instruments is always defined, thereby giving access to a list of parameters common to all instruments of the session e.g. tempo or scale.

An instrument represents a bounded set of music-generating processes (i.e. modules) and a user interface. It also has parameters, outputs, meters and it receives and sends messages. A module is a software component that produces musical data, of audio or control type. It is composed of several parameters and outputs. Modules have a type, possibly from a common set of types (each with a predefined number of parameters and data). This way a module can be copied by another instrument if this instrument handles modules of the same type. For example, a module of type LowPassFilter might have Cutoff Frequency and Q as parameters, whereas a module of type MidiPattern will hold an array of MIDI events. Many instruments have an internal structure that is more complex than a simple set of modules. They may have complex audio graphs with many hierarchical levels. However, for the purpose of user interaction these can usually be flattened as a set of modules and associated parameters .

A parameter is an attribute of a module or instrument that influences its musical production. Parameters can be of various types such as MIDI events, float or integer values, input audio streams, and so on. Parameters can be retrieved, watched, indicated, set and grabbed. These actions are always accessible to the instrument that owns the parameter. Other instruments in the session can request authorization for each of these actions. Concurrent access can be managed in different ways, the simplest being by grabbing a parameter so that only one instrument can access it. An output is a musical attribute that is produced by a module or instrument and that can be retrieved and watched by another instrument without authorization. Outputs can be of the same types as parameters. A meter is a component of an instrument that is not used in the actual sound production, but rather gives an indication on the activity of the instrument e.g. spectrum or loudness. A message is a text, picture or video sent from one instrument to another instrument or group. Messages can be standard (e.g. Start, Stop, Fade Out), defined per session instruments, or dynamically created.

Designing Orchestras using BOEUF

The BOEUF framework allows for the design of orchestras with all collaboration modes described in section 2.1, using the set of components that we defined in section 2.2. Let us illustrate that process with three examples: Sound Bounce, Couacs and Intellectual Improperty. A graphical representation of their design can be seen in Figure 2. Cooperation modes are handled by the access each instrument has to the parameters and outputs of other instruments. In the case of concurrent cooperation, as explained in section 2.1, different strategies of concurrency handling can be defined by each instrument for each parameter. In Couacs concurrent cooperation involves one instrument (player) setting all the parameters of the module of another player at once by shooting this player. Communication modes are implemented through various components. Awareness is enabled by the use of meters and watched parameters and outputs. In Couacs, the players' avatars are scaled using the loudness meters of each instrument, allowing musicians to perceive who is playing what. Indications rely on the use of messages and indicated parameters. In Intellectual Improperty, messages are sent by the conductor to the musicians. Exchanges are performed by copying modules or sets of parameters. For example, in Sound Bounce, modules representing sound balls are sent from one instrument to another. Organisation modes are mostly implemented through the group and instrument components. Nomination of the conductor in Intellectual Improperty, corresponds to the definition of a group for the conductor only. The conductor can define (nominate) groups of musicians (instruments) to send indications (messages) to. An example of selection can be found in Sound Bounce, where each smartphone (instrument) selects another instrument for the exchange of sound balls (modules).

Conclusion

In this paper we introduced BOEUF, a framework for designing, building, and managing orchestras of DMIs. Our next step is to implement this framework as a software library which can be integrated into DMIs. This will facilitate the modes of collaboration described here, and empower the community to more easily explore new musical collaboration possibilities.

Another research area raised by BOEUF has to do with questions of usability when accessing BOEUF components during performance in order to enable or enact the various modes of collaboration. For example, we would like to studyhow to best enable a musician to grant another musician complementary or concurrent access to a parameter on their own instrument.

Fig. 2 .

 2 Fig. 2. Designing orchestras using our framework. Dashed lines are used for Communication, solid lines for Cooperation, and dotted lines for Organisation.

  . Anahata and Les complementaires are both

	Orchestra/Instrument	Cooperation	Communication	Organisation
		Independent Complementary Concurrent Awareness Indications Exchanges Nomination Grouping Selection

for BOEUF OrchEstras Unification Framework, boeuf meaning jam session in French

This project was partially funded through the Marie Curie FP7 framework (Grant Agreement PIEF-GA-2012-330770).


tion, called the wavelet transform, which basically has the same structure [Gr- , the frequency translations being replaced by dilations (the dilation parameter being interpreted as the quotient of a reference frequency by a frequency shift). Both methods have been applied to signal analysis, with comparable performance levels, in different contexts.

The mathematical theory of continuous wavelet analysis and Gabor analysis, and their discretizations (frames), has been mainly developed in [START_REF] Daubechies | The wavelet transform, time-frequency localisation and signal analysis[END_REF]. In particular, it has been shown that there exists an important difference between these two representations, at the level of discretization. Indeed, while it is convenient to sample uniformly the time and frequency parameters in the Gabor case, the time and scale parameters have to be sampled on a dyadic lattice in the wavelet case. This is actually the consequence of a strategic difference: Gabor functions perform an analysis with constant absolute bandwidth ∆ν, while wavelets perform an analysis with constant relative bandwidth ∆ν ν . We consider here their application to specific problems of signal analysis, namely the problems of extraction and characterisation of signals modulated in amplitude and frequency, in the asymptotic (i.e. high frequency) limit. More specificaly, we are interested not only in the asymptotic signals themselves, for which standard techniques based on Hilbert transform in general give good results, but rather on composite signals of asymptotic type, which can be defined as sums of asymptotic signals. Note that a composite signal of asymptotic type has a priori no reason to be itself asymptotic; moreover, the Hilbert transform is in general unable to give precise informations about the characteristics of the components separately, like for instance frequency and amplitude modulation laws. Nevertheless, the first tests of the algorithms presented here will be performed on asymptotic signals.

The plan of the paper is as follows: After briefly recalling in section II the basic definitions and properties of wavelet and Gabor analysis, and standard results on frequency and amplitude modulated and asymptotic signals in section III, we describe in section IV the behaviour of wavelet and Gabor analysis in the asymptotic limit, for asymptotic analysing functions. We show in section V how, for Gaussian analysing functions, the asymptotic assumption on the analysing function can be relaxed. Algorithms are proposed, which allow the extraction of instantaneous frequencies and spectral lines, and in some cases the separation of frequency-modulated components is a signal. Numerical results are presented in section VI, showing in particular that the study of the phases of wavelet and Gabor coefficients allow the construction of algorithms much more precise than those based on the energy localisation criteria. Applications to matched filtering are discussed in section VII, and section VIII is devoted to conclusions.

II: CONTINUOUS WAVELET AND GABOR ANALYSIS:

We describe in this section the basic properties of continuous wavelet and Gabor analysis, in the form we will use in this paper. For more details on such techniques, we refer to [START_REF] Daubechies | The wavelet transform, time-frequency localisation and signal analysis[END_REF], [Gr-Mo], , [Gr- and [Gr-KM-Mo].

II-1: Wavelet analysis:

We will restrict here to the wavelet analysis of the real Hardy space:

H 2 (R) = f ∈ L 2 (R); f (ω) = 0 f or ω ≤ 0
which is sufficient for our purpose. The analysis of L 2 (R) follows from simple modifications.

Start from two functions g, h ∈ H 2 (R), such that the following admissibility condition holds:

0 < |c g,h | = +∞ 0 ĝ(ω) ĥ(ω) ω dω < ∞ (II -1)
where our convention for Fourier transform is the following:

f (ω) =< f, ε ω > L 2 = +∞ -∞ f (t)e -iωt dt
with ε λ (t) = e iλt . g (resp. h) is the analysing wavelet (resp. the reconstructing wavelet).

To these functions are associated the wavelets, i.e. the functions g (b,a) and h (b,a) , defined by:

g (b,a) (t) = 1 a g t -b a (II -2)
(with the same definition for h (b,a) ), where the translation and dilation parameters run over the Poincaré half plane: H = {(b, a) ∈ R 2 ; a > 0} (notice that we have chosen here a different normalisation for the wavelets than that of [Gr-Mo]). The admissibility condition then implies the following resolution of the identity: to each s ∈ H 2 (R) one can associate the following family of coefficients (wavelet coefficients):

T s (b, a) =< s, g (b,a) > L 2 (II -3) which are simply the L 2 products of s(t) by the dilated and translated wavelets (when there is no possible confusion, the subscript s of T s will be omitted). One then has a reconstruction formula:

s = 1 c g,h +∞ 0 +∞ -∞ T s (b, a) h (b,a)
da db a (II -4) which expresses s(t) as a sum of dilated and translated wavelets h (b,a) , the coefficients of the decomposition being the corresponding wavelet coefficients (the equality being understood in the L 2 sense). In other words, the wavelets are used exactly as if they formed a biorthogonal basis, although this is far from being the case, the information provided by T s being highly redundant. It may also be convenient in practice to use a slightly modified version of wavelet analysis, which is simply to analyse and reconstruct a λ-frequency-shifted version of s(t) instead of s itself, for some real number λ. One then has a the following decomposition of s(t):

s = 1 c g,h +∞ 0 +∞ -∞ T λ s (b, a) h (b,a) ε λ da db a (II -5)
which expresses ε -λ s(t) as a sum of dilated and translated wavelets h (b,a) , the coefficients of the decomposition being the corresponding modified wavelet coefficients.

T λ s (b, a) =< s, ε λ g (b,a) > L 2 (II -6)

For simplicity, we will write T s = T 0 s . It is often convenient to work with wavelets having good localisation properties in both the direct space and the Fourier space; in general g can be chosen so that both g and ĝ have exponential decay at infinity (and/or the same for h and ĥ). This moreover allows to interpret the dilation parameter a as the inverse of a frequency shift parameter (up to a multiplicative constant, characteristic of the analysing wavelet), and then to make the connexion with Gabor coefficients.

Note that it is very useful for numerical applications to use a simplified version of the reconstruction formula:

s(t) = 1 k g +∞ 0 T λ s (t, a) ε λ (t) da a (II -7)
which expresses s(t) as a sum over frequencies, with logarithmic measure. This is equivalent to take formally h = δ in eq. (II -5), the admissibility condition reading:

0 < |k g | = +∞ 0 ĝ(ω) ω dω < ∞ (II -8)
Assuming that 0 < c g := c g,g < ∞, one then has that T λ s ∈ L 2 (H, da db a ). However, the wavelet coefficients do not span the whole L 2 (H, da db a ), but rather a subspace of it, called a reproducing kernel space:

H g = F ∈ L 2 (H, da db a ); K g,h .F = F (II -9)
with reproducing kernel given by:

K g,h (b, a; b ′ , a ′ ) = 1 c g,h < h (b ′ ,a ′ ) , g (b,a) > (II -10)
with respect to the measure da db a on the Poincaré half plane. It is worth noticing that K g,h is a projection operator, which is an orthogonal projector if and only if h = g. Notice that while the kernel depends on both g and h, the reproducing kernel space only depends on the analysing wavelet g. The reproducing kernel property can in particular be used to recover wavelet coefficients by interpolation from a discretised version of the wavelet transform.

To understand the information provided by the wavelet coefficients, a standard technique [Gr-KM-Mo] is to code them into two images, representing respectively T λ s (b, a) and Ψ λ (b, a) = Arg T λ s (b, a) . Still assuming that 0 < c g := c g,g < ∞, thanks to the Plancherel equality:

||s|| 2 = 1 c g,g T λ s (b, a) 2 da db a (II -11) T λ s (b, a)
2 is interpreted as an energy density in the Poincaré half plane, and is often refered to as scalogram. We refer to [Gr-KM-Mo] for a description of the usual conventions and interpretations of the images. We will mainly focus here on the study of the phase Ψ λ (b, a) of the wavelet transform T λ .

II-2: Gabor analysis:

Let us describe now the Gabor analysis of L 2 (R) (and not H 2 (R) as in the wavelet case, since H 2 (R) is not stable by frequency translations), in a similar way than wavelet analysis (note that Gabor [Ga] was only concerned with Gaussian analysing and reconstructing functions; we use the term Gabor analysis for simplicity). Let g, h ∈ L 2 (R), and consider the time and frequency shifted copies of them, denoted by:

g (b,ω) (t) = e iω(t-b) g(t -b) (II -12)
and the same notation for h (b,ω) . One then has the following resolution of the identity: To each s ∈ L 2 (R), one associates the corresponding Gabor coefficients:

G s (b, ω) =< s, g (b,ω) > L 2 (II -13)
from which one has the reconstruction formula, the integral converging strongly:

s = 1 2π < h, g > +∞ -∞ +∞ -∞ G(b, ω) h (b,ω) dω db (II -14)
which expresses s(t) as a decomposition over time and frequency shifted copies of a basic function. g(t) (resp. h(t)) is called the analysing (resp. reconstructing) Gabor function, and the g (b,ω) and h (b,ω) the Gabor functions. Note that contrary to the case of the dilation parameter of wavelets, the frequency shift parameter of Gabor functions run over the whole real axis. As in the previous section, the Gabor coefficients do not span the whole L 2 (R 2 , db dω), but rather some reproducing kernel subspace of it, also denoted by H g and defined by eq. (II -9), the reproducing kernel being now given by:

K g,h (b, ω; b ′ , ω ′ ) = 1 2π < h, g > < h (b ′ ,ω ′ ) , g (b,ω) > (II -15)
and still defining a projection operator, orthonormal if and only if h = g. One usually represents the Gabor coefficients by two images, namely the modulus and the argument of the coefficients. As in the wavelet case, a Plancherel equality allows the interpretation of |G(b, ω)| 2 as an energy density, usually called spectrogram, in the time-frequency plane. The argument of G(b, ω) represents a complementary information that we will study here.

III: ASYMPTOTIC SIGNAL AND EXPONENTIAL MODEL:

Throughout sections III, IV and V, we will describe estimations techniques for wavelet and Gabor coefficients in the asymptotic limit. More specifically, we will derive asymptotic developments for such coefficients. As a convention, we will use the symbol x ≈ y to specify that the expression x is approximated by the first order term y in its asymptotic development.

Let us first briefly describe here the basic definitions and properties of asymptotic signals.

III-1: The Canonical Representation of a Real Signal:

An arbitrary real signal s(t) can always be represented in terms of instantaneous modulus and argument, in the form:

s(t) = A(t) cos (φ(t)) (III -1)
Obviously, such a representation is far from being unique, in the sense that to a given real signal s(t) one can associate an infinity of pairs (A, φ) such that eq. (III -1) holds. Several authors (see e.g. [Pi-Ma]) have noticed that among these pairs, it is convenient to specify a particular one, called the canonical pair, defined as follows: consider the Hilbert transform H, defined by:

H = -i F -1 .ε.F (III -2)
where F is the Fourier transform, and ε is defined by its action on functions, according to:

(ε.f ) (ω) = Sgn(ω)f (ω) (III -3)
For instance, if α ∈ R, one obviously has that:

[H. cos] (αt) = sin(αt) (III -4)
By definition, the analytic signal Z s (t) associated with s(t) is obtained by a linear filtering of s(t) cancelling its negative frequencies, which can be expressed as:

Z s (t) = [1 + iH] (t) (III -5)
Z s is completely characterised by the pair (A s , φ s ) by setting:

Z s (t) = A s (t) exp(iφ s (t)) (III -6)
if one assumes that A s is non-negative, and that φ s takes its values in the interval [0, 2π[. (A s , φ s ) is called the canonical pair associated to s(t). Obviously,

Re(Z s ) = s (III -7)
can be written as:

s(t) = A s (t) cos(φ s (t)) (III -8)
which defines its canonical representation.

The canonical representation of s(t) allows the introduction of the instantaneous frequency ν s (t) of s(t), defined by:

ν s (t) = 1 2π dφ s dt (III -9)
Note that while the definition of the instantaneous frequency always makes sense mathematically, its physical significance can be doubtful in some particular situations, for instance when s(t) is not oscillating enough, i.e. if φ s varies slowly compared to A s , or when s(t) is the sum of two sine functions.

Several authors have studied the problem of characterizing the possible canonical pairs (A s , φ s ), in general by using spectral methods. Some sufficient conditions are derived in [Pi-Ma].

III-2: Asymptotic Analytic Signal, and Exponential Model:

Let s ∈ L 2 (R) be a real finite energy signal, of the form (III -1), where A(t) ≥ 0 and φ(t) ∈ [0, 2π[ for all t ∈ R. Assume moreover that s(t) is asymptotic, which essentially means that (see [Co], [Di]):

dφ dt ≫ 1 A dA dt (III -10)
and that the signal s(t) is locally monochromatic. An explicit form for the associated analytic signal can be obtained by using the stationary phase method (see appendix). Set:

s ± (t) = 1 2 A(t)e ±iφ(t) (III -11)
and evaluate the corresponding Fourier transforms s ± (ω). Denoting by t ± (ω) the corresponding stationary points of φ, defined by:

φ ′ (t ± ) = ± ω (III -12)
then:

s ± (ω) ≈ π 2 A(t ± )e ±iφ(t ± ) |φ ′′ (t ± )| e -iωt ± e iδ ± (III -13)
where :

δ ± = π 4 Sgn [φ ′′ (t ± )] (III -14)
Consider the analytic signals:

Z ± (ω) = 2 H(ω) s ± (ω) (III -15)
H being the Heaviside step function. Let Z(t) = Z + (t) + Z -(t), and:

Z + (t) = 1 2π Z + (ω)e iωt dω = 1 √ 2π ∞ 0 A(t + ) e i[φ(t + )+ω(t-t + )+δ + ] |φ ′′ (t + )| dω (III -16)
Using again the stationary phase method to evaluate the integral, the essential contribution to Z + (t) is given by the stationary points ω + of the argument of the integrand, such that:

∂t + ∂ω [φ ′ (t + ) -ω + ] + t -t + = 0 (III -17)
i.e. the points ω ∈ R + such that:

t + (ω) = t (III -18)
The set of such points is non-empty only if φ ′ (t) > 0. Since moreover:

∂t + ∂ω φ ′′ (t + ) = 1 (III -19)
eq. (III -16) yields:

Z + (t) = H [φ ′ (t)] A(t)e iφ(t) (III -20)
The same argument gives:

Z -(t) = H [-φ ′ (t)] A(t)e iφ(t) (III -21)
which shows that: Z(t) ≈ A(t)e iφ(t) (III -22)

In the asymptotic limit, the analytic signal of A(t) cos(φ(t)) is then close to the exponential model A(t) exp(iφ(t)). This property will be used in the following, in particular for numerical applications, where computations will be performed on real signals.

IV: THE CASE OF ASYMPTOTIC ANALYSING FUNCTION:

We now turn to the description of the wavelet (resp. Gabor) coefficients in the asymptotic limit. Assuming that the signal and the wavelets (resp. Gabor functions) are all asymptotic allows to get an approximate expression for the corresponding coefficients, by means of the evaluation of the scalar products by the stationary phase method. As we will see, this allows the specification of some particular sets of curves in the time-scale halfplane (resp. the time-frequency plane), namely the ridge and the wavelet curves (resp. the Gabor curves). The ridge has an interesting property: it describes the frequency modulation law of the signal.

One is then led to a simple geometrical description of the transform (the so-called factorisation property): essentially, the restriction of the transform to the ridge (the socalled skeleton) reproduces the analytic signal of the analyzed signal, while the restriction to a given wavelet curve (resp. Gabor curve) reproduces the corresponding wavelets (resp. Gabor functions).

We first describe here the computation of the wavelet (resp. Gabor) coefficients in the asymptotic limit, and introduce the notions of ridge and wavelet (resp. Gabor) curves; we then derive simple formulas, allowing to build algorithms for the extraction of the ridge.

IV-1: Wavelet Analysis:

Let s ∈ L 2 (R) be an asymptotic locally monochromatic signal, and let Z s (t) be its analytic signal. Set:

Z s (t) = A s (t) exp (iφ s (t)) (IV -1)
Let now g ∈ H 2 (R) be an analytic analysing wavelet, one then has:

T λ (b, a) =< s, ε λ g (b,a) >= 1 2 < Z s , ε λ g (b,a) > (IV -2)
Notice that the analysing wavelet being assumed to be progressive, it automatically removes in a smooth way the negative frequencies of s(t), and directly works with the analytic signal Z s (t), as expressed by eq. (IV -2). The smoothness of such a filtering (consequence of the smothness of ĝ) will imply a greater numerical precision than the one obtained by direct techniques involving Hilbert transform. Assume now that the analysing wavelet g(t) is itself asymptotic, and set:

g(t) = A g (t) exp (iφ g (t)) (IV -3)
then all the dilated and translated wavelets g (b,a) (t) are asymptotic too (since they are of constant shape), and T λ (b, a) takes the form of a rapidly oscillating integral:

T λ (b, a) = 1 2a +∞ -∞ A s (t) A g t -b a e i[φ s (t)-λt-φ g ( t-b a )] dt (IV -4)
the essential contribution to which is provided by the stationary points (see appendix) of the argument of the integrand, i.e. the points t s such that:

φ ′ s (t s ) = λ + 1 a φ ′ g t -b a (IV -5)
Note that any stationary point t s is a function t s (b, a) on the half plane H. Consider (b, a) ∈ H, and let t

(1) s , t (2) s ,... t (N ) s
be the associated stationary points, assumed to be in a finite number. Then the first order approximation to T λ (b, a) decomposes into the finite sum of the contributions T λ (t (n) s ) of these stationary points:

T λ (b, a) ≈ n=N n=1 T λ (t (n) s ) (IV -6)
Set for simplicity:

Φ (b,a) (t) = φ s (t) -λt -φ g t -b a (IV -7)
If there exists an integral number k ≥ 1 such that

d k Φ (b,a) dt k (t (n) s ) = 0 (IV -8) and d k+1 Φ (b,a) dt k+1 (t (n) s ) = 0 (IV -9)
then the stationary phase argument provides an explicit formula for the contribution of the stationary point t

(n) s : T λ (t (n) s ) = C k Z s (t (n) s ) g t (n) s -b a e -iλt (n) s Corr k (b, a) (IV -10)
where:

Corr k (b, a) = a d k+1 Φ (b,a) dt k+1 (t (n) s ) 1 k+1 e -i π 2(k+1) Sgn d k+1 Φ (b,a) dt k+1 (t (n) s ) (IV -11) for some constant C k . In particular, C 1 = π 2 and C 2 = 1 2 Γ 4 3 6 1 3
. It is interesting to notice that such an expression involves the evaluation of the analytic signal only at the stationary points. This property will be crucial in the following.

Using the same techniques, one can also get an explicit expression for the first corrective term in the asymptotic development of T λ (b, a), involving the first derivatives of the amplitudes of the signal and the translated and dilated wavelets. This term is easily seen to vanish in the case k=1, which shows the strongness of the approximation. The next corrective term (involving second derivatives of the amplitudes) is in general nonzero.

From now on, let us assume that to each point in the domain of study Ω ⊂ H is associated a unique stationary point t s (b, a), such that moreover Φ ′′ (b,a) (t s ) = 0 (i.e.k = 1). This assumption essentially means that the components of the considered composite signal of asymptotic type do not interact, and that one can always restrict to an Ω domain in which the wavelet coefficients of all but one component are negligible. T λ (b, a) thus takes the simple form:

T λ (b, a) ≈ π 2 e i π 4 Sgn[Φ ′′ (b,a) (t s )] Φ ′′ (b,a) (t s ) Z s (t s ) e -iλt s 1 a g t s -b a (IV -12)
On such an expression, it clearly appears that some particular sets of points will play an important role in the understanding of the wavelet transform. These sets of points are basically the points (b, a) ∈ Ω for which t s (b, a) = b (ridge of the transform) or t s (b, a) = C st (wavelet curves).

The Ridge of the Transform:

The ridge is defined to be the set of points (b, a) ∈ Ω such that t s (b, a) = b. It immediatly follows from the definition of stationary points that on the ridge:

a = a r (b) = φ ′ g (0) φ ′ s (b) -λ (IV -13)
so that the ridge is a curve R = {(b, a) ∈ Ω; a = a r (b)} ⊂ Ω in the domain Ω, from which one easily recovers the frequency modulation law. Moreover, eq. (IV -12) simplifies on the ridge to:

T λ (b, a r (b)) ≈ π 2 e i π 4 Sgn[Φ ′′ (b,a r (b)) (b)] a r (b) 2 Φ ′′ (b,a r (b)) (b) e -iλb g(0) Z s (b) (IV -14)
The restriction of the wavelet transform of s(t) to the associated ridge gives, up to a corrective function of b, the analytic signal Z s (t) of s(t). To evaluate this corrective function, note that owing to eq. (IV -13), for any (b, a) ∈ Ω:

Φ ′′ (b,a) (t) = - a ′ r (t)φ ′ g (0) a r (t) 2 - 1 a 2 φ ′′ g t -b a (IV -15)
and then:

Φ (b,a r (b)) (b) = - 1 a r (b) 2 a ′ r (b)φ ′ g (0) + φ ′′ g (0) (IV -16)
The correction function then depends on the analysing wavelet, and on the analysed signal only through the associated ridge. Once the ridge is known, one is then able to recover directly A s (b) from the wavelet coefficients T λ (b, a).

The restriction of the wavelet transform to its ridge is called the skeleton of the transform. The knowledge of the ridge and the skeleton on Ω is sufficient to characterize the transform on Ω. It is then natural to introduce the corresponding skeleton operator ), defined by:

F R 0 : L 2 (R) → L 2 (R
F R 0 .s (t) = 2 π a 2 Φ ′′ (b,a r (b)) (b) g(0) e iλb e -i π 4 Sgn[Φ ′′ (b,a r (b) (b)| T λ (b, a r (b)) (IV -17)
and whose action is to select the content of the signal localised around the ridge.

Remark: The phase shift ± π 4 is constant here, since we have assumed that for any (b, a) ∈ Ω, Φ ′′ (b,a) (t s ) = 0, so that the sign of Φ ′′ (b,a) (t s ) does not change on Ω.

The wavelet curves:

Consider (b 0 , a r (b 0 )) ∈ Ω. The wavelet curve through (b 0 , a r (b 0 )) is defined to be the connected component of (b 0 , a r (b 0 )) of the set of points (b, a) ∈ Ω such that t s (b, a) = b 0 . Let (b, a)
∈ Ω be such a point; this implies that:

φ ′ g b 0 -b a = a a r (b 0 ) φ ′ g (0) (IV -18)
The wavelet curves are then uniquely determined by the analysing wavelet. For instance, in the case of fixed-frequency wavelets, φ ′ g (0) = ω 0 , and the wavelet curve through (b 0 , a r (b 0 )) is the line a = a r (b 0 ); this property is a characteristic property of the fixed-frequency analysing wavelets.

Restricting the wavelet transform to a given wavelet curve yields:

T λ (b, a) ≈ π 2 e i π 4 Sgn[Φ ′′ (b,a) (b 0 )] a 2 Φ ′′ (b,a) (b 0 ) Z s (b 0 )e -iλb 0 g b 0 -b a (IV -19)
so that the behaviour of the wavelet transform restricted to a given wavelet curve is governed by that of the corresponding dilated and translated wavelets (the dilation parameter a being a function of the translation parameter b). Once more, T λ (b, a) is obtained from

g b 0 -b a
via a perturbation by the correction function a 2 Φ ′′ (b,a) (b 0 ) . It is worth noticing that one now has a nice geometrical picture of the wavelet transform in the time-scale half plane: when restricted to the ridge, it is essentially determined by the analytic signal of the analysed signal, and when restricted to a given wavelet curve, it is essentially determined by the wavelets, in both cases up to a corrective function, completely determined by the ridge and the analysing wavelet. This is the so-called factorisation property of the transform, that will be frequently used in the sequel.

Ridge extraction:

The next question is that of the extraction of the ridge from the T λ (b, a) coefficients. The most usual and natural information comes from the squared modulus of the coefficients, interpreted as an energy density in the time-scale half-plane H. The analysing wavelet being assumed to be maximum at t = 0, eq. (IV -12) shows that T λ (b, a) is locally maximum at b = t s , if one neglects the influence of the corrective function in the denominator of eq. (IV -14). Taking this corrective function into account shows that the maximality of T λ (b, a) at b = t s only holds in an approximate sense, which is exact if and only if the ridge a r (b) is a linear function of b, i.e. for hyperbolically frequency modulated signals, and φ ′′ g (0) = 0. The ridge of the transform can be extracted in a much more precise way from the phase of the wavelet coefficients. Let:

Ψ λ (b, a) = arg T λ (b, a) (IV -20)
Then by definition of t s (b, a):

∂Ψ λ (b, a) ∂a = - t s -b a 2 Φ (b,a) ′ t s -b a (IV -21)
and then vanishes on the intersection of the lines b = C st and the ridge. This gives an easy way to extract numerically the ridge of the wavelet transform.

Moreover, the argument of T λ (b, a), restricted to a given wavelet curve, has the following other property:

∂Ψ λ (b, a) ∂b t s (b,a)=b 0 = 1 a φ ′ g t s -b a + ∂a ∂b t s -b a 2 φ ′ g t s -b a (IV -22) = φ ′ g (0)
a on the intersection with the ridge (IV -23)

Hence the instantaneous frequency of the restriction of the wavelet transform to a fixed wavelet curve equals the central frequency of the dilated wavelet at the intersection with the ridge. This provides another algorithm for the ridge extraction, which will be discussed in section VI.

Remark: Note that in the case of a fixed-frequency analysing wavelet, such a property does not allow the extraction of the ridge. Indeed, in that case, t s (b, a) = b 0 is equivalent to a = a r (b 0 ), and one has that

∂Ψ λ ∂b (b, a) = φ ′ g (0) a
everywhere. This in particular illustrates the fact that the extraction algorithms described here are based on necessary conditions rather than necessary and sufficient conditions, and that the ridge is a subset of the set of points in Ω such that (V I -23) holds, and not the whole set of such points.

IV-2: Gabor Analysis:

We keep the same notations than in the section IV-1, and we will consider only the points (b, ω) ∈ Ω ⊂ H such that g (b,ω) ∈ H 2 (R), or at least that the Fourier transform of g (b,ω) is numerically negligible at negative frequencies. Under such an assumption, one has that:

G(b, ω) =< s, g (b,ω) >= 1 2 < Z s , g (b,ω) > (IV -25)
Assuming that all the g (b,ω) functions under consideration are asymptotic, the oscillating integral defining G(b, ω):

G(b, ω) = +∞ -∞ A s (t) A g (t -b) e i[φ s (t)-φ g (t-b)-ω(t-b)] dt (IV -26)
can be approximated by the stationary phase method: the essential contribution to eq. (IV -26) is provided by the stationary points of the phase of the integrand, i.e. by the points t s (b, ω) such that:

φ ′ s (t s ) = φ g (t s -b) + ω (IV -27)
As in the previous case, we restrict to a region Ω ⊂ H such that to each (b, ω) ∈ Ω is associated a unique stationary point t s (b, ω), which is a first order stationary point. One then has the following approximation for the Gabor transform:

G(b, ω) ≈ π 2 e i π 4 Sgn[Φ ′′ (b,ω) (t s )] Φ ′′ (b,ω) (t s ) Z s (t s )e -iω(t s -b) g(t s -b) (IV -28)
where one has set:

Φ (b,ω) (t) = φ s (t) -φ g (t -b) -ω(t -b) (IV -29)
As in the case of the wavelet transform, one can introduce the notions of ridge and skeleton of the transform, and that of Gabor curves.

The ridge of the transform:

The ridge is defined to be the set of points (b, ω) ∈ Ω such that t s (b, ω) = b. By definition of the stationary points, one has that, on the ridge:

ω = ω r (b) = φ ′ s (b) -φ ′ g (0) (IV -30)
so that the ridge gives a direct information on the instantaneous frequency modulation law of the signal. Moreover, eq. (IV -28) reduces to:

G(b, ω r (b)) ≈ π 2 e i π 4 Sgn[Φ ′′ (b,ω r (b)) (b)] Φ ′′ (b,ω r (b)) (b) g(0) Z s (b) (IV -31)
Note that:

Φ ′′ (b,ω r (b)) (b) = ω ′ r (b) -φ ′′ g (0) (IV -32)
The correction function appearing in the denominator of eq. (IV -31) is then completely determined by the analysing function and the ridge of the transform, and the Gabor analogue of the skeleton operator (see eq. (IV -20)) can still be introduced.

The restriction of the transform to its ridge is called the skeleton of the Gabor transform.

The Gabor curves:

Let (b 0 , ω r (b 0 )) ∈ Ω be a point on the ridge of the Gabor transform of s(t). The Gabor curve through (b 0 , ω r (b 0 )) ∈ Ω is the connected component of (b 0 , ω r (b 0 )) ∈ Ω of the set of points (b, ω) ∈ Ω such that t s (b, ω) = b 0 . Let (b, ω)
∈ Ω be such a point: this implies that:

φ ′ g (b 0 -b) = φ ′ g (0) + ω(b 0 -b) (IV -33)
which means that such curves are completely determined by the analysing function. The restriction of G(b, ω) to the Gabor curve associated with b 0 reads:

G(b, ω) ≈ π 2 e i π 4 Sgn[Φ ′′ (b,ω) (b 0 )] Φ ′′ (b,ω) (b 0 ) Z s (b 0 )e -iω(b-b 0 ) g(b 0 -b) (IV -34)
which shows that the behaviour of the Gabor coefficients restricted to a Gabor curve is essentially governed by that of the Gabor function e -iω(b-b 0 ) g(b 0 -b).

Ridge extraction:

The ridge can be extracted from the G(b, ω) coefficients in an approximate way by the study of the energy localisation in the time-frequency plane, i.e. by the determination of the local maxima of the modulus |G(b, ω)|. But as in the case of wavelet analysis, the study of:

Ψ(b, ω) = arg (G(b, ω)) (IV -35)
allows a much more precise determination of the ridge. Indeed, a straightforward calculation yields: ∂Ψ(b, ω) ∂ω ∂Ψ(b, ω) ∂b

t s (b,ω)=b 0 = ω + φ ′ g (b 0 -b) - ∂ω ∂b (b 0 -b) (IV -38) = ω + φ ′ g (0)
on the intersection with the ridge (IV -39)

The right hand side is nothing but the proper frequency of the Gabor function g (b,ω) (t). This gives another way of extracting the instantaneous frequency, which will be discussed in the next chapter.

V: NON-ASYMPTOTIC ANALYSING FUNCTION:

It must be noticed that the asymptotic assumption made on the analysing functions in the last section is not very natural either in the wavelet or in the Gabor philosophy, in the sense that the frequency localisation of these functions is enforced, while the time localisation is destroyed. This has moreover a bad consequence on the numerical point of view, in the sense that the computational time gets increased.

In the case where the analysing functions have a Gaussian amplitude, it is however possible to slightly generalize the stationary phase calculus to non asymptotic wavelets (or Gabor curves), and still asymptotic signals. The notions of ridge and wavelet and Gabor curves can still be introduced, and have the same intuitive meaning. The only changes concern the explicit expressions of the coefficients, and the ridge extraction algorithms.

V-1: Wavelet Analysis:

Consider an analytic analysing wavelet g ∈ H 2 (R), such that its amplitude A g (t) is a Gaussian function, which can be supposed without loss of generality to be normal. By a slight modification of the stationary phase principle (see appendix), eq. (IV -4) can now be approximated by:

T λ (b, a) ≈ A s (t s )e iΦ (b,a) (t s ) 2a +∞ -∞ A g t -b a e i (t-t s ) 2 2 Φ ′′ (b,a) (t s ) dt (V -1)
which, up to a rotation in the complex plane, reduces to the computation of a Gaussian integral. Here, t s = t s (b, a) is still defined by eq. (IV -5). A straightforward computation then yields:

T λ (b, a) = T λ (b, a) e iΨ λ (b,a) (V -2 -a) T λ (b, a) ≈ π 2 A s (t s )e -1 2 (t s -b) 2 Φ ′′ (b,a) (t s ) 2 Φ ′′ (b,a) (t s ) 2 + 1 a 4 1 + a 4 Φ ′′ (b,a) (t s ) 2 1 4 (V -2 -b) Ψ λ (b, a) ≈ φ s (t s ) -λt s -φ g t s -b a + 1 2 (t s -b) 2 Φ ′′ (b,a) (t s ) 1 + a 4 Φ ′′ (b,a) (t s ) 2 + 1 2 arctan a 2 Φ ′′ (b,a) (t s ) (V -2 -c
) Contrary to the case studied in section IV, the first corrective term in the asymptotic development, proportional to A ′ s (t s ) is not identically zero in Ω. This term reads:

T λ (1) (b, a) = A ′ s (t s ) 2a e iΦ (b,a) (t s ) +∞ -∞ (t -t s )A g t -b a e i (t-t s ) 2 2 Φ ′′ (b,a) (t s ) dt (V -3)
and can be shown to equal:

T λ (1) (b, a) = T λ (1) (b, a) e iΨ λ (1) (b,a) (V -4 -a) T λ (1) (b, a) = π 2 A ′ s (t s )e -1 2 (t s -b) 2 Φ ′′ (b,a) (t s ) 2 Φ ′′ (b,a) (t s ) 2 + 1 a 4 1 + a 4 Φ ′′ (b,a) (t s ) 2 3 4 |t s -b| (V -4 -b) Ψ λ (1) (b, a) = φ s (t s ) -λt s -φ g t s -b a + 1 2 (t s -b) 2 Φ ′′ (b,a) (t s ) 1 + a 4 Φ ′′ (b,a) (t s ) 2 + 3 2 arctan a 2 Φ ′′ (b,a) (t s ) (V -4 -c) Note that T λ
(1) (b, a) is proportional to |t s -b|, which is not the case of the second corrective term T λ

(2) (b, a) of the asymptotic development.

Ridge and skeleton:

Define the ridge of the transform to be the set of points (b, a) ∈ Ω such that t s (b, a) = b. It immediately follows from the definition of stationary points that eq. (IV -14) is still valid, so that the knowledge of the ridge of the wavelet transform is enough to give a direct access to the frequency modulation law.

The skeleton of the transform, i.e. the restriction of the wavelet transform to the associated ridge has a slightly different form than that given by eq. (IV -15):

T λ (b, a r (b)) ≈ π 2 e i 2 arctan[a r (b) 2 Φ ′′ (b,a r (b)) (b)] 1 + a 4 Φ ′′ (b,a) (t s ) 2 1 4 e -iλb-iφ g (0) Z s (b) |t s -b| (V -5) while from eq. (V -4 -b): T λ (1) (b, a r (b)) = 0 (V -6)
Finally, eq. (IV -16) still holds. Thus, the skeleton of the wavelet transform of the asymptotic signal s ∈ L 2 (R) gives, up to a corrective function which is completely determined by the analysing wavelet and the ridge itself,the analytic signal Z s (t) of s(t). Eq. (V -6) also indicates that the first corrective term to this approximation vanishes identically.

Wavelet curves and ridge extraction:

The extraction of the ridge of the transform from the zeroes of ∂Ψ λ ∂a b=b 0 is no longer valid. Indeed, it is not difficult to check that ∂Ψ λ ∂a b=b 0 involves a (in general nonzero)

term proportional to ∂t s ∂a b=b 0 . Consider now the wavelet curves, as defined in section IV-1, i.e. essentially as the sets of points (b, a) ∈ Ω with constant stationary point. The restriction of the wavelet transform to such a curve is essentially governed by the corresponding translated and dilated wavelets (its expression is obtained by setting t s = b 0 in eq. (V -4), for some fixed b 0 . Note that the phase of T λ (b, a) is now perturbed by the contribution of the modulus of the wavelet.

Nevertheless, it can be a useful quantity for the extraction of the ridge, since on a fixed wavelet curve:

dΨ λ db t s (b,a)=b 0 = 1 a φ ′ g (0) + 1 2 φ ′′′ g (0) 1 + φ g (0)a ′ r (b) + φ ′′ g (0) 2 (V -7)
on the intersection with the ridge, so that for any analysing wavelet such that its central frequency is an inflection point (i.e. φ ′′′ g (0) = 0), eq. (IV -21) is still valid, providing an algorithm for the extraction of the ridge.

V-2: Gabor Analysis:

Assume that the window is a function of Gaussian type, once more assumed to be normal. The approximation:

G(b, ω) ≈ 1 2 A s (t s )e iΦ (b,ω) (t s ) +∞ -∞ e -1 2 (t-b) 2 e i(t-t s ) 2 Φ ′′ (b,ω) (t s ) dt (V -8) leads to: G(b, ω) = |G(b, ω)| e iΨ(b,ω) (V -9 -a) |G(b, ω)| ≈ π 2 A s (t s )e -1 2 (t s -b) 2 Φ ( b,ω) ′′ (t s ) 2 1+Φ ′′ (b,ω) (t s ) 2 1 + Φ ′′ (b,ω) (t s ) 2 1 4 (V -9 -b) Ψ(b, ω) ≈ Φ (b,ω) (t s ) + 1 2 (t s -b) 2 Φ ′′ (b,ω) (t s ) 1 + Φ ′′ (b,ω) (t s ) 2 + 1 2 arctan Φ ′′ (b,ω) (t s ) 2 (V -9 -c)
Moreover, one easily checks that, as in the wavelet case, the first corrective term in the asymptotic development is proportional to |t s -b|.

Ridge and skeleton:

Still defining the ridge of the Gabor transform by t s (b, ω) = b, eq. (IV -30) remains valid. The skeleton of the transform is now given by:

G(b, ω r (b)) ≈ π 2 e i 2 arctan Φ ′′ (b,ω r (b)) (b) 1 + Φ ′′ (b,ω r (b)) (b) 2 1 4 g(0) Z s (b) (V -10)
and reproduces the analytic signal Z s (t), up to the corrective function:

π 2 e i 2 arctan Φ ′′ (b,ω r (b)) (b) 1 + Φ ′′ (b,ω r (b)) (b) 2 1 4 g(0) (V -11)
which is completely determined by the ridge.

Gabor curves and ridge extraction:

Defining the Gabor curves by eq. (IV -36), it is not difficult to see that on a given Gabor curve, eq. (IV -39) becomes:

dΨ db t s (b,ω)=b 0 = ω + φ ′ g (0) + 1 2 φ ′′′ g (0) 1 + ω ′ r (b) + φ ′′ g (0) 2 (V -12)
on the intersection with the ridge, leading to a simple algorithm for the extraction of the ridge in cases where φ ′′′ g (0) = 0.

VI: RIDGE AND SKELETON EXTRACTIONS:

We now turn to numerical applications of the asymptotic formulas described in the previous sections. As stressed in the introduction, we are interested in the study of asymptotic signals, but also of composite signals of asymptotic type, defined as linear combinations of asymptotic signals. Clearly, such signals need not be asymptotic themselves, because of interference phenomena. We will consider in the same footing wavelet and Gabor analysis, which have been shown to have the same behaviour. In the wavelet case, we will only describe the λ = 0 case for simplicity, keeping the possibility of introducing the λ parameter to enforce the asymptotic character of the signal.

Throughout this section, we will only use for numerical computations the so-called Morlet wavelet:

g(t) = e -t 2 2 e i ω 0 t (V I -1)
Although the Morlet wavelet is not rigorously admissible, ω 0 can be chosen sufficiently large so that ĝ(0) is numerically negligible. Under such an assumption, g(t) can be considered as a numerically progressive wavelet, so that eqs. (IV -2) and (IV -25) hold. Notice that in that case, the wavelet curves are nothing but the a = C st curves, and the Gabor curves are the ω = C st curves, which makes the numerical implementation simpler.

VI-1: An algorithm for ridge extraction: non-interacting ridges:

Let us first consider the case of non-interacting components, that is the case where the asymptotic components of the signal live in different regions of the time-frequency half-plane. Because of the good localisation properties of both g and ĝ, we are then led to the case of asymptotic signals. We will say that in such a case, the individual ridges are non-interacting. Our aim is to extract numerically an individual ridge from the (wavelet or Gabor) transform, but also to avoid computing the whole transform by restraining the computation to a neighborhood of the ridge. Let us first give a simple example:

An intuitive view of the problem:

Let us look at the simplest case of a monochromatic signal: let

s(t) = exp(iω s t) (V I -2)
Its wavelet transform with respect to the progressive analysing wavelet g reads:

T (b, a) = |T (b, a)| e iΨ(b,a) = ḡ(aω s ) exp(iω s b) (V I -3)
Clearly, the ridge of the wavelet transform of such a signal is given by:

a r = ω 0 ω s (V I -4)
Moreover, a r can obviously be reached iteratively as follows: let a 0 be a primary value of the dilation parameter. Compute:

ω 1 = ∂Ψ(b, a 0 ) ∂b (V I -5)
and set:

a 1 = ω 0 ω 1 (V I -6)
Clearly a 1 = a r , and:

ω 0 ∂Ψ(b,a r ) ∂b = a r (V I -7)
The ridge is then obtained as a fixed point of an iterative scheme in such a particular case.

The algorithm:

Let us now generalise the approach described above, focusing on the algorithmic aspect. Consider the case of discrete sequences: let T e = 1 ν e be the sampling period, consider the wavelet transform T a (k) = T (kT e , a) of s k = s(kT e ) for a fixed value of the dilation parameter, and Ψ a (k) its argument. Denote by D b an estimator of a discrete differentiation with respect to b. Then, the following algorithm allows the construction of a trajectory a r (kT e ) in the Ω domain such that:

D b .Ψ a (k) = ω 0 a (V I -8)
The algorithm can be described in a similar way than above: let a 0 (t 0 ) be an initial estimate of a r (t 0 ). Then the estimates for a r (t 0 + kT e ) are given by:

a i+1 (t 0 + kT e ) = ω 0 D b .Ψ a i (t 0 + kT e ) (V I -9) a 0 (t 0 + (k + 1)T e ) = a c (t 0 + kT e ) (V I -10)
where a c is the "converged solution", i.e. the value a c = a j of the dilation parameter for which one has, for some fixed positive value of the required precision ǫ (arbitrarily small fixed positive number):

a j+1 -a j a j < ǫ (V I -11)
In practice, such an algorithm converges very rapidly, two or three iterations being necessary for each signal sample. This in particular means that one does not have to compute the whole wavelet transform to get an estimate of its ridge.

Numerical results:

Since the signals studied here are assumed to have non-interacting ridges, we first give the example of an academic signal s(t) = exp(iω s t + iα s ln(1 + βt)), hyperbolically modulated in frequency (it is clearly asymptotic since its amplitude does not vary). Numerical tests can also be performed on real signals, like for instance a sonar signal emitted by a bat. In Fig. VI-2-a and VI-2-b, its frequency modulation law is seen to be close to an hyperbolic one. Notice the existence in Fig. VI-2-a of an auxiliary ridge at higher frequencies, which very presumably results from an echo phenomenon during the experiment. It is quite interesting to compare the results obtained with wavelets to those obtained with standard Hilbert transform-based techniques, in the case where the previous signal is embedded in noise (here a uniform computer noise, with amplitude comparable to that of the signal itself). In that case, the standard method (Fig. VI-3-a) fails to extract correctly the frequency and amplitude modulation laws, because they use the whole timefrequency content of the input signal, while the good localisation properties of both g and ĝ allow to take into account only the subset of the time-frequency plane where the asymptotic signal is important, and thus yield much more precise results (Fig. VI-3-b), in particular for the instantaneous frequency. Moreover, wavelet and Gabor transform techniques also allow the separation of the auxiliary component of the signal, which is completely impossible with the standard techniques. Notice also that the Hilbert transform method is able to yield a precise result in that case only because the second component of the signal is much smaller in amplitude than the first one.

VI-2: Interacting ridges:

Consider now the cases where the input signal s(t) is the sum of say two asymptotic signals, and that the corresponding wavelet (or Gabor) transforms are non-negligible in intersecting domains Ω 1 and Ω 2 . This is what happens for instance in Fig. VI-4, which represents the wavelet transform of two linearly modulated cosine functions, with the same modulation rate α. Interference phenomena cause here a strong perturbation of the phase of the wavelet transform, so that the ridge extraction algorithm provided by eq. (IV -21) is no longer exact. Indeed, in Fig. VI-4, small "bubbles" regularly appear on the ridge plot (we will understand the origin of such bubbles a little later). In some particular cases, the interferences can be made negligible by a suitable choice of the analysing wavelet. For instance, with the same signal than the one studied in fig. VI-4, taking a Morlet wavelet with a much greater value of ω 0 will yield a family of wavelets with a much better localisation in the Fourier space, so that they will be able to separate the two components. This is the case in Fig. VI-5, where the proper frequency of the Morlet wavelet has been chosen large enough to avoid interferences between the two components, and allows their separation and their extraction. Unfortunately, such a simple trick does not always work, in particular in the case of complex signals. For instance, Fig. VI-6-a and VI-6-b represent the phase of the Gabor transform of a trumpet sound and a female voice respectively, and the associated ridges. The ridge structure is in such cases so complex that it cannot be analysed by the methods described here, even with a particular choice of the analysing functions. This is a strong limitation of the asymptotic approximations-based algorithms described in this paper.

In some very simple and particular example, the existence of the bubbles can be identified as coming from interactions between two or more modulated component. Consider then a signal made of two spectral lines with equal amplitudes (set to 1) and pulsations ω 1 and ω 2 : s(t) = exp(iω 1 t) + exp(iω 2 t), to be analysed with Gabor analysis, with analysing function g(t) = exp(-t 2 2σ ) exp(iωt) (fixed frequency analysing function). A straightforward computation shows that b-derivative of the phase of G(b, ω) is given by:

∂Ψ(b, ω) ∂b = ω * + δ tanh(s) 1 + tan 2 (bδ) 1 + tan 2 (bδ) tanh 2 (s) (V I -12)
where:

ω * = ω 2 + ω 1 2 (V I -13) δ = ω 2 -ω 1 2 (V I -14) and s = σ 2 (ω -ω * )δ (V I -15)
The ridge equation is then a fixed point equation on s:

s = σ 2 δ 2 tanh(s) 1 + tan 2 (bδ) 1 + tan 2 (bδ) tanh 2 (s) (V I -16)
Depending on the values of b and of the parameters, the number of solutions of such a fixed point equation determines the existence of bubbles. For large values of σδ, i.e. for non-interacting ridges, the fixed points are close to ω 1 and ω 2 . Otherwise, consider the case of small values of σδ; then whenever tan 2 (bδ) is large enough, the fixed point equation has three solutions (s = 0 and two symmetric ones), while for small values of tan 2 (bδ), only s = 0 is possible, so that a bubble may be formed. This is illustrated in Fig. VI-7.

For real signals, the bubble structure is in general more complex, as illustrated in Fig.

VI-8-a and VI-8-b, showing the ridges of the Gabor transform of a clarinet sound (at different times), containing three frequency modulated components. In that example, the interaction between the different ridges increases with time, and bubbles appear.

VI-3: Spectral lines estimation:

We now restrict our interest to special kinds of signals, that is to asymptotic signals that are sums of spectral lines:

s(t) = sl i (t) (V I -17)
where a spectral line is of the form:

sl i (t) = A i (t) exp [iω i t + iδ i ] (V I -18)
Such signals are quite interesting for instance from the acoustical point of view, since they allow to describe many vibrating systems. We will discuss here a method of estimation of the frequencies ω i and the amplitudes A i (t), under the assumptions that all the A i (t) are analytic and vary slowly in the domain of study.

The case of one spectral line:

Consider first the case of a monochromatic signal of constant amplitude A and pulsation ω s . Then clearly T (b, a) = A exp(iω s b)ĝ(aω s ). Assuming that ĝ is real, which is indeed the case for the Morlet wavelet, clearly ∂Ψ(b,a) ∂b = ω s . Assuming now that A s (t) is an analytic function of t:

sl(t) = A s (t) exp(iω s t) (V I -19)
and since the analysing wavelet is assumed to be maximal at t = 0, consider the following asymptotic development of T sl (b, a):

T sl (b, a) = A s (b)e iω s b ḡ(aω s ) + ∞ k=1 (-ia) k k! d k ḡ dω k (aω s ) d k A s (b) db k (V I -20)
Then, since ĝ′ (ω 0 ) = 0, setting:

a s = ω 0 ω s (V I -21)
one has that:

T sl (b, a s ) ≈ A s (b)e iω s b (V I -22)
up to a second order error.

Remark: Notice that such a simple calculus is not the consequence of the stationary phase estimations described in sections IV and V, which in such a particular case fail to yield approximate expressions for the wavelet transform (since we use the Morlet wavelet, the stationary points of the argument of the integrand are infinite order stationary points). Nevertheless, the philosophy of the technique is the same.

One can then use the same technique than the one used for ridge extraction to determine the value of ω s : compute the wavelet transform of sl(t) for a fixed value a 0 of the dilation parameter, and the derivative of its phase with respect to the translation parameter b: the resulting frequency ω 1 will yield a corresponding value a 1 = ω 0 ω 1 ; then compute T (b, a 1 ) and iterate the algorithm. The fixed point of the iteration will yield an estimate for ω s , the error being of second order in the derivatives of A s .

The case of two spectral lines:

Assume now that s(t) is the sum of two spectral lines, and that the corresponding amplitudes are regular and slowly varying enough so that the preceding approximations are still valid. Clearly, if the two frequencies are sufficiently far away from each other, the fast decay of ĝ will allow to treat them independently, and to separate them. If not, one still can choose a sufficiently large ω 0 to enforce the resolution of the wavelets in the Fourier space; but as stressed in section VI-2, such a trick is not so convenient for numerical computations, since it involves convolutions of the signal by wavelets with a large support. Moreover, increasing the central frequency of the analysing wavelet increase in the same time the values of the derivatives of ĝ appearing in eq. (V I -20), so that the asymptotic assumptions are no longer valid. It is actually possible to refine the method, by cancelling the interferences between the wavelet transforms of the two components. Consider for simplicity the case where the amplitudes of the two components are constant:

s(t) = A 1 cos(ω 1 t) + A 2 cos(ω 2 t + δ 2 ) (V I -23) Then: T (b, a) = A 1 ḡ(aω 1 )e iω 1 b 1 + C(a)e i(ω 2 -ω 1 )b+iδ 2 (V I -24)
where

C(a) = A 2 ḡ(aω 2 ) A 1 ḡ(aω 1 ) (V I -25)
Let ε < 2π ω 2 -ω 1 be a non-negative real number, and let:

T = n 2π ω 2 -ω 1 + ε (V I -26)
for some positive integral number n. Introduce:

Ω a = 1 T b 0 +T b 0 ∂Ψ(b, a) ∂b db (V I -27) A straightforward computation shows that if C(a) < 1, Ω a = ω 1 + O( ε n ), and if C(a) > 1, then Ω a = ω 2 + O( ε n ).
Obviously, the value of C(a) depends only on the value of the dilation parameter one started with. Starting from a 0 , and setting a 1 = ω 0 Ω a , a simple analysis (see e.g. [START_REF] Seip | Some remarks on a method for detecting spectral lines in a signal[END_REF]) shows that Ω a 1 is a fixed point of Ω. One then has an iterative method to estimate the frequencies of the spectral lines. Notice that in that case, the error is proportional to ε n , and that the signal being assumed asymptotic, n can be made large enough to allow to neglect the error. A more detailed analysis of the convergence properties of the algorithm in the case of two spectral lines can be found in [Gu].

The general algorithm:

The spectral lines detection algorithm is the straightforward generalisation of the above discussion. Let s(t) be the input signal, containing a certain number of spectral lines ω 1 ,...ω N . Let a 0 be a starting dilation parameter, and T a (large) real number. Set:

a i+1 = ω 0 Ω a i (V I -28)
Then the {a i } sequence converges numerically very rapidly to a fixed point a ∞ , such that:

a ∞ ≈ ω 0 ω k (V I -29)
for some k = 1, ...N . Moreover, different starting values of the dilation parameter a 0 will yield different spectral lines. It must also be remarked that once a spectral line has been detected, it can be extracted by means of the factorisation property: indeed, consider for simplicity the case of two spectral lines. Then:

1 ḡ(ω 0 ) T (b, ω 0 ω 1 ) ≈ A 1 (b)e iω 1 b + K(b)e iω 2 b (V I -30)
where:

K(b) = A 2 (b) ḡ( ω 0 ω 2 ω 1 ) ḡ(ω 0 ) < A 2 (b) (V I -31)
so that the relative amplitude of the second component has decreased. One can iterate such a procedure by considering the left hand side of eq. (V I -30) as a new input signal, and then neglect the second component after a certain number of iterations. In such a case, the first component is completely characterised, and can eventually be removed from the s(t) signal. Examples of such spectral lines substractions will be shown in the next section.

Separation of the contribution of each spectral line:

In the case of the two frequencies signal, we saw that it is possible to exactly determine the values of the frequencies. Unfortunately, the modulus of the restriction of the wavelet transform to the values of the dilation parameter corresponding to those frequencies are oscillating functions. Nevertheless, for a fixed number of spectral lines of frequencies ω n , n = 1, ..N , it is possible to 'disentangle' the components and then determine the amplitude modulation law of each component independently, without changing the time and frequency supports of the wavelets. For this purpose, one may assume that N and the ω n have previously been determined (this is possible to do thanks to the spectral lines estimation algorithm and the factorisation property described before). Linear transformations on the restrictions of the wavelet transform to these values of the dilation parameter lead to a linear system, whose solution yields the amplitude modulation law of all spectral lines (or more precisely the first order term in their asymptotic development).

Numerical results:

The first application of wavelet based spectral lines estimations was an application to NMR spectroscopy, where a spectral line corresponds to an excitation degree of freedom of a molecule. One then has to determine exactly the value of the frequency each spectral lines, and to separate the spectral lines of the studied molecules from those of the surrounding

F R τ .s (t) = 2 π a r (t) 2 Φ ′′ (t,a r (t)) (t) g(0) e -i π 4 Sgn[Φ ′′ (t,a r (t)) (t)] T s(t + τ, a r (t)) (V II -2)
which is the τ -shifted analogue of the skeleton operator introduced in eq. (IV -17). At τ = 0, F R τ associates with s its wavelet transform restricted to the R ridge, and corrected by a (local) normalisation function and phase shift, needed to compensate the terms introduced by the stationary phase evaluation of the integrals. For τ = 0, the wavelet transform is now restricted to a τ -shifted copy of the R ridge:

R τ = {(b, a) ∈ Ω; a = a r (b -τ )} (V II -3)
Clearly, when s(t) = s(t -τ 0 ), R τ coincides with the ridge of s, so that both functions:

SC(τ ) =< F R τ .s, Z s > (V II -4) SE(τ ) = ||F R τ .s|| 2 (V II -5)
are maximal at τ = τ 0 (to see this, just take for F R τ .s its expression provided by the stationary phase argument in section IV, and evaluate the scalar products). SC is called the skeleton correlation function, and SE is the skeleton energy function. They are to be compared with the autocorrelation function C, output of the matched receiver:

C(τ ) = Z s (t)Z s (t -τ )dt (V II -6)
Consider now s(t) of the form:

s(t) = N i=1 A i s(t -τ i ) + r(t) (V II -7)
for some real numbers A i and τ i , and some deterministic or random r(t). s(t) contains N shifted copies of s(t), and the corresponding ridge is a multicomponent ridge, which contains in particular N copies of the ridge of s(t). Assume moreover that the τ i parameters are large enough to ensure that these components of the ridge are non-interacting, in the sense defined in section VI. Then SC(τ ) is still locally maximal (more precisely it exhibits sharp peaks) at τ = τ i , i = 1, ..N , at least as soon as the wavelet coefficients T r (b + τ, a r (b)) are small compared with all the A i T s (b, a r (b)). This in a first step allows the determination of the time delay parameters τ i , and in a second step the determination of the A i parameters, by using the factorisation property described in section VI. In the opposite case, additional local maxima appear; one says that the false alarm probability becomes more important. Note that the local maxima of SC(τ ) correspond to the cases where the ridges of s(t) coincide with the theoretical ridge of s(t). This explains the name given to the algorithm. θ i . One then gets modified wavelet coefficients, which no longer satisfy the reproducing kernel equation. This in particular means that one does not have any control a priori on the signal reconstructed from these modified coefficients. Nevertheless, a ridge extraction gives a vertical ridge like in eq. (V II -9). Another modification is then necessary.

In a third step, to achieve coherence between all the voices, the output of each voice is multiplied by a corrective function:

Corr k (t) = 2 π a 2 k Φ ′′ (t-θ k ,a k ) (t s ) g(0) (V II -13)
and a phase shift:

δφ k = -2πφ s (t s (t -θ k , a k )) + Cst (V II -14)
which are completely determined by the a-priori known ridge. Finally, the modified outputs of all voices are summed up in a coherent way:

W (t) = N i=1 Corr k (t)e -δφ k T s (t -θ k , a k ) (V II -15)
or in a partially coherent way:

E(t) = N i=1 |Corr k (t)T s (t -θ k , a k )| 2 (V II -16)
The generic structure of the algorithm is summarized in figs. VII-3-a and VII-3-b.. It is then not difficult to show that if the receiver is matched to the frequency modulation law of the incoming signal, i.e. if the incoming signal is hyperbolically frequency modulated, both W (t) and E(t) are maximal at t = T . If the input signal is now of the form given by eq. (V II -7), the output signal will exhibit sharp peaks at all the t = T +τ i , yielding precise estimates for the τ i 's.

Implementation and numerical results:

The main point in the practical implementation of the algorithm is the choice of the analysing wavelet, and the determination of the number of different dilation parameters to be taken into account. For a fixed analysing wavelet, the number of voices is determined by the constant relative bandwidth constraint, and by the specification of the constant overlap of the corresponding dilated wavelets in the Fourier space: the Fourier transforms of the dilated wavelets intersect when they take a fixed value, set to ǫ times the maximal value ĝ(ω 0 ), for some fixed ǫ < 1.

Let us consider the simple example of the hyperbolically frequency modulated signal, with instantaneous frequency:

ν(t) = ν s 1 + αt t ∈ [0, T ] (V II -17)
and set ν e = ν(T ). The values of the considered dilation parameters a k are then completely determined (see [E-To] for more details):

a k = ν 0 ν s 1 + 1 2Q ǫ k+1 1 -1 2Q ǫ k-1 (V II -18)
Here Q ǫ = ω 0 -2 ln(ǫ) is a surtension factor. Figs. VII-4-a and VII-4-b represent respectively E(t) and |W (t)| 2 in the case N = 12 and ω 0 = 15, exhibiting sharp maxima at t = T . In the presence of a Gaussian noise (Figs. VII-5-a and VII-5-b), the signal to noise ratio at the output W (t) enables an accurate time estimation, and may be considered as an unbiased time estimator. Various attempts [E-To] indicate that the improvement of signal to noise ratio due to the coherent receiver is near BT , B being the bandwidth of the received signal.

Remark: ridge rectification and bat's sonar receiver:

Such results may appear as common ones for signal processing. Nevertheless, the model described here is in a strong connexion with biological data about bats. Bats use an ultrasonic sonar system to detect obstacles during cruise, to locate and pursuit their preys [Si]. Among various species of bats, Vespertillionids emit wideband signals, the amplitude and the frequency of which are modulated in such a way that they can be considered as asymptotic (see fig. VI-2-a). The acoustic receiver of bats is built up with nervous cells, the response of which is similar to that of tuned circuits. The relative bandwidth of these cells is constant, like that of wavelets. These cells may then be represented as a "wavelet receiver", and wavelet representation is then well suited to model such a receiver. Recent results by J. Simmons [Si] have moreover pointed out two main facts, which enforce the connexion with wavelet representation:

-Auditory representations of echoes is displayed by nervous cells as a spectrogram (or scalogram) due to the Corti's cells acting as constant relative bandwidth filters.

-The outputs of such tuned cells are the processed parts of the incoming echoe relative to the relative bandwidth of the filters. These outputs (nervous response) are called neurograms, and appear to be deeply related to scalograms.

VIII: CONCLUSIONS:

We have reviewed here the main properties of wavelet and Gabor representations when applied to asymptotic signals or composite signals of asymptotic type. The main ingredient of the applications we described is that time-frequency representations allow to separate components of signals located in different regions of the time-frequency plane, and then to synthetise them independently.

It is then natural to see that the techniques described above give remarkable results in the case of asymptotic signals, and in the case of composite signals of asymptotic type when the components do not interact in the time-frequency plane; nevertheless, it is still difficult to isolate components that strongly interact in the time-frequency plane, as for instance in the case of the bubbles described in section VI. In such cases, new techniques are needed, and there is some hope that wavelet packets [Co-Me-Qu-Wi], [START_REF] Torrésani | Time-Frequency Representations: wavelet packets and optimal decomposition[END_REF] or even more general decompositions and adaptive analysis schemes could yield better results.

Another point that has to be stressed is that in the analysis described here, which is essentially concerned with a careful study of the phase of the time-frequency representation, there is absolutely no difference between wavelet and Gabor analysis (of course the explicit expressions and the algorithms are slightly different, but the general procedure is completely identical). The only difference that appears between the two approaches lies in the behaviour of the squared modulus of the representation, i.e. the time-frequency energy density. This is the consequence of a well known property of time-frequency representations, that is that a given representation is optimal (in the sense of energy localisation) for only one type of frequency modulation law.

Finally, notice that the computations of sections IV and V can trivially be generalised to the two-dimensional wavelet transform with rotations, described un [Mu]. The only difference is that the notion of wavelet curve has to be replaced by that of wavelet surface.

The second term of the asymptotic development, i.e. the term proportional to the first derivative of M (b,a) is easy to estimate too:

T (1) (b, a) = M ′ (b,a) (t s )e i Φ (b,a) (t s ) 2a ∞ -∞
(t -t s )e i 1 2 (t-t s )2 Φ ′′ (b,a) (t s ) dt = 0 (A -9) which proves the strongness of the approximation.

A-2: s(t) asymptotic and |g(t)| Gaussian :

In the particular case of an analysing wavelet with a Gaussian amplitude: (A -12 -c) It must be noticed here that the next term in the asymptotic development, proportional to the derivative of A s (t) is not identically zero in that case. Its computation still reduces to a Gaussian integration, and yields:

g(t) = e -

  Gabor analogue of eq. (IV -18). The ridge of the Gabor transform is then easy to determine numerically. |G(b, ω)| is obtained from |g(b 0 -b)| by a perturbation by Φ ′′ (b,ω) (b 0 ) , and arg G(b, ω) has the following property, when restricted to a fixed Gabor curve:

  Fig. VI-1-a represents the phase and the modulus of T s , and the ridge a = a r (b) is represented on the phase plot. Fig. VI-1-b the frequency and amplitude modulation laws of such a signal, computed via wavelet transform. The results are in excellent agreement with the input data.

  medium. This is what is exhibited in figs. VI-9. Fig. VI-9-a represents the Fourier spectrum of the original NMR signal, where the most important spectral line corresponds to the sorrounding medium, and where one is interested in the extraction and separation of the other ones, which are negligible in amplitude. Fig. VI-9-b represents the Fourier spectrum of the same signal, after the removal of the most important spectral line. Spectral lines estimations are also useful for noise filtering. Fig. VI-10-a represents the Fourier spectrum of a sum of spectral lines imbedded in noise, and fig. VI-10-b shows the sum of the extracted spectral lines. Let us finally show an example where the algorithms are used for modelisation in the aim of resynthesis. Figs. VI-11 show the extraction of 3 harmonics of the 16 main harmonics of the syllabe ′ la ′ . The resynthesis is obtained by summing up all the harmonics [Gu]. transform:

  t 2 2 e iφ g(t) (A -10) the asymptotic assumption on g(t) can be dropped as follows: instead of considering the Taylor series development of M (b,a) , consider the development of A s (t) only. This also means that one assumes the variations of A s to be slow compared with those of the Gaussian amplitude of the wavelets. The first term in the asymptotic development of T (b, a) then reads:T (0) (b, a) = 1 2a A s (t s )e iΦ (b,a) (t s ) -t s ) 2 Φ ′′ (b,a) (t s ) dt (A -11)which then once more reduces to a Gaussian integration, and yields:T (0) (b, a) = T (0) (b, a) e iΨ (0) (b,a) (A -12 -(b, a) = φ s (t s ) -φ g t s --b) 2 Φ ′′ (b,a) (t s ) 1 + a 4 Φ ′′(b,a) (t s ) 2 + 1 2 arctan a 2 Φ ′′ (b,a) (t s )

T

  FIGURE CAPTIONS:

Fig

  Fig. VII-2-a: Same than Fig. VII-1-a, with an additional uniform noise with amplitude ten times larger than the maximal amplitude of the asymptotic signal. Fig. VII-2-b: Same than Fig. VII-1-b, with an additional uniform noise with amplitude ten times larger than the maximal amplitude of the asymptotic signal. Fig. VII-2-c: Same than Fig. VII-1-c, with an additional uniform noise with amplitude ten times larger than the maximal amplitude of the asymptotic signal. Fig. VII-3-a: Definition of the time delays θ k . Fig. VII-3-b: Generic structure of the ridge rectification algorithm. Fig. VII-4-a: E(t) function in the case of an hyperbolically frequency modulated signal, and a receiver with N = 12. Fig. VII-4-b: |W (t)| 2 function in the case of an hyperbolically frequency modulated signal, and a receiver with N = 12. Fig. VII-5-a: Same than Fig VII-4-a, with an additional centered Gaussian noise, with input signal to noise ratio equal to 1/5. Fig. VII-5-b: Same than Fig VII-4-b, with an additional centered Gaussian noise, with input signal to noise ratio equal to 1/5.

arctan a 2 Φ ′′ (b,a) (t s ) (A -13 -c) Throughout this paper, the symbol ≈ will mean that one has replaced an expression by the first order term of its asymptotic development.

VII: APPLICATIONS TO MATCHED FILTERING:

As stressed in section II, the computation of the wavelet (resp. Gabor) coefficients T (b, a) (resp. G(b, ω)) of s(t) can be thought of as band-pass filtering of the signal, which localize it into cells (in the time-frequency plane) of constant relative bandwidth ∆ν ν (resp. constant absolute bandwidth ∆ν). It is then natural to combine such filters to build by hand other filters, adapted to special kinds of signals. This is the purpose of matched filtering, the matched filter being defined as the linear filter which optimizes the signal to noise ratio of the output (see e.g. [Pap], [Bl]). Notice that in such an approach, it is also natural to use the simplified reconstruction formula (II -7) to recover the signal as the sum over all the filtered components.

We consider here the case of asymptotic locally monochromatic signals (for which approximate expressions for the wavelet and Gabor transforms are available), and we restrict to the special case of wavelet analysis; analogous algorithms can be derived in the framework of Gabor analysis.

We will describe two algorithms, called the ridge coincidence algorithm and the ridge rectification algorithm respectively, whose purpose is detection and time-location of a-priori well known asymptotic locally monochromatic components s(t) in a signal, imbedded in an input signal s(t). These two algorithms model filters close to the matched filter, in the sense that the output signal exhibits a sharp peak at the estimated date, like the autocorrelation function of s(t).

Both algorithms basically follow the same philosophy. They first use the wavelets to localise the signal into time-frequency cells of constant relative bandwidth, and the contributions of such cells are summed in different ways. Finally, the same operations are performed with time-shifted copies of the cells. The output signal is in both cases different than the autocorrelation function of s(t), but basically has the same characteristics.

Let us stress that one can build in the same way many slightly different versions of such algorithms; our goal in this section is not to present algorithms in a final form, but rather to illustrate how wavelet analysis can be used to get approximations of the matched filter. Let us also point out the simple geometrical interpretation of these algorithms in the time-frequency half-plane.

VII-1: The Ridge Coincidence Algorithm:

Let s ∈ L 2 (R) be an asymptotic locally monochromatic signal of the form given by eq. (III -8), and let g ∈ H 2 (R) be an analytic analysing wavelet, written in the form given by eq. (IV -3). Denote by:

the ridge of the wavelet transform T (b, a) of s(t). It is convenient to introduce here the τ -shifted skeleton operator

and T s is its wavelet

Numerical results:

As an example, consider an hyperbolically frequency modulated signal, starting at τ 0 = 0. Figs. VII-1-a to VII-1-c represent respectively the C(τ ), SE(τ ) and SC(τ ) functions for such a signal. Notice that SC(τ ) is sharply peaked at τ = 0. In the case of a noisy signal, up to some critical value of the signal to noise ratio, time estimation via the skeleton correlation is still possible, as exhibited in Figs. VII-2-a to VII-2-c (here the asymptotic signal is the same than previously, and the noise is a uniform random signal, whose amplitude equals 10 times the maximal amplitude s(t)). The performances are then comparable to those of the matched filter.

VII-2: The Ridge Rectification Algorithm:

We now describe another algorithm dedicated to matched filtering. The main idea of this algorithm is to transform a signal with a given ridge into another signal, with a ridge as close as possible to a vertical ridge:

for some fixed b 0 . Then, use the simplified reconstruction formula to construct the output of the algorithm from the modified wavelet coefficients. We describe here the generic structure of the algorithm in the particular case of hyperbolic frequency modulation laws, for which it is known that representations associated with the affine group are optimal in terms of energy localisation in the time-scale plane (see e.g. [Be-Be]).

The ridge rectification algorithm goes as follows:

In a first step, pick a fixed number of dilation parameters a i , i = 1, ..N , chosen in such a way that the corresponding relative bandwidth

is constant. To these dilation parameters are then associated the proper frequencies of the corresponding dilated wavelets:

and then the time delays t i such that:

. Fix a reference time T , and set:

In a second step, compute the wavelet transform of the input signal for these fixed values of the dilation parameter, and time-translate the i-voice by a time-delay equal to

APPENDIX: THE STATIONARY PHASE ARGUMENT:

Let us describe briefly here the approximation techniques we use in this paper. We derive in a simple way the first terms of the asymptotic developments of T (b, a) we start from in sections IV and V (the approximations for G(b, ω) and T λ (b, a) go exactly in the same way).

From eq. (IV -4), one has the following expression for the wavelet transform of s ∈ L 2 (R):

where:

From now on, we will assume that both g(t) and s(t) are analytic functions, so that their series development make sense.

A-1: g(t) and s(t) asymptotic: Assume that both s(t) and g(t) are asymptotic signals in the sense defined in section III. Then T (b, a) is defined by a rapidly oscillating integral, so that asymptotic techniques can be used (see e.g. [Co]). Moreover, because of the regularity and the fast decay of g(t), the integration bounds do not contribute to the asymptotic development of T (b, a). Let t s = t s (b, a) be such that: Φ ′ (b,a) (t s ) = 0 (A -4)

Then, by general arguments, the essential contribution to (A -1) is provided by the neighborhood of such points. Moreover, assuming for simplicity that t s is the unique point such that (A -4) holds and that M (b,a) is non negligible, then the first term in the asymptotic development of T (b, a) reads:

A straightforward Gaussian integration thus yields: (A -7)

where we have assumed that: Φ ′′ (b,a) (t s ) = 0 (A -8)