Matthieu Quantin

Benjamin Hervy

Florent Laroche

Jean-Louis Kerouanton

MASS CUSTOMIZATION FOR CULTURAL HERITAGE 3D MODELS

Keywords: Industrial heritage, knowledge management, 3D modeling, information visualization, levels of detail

Digital technologies in the last twenty years have offered cultural heritage (CH) new possibilities in conservation and promotion. 3D digitization has especially become more and more affordable and efficient. This leads to massive digitization projects and increasing amount of CH digital data. As an engineering team working on industrial techniques for reverse engineering, we are deeply affected by this effect. In this paper we propose a way to combine semantic information on top of the acquisition and modeling steps in order to manage heterogeneous historical data. We illustrate our approach with a use case composed of three overlapping historical objects related to Nantes' harbour history.

I. Introduction

Considering the amount of overlapping 3D digitized models of cultural heritage (CH) objects, connecting their references would enrich every one of them and optimize historical research process. Focused on technical and industrial heritage, this paper describes a process from mass digitization to mass customization. Indeed, modeling the context make digitized artifact become unique and customized rather than a 3D mesh among others. The context also provides a way to explore the system the technical objects belong to. Three variables are considered: temporal, spatial and thematic. Each one includes several levels of detail:

• The temporal dimension can be related to a short period of time (machinery production cycle), an intermediate period (machinery lifecycle in the factory) or even a long period (the whole lifecycle including the functioning period and its lifecycle as an heritage object).

• The spatial dimension can be related to both the geometric shape and the geographical area. Although a CH object can have no specific position (e.g. vehicles, steam engine or any generic technical object), it is linked to one or more geographical environments.

• The thematic dimension is not finite and can be related to many research domains: social, economics, mechanics, aesthetic, etc. This aspect is significant to understand different points of view on the object.

The main issues lie in the management of historical knowledge and the access to information. Many projects such as [START_REF] Karmacharya | ArchaeoKM: Managing Archaeological data through Archaeological Knowledge[END_REF] deal with the management of archeological data and try to provide semantic access to geographical information. In our case, we deal with narrower scales and focus on the analysis of an object and its context. We mostly work on the object as part of a technical system. Our approach assumes that historical knowledge should not depend on physical representation (such as a 3D model) quality [START_REF] Laroche | Advanced Industrial Archaeology: A new reverse-engineering process for contextualizing and digitizing ancient technical objects[END_REF]. Nevertheless, it can enhance the information system with the object geometrical features. This paper discusses the idea of making connexions between the physical object and a knowledge database. This connection allows identification of lacks of information, and cross navigation through the three variables mentioned above. In addition, semantic keyword navigation is possible.

We study Nantes' (french city) harbour history through three different historical objects (see Fig. 1). This use case leads us to those considerations.

Based on this use case of conservation and promotion of cultural and industrial heritage, we wanted to make connexions between the different pieces of information, gathering the knowledge capitalized regarding the parts of the use case. Theses connexions make possible to manage artifacts' physical representations (3D, dynamic simulations) based on available knowledge (and not the contrary). Combining both historical sources and produced knowledge (under the form of 3D models or written analysis) would allow researchers and broad public to explore different scales and dimensions of the objects involved. Because we believe 3D representations should create at least as many questions as it answers, our proposal should be a way to identify historical research hypothesis.

II. Methodology: Multi-dimension knowledge modeling process

3D digitization has especially become more and more affordable and efficient. This leads to massive digitization projects and increasing amount of cultural heritage (CH) digital data [START_REF] Addison | Emerging trends in virtual heritage[END_REF]. Every project in CH conservation and promotion we are involved into implies the creation of new historical knowledge under the digital form. Many of them have implicit links: geography, time, thematics. Therefore, we want to manage CH objects uniqueness by creating connexions in order to access information through different dimensions. This is what we call "mass customization for cultural heritage objects".

We state here that every historical object can be represented by a combination of the three axis of Fig. 2 according to common practices in industrial archaeology [START_REF] Martin | Industrial Archaeology[END_REF], [START_REF] Trinder | Industrial archaeology: a discipline?[END_REF].

Therefore, we propose a methodology based on the creation of relationships between descriptive pages dedicated to cultural objects and these 3 axis. First, we create a descriptive page. This descriptive page is a classical historian's production composed of full text document with historical sources. Thus, this work is subjective so we aim to confront it to other ones. Based on the content, we link this page with other ones through the identification of common keywords. This linking step can be semi-automatic thanks to text mining techniques such as named entity recognition. Then, we attach three variables for the positioning of the descriptive page in the Space-Time-Thematic cube. The thematic dimension can be deduced from the previous step of keywords identification. The spatial positioning can be divided into three parts :

1) Mandatory association between the descriptive page and geographic area(s). Whenever this kind of information is not avalaible, the information system has to highlight this page as incomplete. This would allow users to identify potential research axis. 2) Potential association between a 3D part and other spatial information (map, 3D part, etc.) 3) Potential association between descriptive page and 3D part (when available)

This methodology implies that every descriptive page can be defined on Fig. 2 by a group of blocks. Considering our three use cases, they gather a group of descriptive pages (blocks). On Fig. 2, we represent 3 volumes corresponding to an average of all these blocks for each use case.

Our approach needs to be adaptative to different historical data (physic and semantic) acquisition processes. The objective is to model information in order to combine the three dimensions represented on Fig. 2. 3D physical information is not attached to any thematic but to space and time. Our methodology allows to switch from mass digitization to mass customization of virtual artifacts. Virtual artifacts can be both physical or semantic/abstract representation about studied objects.

Once the user accesses a descriptive page, the system provides links according to neighbouring blocks on Fig. 2 and keywords cooccurrences. This way, the user is able to navigate all around his first search. He can explore combined, dynamic 3D environments and semantic networks. The cube is not an user interface, but define what we call the "neighbourhood" of a piece of information.

III. Use case presentation: heterogeneous data

In addition to the 3 dimensions mentioned in section II, we also have to deal with heterogeneous historical materials 1 . Three industrial archaeology studies compose our use case and illustrate this double issue. Yet, our proposition aims to be applied to any other cultural heritage study of this kind. a) "Nantes1900" project: a museum project focusing on industrial heritage through the case of a 100 years old scale model. This scale model is part of Nantes (France) history museum. In this project, we have a delimited spatial area which is the harbour of a city. We also have many heterogeneous historical documents related to hundreds of points of interest. The 3D model of the physical scale model comes from scanner digitization (Creaform HandyScan). It is composed of 100 million points. In addition, we gathered a corpus of 1000 historical documents (postal cards, archives documents, paints, photographies, books, etc.) and hundreds of notes written by historians. Based on this corpus, we created a knowledge database to handle these information. The global level of detail is related to the city scale approximately between 1899 and 1920. The database stores information about many thematics related to industrial activity: social ("labour conditions in the early 1900's", people biographies), economics (companies' chronicle), technology (shipyards, steam engines), architecture and urbanism (urban network), etc. b) "Halles Alstom": an industrial complex being restorated and based in the harbour of Nantes. This historical "object" is part of the scale model mentioned above. In this case, we have a detailed historical analysis that lead to a monography (master thesis). We also have specific historical documents such as 3D models coming from photogrammetry (using Agisoft Photoscan). This analysis focuses on an area related to a factory scale, but with a broader time period . Compared to the Nantes1900 project, only some thematics are covered: technology (mechanics, materials), architecture and urbanism. c) Nantes transporter bridge: a bridge present on the harbour area from 1903 to 1958. Based on historical and technical archives, we created a functional CAD 3D model. This 3D model can simulate the transporter bridge during its functioning cycle. On this object, the different dimensions taken into account are smaller: the geometry is related to a single industrial element, and the thematic dimensions focus on its technical and aesthetic aspect.

IV. Concrete application

A. Results of semantic integration for 3D CH models

In order to illustrate our methodology described in section II, we present here seven descriptive pages (among the hundreds we have for our three use cases) and the way they connect in the system. Transporter Bridge: 2 descriptive pages: the windowed cabin and the category of cable-stayed bridges.

Fig. 4 represents an image of the whole proposed system simplified with the example of the 7 previously mentioned descriptive pages. Historical studies focus on specific subject leading to a particular research context. Such context is marked out with temporal, spatial and thematic boundaries. When the user selects a descriptive page, the system provides a set of connexions with other pages based on keywords cooccurences. Resulting network has spatial, temporal, thematic and semantic features according to the cube Fig. 2. In addition, some of the network's elements can be connected to 3D models. This model is a way to navigate between different levels of detail and dimensions.

In our example above, one possible keyword query returns 3 descriptive pages ("1943 bombing on Halles Alstom","Halles Alstom" and "Nantes transporter bridge"). Then, the user can choose one of them (e.g. "Nantes transporter bridge"). Finally, the system provides all the related content (text analysis, historical sources and 3D models), and the context of the user selection. This context includes semantic linked pages ("windowed cabin", "cable-stayed bridges") and a spatiotemporal-thematic neighbourood. Neighbourhood criteria are defined by the user (e.g. "5 to 10 years" and/or "less than 500 meters").

B. Proposal of virtual reality based interface

We propose several solutions to explore this heterogeneous data. First, the user can navigate in the virtual scene: it allows to visualize semantic and spatio-temporal information through different points of view (place, date, thematics).

The user can select several descriptive pages as a group. Then, he can choose to highlight the entire group to quickly see the locations of items of interest.

Each descriptive page (Figure 5) can show the evolution of the item selected through time. The user can choose a period of time within the available data and check the details relevant to this age. It can be a different picture, 3D model or a different description for example. If we combine all these data types, the interface and the visualization can become messy and non-efficient. The solution would be to use filters and layers to show only what the historian is interested in at the moment. For example, our prototype proposes a cursor to adjust the vizualisation area around the user, allowing to manage the number of semantic details.

V. Conclusion

This paper proposes a way to manage structured historical information combining 3D models, spatio-temporal data and semantic information. Descriptive pages as defined in this paper are the core of this framework because they fit with historian working methods. The use cases we mentioned highlight the usefullness of our system. They show there is a significant interest in navigating through different dimensions (see Fig. 2), and levels of detail, regardless of a 3D model representation. The system allows the historians to create and visualise connexions, based on existing knowledge elements, in order to identify research issues among other things. 3D models finally come as additive tool.

Our approach presents some limits due to the complexity of historical information and 3D models post-processing (which currently is mostly manual). Next step is to improve the linking process so that it would be transparent for historians. We need to find a compromise between automatic and manual processing because of the necessary validation of historical information. Indeed we have to avoid a blackbox system.

Semantically enriched 3D cultural heritage models [START_REF] Maass | Seamless integration of labels into interactive virtual 3D environments using parameterized hulls[END_REF] could be a way to improve visualisation of information within the 3D model. It would allow end users to navigate from the 3D model to the semantic network and vice versa. Another improvement lies in the possibility for users to specialize information models both for semantic (specific ontology) and spatio-temporal data (for example CityGML urban models [START_REF] Costamagna | CityGML for Architectural Heritage[END_REF]).

Contextualizing historical objects is a common issue. Dealing with tied historical use cases, our approach aims at providing a common framework to combine them. This makes crossings possible between knowledge elements related to these historical objects. Thus, each object helps with contextualizing the others.

Fig. 1 :

 1 Fig. 1: Geographical positioning of three studied objects between 1899 and 1913: the global area represents a part of Nantes harbour (C), and the two other objects (highlighted): Transporter bridge (A) and "Halles Alstom" (B).

Fig. 2 :

 2 Fig. 2: Schematic representation of the Space-Time-Thematic cube. Three historical objects are represented in different colors to illustrate possible cross referencing: Nantes1900 in red, Halles Alstom in green and Transporter bridge in blue.

Fig. 3 :

 3 Fig. 3: Use cases' 3D models with different acquisition techniques : 3D scan, photogrammetry, and 3D reconstruction (CAD model).

Fig. 4 :

 4 Fig. 4: Schematic representation of the relationships between 3D models, semantic information and spatio-temporal data.

Fig. 5 :

 5 Fig. 5: Visualizing connexions between items in the virtual scene based on database information.

Heterogeneity concerns sources types, analysis nature, levels of detail, object nature, temporalities of the object


G Beylkin

B Torrésani

IMPLEMENTATION OF OPERATORS VIA FILTER BANKS; Hardy Wavelets and Autocorrelation Shell

We consider implementation of operators via filter banks in the framework of the Multiresolution Analysis. Our method is particularly efficient for convolution operators. Although our method of applying operators to functions may be used with any wavelet basis with a sufficient number of vanishing moments, we distinguish two particular settings, namely, orthogonal bases, and the autocorrelation shell.

We apply our method to evaluate the Hilbert transform of signals and derive a fast algorithm capable of achieving any given accuracy. We consider the case where the wavelet is the autocorrelation function of another wavelet associated with an orthonormal basis and where our method provides a fast algorithm for the computation of the modulus and the phase of signals. Moreover, the resulting wavelet may be viewed as being (approximately, but with any given accuracy) in the Hardy space H 2 (IR).

I INTRODUCTION

In this paper we introduce a method for design of digital filters and consider their implementation and application via filter banks. The design of digital filters is always a trade-off between accuracy and efficiency. For a number of operators this trade-off obtained via traditional filter design techniques is not adequate, especially if high precision is required. As examples, consider the Hilbert transform or operators of fractional differentiation where an accurate traditional implementation over a wide band necessarily implies a long filter.

Signal processing is not the only field where fast and accurate implementation of such operators is of interest. In numerical analysis Fast Multipole Method (FMM) [START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF][START_REF] Greengard | A fast algorithm for particle simulations[END_REF][START_REF] Carrier | A fast adaptive multipole algorithm for particle simulations[END_REF] has been developed to address this problem. Although this method proved to be efficient in numerical analysis, it did not so far find its way into the signal processing. A possible explanation for this may be the traditional reliance on filtering operations in the signal processing community. In fact, multiresolution techniques evolved in signal processing as subband coding techniques [START_REF] Esteban | Application of quadrature mirror filters to split band voice coding systems[END_REF][START_REF] Smith | Exact reconstruction techniques for tree-structured subband coders[END_REF]. The original motivation for subband coding was essentially optimal representation and compression of signals. The introduction of the orthonormal bases of wavelets [START_REF] Stromberg | A Modified Franklin System and Higher-Order Spline Systems on R n as Unconditional Bases for Hardy Spaces[END_REF][START_REF] Meyer | Ondelettes et Opérateurs[END_REF] and the concept of multiresolution analysis [START_REF] Mallat | Multiresolution approximation and wavelets[END_REF][START_REF] Meyer | Ondelettes et fonctions splines[END_REF] have led to the development of a broader concept of harmonic analysis of signals where subband coding became a natural way of representing and analyzing signals. These notions also migrated to numerical analysis, where they were applied to the problem of efficient representation and application of operators [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF]. In particular, it was shown in [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] how to use wavelet bases to "almost diagonalize" certain classes of operators, for example, pseudodifferential and Calderón-Zygmund operators. For signal processing applications the approach of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] is also of interest. Although designed for numerical purposes, algorithms of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] explicitly use quadrature mirror filters (QMFs) with the "exact reproduction property" and may be viewed as a link between signal processing techniques and traditional numerical computing.

