G Brogna

J Antoni

Q Leclere

O Sauvage

Engine noise separation through Gibbs sampling in a hierarchical Bayesian model

Keywords: Referenced source separation, Bayesian computation, Hierarchical modelling, Engine noise

). As opposed to them, the present algorithm autonomously chooses one or the other depending on the characteristics of the analysed signals. Even if the development context is the separation of the sources in an IC engine, the presented method is general and can be applied to any source separation problem.

Introduction

Of major concern with NVH characteristics of car vehicles is the control and harmonization of levels of multiple noise sources. One source of primary interest is the engine noise, particularly in its diesel variant. Its sound quality remains one of the most discriminating characteristics in front of the gasoline engine. As a consequence, diesel engines have been the object of several studies aiming at analysing the reasons for this poor acoustic quality [START_REF] Hastings | Sound quality of diesel engines[END_REF][START_REF] Groemping | Customer driven diesel vehicle sound quality[END_REF][START_REF] Alt | Diesel engine cold start noise improvement[END_REF][START_REF] Renard | Extraction of vibration sources in diesel engines[END_REF]. Several aspects have been studied, above all the combustion whose irregularities are impacting the perception of the diesel engine sound [START_REF] Lowet | Development of a metric to quantify diesel engines irregularities[END_REF][START_REF] Sauvage | Diesel knock noise from combustion phenomenon to perceived signals[END_REF]. Experimental results on the equipment isolated from the working engine are accurate, but not truly indicative of the behaviour of an engine in actual working conditions. For this reason, signal separation methods have become largely applied in the automotive industry: they allow the separation of the contribution of each equipment from the overall noise of a regular engine. The separation is a difficult task since the noises are highly overlapping in both the time and the frequency domains, especially in the case of impulsive events. The coherence based filtering methods, such as the Wiener filter, present two major limitations:

• they have been developed considering the existence of one reference for each phenomenon whose noise has to be extracted. This means that the reference entirely observes the contribution of the source. Actually, such a unique reference is available only for very few noise sources. This yields to the need to use references that are only partly coherent to a single source, but which may be multiple;

• the reference also needs to be "pure", meaning it is not polluted by any additive noise (even uncoherent with other sources). If multiple references are considered for a single source, a subspace filtering method [START_REF] Sharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF] can be applied in order to denoise them. However, the latter assumes that the Signal to Noise Ration (SNR) of all references is the same. In practice, if the references have different nature the hypothesis of equal SNR cannot be assumed.

Some methods for source extraction succeed in the separation even in operating conditions [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF][START_REF] Xianhua | Blind separation of internal combustion engine vibration signal by a deflation method[END_REF][START_REF] Hirano | Using multiple regression analysis to estimate the contribution of engine-radiated noise components[END_REF]. For instance, Antoni et al. presented a general separation method based on the Wiener filter in Ref. [START_REF] Antoni | Novel source separation method for measuring combustion noise and structural loss[END_REF]. El Badaoui et al. [START_REF] Badaoui | Servière Separation of combustion noise and piston-slap in diesel engine, part I: separation of combustion noise and piston-slap in diesel engine by cyclic Wiener filtering[END_REF] used the Wiener filtering method in the context of cyclo-stationary processes [START_REF] Antoni | Apports de l'échantillonnage angulaire et de la cyclostationnarité au diagnostic par analyse vibratoire des moteurs thérmiques[END_REF][START_REF] Antoni | Cyclostationarity by examples[END_REF]. Afterwards, Pruvost et al. developed an improved version of the filter in Ref. [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF]. Finally, this approach is extended to the cyclo-non-stationary regime in [START_REF] Antoni | Separation of combustion noise in IC engines under cyclo-non-stationary regime[END_REF]. These methods answer well to the case of a single but not pure reference. In order to make it pure, its random part -obtained after substraction of the synchronous average -is isolated and used for the computation of the Wiener filter instead of the whole signal. As a consequence, another definition of the Wiener filter arises in the cyclo-stationary context. This new definition is proved to be better than the former [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF] when the random part of the reference is more energetic than the deterministic part. This is not always the case in an operating engine.

A filter with a larger application field would be a linear combination of the classic Wiener filter and the one presented above, weighting them by the mean squared error of their estimation. However, the calculation of the mean squared errors assumes the knowledge of the noise characteristics, which cannot be estimated until after the estimation of the filters and contributions. The method proposed in this work aims at building a hierarchical model in order to estimate the noise characteristics while estimating the separation filters. Moreover, it generalizes the extraction to the case of multiple polluted references with different SNRs. The Bayesian approach seems the most natural to this aim: a Bayesian hierarchical model of the system is proposed.

In the Bayesian context, the posterior joint probability density function is often fairly complex and a direct sampling of it is not possible. Simple Monte Carlo methods, such as the importance sampling, may fail or be slow when the posterior distribution is concentrated with respect to the prior. In this cases using Markov Chain Monte Carlo (MCMC) methods may be preferable. These methods involve the construction of an ergodic Markov chain whose stationary distribution corresponds to the target posterior probability density function. This way, after having attained the convergence of the chain, the samples from the Markov process are samples from the target distribution as well. Several applications of MCMC methods in different fields of science, such as biology [START_REF] Ahmadian | Efficient Markov Chain Monte Carlo Methods for Decoding Neural Spike Trains[END_REF], mechanics [START_REF] Green | Bayesian system identification of dynamical systems using large sets of training data: A MCMC solution[END_REF] and space research [START_REF] Huang | MCMC-Particle-based group tracking of space objects within Bayesian framework[END_REF], have focused on their advantages. The proposed approach will infer the unknown elements through the Gibbs sampler, a popular MCMC algorithm, typically well adapted to hierarchical Bayesian models. First of all, the Gibbs sampler samples from the joint posterior distribution using the conditional posterior distributions of the parameters: this makes possible the use of far simpler probability density functions with remarkable properties (such as Normal laws). Moreover, it keeps the specificities of a Bayesian model. For instance, as opposed to optimisation techniques, the Bayesian credible interval on the results can be easily obtained. The convergence of the application will be verified through parallel chains and the Gelman & Robin convergence test [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF].

The first section of this work defines the problem of source separation and its solution using the Wiener filter [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF]. The limits of this method come when the references are not pure: they are pointed out using signals measured on a diesel engine. Therefore, a new model featuring latent variables is described in the same section. The hierarchical Bayesian approach to the problem and the application of the Gibbs sampler are the subject of the second section. It is shown in particular that the obtained filter is a linear combination of the Wiener filters cited above. The third section puts the method in practice in a numerical simulation with synthesized engine signals and its results are compared to the ones provided by a Wiener filter. The use of synthesized signals allows the comparison of the methods to known solutions. The application to real engine signals allows to actually show that the obtained filter can be interpreted as a combination of the Wiener filters: this is addressed in the fourth section. Finally, some perspectives on the method will be provided in the conclusion. The engine application is just one among others, the presented method updates the existing ones and can be applied to a large scope of separation problems.

2 Problem statement

Signal modelling

Notational conventions: in this paper, a lower-case character stands for a scalar value, a lower-case bold character stands for a vector and an upper-case bold character stands for a matrix. The problem is first formulated in the time domain -with respect to the time variable tand then solved in the frequency domain -with respect to the frequency variable f . The problem consists in the separation of referenced sources. Be an overall measured signal named d(t), composed of the K contributions that have to be separated, x k (t), k = 1, ..., K, and an overall residual noise n(t) independent of the x k (t)'s:

d(t) = K k=1 x k (t) + n(t). (1)
For each contribution, one reference r k (t) is available and the relation between the reference and the contribution through a theoretical filter hk (t) can be written as follows:

x k (t) = (r k * hk )(t) (2)
where the operator * stands for convolution. Figure 1 shows a schematic view of the model chosen for the overall measured signal d(t).

x 1

x 2 1) and (2).

x k d n r 1 r 2 r k h-1 k h-1

Wiener filters and their limits

The Wiener filter method as explained for instance in Ref. [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF] is able to solve such a problem and provide an estimate for the impulse response hk , noted h k (t), if the references are pure and uncorrelated. A reference is pure if: 1) its signal contains information on the considered source alone and 2) if the information is enough to reconstruct the source perfectly through a linear filter, i.e. the coherence function between the reference and the source is everywhere equal to one. If these conditions are respected, then the method can be applied separately to each contribution. In practice in a Diesel engine this is rarely the case: all signals from the same engine are rather correlated, being all generated by the same rotating machine. The solution proposed in the context of cyclostationarity is to extract and consider just their random parts. For instance, figure 2 shows the cylinder pressure in blue: this reference is typically pure, but it is correlated to other engine signals. Its random part is obtained dividing the signal into cycles and computing the mean along the cycles: this yields an estimate of the first order cyclostationary (CS1) part of the signal, also called the synchronous average (SA). Then, removing the SA from each cycle of the signal leaves its random part (Fig. 3). Thus, two possible Wiener filters can be defined: the one obtained from raw signals will be called h k (t) (Fig. 4a), the one obtained from the random part of the signals (removing the synchronous average) will be called h c k (t) (Fig. 4b). The random parts being usually not correlated to the other signals by the mere fact of originating from the same machine, h c k (t) boasts a smaller bias error. On the other hand, this manipulation decreases the SNR, increasing the random error of the estimator. Therefore, the two filters coexists, each one being better than the other in some specific cases. The advantages of h c k (t) over h k (t) have been pointed out in Ref. [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF][START_REF] Leclère | Source separation in diesel engines using Wiener filtering: physical interpretations[END_REF] through an application to the combustion noise. Actually, the combustion noise references (the cylinder pressure signals obtained from in-cylinder pressure probes) are well suited to this kind of application since their deterministic part is negligible from a frequency low enough (Fig. 5). This means that using the random part of the references does not decrease the SNR notably, but indeed decorrelates the signals leading to a better estimation of the filter. For all other contributions, the limits of the method are reached. For other kinds of sensors, such as accelerometers, the deterministic and random parts of the signal keep comparable values over all the frequency domain of interest (Fig. 6). As a consequence, removing the deterministic part of the reference decreases the SNR, not necessarily yielding to a gain in the application of this approach. Moreover, the references must be pure. The random part of the cylinder pressure is a notably pure reference for the combustion noise. On the other hand, the references coming from classical sensors such as accelerometers and microphones are hardly pure. For instance, an accelerometer placed on one injector sees all the combustion events, not just the ones from the corresponding cylinder (see the accelerometer signal in Fig. 2). Usually a correct windowing strategy allows the interesting part of the signal to be isolated, for instance windowing the references around the combustion events of the targeted cylinder. Nevertheless, this approach may be limited when a less impulsive source is present (e.g. the high pressure pump in a diesel engine). One more general solution to this issue is to apply a waterfall approach (Fig. 7). In a multi-source problem, the waterfall approach is the following:

Time

Physical quantity

1. separate the source whose reference is the purest; 2. find its estimated contribution xk (t) to the overall measured signal;

3. remove its contribution to the references of the remaining sources, in order to have cleaned references;

4. go back to the beginning and repeat using the following purest reference.

Figure 7 summarizes the procedure. The measured signal d(t) is used with the purest reference r 1 (t) to compute a Wiener filter h 1 (t) (or h c 1 (t)). This filter yields the first source contribution x1 (t). After removing this contribution from d(t), the remaining information is used together with the second reference r 2 (t) to compute the filter h 2 (t) (or h c 2 (t)), yielding then the second contribution x2 (t). The process is reiterated for each reference, from the purest to the noisiest.

d(t) h 1 (t) r 1 (t) * x1 (t) r 2 (t) * x2 (t) h 2 (t) - - d(t) -x1 (t) d(t) -x1 (t) -x2 (t)
1 st contribution 2 nd contribution 3 rd contribution Figure 7: Explanatory diagram for the waterfall strategy.

The choice between h k (t) and h c k (t) is of utmost importance for the waterfall approach and it leads to different estimations of the contributions. For instance, figure 8 shows the spectra of the injectors contribution x2 (t) to the noise emission of an engine d(t), computed after a first waterfall step yielding the combustion contribution x1 (t). One curve is obtained using in the first step an h 1 (t) Wiener filter, whereas the other derives from the application of an h c 1 (t) filter: they are sensibly different. Actually, in an engine application, the difference in meaning between the two filters has been pointed out by Leclère et al. [START_REF] Leclère | Source separation in diesel engines using Wiener filtering: physical interpretations[END_REF]. If the first source considered is the combustion (k = 1), h 1 (t) extracts the combustion noise and the coherent mechanical noise, while h c 1 (t) extracts the combustion noise alone. As a consequence:

• if h c 1 (t)
is used for the first step of the waterfall, then only the combustion contribution is removed from the remaining source references. The mechanical noise of other moving parts is left polluting them. Thus, the contribution extracted using these "partially" cleaned references originates not only from the corresponding source, but contains also some mechanical noise coherent for instance to the combustion itself. This is why in figure 8 the blue spectrum is more energetic than the red one;

• if h 1 (t) is used for the waterfall, then the combustion noise and the coherent mechanical noise are completely removed from the references. This yields an "over-cleaning" of the references. Indeed, any mechanical source in an engine is composed of a de-terministic and a random part, so that its reference should be representative of both. The deterministic part, which repeats itself at each cycle, is strongly linked to the operation of the machine itself. This means that the deterministic parts of all sources (and references) are rather coherent simply from the fact of belonging to the same rotating machine. This coherent part is what has been called above "coherent mechanical noise". Therefore, removing it from a reference means not only removing the pollution from the other sources, but also removing information on the targeted source. For instance in engine applications, it is reasonable to assume that the considered source contribution has a periodic component to it since it follows the rotating process. However, the contribution spectrum obtained by h 1 (t) waterfall hardly shows this behaviour (Fig. 8): it seems that the harmonic part in the reference has been cleaned too in the waterfall process.

One of the limits of the Wiener filter method begins to stand out: in the mono-source as in the multi-source cases, two definitions for the Wiener filter coexist, both valid inside their field of application, but none generally accurate, posing a problem of choice.

Figure 8: Spectrum of the injectors contribution to the engine noise. Solid red curve: result obtained using the waterfall strategy with an h 1 (t) filter. Dashed blue curve: result obtained using the waterfall strategy with an h c 1 (t) filter.

Finally, note that what stated above can be applied to the case in which multiple references are available for each source. It can be said that a group of references is pure if it allows the reconstruction of the source contribution through a linear combination of linear filters. Therefore, it is not pure if:

1. the combination of the references does not entirely observe the source (for instance a frequency range is not available);

2. each reference is polluted by an additive noise.

