René A Carmona

Wen L Hwang

Bruno Torrésani

MULTI-RIDGE DETECTION AND TIME-FREQUENCY RECONSTRUCTION

The ridges of the wavelet transform, the Gabor transform or any time-frequency representation of a signal contain crucial information on the characteristics of the signal. Indeed they mark the regions of the time-frequency plane where the signal concentrates most of its energy. We introduce a new algorithm to detect and identify these ridges.

The procedure is based on an original form of Markov Chain Monte Carlo algorithm specially adapted to the present situation. We show that this detection algorithm is especially useful for noisy signals with multi-ridge transforms. It is a common practice among practitioners to reconstruct a signal from the skeleton of a transform of the signal (i.e. the restriction of the transform to the ridges). After reviewing several known procedures we introduce a new reconstruction algorithm and we illustrate its efficiency on speech signals.

I. Introduction

A wide class of signals may be conveniently described in terms of time-dependent amplitude and frequency (see for example [START_REF] Flandrin | Temps-Fréquence[END_REF] or [START_REF] Boashash | Estimating and Interpreting the Instantaneous Frequency of a Signal[END_REF]) or sums of such amplitude and frequency modulated components (like in speech analysis and synthesis applications [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF], [START_REF] Quatieri | Phase coherence in speech reconstruction for enhancement and coding applications[END_REF], [START_REF] Mcaulay | Low Rate Speech Coding Based on the Sinusoidal Model[END_REF]). But the main problem remains the numerical evaluation of these time-dependent characteristics. Time-frequency representations [START_REF] Flandrin | Temps-Fréquence[END_REF] offer a convenient set up and the problem of the estimation of the local amplitude and frequency is well understood in the noise free case with only one component, the situation gets much more complex in the presence of noise and/or of several components. We proposed in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF] a new constrained optimization approach for processing one-component signals in very noisy situations. At high signal to noise ratios (SNR's for short), one-component signals can be analyzed by means of their instantaneous amplitude and frequency for which the estimation theory is well developed [START_REF] Boashash | Estimating and Interpreting the Instantaneous Frequency of a Signal[END_REF].

The nonparametric approach used in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF] was based on the solution of variational problems. This choice was motivated by the desire to handle low SNR's. We also note that bilinear representations such as the (generalized) Wigner distributions [START_REF] Flandrin | Temps-Fréquence[END_REF], [START_REF] Boashash | Polynomial Wigner-Ville Distributions and Time-Varying Higher Order Spectra[END_REF] can yield extremely precise results in the one-component case, but may completely fail in the multicomponent situation because of the presence of interference terms. We describe here a new approach capable of handling multiple component signals. This new approach is based on a new Markov Chain Monte Carlo (MCMC for short) which uses the energetic distribution provided by a timefrequency representation of the signal. Given the fact that the energy of the signal concentrates around curves in the time-frequency plane which we shall call "ridges", the Markov chain is constructed in such a way that the random walkers (hereafter called "crazy climbers") are attracted by these one-dimensional structures. We also add another step to the analysis. It is called "synthesis step" because it is devoted to the reconstruction of the "part" of the signal that led to the ridge. Throughout this paper, the discussion will be restricted to the cases of the Gabor and wavelet transforms.

The reason is our desire to consider applications to speech signals for which the Gabor transform is better suited. It is important to notice that since our detection algorithm is only a special post-processing of a timefrequency transform, it can be used with other time-frequency energetic representations, such as the family of Wigner distributions for example. On the contrary, the reconstruction algorithm is specific to the considered representations. We develop it for wavelet and Gabor transforms, but the modifications required to extend it to other linbear time-frequency representations are straightforward.

The thrust of the present paper is twofold. First we give a new ridge detection algorithm which can efficiently detect multiple ridges in the modulus of a transform and second we propose a signal reconstruction procedure from the knowledge of the skeleton of its transform on arbitrary points of its ridges. Our detection procedure is based on an original Markov Chain Monte Carlo algorithm. It is designed in such a way that weighted occupation densities draw the ridges on the time-frequency plane. Since its nature is post-processing of a time-frequency energetic representation, is not restricted to linear transforms (in particular it can be used with bilinear representations as discussed in [START_REF] Auger | Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method[END_REF]) and it can draw as many ridges as needed. But most importantly, its robustness to noise is remarkable. The reconstruction procedure is based on the classical idea of spline smoothing as presented in [START_REF] Wahba | Spline Models for Observational Data[END_REF]. It is restricted to linear transforms and it was used in the case of the wavelet transform (without much explanation) in the companion correspondence [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF]. We present it in full details here.

As for the ridge detection, it performs very well in noisy situations. Both components of our work (ridge detection and reconstruction) are illustrated on two specific data sets. The first one is the superposition of a real life sonar bat signal with an artificially generated chirp. It is analyzed with the wavelet transform.

We close this introduction with a short summary of the contents of the paper. After a short section devoted to notation and discretization issues, our detection algorithm is presented in a general setting in Section III.

We also try to emphasize the similarities and the differences with known procedures such as minimization by simulated annealing or the more recent "re-assignment" or "re-allocation" procedures advocated in [START_REF] Auger | Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method[END_REF] or [START_REF] Maes | [END_REF].

Section IV contains a detailed discussion of a first numerical examples and a discussion of the specifics of noisy signals. Section V gives our reconstruction procedure of the original signal from the estimates of the transform on the ridges. Some of the details of the penalization procedure are postponed to the Appendix. Because we choose to illustrate the efficiency of this reconstruction on speech signals, a short Section VI is devoted to the specifics of the sinusoidal model for speech.

All the detection and reconstruction algorithms presented in this paper have been implemented in the Splus environment. The data files and the S-code needed to produce all the numerical results and figures given as examples in this paper are 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].

II. Notation for the Continuous Gabor and Wavelet Transforms

We set the stage by introducing the time-frequency representations which we use to illustrate the crazy climber algorithm. The latter may be used as post processing of any time frequency representation, we shall restrict the present discussion to the cases of the continuous wavelet and Gabor transforms for the sake of simplicity. Indeed, the behavior of these transforms is easy to understand for amplitude and frequency modulated signals.

We work with the complex Hilbert space L 2 (IR) of square-integrable functions. Our convention for the Fourier transform is:

f (ξ) = ∞ -∞ f (x)e -iξx dx (1)
and consequently the Plancherel formula reads

|| f || 2 = 2π||f || 2 .

A. Continuous Transforms

Let ψ(x) be a fixed integrable function such that:

0 < c ψ = ∞ 0 | ψ(ξ)| 2 dξ ξ < ∞ . (2)
Such a function ψ is called an analyzing wavelet and the wavelet transform of a signal f ∈ L 2 (IR) with respect to ψ is defined by:

T f (b, a) = 1 a ∞ -∞ f (x)ψ x -b a dx. (3)
Throughout this paper, we use the Morlet wavelet ψ(x) = e -x 2 /2 e iω0x , but any complex-valued wavelet belonging to the complex Hardy space H 2 (IR), i.e. such that ψ(ξ) = 0 ∀ξ ≤ 0, would do as well. We use the notation:

ψ (b,a) (x) = 1 a ψ x -b a (4)
for the wavelet with scale a and location b. In addition, we shall find convenient to introduce the auxiliary variable ϕ = log(a). With such notation, it follows directly from Taylor's formula that if we consider signals of the form:

f (x) = N k=1 A k (x) cos(φ k (x)) (5)
where the amplitudes A k (x) are continuously differentiable and the phases φ k (x) are twice continuously differentiable, then their wavelet transform can be written in the form: Next we describe the case of the Gabor transform (although Gabor's original representation was discrete, we use a continuous version which we still call the Gabor transform). The Gabor transform of a signal f ∈ L 2 (IR) of finite energy is defined as:

T f (b, a) = 1 2 N k=1 A k (x)e iφ k (b) ψ(aφ ′ k (b)) + r(b, a) (6
G f (b, ω) = ∞ -∞ f (x)g(x -b)e -iω(x-b) dx , (7)
where g(x) is a window function with a good time-frequency localization. We use the Gaussian window , where s is a scale parameter, but other choices such as the Hamming windows (which are very popular in speech processing) would be as convenient. We shall use the notation

g(x) = g s (x) = 1 s √ 2π e -x 2 /2s
g (b,ω) (x) = g(x -b)e iω(x-b) (8)
for the time-frequency atoms used in the definition of the Gabor transform. The same argument as before shows that the continuous Gabor transform of signals of the type (5) may be written in the form:

G f (b, ω) = 1 2 N k=1 A k (x)e iφ k (b) ĝ(φ ′ k (b) -ω) + r(b, ω). (9)
Again the remainder term r(b, ω) depends upon the derivatives of the amplitudes and the local frequencies.

Assuming for simplicity (this is the case for the Gaussian windows as well as for the Hamming windows) that the Fourier transform of the window has fast decay away from the origin of frequencies, we end up again with a Gabor transform square modulus M (b, ω) = |G f (b, ω)| 2 exhibiting a certain number of ridges.

B. Discretizations

In practice, we have only access to discrete data, and the wavelet and Gabor transforms have to be discretized.

The discretization of the Gabor transform is traditionally a regular sampling of the continuous formulas. In practice, we sample the Gabor transform at the same rate as the signal.

Sampling the wavelet transform is somewhat more problematic. Indeed, the classical discretization scheme for wavelet transforms involves the so-called dyadic grid, of the form (kb 0 a j 0 , a j 0 ), j, k ∈ Z Z. Such a choice is not convenient for our purpose, because it is not shift-invariant: the discrete transform of a shifted copy of the signal is not the shifted copy of the original signal's transform. To overcome such a drawback, it is convenient to introduce more redundancy into the transform, and work with discretization grids of the form (kb 0 , a j 0 ), j, k ∈ Z Z. Or, in terms of the auxiliary variable ϕ, we will use a grid of the form (kb 0 , jϕ 0 ), j, k ∈ Z Z, with ϕ 0 = log(a 0 ), i.e. a regular grid again.

III. Drawing the Ridges of a Surface

As we have seen, characterization of the signal's instantaneous frequency by the wavelet or the Gabor transform can be achieved by extracting numerically the ridges as the set of local maxima of the transform's square modulus, denoted hereafter by M (b, ϕ), where ϕ is either the log of the scale variable (wavelet case) of the frequency variable (Gabor case). The problem which we are addressing in this section is to determine the ridges of the surface M (b, ϕ). When the surface is the energy landscape over the time-frequency plane given by a specific transform, a second challenging problem is, for each ridge, to estimate the corresponding component of the signal, and reconstruct it. We shall address this second problem in Section V below.

According to [START_REF] Hall | Ridge Finding from Noisy Data[END_REF], ridges may be defined in a very general setting as curves on a surface z = f (x, y). They are completely characterized by their projection on the (x, y) plane. We also call these curves ridges by a convenient abuse of language. We shall be interested in the special class of ridges which we describe below.

We start with a subset D of the upper (time-frequency) half plane. D will be bounded in the applications but we can think of D as the whole upper half-plane for the purpose of the present discussion. We shall use the notation (b, ϕ) for the points of the domain D. We consider a nonnegative function M (b, ϕ) defined on D.

We define the ridge set R as the set of local maxima in ϕ of the functions ϕ ֒→ M (b, ϕ) when the variable b is held fixed. We assume that the surface M (b, ϕ) is smooth enough so that the ridge set is the finite union of the graphs of smooth functions slowly varying on their respective domains. In other words we assume that:

R = ∪ L ℓ=1 R ℓ (10)
where each R ℓ is the graph of a smooth function [b ℓ,min , b ℓ,max ] ∋ b ֒→ ϕ ℓ (b) defined on a (possibly strict) subset of the domain of the variable b. In the practical applications which we have in mind, the ridge functions ϕ ℓ (b) are slowly varying. But notice that we do not make any assumption on the lengths of the individual ridges R ℓ or even the fact that they could cross 1 .

A. The Crazy Climber Detection Algorithm

The main idea of the "crazy climber" algorithm is as follows. A large number of particles (the climbers) are initially randomly seeded on the domain D at step 0. Then each climber evolves according to a Markov chain on D with a transition mechanism depending upon the local values of the M (b, ϕ) function. This chain is designed to relax to a steady distribution which is essentially concentrated on the ridges. The projection of the motion on the b-axis is the standard symmetric random walk, say with elastic reflection at the boundary points of the interval and consequently, the projection of the steady distribution on the b-axis is the uniform distribution. Vertically, the climbers are then encouraged to "climb on the hills" to reach the ridges by a

Hastings-Metropolis penalization and a temperature schedule similar to the simulated annealing algorithm.

But contrary to simulated annealing, the crazy climber algorithm looks for all the local maxima when the time variable is fixed instead of searching only for the global maxima. In fact, a natural implementation of the simulated annealing would lead to the simulation of one single sample path of a Markov chain in the very complex space of all the ridge candidates (see for example [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF] for an implementation in this spirit) while the crazy climber detection algorithm requires the simulation of many sample paths of a Markov chain in the timefrequency plane (or more precisely its discretized version) which is a state space of a much simpler complexity indeed.

Obviously the implementation involves a discretized version of the time-frequency plane. We assume that the time interval over which the signal is analyzed is discretized into a finite set

{b 0 , b 1 , • • • , b B-1 } with B elements.
We also assume that the values of the frequency variable ϕ are discretized into a finite set {ϕ

0 , ϕ 1 , • • • , ϕ K-1 }.
1 In fact, in the cases under consideration, namely wavelet and Gabor transforms, ridges never cross; even if one generates a signal whose analytical expression is the sum of two frequency modulated components with intersecting frequency modulation curves, the ridges will not reproduce exactly the frequency modulations near their expected intersection. We have no simple explanation for this "experimental" fact. Notice that such a remark is specific to wavelet or Gabor transforms. Indeed, ridges of Wigner-Ville representations do cross.

We thus reduce the analysis of the square modulus of the transform to the analysis of a finite B × K matrix with nonnegative entries.

A.1 Crazy Climbers

At time t = 0 we initialize the positions X α (0) of N climbers on the grid Γ

= {0, • • • , B -1} × {0, • • • , K -1}.
The climbers are labeled by the parameter α = 1, • • • , N . The initial positions are chosen independently of each other, uniformly over the grid Γ. The climbers evolve independently of each other according to the same law. If a climber is at the point (j, k) at time t, i.e. if X α (t) = (j, k), then its position at time t + 1, say

X α (t + 1) = (j ′ , k ′
) is determined according to the following law: j ′ = j -1 with probability 1/2 and j ′ = j + 1 with probability 1/2 (we do not discuss the particular cases j = 0 and j = B -1 involving boundary conditions not to confuse the issue). Then when the climber has decided to move to the left (when j ′ = j -1) or to the right (when j ′ = j + 1) in the horizontal direction, a possible vertical move is considered. As for the horizontal component, the climber tries to move up (i.e. k ′ = k + 1) or down (i.e. k ′ = k -1) with equal probabilities. Again we ignore the boundary conditions for the sake of simplicity. Unlike in the case of the horizontal direction, the move does not always take place. The transition from (j ′ , k) to (j ′ , k ′ ) takes place if the value of the function increases, i.e. if its so-called Delta ∆M = M (j ′ , k ′ ) -M (j ′ , k) is nonnegative. On the other hand, the move does not necessarily take place if the function decreases, i.e. if ∆M < 0. Indeed, in this case the transition is made i.e. X α (t + 1) = (j ′ , k ′ ) with probability exp[∆M/T (t)] and the climber does not move vertically, i.e.

X α (t + 1) = (j ′ , k) with probability 1 -exp[∆M/T (t)].
At each time t we consider two occupation measures. The first one is defined by:

µ (0) t = 1 N N α=1 δ Xα(t) .
It is obtained by putting a mass 1/N at the location of each of the climbers on the grid. In other words, µ

t (A) is the proportion of climbers in A at time t. The second one is the "weighted" occupation measure µ t defined by:

µ t = N α=1 M (X α (t))δ Xα(t)
and obtained by putting a mass equal to the value of the function M at the location of the climber.

We finally consider the corresponding "integrated occupation measures", defined by ergodic averages as follows:

µ 0 I = 1 T T t=1 µ 0 t ( 11
)
and

µ I = 1 T T t=1 µ t ( 12
)
September 21, 1998 DRAFT

The occupation measure µ 0 I is only given here for the sake of completeness. Indeed its main shortcoming is the fact that it assigns nonzero mass to regions without ridges if the lengths of the ridges are smaller than the length of the window. This is due to the very nature of the unrestricted horizontal motions of the climbers. Because the modulus of the denoised versions of the functions M which we use in the applications are essentially zero away from the ridges, the occupation measure µ I gives much better results when it comes to detecting ridges.

Further Remark: Notice that the climbers evolve independently of each other without interaction and with the same distribution so that the computer code generating the motions of the climbers is the same for all the climbers. This is an indication that the algorithm can naturally be parallelized on a SIMDIM machine (such as for example the massively parallel computers MASPAR I & II). We shall not report of such an implementations here.

A.2 A Simple Example

To illustrate the crazy climber algorithm, we first present in Figure 1 a simple example, namely that of the Gabor transform of a sine wave multiplied by a Gaussian envelop (top left). The square modulus of the Gabor transform with a Gaussian window is displayed at the top right of the figure, and the two integrated occupation measures are displayed in the bottom of the figure. We clearly see on the figure the different meanings of the two measures: in particular the weighted occupation measure appears as a shrunk copy of the Gabor square modulus. In addition, note that the weighting of the measure also plays the role of a thresholding for the ridge detection.

A.3 Chaining

The output of the algorithm described above is a measure on the domain D. We identify it with its density which is a function on D. The next step of the algorithm is to identify the various ridges R ℓ . This is done via a chaining procedure which replaces the occupation density by one-dimensional curves. This procedure is based on the following two steps:

1. Thresholding of the density function given as the output of the crazy climbers algorithm. The result is a series of ridges, which are graphs of curves

ϕ = ϕ ℓ (b), b = b ℓ 0 • • • b ℓ k ℓ .

IV. Numerical Results and Examples

The crazy climber algorithm was tested on several signals containing multiple ridges. We restrict the present discussion to two examples, one treated with the wavelet transform and the other with the Gabor transform.

September 21, 1998 DRAFT

A. First Numerical Results

The first signal (displayed at the top of Figure 2) is the sum of a (real) sonar signal emitted by a bat (this particular signal, displayed at the very top of Figure 2, together with noisy versions, was intensively studied in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF]) and a "linear chirp", i.e. a function of the form A(x) cos(φ(x)) with A(x) a Gaussian function and φ(x) a quadratic phase. Figure 2 shows the modulus of the wavelet transform of the signal, together with the three different ridges found by the crazy climber method: the main ridge of the bat signal, the first harmonic component, and the chirp signal. The horizontal axis is the time axis, and the vertical axis corresponds to the logarithm of the scale (i.e. the variable ϕ alluded to above). The modulus is represented with gray levels proportional to the values of the modulus of the wavelet transform. The different ridges are displayed with different gray levels.

The second signal is a speech signal, namely 250 ms of the word /one/ sampled at 8kHz. The signal is displayed at the very top of Figure 5, and the squared modulus of its Gabor transform is the third item of the figure (here we used a Gaussian window of size approximately equal to 16ms (following the lines of [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF]), and we computed the Gabor transform over the range 0Hz-4000Hz (with 100 different values for the frequency, regularly sampled, see the discussion above)). The horizontal axis is the time axis, and the vertical axis is the frequency axis (the conventions for the square modulus and ridge displays are the same as before). The crazy climbers algorithm (run with 500 climbers, and 10000 time steps each) found 18 different ridges, which are displayed at the bottom of Figure 5.

We shall come back to these examples when discussing the reconstruction from ridges, which are displayed on the same figures.

B. The Case of Noisy Signals

In many applications we can assume that we have observations f (x) of an unknown signal f 0 (x) in the presence of an additive noise ǫ(x) with mean zero. In other words we work with the model:

f (x) = f 0 (x) + ǫ(x)
and we assume that the noise is given by a mean zero stationary process with (unknown) covariance:

IE{ǫ(x)ǫ(y)} = Γ(x -y).
The case Γ = I (i.e. Γ(x -y) = δ(x -y)) corresponds to the case of an additive white noise. In some situations, "a-priori" knowledge of the noise is available. For instance it may happen that the power spectrum of the noise is known, or a piece of the signal is known to contain only noise, which gives us the chance to learn about the statistics of this noise. Then the detection algorithm may be improved by "renormalizing" the time-frequency representation, i.e. subtracting what is supposed to be the "typical" contribution of the noise.

This contribution could be chosen to be the expectation IE{M ǫ (b, ϕ)}. Moreover, if an a-priori model for the noise is available, such a quantity may be estimated by Monte-Carlo simulations, or sometimes by a direct computation.

Example: Assume that {ǫ(x)} is a second order mean zero stationary process with power spectrum of the form p(ξ) ∼ σ 2 ξ α , and that we are using with the continuous wavelet transform. In this case, IE{M ǫ (b,

a)} = IE{|T ǫ (b, a)| 2 } ∼ K α σ 2 a -α-1 provided the analyzing wavelet ψ(x) is such that K α = u α | ψ(u)| 2 du < ∞.
In most practical applications, we only have one realization of the noise component and it is impossible to compute directly this expectation. A simple ergodic argument justifies the use of the estimate:

V (ϕ) = 1 B B 0 M (b, ϕ)db. ( 13
)
Then in the penalty term used to define the ϕ motion of the climbers, the squared modulus of the time-frequency transform may be replaced with:

M (b, ϕ) = M (b, ϕ) -V (ϕ) . ( 14
)
The effect of such a modification is to avoid "trapping" the ridge in regions dominated by the noise.

As an illustration, we display in Figure 4 the ridges of the same signal as before, embedded into a Gaussian white noise, with SNR = 0dB. We can see that the main ridge is quite well reconstructed, and that the ridge of the chirp is also recovered, although only a part of it has been detected, namely the most energetic part.

The first harmonic component of the bat signal has not been detected (the corresponding wavelet transform square modulus was too low compared to the typical size of the noise, here of the form K/a for some constant K).

V. Reconstructions from the Ridges

We now address the problem of the reconstruction of a signal from the knowledge of a transform on its ridge(s). This problem was already discussed in [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF] as well as many articles on speech processing (see e.g. [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF]).

We describe here an approach based on a penalization procedure, alluded to in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF]. Even though the ridge detection part of the algorithm was independent of the time-frequency representation, this is not the case for the reconstruction. In particular, our approach is not adapted to bilinear time-frequency representations such as those given by the Wigner-Ville transform. We shall restrict ourselves to the cases of the wavelet and the Gabor transforms. Nevertheless, our approach extends to linear representations such as those obtained from matching pursuit as discussed in [START_REF] Mallat | Matching Pursuit with Time-Frequency Dictionaries[END_REF] or those developed in [START_REF] Torrésani | Time-Frequency Distributions: Wavelet Packets and Optimal Decompositions[END_REF].

A. General Discussion

In order to put our reconstruction algorithm in perspective we first review the procedures currently used.

The methods outlined in subsection V-A.3 will be developed with more details in section V-B.

September 21, 1998 DRAFT

A.1 The Transform Skeleton

The first reconstruction is the simplest one (once the ridges have been estimated). It consists in restricting the transform (whether we are working with the wavelet transform of the Gabor transform) to the ridges. It is motivated by the approximate formulas ( 6) and ( 9). More precisely, using the notation used throughout the paper, this reconstruction is given in the Gabor case by:

f (x) = 2ℜ L ℓ=1 G f (x, ω ℓ (x)) (15)
where the summation in the right hand side is restricted to the ℓ's for which ω ℓ (x) makes sense (in particular, f (x) = 0 when there is no ridge at "time" x), and by a similar formula in the wavelet case. The restriction of a transform to a ridge is sometimes called the "skeleton" of the transform [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF].

This is a very simple scheme and as Figure 3 shows, the results of this naive reconstruction can be extremely good. Its main shortcoming is that it requires the knowledge of the transform at ALL the points of the ridges . This limitation makes it impossible to subsample the ridge (for compression purposes for example).

A.2 Parametric Reconstruction

The main drawback of the above mentioned reconstruction schemes is that they are fundamentally linear, and as such, constrained by the theory of frames of wavelets and/or Gabor functions (see [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF] for a review).

This puts severe limitations on possible subsampling.

For the sake of completeness we quote a non-linear reconstruction scheme that has been successfully used in speech processing in the framework of the so-called sinusoidal model. See [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF] for a review. The main observation is that although the wavelet and the Gabor transforms restricted to ridges are often very oscillatory (and then uneasy to compress), the corresponding amplitudes A ℓ (b) and frequencies ν ℓ (b) = 1 2π φ ′ ℓ (b) are often slowly varying. This is the case for the Gabor transform of speech signals when the window is broad band.

The rationale is then to try to restrict these functions to a specific class which can be characterized by a small number of numerical parameters. The reconstructed signal is then obtained by any parametric estimation of the amplitude and the frequency, and the by integrating the estimate of the frequency to recover the phase

φ ℓ (b) = φ ℓ (b 0 ) + 1 2π b b0 ν ℓ (x)dx ,
and then computing A ℓ (b) cos(φ ℓ (b)). This may be done for all the ridges of the considered time-frequency representation, and the reconstructed signal is then the sum of all the components of the form A ℓ (b) cos(φ ℓ (b)).

A.3 Ridge Penalization

In order to reduce the amount of data necessary for the reconstruction, it is tempting to start by selecting a few sample points from each ridge, and use such points and the corresponding values of the representation to September 21, 1998 DRAFT reconstruct a signal. The reconstructed signal then takes the form of a superposition of elementary waveforms located at the sample points of the ridge. Such waveforms are in general mere perturbations of the wavelets or Gabor functions one starts with.

As we shall see, such solutions may be derived from general principles, namely by minimizing a suitably chosen quadratic functional, and using the values of the time-frequency representation at the sample points of the ridge as linear constraints.

B. The Penalization Approach

We now focus on the penalization approach, and we treat the Gabor case and the wavelet case in the same setting. Our purpose is to present a reconstruction algorithm which produces, from the mere knowledge of the considered time-frequency transform at sample points of the ridges, a very good approximation of the original signal (this reconstruction procedure was implemented and tested in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF] in the case of the wavelet transform of a signal when the latter had only one ridge).

We assume that the ridges can be parametrized by continuous functions

[b ℓ,min , b ℓ,max ] ∋ b ֒→ ϕ ℓ (b) ∈ (0, ∞)
where ℓ ∈ {1, . . . , L} is the ridge label. These ridges are usually constructed as smooth functions resulting from fitting procedures (spline smoothing is an example we are using in practical applications) from the sample points obtained from ridge estimation algorithms such as the crazy climbers algorithm presented in this paper or the snake annealing described in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF].

B.1 Statement of the Problem

We consider a (linear) time-frequency transform of an unknown signal of finite energy f 0 (x), which we denote generically by T f (b, ϕ) = f 0 , e (b,ϕ) , where the time-frequency atoms e (b,ϕ) are either wavelets or Gabor functions. We assume that the values of T f (b, ϕ) are known at sample points (b ℓ,j , ϕ ℓ,j ), with j = 1, . . . n ℓ , ℓ = 1, . . . L, which are regarded as representative of the ridges associated with the (unknown) signal f 0 (x). We use the notation z ℓ,j for the value of the transform of f 0 at the point (b ℓ,j , ϕ ℓ,j ), and e ℓ,j (x) for the corresponding function e (b ℓ,j ,ϕ ℓ,j ) (x). The set of sample points (b ℓ,j , ϕ ℓ,j ) together with the values z ℓ,j constitute what we call the skeleton of the transform of the signal to be reconstructed. Let us denote by L ℓ,j the linear form defined by

L ℓ,j f = f, e ℓ,j . (16)
As we already explained, we use smooth functions b ֒→ ϕ ℓ (b) which fit the sample points and we use the graphs of these functions as our best guesses for the ridges of the modulus of the transform of f 0 . The reconstruction problem is to find a signal f (x) of finite energy whose transform T f (b, ϕ) satisfies

L ℓ,j f = z ℓ,j , ℓ = 1, • • • L, j = 1, . . . , n ℓ (17)
September 21, 1998 DRAFT and has the union R of the graphs of the functions ϕ ℓ (b) as set of ridges. We recall that this last statement means that for each b, the points (b, ϕ ℓ (b)) of the time-frequency plane are the local maxima of the function

ϕ ֒→ |T f (b, ϕ)| 2 .
We are about to show how to construct such a signal. We will also show that it is a very good approximation of the original signal f 0 (x).

B.2 The Penalty Functions

We start with the case of the wavelet transform. Let us assume that we are given a series of L ridges of the form b ֒→ a ℓ (b), ℓ = 1, . . . L. As explained in the introduction, the reconstructed signal is obtained as the function f (x) minimizing a quadratic functional F [f ] = f, Qf with the constraints T f (b ℓ,j , a ℓ,j ) = z ℓ,j . In the case of the wavelet transform we use Q given by the integral operator defined by the kernel:

W(x, y) = δ(x -y) + ǫW 2 (x, y) . ( 18
)
The aim of the term W 2 (x, y) is to enforce smoothness of |T f | on the ridges, i.e. of the functions b → |T f (b, a ℓ (b))I. t is given by:

W 2 (x, y) = ℓ db a ℓ (b) 4 ψ x -b a ℓ (b) ψ y -b a ℓ (b) a ′ ℓ (b) 2 -ω 2 0 +ψ ′ x -b a ℓ (b) ψ ′ y -b a ℓ (b) (x -b)(y -b) a ℓ (b) 2 + 1 + x -2b + y a ℓ (b) +ψ x -b a ℓ (b) ψ ′ y -b a ℓ (b) a ′ ℓ (b) 1 + y -b a ℓ (b) +ψ ′ x -b a ℓ (b) ψ y -b a ℓ (b) a ′ ℓ (b) 1 + x -b a ℓ (b) . (19)
In the case of the Gabor transform the quadratic form Q is given by the integral kernel:

G(x, y) = δ(x -y) + ǫG 2 (x, y) , ( 20
)
where G 2 (x, y) enforces smoothness of |G f | on the ridges, i.e. of the functions b → |G f (b, ω ℓ (b))|. It is given by:

G 2 (x, y) = ℓ g ′ (x -b)g ′ (y -b) + g(x -b)g(y -b) (x -b)(y -b)ω ′ ℓ (b) 2 - (x + y -2b)ω ℓ (b)ω ′′ ℓ (b) cos (ω ℓ (b)(x -y)) db + g ′ (x -b)g(y -b) ω ℓ (b) -(y -b)ω ′ ℓ (b) - g(x -b)g ′ (y -b) ω ℓ (b) -(x -b)ω ′ ℓ (b) sin (ω ℓ (b)(x -y)) db (21)
For the sake of simplicity we have assumed that the window g(x) is real valued.

C. Solution of the Optimization Problem

The reconstruction problem has been reformulated as a problem of minimization of a sesquilinear functional, with linear constraints. This is a classical problem: We outline the solution for the sake of completeness. It may be conveniently reformulated as a minimization problem in the real domain rather than the complex domain by noticing that since we restrict ourselves to real-valued signals, both kernels are sesquilinear (G(x, y) is even real valued), and may be be replaced by their real part (which we shall still denote by the same letter). For each ℓ = 1, . . . L, the n ℓ complex constraints (17) may be replaced by the 2n ℓ real constraints as follows. Set

   ρ ℓ,j (x) = ℜe ℓ,j (x), j = 1, • • • , n ℓ , ρ ℓ,j (x) = ℑe ℓ,j (x), j = n + 1, • • • , 2n ℓ , (22)
and r j = ℜz j and r n+j = ℑz j (here, ℜ and ℑ stand for the real and imaginary parts respectively). Then the new constraints read

R j (f ) = f, ρ ℓ,j = r j , j = 1, . . . , 2n . (23)
Consequently, there exist real numbers λ ℓ,j , j = 1, . . . 2n ℓ , ℓ = 1, . . . L (the Lagrange multipliers of the problem) such that the solution f (x) of the optimization problem is given by:

f (x) = L ℓ=1 2n ℓ j=1 λ ℓ,j ρℓ,j (x) (24)
where the functions ρℓ,j (x) are defined by: ρℓ,j = Q -1 ρ ℓ,j . The Lagrange multipliers are determined by requiring that the constraints (23) be satisfied. In other words, by demanding that the wavelet transform of the function f given in (24) be equal to the z j 's at the sample points (b j , ϕ j ) of the time-scale plane. This gives a system of linear equations from which the Lagrange multipliers λ j 's can be computed. More precisely, if we denote by Λ the column vector consisting of the Lagrange multipliers, and by R the vector of values r ℓ,j , we obtain

Λ = M -1 R , (25)
M ℓ,j;ℓ ′ ,j ′ = Q -1 ρ ℓ,j , ρ ℓ ′ ,j ′ . (26)
Remark: If we denote by N = ℓ n ℓ the total number of ridge samples, we have to deal with N × N matrices, i.e. very large matrices. However, such matrices are in general sparse because of the localization properties of the functions. In particular, if the different ridges of the transform are "well separated", i.e. if they are located in different regions of the time-frequency domain, the matrix M may be approximated by a block diagonal matrix. In other words, the contributions of all ridges may be reconstructed independently. This reduces the complexity from N 3 to Ln 3 , where n is an average number of sample points per ridge. The Lagrange multipliers are then obtained as follows. If we denote by Λ ℓ the column vector consisting of the Lagrange multipliers for fixed ℓ, and by R ℓ the vector of values r ℓ,j for fixed ℓ, we obtain the following expression for the Lagrange multipliers λ ℓ,j with fixed ℓ

Λ ℓ = M -1 ℓ R ℓ , (27)
with a new n ℓ × n ℓ matrix M ℓ whose elements are given by

(M ℓ ) j,j ′ = Q -1 ρ ℓ,j , ρ ℓ,j ′ . (28)

D. The Reconstruction Algorithm

The results of the discussion of this section may be summarized in an algorithmic walk through our solution to the reconstruction problem.

• Determination of a finite set {R ℓ } ℓ=1,•••,L of ridges and, on each of them, of a set of sample points (b 1 , ϕ ℓ,1 ) =

ϕ ℓ (b 1 ), • • • , (b n , ϕ ℓ,n(ℓ) = ϕ ℓ (b n (ℓ)
) on the ridge.

• Construction of smooth estimates b ֒→ ϕ ℓ (b) of the ridges from the sample points.

• Computation of the matrix (i.e. W(x, y) or G(x, y)) of the smoothness penalty along the ridge estimate.

• Computation of the reconstruction time-frequency atoms ẽℓ,j = Q -1 e (b ℓ,j ,ϕ ℓ,j ) localized (in the timefrequency plane) at the ridge sample points.

• Computation of the coefficients λ ℓ,j .

The solution f of the reconstruction problem is then given by formula (24). Numerical examples are discussed below.

E. The Case of Noisy Signals: Smoothing Spline Type Reconstruction

Section V-B addresses the problem of the reconstruction from observations of the wavelet or Gabor transform at a finite sample of points of the time/scale or time/frequency plane. This problem was considered in [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF] when the values of the transform were assumed to be observed faithfully. In this paper we consider the possibility of an additive (possibly colored) noise in the observations of the input signal and the possibility of noise in the computation of the transform of the signal. As before our approach is motivated by the smoothing splines technique as presented in [START_REF] Wahba | Spline Models for Observational Data[END_REF]. The generalization presented in this paper was alluded to as a possible extension to the reconstruction algorithm derived and used in [START_REF] Carmona | Spline Smoothing & Extrema Representation: Variations on a Reconstruction Algorithm of Mallat and Zhong[END_REF] and [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF]. The motivation of [START_REF] Carmona | Spline Smoothing & Extrema Representation: Variations on a Reconstruction Algorithm of Mallat and Zhong[END_REF] was to simplify and shed light on the algorithm introduced by Mallat and Zhong in [START_REF] Mallat | Characterization of Signals from Multiscale Edges[END_REF] to reconstruct a signal from the extrema of its dyadic wavelet transform. The motivation of [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF] was to generalize this approach to the case of the continuous wavelet transform, the role of the extrema of the dyadic wavelet transform being played by the ridges of the continuous wavelet transform. Notice that the estimation of the ridge was taking into account the possible presence of noise while the reconstruction algorithm was assuming that the observations of the transform were correct. The reconstruction which we present now is based on a variational approach involving a penalty on the smoothness of the transform along the estimated ridges. But contrary to [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF], the observations of the transform along the ridges are not brought into the problem in the form of "knowledge" -constraints. Instead they are used to define a second penalty component. This form of the variational problem allows for a delicate balance between the fit of the transform of the solution to the observations and the smoothness of the modulus of the September 21, 1998 DRAFT transform along the ridges. Moreover the generality of the present approach makes it easy to avoid penalizing a finite dimensional space of signals which one can choose apriori.

The purpose of this subsection is to derive the formulas for the reconstruction of the original signal from the observations of the values of the transform at sample points of time/scale or time/frequency plane. As before, we use the notation T f (b, ϕ) for the transform of the signal f . This notation stands for the wavelet transform T f (b, a) as well as for the Gabor transform G f (b, ω).

We assume that we are dealing with the signal+noise model introduced earlier in Subsection IV-B and after computation of the transform of the observations, estimation of the ridges of the transform and sampling these estimates, we end up with a discrete set {(b ℓ,j , ϕ ℓ,j ), ℓ = 1, . . . L, j = 1, . . . n ℓ }, in the transform plane and observations T f (b ℓ,j , ϕ ℓ,j ) of the transform of the unknown signal at these points. We assume that the observations follow the usual linear model:

z ℓ,j = T f (b ℓ,j , ϕ ℓ,j ) + ǫ ′ ℓ,j
where the computational noise terms ǫ ′ ℓ,j are assumed to be identically distributed and uncorrelated between themselves and with the observation noise terms ǫ(x). Hence the final model is of the form:

z ℓ,j = L ℓ,j f 0 + ǫ ℓ,j , ℓ = 1, . . . L, j = 1, . . . , n ℓ (29)
where L ℓ,j is the linear form representing the value of the transform at the point (b ℓ,j , ϕ ℓ,j ) and where:

ǫ ℓ,j = T ǫ (b ℓ,j , ϕ ℓ,j ) + ǫ ′ ℓ,j .
The assumption that the two sources of noise are uncorrelated implies that the covariance matrix Σ of the ǫ ℓ,j

is the sum of the covariance of the T ǫ (b ℓ,j , ϕ ℓ,j ) and the covariance of the ǫ ′ ℓ,j . The latter being of the form σ ′2 I we have: Σ = σ ′2 I + Σ (1) where the entries of the matrix Σ (1) are given by the formula: Σ (1) ℓ,j;ℓ ′ ,j ′ k = e ℓ,j (x) Γ(x -y) e ℓ ′ ,j ′ (x) dxdy.

The reconstruction algorithm is formulated as the solution of the minimization problem:

min f 1 n Σ -1/2 (Z -T f (•, •) 2 + λ Qf, f (30)
where Z denotes the vector of observations z ℓ,j and T f (•, •) denotes the vector of values of the transform of the candidate function f at the points (b ℓ,j , ϕ ℓ,j ) and where the constant λ > 0 is introduced to balance the effects of the two components of the penalty. Theorem 1.3.1 of [START_REF] Wahba | Spline Models for Observational Data[END_REF] implies that the solution is given by: fλ

(x) = n j=1 λ ℓ,j ẽℓ,j (x) = n j=1 λ ℓ,j Q -1 e ℓ,j (x) . (31)
The coefficients λ ℓ,j in (31) are given by

Λ = nλI + Σ -1 Σ -1/2 Z (32)
where the matrix Σ is given by Σj,k = ψj , ψ k .

Remarks:

• Notice that we did not use the full generality of the smoothing spline problem as defined in [START_REF] Wahba | Spline Models for Observational Data[END_REF]. Indeed, we could have chosen a quadratic penalty of the form Q 1/2 P 1 f 2 where P 1 is the projection onto the orthogonal complement of a subspace of finite dimension. In this generality it is possible to avoid penalizing special subspaces of functions (for example, the space of polynomial functions of degree smaller than a fixed number, ....). Since the form of the solution is much more involved and since we did not find an application justifying this level of generality, we decided to use the smoothing spline approach in our simpler context.

The reader interested in this specific feature of the smoothing splines technique can consult [START_REF] Wahba | Spline Models for Observational Data[END_REF].

• The approach presented here was alluded to as a possible extension to the reconstruction algorithm derived and used in [START_REF] Carmona | Spline Smoothing & Extrema Representation: Variations on a Reconstruction Algorithm of Mallat and Zhong[END_REF] and [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF]. The latter corresponds to the case where the knowledge of the wavelet transform of the unknown signal is assumed to be perfect. In other words to the case where both Γ and σ ′2 are assumed to be zero. It is easy to see that, under these extra assumptions, the reconstruction procedure given by the above minimization problem reduces to the minimization of the quadratic form f, Qf under the constraints [START_REF] Ph | Ridge and Skeleton Extraction from the Wavelet Transform[END_REF]. This is the problem which was solved in [START_REF] Carmona | Spline Smoothing & Extrema Representation: Variations on a Reconstruction Algorithm of Mallat and Zhong[END_REF] and [START_REF] Carmona | Characterization of Signals by the Ridges of their Wavelet Transforms[END_REF]. It appears as a particular case of the more general procedure presented here. The advantages of the latter were explained in the introduction.

We shall not reproduce this discussion here.

• Notice that the reconstructed signal appears as a linear function of the observations. Nevertheless, our whole analysis is nonlinear because of the ridge estimation and the sampling of the latter.

F. An Example

Let us return to the wavelet analysis of the bat signal with the additional chirp. We used L = 3 ridges, say R 1 , R 2 and R 3 , and we chose on each ridge estimate a number of samples proportional to the length of the ridge and inversely proportional to the corresponding scale, according to the sampling theory of wavelet transforms, see [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF].

We used the value ǫ = .5 to reconstruct the signal. The result of the reconstruction is given in the second part of Figure 3. The last two plots of Figure 3 give the reconstructions of the two components: the bat signal, reconstructed from two ridges (to be compared with the top plot in Figure 2), and the chirp (the original chirp and the reconstructed one are displayed on the same plot, bottom of Figure 3). As we can see, the agreement is very good (except at the end of the chirp, where the ridge was a bit smaller than the true signal. In addition, we stress that the number of coefficients needed to characterize such a signal (i.e. twice the number of complex constraints) was approximately one fifth of the number of samples. Although compression was not our goal, the method seems to have an definite potential.

VI. Ridges and the Sinusoidal Model for Speech Signal

A popular representation of speech signals is to regard the signal as the output of a slowly time-varying filter excited by a glottal waveform. The filter models the resonant characteristics of the vocal tract. We shall not go into details of speech modeling here (we refer to [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF] and [START_REF] Maes | [END_REF] for a detailed presentation), but we notice that the resulting model for speech signal is of the form given in equation ( 5). Hence it is natural to use a time-frequency representation in order to separate the components of the signal and to express them separately. Since those components are close to harmonic components, the Gabor transform is better suited than the wavelet transform for the description of the speech signals which we want to consider. Indeed, since the wavelet processing may be viewed as a filter bank of constant relative frequency, it is not able to separate the high frequencies components (see nevertheless [START_REF] Maes | [END_REF] for a method to separate the first low frequency components).

Hence, we use the continuous Gabor transform with a Gaussian window of length approximately equal to 160 ms. Our approach is, at least in spirit, similar to that of [START_REF] Mcaulay | Speech Analysis/Synthesis Based on a Sinusoidal Representation[END_REF]. However there is a major difference: since the detection algorithm described in section III returns ridges, i.e. one-dimensional structures, the chaining method required by the Mc Aulay-Quatieri approach is not needed here.

We illustrate this discussion on the example of the /one/ signal displayed at the top of Figure 5. Our results were obtained using approximately 200 ridge samples, i.e. 400 real constraints, while the signal's length is 2048. As can be seen on the top two plots of the figure, the reconstructed signal is very close to the true one.

Of course such a comparison is not significant from the speech processing point of view. However, we stress the fact that the main features of the signal are preserved (in particular the pitch). We also tried the "ultimate test": Listening to the two sounds. They turn out to be almost indistinguishable.

VII. Conclusions

We presented a new technique to detect ridges in a surface. This algorithm is based on the stochastic relaxation of a particle system of a new type. Our detection technique performs extremely well, especially at very low SNR's. It can be used to detect ridges in all the energetic distribution representations of a signal and it is especially useful for multicomponent signals. We also presented a reconstruction procedure from the knowledge of a linear transform (such as the wavelet or the Gabor transform) on the ridges. In the case of the Gabor transform, we showed that it was performing very well on speech signals, even in the presence of significant noise disturbances.

The most important extension to the results presented in this paper would be a real time implementation.

It is relatively easy to find approximations of the reconstruction procedure which would be amenable to on line implementations. It seems more difficult to modify the ridge detection algorithm to accommodate frequent updates.

Acknowledgments: This work was done while the third named author (B.T.) was visiting the Department of Mathematics of the University of California at Irvine. Its warm hospitality should be acknowledged.

Appendix: Derivation of the Reconstruction Kernels

We give here a derivation of formulas ( 18) and [START_REF] Wahba | Spline Models for Observational Data[END_REF] and their counterparts (20) and (21) in the case of the Gabor transform. In both cases the first term aims to enforce the localization of the transform near the ridges. This is achieved by simply minimizing F 1 [f ] = ||f || 2 with the constraints [START_REF] Ph | Ridge and Skeleton Extraction from the Wavelet Transform[END_REF]. This term alone would yield a solution of the form:

f (x) = ℓ,j λ ℓ,j ψ ℓ,j (x)
where the time-frequency atoms e ℓ,j (x) = e (b ℓ,j ,a ℓ,j ) (x) are either the wavelets or the Gabor functions and where the coefficients λ ℓ,j are obtained from the z ℓ,j by multiplication with the inverse of the matrix [ e ℓ,j , e ℓ ′ ,j ′ ] ℓ,j,ℓ ′ ,j ′ .

Numerical tests show that such a solution gives accurate results if the sampling of the ridge is fine enough.

Otherwise, an extra term has to be introduced, in order to enforce the smoothness of |T f | 2 on the ridges. An adequate candidate for such a term could be given by the H 1 -norm of the restriction of the modulus of the transform to the ridges. In the case of the wavelet transform, such a term reads ℓ | d db |T f (b, a ℓ (b))|| 2 db, but unfortunately this does not define a sesquilinear form. However, if we set Ω(b, a) = arg T f (b, a), it follows from the analysis of [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF] that near the ridges number ℓ, we have d db Ω(b, a) ≈ ω0 a ℓ (b) . This suggests to approximate

An explicit computation shows that the kernel of such a form is precisely as in [START_REF] Torrésani | Time-Frequency Distributions: Wavelet Packets and Optimal Decompositions[END_REF] and [START_REF] Wahba | Spline Models for Observational Data[END_REF].

The complete derivation is left to the reader.

September 21, 1998 DRAFT persp.ps

  ) with r(b, a) ∼ O(|A ′ k |, |φ ′′ k |). As a consequence, if the wavelet ψ(x) is localized near a certain value ξ = ω 0 in the frequency domain, the wavelet transform square modulus M (b, a) = |T f (b, a)| 2 is concentrated near the N curves with equations a = a k (b) = ω 0 /φ ′ k (b). These curves of the time-frequency plane are called the ridges of the transform, and our goal is to present efficient algorithms to estimate these ridges.

  2. "Propagation" in the b direction: given a point (b, ϕ) belonging to a given ridge R ℓ , look for the "best neighbor" among (b + ǫ b , ϕ) and (b + ǫ b , ϕ ± ǫ ϕ ) (here ǫ b and ǫ ϕ are parameters fixed in advance); then iterate the process until only values below the threshold can be reached.

2 - ω 2 0

 22 d db T f (b, a ℓ (b)) ≈ d db |T f (b, a ℓ (b)| e iΩ(b,a ℓ (b) + i ω 0 a ℓ (b) |T f (b, a ℓ (b))|e iΩ(b,a ℓ (b)) ,(34)and thus use the sesquilinear functionalF [f ] = f, Qf = ||f || 2 + ǫ ℓ d db T f (b, a ℓ (b)) a ℓ (b) 2 |T f (b, a ℓ (b))| 2 db .

2 db 2 -

 22 The derivation in the Gabor case goes along the same lines. Instead of minimizing||f || 2 + ǫ ℓ d db |G f (b, ω ℓ (b)|we approximate it with the sesquilinear formF [f ] = f, Qf = ||f || 2 + ǫ ℓ d db G f (b, ω ℓ (b)) ω ℓ (b) 2 |G f (b, ω ℓ (b))| 2 db .

Fig. 1 .

 1 Fig. 1. The occupation measures for a simple windowed sine wave: left top, the signal; right top: its Gabor transform modulus; left bottom: unweighted occupation measure; right bottom: weighted occupation measure.

Fig. 2 .

 2 Fig. 2. Bat sonar signal with an additional chirp.

Fig. 3 .

 3 Fig. 3. Reconstruction from the ridges: last plot: full curve: reconstructed chirp; dashed curve: the original chirp.

Fig. 4 .

 4 Fig. 4. Bat's sonar signal and chirp embedded into a Gaussian white noise (SNR=1dB)

Fig. 5 .

 5 Fig. 5. 250 ms of the speech signal /one/ (sampling frequency: 8 kHz)

  2

September 21, 1998 DRAFT


Keywords: Polishing, Surface Roughness, 5-Axis Machine, Abrasion, Free Form Surface

Despite the large improvements in the field of high speed machining, manual polishing operations are still necessary to achieve high quality surfaces named "mirror" for mould and dies or prosthesis parts. This paper deals with the description of a new polishing method performed on a 5axis machine. This method is based on the use of an abrasive discs mounted on an flexible rotating tool describing trochoidal trajectories. Firstly, the influence of each polishing parameter on the surface roughness is highlighted on a conventional polishing machine currently used in material sciences. Based on this knowledge, a specific application concerning the polishing of a stainless medical prosthesis is detailed in order to get a mirror surface roughness which Ra is less than to 50 nm).

I. INTRODUCTION

In the bio-medical domain, polishing is a big part of the plan for procedures of medical prosthesis. Indeed, this type of process is indispensable in order to decrease surface roughness of the part. While prosthesis will be set up on a human body, all work pieces have to get a quasiperfect surface roughness for human safety. Today, the polishing operation is usually performed manually, but this working method is a long time consuming and monotone. Futhermore, in this domain, the demand increases rapidly, so bio-medical companies have to imagine new methods to polish faster than manual polishing.

Automation of polishing process on a 5-axis milling center will increase the productivity and decrease the labor intensity.

Many automatic polishing machines using bound abrasives have been investigated. Indeed, a special polishing machine tool has been developed by G. Wang et al. [START_REF] Wang | Research on polishing process of a special polishing machine tool[END_REF], named serial-parallel hybrid polishing machine tool. They proposed to combine a three parallel axis machine tool with a series of rotational mechanisms for degrees of freedom. This type of hybrid machine-tool gets advantage of the two mechanisms without its disadvantage enabling to automatize the polishing operations directly on a machine-tool.

However, to use this kind of machine, the operator must change the machine between the milling and the polishing stages. This action leads to failure in repositioning the work piece. So it would be better to keep workpiece on the same machine for limiting the duration of the total manufacturing process.

The remedy to this problem was suggested by X. Wu et al. [START_REF] Wu | New polishing technology of free form surface[END_REF]. They developed a new polishing technology of free form surface by grinding center which was studied also by Ryuh et al. [START_REF] Ryuh | An automatic tool changer and integrated software for a robotic die polishing station[END_REF]. This type of machine tool is capable of performing milling and polishing without removing the workpiece.

On this grinding center, an elastic ball tool was used. This paper shows that it is possible to make the milling and the polishing operations on the same machine tool.

C. Tournier et al. [START_REF] Tournier | Automatic polishing process of plastic injection molds on a 5-axis milling center[END_REF] have developed this idea on a 5axis milling center.

For polishing on a machine-tool, scientific research have been done on elastic polishing tools and flexible abrasive tool by J.H. Ahn et al. [START_REF] Ahn | Development of a sensor information integrated expert system for optimizing die polishing[END_REF] , S.S. Cho et al. [START_REF] Cho | Curved surface finishing with flexible abrasive tool[END_REF] and J.P. Huissoon et al. [START_REF] Huissoon | Automated polishing of die steel surfaces[END_REF] .

However, the phenomena of material removal by abrasion are not well known. Therefore, further studies seem mandatory to understand this process.

The first step that this research will examine is the evolution of the surface roughness parameter defined by ISO Standard [START_REF] Cho | Curved surface finishing with flexible abrasive tool[END_REF] (Ra; arithmetic surface roughness) as a function of the polishing time to understand what occurs during this superfinish operation.

In addition, the influences of various process parameters on the quality of the resulting surface have also been studied. This study will help to understand how abrasive machining could be carried out on a stainless steel provided by a bio-medical company.

II. EXPERIMENTAL SET UP

To understand the phenomena which govern the polishing operation, the following experiments will be carried out firstly on flat surfaces before to be applied on free form surfaces. Actually, samples will be polished on an automatic polishing machine shown in Fig. 1 using abrasive papers of different grain sizes, and different diamond pastes.

The polished material is a stainless steel (X4 Cr Ni Mn M N 21-9-4-21-10-3) especially machined for bio-medical application such as prosthesis.

The different steps of the polishing stage are the next:

Experimental method of automatic polishing on 5-axis machine

Virgile Lacharnay1 , Christophe Tournier The polishing of the samples will be stopped every 10 seconds to measure the surface roughness parameters. When the abrasive paper will be worn, it will be changed to a finer paper. This manipulation will allow plotting the evolution of the surface roughness parameters mentioned above as a function of polishing time.

The measurements are performed on a 3D surface roughness device represented on Fig. 2. 3, it is consistent that the Ra value decreases until a limit value as the time of polishing increases. This time limit represents the damage of the abrasive paper meaning that the grains have no more effect on the material removal. When this plateau is reached, the paper must be changed to a finer abrasive paper.

It could be noticed a damage of the surface roughness when the abrasive paper is changed, meaning that the following abrasive paper is going to damage the surface roughness obtained previously before to improve it in a second time. In view of this degradation, it is not necessary to go to the limit of the paper. Thus, it is possible to reduce the time of each polishing stage On the other hand, during the use of the diamond paste, there is a continuous improving of the surface roughness, showing that these pastes are well adapted to the polishing operations. From an industrial point of view, in this case the time is reduced of 30 % which is significant.

The force acting on the paper will be managed directly on the machine-tool by the deformation of the elastic rotating tool.

Influence of the rotational speed (rpm) The same experimental protocol is also used. Two rotational speeds 200 rpm and 400 rpm, with a 20 N polishing force have been tested. The curves underline that more the rotational speed is large more the limit in term of surface roughness is reached. That means the material removal rate is higher. Also, the same final Ra is obtained in both cases.

On this trial, the saving time is even more important compare to Fig. 4. The difference in time between both tests is around 110 seconds, representing a gain of 48 %.

From the industrial point of view, this parameter is easy to adjust on a CNC machine-tool.

IV. INFLUENCE OF THE PRIMARY SURFACE ROUGHNESS BEFORE POLISHING

To observe this influence, two flat samples were manufactured in parallel plans with two different radial pitch (sample 1: 0,5 mm and sample 2: 0,1 mm). The best finished surface takes 1 minute more as compared to the other.

Fig. 6 shows the tool path of the milling tool.

To observe the influence of this primary surface roughness on the evolution of the final surface roughness obtained by polishing, two samples were polished with an effort of 20 N and a 300 rpm rotation speed. On Fig. 7, the sample 1 indicates a longer polishing time to achieve the best surface roughness Ra (500 nm), highlighting that the surface quality before polishing is a significant parameter to take into account. Thus, in this experiment, the sample 1 possesses an initial arithmetical surface roughness Ra=2200 nm and the sample 2 Ra=1250 nm after milling. It is obvious that the paper will have to remove less material for the sample 2 to reach the Ra plateau (500 nm). Due to this observation the gain in time for the polishing operation is 1 min and 30 s.

Althought the milling time is 1 minute more for sample 2, due to the above comment, there is a gain of 30 s for the complete process for these one compared to sample 1.

Furthermore, the form defect is better controlled in milling than in polishing, leading to a best geometrical surface for sample 2. One of the advantages of this method is to combine the milling and polishing operations on the same 5-axis milling center.

The polishing tool path is a combination of a fractal and trochoidal trajectory as shown in Fig. 9. It describes the average trajectory of the polishing tool axis to get the most uniform polishing. For the polishing on 5-axis machine-tool, the flexible tool used and the various papers are represented in Fig. 10.

Six operations of polishing are then programmed. The first stage of experiments consists to polish with the two abrasive papers (600 and 1 200 grains/cm²) and the second stage consists to slide a felt cardboard with three diamond pastes (8, 3 and 1 µm). The last stage is the sweeping of the surface with the felt cardboard without diamond paste.

It is noted that the initial surface roughness of the milling surface is Ra=1060 nm. In consequence the abrasive paper 240 grains/cm² will not be used by referring to Fig. 3, because it is superfluous.

The different tests have been carried out according to a Taguchi design experiment. The outputs of these trials will enable to find a model of the polishing surface roughness according to the parameters mentioned in the preliminary study. Tab. 1 lists the polishing parameters and their levels considered in this study. The feed rate V f , the tilting angle θ and the displacement e are the considered parameters which have a large influence on the surface roughness and the form error of the workpiece. Each factor gets two levels with the hypothesis that there are no interactions between them.

The form error after polishing is measured with a three-dimensional measuring machine. 120 points on the workpiece have been traced to obtain the defect of the surface. The results are written in Tab. 2. On Fig. 12 it is possible to notice that the displacement (e) has the major influence on the form error, more than twice as compared to the tilting angle. Also, it is underlined that the feed rate effect is larger to this tilting angle. This design experiments allows to classify the factors according to their influence on the considered criterion, i.e the form error.

One important issue of this model is to adjust precisely the polishing parameter in order to reach a given form error in order to avoid to manufacture surface with over quality. Then a test has been performed with the best value of the parameters in order to achieve the best quality surface. This workpiece is showing a mirror surface on Fig. 13.

VI. CONCLUSIONS

This paper has highlighted that milling and polishing are possible on the same machine-tool in order to save time and to avoid the dismounting of the workpiece. Firstly preliminary tests have been carried out on automatic polishing machine to observe the relevant parameters during the polishing operation. The polishing force, the spindle frequency, the density of grains/cm² and the initial surface roughness manufactured by milling during the previous stage have been underlined. Based on this knowledge, industrial trials have been performed on a 5-axis milling machine according to a Taguchi design experiment. The results of this campaign have shown the effects of the major process parameter, highlighting the displacement of the flexible tool as the major parameter. A model to predict the form error has been found. Related to this model, the best achievable surface indicates a form error less than 10,47 µm and a surface roughness less than 15 nm (named mirror surface).

Fig. 1

 1 Fig. 1 Automatic polishing machine

Fig. 2

 2 Fig. 2 3D surface roughness device (Wyko)

Fig. 3

 3 Fig. 3 Evolution of Ra according to the time (s) Influence of the polishing force To observe the influence of the polishing pressure, two samples having the same surface roughness have been milled with the same parameters during a semi finish operation. Two different polishing force have been tested; 15 N and 30 N. Then they were polished with a 300 rpm constant rotational speed with the paper 600 grains/cm². The comparison of those results is given on Fig. 4.

Fig. 4

 4 Fig. 4 Effect of the polishing forceOn these two curves which look similar, the time to reach the plateau is quicker for the highest polishing force. In both cases the same final Ra is obtained (250 nm). The gain in time to reach the value is close 60 seconds. In

Fig. 5

 5 Fig. 5 Effect of the polishing force As shown on Fig.5 the rotational speed plays a major role for a polishing operation.The curves underline that more the rotational speed is large more the limit in term of surface roughness is reached. That means the material removal rate is higher. Also, the same final Ra is obtained in both cases.On this trial, the saving time is even more important compare to Fig.4. The difference in time between both tests is around 110 seconds, representing a gain of 48 %.From the industrial point of view, this parameter is easy to adjust on a CNC machine-tool.

Fig. 6 Fig. 7

 67 Fig. 6 Manufacturing of the primary surface roughness by milling

  Fig. 8 a) Knee prosthesis b) Free form surface to manufacture

Fig. 9

 9 Fig. 9 Fractal and trochoidal tool path

Fig. 10

 10 Fig. 10 Description of the polishing tool

(2 3 )

 3 According to those results a model indicating the effect of each parameter is found according to the Taguchi theory. Y i represents the form defect. a 1 (-20), a 2 (20), b 1 (-29), b 2 ,(29) c 1 (-13), c 2 (13), are the effects of the polishing parameters.l is the average value of the form error (72,90 µm).

Fig. 12

 12 Fig. 12 Effect of the polishing parameter on the form error

Fig. 13

 13 Fig. 13 Surface obtained for trial N°1: Form error = 10,47 µm Surface roughness Ra=15 nm

TABLE 2 TAGUCHI

 2 DESIGN EXPERIMENT L 4

LURPA/Department of Mechanical Engineering, ENS, 61 Ave du Président Wilson, 94

CACHAN, France 2 LaBoMaP/Department of Machining, Arts et Métiers ParisTech, 71 250 Cluny, France (gerard.poulachon@ensam.eu)