The method developed in this paper is different from that of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] and may be also easily implemented both in software and hardware. It avoids the construction of non-standard forms which, although quite efficient, are not as simple as the "filter banks approach" described in this paper. Briefly, we decompose a signal into different scales (subbands) and implement operators as subband filters. Let us, for example, consider a convolution operator T . The wavelet ψ may be written as ψ(ξ) = m 1 (ξ/2) φ(ξ/2), where ψ and φ are the Fourier transform of the wavelet and of the scaling function. The 2πperiodic square-integrable function m 1 represents one of the Quadrature Mirror Filters (QMFs).

Our approach is based on an observation that if the wavelet ψ(x) is sufficiently well localized in the Fourier domain, one may write (T ψ)(ξ) ≈ m T (ξ/4) φ(ξ/4), (

where m T is a 2π-periodic function which is computed given the symbol of the operator T . The accuracy of the approximation in (1.1) is controlled by the number of vanishing moments of the wavelet ψ (it might be necessary to consider (1.1) on each scale separately if the symbol of T is not homogeneous). As a result, the operator T is implemented using filters m T (may be different on different scales), where m T plays a role similar to that of the filter m 1 of the QMF pair. The major difference, which is already visible in equation (1.1), is that the m T filter performs a scaling by a factor 4 instead of 2. This has some practical implications, which are discussed throughout the paper. The factor 4 may be replaced by a factor 2 n , n ≥ 2, as a way to improve accuracy. In this way, the procedure for design of these subband filters allows us to attain any desired accuracy. The approach of this paper (as that of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF]) may be traced back to Calderón-Zygmund and Littlewood-Paley approach to harmonic analysis of functions and operators which we apply here to design filters given the symbol of an operator. Our method may be used with any wavelets (associated with quadrature mirror filters) which possess sufficient number of vanishing moments. In cases where the associated scaling function also has vanishing moments (which implies that the corresponding coefficients are well approximated by samples on fine scales), our algorithm leads to a fast method for computing the Hilbert transform (and, thus, modulus and phase) of signals. This is the case for the autocorrelation wavelets derived in [START_REF] Saito | Multiresolution representations using the auto-correlation functions of compactly supported wavelets[END_REF] which we consider here in some detail. Although our approach is quite general, we concentrate on several specific examples, such as the Hilbert transform, operators of differentiation and, more generally, convolution operators. In particular, we consider the Hilbert transform both as an example and an important special case, because of its particular status (it is one of the simplest and most popular examples of Calderón-Zygmund operators) and its relevance in signal processing. In our approach the Hilbert transform (as well as a number of other operators) is completely expressed in terms of filter banks, which makes it easy to handle for a wide variety of scientific communities. In addition to the Hilbert transform, we construct derivative and integration operators including those of fractional order.

Our approach also allows us to consider the following related problem. In signal processing it is often useful to deal with wavelets that belong to the complex Hardy space H 2 (IR), i.e., wavelets such that their Fourier transform is zero for negative frequencies. For instance, such wavelets are considered to be more efficient for the identification of "chirps" (i.e. amplitude and frequency modulated components) in signals. In particular, it is easy to identify the "carrier" frequency and remove it by shifting in the frequency domain if necessary. However, there does not exist any orthonormal multiresolution analysis of H 2 (IR) where the associated wavelet has such a property [START_REF] Auscher | Il n'existe pas de base d'ondelettes régulières dans l'espace de Hardy H 2 (IR)[END_REF]. Nevertheless, as we show in this paper, it is possible to keep the algorithmic structure of multiresolution analysis and use wavelets that approximate the Hilbert transform of a given real-valued function with any given (but finite) accuracy. The sum of the original wavelet and i times its approximate Hilbert transform yields a new wavelet that is approximately in H 2 (IR). As a direct consequence, we obtain a fast algorithm for the computation of the Hilbert transform, and a pyramidal algorithm for discrete wavelet transform with complex analytic (or progressive) wavelet [START_REF] Grossmann | Reading and understanding continuous wavelet transforms[END_REF]. This may be thought of as a starting point to carry on an analysis similar to that developed in [6, ?, 12].

The first part of the paper is devoted to the representation of operators in terms of filter banks. To illustrate our approach, we derive in Section II the function m T in (1.1) for the Hilbert transform. We then present in Section III a general approach of filter banks implementation of convolution operators (in these two sections, we consider compactly supported orthogonal wavelets and obtain O(N ) algorithms). We then turn to the particular case of the autocorrelation of Daubechies' compactly supported wavelets in Section IV. We again consider first the Hilbert transform and then, in Section V, develop approximations of other operators (e.g., operators of fractional differentiation and integration) by our technique.

In the second part of the paper we address the signal processing problems using autocorrelation wavelets. In Section VI we consider computing the Hilbert transform as well as modulus and phase of signals and, in Section VII, address the problem of decomposition of band-pass signals into amplitude and frequency modulated components. The main tools are the same as in the first part of th paper, namely filter banks implementation of Hilbert transform. However, the algorithm we describe is O(N log N ) since we choose to use a translation-invariant version of wavelet transform. Finally, Section VIII is devoted to conclusions.

II THE HILBERT TRANSFORM

As a way of introduction, let us consider our approach for implementing the Hilbert transform,

(Hf )(x) = 1 π p.v. ∞ -∞ f (s)
x -s ds .

(2.1)

The general case and some other operators will be considered in the following sections.

Let us start with the usual MRA (see Appendix). It is well known that the Hilbert transform of the wavelet ψ(x) is still a wavelet. Since we will be intererested in computing coefficients Hf, ψ j k , let us consider the function H * ψ given in the Fourier domain by

(H * ψ)(ξ) = i sgn(ξ) ψ(ξ) = i sgn(ξ)m 1 ξ 2 φ ξ 2 , (2.2)
where H * denotes the adjoint of H. Let us consider the 4π-periodic function

m 2 (ξ) = i k sgn(ξ + 4πk)m 1 (ξ + 4πk)χ [-2π,2π] (ξ + 4πk). (2.3)
Our main observation is as follows: although i sgn(ξ)m 1 (ξ) is not a periodic function, the product i sgn(ξ)m 1 (ξ) φ(ξ) may be well approximated by m 2 (ξ) φ(ξ), due to the fast decay of φ(ξ).

Proposition II.1

1. The Fourier coefficients b ℓ of the function

m 2 (ξ) = 1 √ 2 ℓ b ℓ e iℓξ/2
are given by

b ℓ =      - 1 π k g k k -ℓ 2 for odd ℓ 0 for even ℓ (2.4)
2. The Fourier coefficients b ℓ have the following asymptotics

b 2ℓ-1 ∼ O((2ℓ -1) -M -1 ) (2.5)
Notice that asymptotics (2.5) coincides with the asymptotics expected from the general approach of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF]. Proof:

1. Consider the Fourier coefficients of m 2 (ξ) b ℓ = i √ 2 4π 2π -2π m 1 (ξ)sgn(ξ)e -iℓξ/2 dξ (2.6) = i 4π 2π 0 m 1 (ξ)e -iℓξ/2 -m 1 (-ξ)e iℓξ/2 dξ (2.7) = - 1 2π k g k 2π 0 sin(k - ℓ 2 )ξ dξ (2.8)
which yields (2.4).

2. The decay of the coefficients b 2ℓ-1 is governed by the regularity of m 2 (ξ) at the origin, i.e. by the number of vanishing moments of the wavelet. The asymptotics (2.5) is obtained using Taylor series expansion of (2.4) and taking into account the vanishing moments of the sequence {g ℓ }:

b 2ℓ-1 = - 1 π k g k k -2ℓ-1 2 (2.9) = 2 π k g k 2ℓ -1 ∞ p=0 2k 2ℓ -1 p (2.10) = 2 π k g k 2ℓ -1 ∞ p=M 2k 2ℓ -1 p (2.11) ∼ O((2ℓ -1) -M -1 ) (2.12)
which proves (2.5).

In other words, the Hilbert transform is essentially a local operator on functions which have sufficient number of vanishing moments.

Remarks:

1. We notice that the sequence {b 2ℓ-1 } may be viewed as the Hilbert transform of the sequence {g k }. We note that for sequences the singular behaviour of the Hilbert transform at the origin is avoided by the 1 2 term in the denominator, and the slow decay at infinity is replaced by (2.5).

2. A statement similar to Proposition II.1 may be proved for wavelets with rational m 0 and m 1 . Let

m 1 (ξ) = P (e iξ ) Q(e iξ ) (2.13)
be a rational 2π-periodic function, where both P and Q are polynomials. Set p(ξ) = P (e iξ ), q(ξ) = Q(e iξ ) and consider

p 2 (ξ) = i k sgn(ξ + 4πk)p(ξ + 4πk)χ [-2π,2π] (ξ + 4πk), (2.14)
and the 4π-periodic function

m 2 (ξ) = p 2 (ξ) q(ξ) (2.15)
Since p(ξ) carries all the vanishing moments of m 1 (ξ), the quality of the approximation of i sgn(ξ)m 1 (ξ) φ(ξ) by m 2 (ξ) φ(ξ) is controlled by that of the approximation of i sgn(ξ)p(ξ) φ(ξ) by p 2 (ξ) φ(ξ). Such approximation was discussed in Proposition II.1.

Using Proposition II.1, we obtain from (2.2) the following approximation

(H * ψ)(ξ) ≈ m H * (ξ/4) φ(ξ/4), (2.16)
where m H * is a 2π-periodic function,

