Adrien Meynard
email: adrien.meynard@ens-lyon.fr

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

Synthesis-based time-scale transforms for non-stationary signals

Keywords: Non-stationarity, Time-frequency and time-scale synthesis, Time warping, Maximum a posteriori estimation, Expectation-Maximization, Sharp time-scale representation

Synthesis-based time-scale transforms for non-stationary signals

Adrien Meynard, Bruno

Introduction

In the statistical signal processing or time series analysis literature, the spectral analysis problem aims at estimating the spectral content of a stationary signal or time series, called power spectrum (in other words, the distribution of power as a function of frequency) from a discrete, finite-length sample. The existence of such a power spectrum is guaranteed by the Wiener-Khinchin theorem. In addition, under suitable assumptions, random signals may also be represented in terms of filtering of white noise, where the filter's frequency response is defined by the square root of the power spectrum. This representation is known as Cramér's representation, we refer to [START_REF] Koopmans | The spectral analysis of time series[END_REF] for a thorough presentation. In this context, spectral analysis may be addressed by parametric models (for example AR or ARMA), or nonparametric approaches (see e.g. [START_REF] Koopmans | The spectral analysis of time series[END_REF][START_REF] Stoica | Introduction to Spectral Analysis[END_REF]). In the latter case, one generally resorts to Fourier-based techniques such as periodograms, involving windowing and time averaging.

However, many signals are not stationary, and these tools do not apply anymore. The extension of spectral analysis to non-stationary situations has been considered by several authors (see e.g. [START_REF] Priestley | Non-linear and non-stationary time series analysis[END_REF][START_REF] Dahlhaus | A likelihood approximation for locally stationary processes[END_REF]), who proposed generalizations adapted to specific classes of non-stationary signals.

Among possible approaches, time-frequency and time-scale analysis [START_REF] Flandrin | Explorations in Time-Frequency Analysis[END_REF][START_REF] Boashash | Time-Frequency Signal Analysis and Processing: A Comprehensive Reference[END_REF][START_REF] Carmona | Practical Time-Frequency Analysis[END_REF]] has received considerable attention during the last decades, as it allows one to step away from traditional stationary signal models expressed as linear combinations of sine waves. Linear time-frequency or time-scale transforms have been widely used in applications, to provide visualizations from which non-stationary features can be estimated. They also benefit from a strong, well established mathematical support (see e.g. [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF]). A main aspect of linear transform is the existence, under suitable assumptions, of inverse mappings, which express signals as linear combinations of time-frequency localized waveforms. These transforms are known to be limited by uncertainty principles: even though the classical Heisenberg inequality does not explicitly involve a joint time-frequency domain, the short-time Fourier transform (to quote an example) can be viewed as the output of a filter bank, and the impulse responses of the filters are affected by Heisenberg's inequality. Super-resolution techniques such as reassignment or synchrosqueezing [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF] appear as mere post-processing of one or several linear transforms affected by uncertainty, which attempt to correct the loss of timefrequency resolution by suitable non-linear modifications. An important asset of synchrosqueezing is the preservation of invertibility: reconstruction is possible from synchrosqueezed wavelet or short time Fourier transforms.

Amplitude and Frequency modulated (AM-FM) signal models have received considerable attention as well, often in different contexts. Originally defined in terms of the Hilbert transform (see [START_REF] Picinbono | Some remarks on instantaneous amplitude and frequency of signals[END_REF]), instantaneous amplitude and frequency can also be estimated using such tools. The problem was later on addressed through time-frequency transforms [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF][START_REF] Carmona | Practical Time-Frequency Analysis[END_REF] and super-resolution techniques [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF]. More generally, models involving sums of AM-FM signals have also been thoroughly studied in various contexts, including speech modeling [START_REF] Mcaulay | Speech analysis/synthesis based on a sinusoidal representation[END_REF] and music processing [START_REF] Depalle | Tracking of partials for additive sound synthesis using hidden Markov models[END_REF], and more recently for physiological signal analysis [START_REF] Wu | Instantaneous frequency and wave shape functions (I)[END_REF] and [START_REF] Chui | Signal decomposition and analysis via extraction of frequencies[END_REF]. These models express signals as linear combinations of amplitude and frequency modulated "modes", to be identified. Again there is a vast literature on the subject which we will not review here, let us just mention empirical mode decomposition methods [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], iterative filtering [START_REF] Cicone | Adaptive local iterative filtering for signal decomposition and instantaneous[END_REF][START_REF] Cicone | Iterative filtering as a direct method for the decomposition of nonstationary signals[END_REF] or the frequency extraction approach of [START_REF] Chui | Signal decomposition and analysis via extraction of frequencies[END_REF], which have enjoyed significant interest recently.

However, not every non-stationary signal can be described as AM-FM or a sum of AM-FM components. Many physically relevant signals indeed involve FM-like non-stationarity, but the latter is combined with a time varying spectrum broader than standard AM-FM signals. An example can be found in To the best of our knowledge, the time warping paradigm has been first studied in depth by Clerc and Mallat in [START_REF] Clerc | Estimating deformations of stationary processes[END_REF], who developed estimators based upon small scale limite of wavelet transform. It has been revisited more recently in [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes[END_REF][START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF], in a context of Gaussian random signals, thus avoiding the small scale limit. In [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF], an approximate maximum likelihood approach called JEFAS for the joint estimation of time warping and power spectrum has been developed and studied. JEFAS, which assumes slow variations of the amplitude modulation and time warping function, is essentially a post-processing of a time-scale (wavelet) transform.

In this paper, we develop and study an alternative approach, that may be called a synthesis approach, inspired by a model developed by Turner & Sahani [START_REF] Turner | Time-frequency analysis as probabilistic inference[END_REF], which interprets time-frequency transform in terms of probabilistic inference. The Turner-Sahani was applied to various time-frequency problems. Our model expresses non-stationary signals as the result of a time-frequency or time-scale synthesis from a random time-frequency representation, governed by a prescribed prior distribution, which involves the time warping model. We study the corresponding model, and compute the posterior distribution which is used to construct an estimator for the model parameters. This approach, called JEFAS-S (where the "S" stands for "synthesis"), also outputs an adapted time-scale representation from which signal synthesis is possible. Remarkably enough, JEFAS-S is able to handle time warping models in situations where the time warping function can have much faster variations. In addition, JEFAS-S can also yield time-frequency representations which are much sharper than conventional transforms, and compare well to synchrosqueezed or reassigned transforms, while allowing signal synthesis.

Some of the results presented in this paper were already sketched in [START_REF] Meynard | Time-scale synthesis for locally stationary signals[END_REF]; we refer also to [START_REF] Meynard | Stationnarités brisées : approches à l'analyse et à la synthèse[END_REF], where the AM-TW model and variants were analyzed in details.

This paper is organized as follows. Background material is given in Section 2, including the amplitude modulation and time warping (AM-TW) model, elements of time-frequency and time-scale analysis, and an account of the JE-FAS approach. The main contributions are described in Section 3, which in particular develops the statistical model, the estimation procedure and corresponding algorithm. Numerical results and discussion are given in Section 4, and Section 5 is devoted to conclusions. Some more technical aspects are developed in the Appendix.

Models and time-frequency tools

The time warping model

Amplitude and frequency modulations are standard tools in signal processing. However, amplitude and frequency modulated (AM-FM) signals

y(t) = a(t) cos(γ(t))
may also be understood as amplitude modulated and time warped (AM-TW) signals, time warping being defined as a time-varying clock change. More precisely

y = M α D γ x , (1)
where M α is the amplitude modulation operator (pointwise multiplication with the positive valued function α), and D γ is the time warping operator (composition with the monotone function γ, assumed to be differentiable and strictly increasing), and a = α √ γ ′ :

M α x(t) = α(t)x(t) , D γ x(t) = γ ′ (t)x(γ(t)) . (2)
Locally harmonic signals of the form y(t) = α(t) k cos(kγ(t)) may also be described likewise. We stress that time warping is often more physically relevant than frequency modulation, examples are provided by car engine sounds during acceleration, Doppler effect (which is generated by speed difference), or ECG signals during heart rate variations.

The AM-TW model was extended to the setting of random signals in [START_REF] Clerc | Estimating deformations of stationary processes[END_REF], in the form

Y = M α D γ X , (3)
here X is a second order wide sense stationary random process, with zero mean, and absolutely continuous spectral measure. The corresponding power spectrum will be denoted by S X . In analogy with the spectral estimation problem for a stationary random processes X (estimate the power spectrum S X from one or several realizations of the process), spectral analysis in this context may be described as the problem of joint estimation of α, γ and the underlying power spectrum S X . The problem has been addressed by various authors (see e.g. [START_REF] Priestley | Non-linear and non-stationary time series analysis[END_REF][START_REF] Dahlhaus | A likelihood approximation for locally stationary processes[END_REF]), for different signal models. Since our focus is on AM-TW models, we won't discuss these generalizations further, and limit to time-frequency and time-scale transforms which provide well suited representations for AM-TW models.

Time-frequency and wavelets

Let us first consider the continuous time setting. For the sake of simplicity we limit ourselves to time-scale representation, i.e. wavelet analysis and synthesis, similar developments can be made for more general linear time-frequency transforms. We denote by ψ ∈ L 1 (R) ∩ L 2 (R) the analysis wavelet, and by ψ s the scaled wavelets defined by ψ s (t) = q -s/2 ψ q -s t , s ∈ R , for some reference scale q > 1. The corresponding wavelet transform associates with every x ∈ L 2 (R) the sequence of band-pass filtered signals W x (s, •) defined by the convolution product

W x (s, •) = x * ψ s ,
where ψ s (t) = ψ s (-t). In the continuous time setting, it is well known that such a transform is invertible provided ψ satisfies a suitable admissibility condition [START_REF] Grossmann | Decomposition of Hardy functions into square integrable wavelets of constant[END_REF]. If the admissibility condition is fulfilled, and after suitable normalization of the wavelet ψ, this yields the inversion formula

x = ∞ -∞ ψ s * W x (s, •) ds s . (4)
As a consequence of the admissibility condition, the wavelet ψ must be zeromean, i.e. an oscillatory function. Classical choices include the Morlet wavelet (which has to be modified to fulfill the admissibility condition), derivatives of a Gaussian function and many others. It is actually often relevant to consider socalled analytic wavelets, i.e. complex wavelets whose Fourier transform vanishes for negative frequencies. Among these, popular choices include the so-called Morse wavelets [START_REF] Olhede | Generalized Morse wavelets[END_REF], defined by their Fourier transform ψ(ν) = Kν k e -ν γ for ν > 0, the time-frequency localization of which is controlled by the parameters k ∈ N * and γ ∈ R * + (K is a normalization constant) and the Altes wavelet [START_REF] Altes | Some invariance properties of the wide-band ambiguity function[END_REF]. We will focus here on the so-called sharp wavelet, used in [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF]

ψ(ν) = exp ln(ǫ) δ(ν, ν 0 ) δ(ν 1 , ν 0 ) , ν ∈ R + . (5)
Here ν 0 is the mode of ψ and ν 1 is a cutoff frequency, defined so that ψ(ν 1 ) = ǫ for some prescribed tolerance ǫ < 1. The divergence δ is defined by δ(a, b) =

1 2 a b + b a -1.
The quality factor of ψ, i.e. the center frequency to bandwidth ratio, can be expressed as Q = 1/ C(C + 4), where C = -δ(ν 1 , ν 0 ) ln(2)/ ln(ǫ) > 0.

Using such wavelets, the wavelet transform W x of real valued signals is complex valued. The synthesis formula (4) has to be modified accordingly, to read (for a suitably chosen normalization of ψ)

x = Re ∞ -∞ ψ s * W x (s, •) ds s . (6)
Since we limit to time-scale representation, we follow the standard logarithmic discretization scheme for the scale variable q s , equivalently linear discretization scheme for s. Time is sampled on a regular, scale independent lattice. The discrete wavelet transform used throughout this paper is therefore very close to the so-called stationary wavelet transform (except that scales are not necessarily integer powers of 2).

JEFAS: an analysis-based approach

The current paper is mainly devoted to a synthesis-based approach, however it also relies on an earlier approach called JEFAS (Joint Estimation of Frequency, Amplitude and Spectrum) which we outline here for completeness.

Background

Let Y denote a signal generated by time warping from some unknown signal X, as in [START_REF] Priestley | Non-linear and non-stationary time series analysis[END_REF]. X is assumed to be zero mean and wide sense stationary, with power spectral density S X distributed following a normal law

The approach relies on approximations of the following form [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF]. Assume polynomial decay for the wavelet:

|ψ(t)| ≤ (1 + |t| β ) -1 < ∞ for some β > 2.
Assume further that the power spectrum S X of the underlying stationary signal X is such that ∞ 0 ν 2ρ S X (ν) dν < ∞ with ρ = (β -1)/(β + 2). Then we have

W Y (s, τ ) = α(τ )W X (s + log q (γ ′ (τ )), γ(τ )) + R(s, τ ) = W Y (s, τ ) + R(s, τ ) . (7)
Here W Y is the approximate transform, and R is a remainder. In a nutshell, amplitude modulation and time warping in the time domain amount, up to some accuracy, to amplitude modulation and translation in the time-scale domain.

For completeness, we sketch here the proof of the result, which is given in details in [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF] (a deterministic version can be derived similarly following the lines of [START_REF] Clerc | Estimating deformations of stationary processes[END_REF]). The idea is that whenever the functions α and γ are smooth enough, they can be Taylor-expanded around a time value of interest: given a test function g located near t = τ (i.e. decaying fast enough as a function of

|t -τ |), Taylor expansions near t = τ yield M α g(t) ≈ M α(τ ) g(t) , D γ g(t) ≈ D τ γ g(t) , with D τ γ ∆ = T τ D -log q (γ ′ (τ )) T -γ(τ )
T τ denoting translation by τ . Therefore, we approximate the wavelet transform of Y by

W Y (s, τ ) ≈ M a(τ ) D τ γ X, T τ D s ψ = α(τ ) D -log q (γ ′ (τ )) T -γ(τ ) X, D s ψ
which yields the desired expression since the adjoint operators of D and T are such that D * s = D -s and T * τ = T -τ . Theorem 1 in [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF] provides upper bounds for the remainder term R.

Provided the remainder R in ( 7) can be neglected, W Y may be approximated by W Y , which is a circular (complex) Gaussian random field (see [START_REF] Picinbono | Second-order complex random vectors and normal distributions[END_REF] for details on circular normal distributions, the main aspects are sketched in Appendix A). The distribution of W Y is characterized by its covariance kernel C Y , which may be written in terms of the covariance C X of X, α and γ. This suggests to estimate the unknown parameters using maximum likelihood (ML) approach. Unfortunately, ML is intractable in this case, as it leads to huge nonlinear optimization problems. The JEFAS algorithm of [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF] proposes an ansatz that reduces significantly the computational load: fixed time slices of the wavelet transform W Y are considered independent.

Estimation

Suppose one is given a finite dimensional wavelet transform W ∈ C Ms×Nτ , and denote by W = vec(W) ∈ C MsNτ the vectorized version. W is assumed to result from the sampling of the continuous transform W Y of a time warped signal Y as in (3), on a time-scale grid Λ = M s × N t . The goal is to estimate the set of parameters Θ = (θ 1 , θ 2 , θ 3 ), where

θ 1 , θ 2 , θ 3 ∈ R Nτ are defined by θ 1 (τ n ) = α(τ n ) 2 , θ 2 (τ n ) = log q (γ ′ (τ n )) and θ 3 (τ n ) = γ(τ n ).
Here τ n denotes the n-th time sample, for n = 1 . . . N τ . Parameters corresponding to time τ n are denoted by Θ n .

As alluded to above, consider fixed time slices w n , n = 1 . . . N τ of W as independent random vectors leads to approximating the covariance matrix C W (Θ) by a block diagonal matrix, each block C n (Θ n ) corresponding to a fixed time index n. The negative log likelihood B(Θ) = -L (Θ) takes the form of a Bregman log det divergence [START_REF] Cherian | Jensen-Bregman logdet divergence with application to efficient similarity search for covariance matrices[END_REF] between the covariance matrices WW H and C(Θ) (up to an additive constant that depends on W). As a consequence of the independence assumption, B(Θ) can be approximated by a sum of reduced (i.e. fixed time) terms B r (Θ n ), which can be optimized independently

B(Θ) = Tr WW H C W (Θ) -1 -ln det C W (Θ) -1 (8)
≈ Nτ n=1 B r (Θ n ) , where (9)
B r (Θ n ) = Tr w n w H n C n (Θ n ) -1 -ln det C n (Θ n ) -1 (10)
Notice that the covariance matrix C n (Θ n ) also depends on the unknown power spectrum S X of the underlying stationary signal X, its matrix elements are given as follows: for each pair of scales i, j = 1, . . . M s ,

(C n (Θ n )) ij = θ 1 (τ n ) f (s i + θ 2 (τ n ), s j + θ 2 (τ n )) , ( 11
)
where f is a Hermitian symmetric, positive semi-definite function defined by

f (s, s ′ ) = q (s+s ′ )/2 ∞ 0 S X (ξ) ψ (q s ξ) ψ q s ′ ξ dξ . ( 12
)
The function f is called covariance template of the model. Notice also that the covariance matrix C n (Θ n ) does not depend explicitly on parameters θ 3 (τ n ) anymore. This is a consequence of the independence assumption made on fixed-time slices of the time-scale transform.

Therefore, estimating the parameters requires the knowledge of the power spectrum S X . Given current values of the parameters, an estimate W X of the wavelet transform W X of the underlying stationary process can be obtained by interpolation using [START_REF] Carmona | Practical Time-Frequency Analysis[END_REF], from which an estimate of the power spectrum can be computed by time averaging | W X | 2 (see Chapter 6 in [START_REF] Carmona | Practical Time-Frequency Analysis[END_REF]).

The JEFAS algorithm proceeds by iterative updates of the parameter Θ and power spectrum S X . At convergence, estimates for the time warping γ and amplitude modulation α are finally obtained by interpolation from the parameters Θ. Details of the algorithm can be found in Algorithm 1 in [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF], which also provides statistical performance assessments.

Discussion

We now discuss some pros and cons of JEFAS.

We believe JEFAS provides a relevant approach to extend the spectral estimation problem to non-stationary situations that are well accounted for the the AM-TW (amplitude modulation combined with time warping) model. While an algorithm similar to JEFAS can also be developed for AM-FM models, time warping is probably a more relevant model than frequency modulation in a number of situations of interest.

JEFAS is rooted on an approximate maximum likelihood approach, and therefore benefits from tools of the statistical machinery. This includes among others precision estimates (see [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF][START_REF] Meynard | Stationnarités brisées : approches à l'analyse et à la synthèse[END_REF]), which are useful to assess quantitatively the performances of the approach.

However, JEFAS is mainly a post-processing of a time-frequency transform (here a wavelet transform), and thus inherits from the weaknesses of the latter. These include time-frequency resolution limited by classical time-frequency uncertainty, and a difficulty to handle fast variations of the time-dependent spectrum. In addition, JEFAS is based upon simplifying assumptions, in particular assumptions of slow variations of the time warping function. This increases the difficulty of JEFAS to handle fast non-stationarities.

These weaknesses motivate the development of an alternative approach, that avoids a prior computation of time-frequency transform and explicit slow variations assumptions. This approach is described in the next section.

JEFAS-S: a time warping based signal synthesis model

We now turn to the main contribution of this paper, namely the synthesisbased approach. Instead of starting from a transform, the approach is based on a synthesis model for a non-stationary random signal Y , in the form

Y = Re s∈S ψ s * W s + ε , ( 13
)
where ε is a zero mean, Gaussian white noise with variance σ 2 . For each scale s, W s is a random process, with prescribed prior probability distribution, which will implement the time warping model. Based on the latter, the goal is to estimate W s from a realization y of Y for all scales s under consideration.

Problem formulation in finite dimensional setting

As before, we denote by (τ 1 , . . . τ Nτ ) the vector of times at which the signal is sampled. Signals of interest and wavelets ψ s are then N τ -dimensional vectors, for which we further assume periodic boundary conditions:

x(τ n+Nτ ) = x(τ n ).
Let y ∈ R Nτ be the signal under consideration. We will seek estimates for W s at time samples. We denote by S = {1, 2, . . . M s } the set of considered scales. For each n = 1, . . . N τ , let w n ∈ R Ms denote the vector of wavelet coefficients at time τ n , and let W = w 1 . . . w Nτ be the resulting wavelet transform matrix. For every n = 1, . . . N τ , introduce the matrix

Ψ n =      ψ s1 (τ 1-n ) ψ s2 (τ 1-n ) . . . ψ sM s (τ 1-n ) ψ s1 (τ 2-n ) ψ s2 (τ 2-n ) . . . ψ sM s (τ 2-n ) . . . . . . . . . . . . ψ s1 (τ Nτ -n ) ψ s2 (τ Nτ -n ) . . . ψ sM s (τ Nτ -n)      ∈ C Nτ ×Ms , (14)
then equation ( 13) reads

y = Re Nτ n=1 Ψ n w n + ε = Re (DW) + ε , (15)
where we have introduced the matrices D ∈ C Nτ ×Nτ Ms and W ∈ C Nτ Ms defined by

D = Ψ 1 Ψ 2 . . . Ψ Nτ , W = vec(W) =      w 1 w 2 . . . w Nτ      . (16)

Random model for noise and time-scale coefficients

As mentioned above, the noise ε in ( 13) is modelled as a real Gaussian white noise, with zero mean and variance σ 2 , therefore ε ∼ N (0, σ 2 I Nτ ).

We now introduce the time-scale warping prior model on W. This prior model is designed so as to mimic the behavior of wavelet transform of time warped stationary processes, and involve the following elements 1. Columns w n are assumed independently distributed, according to a circular complex Gaussian law (this was an approximation in JEFAS, and is now a prior model in JEFAS-S). We refer to e.g. [START_REF] Picinbono | Second-order complex random vectors and normal distributions[END_REF] for an extended description of multivariate complex Gaussian distributions, a short account in given in Appendix A. 2. The covariance matrix C n of a given column w n takes a form similar to equations [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF], for some Hermitian symmetric, positive semi-definite covariance template f . The choice for f depends on the application and the a priori shape of the underlying stationary signal. Two covariance templates are investigated in this paper. (a) A natural choice for f is given by [START_REF] Mcaulay | Speech analysis/synthesis based on a sinusoidal representation[END_REF]. The only difference is that since amplitude modulation is not considered here, the parameter θ 1 (τ n ) is not present any more. This prior brings a prescribed regularity of representation along the scale axis. This ensures a good robustness of the representation to errors in estimation of the translation parameter θ 2 (τ n ).

(b) When the underlying stationary signal is a sum of K sine waves, a priori constraints on the sparse structure of the power spectrum can be made. To do this, we simply modify the expression of the covariance template f by replacing S X in [START_REF] Mcaulay | Speech analysis/synthesis based on a sinusoidal representation[END_REF] by

S X,sp (ξ) ∆ = K k=1 a 2 k δ(ξ -ξ k ) , (17)
where a k and ξ k are the amplitudes and frequencies of the sine waves forming the stationary signal. This choice provides sharper time-scale representations.

Remark 1 (Comparison with JEFAS).

1. Unlike JEFAS, amplitude modulation is not included in JEFAS-S. This is motivated by the fact that in JEFAS, amplitude modulation estimation was shown to be poorer than time warping. Nevertheless, there is no theoretical obstruction to the inclusion of amplitude modulation in the a priori model if necessary. The estimation of amplitude modulation is briefly addressed in Section 3.5.5 below. 2. The property of independence of the columns of the time-scale transform is a simplification in JEFAS, which is actually quite questionable, in particular for large scales. In JEFAS-S, this simplification is implemented as a property of the prior probability distributions, which is less constraining: choosing independent prior means that we don't want to assume that columns are a priori informative about each other, and leave it to the observations to tune these dependences. As we shall see, the posterior distribution of columns does not have this independence property any more.

Summarizing, the model parameters are the noise variance σ 2 , the warping parameters θ 2 (τ n ) and the covariance kernel f .

Bayesian estimation 3.3.1. Posterior expectation

Let us first assume that the model parameters are known, we want to estimate the vector W, using the posterior expectation. Since we are in the Gaussian case, the posterior expectation coincides with the maximum a posteriori estimator.

The prior distribution of W is the circular complex Gaussian

p(W) = 1 π Nτ Ms det(Γ pr ) exp - 1 2 W H Γ -1 pr W (18)
(see Appendix A for details) where the prior covariance matrix Γ pr is the block diagonal matrix with blocks C n :

Γ pr = bdiag (C 1 , . . . C Nτ ) ∈ C Nτ Ms×Nτ Ms (19)
As mentioned above, the noise ε ∼ N (0, σ 2 I Nτ ) is white Gaussian and real valued. Therefore the likelihood of observations is given by

p(Y|W) ∼ N (Re(DW), σ 2 I Nτ ) (20)
When both noise and prior are real Gaussian, it is known that the posterior distribution is also real Gaussian. We provide here a similar result with circular complex Gaussian prior, which shows that the posterior distribution is complex Gaussian, non circular. Proposition 1. Let W ∈ C Nτ ×Ms be a circular complex Gaussian random vector, distributed following [START_REF] Cicone | Iterative filtering as a direct method for the decomposition of nonstationary signals[END_REF]. Let y ∈ R Nτ be a random vector whose likelihood p(y|W) is as in [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes[END_REF]. Then the posterior distribution p(W|y) is a complex Gaussian law p(W|y)

∼ CN (µ po , Γ po , R po ) (21)
with

µ po = Γ pr D H C -1 y y (22)
Γ po = Γ pr - 1 4 Γ pr D H C -1 y DΓ pr (23)
R po = -D H D + 4σ 2 Γ -1 pr -1 D T D Γ po (24)
and

C y = σ 2 I + 1 2 Re DΓ pr D H ( 25
)
The proof is given in Appendix B. The expression of the relation matrix R po is provided for the sake of completeness, although it will not be used in the sequel.

Let W = E W|y {W} denote the conditional expectation of W. We have that W = µ po , the posterior mean. The columns w n of the time-frequency representation W can be written as

w n = 1 2 C n Ψ H n C -1 y y , (26)
where the block diagonal structure of Γ pr yields the following form for C y

C y = σ 2 I Nτ + 1 2 Re Nτ n=1 Ψ n C n Ψ H n . (27)
Notice that C y is a matrix of dimension N τ × N τ , which can be large for real world signals, yielding heavy computational load for evaluation, storage and inversion.

Remark 2. As stressed in Remark 1, vectors w n are modeled as independent vectors in the prior model. However, estimated vectors w n , which take into account the actual signal, are not uncorrelated anymore:

E w n w H n ′ = δ nn ′ C n - 1 4 C n Ψ H n C -1 y Ψ n ′ C n ′ ,
which has no reason to vanish for n = n ′ .

Parameter estimation: JEFAS-S

We now turn to the parameter estimation problem. In this section we assume that the covariance template f and the noise variance σ 2 are fixed. The choice of f will be discussed later on, and the noise variance σ 2 will play the role of regularization parameter.

Let us recall the parameters of the AM-TW model, that were introduced in Section 2.3.

2: Θ = (θ 1 , θ 2 , θ 3 ), with θ 1 (τ n ) = α(τ n ) 2 , θ 2 (τ n ) = log q (γ ′ (τ n ))
and θ 3 (τ n ) = γ(τ n ). JEFAS-S only considers the TW model (α = 1) and, as in JEFAS, the time decorrelation in the prior model implies that only θ 2 enters into the estimation procedure.

For the sake of simplicity, and to avoid confusions with the JEFAS parameters, let us denote here by ϑ n = θ 2 (τ n ), and ϑ = (ϑ 1 , . . . ϑ Nτ ) ∈ R Nτ the remaining parameters to be estimated. With these notations, the covariance matrix associated to time τ n takes the form

C n = C(ϑ n ) , with (C(ϑ)) ij = f (s i + ϑ, s j + ϑ) , ∀ϑ . ( 28
)
Here f is a covariance template, as introduced in ( 12), which will control the behavior of the time-frequency representation W.

ϑ can be estimated using an Expectation-Maximization (EM) strategy [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. EM is an iterative algorithm that can be used to search for the maximum likelihood estimators of model parameters. In general, EM is not guaranteed to converge to a global maximum of the likelihood. However, the likelihood is guaranteed to increase at each iteration.

After some initialization, parameters ϑ are updated iteratively. Denoting by ϑ (k-1) the estimate from iteration k -1, a new estimate ϑ (k) is computed as follows • Expectation step: introduce the auxiliary function

Q ϑ, ϑ (k-1) = E (k-1) {ln (p ϑ (y, W))} (29)
where p ϑ (y, W) is the joint probability density function of W and y, which depends on ϑ, and where we have denoted by E (k-1) the expectation with respect to W conditional to y, using parameter value ϑ (k-1) .

• Maximization step: solve

ϑ (k) = arg max ϑ Q ϑ, ϑ (k-1) . ( 30
)
The application of the EM paradigm leads to the following update rules [START_REF] Meynard | Stationnarités brisées : approches à l'analyse et à la synthèse[END_REF].

For the sake of simplicity let us introduce some notations. After iteration k -1, denote by W (k-1) and w the n-th diagonal block of the estimate of the posterior covariance matrix Γ po from iteration k -1 (which is a function of ϑ (k-1) ).

Proposition 2 (EM update rules). With the above notations, the update rules for the time scale representation and parameters read 1. Time scale representation

w (k) n = 1 2 C (k-1) n Ψ H n C (k-1) y -1 y , (31)
2. Parameters:

ϑ (k) = arg min ϑ Nτ n=1 ln det(C(ϑ n )) + w (k-1) H n C(ϑ n ) -1 w (k-1) n + Tr C(ϑ n ) -1 Re Γ (k-1) n ( 32
)
The proof is given in Appendix C.

Remark 3. The parameter update can also be written as

ϑ (k) = arg min ϑ∈R Nτ Nτ n=1 B (k-1) n (ϑ n ) ,
where for all ϑ (33)

B (k-1) n (ϑ) = Tr w (k) n w (k)H n + Γ (k-1) n C(ϑ) -1 -ln det C(ϑ) -1 . ( 34
)
Notice the similarity with equations ( 8)- [START_REF] Picinbono | Some remarks on instantaneous amplitude and frequency of signals[END_REF], involving now Bregman type log det divergences between covariance matrices C(ϑ) and w

(k) n w (k)H n + Γ (k-1) n .

Signal synthesis

After convergence of the estimation algorithm, denote by w n , C n and C y the corresponding estimated quantities. According to the JEFAS-S model [START_REF] Chui | Signal decomposition and analysis via extraction of frequencies[END_REF], the denoised signal given by

y 0 = Re Nτ n=1 Ψ n w n ( 35
)
can be estimated as

y 0 = Re Nτ n=1 Ψ n w n = 1 2 Re Nτ n=1 Ψ n C n Ψ H n C y -1 y . ( 36
)
The quality of the estimation can be assessed by the bias and covariance of the estimator, which are given in the following result. Proposition 3. Let y 0 denote the above estimator of the denoised signal.

1. The estimator y 0 is biased, with bias

b = E { y 0 -y 0 } = -σ 2 C -1 y y 0 . ( 37
)
2. The covariance matrix of the estimation error e = y 0 -y 0 is given by

R e = σ 2 I -σ 2 C -1 y 2 . ( 38
)
Algorithm 1 ( W, ϑ, S X ) = JEFAS-S(y, σ 2 , ̺, N ′ τ , f ) Require: Initial estimates of θ(0) and S (0) X using JEFAS.

1: k ← 1.

2: while stopping criterion (39) = FALSE do

3:
for all n ∈ {1, . . . , N τ } do 4:

Restrict Ψ n , C y ϑ (k-1)
and y to the interval [n-N ′ τ /2, n+N ′ τ /2].

5:

Evaluate w

(k)
n using the corresponding reduced version of (31).

6:

end for 7:

Estimate ϑ (k) by solving [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

8:

Estimate S (k)

X using the wavelet-based estimate.

9:

Update matrices C n and Γ pr using ϑ

(k) , S (k)
X and f . 10:

k ← k + 1. 11: end while Remark 4. Let us denote by γ 2 1 , γ 2 2 , . . . γ 2
Nτ the eigenvalues of the Hermitian, positive semidefinite matrix 1 2 Re(DΓ pr D H ), sorted in increasing order. 1. We notice that the bias may be bounded as

b ≤ 1 + γ 2 1 σ 2 -1 y 0 ,
(γ 2 0 is the smallest eigenvalue). The ratio γ 2 1 /σ 2 may be interpreted as a signal to noise ratio related to the prior model. The larger the signal to noise ratio, the smaller the bias. 2. The variance Tr(R e ) of the estimation error may also be upper bounded, using the submultiplicativity of the trace (i.e. Tr(AB) ≤ Tr(A)Tr(B)), so that

Tr(R e ) ≤ N 2 τ σ 2 1 - 1 N τ Nτ n=1 1 + γ 2 n σ 2 -1 2
.

Algorithm

The algorithm whose pseudo-code is given in Algorithm 1 implements the approach described above. The algorithm is iterative and therefore requires an initialization and a stopping criterion, which are described below. The input includes the input signal y, the noise variance σ 2 and a parameter ̺ controlling the stopping criterion. The algorithm also involves a dimension reduction step, which speeds up some matrix operations. This step depends on a parameter N ′ τ ≪ N τ , which is discussed below.

Initialization

The algorithm requires an initial guess S (0)

X for the power spectrum of the underlying stationary signal and for the parameters ϑ (0) . These can be obtained by running JEFAS on the input signal y.

Remark 5. In general, EM is not guaranteed to converge to a global maximum of the likelihood. However, the likelihood is guaranteed to increase at each iteration. When multiple local minima are expected, the initialization plays a crucial role.

In our experiments, the initialization described above turned out to provide satisfactory results. When JEFAS fails to converge, an alternative initialization can be obtained by computing at each time a local scale defined as a weighted average of the scale, the weights being given by the time-scale transform.

Stopping criterion

EM ensures that the likelihood monotonically increases during iterations. The algorithm stops at iteration k when the condition

L (ϑ (k) ) -L (ϑ (k-1) ) < ̺ ( 39
)
is true. Here ̺ > 0 is a parameter fixed by the user.

Dimension reduction

The algorithm requires solving at each iteration a high dimensional linear system in [START_REF] Oberlin | The second-order wavelet synchrosqueezing transform[END_REF]. However, the large matrix C y turns out to be strongly concentrated around its diagonal. This can be exploited to speed up the evaluation of w n in ( 26) by restricting y on a neighborhood

[n -N ′ τ /2, n + N ′ τ /2] of the time index n with a certain bandwidth N ′ τ ≪ N τ .
The latter has to be matched to the effective support size of wavelets at the largest scale upon consideration.

Power spectrum estimation

JEFAS-S requires the knowledge of the power spectrum S X of the underlying stationary signal, both for initialization and during the iterations.

Concerning the initialization, an estimate can be obtained using the JEFAS approach, if the latter is successful. Otherwise, a very crude estimate can be computed using a Welch periodogram from the input signal y.

During the iterations of Algorithm 1, the power spectrum is re-estimated by using the two-step procedure 1. Correcting the time-scale transform W (k-1) for scale translations ϑ

(k-1) n 2.
Estimating the wavelet spectrum [START_REF] Carmona | Practical Time-Frequency Analysis[END_REF], i.e. the time-average of the corrected time-scale representation.

Amplitude modulation estimation

If amplitude modulation is taken into account, the a priori covariance template takes the form [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF]. The whole procedure described in the above remains true. It is only necessary to add a step for the update of the estimate of the amplitude modulation parameter between lines 7 and 8 of Algorithm 1.

Solving problem [START_REF] Li | Adaptive short-time Fourier transform and synchrosqueezing transform for non-stationary signal separa[END_REF] with respect to the amplitude modulation parameter yields a closed-form expression (see eq. ( 17) in [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF]). Therefore, considering the full AM-TW model instead of the TW model alone is computationally inexpensive. As mentioned in Remark 1, amplitude modulation is not considered here due to the poor performance of the AM estimator, which can be assessed in terms of the variance of the estimator. In [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF] a closed-form expression was given for the corresponding Cramér-Rao lower bound, which was shown to be quite large in general, impairing the quality of estimates.

Numerical results

The numerical results presented in the current section are reproducible. The MATLAB code is available online at https://github.com/AdMeynard/JEFAS.

Illustration on wideband non-stationary signal

We implement JEFAS-S on an AM-TW synthesized signal. The underlying stationary signal is formed by a stationary Gaussian random vector of length N τ = 2048. The power spectrum was chosen as the output of JEFAS applied to the audio wind recording whose scalogram is displayed in Fig. 1. The time warping function was also derived from the output of JEFAS on the real signal. Amplitude modulation was not taken into consideration. Eventually, a Gaussian white noise of variance σ 2 = 0.04 was added to the observations. The stationary and TW signals, and their respective wavelet transforms are shown in Fig. 2.

As reported in [START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF], slow varying AM-TW signals are satisfactorily processed using JEFAS. Indeed, time warping is well described in the wavelet domain as a time-dependent scale shift of the coefficients. It is therefore estimable via JE-FAS. The time warping function is accurately estimated after 6 iterations: the mean square error is 0.00155. The initialization of JEFAS-S with the JEFAS output results in a fast convergence in 6 iterations. The mean square error is 0.00150. It is incrementally improved at the cost of a significant increase in computing time: 33 s with JEFAS against 157 s with JEFAS-S (experiments performed on a laptop equipped with an Intel Core i9 CPU). However, in addition to the estimates of the parameters characterizing the AM-TW model, JEFAS-S provides a time scale representation and a synthesis operator for the noise-free signal. Both quantities are displayed in Fig. 3. Visual inspection shows that the time-scale representation provided by JEFAS-S displays more clearly the spectral content coming from the AM-TW signal than the wavelet transform. For example, the frequency bands located around 200 Hz and 600 Hz appear clearly, which was not the case in Fig 2 . Moreover, the synthesis operator makes it possible to recover the meaningful part of the signal by filtering out the background noise.

To precisely quantify the ability of JEFAS-S to achieve denoising on AM-TW signals, we apply the synthesis operator to this AM-WP signal with various noise levels (σ 2 ranging from 0.0005 to 1). The input signal-to-noise ratio (SNR) is defined as

SNR I = 10 log 10 S X (ξ)dξ σ 2 . ( 40
)
The reconstruction formula from the time-scale representation, expressed in (35), provides a supposedly denoised signal. Hence, we measure the output SNR given by

SNR O = 10 log 10 E{ y 0 2 } E { ỹ0 -y 0 2 } . (41)
In this example, the quality of denoising only depends on the performance of the synthesis estimator determined by its bias (37) and its variance (38). Indeed, we have E ỹ0 -ỹ0 2 = Tr(R e ) + b T b. The output SNR is plotted against the input SNR in Fig. 4. Numerical results show that for input SNRs smaller than 19 dB, JEFAS properly performs denoising. When the input SNR is higher than 19 dB, the intrinsic limitations of the Bayesian estimator become preponderant. The synthesized signal is thus degraded. However, the smaller the input SNR, the better the denoising. The a priori knowledge of the non-stationarity class to which the signal belongs makes the denoising operation adaptive, which can be seen as a form of non-stationary filtering.

The subsequent illustration shows that JEFAS-S is also able to handle AM-TW signals on which JEFAS fails.

Illustration on a narrowband non-stationary signal

We illustrate here the ability of JEFAS-S to handle non-stationary TW signals with fast varying spectral content.

To this end, a non-stationary signal is synthesized from a stationary signal x comprised of two sine waves with close frequencies. We have

x[n] = cos 2πf 1 n N τ + ϕ 1 + 2 cos 2πf 2 n N τ + ϕ 2 , (42)
where N τ = 1024, f 1 = 184 and f 2 = 154. The phases ϕ 1 and ϕ 2 are independent random scalars drawn from a uniform distribution in the interval [0, 2π).

We then apply the time warping operator. The warping function γ is derived from the electrocardiogram of a patient with irregular heartbeat. Its derivative is obtained evaluating the R-peak to R-peak intervals, followed by a cubic spline interpolation. This fast varying function is shown in Fig. 6. Eventually, a Gaussian white noise of variance σ 2 = 2.5 × 10 -3 is added to the observations. The stationary and TW signals, and their respective wavelet transforms are displayed in Fig. 5.

JEFAS and JEFAS-S

The wavelet transform of the TW signal is highly disturbed by interference between the two components. Indeed, on the one hand, the ability to discriminate two close sine waves requires the analysis wavelet to be sharply localized in frequency, but then interference appears along the time axis. On the other hand, the fast variations of the warping function are only discernible if the wavelet is sharp in time, but interference then appears along the scale axis. The impossibility to avoid interference is a consequence of the uncertainty principle. Hence, the approximated behavior of the wavelet transform JEFAS is based on is not valid. On this example JEFAS turns out to converge in 3 iterations, but the algorithm does not converge satisfactorily towards the ground truth warping function. This is shown in the middle of Fig. 6. The rapid variations of γ ′ are strongly smoothed by the JEFAS estimate. By setting ̺ = 0.5, JEFAS-S converges in 111 iterations. The estimated γ ′ function is represented in the bottom of Fig. 6. The rapid variations are now accurately reproduced. The accuracy of the estimate is markedly improved: the mean square error goes from 0.0261 for JEFAS to 0.0069 for JEFAS-S. The timescale representation estimated with JEFAS-S is depicted in the top of Fig. 7. Although the presence of two components is not visible, we clearly visualize the rapid variations of the time warping function. Unlike the wavelet transform, interference is absent here.

Sharpened time-scale representation

Once the parameters of the model have been estimated, we can take advantage of the expression [START_REF] Olhede | Generalized Morse wavelets[END_REF] to adapt the form of the time-scale representation to new a priori constraints. To do this, we simply change the expression of the covariance template f to its sharp form provided by [START_REF] Cicone | Adaptive local iterative filtering for signal decomposition and instantaneous[END_REF]. Its implementation requires the estimates of the number of sine waves forming the stationary signal, their amplitudes and their frequencies. These estimates are obtained using the local maxima of the spectrum estimated by JEFAS-S. In this example, we find two local maxima, i.e K = 2. By construction, the wavelet-based spectrum estimator is smooth. The estimation of K is therefore robust to noise. The estimated frequencies are ξ1 = 155 and ξ2 = 184. The resulting time-scale transform is shown in the middle Fig. 7, and as expected turns out to be much sharper than the time-scale representation given on the top image. In addition, the two components forming the signal appear quite distinctly, the synthesis was indeed able to separate these two components.

Besides, the comparison with the first-order synchrosqueezed wavelet transform, depicted in the third panel of Fig. 7, shows that our sharpened representa- tion presents a much better concentration of the time-scale coefficients along the instantaneous frequencies, and does not suffer from the smoothing effect which is visible on the synchrosqueezed scalogram. The second-order synchrosqueezed wavelet transform [START_REF] Oberlin | The second-order wavelet synchrosqueezing transform[END_REF], displayed in the fourth panel of Fig. 7, accounts for the rapid variations of the instantaneous frequency to allow for better localization in the time-frequency plane, but does not dissociate the two frequency components either. In addition, it tends to oversmooth the time warping function in comparison with JEFAS-S. State-of-the-art synchrosqueezed transforms are able to provide sharp representations for signals containing multiple components with unnecessarily related instantaneous frequencies (see e.g. [START_REF] Li | Adaptive short-time Fourier transform and synchrosqueezing transform for non-stationary signal separa[END_REF]). In particular, methods that handle crossover instantaneous frequency curves are now available [START_REF] Chui | Time-scale-chirp rate operator for recovery of non-stationary signal components with crossover instanta[END_REF]. This source separation task is not tractable with our method, as the model assumes simultaneous evolution of the instantaneous frequencies. Our approach, as a by-product of JEFAS-S, is helpful when a single non-stationarity source is present.

Conclusions

We have developed in this paper a new approach for time-scale representation of a class of non-stationary signals, which is able to handle fast varying nonstationarities. In comparison with existing approaches, a main originality is that JEFAS-S is based upon explicit probabilistic modeling (in the spirit of [START_REF] Turner | Time-frequency analysis as probabilistic inference[END_REF]) of a time-scale synthesis representation, which allows to bypass (to a certain extent) constraints due to uncertainty principles.

JEFAS-S can handle both locally broadband and narrowband signals. In the latter case, i.e. for signals generated by time warping of sums of sine waves, JEFAS-S is also able to yield extremely sharp time-scale representations, sharper and much less biased than reassigned or synchroszqueezed spectrograms. Since JEFAS-S is a synthesis based approach, we also point out that the signal can be synthesized from the sharpened representation (as well as a stationary signal constructed from a corrected time-scale transform).

The approach is limited so far to non-stationarities generated by Time Warping (TW), but could be extended to other types, for example by adding amplitude modulation, or introducing frequency modulation or more general frequency nonlinear transformation (possibly using a more general time-frequency representation within which the nonlinearity could be described by a translation on a nonlinearly modified frequency axis). 4. The distribution of a complex Gaussian random vector with mean µ, covariance matrix Γ and relation matrix R is denoted by CN (µ, Γ, R), and takes the form

p(Z) = K exp - 1 2 Z -µ Z -µ T Γ R R Γ -1 Z -µ Z -µ , (A.3)
for some constant K. The inverse of the above block covariance-relation matrix takes the form

Γ R R Γ -1 = A B B A , with    A = Γ -RΓ -1 R -1 B = -ARΓ -1
. (A.4)

5. The complex Gaussian random vector W is called circularly symmetric, or circular, if for every (deterministic) ϕ ∈ (-π, π], e iϕ W and W have the same distribution. In such a case, µ = 0 and R = 0, one then writes W ∼ CN c (0, Γ). the latter expression following from Woodbury's formula, which states that for every conformable matrices A, U, C, V , one has (A + U CV ) -1 = A -1 -A -1 U (C -1 + V A -1 U ) -1 V A -1 . This proves [START_REF] Turner | Time-frequency analysis as probabilistic inference[END_REF].

We now turn to the expression of the posterior covariance matrix. Equating the expressions in (A.4) with the expressions obtained above yields

Γ po -R po Γ -1 po R po = 4σ 2 D H D + 4σ 2 Γ -1 pr -1 R po = -D H D + 4σ 2 Γ -1 pr -1 D T D Γ po .
Inserting the expressions of R po Γ

Figure 1 :

 1 Figure 1: Scalogram of a 10-second-long wind recording.

Figure 1 ,

 1 Figure 1, which represents the so-called scalogram (squared modulus of wavelet transform) of a 10-second-long wind recording. Besides the larger local bandwidth already mentioned, the time variations of the local spectrum are quite often poorly accounted by amplitude and frequency modulation, and time warping (i.e. time dependent clock change) is often more physically realistic. For example, the sound of a accelerating car engine is suitably described by time dependent clock change, as well as stationary signals deformed by Doppler effect. One may also argue that time warping is a convenient model to describe non-stationary ECG signals, which originate from heart beats whose rhythm vary as a function of time.To the best of our knowledge, the time warping paradigm has been first studied in depth by Clerc and Mallat in[START_REF] Clerc | Estimating deformations of stationary processes[END_REF], who developed estimators based upon small scale limite of wavelet transform. It has been revisited more recently in[START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes[END_REF][START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF], in a context of Gaussian random signals, thus avoiding the small scale limit. In[START_REF] Meynard | Spectral analysis for nonstationary audio[END_REF], an approximate maximum likelihood approach called JEFAS for the joint estimation of time warping and power spectrum has been developed and studied. JEFAS, which assumes slow variations of the amplitude modulation and time warping function, is essentially a post-processing of a time-scale (wavelet) transform.In this paper, we develop and study an alternative approach, that may be called a synthesis approach, inspired by a model developed by Turner & Sahani[START_REF] Turner | Time-frequency analysis as probabilistic inference[END_REF], which interprets time-frequency transform in terms of probabilistic inference. The Turner-Sahani was applied to various time-frequency problems. Our model expresses non-stationary signals as the result of a time-frequency or time-scale synthesis from a random time-frequency representation, governed by a prescribed prior distribution, which involves the time warping model. We study the corresponding model, and compute the posterior distribution which is used to construct an estimator for the model parameters. This approach, called JEFAS-S (where the "S" stands for "synthesis"), also outputs an adapted time-scale representation from which signal synthesis is possible. Remarkably

  for the time-scale representation and its columns, and by C (k-1) n the current estimate of the n-th block of C y . Let also C (k-1) y denote the corresponding estimate of C y (thus depending on ϑ (k-1) ). Similarly, we denote by Γ (k-1) n

Figure 2 :

 2 Figure 2: Synthetic stationary signal (upper left) and corresponding TW signal (upper right), derived from an audio recording of the wind. Their scalograms are below.

Figure 3 :

 3 Figure 3: JEFAS-S time-scale representation (top), and a segment of the derived synthesized signal (bottom, red) compared to the observations (black) and the noise-free signal (blue).

Figure 4 :

 4 Figure 4: Output SNR given by the synthesis estimator, plotted against the input SNR. The dotted line marks the separation between denoising (upper part) and signal degradation (lower part).

Figure 5 :

 5 Figure 5: Synthetic stationary signal (upper left) and corresponding TW signal (upper right), derived from heartbeat measurements. Their scalograms are below.

Figure 6 :

 6 Figure 6: Derivative of the warping function, derived from real heartbeat measurements (top). The resulting JEFAS (middle) and JEFAS-S (bottom) estimates are superimposed.

Figure 7 :

 7 Figure 7: JEFAS-S time-scale representation (first panel), and its sharpened version (second panel). The first (third panel) and second-order (fourth panel) synchrosqueezed wavelet transforms are below for comparison.

3 .

 3 The mean and covariance matrix of the associated real random vector Z are given by + R) Im(-Γ + R)Im(Γ + R) Re(Γ -R) . (A.2)

2 Γ 1 = σ 2 I + 1 2 Re

 212 Setting z = Re(Dµ po ), then 2σ 2 z = Re DΓ pr D H (y -z) which yields z = 2σ 2 I + Re(DΓ pr D H ) -1 Re DΓ pr D H y and therefore µ po = 1 pr D H C -1 y y , where C y = σ 2 I -2σ 2 I + Re(DΓ pr D H ) -1 Re DΓ pr D H -DΓ pr D H

- 1 po 1 . 1 =Γ 1 pr + D T D - 1 D

 11111 and R -1 po into the first of the two above equations yieldsΓ po = 4σ 2 I + D H D -DObserve that by Woodbury's formula we have thatI -D D H D + 4σ 2 Γ I -D D T D -D H D + 4σ 2 Γ so that, setting C = (4σ 2 ) -1 I + D D H D + 4σ 2 Γ po = Γ -1 pr + D H CD -1 = Γ pr -Γ pr D H C -1 + DΓ pr D H -1 DΓ pr . Now C -1 = 4σ 2 I -D -D H D -4σ 2 Γ -T = 4σ 2 I + D Γ pr D T ,therefore we obtainΓ po = Γ pr -Γ pr D H 4σ 2 I + DΓ pr D H + D Γ pr D T -1 DΓ prwhich yields[START_REF] Meynard | Time-scale synthesis for locally stationary signals[END_REF].

Appendix A. Complex Gaussian distributions

We recall here basic facts on multivariate complex Gaussian distributions, referring to [START_REF] Picinbono | Second-order complex random vectors and normal distributions[END_REF] for more details.

1. A complex random vector W ∈ C N follows a complex Gaussian distribution if the real random vector

follows a real Gaussian distribution N (µ Z , Γ Z ). 2. A complex Gaussian distribution is characterized by its mean µ ∈ R N , and its covariance and relation matrices

) the superscript H denoting Hermitian conjugation and the superscript T denoting transposition. Γ is then Hermitian symmetric, and R is symmetric.

Appendix B. Proof of Proposition 1

Being the composition of complex Gaussian distributions, the posterior distribution of W given y is complex Gaussian, non-symmetric, and its probability density function takes the form

Therefore, using the notations of (A.3) for the inverse covariance-relation matrix, with µ = µ po , Γ = Γ po and R = R po , we obtain

Using ( 18) and ( 20) together with Bayes' formula yields

By identification with the previous expression, we obtain the parameters

After some algebra, we obtain that µ po is a solution of the equation

Appendix C. Proof of Proposition 2

We first notice that according to [START_REF] Olhede | Generalized Morse wavelets[END_REF], the estimate for the time-scale representation given ϑ (k-1) reads

which is the desired expression. EM proceeds by alternating the expectation and maximization steps, which are described below.

Appendix C.1. Expectation

The goal is to evaluate the quantity defined in [START_REF] Cherian | Jensen-Bregman logdet divergence with application to efficient similarity search for covariance matrices[END_REF], which may be written as

According to [START_REF] Omer | Time-frequency and time-scale analysis of deformed stationary processes[END_REF], the first term does not depend on ϑ and can be ignored in the optimization. The second term reads

Focusing on the last term, setting X = Re(W), Y = Im(W) and Z = (X T , Y T ) T , we have

Now, since at iteration k and conditionally to y, Z|y ∼ CN m (k-1) , G (k-1) , with

and where we have set for simplicity

.

Therefore we have that

po .

Putting things together we finally obtain which is the desired expression.


Issraa Shahine

Jean-Yves Mevellec

Mireille Richard-Plouet

Bernard Humbert

Pierre-Yves Tessier

Nanoporous gold stacked layers as substrates for SERS detection in liquid or gases with ultra-low detection limit and long term stability

Surface enhanced Raman scattering (SERS) substrates constituted of stacked ultrathin nanoporous gold layers are used to detect very low concentrations of molecules in liquids or gases. The SERS substrates are obtained by copper chemical etching of alternative copper and gold stacked nanolayers. This process is a reliable method for fabricating uniform and reproducible SERS substrates, with a robust SERS response at extremely low detection limits.

By optimizing fabrication conditions combined with a thorough analysis of the SERS mappings and using 2,2-bipyridine (BP) as probe molecules, the detection threshold corresponding to a detectable SERS response reaches a BP concentration of 10 -18 mol.L -1 in water. An additional Raman mechanism is also highlighted by µ-Surface Enhancement Spatially Offset Raman spectroscopy (µ-SESORS): gold ligaments, inside nanoporous layers, act as waveguides for the incident light, leading to a significant increase in the size of the active SERS area. Moreover, these SERS substrates could be stored during several days without significant decrease of their properties, and they can be reactivated and reusable after SERS analysis. The ability to detect low BP vapor pressure is also demonstrated.

INTRODUCTION

Surface enhanced Raman scattering (SERS) is widely used for the detection of low concentration analytes. This sensing technique covers a large field of applications in chemistry, [START_REF] Schlücker | Surface-Enhanced Raman Spectroscopy: Concepts and Chemical Applications[END_REF][START_REF] Fan | A Review on Recent Advances in the Applications of Surface-Enhanced Raman Scattering in Analytical Chemistry[END_REF][START_REF] Langer | Present and Future of Surface Enhanced Raman Scattering[END_REF][START_REF] Pérez-Jiménez | Surface-Enhanced Raman Spectroscopy: Benefits, Trade-Offs and Future Developments[END_REF] bioanalysis 5- 7 and biomedical spectroscopy. [START_REF] Kneipp | Surface-Enhanced Raman Scattering and Biophysics[END_REF][START_REF] Qian | Single-Molecule and Single-Nanoparticle SERS: From Fundamental Mechanisms to Biomedical Applications[END_REF][START_REF] Tanwar | Advancing Raman Spectroscopy from Research to Clinic: Translational Potential and Challenges[END_REF][START_REF] Nicolson | Spatially Offset Raman Spectroscopy for Biomedical Applications[END_REF] SERS has also been studied for its ability to the detection of proteins, viruses, and microorganisms in immunoassays. [START_REF] Porter | SERS as a Bioassay Platform: Fundamentals, Design, and Applications[END_REF][START_REF] Bruzas | Advances in Surface-Enhanced Raman Spectroscopy (SERS) Substrates for Lipid and Protein Characterization: Sensing and Beyond[END_REF][START_REF] Zheng | Label-Free SERS in Biological and Biomedical Applications: Recent Progress, Current Challenges and Opportunities[END_REF][START_REF] Cui | Diagnostic Methods and Potential Portable Biosensors for Coronavirus Disease 2019[END_REF] SERS spectroscopy is based on the dramatic increase of the Raman scattering signals by the molecules adsorbed on the surfaces of noble metal nanostructures, that exhibit surface plasmon resonance when exposed to a light of an appropriate wavelength. Due to the strong near-field amplification properties of the plasmonic nanostructures combined with a direct chemical interaction between the adsorbed chemical species on the metal surface, a huge electromagnetic enhancement of the local field is produced at both the excitation wavelength and the Raman wavelength. And thus the Raman signal of the molecules to be detected is significantly amplified. [START_REF] Langer | Present and Future of Surface Enhanced Raman Scattering[END_REF][START_REF] Pérez-Jiménez | Surface-Enhanced Raman Spectroscopy: Benefits, Trade-Offs and Future Developments[END_REF][START_REF] Kumar | Hot Spots in Ag Core-Au Shell Nanoparticles Potent for Surface-Enhanced Raman Scattering Studies of Biomolecules[END_REF][START_REF] Lal | Nano-Optics from Sensing to Waveguiding[END_REF][START_REF] Willets | Localized Surface Plasmon Resonance Spectroscopy and Sensing[END_REF][START_REF] Stiles | Surface-Enhanced Raman Spectroscopy[END_REF] Moreover, localized surface plasmon resonance can be produced in the interstices between metal nanostructures, or in nanogaps between aggregated nanoparticles called "hot spots", which contribute greatly to the creation of a strong electromagnetic field. [START_REF] Pérez-Jiménez | Surface-Enhanced Raman Spectroscopy: Benefits, Trade-Offs and Future Developments[END_REF][START_REF] Zong | Surface-Enhanced Raman Spectroscopy for Bioanalysis: Reliability and Challenges[END_REF] A significant enhancement of the Raman signal is then generated when the molecules to be detected are located in these hot spots. Therefore, hot spots play a dominant role in the detection of ultra-low molecule concentration or even single molecule with an enhancement factor of the emitted Raman signal that could reach 10 10 or more. [START_REF] Pérez-Jiménez | Surface-Enhanced Raman Spectroscopy: Benefits, Trade-Offs and Future Developments[END_REF][START_REF] Kneipp | Surface-Enhanced Raman Scattering and Biophysics[END_REF][START_REF] Qian | Single-Molecule and Single-Nanoparticle SERS: From Fundamental Mechanisms to Biomedical Applications[END_REF] Considerable efforts have been devoted towards the fabrication of efficient SERS substrates, aiming to detect very low concentrations of molecules. A wide range of materials could be used as SERS substrates based on aggregation of nanoparticles or on metal nanostructures fabricated directly on solid substrates. [START_REF] Fan | A Review on Recent Advances in the Applications of Surface-Enhanced Raman Scattering in Analytical Chemistry[END_REF][START_REF] Langer | Present and Future of Surface Enhanced Raman Scattering[END_REF][START_REF] Porter | SERS as a Bioassay Platform: Fundamentals, Design, and Applications[END_REF][START_REF] Le Ru | Surface Enhanced Raman Scattering Enhancement Factors: A Comprehensive Study[END_REF][START_REF] Mosier-Boss | Review of SERS Substrates for Chemical Sensing[END_REF][START_REF] Rao | Nanoporous Gold Nanoparticles and Au/Al2O3 Hybrid Nanoparticles with Large Tunability of Plasmonic Properties[END_REF][START_REF] Camelio | Linear Chains of Ag Nanoparticles Embedded in Dielectric Films for SERS Applications in Analytical Chemistry[END_REF][START_REF] Jonker | A Wafer-Scale Fabrication Method for Three-Dimensional Plasmonic Hollow Nanopillars[END_REF][START_REF] Wu | DNA Structure-Stabilized Liquid-Liquid Self-Assembled Ordered Au Nanoparticle Interface for Sensitive Detection of MiRNA 155[END_REF][START_REF] Xu | SERS Approach to Probe the Adsorption Process of Trace Volatile Benzaldehyde on Layered Double Hydroxide Material[END_REF] Particularly, metal nanoporous films have attracted much attention due to their remarkable morphology characterized by three-dimensional networks of interconnected ligaments at the nanoscale. [START_REF] Zhang | Effect of Residual Silver on Surface-Enhanced Raman Scattering of Dealloyed Nanoporous Gold[END_REF][START_REF] Qian | Surface Enhanced Raman Scattering of Nanoporous Gold: Smaller Pore Sizes Stronger Enhancements[END_REF][START_REF] Chen | Nanoporous Copper with Tunable Nanoporosity for SERS Applications[END_REF][START_REF] Lang | Characteristic Length and Temperature Dependence of Surface Enhanced Raman Scattering of Nanoporous Gold[END_REF][START_REF] Liu | Single Molecule Detection from a Large-Scale SERS-Active Au79Ag21[END_REF][START_REF] Chao | Gold-Rich Ligament Nanostructure by Dealloying Au-Based Metallic Glass Ribbon for Surface-Enhanced Raman Scattering[END_REF][START_REF] Gan | Nanoporous Metals Processed by Dealloying and Their Applications[END_REF][START_REF] Ruffino | Nanoporous Gold-Based Sensing[END_REF][START_REF] Zalm | Review-Recent Advances in the Development of Nanoporous Au for Sensing Applications[END_REF][START_REF] Cialone | Nanoporous FePd Alloy as Multifunctional Ferromagnetic SERS-Active Substrate[END_REF][START_REF] Xiao | Nanoporous Gold: A Review and Potentials in Biotechnological and Biomedical Applications[END_REF] Their porosity gives them a high specific surface area to adsorb molecules and improve their performances as sensors. SERS sensitivity of these films could be attributed to the presence of hot spots located in nanogaps between ligaments. A keypoint is to find new synthesis strategies to change the nature of plasmonic couplings, in particular between ligaments, in order to increase the local electric fields and/or the density of hot spots in these nanoporous films.

Recently, we have shown that it is possible to easily synthesize new types of gold nanoporous lamellar films with a strong anisotropy of their morphology along the perpendicular axis to the substrate surface. [START_REF] Chauvin | Lamellar Nanoporous Gold Thin Films with Tunable Porosity for Ultrasensitive SERS Detection in Liquid and Gas Phase[END_REF] The synthesis method is based on the copper etching by nitric acid in stacked nanolayers of gold and copper. During the process, copper, which is the less noble metal, is dissolved to create porosity between gold nanolayers. At the same time, there is a diffusion of gold, the noble metal, at the electrolyte/solid interfaces which leads to the creation of ligaments.

The final morphology of nanoporous films could be considered as stacked layers of ultrathin nanoporous gold films, and this very specific anisotropic morphology is due to the specific initial Au/Cu lamellar design of the film before its etching in nitric acid. In this method, the initial thickness of the gold nanolayers before copper etching plays a key role in the final morphology of the gold nanoporous film. This morphology directly impacts the SERS behavior of the gold nanoporous film. It is attributed to the critical distance between neighboring gold ligaments according to the morphology which is one of the main parameters for the creation of "hot spots". This first preliminary study on this type of substrates revealed sensitivities to bipyridine (BP) molecules of the order of 10 -12 mol.L -1 . [START_REF] Chauvin | Lamellar Nanoporous Gold Thin Films with Tunable Porosity for Ultrasensitive SERS Detection in Liquid and Gas Phase[END_REF] The objective of this present paper is to show that a very detailed statistical analysis based on SERS mappings of the studied Au nanoporous film leads to the possibility to reach ultimate molecules detection threshold much lower than 10 -12 mol.L -1 . Another goal is to demonstrate that this process is a reliable method for fabricating uniform, highly reproducible and efficient reusable SERS substrates, with a robust SERS response at extremely low detection limits. Moreover, µ-Surface Enhancement Spatially Offset Raman Spectroscopy (µ-SESORS) characterization is carried out to study the possibility that nanoporous layers act as waveguides for the light leading to an increase of the actual SERS tested area.

EXPERIMENTAL METHODS

Synthesis of nanoporous Au thin films and BP adsorption.

The thin film deposition is carried out on monocrystalline (100) silicon substrates using a gold target (purity, 99.99 %) and a copper target (purity, 99.99 %). Before the Cu/Au multilayers deposition, a chromium film (120 nm thick) is primarily deposited followed by a 50 nm thick protecting film of pure Au. Before film deposition, the base pressure is lower than 4×10 -5 Pa. The Ar plasma pressure is regulated at 0.5 Pa, and the electrical power applied to the Cu and Au targets is fixed to 150 W and 30 W respectively. Concentrated nitric acid (65%, 16 mol.L -1 Carlo Erba Reagents) is used for the copper etching of the stacked multilayers at room temperature for 300 min. After etching, the samples are directly dipped into deionised (DI) water for 2 min in order to stop the etching reaction. The immersion in DI water is repeated for additional 5 min. The second washing process ensures the removal of the residual nitric acid possibly remaining in the samples. After etching, except for the storage study, the washed samples are immediately immersed in 10 ml batches of DI water solutions of 2,2-BP with concentrations varying from 10 -18 mol.L -1 up to 10 -4 mol.L -1 . The incubation of the gold substrates in BP solutions lasts for 24 hours. After this adsorption process, the samples are dried and immediately tested by Raman measurement.

For the detection of the BP gas, the etched, washed in DI water and dried SERS substrate is characterized by Raman measurements to establish a blank reference. Then, solid grains from BP powder are put at about 1 cm apart from the sample being under SERS microscope to follow the SERS spectrum evolution with time.

Characterization of nanoporous Au thin films. Scanning electron microscopy (SEM)

imaging is carried out on a JEOL JSM 7600 F microscope at 5 kV. The XPS measurements are performed using an AXIS Nova from Kratos having a monochromatic Al Kα X-ray source (1486.6 eV). The spectra are recorded with a pass energy of 20 eV and are decomposed using CasaXPS. [START_REF] Biesinger | Advanced Analysis of Copper X-Ray Photoelectron Spectra[END_REF] The spectra are calibrated by setting adventitious carbon C 1s component at 284.7 eV. Extinction spectra are acquired using a Perkin lambda 1050 equipped with a 150 mm integrating sphere. Raman measurements are performed with a Renishaw InVia Raman Microscope using 633 nm excitation laser line. The spectral resolution is 3.5 cm -1 given by an 1800 grooves/mm grating. The setup is equipped with an automated XYZ table, where the displacement motors are able to generate 100 nm steps with a spatial precision verified on an AFM grid. The SERS spectra are recorded with an acquisition exposure time of one second and a laser beam of 50 µW focused by a 50x Olympus objective (numerical aperture of 0.9) on around 4 µm 2 . SERS mapping spectra have been collected in an area of 60 × 60 μm 2 under the same conditions of the single measurements with a spatial gap of 5 μm between each spectrum.

Therefore, 169 Raman spectra (13 x 13) are obtained for each mapping.

In order to obtain the different spectral sources constituting the experimental spectra and to follow their evolutions with the bipyridine concentration, we adopt a curve resolution method using all the spectral data recorded from 10 -12 mol.L -1 to 10 -4 mol.L -1 , corresponding to a total number of 845 spectra. We use the Bayesian Positive Source Separation (BPSS) approach that was developed [START_REF] Carteret | Polymorphism Studied by Lattice Phonon Raman Spectroscopy and Statistical Mixture Analysis Method. Application to Calcium Carbonate Polymorphs during Batch Crystallization[END_REF] 37 in the team for spectral mixture analysis, to provide an estimation of the unknown component spectra and the relative concentrations of the underlying species. Indeed, where spectroscopic techniques fail to resolve unknown or unstable individual constituents contained in multi-component mixtures, vibrational spectroscopies combined with BPSS treatment may offer a promising and efficient approach. [START_REF] Carteret | Polymorphism Studied by Lattice Phonon Raman Spectroscopy and Statistical Mixture Analysis Method. Application to Calcium Carbonate Polymorphs during Batch Crystallization[END_REF] Samples with an Au layer thickness of 4, 8, 10, 12, 16 or 20 nm are prepared. The Cu layer thickness is set to 90 nm. Before the Cu/Au multilayers deposition, a chromium film used as adhesion layer is primarily deposited followed by a protecting film of pure Au. This latter layer prevents the oxidation of the Cr layer during the Cu etching by nitric acid. After Cu etching, the SEM imaging of samples with initial Au thicknesses from 4 to 20 nm is performed in order to compare the morphology of the ligaments (Figure 2). The SEM images reveal the formation of Au ligaments for the samples with an initial gold thickness of 4 to 12 nm. Well defined interconnected ligament structures are noticed for 10 and 12 nm Au initial thickness sample and the width of the ligaments is in the range of 40 nm. Above 12 nm Au thickness, the etching process does not lead to the creation of ligaments. In the cross section view SEM images, the typical stacked morphology clearly appears especially for the 16 and 20 nm initial gold thickness whereas for 10 and 12 nm gold thickness, on the plan SEM view, at least 2 ligament stacked layers are visible. The etching process in the liquid phase may lead to the delamination of the last Au surface layers from the substrates, which decreases the number of visible layers.

RESULTS AND DISCUSSION

Effect

After 24 hours incubation in 10 mL of 10 -4 mol.L -1 BP solutions, separately, the different nanoporous Au samples (washed and dried samples) are tested by SERS measurements. Each Raman spectrum is recorded under the same conditions. In these preliminary measurements, a mapping has not been performed, and the illustrated Raman spectra (Figure 3) are the average of 11 simple measurements of the surface. The distinctive Raman fingerprint of BP molecules with the typical features of the active modes is observed in Figure 3(a). Figure 3(b) shows the zoom in the characteristic peak in the region corresponding to the breathing mode of the pyridine rings. [START_REF] Forato | Functionalized Core-Shell Ag@TiO2 Nanoparticles for Enhanced Raman Spectroscopy: A Sensitive Detection Method for Cu(II) Ions[END_REF] This peak will be referenced as BM peak (breathing mode peak). We have followed the integrated intensity of this Raman peak to evaluate the SERS performance of the Au stacked layers as a function of BP concentration because this intensity is not very sensitive to the environment evolution. [START_REF] Chauvin | Lamellar Nanoporous Gold Thin Films with Tunable Porosity for Ultrasensitive SERS Detection in Liquid and Gas Phase[END_REF] The 10 nm initial Au thickness reveals the highest SERS signal for the detection of BP molecules. This behavior is attributed to the specific nanoporous gold architecture, facilitating the diffusion of molecules through the porous network, and also leading to an increase in the number of molecules located in the vicinity of the Au ligaments. For thinner or thicker Au initial layers, this specific ligament morphology is less pronounced or even disappeared leading to a reduction of the intensity of the Raman signal.

Due to the highest SERS signal measured in this preliminary experiment, the 10 nm at initial state is the selected Au thickness in all the upcoming experiments of the study. These samples will be denoted as SERS substrates.

Towards the verification of the homogeneity and the reproducibility of the SERS substrates, different Raman spectra and SERS mappings are measured over different zones on the same substrate (Supporting information, Figure S1) that is incubated in 10 -4 mol.L -1 BP solution, as well as on different similar substrates (Supporting information, Figure S2). BP solutions of 10 -4 mol.L -1 concentrations are considered as a starting point and a reference solution for all other used BP concentrations.

Storage and time stability. The stability and storage conditions of SERS substrates as a function of time are studied because they will determine an actual use as sensors in practical applications. In this study, we have observed that the time between the Au/Cu deposition and the Cu etching has no effect on the SERS behavior of the Au nanoporous film: indeed, after storage for more than a month, the films could be etched and present very similar SERS results equivalent to fresh samples. This result could be explained by the passivation of the surface by the last layer of gold which prevents the Cu layers to be oxidized. Moreover, at room temperature, we assume that the Cu diffusion remains limited and does not change significantly the composition of the stacked layers.

The step of BP adsorption over the substrates has to be done immediately after the etching process in order to avoid the quenching of the Raman signals when the etched samples are kept in ambient atmosphere before BP adsorption (Supporting information, Figure S3). A solution to avoid this problem is to store the etched samples in DI water before incubation in BP solution.

Figure 4 shows the effect of storing the etched samples in DI water before incubation in BP solutions. Herein, after deposition, the samples are etched, and then stored in DI water for different times (3 to 11 days). Then, adsorption of BP solutions of 10 -4 mol.L -1 concentration is performed. After this adsorption process, the samples are dried and immediately tested by Raman measurement. A fresh sample referred as "0 day" is also used as a reference. This sample is washed in DI water just after etching, then directly incubated in BP solution for 24 hours, to finally being dried and measured.

The SERS response is similar for all substrates (Figure 4 The XPS measurements of the Au surface when stored in water compared of that not-stored have been performed. The extracted results are shown in Table 1. The remaining Cu concentration on the surface is very low after etching, and we observe that the ratio of Cu/Au atomic concentrations decreases after storage in water. The position of the Cu 2p XPS peak at 932 eV and the Auger LMM spectrum indicate that after etching Cu is oxidized as Cu+ (Supporting information, Figure S4). [START_REF] Biesinger | Advanced Analysis of Copper X-Ray Photoelectron Spectra[END_REF] This result indicates that residual Cu2O may remain at the surface of the Au ligament. The storage in water during several days leads to the removing of these copper species since the intensity of the typical peaks strongly decreased. An expected typical surface contamination by oxygen and carbon is detected and we can also detect traces of nitrogen. The position of the N 1s XPS peak is close to 398 eV, indicating that this presence cannot be attributed to residual nitrite after Cu etching by nitric acid. The step of rinsing in water just after the etching has removed the residue of copper nitrite from the surface of the sample.

For more investigations, visible extinction spectra of dried Au samples with and without storing in DI water are measured (Figure 5 (a)). For the as-prepared fresh Au samples, there is a strong absorption band peak at about 485 nm, which is referred to the surface plasmon resonance band of Au ligaments. [START_REF] Lang | Localized Surface Plasmon Resonance of Nanoporous Gold[END_REF] After the as-prepared Au sample is immersed in DI water for 11 days, the maximum absorption peak red shifts towards 500 nm. In the SEM images of the surface presented on Figure 5(b&c), we can notice that the ligament morphology is not significantly affected by storage in water. Therefore, we assume that the shift of the plasmon band cannot be attributed to a strong evolution of the Au ligament dimensions. For each point of the Raman mapping, we have integrated the spectral peaks characteristic of BP signature. These integrated values are reported in function of the (x,y) positions of the scanning area. These SERS mappings for BP concentrations in the 10 -8 mol.L-1 to 10 -4 mol.L -1 range are shown in Figure 6 with a linear color scale for intensity between 0 and 15000 a.u. The BP SERS signal is well detected all over the tested area for the 10 -4 and 10 -6 mol.L -1 concentrations and at 10 -8 mol.L -1 low intensities are yet detected in some pixels of the map but masked by the color scale used in the Figure 6. If we decrease the maximum intensity of the colored intensity scale bar by a factor of 10 (from 15000 to 1500), it is possible to observe the repartition of the BP SERS signal on the mapping for the concentrations between 10 -18 and 10 -8 mol.L -1 (Figure 7). Notable SERS signals are observed at BP concentrations lower than 10 -12 mol.L -1 . The number of pixels, where a strong BP SERS signal is detected, decreases as the BP concentration goes down. To give a quantitative representation of this decrease, the average intensity of the mapping as a function of the BP concentration is given in Figure 8 for all the BP concentration ranges which have been studied. Herein, three regions appear: an increase of the intensity for a BP concentration varying from 10 -18 to 10 -12 mol.L -1 , a plateau between 10 -10 and 10 -8 mol.L -1 , and finally an expected increase of the BP SERS intensity with increasing BP concentration from 10 - This behavior comes from the fact that the most enhancing SERS hot spots are occupied by BP molecules even with low concentration because of their high chemical potential of reactivity.

Once these sites are saturated, the BP adsorption continues but with an enhancement SERS factor less efficient and then not displayable on the spectra. The previous quantitative analysis has been checked by driving a statistical analysis without a priori approach on all spectra. By decreasing the BP concentration, Raman signal specific for BP is diluted in a noisy background. Therefore, an advanced mixture analysis algorithm based on Bayesian positive source separation (BPSS) has been applied that allows the recovery of the BP Raman spectra as the main source of the SERS spectra (Figure 9). 39 Figure 9(a) displays the results of the BPSS method obtained from the 845 spectra collected for the 10 -12 to 10 -4 mol.L -1 concentrations. It is worth noting that the spectra consist of three sources whose behaviors as a function of solution concentration are different. The first source corresponds to the Raman spectrum of the BP molecules (Figure 9(a)) and it is very similar to the "usual" SERS signal obtained with colloidal silver/bipyridine system. Accordingly, the weight of this source increases with the concentration but not linearly (Figure 9(b)). In contrast, the weights of the second and third sources are quasi-constant, which suggests that it corresponds to the 'background' contributions arising from the substrate or other pollutions. The spectra of the two other sources are given in supporting information (Figure S6). The weight of the BP source versus concentrations (Figure 9(b)) is clearly similar to the first analysis (Figure 8) between 10 -12 to 10 -4 mol.L -1 . Therefore, this approach without a priori validates the first analysis displayed in Figure 7 and Figure 8. For lower concentrations than 10 -12 mol.L -1 , the BPSS becomes a no-valid approach because the statistical abundancy of the source 1 (BP source) extracted previously is too weak.

Even at ultra-low BP concentrations, such as 10 -18 mol.L -1 , a significant BP SERS signal is detected and hot spots are visible on few pixels of the mapping. At this very low BP concentration level, we can estimate the density of the adsorbed molecules per unit area of the substrate. An estimation could be done taking into consideration the volume of the BP solution used for the adsorption tests (10 mL), and the surface area of the samples (8 x 8 mm 2 ), with an assumption that all the molecules present in the solution are uniformly adsorbed on the gold surface. For a concentration as low as 10 -18 mol•L -1 the quantities of molecules on the active surface could be estimated to about 102 molecules/mm 2 . Taking into account that the laser spot size is of 3 µm diameter, the total surface tested by the mapping (13 x 13 pixels) remains very low (about 10 -3 mm 2 ). In this case, the probability of the presence of a molecule on the mapping surface is therefore very low (around 0.1). This low probability means that the possibility to detect BP molecules at such a level of ultra-low concentration of 10 -18 mol.L -1 is not completely assured, and is subjected to chance. Moreover, we can observe on the mappings (Figure 8) that the presence of high BP SERS intensity on pixels, attributed to the hot spots, is generally linked to the high BP SERS intensity in the neighboring pixels. We therefore suspect that there is an additional Raman mechanism that induces the signal propagation in the layer beyond the laser excitation zone, and contributes to an increase of the actual SERS tested area. To verify this possibility, we have realized spatially offset Raman spectroscopy. First, the µ-SESOR test is performed on the surface of a silicon wafer. It shows that the signal could not be collected for an offset of 6 µm and above (Supporting information, Figure S7). In µ-SESOR experiments performed on the stacked gold nanoporous layer samples, the intensity of the BM peak remains significant up to 15 µm and a low signal is still detected at 25 µm. This result indicates that the Raman signal can be collected from an area much larger than the area defined by the excitation laser beam. Moreover, compared to a single nanoporous layer, the fact that the samples are constituted of several stacked nanoporous layers could also contribute to the increase of the sites for the creation of hot spots. Finally, these two combined effects, i.e. the increase of the tested area due to the propagation of the electronic field along the ligaments parallel to the substrate surface, and the increase of the density of hot spots, may therefore lead to the possibility to detect ultra-low concentrations of BP molecules by SERS. In order to verify that the BM BP peak does not correspond to the Au surface itself, the evolution of the Raman signal of the blank sample without the presence of BP in the room is examined with time. For that, the Raman of the blank sample is repeated after 24 hours, and then a Raman spectrum is taken at the same area for 3 consecutive times each 10 minutes where no BM peak is found. This very interesting result requires also a comment directly related to the measurement performed with BP solution at very low concentration in the previous section. For instance, at a BP concentration of 10 -18 mol.L -1 , the concentration of BP molecules in water is about 600 mL -1 .

Gas

It means that the slightest contamination of vials or air by BP molecules can lead to a very significant error in the value of the real concentration of the BP solution. The protocol for diluting and preparing BP solutions must therefore be carried out under extremely rigorous conditions of cleanliness and air purity control. Tests at concentration lower than 10 -18 mol.L -1 are hardly possible because of the risk of this remaining pollution. This remark is confirmed by tests of BP solution that we have performed at concentration of 10 -20 mol.L -1 for which SERS mapping can show some intense pixel in the BP SERS mapping that we attribute to this residual BP pollution. However, this experiment highlights the ability of the nanoporous Au substrates to detect the presence of BP molecules in room air.

Recovery of Au sample substrates.

In contrary of storing the samples before Raman measurements, Au samples could be reversibly restored. This means that, after BP adsorption on the substrates, the BP molecules could be removed by means of acetic and nitric acid, and then, re-adsorption of the same BP solution is carried on. The results are shown in Figure 14.

The as-prepared sample has intense peaks at 1011 cm -1 and 1058 cm -1 (black colored spectrum in Figure 14(a&b)). This sample is stored at room temperature for 6 months (6 months aged), the corresponding peaks are shifted into 1032 cm -1 and 1060 cm -1 (red spectrum in Figure 14(a &c)). After the immersion of the aged sample in acetic acid for 15 minutes, a shoulder at 1013 cm -1 appeared (blue spectrum in Figure 14(a&c)), which means that the aged sample is polluted due to storage, preventing the BM peak corresponding to BP molecules to be observed. After the immersion in acetic acid for additional 15 min, the sample showed no peaks (green spectrum in Figure 14(a &c)). This means that the sample is totally recovered with removing all the BP adsorbed on the surface. The sample is then immersed in 16 mol.L -1 HNO3 concentrated acid to insure the removal of all the BP residuals possibly remaining on the surface (violet spectrum in Figure 14(a&c)), where no BP peaks are observed. Therefore, the SERS substrate is totally cleaned and recovered to be re-used for another purposes. Finally, the sample are incubated in BP solution in order to evaluate the ability of our recovered surface for another detection. BP peaks re-appeared at 1013 cm -1 and 1060 cm -1 although in lower intensities (orange spectrum in Figure 14(a&b)). Therefore, the sensitivity to BP detection is recovered and sample could be used another time as SERS substrate.

CONCLUSIONS

SERS substrates using multi-stacked layers of nanoporous gold ligaments can be used for molecules detection in liquid or gaseous media. Our method using a thorough analysis of the SERS mapping allows to reach extremely low detection limits for BP molecules reaching concentrations as low as 10 -18 mol.L -1 in water solutions. An additional Raman mechanism has been revealed by µ-SESORS that may appear in the very specific morphology of the nanoporous gold stacked layer. In this mechanism, oriented Au ligaments parallel to the surface play the role of waveguide for the incident light, thus widening the effective area for the SERS. Combined with this effect, the intrinsic nature of the stacked Au nanoporous structures leads to an increase of the density of hot spots.

This study has also shown the good homogeneity of the SERS surfaces and the reproducibility of the process to prepare the SERS substrates. These type of SERS substrates can be stored before their use and present an ability to recover their properties after molecules adsorption and Raman tests.

All these results open the way to the development of real SERS sensor for the detection of molecules at very low concentrations in liquids or gases. The SERS signal, at low concentration, appears to be not linear dependent of the value of the molecule concentration. Therefore, it is clear that application of these SERS substrates has to be dedicated to detection of the presence or not of molecules at low concentration in the media. This point represents the most promising perspective. The next research will be focused on the test of these SERS substrates for other type of molecules representing an interest for application such as detection of neurotoxic gases.

ASSOCIATED CONTENT Supporting Information. Raman spectra and Raman mapping for: homogeneity and reproducibility evaluation, effect of time between etching and BP adsorption, BP active modes in Raman spectra. XPS spectra. Raman spectra for source 2 and 3 deduced from BPSS approach. µ-SESOR on Si

  of the initial gold thickness. The representation of the obtained stacked Au-Cu layers is shown in Figure 1. Films are deposited on silicon substrate by DC magnetron sputtering in a pure argon plasma. The deposition of each Cu and Au layer is repeated alternatively five times allowing the growth of the desired five Cu/Au stacked layers.

Figure 1 .

 1 Figure 1. Schematic representation (a) and SEM image (b) of the cross-section with the five

Figure 2 .

 2 Figure 2. Cross section and plan view SEM images showing Au samples with different

Figure 3 .

 3 Figure 3. Raman spectra on Au samples incubated in 10 -4 mol.L -1 BP solutions for 24 hrs after

  (a)). By examining the BM peak integrated intensity (Figure4(b)), the fresh samples without storage (0 days) shows the highest intensity. Storing the samples for 3 days prior to BP adsorption results in a low decrease of the intensity and this way appears to be very useful and efficient to keep pristine Au nanoporous SERS substrates until their use for BP adsorption and SERS analysis. The intensity is stabilized at around 65 % of the initial value for further storage in DI water after 8 and 11 days. We suppose that the decrease in the Raman signal results in a competition between water and BP molecules for their adsorption on Au ligament surface.

Figure 4 .

 4 Figure 4. Raman spectra on samples with a thickness of 10 nm of gold and 90 nm of copper at

Figure 5 .

 5 Figure 5. (a) Extinction spectra for the as-prepared fresh Au substrate and Au substrate stored

Figure 6 .

 6 Figure 6. 13x13 SERS mapping of the samples after adsorption of BP solutions of

Figure 7 .

 7 Figure 7. 13x13 SERS mapping of the samples after adsorption of BP solutions of

Figure 8 .

 8 Figure 8. Evolution of the average integrated intensity of the BM peak of bipyridine deduced

Figure 9 .

 9 Figure 9. (a) Source 1 spectrum for concentrations from 10-12 to 10-4 mol.L-1 extracted from

µ-

  Surface Enhancement Spatially offset Raman (µ-SESOR) spectroscopy. In µ-SESOR, the collection and excitation areas are disjoined by shifting the excitation area and remaining fixed the collection area. The shift between the laser excitation (bright spot on Figure 10(a) and red point on Figure 10(c)), and the collection area (crosshair Figure 10(a) over the sample and blue point Figure 10(c)), is obtained by the lateral moving of the laser beam in front of the lens.

Figure 10 .

 10 Figure 10. (a,b) Images under the Raman microscope clarifying the separation of the

Figure 11 .

 11 Figure 11. µ-SESOR experiments performed on the stacked gold nanoporous layer samples

Figure 12 .

 12 Figure 12. Schematic illustration of the mapping taken for SERS and µ-SESOR experiments.

  detection. The sensitivity of the substrates is tested by measuring the ability to detect an enhanced Raman signal of BP gas vapor under ambient air. In this experiment, the Raman signal of a blank sample which corresponds to pure nanoporous gold is measured. During the SERS measurement, small grains of BP powder are placed at about 1 cm space near the sample, and the same mapping is taken every 40 minutes. The results are shown in Figure13. Before the introduction of the BP grains in the room, the blank sample shows only a noise Raman signal in the range of the BP BM peak. At 40 minutes after the introduction of the BP grain near the substrate, a shoulder begins to appear on the spectrum orresponding to BM peak. Raman spectrum is measured each 40 minutes, to reach 225 minutes. The intensity of the BM BP peak is highly increased with time of exposition to the presence of BP grain.

Figure 13 .

 13 Figure 13. Raman spectra on a sample with a thickness of 10 nm of gold and 90 nm of copper

Figure 14 .

 14 Figure 14. Average Raman spectra of a sample with a thickness of 10 nm of gold and 90 nm

Table 1 .

 1 Atomic concentrations extracted from XPS data for the non-stored (0 days) fresh sample and 11 days water-stored Au sample, in addition to the ratio of the atomic concentration of Cu/Au before and after water-storage.

		C (at.%) N (at.%) O (at.%) Cu (at.%) Au (at.%)	Cu/Au
	Fresh sample (0 days)	34.6	1.6	14.0	2.7	47.1	0.06
	Water-Stored sample (11 days)	48.0	4.7	19.1	0.3	27.9	0.01

to 10 -4 mol.L -1 . The presence of a plateau for BP concentration has been previously observed.[START_REF] Zalm | Review-Recent Advances in the Development of Nanoporous Au for Sensing Applications[END_REF]

AUTHOR INFORMATION Corresponding Author

Pierre-yves.tessier@cnrs-imn.fr

ACKNOWLEDGMENT