Increasing the number of references for the same source helps to overcome the first limit, however it does not necessarily help to overcome the second. In order to do so, the references can be denoised through a decomposition on signal and noise subspaces [START_REF] Sharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF]. Nevertheless, this method assumes a spatially white noise (homogeneous among the references), which is a hypothesis that does not stand if the references are of different natures. After all that is said, the need for a method that computes a filter adapted to the characteristics of the considered signals and that manages to take into account several correlated references with heterogeneous SNRs is evident. This is the aim of the method proposed in this paper.

Consideration of correlated references

The engine signals will be treated in the framework of cyclostationarity and cycloergodicity. In a broad definition a cyclostationary signal is one that exhibits some hidden periodicity of its energy flow [START_REF] Antoni | Cyclostationarity by examples[END_REF]. Each one-cycle-long portion of the signal can be considered as a realization of a stochastic process governed by constant characteristics (in our case: the Fourier coefficients corresponding to one frequency form a stationary and ergodic process along the cycles). The Wold isomorphism [START_REF] Antoni | El Badaoui Cyclostationary modelling of rotating machine vibration signals[END_REF] allows the model describing the mixing of the sources in a frequency band [f l , f u ] to be written as:

r i (f ) = G(f )z i (f ) + n ri (f ) (3)
d i (f ) = h t (f )z i (f ) + d 0 (f ) + n di (f ) (4)
where • t stands for the transpose operator, i ∈ {1, ..., I obs } with I obs the number of considered cycles and r i (f ) ∈ C Nr is the column vector of the N r measured references at frequency f . As stated in Eq. ( 3), these correlated references are considered as a linear composition of N z N r uncorrelated latent variables stored in the vector z i (f ) ∈ C Nz and an additive noise n ri (f ) ∈ C Nr . The matrix G(f ) ∈ C Nr×Nz describes the relation between the latent references z i and the measured ones r i at a frequency f . Equation ( 4) describes how the sources contribute to the global output signal d i (f ) ∈ C. It is considered composed of a deterministic part independent from the sources d 0 (f ) ∈ C, a series of contributions from the sources of interest and a zero mean additive noise n di (f ) ∈ C. The column vector h(f ) ∈ C Nz contains the transfer functions for each latent variable at a frequency f . The product between the transfer function and the latent variable yields to the contribution of the considered source to the global noise.

From here on, the dependence on the frequency will be omitted to simplify the notation.

Physical meaning of the latent variables

As the model has been presented in the previous section, the inputs of the method are the measured references r i and the global response

d i .
The outputs are all the other terms of the equations. In order to ensure the identifiability, the matrix G is assumed to have a triangular inferior shape. Actually, if this matrix were full, the problem stated by Eq. ( 3) would resume to an eigenvalue decomposition, without any direct way to know the physical meaning of the variables. On the other hand, the triangular inferior shape ensure the physical meaning of the latent references to be inferred from the measured ones. This implies that the measured references must be wisely chosen in order for the problem to be solved keeping a physical meaning. Take for instance a cycle i, let r 1 to r 6 be the measured references, which the user considers well describing the sources he/she would like to extract from the overall engine noise. Let these sources be A, B and C, thus let z ∈ C 3 the column vector of the corresponding latent variables. The link between the sources and the latent variable is ensured by the measured references and the matrix G ∈ C 6×3 . Equation (3) becomes:

        r 1 r 2 r 3 r 4 r 5 r 6         =         G 11 0 0 G 21 G 22 0 0 G 32 G 33 G 41 0 G 43 G 51 G 52 G 53 G 61 G 62 G 63           z 1 z 2 z 3   +         n r1 n r2 n r3 n r4 n r5 n r6         . ( 5
)
As it can be seen the first row implies that the physical meaning of the first latent variable is directly linked to the first measured reference alone. If the measured reference r 1 is chosen wisely as to be pure and contain information from the source A alone, it can be implied that the latent variable z 1 describes this source too. Having determined the first latent reference, the second row allows the determination of z 2 even if the measured reference r 2 contains information about both source A and B. Finally, all other references are not required to be pure, since, having determined the first two latent references, z 3 as a latent pure reference on the source C can be determined as well.

It can be remarked that Eq. ( 5) describes the spatial decomposition [START_REF] Sharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF] on the signal subspace, generated by the virtual references z i , and the noise subspace, generated by the noises n ri . Nevertheless, the method proposed in the present paper does not need to assume the homogeneity of the noises, i.e. with identical variances. As a consequence, it can be applied to references having different natures.

3 Inference of parameters through Gibbs sampling

Principles

This section presents the way of solving the problem through the implementation of Gibbs sampling in a hierarchical Bayesian model. Let x be the vector storing the unknown parameters and Y the measurements, the application of the Bayes' rule allows obtaining the joint posterior probability density function (pdf) of x as:

[x|Y ] ∝ [Y |x] • [x]. (6)
In Eq. ( 6) the following elements can be recognised:

• [x|Y ]: the posterior joint pdf of the unknown parameters, i.e. the joint pdf of the parameters given the measurements Y .

• [Y |x]: the likelihood function, which stands for the probability of observing the data Y given a set of unknown parameters x.

• [x]: the prior pdf of the unknown parameters, which reflects what we know about the parameters before having collected the data Y .

The Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF] is a procedure which, among the MCMC methods, allows Bayesian inference on the unknown parameters using conditional pdfs of one parameter given the others, instead of the joint pdf of all parameters. In this procedure, we will cycle through the parameters drawing each one in its conditional pdf given all other parameters at their most recent values. Each sample of a parameter is a step of a Markov chain whose long-run distribution represents the marginal pdf of the parameter. This method is particularly well suited in a hierarchical model [START_REF] Bolstad | Understanding computational Bayesian statistics[END_REF].

For the application to the separation problem, each term in equations ( 3) and ( 4) is considered as a random variable. The hypothesis on the stationarity and ergodicity of the signals hopping from one cycle to the other yields to the real and imaginary parts of the random variables to be independent, which facilitates considerably the calculations. The noises n di and n ri are included in the measurements, which are then considered issued from circularly symmetric complex normal laws, given the other parameters. Thus, the likelihood functions are:

(3) → [r ji |g j , z i ] ∼ N C (g j z i , σ 2 rj ), ∀i = 1, ..., I obs , ∀j = 1, ..., N r (7) (4) → [d i |h, z i , d 0 ] ∼ N C (h t z i + d 0 , σ 2 d ), ∀i = 1, ..., I obs (8)
where g j is the j th row of the matrix G and σ 2 d , σ 2 rj are the variances of the zero-mean Gaussian noises in equations ( 3) and ( 4). Each unknown parameter is considered following a probability density distribution whose parameters are considered as random variables as well. Figure 9 shows the associated hierarchical system. Hereafter, there are prior pdfs that have to be chosen for each random variable, where InvGamma(a, b) stands for an Inverse-Gamma law of shape and rate parameters respectively a and b:

[d 0 ] ∼ N C (0, σ 2 d0 ) [σ 2 d0 ] ∼ InvGamma(a 0 , b 0 ) [σ 2 d ] ∼ InvGamma(a d , b d ) [σ 2 rj ] ∼ InvGamma(a r , b r ), ∀j = 1, ..., N r [z i ] ∼ N C (0, I), ∀i = 1, ..., I obs [h t ] ∼ N C (0, diag(σ 2 h1 , ..., σ 2 h Nz )) [σ 2 h k ] ∼ InvGamma(a h k , b h k ), ∀k = 1, ..., N z [g j ] ∼ N C (0, diag(σ 2 Gj1 , ..., σ 2 G jNz )), ∀j = 1, ..., N r [σ 2 G jk ] ∼ InvGamma(a G jk , b G jk ), ∀j = 1, ..., N r , ∀k = 1, ..., N z . (9)
σ 2 In accordance with the Bayesian method, the prior laws contain the subjective view or the actual knowledge of the user on the parameters. They are usually chosen quite flat in order to reflect the lack of knowledge on the parameters before any measurement. The hyperparameters (a, b) are chosen by the user, but the higher they are in the hierarchical model, the less influence they have on the posterior pdf of the random variables [START_REF] Bolstad | Understanding computational Bayesian statistics[END_REF]. In this case of application, the following information has been assumed known on the priors:

d 0 a 0 b 0 d 0 σ 2 d a d b d σ 2 h k a h k b h k h k z i σ 2 r j a r b r σ 2 G jk a G jk b G jk G jk r ji d i
• the variances are supposed to follow Inverse Gamma laws since a variance can only have positive values; • all the other parameters are supposed to follow circularly symmetric complex Normal laws; • elements of the Gaussian random vector h are a priori uncorrelated; • h is heteroskedastic, i.e. each element of h can have a different variance. Each variance follows an inverse gamma distribution; • elements of G are a priori uncorrelated; • G is heteroskedastic. Each variance follows an inverse gamma distribution; • z i follows a multivariate circularly symmetric complex normal distribution whose covariance matrix is the identity matrix I. This means that the latent variables are supposed uncorrelated. An identity covariance matrix is imposed in order for the variances of h and G to be identifiable following Eqs. ( 3) and (4);

• in order to obtain an estimation of z i for each cycle and rebuild the time domain signal, the i observations are considered independent realizations of a stochastic process. The Bayesian inference for z i is thus realized for each cycle separately.

Finally it can be noted that Inverse Gamma distributions have been chosen for the variances and Gaussian ones for the means, since these two distributions are conjugate priors for the likelihood1 . This is a sensible choice in order to greatly improve the performances of the following Gibbs sampler [START_REF] Bolstad | Understanding computational Bayesian statistics[END_REF].

Proposed Gibbs sampler

In what follows the operator • H is the Hermitian transpose and rest stands for all the random variables, except the random variables whose pdf is being expressed. The calculations are detailed in Appendix B.

The application of the procedure follows 11 steps.

Initialize the values of

d 0 , σ 2 d0 , σ 2 d , σ 2 rj , h k , σ 2 h k , G jk and σ 2 G jk . 2. ∀i = 1, ..., I obs , draw a sample from [z i |rest] ∼ N (µ zi , Ω zi ) with µ zi = Ω zi (h * (d i -d 0 )/σ 2 d + G H   σ -2 r 1 0 0 0 . . . 0 0 0 σ -2 r Nr   r i ) Ω zi = (h * h t /σ 2 d + G H   σ -2 r 1 0 0 0 . . . 0 0 0 σ -2 r Nr   G + I) -1 .
where • * stands for the complex conjugate and • t for the matrix transpose.

Draw a sample from [σ 2

d |rest] ∼ InvGamma(I obs + a d ,

I obs i=1 |d i - h t z i -d 0 | 2 + b d ). 4. ∀j = 1, ..., N r , draw a sample from [σ 2 rj |rest] ∼ InvGamma(I obs + a r , I obs i=1 |r ji -g j z i | 2 + b r ). 5. ∀k = 1, ..., N z , draw a sample from [σ 2 h k |rest] ∼ InvGamma(1 + a h k , |h k | 2 + b h k ). 6. ∀j = 1, ..., N r and ∀k = 1, ..., N z , draw a sample from [σ 2 G jk |rest] ∼ InvGamma(1 + a G jk , |G jk | 2 + b G jk ). 7. Draw a sample from [h t |rest] ∼ N (µ h , Ω h ) with µ h = I obs i=1 (d i -d 0 )z H i Ω h /σ 2 d Ω h = ( I obs i=1 z i z H i /σ 2 d +    σ -2 h 1 0 0 0 . . . 0 0 0 σ -2 h Nz   ) -1 .
8. ∀j = 1, ..., N r , draw a sample from [g j |rest] ∼ N (µ gj , Ω gj ) with

µ gj = I obs i=1 r ji z H i Ω gj /σ 2 rj Ω gj = ( I obs i=1 z i z H i /σ 2 rj +    σ -2 G j1 0 0 0 . . . 0 0 0 σ -2 G jNz   ) -1 .
9. Draw a sample from [d 0 |rest] ∼ N (µ d0 , ν d0 ) with

µ d0 = ν d0 I obs i=1 (d i -h t z i )/σ 2 d ν d0 = (σ -2 d0 + I obs σ -2 d ) -1 .
10. Draw a sample from [σ 2 d0 |rest] ∼ InvGamma(1 + a d0 , |d 0 | 2 + b d0 ). 11. Go back to step 2 and repeat the process until a sufficiently large sample is collected after convergence.

The sequence of samples produced for one of the random variables forms a Markov chain, whose distribution after convergence is the posterior marginal distribution of the random variable. The first draws will not be from the target distribution: the Markov chain needs a burn-in phase after which its distribution is deemed converged. The length of the burn-in phase has to be chosen considering the convergence rate of the Markov chains.

In this application, the convergence has been inspected through the Gelman & Rubin method [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF][START_REF] Bolstad | Understanding computational Bayesian statistics[END_REF]. Several chains are created starting from different initialization points. Then the inter-chain and the inchain variances are compared and the information is condensed in a convergence indicator.

For any parameter, a random sample from its prior law can be chosen as initialization point: the Markov Chain Monte Carlo methods guarantee the convergence towards the posterior marginal distribution for any starting point given a large enough number of iterations. However, it could be useful to optimize the convergence in order to reduce the number of iterations and the calculation time. One way of doing this is to initialize the parameters with the results from the application of a classical Wiener filter. As a consequence, the initialization points will already be near the main mode of the posterior marginal distributions, easing the convergence towards it. An example of the application of the Gelman & Rubin method starting from initialization points scattered around the Wiener filter results is produced in the following section.

It is worth mentioning that, as it can be seen in the algorithm steps, the more the observations the less the prior laws influence the target distributions.

Interpretation of the Gibbs filter as a linear combination of Wiener filters

As already stated in sub-section 2.1.1, one of the issues in using Wiener filters is the need to choose between h k (f ) and h c k (f ) without knowing before hand which one estimates best the signal of interest.

In the case of one source and one reference, it can be easily proven that the mean of the filter µ h after convergence (i.e. the expectation of the posterior law of the filter) is a linear combination of the Wiener filters h(f ) and h c (f ). In what follows each term depends on the frequency, but for the sake of simplicity this dependence is left out. Just note that the analysis proposed can be done at each frequency. After development (Appendix C), the following expression can be obtained:

µ h = k h c h c + k h h = αCh c + (1 -α)h, (10)
α = d0 d0+x+ σ d / √ I obs C = I obs i=1 |z c i | 2 I obs i=1 |zi| 2 . ( 11
)
In the former equations x is the deterministic part of the source contribution, while z c i stands for the random part of the reference. The term σ d / √ I obs , with a standard Gaussian random variable, stands for the measurement noise. Figures 10 and11 help the interpretation of the obtained coefficients:

• if the whole reference signal z i is large with respect to its random part z c i , then C tends to zero and the filter obtained through the proposed method tends to h. This corresponds well to the Wiener filtering strategy of favouring the h filter when the random part of the references is too small: removing the synchronous average from the references would lead to a small SNR and a poor h c filter estimation. It is interesting to point out that, even if C is null, the coefficient k h is not necessarily equal to 1. This means that k h can be considered as a corrective term, thus according to the remark of Leclère who noted that the filter h is always bigger than h c and it can thus be thought to over-estimate the contribution [START_REF] Leclère | Source separation in diesel engines using Wiener filtering: physical interpretations[END_REF];

• if the reference signal z i consists mainly in its random part z c i then C tends to 1 and the h c filter is acceptable. Normally, in this case, if the reference has been wisely chosen, the aimed contribution too consists mainly of its random part, then α tends towards 1 and the h c filter is preferred. This also is consistent with the strategy based on Wiener filter;

• looking more thoroughly at the coefficient α, for an easier understanding it can be rewritten as:

α = 1 1 + x/d 0 + √ I obs σ d d0 . ( 12
)
If the deterministic part of the signal noise (d 0 ) is large with respect to the variance of the noise σ d and to the deterministic part of the contribution x, then α tends to 1 and the h c filter is preferred (Fig. 10). This means that if d 0 is important, then the filter should be computed from the random parts of the signals, which is in accord with the Wiener filter strategy. The same way, increasing the number of observations I obs yields to an increase in the coefficient α tipping the scales in favour of h c (Fig. 11);

• on the other hand, the coefficient α tends to zero if x d0 and/or σ d d0 are large, yielding more weight to the h filter (Fig. 10). This is the case when the deterministic part of the source contribution composes most of the measured signal d i and it is thus bigger than the deterministic part of the noise (d 0 ). It holds too when the variance of the zero mean Gaussian additive noise is far bigger than its average d 0 (which represents the deterministic part of the noise), meaning that the noise is essentially random and not correlated with the reference signal and contribution. In both cases, there is no advantage in computing the h c filter and the use of the h Wiener filter is preferred, which is in accord with what is pointed out in this analysis.

z c i z i ).
The former analysis shows that the proposed Gibbs filter can be considered as a linear combination of the two Wiener filters previously introduced in sub-section 2.1.1. Whereas the adoption of a Wiener filtering strategy yields to the choice of one filter or the other, the proposed method blends both a function of the characteristics of the processed signals.

z c i z i ).

Experimental results on simulated signals

In this section, the algorithm explained in section 3.2 has been applied to signals simulating an engine outputs. The aim is to verify the algorithm and compare it to the Wiener filter separation. As a consequence, the application to a global noise whose components are known is compulsory. In order to do so, synthesized signals which simulate the references and the total noise are obtained through the application of Eqs. ( 3) and ( 4) to known latent variables z, spectrofilters h and matrix G. These, jointly with their respective variances, are shaped after the results of the algorithm being previously applied to actual engine signals. The data d i and r i are composed of the Fourier coefficients of respectively the global noise and the measured references corresponding to one single frequency and to the i = 1, ..., I obs observations. In what follows, I obs = 400 engine cycles, the measured references are N r = 8 and the latent ones are supposed to be N z = 3. The measured references have been chosen so that matrix G has the following form:

G =             G 11 0 0 0 G 22 0 G 31 G 32 0 G 41 G 42 0 G 51 G 52 0 G 61 0 G 63 G 71 0 G 73 G 81 0 G 83             . ( 13
)
Each measured reference, i.e. each term in the vector r i has a dedicated variance σ 2 rj whose prior value is fixed considering the nature of the measured quantity (e.g. microphone or accelerometer). The results presented in what follows concern the convergence of the algorithm and its comparison to the separation based on the Wiener filter method. Figure 12 shows the convergence of σ 2 d over 10 4 iterations of the Markov chain. The convergence is rather fast and the zoom in Fig. 13 shows how the chains, whose initialization point is not the same, converge toward the same distribution. Given a parameter, the initialization points for the chains are scattered around the result returned by the classical Wiener filter: one initialization point corresponds to the "Wiener" solution, another to a tenth of its value and the remaining chains are initialized randomly drawing from a standard uniform distribution defined between 0 and twice the Wiener estimation of the parameter. The number of chains used is N c = 4 and the Gelman & Rubin convergence indicator R is equal to 1.0001. Knowing that the asymptotic value of this indicator for the number of iterations which tends to infinite is 1 and knowing that values inferior to 1.1 can be considered acceptable, it can be concluded that the chains have well converged. Figure 14 shows the convergence of the three contributions separated from the global engine noise. As it can be seen, the convergence is not toward a value but rather toward a distribution, with its mean and variance. Figure 15 makes the comparison between the relative error of the Wiener filter, based on a mean-square optimization, and the one of the present method. For each observation i, the mean of the Markov chains from the Gibbs sampling and the value result of the Wiener filtering are normalized by the theoretical contribution corresponding to the measure d i . It is interesting to remark that the algorithm does not exactly converge toward the theoretical value. This is expected, since the measure d i is itself just a realization of a stochastic process, so there is no reason for it to correspond to the mean of the Markov chain, whose distribution describes entirely the stochastic process. However, the 95% credible interval allows to verify that this measure is usually included in it. In order to compare the results of the proposed method to the ones presented in [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF], an indicator e k is introduced. For a contribution, it is obtained as the sum along the engine cycles of the absolute value of the differences between the known contribution and its estimation:

Gibbs sampler iterations

e k = I obs i=1 |x ki -h k z ki | I obs i=1 |x ki | , ∀k = 1, ..., N z , (14)
where x ki is the known k th contribution for the engine cycle i. It is important to note that this indicator might be frequency dependent according to the frequency range studied. In this section, one single Fourier coefficient has been synthesized (i.e. one single frequency analysed), so that e k is a scalar value. Table 1 shows the result. As it can be seen, the error of the proposed method is less than the one obtained through the application of the Wiener filter h c . This is expected since the presence of correlated references implies an additional waterfall step in the Wiener filter strategy: this step, as explained in paragraph 2.1.1, may remove from the references information useful for the correct identification of a source. On the other hand, the proposed approach considers the possible correlation between the measured references directly in the formulation and extracts from them some uncorrelated latent references in a nonsupervised way. These latent references are then used for the source separation. Moreover, as Fig. 15 shows, the result from the present method can be always followed by the 95% credible interval which conveys some information on its dispersion. The bold black line is the reference (the measure); in dashed blue is the error of the Wiener filter h c ; in solid red is the error of the proposed algorithm, whose result is followed also by the credibility interval (dashed red). Zoom between the 100 th and 150 th engine cycles in the figures at the right.

Experimental results on measured signals

In this section the proposed algorithm is applied to measured signals issued from a diesel engine in order to extract the contributions to the sound pressure of the combustion, the injection and high pressure pump. The presented results concern the contributions of one of the cylinders of the engine, being the signals windowed around the Top Dead Center (TDC)2 of the considered cylinder. The number of observed engine cycles is I obs = 201 and, after having applied the Fourier Transform, the evolution of the Fourier coefficients for one frequency is observed along the cycles. Therefore:

• d i is the evolution along engine cycles of one Fourier coefficient of the overall engine noise measured by a microphone;