m H * (ξ) = m 2 (2ξ)m 0 (ξ).
(2.17)

The coefficients dj k = Hf, ψ j k may then be computed as

Hf, ψ j k ≈ 2 j/2 f (ξ)e ik2 j ξ m 2 (2 j-1 ξ) φ (2 j-1 ξ)dξ (2.18) = ℓ b 2ℓ-1 f (ξ)e i(2k-ℓ-1/2)2 j-1 ξ φ (2 j-1 ξ)dξ (2.19) = ℓ b 2ℓ-1 s j-1 2k-ℓ+1/2 , (2.20)
since the coefficients b 2ℓ+1 are real.

Notice that as a consequence of equation (2.20), we now need "half-integer samples" of the coefficients s j-1 . In view of equation (2.16), this may be interpreted as follows. The shift by 1/2 actually amounts to a switch from V j to V j-1 , as sugested by the following lemma:

Lemma II.1 Let f (x) ∈ V 0 or f (x) ∈ W 0 , and set f 1 2 (x) = f x + 1 2 . Then f 1 2 (x) ∈ V -1 .
The proof of Lemma II.1 is very simple. Observe that f (x) ∈ V 0 is equivalent to f (ξ) = q(ξ) φ(ξ), where q(ξ) is a 2π-periodic square-integrable function, and that translation by 1/2 is equivalent to multiplication by exp{iξ/2} in the Fourier space, so that f 1 2 (ξ) = exp{iξ/2}m 0 (ξ/2) φ(ξ/2) which, in turn, implies that f 1 2 (x) ∈ V -1 . The argument for W 0 is the same.

From Lemma II.1 follows a simple algorithm to compute (2.20). In order to obtain "half-integer samples" on all scales except the finest scale j = 0, it is simply sufficient to avoid subsampling at the first step of the algorithm 1 .

At scale j = 0, we need to use interpolation to obtain "half-integer samples" s 0 k+1/2 , using the assumption that f ∈ V 0 . We have

s 0 k+1/2 = f, φ 0 k+ 1 2 = ℓ s 0 ℓ φ 0 ℓ , φ 0 k+ 1 2 , (2.21)
where coefficients φ 0 ℓ , φ 0 k+ 1 2 are easily obtained using autocorrelation function of scaling function described in Appendix.

Summarizing, we obtain the following O(N ) scheme for computing the coefficients sj k of the Hilbert transform of a function f (x) on subspace V j , j ≥ 0, assuming that projection of f (x) onto an approximation space, say V 0 is known.

1. Compute the coefficients s j k , j ≥ 1 and k ∈ Z Z/2 using equation (9.9) in Appendix and the Hilbert differences dj k = Hf, ψ j k via the pyramidal algorithm given in equation (2.20) on all scales.

2. Use the usual filters for the reconstruction on all scales: replace equation (9.10) by the following

sj k = ℓ h 2ℓ-k sj+1 ℓ + g 2ℓ-k dj+1 ℓ (2.22)
1 Notice that in the particular case of the Hilbert transform, we would then use the even (non-subsampled) scaling function coefficients to compute the difference coefficients d j k , and the odd ones for the Hilbert transform coefficients dj k .

We note that the computational cost is a factor of two compared with the usual wavelet transform.

Remark: Recently Auscher [START_REF] Auscher | Il n'existe pas de base d'ondelettes régulières dans l'espace de Hardy H 2 (IR)[END_REF] and, independently, Lemarié [START_REF] Lemarié-Rieusset | Sur l'existence des analyses multirésolution en thórie des ondelettes[END_REF] have shown (as a part of a more general result) that given a Multiresolution Analysis (MRA), it is possible to associate another MRA with the Hilbert transform of the associated wavelet ψ(x). We notice that our construction is different in the sense that we never need to consider the scaling function associated with the new MRA. Also, we always derive approximate formulas since our goal is to develop efficient approximations suitable for numerical implementations.

III IMPLEMENTATION OF OPERATORS VIA FILTER BANKS

Let us now turn to the more general case of filter bank implementation of linear operators. We show that the approach we developed for the Hilbert transform may be generalized to convolution operators with non-oscillatory kernels. We also analyze the connection with the BCR approach [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] and express the filter bank implementation as an approximation (with controlled accuracy) to the non-standard and standard forms (NS-form and S-form) approaches of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF], where we take into account only a few blocks of corresponding representations.

III.1 Convolution operators

Let us consider a more general convolution operator with symbol a(ξ),

ĝ(ξ) = a(ξ) f (ξ), (3.1)
and compute coefficients dj k of the projection of g on W j ,

dj k = g(x)ψ j k (x)dx = 2 j/2 2π ĝ(ξ) ψ(2 j ξ) e ik2 j ξ dξ. (3.2)
We look for the coefficients bj

l such that dj k ≈ ν bj ν s j-1 ν+2k , (3.3)
where

s j-1 ν = 2 (j-1)/2 2π f (ξ) φ(2 j-1 ξ) e iν2 j-1 ξ dξ, (3.4)
and index ν is not necessarily an integer. Typically, we will set ν to be half-integer and, if more precision is needed, we will demonstrate that ν can be taken to be in 2 -N Z Z. Using (3.1), we write (3.2) as

dj k = 2 j/2 2π a(ξ) f (ξ) m 1 (2 j-1 ξ) φ(2 j-1 ξ) e ik2 j ξ dξ, (3.5)
and note that it is sufficient to find an approximation

a(ξ) m 1 (2 j-1 ξ) φ(2 j-1 ξ) e ik2 j ξ ≈ 2 -1/2 ν bj ν-2k φ(2 j-1 ξ) e iν2 j-1 ξ , (3.6) or a(ξ) m 1 (2 j-1 ξ) ≈ 2 -1/2 ν bj ν e iν2 j-1 ξ . (3.7)
over the essential support of the function φ(2 j-1 ξ). Let us replace 2 j-1 ξ by ξ in (3.7) and approximate a(2 -j+1 ξ) m 1 (ξ) by the following 4π-periodic function (i.e., restrict ν to half integers,

ν = n/2) m j # (ξ) = k a(2 -j+1 (ξ + 4πk)) m 1 (ξ + 4πk) χ [-2π,2π] (ξ + 4πk). (3.8)
Due to M -1 vanishing derivatives of m 1 (ξ) at points 2πk, k ∈ Z Z, at these points the "break" in the function due to the periodization occurs only in the higher derivatives. Since derivatives of m 1 (ξ) vanish at ξ = 0, we can consider symbols that have a singularity at ξ = 0, e.g., the Hilbert transform and fractional derivatives.

The Fourier coefficients of m j # may be found by computing

b j n = √ 2 4π 2π -2π a(2 -j+1 ξ) m 1 (ξ) e -inξ/2 dξ , (3.9)
since they are related to the bj ν coefficients by a complex conjugation and a redefinition of the index ν, and we have

m j # (ξ) = 1 √ 2 n b j n e inξ/2 . (3.10)
The coefficients b j n have a fast asymptotic decay. Let us consider two cases, first where the symbol a(ξ) has at least M continuous derivatives (M is the number of vanishing moments of the basis) and, second, where a(ξ) has a singularity at ξ = 0 but has at least M continuous derivatives elsewhere. In the first case we simply integrate (3.9) by parts M -1 times using 2i n

m d m dξ m e -inξ/2 = e -inξ/2 , (3.11)
and notice that the boundary terms vanish so that

b j n = O(n -M +1 ). (3.12)
In the second case we split the integral into two over [-2π, 0] and [0, 2π], and then integrate by parts to obtain again (3.12). Summarizing the results of this section, we show that the action of a convolution operator T with symbol a(ξ) on a function f (x) may be obtained as follows:

T f (x) = j k dj k ψ j k (x) = j k n b j n s j-1 2k-n/2 ψ j k (x) (3.13)
where the coefficients b j n are given in (3.9). This implies that in order to evaluate T in the wavelet basis, we compute

dj k = n b j n s j-1 2k-n/2 (3.14)
Again, we notice that "half-integer samples" of the coefficients s j-1 are needed, and refer to the discussion in Section II.

The filter B j = {b j n } defined in (3.9) may be used similar to the filter G in QMF pair. Namely, a signal f is decomposed using the filter pair H and B j and then reconstructed with the usual QMF pair H and G to yield the desired result. Notice that the filters B j depend on the scale. If the symbol a(ξ) is homogeneous of degree m, then b j n = 2 -jm b 0 n . Remark: It is easy to see the similarities with the decomposition into a biorthogonal basis. We note, however, that there is a single MRA in our approach.

III.2 Time-dependent symbols

A number of interesting questions arises if we consider a more general class of symbols of pseudodifferential operators,

T σ f (x) = 1 2π I R σ(ξ, x) f (ξ)e iξx dξ, (3.15)
where symbol σ(ξ, x) ∈ S ′ (IR 2 ). The operator T σ may be expressed as an integral operator of the form

T σ f (x) = I R K(x, y)f (y)dy . (3.16)
Here K(x, y) ∈ S ′ (IR 2 ) is the distribution kernel of T σ , given by

K(x, y) = F -1 1 σ (x -y, x) = L(x -y, x) , (3.17)
where F 1 denotes Fourier transform with respect to the first variable. To develop our approach, we need to specify further the symbol class we are working with. We restrict to the class of the so-called Calderón-Zygmund kernels, i.e. kernels K(x, y) such that

∂ α x ∂ β y K(x, y) ≤ C α,β |x -y| 1+α+β . Let f (x) ∈ L 2 (IR), and let us compute the projection of T σ f onto W j , T σ f, ψ j k = 2 -j/2 I R×I R L(x -y, x)f (y)ψ 2 -j x -k dxdy (3.18)
Let us focus on the integral with respect to x first. We write

L(x -y, x)f (y)ψ 2 -j x -k dx = L x -y, k2 j ψ 2 -j x -k dx + R(y; j, k) (3.19)
where R(y; j, k) is some remainder. It follows from general arguments involving the vanishing moments of

ψ(x) that |R(y; j, k)| = O 2 M (j-1/2) .
From now on, we assume that the remainder may be neglected, i.e., that we are at a sufficiently fine scale. Assuming that we may change the order of summation in (3.18), we arrive at an approximation

T σ f, ψ j k ≈ 1 2π σ(ξ, k2 j ) f (ξ)e ik2 j ξ m 1 (2 j-1 ξ) φ (2 j-1 ξ) dξ . (3.20)
Repeating considerations of Section III.1, we construct the 4π-periodic function

m σ (ξ; k, j) = n∈Z Z m 1 (ξ + 4πn) σ (2 -j+1 (ξ + 4πn), k2 j ) χ [-2π,2π] (ξ + 4πn) . (3.21) Setting m σ (ξ; k, j) = 1 √ 2 ℓ b j k,ℓ e iℓξ/2 , (3.22)
we obtain the Fourier coefficients of m σ (ξ; k, j)

b j k,ℓ = n g n 1 4π 2π -2π σ (2 -j+1 ξ, k2 j ) e i(n-ℓ/2)ξ dξ . (3.23)
It is clear that the coefficients b j k,ℓ have the expected asymptotic behavior as ℓ → ∞,

b j k,ℓ = O ℓ -L-1 . (3.24)
Finally, we obtain the following algorithm for computing the wavelet coefficients in (3.18),

T σ f, ψ j k = ℓ b j k,ℓ s j-1 ℓ . (3.25)
This expression is similar to (3.14), except that the sum is no longer a convolution. Thus, strictly speaking, the algorithm in (3.25) is not a filter bank, since filter bank algorithms are usually understood to consist of convolutions.

III.3 Connection with BCR approach

It is reasonable to expect that a subclass of Calderón-Zygmund operators (see e.g., vol. 2 of [START_REF] Meyer | Ondelettes et Opérateurs[END_REF]) may be implemented numerically via filter banks. Let us consider the class of symbols S 0 1,1 , where σ ∈ S 0 1,1 satisfies

|∂ α ξ ∂ β x σ(ξ, x)| ≤ C(α, β)(1 + |ξ|) β-α . (3.26)
It was shown in [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] that in wavelet bases operators of this class may be represented by sparse matrices. All information is contained in the following set of coefficients

     α j kℓ = T ψ j k , ψ j ℓ β j kℓ = T φ j k , ψ j ℓ γ j kℓ = T ψ j k , φ j ℓ , (3.27)
which gives rise to the "NS-form", an alternative to the "S-form" consisting of the elements T ψ j k , ψ j ′ k ′ (see [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] for more details).

To explain the relation of the filter bank approach to that using NS-form, let us consider wavelets with good localization in the Fourier domain (e.g., Battle-Lemarié wavelets), so that for a given precision we need to consider "interaction" between scales which are immediate neighbors. In this case we may consider the simplified S-form where only interaction between neighboring scales is taken into account. Thus, for a given subspace W j , only its mappings from subspaces W j+1 , W j and W j-1 are significant, and these are subspaces of V j-2 . In this approximation we then consider the mapping V j-2 → W j which is exactly the one considered in

T ψ(ξ) = m T (ξ/4) φ(ξ/4) , (3.28)
where

m T (ξ) = m # (2ξ)m 0 (ξ) (3.29)
and m # is defined in (3.8) (we suppress the scale index j).

If more accuracy is required, one may consider mappings between more scales, e.g., W j+2 , W j+1 , W j , W j-1 and W j-2 , which amounts to considering the mapping V j-3 → W j . This corresponds to an approximation

T ψ(ξ) = mT (ξ/8) φ(ξ/8) , (3.30)
where

mT (ξ) = m# (4ξ)m 0 (2ξ)m 0 (ξ) (3.31)
and m# is defined similar to m # , except that the 4π-periodization in (3.8) is replaced with the 8π-periodization.

Let us then consider again a convolution operator T with symbol a(ξ). Following [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF], we consider

P 0 T P 0 = P J T P J + J j=1 [Q j T Q j + Q j T P j + P j T Q j ] (3.32) = P J T P J + J j=1 [Q j T P j-1 + P j T Q j ] . (3.33)
The action of first term Q j T P j-1 on f (x) ∈ L 2 (IR) (putting together P j and Q j is motivated by Lemma II.1) may be evaluated with the same type of approximation as discussed before,

Q j T P j-1 f (x) = k,ℓ s j-1 ℓ φ j-1 ℓ , T * ψ j k ψ j k (x) . (3.34)
Using notation of the previous section, we have

T φ j-1 ℓ , ψ j k = 1 2π 2 j-1/2 e i(2k-ℓ)2 j-1 ξ a(ξ) m 1 (2 j-1 ξ) φ 2 j-1 ξ 2 dξ (3.35) ≈ 1 2π 2 j-1/2 e i(2k-ℓ)2 j-1 ξ m j # (2 j-1 ξ) φ 2 j-1 ξ 2 dξ (3.36) = n b j n 1 2π e i(2k-ℓ-n/2)ξ φ(ξ) 2 dξ . (3.37)
Finally, using the definition of coefficients a 2ℓ-1 in Appendix, we obtain for the coefficient of

ψ j k (x) in (3.34), ℓ T φ j-1 ℓ , ψ j k s j-1 ℓ = ℓ b j 2(2k-ℓ) s j-1 ℓ + 1 2 ℓ,n a 2n-1 b j 2(2ℓ-k)+(2n-1) + b j 2(2ℓ-k)-(2n-1) s j-1 ℓ .
(3.38) Comparing this expression with equation (3.13) and interchanging the order of summation, we recognize here the same structure as that described in Section III.1. We may interpret the summation in (3.38) as interpolation to obtain the half-integer translates coefficients s j-1 k . Coefficients b k are given in (3.9).

Let us now turn to the P j T Q j term. We notice that in order to describe this term, it is sufficient consider only Q j+1 T Q j ,

P j T Q j f ≈ Q j+1 T Q j f ,
due to considerations above. This term represents mapping from scale j to scale j + 1, and

Q j+1 T Q j f (x) = k,ℓ d j k T ψ j k , ψ j+1 ℓ ψ j+1 ℓ (x) = ℓ q j ℓ ψ j+1 ℓ (x) . (3.39)
As before, we approximate the coefficient ω j k-2ℓ = T ψ j k , ψ j+1 ℓ as follows

ω j k-2ℓ = 1 2π 2 j+1/2 e i(2ℓ-k)2 j ξ a(ξ) m 0 (2 j-1 ξ) m 1 (2 j ξ) m 1 (2 j-1 ξ) φ 2 j-1 ξ 2 dξ ≈ 1 2π 2 j+1/2 e i(4ℓ-2k)2 j-1 ξ a(ξ) m j+1 # (2 j ξ) m 0 (2 j-1 ξ) m 1 (2 j-1 ξ) φ 2 j-1 ξ 2 dξ ≈ n b j+1 n c n+2k-4ℓ , (3.40)
where we have set for simplicity

m 0 (ξ) m 1 (ξ) = 1 2 n c n e inξ . (3.41)
Again, we obtain for the coefficients q j ℓ in (3.39) a filter bank type relation,

q j ℓ = k ω j k-2ℓ d j k . (3.42)
The results of this section may be summarized as follows. From equations (3.13), (3.38) and (3.42) we have

dj k = n b j n s j-1 2k-n/2 = ℓ b j 2(2k-ℓ) s j-1 ℓ + 1 2 ℓ,n a 2n-1 b j 2(2ℓ-k)+(2n-1) + b j 2(2ℓ-k)-(2n-1) s j-1 ℓ + k ω j-1 k-2ℓ d j-1 k , (3.43)
which is a filter bank type representation obtained using elements of the NS-form. The terms in (3.43) may be interpreted as follows: the first two terms of the r.h.s. of (3.43) may be viewed as representing the "half-integer samples" discussed in Section II via interpolation within V j-1 and the third term is an element of W j-1 so that (3.43) represents the mapping V j-2 → W j .

IV HILBERT TRANSFORM OF AUTOCORRELATION WAVELETS

In Section II we have described an approximation of the action of the Hilbert transform on wavelets. These approximate filters have to be applied to the coefficients of the function on subspaces V j . Since the discrete Hilbert transform is usually defined directly on the samples of the function, it is advantageous to require that the coefficients s j k are (at least approximately) the values of the function. This requirement may be satisfied by considering interpolating scaling functions. Examples of such scaling functions proposed in [START_REF] Saito | Multiresolution representations using the auto-correlation functions of compactly supported wavelets[END_REF] are obtained as autocorrelations of the usual compactly supported scaling functions. The properties of such autocorrelation wavelets and scaling functions are described in Appendix. We note that by using symmetric interpolating wavelets in this section, we give up orthogonality of the basis.

In view of the applications we consider further in the paper, we will elaborate on the case of the so-called dyadic wavelet transform (which is redundant with respect to the translation variable, see Appendix). As a result we obtain an O(N log N ) algorithm for decomposition and computing the Hilbert transform. Let us make clear that the redundancy may be avoided (yielding an O(N ) algorithm) if we follow considerations of Section II.

IV.1 Representation of the Hilbert transform

Let us consider m1 (ξ) = i sgn(ξ)|m 1 (ξ)| 2 (4.1)
and denote by mc

1 (ξ) its restriction to the interval [-2π, 2π]. Let m 2 (ξ) = ∞ k=-∞ mc 1 (ξ + 4πk), (4.2)
and consider its Fourier series,

m 2 (ξ) = 2i ∞ k=1 b k sin kξ 2 , (4.3)
where m 2 is a 4π-periodic function. The adverse effect of the restriction to [-2π, 2π] and of the 4π-periodization in (4.2) is weakened by the fact that the problematic point in multiplying by i sgn(ξ) is at the origin ξ = 0, where m 1 has a zero of order L. Therefore, the sequence b k in (4.3) has fast decay and, since Φ(ξ) is concentrated around ξ = 0, we may expect the product m 2 (ξ) Φ(ξ) to be a good approximation for sgn(ξ)|m 1 (ξ)| 2 Φ(ξ). Let Tj f = T j (Hf ) denote the jth scale of dyadic wavelet transform of the Hilbert transform of f . We prove Theorem IV.1 Let Ψ be the autocorrelation of the Daubechies' compactly supported wavelet with L/2 vanishing moments. Then the coefficients of Tj f = T j (Hf ) may be approximated by those obtained from the following pyramidal algorithm:

W j f (n) = ∞ k=1 b 2k-1 [S j f (n + k2 j-1 -2 j-2 ) -S j f (n -k2 j-1 + 2 j-2 )] , (4.4)
where the sequence b k is given by

b k =        0 , f or k = 2m, -1 (2m -1)π 1 - L/2 ℓ=1 a 2ℓ-1 1 1 -4( 2ℓ-1 2m-1 ) 2 , f or k = 2m -1, (4.5)
and decays at infinity as

b 2k-1 = O((2k -1) -L-1 ) . (4.6)
In addition, the exact and approximate Hilbert transform coefficients satisfy

|| Tj f -W j f || ≤ 2K 2 (1 + 2π) -2αL ||f || , (4.7)
where α is the Hölder regularity of the scaling function φ(x).

The theorem is proved below, and the proof follows the lines of that of Proposition II.1. We detail it for completeness.

IV.1.1 Computation of coefficients b 2k-1 and their behavior for large k

We have

b k = 1 2π 1 2i 2π -2π i sgn(ξ) sin kξ 2 |m 1 (ξ)| 2 dξ = 1 8π 2π -2π sin kξ 2 sgn(ξ)dξ - 1 8π L/2 ℓ=1 a 2ℓ-1 2π -2π sin kξ 2 cos ((2ℓ -1)ξ)sgn(ξ)dξ . (4.8)
The computation of the integrals in (4.8) yields (4.5).

If

k is large enough, 2k -1 > 2(2ℓ -1), then 1 1-4( 2ℓ-1 2k-1 )
2 may be replaced by its Taylor series, namely, 1

1 -4( 2ℓ-1 2k-1 ) 2 = ∞ p=0 2 2ℓ -1 2k -1 2p (4.9)
According to Lemma IX.2 of Appendix, the sequence {a 2ℓ-1 } has L-1 vanishing even moments, and we have

b 2k-1 = -1 π(2k -1) 2 ∞ p=L (2k -1) -2p L/2 ℓ=1 a 2ℓ-1 2 2p (2ℓ -1) 2p = -1 π (2k -1) -L-1 L/2 ℓ=1 a 2ℓ-1 [2(2ℓ -1)] L 1 -(2 2ℓ-1 2k-1 ) 2 = O((2k -1) -L-1 ) . (4.10)

IV.1.2 Pyramidal algorithm

Let us consider the approximate wavelet transform of Hf ,

W j f (n + w) = 1 2π f (ξ)e iξ(n+w) m 2 (2 j-1 ξ) Φ(2 j-1 ξ) dξ = ∞ k=1 b 2k-1 [S j-1 f (n + w + k2 j-1 -2 j-2 ) -S j-1 f (n + w -k2 j-1 + 2 j-2 )] .
(4.11) By setting w = 0, we arrive at (4.4).

Remarks:

1. As in the orthogonal case, we note that for j = 1 (and only in this case, since we are now using the dyadic wavelet transform) we need samples of S 0 f for half-integer n. Thus, an additional interpolation procedure is required for the first step of the algorithm, j = 1. An alternative is to set w = 2 j-2 , which would yield half integer samples of the Hilbert transform coefficients. We will come back to this point later on.

2. In order to compute the wavelet coefficients of f (x) with respect to the Hilbert transform of Ψ(x), we just need to use (4.4), since f, HΨ jk = -Hf, Ψ jk .

IV.1.3 Accuracy estimate

Let us introduce the following function

Θ(ξ) = m 2 ξ 2 Φ ξ 2 ≈ H * Ψ(ξ) (4.12)
where m 2 is the 4π-periodic function given in equation (4.2). The function Θ is also a wavelet, and has, in fact, the same number of vanishing moments as Ψ.

Considering f ∈ L 2 (IR), we fix j and compute

|| Tj f -W j f || 2 ≤ || Θ j -HΨ j || 2 ∞ ||f || 2 , (4.13)
and

| Θ j (ξ) -HΨ j (ξ)| 2 = |m 2 (2 j-1 ξ) -m1 (2 j-1 ξ)| 2 | φ(2 j-1 ξ)| 4 . (4.14)
Since we know that for some positive

K [5] | φ(ξ)| ≤ K[1 + |ξ|] -αL , (4.15)
and that m 2 -m1 = 0 inside the interval [-2π, 2π], it follows that

|| Θ j -HΨ j || 2 ∞ ≤ 4K 4 (1 + 2π) -4αL . (4.16)
This completes the proof of Theorem (IV.1).

IV.2 Numerical examples

The coefficients b 2k-1 are easy to compute numerically using equation (4.5). We present in Figure 1 a plot of the approximate filter m 2 (ξ) (up to a factor i) for the case of the autocorrelation of Daubechies' wavelet with 5 vanishing moments. As expected, we observe the sign flip in [-2π, 0]. We also provide tables of the twenty top coefficients b 2k-1 in (4.5) for the autocorrelation of Daubechies' compactly supported wavelets with L = 2, 4, 6, ...12 (the numerical values of the a 2ℓ-1 coefficients are listed in [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF]). The coefficients b 2k-1 have been computed using Mathematica TM . To check numerically the accuracy of our approach, we compare Θ(ξ) and Ψ(ξ) for positive values of ξ. The difference is plotted in Figure 2 for L = 2 and L = 5, respectively. In both cases only the twenty top coefficients b 2ℓ-1 have been considered for the evaluation of m 2 (ξ). Table 3 contains a numerical estimate for the absolute value of the error max ξ∈[0,π] Θ -Ψ , computed for L = 2 to L = 12.

n L = 2 L = 4 L =

IV.3 Improving the accuracy

It is clear from the estimate (4. [START_REF] Meyer | Ondelettes et Opérateurs[END_REF]) and examples in Figure 2 that the accuracy may be controlled by increasing the number of vanishing moments of wavelets. Also, the accuracy estimate may be improved by considering the restriction of m1 to a larger interval. Let mc 1 denote the restriction of m1 to the interval [-2 n π, 2 n π], and set

m 2 (ξ) = ∞ k=-∞ mc 1 (ξ + 2 n πk). ( 4
b k = 1 2 n+2 π 2 n π -2 n π sin (2 -n kξ) sgn(ξ)dξ - 1 2 n+2 π L/2 ℓ=1 a 2ℓ-1 2π -2π sin (2 -n kξ) cos ((2ℓ -1)ξ) sgn(ξ)dξ . (4.18)
The same computation as before yields

m 2 (ξ) = 2 i ∞ k=1 b 2k-1 sin ((2k -1)2 -n ξ), (4.19)
where

b 2k-1 = -1 (2k -1)π 1 - L/2 ℓ=1 a 2ℓ-1 1 1 -2 2n ( 2ℓ-1 2k-1 ) 2 , (4.20) and b 2k-1 = O((2k -1) -L-1 ) (4.21) for 2k -1 > 2 n (2ℓ -1).
Let us consider the function Θ defined by Θ(ξ) = m 2 ( ξ 2 ) Φ( ξ 2 ). For f ∈ L 2 (IR) and fixed j we have

W j f (m) = 1 2π f (ξ)e iξm Θ(2 j ξ) dξ = ∞ k=1 b 2k-1 [S j-1 f (m + k2 j-n -2 j-n-1 ) -S j-1 f (m -k2 j-n + 2 j-n-1 )] . (4.22)
The accuracy estimate is derived essentially as before. The only change is that 2π has to be replaced by 2 n π. We then obtain

|| Tj f -W j f || 2 ≤ 4K 4 (1 + 2 n π) -4αL ||f || 2 , (4.23)
as a generalization of (4.13) and (4.16).

V OTHER EXAMPLES

A particularly simple implementation (similar to that for the Hilbert transform) is possible for the convolution operators with non-oscillatory kernel considered in [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF]. The coefficients of filters ( similar to m 2 ) might be scale dependent. If the operator is homogeneous of some degree, then filters on different scales will differ only by a scaling factor. This section is devoted to the study of the action of various operators of differentiation and integration inluding those of fractional order. Throughout the section, we will only use the autocorrelation wavelets that we used for the Hilbert transform, leaving to the reader the (straightforward) computations in the orthogonal case, as described in Section II.

V.1 Derivative operators

The Hilbert transform is homogeneous of degree zero and, therefore, the operator (and, thus, m 2 ) is the same on all scales. Since derivative operators are homogeneous, the filter will be the same on all scales except for a scaling factor,

d n dx n Ψ jk = 2 -nj d n Ψ dx n jk . (5.1)
Thus, it is sufficient to evaluate the derivative operator on the function ψ. Our analysis is based on an approximation of ξ|m 1 (ξ)| 2 by the 4π-periodization of its restriction to [-2π, 2π]. Let us consider (for simplicity) the case n = 1, and the 4π-periodic function

m 3 (ξ) = ∞ k=-∞ (ξ + 4πk)|m 1 (ξ + 4πk)| 2 χ [-2π,2π] (ξ + 4πk) . (5.2)
Since m 3 is an odd function, we have

m 3 (ξ) = ∞ k=1 δ k sin kξ 2 . (5.3)
An explicit computation of the Fourier coefficients δ k yields

δ k =              - 2(-1) k k 1 - L/2 ℓ=1,ℓ =ℓ 0 a 2ℓ-1 1 1 -( 2(2ℓ-1) k ) 2 + 1 2k a k/2 if k = 2(2ℓ 0 -1) - 2(-1) k k 1 - L/2 ℓ=1 a 2ℓ-1 1 1 -( 2(2ℓ-1) k ) 2 otherwise.
(5.4)

As in the case of the Hilbert transform, the behavior of the Fourier coefficients is governed by the number of vanishing moments of the wavelet. Using Lemma IX.2 and the Taylor series expansion of δ k for k > 2(2L -1), we estimate δ k as

δ k = -2(-1) k k -L-1 L/2 ℓ=1 a 2l-1 (2(2ℓ -1)) L 1 -( 2(2ℓ-1) k ) 2 = O(k -L-1 ) (5.5)
For a given precision, the series (5.5) may be truncated. As an example, we provide in Tables 4 and 5 coefficients δ k for the autocorrelation of the Daubechies wavelet with 4, 5 and 6 vanishing moments. The shape of filter m 3 in the Fourier domain is shown in Figure 3. Table 5: Approximate coefficients δ k in (5.5) for the derivative of the autocorrelation wavelets with 4, 5 and 6 vanishing moments (continued).

-

V.2 Second order derivative

For the second order derivative we have to consider the Fourier expansion of m 4 (ξ), the 4πperiodization of

-ξ 2 |m 1 (ξ)| 2 χ [-2π,2π] (ξ) (5.6)
Since m 4 (ξ) is even, we have

m 4 (ξ) = ∞ k=0 ∆ k cos kξ 2 , (5.7)
where

                             ∆ 0 = - 2π 2 3 + a 2ℓ-1 2ℓ -1 ∆ k = - 8(-1) k k 2 + (-1) k ℓ =ℓ 0 a 2ℓ-1 1 (2ℓ -1 -k 2 ) 2 + 1 (2ℓ -1 + k 2 ) 2 + a 2ℓ 0 -1 4(-1) k k 2 + 2π 2 3 if k = 2(2ℓ 0 -1) = - 8(-1) k k 2 + (-1) k a 2ℓ-1 1 (2ℓ -1 -k 2 ) 2 + 1 (2ℓ -1 + k 2 ) 2
elsewhere.

(5.8) The derivation of (5.8) is straightforward and is similar to that of (5.5). The graph of -m 4 (ξ) is shown in Figure 4, together with that of ξ 2 .

V.3 Derivative of the Hilbert transform

The symbol of the composition of the Hilbert transform and differentiation is σ(ξ) = -i|ξ| and this operator appears prominently in the inversion of the Radon transform on the plane, e.g. in X-ray tomography. Setting the evaluation of the Fourier coefficients yields

m 5 (ξ) = σ(ξ)|m 1 (ξ)| 2 χ [-2π,2π] (ξ) = -i ∞ 0 µ k cos kξ 2 , (5.9)
                           µ 0 = π 2 µ 4k = 0 µ 2(2k-1) = - π 2 a 2k-1 µ 2k-1 = - 1 π      2 2k -1 2 - a 2ℓ-1 (2ℓ -1) 2 + 2k-1 2 2 (2ℓ -1) 2 -2k-1 2 2 2      . (5.10)
The graph of i m 5 (ξ) is shown in Figure 5, together with that of |ξ|.

V.4 Fractional derivatives

Contrary to what the name suggests, it is better to view fractional derivatives as integral operators since they are non-local in a way similar to the Hilbert transform. This non-local behavior is manifested in the Fourier domain by the action of the operator that "breaks" the function at ξ = 0. However, in the wavelet representations such operators are approximately local if applied to functions which do not have projections on the coarsest subspace. In particular, band-limited signals are an example of such class of functions. The projection on the coarsest subspace has to be treated separately (if necessary) and requires very few operations since the function is represented by a small number of samples.

Using the following definition of fractional derivatives,

(∂ α x f ) (x) = +∞ -∞ (x -y) -α-1 + Γ(-α) f (y)dy , (5.11)
where α = 1, 2 . . . (if α < 0, then (5.11) defines fractional anti-derivatives), we find its representation in the Fourier domain as

a(ξ) = e -iαπ/2 ξ α + + e iαπ/2 ξ α -, (5.12)
where ξ α + = ξ α for ξ > 0 and is zero otherwise, and ξ α -= |ξ| α for ξ < 0 and is zero otherwise.

Since ∂ α x Ψ jk = 2 -αj (∂ α x Ψ) jk , (5.13) -6 -4 -2 2 4 6 -3 -2 -1 1 2 3
Figure 6: Modulus of the approximate filter m 6 (ξ) in (5.14) for the derivative of order 1/2 of the Hilbert transform of autocorrelation of Daubechies' wavelet with 6 vanishing moments (dashed line: the symbol ξ 1/2 ).

it is sufficient to evaluate the operator on the function Ψ. We have, as before,

m 6 (ξ) = ∞ k=-∞ a(ξ + 4πk) |m 1 (ξ + 4πk)| 2 χ [-2π,2π] (ξ + 4πk). (5.14)
The Fourier coefficients of m 6 (ξ) are given by

γ ℓ = 1 4π 2π -2π a(ξ) |m 1 (ξ)| 2 e -iℓξ/2 dξ = 1 2π Re 2π 0 |m 1 (ξ)| 2 e -iαπ/2 ξ α e -iℓξ/2 dξ. (5.15)
Setting

u α k = 1 4π 2π 0 ξ α cos kξ + απ 2 dξ , (5.16)
we obtain

γ ℓ = u α ℓ -1 2 k a 2k-1 u α ℓ+2(2k-1) + u α ℓ-2(2k-1) .
(5.17)

Again, the decay of the γ ℓ coefficients is governed by the regularity of m 6 . Since the 2π-periodic function |m 1 (ξ)| 2 vanishes at ξ = 0 and ξ = 2π together with its derivatives of order up to L -1, one directly obtains the asymptotics of γ ℓ ,

γ ℓ = O(ℓ -L-1 ) (5.18)
As an example, we display in Figure 6 the derivative of order α = 1/2 using the autocorrelation of Daubechies' wavelets with 6 vanishing moments. We note that m 6 (ξ) is complex-valued.

V.5 Integration operators

Let us now consider integration operator with symbol

σ(ξ) = 1 iξ . (5.19)
The same procedure as before yields For the coefficients λ k , we obtain

m 7 (ξ) = 1 iξ |m 1 (ξ)| 2 = k λ k sin kξ 2 . (5.20)
λ k = -i 2π Si(kπ) - 1 2 ℓ a 2ℓ-1 [Si((k + 2(2ℓ -1))π) -Si((k -2(2ℓ -1))π)] , (5.21)
where

Si(x) = x 0 sin(y) y dy .
The graph of i m 7 (ξ) is shown in Figure 7, together with that of 1/ξ.

VI THE HILBERT TRANSFORM OF SIGNALS

We now turn to signal processing problems. The purpose of this Section is to illustrate one of the applications of our method for computing the Hilbert transform of a signal. As we shall see, it is interesting to work in a context in which the scaling function has vanishing moments, since the samples of the signal may then be identified (within a certain accuracy) with the coefficients of its projection onto some V j space2 . For this reason we shall use autocorrelation wavelets (other choices such as high order Battle-Lemarié wavelets or Coiflets, whose scaling function also possess vanishing moments, would do the job as well). The autocorrelation wavelets also offer an advantage of a trivial reconstruction formula (simple summation of wavelet coefficients over scales, see equation (9.31) in Appendix, the price to pay being an O(N log(N )) complexity).

VI.1 Band-pass signals

Let f ∈ C r (IR), r > L, and assume also that Hf ∈ C r ′ (IR), with r ′ > L. Then, according to (9.31), we have the following wavelet decompositions (we refer to (9.28) and (9.29) in Appendix for the description of our notation),

f (k) = S j 0 f (k) + O(2 j 0 (L-1) ) (6.1) = S J f (k) + J j=j 0 +1 T j f (k) + O(2 j 0 (L-1) ) (6.2) [Hf ](k) = S j 0 [Hf ](k) + O(2 j 0 (L-1) ) (6.3) = S J [Hf ](k) + J j=j 0 +1 T j [Hf ](k) + O(2 j 0 (L-1) ) (6.4)
But the Hilbert transform is anti self-adjoint,

T j [Hf ](k) = Hf, ψ jk = -f, [Hψ] jk , (6.5)
so that for -j 0 large enough we have

[Hf ](k) ≈ S J [Hf ](k) + J j=j 0 +1 W j f (k) , (6.6)
where W j f is defined in (4.4). We then obtain the following

Theorem VI.1 Let f ∈ C r (IR) be such that Hf ∈ C r ′ (IR), with r, r ′ > L. Then [Hf ](n) = S J [Hf ](n) + W j f (n) + O((1 + 2π) -αL ) . (6.7)
Therefore, as long as for a sufficiently sparse scale J the low-pass component S J f (k) of a signal f (k) can be neglected, the algorithm in (6.7) provides a good approximation of the Hilbert transform of f .

VI.2 Examples

Speech signal (or at least voiced speech) is an example of signals that may be modeled as superpositions of amplitude and frequency modulated components (see e.g. [START_REF] Mcaulay | Low rate speech coding based on the sinusoidal model[END_REF][START_REF] Maes | The Wavelet Transform in Signal Processing, with Applications to the Extraction of Speech Modulation Model Features[END_REF]). Although wavelet decompositions do not seem optimal for applying the analysis we have in mind to speech, let us consider it as an illustration.

In Figure 8 we show a half a second example of sound /one two/ sampled at 8kHz. In Figures 9 and 10 we show the "band-pass" component of the signal j T j f (n) and the corresponding approximate Hilbert transform j W j f (n), respectively. Figure 11 is a zoom of Figures 9 and 10 in which the real and imaginary parts of the reconstructed analytic signal,

Z f (n) = j (T j f (n) + iW j f (n)) , (6.8)
are represented. Finally, in Figure 12 we shows the squared modulus |Z f (n)| 2 , i.e. the square of the instantaneous amplitude of the signal (in the sense of Ville [START_REF] Ville | Théorie et Applications de la Notion de Signal Analytique[END_REF]). Notice that the instantaneous amplitude still has a lot of oscillations characteristic of the presence of many additive components in the signal within the considered frequency band.

VII ON THE REPRESENTATION OF SIGNALS BY LOCAL PHASES AND AMPLITUDES

It is well-known that an arbitrary continuous-time signal may be represented (e.g. using the method described in the previous Section) in terms of its local phase or its phase derivative, the instantaneous frequency, and local amplitude (following the pioneering work of J. Ville [START_REF] Ville | Théorie et Applications de la Notion de Signal Analytique[END_REF]).

More precisely, writing

Z f (x) = f (x) + i[Hf ](x) , (7.9)
we obtain an analytic function (the so-called analytic signal) which may be associated with the so-called canonical pair,

A f (x) = |Z f (x)| instantaneous amplitude , ω f (x) = arg Z f (x) instantaneous phase. (7.10)
The instantaneous frequency is then defined as

ν f (x) = 1 2π ω ′ f (x) . (7.11)
The purpose of such representation is to obtain the local phase and amplitude in the hope that they are much less oscillatory that the original signal (and then more easily compressible if the target application is compression). Moreover, in such a case, the instantaneous amplitude and frequency are often intimately connected with physical quantities. In general, however, the instantaneous frequency and amplitude may be as complicated as the signal itself. This is particularly clear in the example shown in the previous section (see Figure 12), where the global amplitude of the considered speech signal has fast oscillations.

An explanation of this fact is as follows. In the speech signal, a given phoneme may often be modelled as a superposition of short chirps, each having its own instantaneous frequency. It is then clear that there is no natural way of assigning a unique instantaneous frequency to such a phoneme, since the instantaneous frequency oscillates fast due to the interferences between chirps. An adequate description of the speech signal thus has to take into account this "multicomponent" character of the signal.

It is natural to expect that by splitting the frequency band, the amplitudes of the subbands will have slower oscillations, so that the representation of the subbands in terms of local phase and amplitude becomes useful. Moreover, if the considered subband "contains" one and only one of the chirps of the phoneme, standard approximations (see for instance [6, ?]) show that analytic signal provides a good approximation of the behavior of the component.

It turns out that in the discrete case such representation is quite easy to obtain from our approximate Hilbert transform algorithm. Indeed, the main aspect of our approach is to derive approximate expressions for the Hilbert transform of the wavelet Ψ(x), together with a fast algorithm for the computation of the corresponding coefficients.

As a by-product, our method yields a decomposition of band-pass signals as

f (n) = Re j Z j f (n) , (7.12)
where

Z j f (n) = T j f (n) + iW j f (n) (7.13)
may be thought of as a "discrete analytic subband" of the signal. Again, it is easy to obtain from such analytic subbands the local amplitudes and frequencies,

A j f (n) = |Z j f (n)| , ν j f (n) = 1 2π T j f (n)W j f ′ (n) -T j f ′ (n)W j f (n) A j f (n) 2 . (7.14)
Expression in (??) and (7.14) is a discrete approximation of the continuous expression for the instantaneous frequency of an analytic signal. In particular, it involves the derivatives of T j f and W j f which may be evaluated using the representation of the derivative in bases of compactly supported wavelets derived in [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF].

As an example in Figure 13, we illustrate the representation of the scale decomposition of the speech signal /one two/ sampled at 8kHz. Plots represent coefficients T j f (n) and the corresponding local amplitudes A j f (n), respectively. In computations, we used autocorrelations of Daubechies' wavelet with 9 vanishing moments. It turns out that each one of the T j f signals has a much simpler structure than the original signal itself, so that its local amplitude is a much more natural object than the global one. It is reasonable to expect that this method could be used as a method for compression (a non-linear compression scheme). Indeed, the local amplitudes and frequencies being slowly varying, are easier to compress. We also see a potential for feature extraction (for example speaker identification in speech processing, along the lines of [START_REF] Maes | The Wavelet Transform in Signal Processing, with Applications to the Extraction of Speech Modulation Model Features[END_REF]).

Let us stress that the results of this paper may be generalized to wavelet packet decompositions (see e.g. [START_REF] Wickerhauser | Adapted Wavelet Analysis from Theory to Software[END_REF]), where the wavelet basis appears as a particular case of a family (library) of orthonormal basis decompositions generated from a pair of quadrature mirror filters. For this purpose we have a triplet of filters m 0 (ξ), m 1 (ξ), m 2 (ξ) from which we may generate all wavelet packets and their (approximate) Hilbert transform. This permits us to look for the decomposition that is optimal in terms of information cost (within the above-mentioned non-linear compression scheme) for a given signal. We plan to address this problem elsewhere.

VIII CONCLUSIONS

We have described here a method for approximating the action of a class of operators on wavelets, and obtained several fast effective algorithms for the numerical evaluation of such operators in the form of filter banks. Such algorithms are easy to implement in both software and hardware. The operators under consideration are essentially convolution operators, i.e. operators characterized by a multiplyer in the Fourier domain, and we described a possible extension to some classes of pseudodifferential operators. Our method has to be thought of as an alternative to the NS-form approach in [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF]. Indeed, we demonstrate that we obtain the same results with comparable accuracy.

Our construction is illustrated using the "autocorrelation wavelets" described in [START_REF] Saito | Multiresolution representations using the auto-correlation functions of compactly supported wavelets[END_REF], but may be applied to any wavelet decomposition associated with quadrature mirror filters with several vanishing moments. It is worth noticing that the underlying approximation are quite close to those made in [START_REF] Muschietti | Pyramidal algorithms for littlewood-paley decompositions[END_REF] for finding approximate fast wavelet transform algorithms.

In the case of the Hilbert transform, the results of Section VII clearly indicate that the representation of signals by local amplitude and phase is in general not appropriate, since it yields a fast oscillating amplitude (and frequency). However, it appears likely that this difficulty may be avoided by associating amplitude and phase to the subbands of the signal. Such method may be developed in the framework of our approach, which we plan to address separately in regards to speech processing applications.

IX Appendix: Wavelet bases and Multiresolution Analysis

IX.1 Multiresolution Analysis

Let us introduce our notation. We will work in the context of multiresolution analysis [START_REF] Mallat | Multiresolution approximation and wavelets[END_REF] and [START_REF] Meyer | Ondelettes et fonctions splines[END_REF],

. . . ⊂ V 2 ⊂ V 1 ⊂ V 0 ⊂ V -1 ⊂ V -2 ⊂ . . . ⊂ L 2 (IR) , (9.1)
with the scaling function φ(x) and wavelet ψ(x), and define subspaces W j as orthogonal complements of V j in V j-1 , V j-1 = V j ⊕ W j . (9.2)

We set ψ j k (x) = 2 -j/2 ψ(2 -j x -k) φ j k (x) = 2 -j/2 φ(2 -j x -k) which may be viewed as upsampling followed by convolution.

(9.
In the signal processing part of the paper we use the dyadic wavelet transform, a "translation-invariant" version of multiresolution decompositions. In such a case, one considers all the integral translates of the wavelet and the scaling function, i.e. ψ jk (x) = 2 -j/2 ψ 2 -j (x -k) φ jk (x) = 2 -j/2 φ 2 -j (x -k) .

(9.11)

The corresponding coefficients of f ∈ L 2 (IR) are denoted by

d jk = f, ψ jk s jk = f, φ jk , (9.12)
and their numerical evaluation may be realized via an O(N log(N )) algorithm similar to (9.9), We denote by M be the number of vanishing moments of ψ(x), i.e., I R

x m ψ(x)dx = 0, m = 0, 1, .., M -1, (

and consider compactly supported wavelets and associated quadrature mirror filters. The filter coefficients G inherit the vanishing moments from the wavelet, namely, ℓ ℓ m g ℓ = 0, m = 0, 1, .., M -1. (9.15)

IX.2 Autocorrelation Shell

We also develop the algorithm for the Hilbert transform in the context of the autocorrelation shell described in [START_REF] Saito | Multiresolution representations using the auto-correlation functions of compactly supported wavelets[END_REF], where the wavelet and scaling functions are autocorrelations of the wavelet and scaling function associated with I. Daubechies' wavelets (see e.g. [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]). Let m 0 (ξ) and m 1 (ξ) be the 2π-periodic square-integrable quadrature mirror filters associated with an orthonormal basis of compactly supported wavelets. Let us denote by φ and ψ associated scaling function and wavelet, respectively. Let Φ(x) and Ψ(x) denote the autocorrelation functions of φ(x) and ψ(x). Clearly, Ψ(ξ) ∼ ξ L as ξ → 0, so that Ψ(x) has L vanishing moments. In addition, it is easy to see that the scaling function also has vanishing moments, x m Φ(x)dx = 0, m = 1, ..., L -1 .

(9.25)

We will need in the sequel the following properties of the sequence {a 2ℓ-1 }:

Lemma IX.2 1.

L/2 1 a 2ℓ-1 = 1.

2.

Figure 2 :

 2 Figure 2: Error of the approximation Θ(ξ)-Ψ(ξ) for the autocorrelation of Daubechies' wavelet with M = 2 and M = 5.

Figure 3 :Figure 4 :

 34 Figure 3: Approximate filter m 3 (ξ) in (5.2) for the derivative of the autocorrelation of the Daubechies wavelet with 6 vanishing moments

Figure 5 :

 5 Figure 5: Approximate filter im 5 (ξ) in (5.10) and the symbol |ξ| for autocorrelation of the Daubechies wavelet with 5 vanishing moments.

Figure 7 :

 7 Figure 7: Approximate filter im 7 (ξ) in (5.20) and the symbol 1/ξ for the primitive of the autocorrelation of the Daubechies wavelet with 5 vanishing moments.

Figure 8 :

 8 Figure 8: An example of sound /one two/ sampled at 8kHz.

Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 9101112 Figure 9: Real part of the reconstructed signal in (6.8).

Figure 13 :

 13 Figure 13: Scale decomposition of the /one two/ signal

3 ) and expand any f ∈ L 2 ( 5 )

 325 As usual, there exist 2π-periodic trigonometric polynomials m 0 (ξ) and m 1 (ξ) such thatφ(2ξ) = m 0 (ξ) φ(ξ) ψ(2ξ) = m 1 (ξ) φ(ξ),(9.6)where φ and ψ are the Fourier transforms of φ and ψ, e.g.,φ(ξ) = I R φ(x) e -iξx dξ . (9.7)Filters m 0 and m 1 satisfy the "exact reconstruction condition",|m 0 (ξ)| 2 + |m 1 (ξ)| 2 = 1 .(9.8)Let us denote by H = {h ℓ } L-1 ℓ=0 and G = {g ℓ } L-1 ℓ=0 the associated quadrature mirror filters. A direct consequence of equation (9.6) is the pyramidal algorithm for the computation of the coefficients requires O(N ) operations and may be viewed as a convolution followed by decimation (or downsampling). The same filters are used in the reconstruction algorithm (again an O(N )

h

  ℓ s j-1,k-2 j-1 ℓ d jk = ℓ g ℓ s j-1,k-2 j-1 ℓ . (9.13)

1 2 L/ 2 ℓ=1aa

 122 Both Φ(x) and Ψ(x) are interpolating functions, i.e.Φ(n) = Ψ(n) = δ n,0 .(9.18)Thus, we have for functions on the subspace spanned by {Φ(2j x -k)} k∈Z Z , f (x) = k f (2 -j k)Φ(2 j x -k).(9.19)If {h k } L-1 k=0 are the Fourier coefficients of m 0 (ξ), then one has the following two-scale difference equations:Φ(x) = Φ(2x) + 1 2ℓ-1 [Φ(2x -2l + 1) + Φ(2x + 2ℓ -1)] , 2ℓ-1 [Φ(2x -2l + 1) + Φ(2x + 2ℓ -1)] ,

L/ 2 1 2 L/2 1 a

 221 (2ℓ -1) 2m a 2ℓ-1 = 0, m = 1, ...L/2 -1.Proof:Since the function Ψ(x) has L vanishing moments, we have0 = Ψ(x)dx = Φ(2x)dx -1 2ℓ-1 [Φ(2x -2l + 1) + Φ(2x + 2ℓ -1)] dx , (9.26)which implies the first property. On the other hand, using (9.25) for 2, 4, ...2m -2, m < L/2,

Table 1 :

 1 Coefficients b 2ℓ-1 in (4.5), for L = 2 to 6.

	6

Table 2 :

 2 Coefficients b 2ℓ-1 in (4.5), for L = 8 to 12.

	18

Table 3 :

 3 Error max ξ∈[0,π] Θ -Ψ as a function of the number of vanishing moments, for L = 2 to 12.Here m 2 is a 2 n+1 π-periodic function. Let b k denote the Fourier coefficients of m 2 ,

Table 4 :

 4 Approximate coefficients δ k in (5.5) for the derivative of the autocorrelation wavelets with 4, 5 and 6 vanishing moments .

	n	L = 8	L = 10	L = 12
	1	2.784770784770785	2.793392366147784	2.799175787195709
	2	-0.67291259765625	-0.6661834716796875	-0.6615571975708008
	3	-0.819050230814937	-0.842451665981078	-0.858413908073351
	4	0.3888888888888888	0.4090909090909091	0.4230769230769232
	5	-0.3810738968633705	-0.4129362628218463	-0.4354273446748881
	6	0.08673095703125	0.0863571166992188	0.0855073928833008
	7	0.1838684191625368	0.2166304283951343	0.2410324766511379
	8	-0.0883838383838384	-0.1101398601398601	-0.1269230769230769
	9	0.0828421514920371	0.1104562019893828	0.1326570763346433
		-0.02213134765624998	-0.02583160400390625	-0.02801492214202881
		-0.03276200923259743	-0.05229010340572205	-0.06972013787429609
		0.0135975135975136	0.02447552447552448	0.034841628959276
		-0.01062078865282526	-0.02220087434526055	-0.03411181962890829
		0.003513881138392858	0.006190844944545201	0.00838465350014823
		0.002485393803852205	0.00813401608533453	0.015200246275297
		-0.0006798756798756801	-0.003239407651172366	-0.006876637294593952
		0.0002619928889691353	0.002425305600742754	0.006019732875689713
		-0.0001186794704861111	-0.00081475575764974	-0.001737650235493968
		0.00005945918191033914	-0.0005206695929446174	-0.002048813682806388
		-0.00003199414964121239	0.0001363961116283008	0.0007859014050964675
		0.00001817909841825415	-0.00005034211869672657	-0.000566641522074894
		-0.00001078904277146913	0.00002184781161221591	0.0001923626119440171
		6.637629002919957 • 10 -6	-0.00001049069169242624	0.0001128332173139814
		-4.209756531740421 • 10 -6	5.412544112237683 • 10 -6	-0.00002847468859045373
		2.740848684599229 • 10 -6	-2.95028230500094 • 10 -6	0.00001012721018176066

An alternative would be to use spline wavelets and the associated Lagrange interpolation to obtain the connection between approximation coefficients and samples, or to use the more general algorithms developed in[START_REF] Delyon | On the computation of wavelet coefficients[END_REF].

Acknowledgments.

The research of G.B. was partially supported by ARPA grant F49620-93-1-0474 and ONR grant N00014-91-J4037.

) 2m a 2ℓ-1 Φ(2x)dx , (9.27) which implies the second property.

If f ∈ L 2 (IR), we will denote by T j f and S j f the dyadic wavelet transform of f and the scaling function transform of f ,

We have

).

The two-scale difference equations imply that the computation of the S j f (n) and T j f (n) coefficients can be realized through the following pyramidal algorithm,

(9.30)

Moreover, the perfect reconstruction formula (9.8) yields the following simple reconstruction formula from the autocorrelation wavelet coefficients,

T j f (n) . (9.31)


René A Carmona

Wen Liang Hwang

Bruno Torresani

Bruno Torrésani

CHARACTERIZATION OF SIGNALS BY THE RIDGES OF THEIR WAVELET TRANSFORMS

come   L'archive ouverte pluridisciplinaire

Introduction

The characterization and the separation of amplitude and frequency modulated signals is a classical problem of signal analysis and signal processing. Applications can be found in many situations, such as for instance radar/sonar detection and speech processing [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF]. Many methods have been proposed in the past few years to analyze the time-frequency localization of signals. The most noticeable are the family of bilinear representations such as the Wigner representation and its generalizations (see [START_REF] Boashash | Estimating and Interpreting the Instantaneous Frequency of a Signal[END_REF][START_REF] Flandrin | Temps-Fréquence[END_REF] for a review) and the linear representations such as the wavelet and Gabor transforms.

In 1990 by the Marseille group proposed a new algorithm (see [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF] for a survey) based on the study of the phase of the wavelet (or Gabor) transform. The present work is an attempt to extend the latter to noisy situations. The main thrust of this correspondence is to use the localization properties of the modulus of the transform (which is generally more robust than the phase, eventhough the latter provides more precise estimates [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF]). In the case of frequency-modulated signals, the wavelet transform is "concentrated" in the neighborhood of curves (the ridges of the transform). We develop a scheme in which these curves are searched as such, in a (high dimensional) space of ridges, via a stochastic relaxation procedure. This alternate characterization of the ridges is better suited to the needs of noisy signal analyses. We also propose a stable method for signal reconstruction from the numerically computed ridges. This method is also based on an L 2 -minimization procedure.

For the sake of simplicity, our discussion is restricted to the case of the wavelet transform. But since our algorithms deal only with post-processing of time-frequency transforms, they can be extended to any time-frequency energetic representations. The case of the Gabor transform will be considered in the companion paper [START_REF] Carmona | Multiridge Detection and Time-Frequency Reconstruction[END_REF] where still another stochastic search algorithm, adapted to different situations, will be introduced.

We close this introduction with a short summary of the contents of the paper. Section 2 below is devoted to the statement of the problem and the definition of the ridges, Section 3 presents the main features of the variational problems which we propose and solve to estimate the ridges. We also give a Bayesian interpretation of this approach and we describe how one can modify the penalization functional in order to accommodate the presence of an additive noise in the signal. This section ends with a discussion of an example of a bat sound signal which we embed in noise. Section 4 is devoted to a quick account of the reconstruction problem. Namely, given a set of points in the time-scale domain, find the signal most likely to have a ridge going through these points. As before we outline the mathematical derivations and we illustrate the efficiency of the method on a numerical example.

Ridges

The first goal of this section is to set up an abstract formalism for the mathematical definition of the ridges of functions of two variables. Then, we propose two Monte Carlo algorithms to detect and identify these ridges. Let ψ ∈ L 1 (IR) be such that 0 < c ψ = ∞ 0 | ψ(ξ)| 2 dξ/ξ < ∞, i.e. fulfills the wavelet admissibility condition. The corresponding wavelet transform of f (x) is given by:

T f (b, a) = f, ψ (b,a) = 1 a ∞ -∞ f (x)ψ x -b a dx = e -ϕ ∞ -∞ f (x)ψ (e -ϕ (x -b))dx , (1)
where we have introduced the auxiliary variable ϕ = log(a). We are mostly interested in frequency modulated signals which can be written as the sum of finitely many components of the form:

f (x) = A(x) cos φ(x) , (2)
but for the purpose of the present correspondence we shall restrict ourselves to monocomponent signals.

See [START_REF] Carmona | Multiridge Detection and Time-Frequency Reconstruction[END_REF] for a detailed analysis of the multicomponent case. It is convenient to use the so-called "progressive wavelets", i.e. with vanishing negative frequencies. If ψ(x) is such a wavelet, then the wavelet coefficients of f (x) are given by:

T f (b, a) = f, ψ (b,a) = 1 2 Z f , ψ (b,a)
where Z f (x) is the "analytic signal" of f (x) given by: Z f (x) = 1 π P. f (x+y) dy y , where P denotes principal value integral. It is well-known [1] that a signal of the form (2) with A(x) and φ ′ (x) slowly varying gives Z f (x) ≈ A(x) exp{iφ(x)}. If we assume that the Fourier transform ψ(ξ) is peaked near a particular value ξ = ω 0 of the frequency, like for instance the Morlet wavelet given in the Fourier domain by ψ(ξ) = exp{-(ξ -2π) 2 /2}, it follows from standard arguments [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF] that the wavelet transform may be approximated as:

T f (b, a) ≈ 1 2 A(b) exp{iφ(b)} ψ(aφ ′ (b)) + O(|A ′ |/|A|, |φφ ′′ |/|φ ′ | 2 ) . (3)
From the localization properties of the wavelet in the Fourier domain one can see that the modulus |T f | of the wavelet transform is essentially maximum in the neighborhood of a curve a = a r (b) = exp ϕ(b), the ridge of the wavelet transform, related to the instantaneous frequency of the signal by a r (b) = e ϕr(b) = ω 0 /φ ′ (b). In [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF], the phase coherence of the wavelet transform was used to get a numerical estimate of the ridge. Since the phase can be somewhat difficult to control in noisy situations, we shall mainly focus here on the localization of the maxima of the modulus of the wavelet transform.

Ridge Detection: Variational Approaches

The purpose of this section is to give two examples of ridge detection algorithms both derived from variational problems. In both cases the ridge is searched in a high-dimensional space of curves and the ridge estimate appears as the graph of the argument of the minimization of a suitable penalty function.

Unlike the methods in [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF], the penalty function is mainly on the square modulus of the wavelet transform.

The loss of accuracy is weak, since the signals for which the methods are designed are supposed to have slowly varying frequencies. In the first case, the ridge is the graph of a function b → ϕ(b) while it is the graph of a parametric curve in the second case. The results of this section can be used beyond the single component case (the wavelet transform has a single ridge) when the ridges can be separated by a preprocessing localization procedure and then analyzed separately.

A Direct Search Algorithm:

We first assume that the ridge of the wavelet transform of the signal f can be parametrized by a function b ֒→ ϕ(b) defined for all the values of b. For the sake of the present discussion we denote by Φ the space of all the twice differentiable functions with square integrable derivatives. We then define the penalty function F f on the set Φ of ridge candidates ϕ by:

F f (ϕ) = -|T f (b, e ϕ(b) )| 2 db + λϕ ′ (b) 2 + µϕ ′′ (b) 2 db (4)
Such a penalty function clearly implements the two following features: the smoothness of the ridge and the localization in the time-frequency plane (for λ = µ = 0, minimizing F f (ϕ) is equivalent to searching maxima of |T f | 2 in the a direction). Our estimate of the unknown ridge of the wavelet transform of the signal f will be the function ϕ(b) which minimizes F f (ϕ). The Euler equation associated with this minimization problem can easily be obtained, and once discretized into a finite difference equations, solved numerically. However, such an approach is efficient only for weak noise. The presence of a strong noise component implies the existence of many local extrema in which the algorithm may get trapped.

We need a procedure which can jump over the local extrema to reach the global one(s). A natural candidate for this is the simulated annealing algorithm [START_REF] Van Laarhoven | Simulated Annealing: Theory and Applications[END_REF].

Snake Penalization

We now consider a ridge as a parametrized curve r :

s ∈ [0, 1] → r(s) = [ρ 1 (s), ρ 2 (s)
] in the time-scale plane. The ridge then takes the form of a "snake" (see [START_REF] Kass | Snakes: Active Contour Models[END_REF] for a description of the method in an image processing context). We use a cost function which takes into account the modulus of the wavelet transform, as well as additional terms needed in order to ensure the smoothness of the ridge (both in the b and a directions). We set:

F f (r) = -|T f (ρ 1 (s), ρ 2 (s)) | 2 ds + λ a ρ ′ 2 (s) 2 + µ a ρ ′′ 2 (s) 2 + λ b ρ ′ 1 (s) 2 + µ b ρ ′′ 1 (s) 2 ds , (5)
where λ a , λ b , µ a and µ b are positive constants. In the "snake terminology" of [START_REF] Kass | Snakes: Active Contour Models[END_REF] the second term is the "internal energy" of the snake. Its role is to control the smoothness and the rigidity of the snake. The first term is the "external energy" of the snake. It accounts for the interaction of the snake with the wavelet transform modulus. For the reasons mentioned in the previous section, we turn to stochastic optimization techniques (see [START_REF] Kass | Snakes: Active Contour Models[END_REF] for a direct solution of the corresponding Euler equations)

for the numerical solution of such a minimization problem.

Remark: In many applications the signal f (x) is the sum of a pure component f 0 (x) and a noise component n(x). When some information on the noise is available, it may be included into the penalty function (see e.g [3,[START_REF] Carmona | Multiridge Detection and Time-Frequency Reconstruction[END_REF][START_REF] Carmona | Practical Time/Frequency Analysis: Wavelet and Gabor Transforms with an Implementation in S. (book under review) Swave: a Public Domain Package of S Functions for Time-Frequency Analysis[END_REF] for more details on this point).

Bayesian Interpretation:

Both ridge detection procedures have a Bayesian interpretation.

Let us present it in the case of the direct search. Consider the prior probability measure defined formally by µ prior (dϕ

) = Z -1 1 exp{-µ a |ϕ ′ (b)| 2 + λ a |ϕ ′′ (b)| 2
db} ′′ dϕ ′′ and the conditional probability

µ ϕ (df ) = Z -1 2 e |T f (b,ϕ(b))| 2 ′′ df ′′ .
which gives the probalility, conditioned by ϕ, that the signal is in the infinitesimal ′′ df ′′ in the space of finite energy signals. Then according to Bayes' rule, the conditional probability knowing the signal is given by µ posterior (dϕ|f ) = exp -F [ϕ] ′′ dϕ ′′ /Z for some constant Z.

Maximizing µ posterior (dϕ|f ) is equivalent to minimizing (4).

Cost Minimization by Simulated Annealing:

We included the ridge detection procedures described above in a package of S functions made available on the Internet [START_REF] Carmona | Practical Time/Frequency Analysis: Wavelet and Gabor Transforms with an Implementation in S. (book under review) Swave: a Public Domain Package of S Functions for Time-Frequency Analysis[END_REF]. The implementation was done by solving the variational problems by simulated annealing (see [START_REF] Van Laarhoven | Simulated Annealing: Theory and Applications[END_REF] for background on this combinatorial optimization technique). The details are spelled out in [3] and the book [START_REF] Carmona | Practical Time/Frequency Analysis: Wavelet and Gabor Transforms with an Implementation in S. (book under review) Swave: a Public Domain Package of S Functions for Time-Frequency Analysis[END_REF].

Examples:

Numerical experiments have been made on various types of academic and real signals.

We illustrate the method described above with a (real) sonar signal emitted by certain species of bats.

The signal is frequency modulated, with approximately hyperbolic instantaneous frequency. The wavelet transform (w.r.t. Morlet's wavelet) was computed for frequencies ranging from ν s /16 to ν s /2, with ν s the sampling frequency, in geometric progression (i.e. of the form a = 2a n/20 0

, n = 0, . . . 59). Figure 1 shows the wavelet transform of the signal (left) and the wavelet transform of the same signal with additive Gaussian white noise, with input SN R = -5dB. Superimposed are the ridges estimated with the direct search procedure. Both transforms are coded with gray levels proportional to their modulus square. In Figure 2, we show the ridge estimated with the snake procedure (left, notice that the boundaries have been fairly well reproduced), and a comparison of ridge estimations in various situations (right).

Reconstruction from the Skeleton on a Ridge

We present in this section a new algorithmic reconstruction of a signal from the knowledge of sample values of its wavelet transform on the ridges of its modulus. For the sake of simplicity we restrict ourselves to the case of a single ridge. See nevertheless [START_REF] Carmona | Multiridge Detection and Time-Frequency Reconstruction[END_REF] for the analysis (in the case of the Gabor transform) of the more general case of finitely many arbitrary ridges. Let us focus on ridges given in the form b → ϕ(b). In practical applications one only knows sample points (b 1 , a 1 ), . . . , (b n , a n ) and the smooth function b ֒→ ϕ(b) which we use in lieu of the true (unknown) ridge function is merely a guess which one constructs from the sample points. We use a smoothing spline (but any other kind of nonlinear regression curve would do as well). From now on, ϕ(b) is a smooth ridge function which is constructed from the n sample data points.

Statement of the Problem:

We are concerned with the implementation of the folk belief that a signal can be characterized by the values of the restriction of its wavelet transform to its ridges.

Illustrations can be found in [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF], where it is shown that in the case of signals of the form (2) the restriction of the wavelet transform to its ridge of the wavelet transform behaves as A(x) exp[iφ(x)] (see also [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF] for similar remarks for the Gabor transform in the context of speech, yielding good quality reconstruction with high compression rate). Such an approach can be used in non-noisy situations, but it does fail in the presence of a significant noise component. We assume that the values of the wavelet transform, say z j , are known at sample points (b j , a j ). The set of sample points together with the values z j constitutes the wavelet transform skeleton of the signal to be reconstructed. We look for a signal f (x) of finite energy whose wavelet transform has the graph of the function a r (b) as ridge and satisfies:

T f (b j , a j ) = z j , j = 1, . . . , n . (6)
4.2 The Penalization Approach: We use a strategy which was successfully used by Mallat and Zhong to reconstruct a signal from the extrema of its dyadic wavelet transform [START_REF] Mallat | Characterization of Signals from Multiscale Edges[END_REF][START_REF] Carmona | Spline Smoothing & Extrema Representation: Variations on a Reconstruction Algorithm of Mallat and Zhong[END_REF]. In the present setting, we look for a signal f (x) which satisfies the constraints (6) while the L 2 -norm in the scale variable a of the modulus is kept to a minimum for each b. This may be achieved by minimizing

F 1 (f ) = 1 c ψ db da |a| |T f (b, a)| 2 (7)
(note that when the integration is performed on the whole half-plane, F 1 (f ) = ||f || 2 by the energy conservation formula). Since the cost function F 1 (f ) is a quadratic form in the unknown function f , the solution is easily computed by means of Lagrange multipliers. A solution can be constructed as a linear combination of the wavelets ψ (b j ,a j ) at the sample points of the ridge, the coefficients being given by the solution of a n×n linear system. This solution is not completely satisfactory, especially when the number of sample points is small. It ignores the empirical fact that (in most of the practical cases) the restriction to the penalty. Unfortunately, F2 (f ) is not quadratic in f . In order to remedy this problem we remark that according to the analysis of [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF], d db Ω(b, a r (b)) ≈ F ′ (b) = ω 0 /a r (b) , where Ω f (b, a) = arg T f (b, a). Then we replace F2 (f ) by a quadratic form which gives a good approximation of it, and F with

F (f ) = 1 c ψ db da |a| |T f (b, a)| 2 + ǫ bmax b min | d db T f (b, ϕ(b))| 2 - ω 2 0 ϕ(b) 2 |T f (b, ϕ(b))| 2 db = Qf, f . (8)
4.3 Solution of the Optimization Problem: The constrained minimization problem can be solved using Lagrange multipliers. The solution is given by

f (x) = 2n j=1 λ j Q -1 ψ j (x) ( 9
)
where Q is the operator (matrix after discretization of the problem) defined in [START_REF] Kass | Snakes: Active Contour Models[END_REF] and the functions ψ j are defined by: ψ j (x) = a -1 j ψ ((x -b j )/a j ) , j = 1, . . . , n, The Lagrange multipliers are determined by imposing the constraints [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF]. This gives a system of (2n) × (2n) linear real equations from which the Lagrange multipliers λ j 's can be computed.

Examples:

To illustrate the reconstruction procedure, we selected a subset of n = 500 consecutive samples from the bat signal (Figure 3, left). We used 40 sample points on the estimate of the ridge and the value ǫ = .5 to reconstruct the signal. The result of the reconstruction is given at the bottom of Figure 3. As may be seen, the reconstruction is of extremely good quality. An analysis (not presented here) of the modulus of the wavelet transform of the reconstructed signal shows that, because we chose a ridge estimate which ignored the existence of a secondary ridge, the latter is not present in the reconstruction. Further results (see [START_REF] Carmona | Multiridge Detection and Time-Frequency Reconstruction[END_REF][START_REF] Carmona | Practical Time/Frequency Analysis: Wavelet and Gabor Transforms with an Implementation in S. (book under review) Swave: a Public Domain Package of S Functions for Time-Frequency Analysis[END_REF]) confirm the quality of the reconstruction method. This justifies a posteriori the approximation we made in the definition of the quadratic penalty function.

Conclusions

We presented a new approach to the problem of ridge detection in an energetic distribution of a signal. Our approach is based on the minimization of a penalty function on the set of all possible ridge candidates. The penalty function takes into account a-priori information on the signal (namely the time-frequency representation of the signal which is essentially localized around a curve, this curve is smooth,...) and possibly on the noise (through an a-priori noise model, or simulations). The minimization is achieved through Monte-Carlo type methods. We also proposed a new synthesis procedure which requires only a small number of values of the transform on the ridge and which is very robust to noise. We have focused here on the case where the time-frequency representation is given by the square modulus of the wavelet transform (the scale variable being interpreted as an inverse frequency variable). Any time-frequency representation can be used as well. The case of the the square modulus of the Gabor transform will be considered in a forthcoming publication. [START_REF] Carmona | Multiridge Detection and Time-Frequency Reconstruction[END_REF].

  of the modulus |T f (b, a)| to the ridge, i.e. the function b ֒→ |T f (b, a r (b))|, is smooth and slowly varying. In order to force the solution of the constrained optimization problem to respect this requirement, we introduce the extra term F2 (f ) = bmax b min d db |T f (b, a r (b))| 2 db, and consider the minimization of the cost function F (f ) = F 1 (f )+ǫ F2 (f ) where the free parameter ǫ > 0 is chosen to balance the two contributions

Figure 1 :

 1 Figure 1: Left) Intensity plot of the modulus square of the wavelet transform of the bat signal. The ridge superimposed. Right) Ridge estimate, annealing method; SN R = -5dB.

Figure 2 :Figure 3 :

 23 Figure 2: Left) Ridge estimate, snake method; SN R = -5dB. Right) Comparison of several ridge estimates for the bat sonar signal. Notice that the scale is now increasing upward as opposed to all the other plots where it is increasing downward.

* Partially supported by ONR N00014-91-1010 † Supported by NSF IBN 9405146 0