• among the measured references r i , the first one (i = 1) derives from the cylinder pressure, the second one (i = 2) from an electric information on the injectors and the others (i = 3, . . . , 8) from accelerometers placed on the engine (Fig. 16). In particular, the combustion is considered as the dominant source and as such it is "visible" in all the accelerometer signals. Moreover, each accelerometer is sensitive to the mechanical part on which it is placed. As a consequence, the accelerometers placed on the injectors (i = 3, . . . , 5) are sensitive to the combustion and the injection, whereas the ones placed on the high pressure pump (i = 6, . . . , 8) are sensitive to the combustion and the pump movement. Thus, the assumed topology of the matrix G is:

G =             G 11 0 0 0 G 22 0 G 31 G 32 0 G 41 G 42 0 G 51 G 52 0 G 61 0 G 63 G 71 0 G 73 G 81 0 G 83             ; (15)
• such a topology yields to z i (i = 1, 2, 3) components being virtual pure references for the combustion, the injection and the high pressure pump, respectively. As it has been done for the synthesized signals, each measured reference has been assigned a different variance σ 2 ri whose prior value is chosen taking into consideration the physical nature of the signal (heteroskedastic mode). The convergence is again rather fast. It is checked through the Gelman & Rubin method using four parallel Markov chains. For each random variable, the initial values of the chains are randomly scattered around the results from a classical Wiener filter separation realized beforehand. The convergence indicator after 10 4 iterations is R = 1.0001, therefore the chains can be considered to have converged. As opposed to the application on synthesized signals (section 4), in this case the separation results cannot be compared to a reference solution, but only compared to the Wiener filtering results. What is most interesting to analyse is the combustion contribution. Its reference being very pure, the Wiener filtering is expected to yield a good estimation of the contribution. The Gibbs extracted contribution is seen to converge towards it (Fig. 18(a)). What is more, the purity of the reference is clearly visible in the quality of the estimation, as hinted by the narrow 95% credible interval. As opposed to the cylinder pressure, the other references are more polluted and this is testified by the larger credibility intervals of the estimations (Fig. 18(b-c)). The proposed algorithm has been thought for just such cases: combining more polluted references to find a virtually pure reference and proposing an estimator jointly with its credible interval giving some information on the accuracy of the task achieved. From what has been said in the previous paragraph it follows that the Gibbs combustion extraction yields similar results to the Wiener extraction. After having applied the same algorithm to the whole frequency range of the measured signals, it is interesting to note whether the Gibbs filter has converged toward h or toward h c (see discussion in subsection 3.3). It can be remarked that the Gibbs result converges generally toward the filter h c (Fig. 19). According to Leclère [START_REF] Leclère | Source separation in diesel engines using Wiener filtering: physical interpretations[END_REF] and as explained in section 2.1.1, indeed the h c Wiener filter boasts a smaller bias error than the h one, but it is really more convenient only if the reference random part is important enough to yield an acceptable SNR and thus a small random error. In the present application it seems to be the case from f 0.046 onward. On the other hand, in the low frequency range the random part of the reference is too weak and the estimation error of the h c filter is too important for it to be acceptable. It can be seen that, in such cases the Gibbs estimator converges toward the Wiener filter h (Fig. 20). The credibility intervals prove that the two filters are significantly different.

Conclusion

A hierarchical Bayesian separation method has been proposed and validated on synthesized and measured signals. The main advantages of this approach lay in (1) achieving the separation task when the available references are correlated and noisy, (2) using a Bayesian framework which yields to the consideration of the prior knowledge, (3) accompanying the results with credible intervals which testify for the quality of the results, (4) allowing a non-supervised choice of the filter as a function of the signals characteristics. Using this approach, some points need particular care. First of all, the more engine cycles are observed, the better the results will be, since the influence of the prior laws will be reduced. However, this yields to long signal acquisitions during which the cyclostationarity and cycloergodicity hypotheses have to be verified all along. In the engine case, this means that one stabilized operating condition is analysed in a run of the algorithm. Moreover, the a priori given shape of the matrix G imposes constraints on the measured references. Here also, expertise and prior knowledge are compulsory in order to safely assume the source contribution in each reference. Finally, depending on the analysed signals, the posterior probability distributions could be multi-modal. In this case, more advanced sam-pling techniques could be used in order to improve the chain mixing properties [START_REF] Ahmadian | Efficient Markov Chain Monte Carlo Methods for Decoding Neural Spike Trains[END_REF]. The proposed approach proved itself worthy of interest for the source separation in the particularly difficult case of IC engines, which are characterized by multiple sources overlapped in the time and the frequency domains. Nevertheless, the approach is general and it can be used in all application where correlated references are available.

a Multivariate Gaussian distribution with parameters:

µ zi = Ω zi (h * (d i -d 0 )/σ 2 d + G H   σ -2 r 1 0 0 0 . . . 0 0 0 σ -2 r Nr   r i ) Ω zi = (h * h t /σ 2 d + G H   σ -2 r 1 0 0 0 . . . 0 0 0 σ -2 r Nr   G + I) -1 . (B.4) 2. Posterior law of σ 2 d : Likelihood : [d i |h, z i , d 0 , σ 2 d ] ∼ N di (h t z i + d 0 , σ 2 d ) Prior : [σ 2 d ] ∼ InvGamma(a d , b d ) (B.5) [σ 2 d |rest] ∝ I obs i=1 [d i |h, z i , d 0 , σ 2 d ][σ 2 d ] ∝ I obs i=1 N di (h t z i + d 0 , σ 2 d ) • InvGamma(a d , b d ) ∝ 1 (σ 2 d ) Iobs exp - 1 σ 2 d I obs i=1 |d i -d 0 -h t z i | 2 1 (σ 2 d ) a d +1 exp - b d σ 2 d ∝ 1 (σ 2 d ) I obs +a d +1 exp - 1 σ 2 d ( I obs i=1 |d i -d 0 -h t z i | 2 + b d ) (B.
I obs i=1 |d i -d 0 -h t z i | 2 + b d . 3. Posterior law of σ 2 rj , ∀j = 1, . . . , N r : Likelihood : [r ji |g j , z i , σ 2 rj ] ∼ N rji (g j z i , σ 2 rj ) Prior : [σ 2 rj ] ∼ InvGamma(a r , b r ) (B.7) [σ 2 rj |rest] ∝ I obs i=1 [r ji |g j , z i , σ 2 rj ][σ 2 rj ] ∝ I obs i=1 N rji (g j z i , σ 2 rj ) • InvGamma(a r , b r ) ∝ 1 (σ 2 rj ) Iobs exp - 1 σ 2 rj I obs i=1 |r ji -g j z i | 2 1 (σ 2 rj ) ar+1 exp - b r σ 2 rj ∝ 1 (σ 2 rj ) I obs +ar+1 exp - 1 σ 2 rj ( I obs i=1 |r ji -g j z i | 2 + b r ) (B.8)
This corresponds to an Inverse-Gamma distribution with parameters a r = I obs + a r and b r =

I obs i=1 |r ji -g j z i | 2 + b r . 4. Posterior law of σ 2 h k , ∀k = 1, . . . , N z : Likelihood : [h k |σ 2 h k ] ∼ N h k (0, σ 2 h k ) Prior : [σ 2 h k ] ∼ InvGamma(a h k , b h k ) (B.9) [σ 2 h k |rest] ∝ [h k |σ 2 h k ][σ 2 h k ] ∝ N h k (0, σ 2 h k ) • InvGamma(a h k , b h k ) ∝ 1 σ 2 h k exp - 1 σ 2 h k |h k | 2 1 (σ 2 h k ) a h k +1 exp - b h k σ 2 h k ∝ 1 (σ 2 h k ) 1+a h k +1 exp - 1 σ 2 h k (|h k | 2 + b h k ) (B.10)
This corresponds to an Inverse-Gamma distribution with parameters

a h k = a h k + 1 and b h k = |h k | 2 + b h k . 5. Posterior law of σ 2 G jk , ∀j = 1, . . . , N r ∀k = 1, . . . , N z : Likelihood : [G jk |σ 2 G jk ] ∼ N G jk (0, σ 2 G jk ) Prior : [σ 2 G jk ] ∼ InvGamma(a G jk , b G jk ) (B.11) [σ 2 G jk |rest] ∝ [G jk |σ 2 G jk ][σ 2 G jk ] ∝ N G jk (0, σ 2 G jk ) • InvGamma(a G jk , b G jk ) ∝ 1 σ 2 G jk exp - 1 σ 2 G jk |G jk | 2 1 (σ 2 G jk ) a G jk +1 exp - b G jk σ 2 G jk ∝ 1 (σ 2 G jk ) 1+a G jk +1 exp - 1 σ 2 G jk (|G jk | 2 + b G jk ) (B.12) This corresponds to an Inverse-Gamma distribution with param- eters a G jk = a G jk + 1 and b G jk = |G jk | 2 + b G jk . 6. Posterior law of h t is: Likelihood : [d i |h, z i , d 0 , σ 2 d ] ∼ N di (h t z i + d 0 , σ 2 d ) Prior : [h t ] ∼ N h t (0, diag(σ 2 h1 , . . . , σ 2 h Nz )) (B.13) [h t |rest] ∝ I obs i=1 [d i |h, z i , d 0 , σ 2 d ][h t ] ∝ I obs i=1 N di (h t z i + d 0 , σ 2 d ) • N h t (0, diag(σ 2 h1 , . . . , σ 2 h Nz )) ∝ exp - 1 σ 2 d I obs i=1 ((d i -d 0 ) -h t z i )((d i -d 0 ) -h t z i ) H • exp -h t diag( 1 σ 2 h1 , . . . , 1 σ 2 h Nz )(h t ) H ∝ exp -h t I obs i=1 z i z H i σ 2 d + diag( 1 σ 2 h1 , . . . , 1 σ 2 h Nz ) (h t ) H • exp 1 σ 2 d I obs i=1 (d i -d 0 )z H i (h t ) H + h t σ 2 d I obs i=1 z i (d i -d 0 ) * (B .14)
A Multivariate Gaussian distribution for a random variable h t with mean m and covariance matrix S would lead to:

[h t |m, S] ∝ exp -h t S -1 (h t ) H + h t S -1 m H + mS -1 (h t ) H (B.
µ h = I obs i=1 (d i -d 0 )z H i Ω h /σ 2 d Ω h = ( I obs i=1 z i z H i /σ 2 d +    σ -2 h 1 0 0 0 . . . 0 0 0 σ -2 h Nz   ) -1 . (B.16) 7.
Posterior law of g j , ∀j = 1, . . . , N r is:

Likelihood : [r ji |g j , z i , σ 2 rj ] ∼ N rji (g j z i , σ 2 rj ) Prior : [g j ] ∼ N gj (0, diag(σ 2 Gj1 , . . . , σ 2 G jNz )) (B.17) 36 [g j |rest] ∝ I obs i=1 [r ji |g j , z i , σ 2 rj ][g j ] ∝ I obs i=1 N rji (g j z i , σ 2 rj ) • N gj (0, diag(σ 2 Gj1 , . . . , σ 2 G jNz )) ∝ exp - 1 σ 2 rj I obs i=1
(r j i -g j z i )(r j i -g j z i )

C Expressing the Gibbs filter as a combination of Wiener filters

For the sake of clarity the development corresponds to the case where N z is equal to one. Moreover, in order to compare the Bayesian filter to the Wiener filters, the variances σ h k are reduced to zero to simulate the deterministic result characteristic of the Wiener approach. Under these hypothesis, the variance of the Gibbs estimate for the filter becomes

ω h = σ 2 d I obs i=1 |z i | 2 (C.1)
and the expected value of the filter can thus be written as:

µ h = I obs i=1 (d i -d 0 )z * i I obs i=1 |z i | 2 (C.2)
where • * stands for the complex conjugate. Some manipulation is necessary to make the Wiener filters appear. To this aim, their formulation is reminded here after: where x i is the contribution of the source, x its deterministic part and the mean of the random noise I obs i=1 n i follows a complex Gaussian distribution defined by σ d / √ I obs with ∼ N C (0, 1) being a standard circular complex Gaussian random variable. This yields equations [START_REF] Hirano | Using multiple regression analysis to estimate the contribution of engine-radiated noise components[END_REF] and [START_REF] Antoni | Novel source separation method for measuring combustion noise and structural loss[END_REF].

Figure 1 :

 1 Figure 1: Explanatory diagram of the signal model as defined by equations (1) and (2).

Figure 2 :Figure 3 :

 23 Figure 2: Signals provided by a cylinder pressure probe (solid red line) and by an accelerometer (dashed blue line) on the corresponding injector.

Figure 4 :

 4 Figure 4: Explanatory diagram of the Wiener filter h k (t).

Figure 5 :Figure 6 :

 56 Figure 5: Spectra of the random (dashed red line) and deterministic part (solid blue line) of a cylinder pressure signal.

Figure 9 :

 9 Figure 9: Hierarchical Bayesian model: data (boxes), parameters (ellipses), hyper-parameters (diamonds).

Figure 10 :

 10 Figure 10: Sensitivity of k h and k h c to x d 0 (a) and σ d d 0 (b) for several values of the reference random part (ratios

Figure 11 :

 11 Figure 11: Sensitivity of k h and k h c to the observation number I obs for values of the reference random part (ratios

Figure 13 :

 13 Figure 13: Convergence of a Markov chain toward the distribution describing σ 2 d (zoom on the burn-in phase).

Figure 14 :

 14 Figure 14: Convergence of a Markov chain toward the distribution of the contributions xki = h k z ki , ∀k = 1, ..., N z and i = {200, 400}. The convergence of the first, second and third source are shown respectively in (a-b),(c-d), (e-f). This for their values at the engine cycles i =200 (a-c-e) and i = 400 (b-d-f). In each figure, the bold solid line is the target result, the bold dashed line is the result obtained by Wiener filter h c and the thin coloured lines are four Markov chains.

Figure 15 :

 15 Figure 15: Estimation error of the algorithms ∀i = 1, ..., I obs . The first source is in (a), the second in (b) and the third in (c). The bold black line is the reference (the measure); in dashed blue is the error of the Wiener filter h c ; in solid red is the error of the proposed algorithm, whose result is followed also by the credibility interval (dashed red). Zoom between the 100 th and 150 th engine cycles in the figures at the right.

Figure 16 :

 16 Figure 16: Sketch of the measurement protocol for one cylinder.

Figure 17 :

 17 Figure 17: Convergence of a Markov chain toward the distribution describing σ 2 d (zoom on the burn-in phase).

Figure 18 :

 18 Figure 18: Comparison between the Wiener estimation based on h c (blue lines) and the Bayesian estimation (red solid lines) of the contributions h k z ki , ∀k = 1, ..., N z . (a), (b) and (c) show respectively the first, second and third contributions. The Bayesian result is followed by the 95% credibility interval (red dashed lines).

Figure 19 :

 19 Figure 19: Convergence of the Gibbs combustion filter (red solid line) towards the h c Wiener filter (green dashed line -plus markers) in a high frequency interval, for the real part (a) and the imaginary part (b). It is shown that the h Wiener filter (blue dashed line -round markers) is out of the 95% credibility interval of the Gibbs filter (red dashed lines). The frequency scale is normalized by the sampling frequency.

Figure 20 :

 20 Figure 20: Convergence of the Gibbs combustion filter (red solid line) towards h Wiener filter (blue dashed line -round markers) in a low frequency interval, for the real part (a) and the imaginary part (b). It is shown that the h c Wiener filter (green dashed line -plus markers) is generally out of the 95% credibility interval of the Gibbs filter (red dashed lines). The frequency scale is normalized by the sampling frequency.

6 )

 6 This corresponds to an Inverse-Gamma distribution with parameters a d = I obs + a d and b d =

8 .| 2 σ 2 d0) 1+a d 0 +1 exp - 1 σ 2 d0(

 822 Posterior law of d 0 is:Likelihood : [d i |h, z i , d 0 , σ 2 d ] ∼ N di (h t z i + d 0 , σ 2 d ) Prior : [d 0 ] ∼ N d0 (0, σ d0 ) (B.21) [d 0 |rest] ∝ I obs i=1 [d i |h, z i , d 0 , σ 2 d ][d 0 ] ∝ I obs i=1 N di (h t z i + d 0 , σ 2 d ) • N d0 (0, σ d0 ) i -h t z i ) -d 0 ) H ((d i -h t z i ) -d 0 ) • exp -|d 0 -h t z i ) (B.22)A Gaussian distribution for a random variable d 0 with mean m and variance s 2 would lead to:[d 0 |m, s 2 ] ∝ exp -|d 0 .22) and (B.23) have the same form and by identification it can be recognised that the posterior distribution of d 0 is a Gaussian distribution with parameters:µ d0 = ν d0 I obs i=1 (d i -h t z i )/σ 2 d ν d0 = (σ -2 d0 + I obs σ -2 d ) -1 . (B.24) 9. Posterior law of σ 2 d0 : Likelihood : [d 0 |σ 2 d0 ] ∼ N d0 (0, σ 2 d0 ) Prior : [σ 2 d0 ] ∼ InvGamma(a d0 , b d0 ) (B.25) [σ 2 d0 |rest] ∝ [d 0 |σ 2 d0 ][σ 2 d0 ] ∝ N d0 (0, σ 2 d0 ) • InvGamma(a d0 , b d0 ) |d 0 | 2 + b d0 ) (B.26)This corresponds to an Inverse-Gamma distribution with parameters a d0 = a d0 + 1 and b d0 = |d 0 | 2 + b d0 .

2 ( 2 ( 2 ( 2 ( 2 (i=1 |z i -z| 2 I obs i=1 |z i | 2 I|zi| 2 ,

 22222222 i -d)(z i -z) * I obs i=1 |z i -z| 2 = I obs i=1 (d i -d)z * i I obs i=1 |z i -z| 2 (C.4)where the synchronous average (deterministic part) of the terms is noted •. Equation (C.2) can then be manipulated as follows:µ h = I obs i=1 (d i -d 0 )z * i I obs i=1 |z i | d i -d) -d 0 + d)z * i I obs i=1 |z i | α + (1 -α))(d i -d) -d 0 + d)z * i I obs i=1 |z i | i -d)z * i I obs i=1 |z i | 2 + I obs i=1 ((1 -α)(d i -d) -d 0 + d)z * i I obs i=1 |z i | i=1 (d i -d)z * i I obs i=1 |z i | 2 + (1 -α) I obs i=1 d i z * i I obs i=1 |z i | 2 + I obs i=1 (αd -d 0 )z * i I obs i=1 |z i | obs i=1 (d id)z * i I obs i=1 |z i -z| 2 + (1 -α) I obs i=1 d i z * i I obs i=1 |z i | 2 + I obs i=1 (αd -d 0 )z * theWiener filter h c can be recognised in the first term, whereas the h filter in the second. Therefore, equation (C.2) can be rewritten as:µ h = αCh c + (1 -α)h + I obs i=1 (αd -d 0 )z * i I obs i=1 |z i | 2 (C.11)Finally, the arbitrary value α can be chosen in order to eliminate the third term: 0 + x + σ d / √ I obs (C.12)

Table 1 :

 1 Comparison between Wiener and Gibbs errors.

	source contribution Wiener Gibbs
	1	0.300	0.227
	2	0.624	0.077
	3	0.106	0.082

  [START_REF] Antoni | El Badaoui Cyclostationary modelling of rotating machine vibration signals[END_REF] Equations (B.14) and (B.15) have the same form and by identification it can be recognised that the posterior distribution of h t is a Multivariate Gaussian distribution with parameters:

  H • exp -g j diag( 1 σ 2 Multivariate Gaussian distribution for a random row vector g with mean m and covariance matrix S would lead to:[g|m, S] ∝ exp -gS -1 g H + gS -1 m H + mS -1 g H (B.19)Equations (B.18) and (B.[START_REF] Huang | MCMC-Particle-based group tracking of space objects within Bayesian framework[END_REF]) have the same form and by identification it can be recognised that the posterior distribution of g j is a Multivariate Gaussian distribution with parameters:

													Gj1	, . . . ,	1 G jNz σ 2	)g H j
	∝ exp -g j (	I obs i=1	z i z H i σ 2 rj	+ diag(	1 σ 2 Gj1	, . . . ,	1 G jNz σ 2	))g H j +	1 σ 2 rj	I obs i=1	r ji z H i g H j +	g j rj σ 2	I obs i=1	z i r * ji
										(B.18)
	A µ gj =	I obs	r ji z H i Ω gj /σ 2
					i=1

rj

Ω gj = (

In Bayesian probability, a prior is called a conjugate prior for a certain likelihood, if prior and posterior are in the same distribution family. For instance, for a Gaussian likelihood, an Inverse Gamma prior for the variance yields an Inverse Gamma posterior as well.

In an internal combustion engine, the TDC denotes the position of a piston in which it is furthest from the crankshaft.

Acknowledgments

This work was funded by Groupe PSA and performed within the framework of the OpenLab Vibro-Acoustic-Tribology@Lyon and the LabEx CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007).

A Usual probability density functions

Multivariate real normal distribution

The multivariate normal distribution of a N -dimensional random vector x ∈ R N ×1 of mean µ and covariance matrix Σ is noted:

In the non-degenerate case (i.e. Σ is definite positive), the PDF of the multivariate real normal law is written:

where det • stand for the determinant of a matrix.

Inverse-Gamma distribution

The Inverse-Gamma distribution can be parametrized using the shape parameter α and the rate parameter β. The Inverse-Gamma distribution of a random variable x is noted:

The corresponding PDF is written as:

Multivariate complex normal distribution

If the complex normal random vector x is proper, than the PDF is written:

) where µ x and Σ xx are respectively the mean and the covariance of x.

B Posterior laws calculation

In what follows, the details about the computation of the Bayesian posterior distributions are proposed. The basis of the posterior computation is the Bayes theorem as presented in [START_REF] Sauvage | Diesel knock noise from combustion phenomenon to perceived signals[END_REF] and [START_REF] Sharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF]. The posteriors are obtained as follows:

1. Posterior law of z i , ∀i = 1, . . . , I obs is:

A Multivariate Gaussian distribution for a random variable z with mean m and covariance matrix S would lead to: