Julie Drouet

Quentin Leclère
email: quentin.leclere@insa-lyon.fr.

Etienne Parizet

Experimental modeling of Wiener filters estimated on an operating diesel engine

Keywords: Wiener filter, diesel engine combustion noise, source separation, modal estimation

Sound source separation in diesel engines can be implemented using a Wiener filter, or spectrofilter, that can extract the combustion contribution in the overall noise. In this study this filter characterizes the transfer function between a cylinder pressure and a measurement point. An engine is characterized by several filters (one for each cylinder) which are estimated for many operating conditions (engine speed and load). A purpose of this work is to obtain an averaged spectrofilter allowing the synthesis of combustion noise in all operating conditions. This synthesis should be accurate enough to be used in perceptive studies. In order to refine the spectrofilter estimation in the medium frequency band, this paper consists in taking advantage of the multitude of information given by the estimations from different operating conditions. To do this, an experimental model is adopted so modal parameters are extracted from the great number of measured filters. Different procedures as the ESPRIT method or the LSCE method (modal analysis) are used to decompose the impulse responses on a complex exponential basis. The spectrofilters estimated from different operating conditions are analyzed and compared in this reduced basis, in order to identify the underlying structural parameters. These parameters are compared to the results of an experimental characterization of the stopped engine. The accuracy of the synthesis (number of components of the filter) is an important issue because these filters will be used in perceptive applications, extracting combustion noises.

Introduction

Noise from diesel engines is a major source of annoyance for car passengers and people outside the vehicle. This has led to many studies aiming at reducing this noise, which is partly due to the high pressure combustion in the cylinders, but also to mechanical sources involved in the engine operation (distribution, belts, impacts due to clearances in mechanical joints). It is therefore very important to identify the contribution of potential sources, in order to implement efficient noise reduction solutions. Several studies are dedicated to this issue. Firsts results are obtained from the comparison between different operating conditions: it is shown in [START_REF] Priede | In search of origins of engine noise -an historical review[END_REF] how mechanical and combustion contributions to the overall noise are varying with operation. Under a certain combustion pressure level, called the critical cylinder pressure level, the noise level does not vary: it means that the mechanical noise dominates the engine noise. Above this limit, the noise level increase linearly in function of the combustion pressure level: the engine noise is thus mainly due to combustion. Another approach, studied in [START_REF] Leipold | Noise, emissions and performance of the Diesel engine-a comparison between di and idi combustion systems[END_REF], is based on the comparison between the engine operating by itself and driven with an eletric motor. The difference between noise levels are considered as the contribution of the combustion. This study shows the benefits of the indirect injection technology over the standard direct injection, concerning the combustion noise reduction. Of course these methods are based on the rough assumption that the mechanical sources are not affected by changes in the loading conditions. Several papers are studying the possibility to use signal processing tools to separate the contribution of combustion in the overall engine noise [START_REF] Chung | Measurement of frequency reponses and the multiple coherence function of the noise-generation system of a Diesel engine[END_REF][START_REF] Alfredson | The partial coherence technique for source identification on a Diesel engine[END_REF][START_REF] Yawata | Identification of internal noise sources in Diesel engines[END_REF]. The idea is to use the coherence between the acoustic pressure outside the engine and cylinder pressure signals, measured using pressure sensors introduced in the cylinders through the cylinder head. However, there is a strong difficulty due to the fact that mechanical excitations are strongly coherent with the cylinder pressure signals, making it difficult to obtain a valuable result. However, this difficulty can be overcome by introducing some randomness in the injection, as shown in [START_REF] Albright | Conditioned source analysis, a technique for multiple input system identification with application to combustion energy separation in piston engines[END_REF], but this approach is somewhat intrusive, and not without effect on the engine operation. The characterization of the transfer between the cylinder pressure and the acoustic pressure outside the engine is not new. Pioneer studies [START_REF] Austen | Origins of diesel engin noise[END_REF] are defining the structural attenuation as the difference between engine noise and cylinder pressure spectra, in decibels (dB). It is shown that this structural attenuation is relatively constant when changing operating conditions, and that structural attenuations of different engines are quite similar. The outcome of these works has been the commercialization of a measurement system called the Combustion Noise Meter [START_REF] Reinhart | An evaluation of the lucas combustion noise meter on cummins 'b' series engines[END_REF], which consists in filtering the cylinder pressure by a typical structural response plus an A-weighting curve, to finally give the contribution of the combustion to the overal noise. But this approach does not use any microphone, and the result only gives a relative indication of the com-bustion noise level, allowing to estimate for instance the effect of some changes in the injection map on the combustion noise. The structural attenuation, as defined here, is a ratio between two power spectral densities, often studied in octave or third octave bands. It gives an estimation of the modulus of the linear filter between the cylinder pressure and the acoustic pressure, averaged in frequency bands. Thus, it does not allow the direct synthesis of the acoustic pressure from cylinder pressures, which is required to consider the separation of mechanical and combustion contributions in the time domain. Several works are reported in the literature, addressing the issue of diesel engines source separation in the time domain [START_REF] Badaoui | 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][START_REF] Wang | Optimization of combustion noise of modern diesel engines for passengers cars[END_REF]. The most convincing studies belong to the family of supervised source separation approaches: the source (cylinder pressures) is measured, and its contribution to the measured output (microphone) is to be estimated. It has been shown quite recently [START_REF] Antoni | Effective vibration analysis of ic engines using cyclostationarity. part ii -new results on the reconstruction of the cylinder pressures[END_REF][START_REF] Pruvost | Extraction du bruit de combustion d'un moteur Diesel[END_REF][START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter. mechanical systems and signal processing[END_REF] how to implement cyclostationarity tools in this context. Suppressing the periodic parts of signals overcomes the difficulty of high correlation between combustion and mechanical excitations [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter. mechanical systems and signal processing[END_REF]. One result of this operation is that filters obtained from random parts of input and output signals are much more stable to changes in operating conditions than filters obtained with raw signals.

The aim of this paper is to show how the filters estimated with these approaches can be used to analyse the structural response of the engine in operation. The originality of this work mainly concerns two aspects. The first one is that the spectrofilters resulting from the source separation algorithm are compared to transfer functions measured with an impact hammer, with the engine stopped. It is shown how this comparison can be used to analyze the transmission paths of the combustion noise. The second original contribution is the application of modal analysis approaches (LSCE [START_REF] Brown | Parameter estimation techniques for modal analysis[END_REF] and ESPRIT [START_REF] Roy | Esprit -estimation of signal parameters via rotational invariance techniques[END_REF][START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoïdes modulées[END_REF]), aiming to compare filters obtained at different operating conditions in a modal framework. The first part of this work is dedicated to some general considerations about engine noise sources, and to the transmission paths of the combustion noise. The second part presents spectrofilters obtained from measurements in operating conditions. These spectrofilters are compared to standard impact hammer FRF measurements on the stopped engine. Then, the basic theory of ESPRIT and LSCE is exposed in the third section, highlighting the similarities and differences between these two methods. They are also compared in practice with an application on a synthetic signal. Finally, in the last section, ESPRIT is carried out on the spectrofilters measured in operation, pointing out the limits of the approach.

Diesel engine noise sources and combustion noise transmission paths

Diesel engines are complex mechanical systems, in which many sources are contributing to the overall noise. There are several ways of classifying engine noise sources. The classification used in this work, shown in Fig. 1, is loosely based on the literature [START_REF] Russel | Diesel engine noise : control at source[END_REF][START_REF] Baxa | Noise Control In Internal Combustion Engines[END_REF]. The internal sources are separated into three categories: the combustion, generating high pressure pulses on internal faces of the combustion chamber, the load-independent mechanical sources, whose behavior is affected only by the rotation speed, and the load-dependent mechanical sources, whose amplitude generally increase with the load. For instance inertia loads and distribution (camshaft, valves) belong to the class of load-independent sources, while the injection system (injectors, injection circuit, high pressure fuel pump) is a load-dependent source. This classification has the advantage to be closely linked to the physical phenomena, but it has also limitations. Some difficulties appear for instance when considering the excitation generated by the piston on the cylinder block, resulting from a combination of the three categories. This work is focusing on the combustion noise, resulting from the high pressure pulses in the combustion chamber generated by the auto-ignition of the fuelair mixture. This pressure pulses are distributed on the combustion chamber walls: the cylinder head, the very top of the cylinder and the piston. The loads directly applied to the cylinder block (the cylinder head and the very top of the cylinder) are responsible of the so-called direct gaz path (see Fig. 2). The loads applied to the piston are transmitted transversely to the cylinders (indirect transverse path) and to the crankshaft bearings through the connecting rod and crankshaft (indirect vertical path). The loads at the crankshaft bearings are mainly vertical loads, but the elasticity of the whole mechanism generates also moments inducing the bending deflection of bearings (see [START_REF] Priede | Characteristics of exciting forces and structural response of turbocharged Diesel engines[END_REF][START_REF] Leclere | Indirect measurement of main bearing loads in an operating diesel engine[END_REF]).

The characterization and ranking of these transmission paths is of prime interest to be able to conceive efficient noise reduction solutions. An interest of this work is to show how the filters obtained from the source separation algorithm can be used by mechanical engineers to get useful information about this issue.

2 Determination of spectrofilters

Experimental setup

Measurements have been realized on a 1.9 L dci diesel engine. Two accelerometers (PBC) are mounted vertically, one inside the engine on the external bearing cap crankshaft ♯1 (distribution side), and the other one on the cylinder head in front of cylinder ♯1 (distribution side too). A microphone (BK) was placed 36 cm above the engine (see Fig. 3).

Firstly, according to the method described in [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter. mechanical systems and signal processing[END_REF], Wiener filters have been estimated on the operating engine, from additional measurement of cylinder pressure and crankshaft angle. Signals were centered by removing the cyclic average, and windowed with a Tukey window (window's length equal to approximately a quarter of engine cycle). The cyclic average is synchronized for each cylinder to the corresponding TDC, using the angle encoder signal. The total length of each recording is 30seconds, that is to say about 200 and 350 engine cycles at 810 and 1300 rpm. The whole computational details are given in [START_REF] Drouet | Séparation des sources de bruit des moteurs Diesel : Application en hiérarchisation de source et qualité sonore[END_REF]. H1 estimators were then computed between each couple (accelerometer -microphone) / cylinder pressure.

Secondly, to get a set of response functions determined by conventional methods an impact hammer study has been conducted on the (still hot) engine with the same response transducers. The impact was realized vertically on the cylinder head, just above the cylinder ♯1 (see Fig. 3). Even if the combustion excitation could not be reproduced, the structure was the same, and could emphasize similar modal parameters.

Analysis of measured response functions

Combustion excitations are equivalent to spatially distributed impacts on the cylinder head and the top of the cylinder (later on called gas excitation path), and on the cylinder and bearing via moving parts (piston-crankshaft path). The hammer excitation point was closed to the piston so that it can be hypothesized that the same paths contributed to the measurements in the stopped conditions. The spectrofilter between the external accelerometer and the cylinder ♯1 is drawn in Fig. 4 in addition to the response function between the same accelerometer and the impact hammer. The modulus of the spectrofilters was divided by the cylinder section, in order to obtain a response function in ms -2 N -1 .

Fig. 4. FRF magnitude for the cylinder head accelerometer, for different operating points (810 rpm (dash-dot blue), 1300rpm with 92Nm (dash red) and 150Nm (solid green)) and for the impact hammer (engine stopped) (solid black)

Identified FRFs (frequency response function) for operating and stopped engine are similar between 1 kHz and 3 kHz, which implies that the cylinder head impact is similar to the combustion excitation, emphasizing the predominance of the gas excitation path (which could be expected for an accelerometer placed on the cylinder head).

Fig. 5. FRF magnitude for the bearing cap accelerometer for different operating points (810 rpm (dash-dot blue), 1300rpm with 92Nm (dash red) and 150Nm (solid green)) and for the impact hammer (engine stopped)(solid black)

For the internal accelerometer placed on the bearing cap (Fig. 5), FRF obtained from the impact strongly underestimates the spectrofilter between 1 and 2 kHz. This can be explained by the deformation of the bearing [START_REF] Chiollaz | Engine noise characterisation with wigner-ville tiemfrequency analysis[END_REF] which is submitted to a high constraint when the engine is running. This is not reproduced by the impact on the cylinder head on stopped engine. The important peak observed around 1.5 kHz on the spectrofilters does not appear on the impact response. This peak can be attributed to the piston-crankshaft path of the combustion noise, and is probably due to a modal behavior of the whole path. This observation agrees with results of the literature [START_REF] Russel | Diesel engine noise : control at source[END_REF].

Fig. 6. FRF magnitude for microphone for different operating points (810 rpm (dash-dot blue), 1300rpm with 92Nm (dash red) and 150Nm (solid green)) and for the impact hammer (engine stopped)(solid black)

For the microphone measurements (Fig. 6), the propagation phenomenon is more complex, because of the acoustic radiation and propagation. However, similar observations can be made : for frequencies up to 1.2 kHz, the impact overestimates the response because of the cancellation of gas-excitation and piston-crankshaft paths happening during operation. Above 1.8 kHz the level of the impact response is similar to the spectrofilters ones, which means that the cancellation effect is no more present and that the overall radiated sound level does not dependent on the excitation location. The same analysis for calculated response functions can be found in the literature [START_REF] Wang | Optimization of combustion noise of modern diesel engines for passengers cars[END_REF].

3 Presentation of the ESPRIT and LSCE methods

Theory

Both methods decompose the measured signal on a complex exponential basis :

s(t) = ℜe K k=1 r k e λ k t (1)
When the signal is sampled, this can be written as:

s n = ℜe K k=1 r k z n-1 k , n ∈ [1...N ] (2)
With

z k = e λ k ∆t (3)
λ k = -δ k + j ω dk = -α k ω uk + j ω uk 1 -α 2 k ( 4
)
δ k is the damping factor, α k the damping ratio, ω uk the undamped natural frequency and ω dk the damped natural frequency of mode k, the complex coefficients r k are called the residues.

The goal of these methods is to identify the poles z k and then compute the residues r k to assess the signal s n . Once the poles are identified the damping ratio and the undamped natural frequency can easily be determined.

ESPRIT is a high resolution modal analysis method [START_REF] Roy | Esprit -estimation of signal parameters via rotational invariance techniques[END_REF][START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoïdes modulées[END_REF][START_REF] Le Carrou | Vibro-acoustique de la harpe de concert[END_REF] based on particular properties of the signal autocorrelation matrix. It is widely used in musical applications to extract modal components with very close frequencies, as developed in [START_REF] Ege | High-resolution modal analysis[END_REF]. The LSCE method is a well known SDOF parameter estimation algorithm [START_REF] Brown | Parameter estimation techniques for modal analysis[END_REF], available in standard modal analysis softwares [26].

Pole identification

The first step consists in identifying poles from a hypothesis on the number of components K. The signal decomposition (1) implies couples (r k , λ k ) and (r * k , λ * k ) represent the same component. Thus searching for 2K complex poles will allow the identification of K different components.

Both methods organize the signal in a Hankel matrix. For the ESPRIT method, the number of lines, denoted as n, depends on the analyzed signal length N . It must be chosen between N/3 and 2N/3 according to [START_REF] Ege | High-resolution modal analysis[END_REF]. In the LSCE method, the number of columns of the Hankel matrix is equal to 2K.

X esprit =        x 1 • • • x N -n-1 . . . . . . . . . x n • • • x N -1        (5)
X lsce =        x 1 • • • x 2K . . . . . . . . . x N -2K • • • x N -1        (6)
There are no other direct similarities between the two methods. ESPRIT uses auto-correlation matrix R XX = 1 l XX H , which is decomposed in eigenvalues. The eigenvectors w i corresponding to the 2K largest eigenvalues are used to construct a Φ matrix. Let W be the matrix composed of the w i vectors in columns. W ↑ is obtained from W by removing it first line and W ↓ by removing the last line of W . Φ is defined as :

Φ = W + ↓ W ↑ (7)
where W + ↓ represents the pseudo-inverse of W ↓ . The eigenvalues of Φ are the poles z k , from which modal parameters are deduced from equations (3) and (4).

The LSCE method is based on the assumption that the studied system can be described by a linear differential equation with constant coefficients of order 2K.

A 2k d dt 2K x(t) + • • • + A 0 x(t) = 0 (8)
For a discretized signal, this is equivalent to :

2K+1 k=1 a k x n+k-1 = 0 (9)
a 2K+1 is arbitrarily fixed to 1 to identify coefficients a k (k ∈ [1 . . . 2K]). Thus each x n (n ∈ [2K + 1 . . . N ]) value is a function of the 2K preceding values. This leads to the follow system :

X lsce        a 1 . . . a 2K        = -        x 2K+1 . . . x N        (10)
This system has to be inverted. The unicity of the solution implies K < 4N , so that the system is over-determined and can be resolved by pseudo-inversion.

       a 1 . . . a 2K        = -X + lsce        x 2K+1 . . . x N        (11)
Each component of discrete signal x n , is written as rz n-1 , and has to be a solution of [START_REF] Reinhart | An evaluation of the lucas combustion noise meter on cummins 'b' series engines[END_REF]. This leads to the polynomial equation [START_REF] Antoni | Effective vibration analysis of ic engines using cyclostationarity. part ii -new results on the reconstruction of the cylinder pressures[END_REF], the roots of which are the poles to be identified.

2K+1 k=1 a k z k-1 n = 0 (12)
Only complex roots with a positive imaginary part are retained.

Identification of the residues

The identification of residues is made by a pseudo-inversion of the system constructed from equation ( 2) :

[r k ] = V + x (13)
where r and x represent the residues and the signal. V is the Vandermonde matrix, determined from the K identified poles.

V =               1 1 • • • 1 z 1 z 2 • • • z K z 2 1 z 2 2 • • • z 2 K . . . . . . . . . . . . z N -1 1 z N -1 2 • • • z N -1 K               (14)
The signal reconstruction s is obtained from the following equation :

s = 2ℜe([r k ]V ) (15)

Application on a synthetic signal

As the modal components of the spectrofilters are unknown, the properties of the methods are evaluated using a synthetic signal, illustrated in Fig. 7.

The undamped natural frequencies of this signal are defined between 500 to 5000 Hz, regularly spaced of 150 Hz, giving a whole number of 31 modes. The damping ratio decreases proportionally to the inverse of the natural frequencies and varies between 0.07 and 0.007. There is no noise added to the synthetic signal. The choice of the number of extracted components K is crucial for both methods. If it is lower (or higher) than the real number of signal components the computation will be normally called under (or over)-estimated. When these numbers are exactly the same (ideal conditions) the methods should both identify the true components of the signal.

The under-estimation of the modes theoretically brings to false identification.

True and identified poles are drawn in Fig. 8 in the complex plane. The angle between 0 and π represents the damped natural frequencies varying between 0 and F s/2. The modulus represents the damping : close to unity when the damping is low (<<1); lower than unity if the damping is high; and greater than unity when the damping is negative (non physical increasing exponential envelope). ESPRIT succeeds to extract twenty poles in the correct range of frequency and damping. Poles identified with LSCE are strongly different from the true ones. This difference can be explained by Fig. 9. ESPRIT identifies poles only in the frequency range of the spectrofilter, whereas LSCE searches poles over the whole frequency range between 0 and F s/2. Fig. 9 is another representation of these observations giving an information about the modal density of filters averaged by 500 Hz frequency bands. When the number of components is over-estimated, the true components are identified by both methods, in addition to other wrong components to satisfy the expected K value. Nonetheless, below 105 required modes LSCE cannot extract all the true components. To understand this phenomenon the extracted poles, in ideal conditions (K = 31) for the two methods, and for K = 105 for LSCE, are plotted in Fig. 10. The necessity of LSCE to be carried out with a strong over-estimation of K is explained by the fact that it finds modes in the overall frequency domain, whereas ESPRIT works only in frequency bands of high energy. The number of modes by 500 Hz frequency bands is drawn in Fig. 11, for true and identified poles. The modal density of ESPRIT is strictly equal to the true one, the poles being exactly recognized. For LSCE, between 500 to 5000 Hz, the modal density is clearly underestimated for K = 31, and wrong modes are determined outside this frequency range. The correct modal density is found for K = 105, i.e 31 in the frequency range 500-5000 Hz, but 74 "computational" modes are also identified in high frequency, almost regularly spaced between 5000 and F s/2. These observations bring out that signals should be resampled in function of the studied frequency range, to avoid LSCE to identify non physical modes in the high frequency range. Because the aim of the study is the information reduction of the spectrofilter, in the following of the paper the components number are voluntary underestimated. Considering these results on synthetic data, the ESPRIT method is preferred to LSCE for the modal analysis of measurements.

Sensitivity to additive noise

In order to assess the robustness of the ESPRIT method, the identification of poles in ideal conditions (K = 31) is carried out with an additive noise of 1 and 10 %, (signal to noise ratios equals to respectively 40 and 20dB). The poles identified by ESPRIT in noisy conditions are shown in Fig. 12. The results for 1% noise are satisfying : all the poles are correctly identified. For a noise level of 10 %, some poles (6 over 31, about 20 %) are identified at a wrong place in the complex plane. These observations illustrate the potential effects of uncertainties on the extraction of the poles. However, the filters obtained in operation are assessed on a large number of signal realizations (engine cycles), and the estimation uncertainties are expected to be relatively controlled.

Modal analysis of measurements

In the following ESPRIT method has been applied to the measured FRFs ( see part 1). Results are given for the Wiener filters relating the cylinder pressure ♯1 and the pressure measured above the engine. The modal density of the structure is rather high (about 800 modes between 500 and 5000 Hz, according to a finite element model presented in [START_REF] Leclere | Etude et dévelopement de la mesure indirecte d'efforts. Application à l'identification des sources internes du moteur[END_REF]). The number of extracted components K was set to 150 to minimize the estimation error of the reconstructed signal, so it is under-estimated as expected.

Fig. 13 shows that the signal reconstruction is correct for frequencies higher than 1.5 kHz. The discrepancy appearing in the low frequency range can be due to a lower energy in this range : modal parameters are mainly identified in high energy areas. This phenomenon can be corrected with an analysis conducted by frequency bands [START_REF] Laroche | The use of the matrix pencil method for spectrum analysis of musical signals[END_REF] or using spectral whitening [START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoïdes modulées[END_REF]. For the spectrofilters a spectral whitening is preferred to a frequency band analysis, because the number of dominant modes in each band cannot be easily estimated. In this study, whitening has been applied between 200 and 5000 Hz, in 600 Hz wide bands. In each band, the filtered signals are weighted by the inverse of their RMS values, and then summed to obtain the whitened spectrofilter. The poles are extracted from this whitened filter, but the signal reconstruction (after the identification of residues) is done on the original signal. For the frequency response obtained with the impact hammer, whitening is not necessary because the dynamic range (about 15dB) is much lower than the one obtained in operation (about 30dB). It can be seen in Fig. 14 than ESPRIT estimation fits well the measured FRF on the whole frequency range. Fig. 15 shows the damping ratios averaged in 500Hz wide frequency bands. These ratios decrease with the frequency. Slopes are similar in all cases, but the estimated damping ratios in operating conditions seem to increase with the rotational speed of the engine. That may be due to the windowing applied during the spectrofilter calculation. For 1300 rpm, the temporal window length, corresponding to about one stroke [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter. mechanical systems and signal processing[END_REF], is short as compared to the impulse response length of the engine block. This windowing leads to an overestimation of the damping. This phenomenon is illustrated in Fig. 16.

Windowing shortens the duration of the Wiener filter impulse response when the rotational speed increases, which means that the damping ratios will be artificially increased as compared to the damping ratio estimated from the impact. The lengths of spectrofilters are well correlated with the lengths of the one-stroke windows, confirming the effect of the windowing on the overestimation of the damping, which is about 30% at idle (810rpm), and 45% at 1300 rpm.

Effects of the spectrofilter estimation errors on the synthesized combustion noise level

In the previous section, it has been shown that the length of the window which is used in the spectrofilter estimation has a significant effect on the identified structural parameters. The damping can be greatly overestimated when the window length is shortened, due to the increase of the engine rotation speed. The consequence is that the levels of the combustion noise synthesized using such filters is expected to be underestimated. An additional computation was conducted in order to confirm those expectations. The combustion noise was synthesized for a high engine speed (2800rpm), using spectrofilters between the cylinder pressures and the microphone point estimated at 3 different speeds : idle (810 rpm), 1800rpm, and 2800rpm. The combustion noise is synthesized simply by convolving 30 seconds of the 4 cylinder pressure signals measured at 2800rpm by the inverse FFT of the different spectrofilters. The mechanical contribution is obtained by subtracting the combustion contribution from the measured acoustic pressure signal. The noise levels in dB(A) are then averaged over 30 seconds. The full processing time is less than one minute using a standard PC. Note that for the last speed, the same recording was used for the spectrofilters estimation and for the synthesis. The results are given in table 1. 1 Noise levels of combustion and mechanical noise at 2800rpm, full load, using spectrofilters extracted at idle (810 rpm), 1800rpm, and 2800rpm.

Filters

First of all, it can be noted that, for this relatively high engine speed, the combustion noise is relatively low as compared to the mechanical noise. Furthermore the total noise level is not equal to the quadratic sum of combustion and mechanical noise levels, both components being partially correlated. The difference between the combustion contribution and the overall noise level decreases when using the spectrofilters obtained at lower speeds. This difference is about 10 dB(A) when using the filters estimated at 2800rpm, but 8.5 and 6 dB(A) when using the filters estimated at 1800 and 810 rpm, respectively. These results confirm what was expected : the overestimation of the filter's damping has a significant effect on the resulting level of combustion noise. However, increase of the combustion noise level does not imply the mechanical noise level decreases: this imply that this level remains approximately un-changed whatever the spectrofilter. A decrease of the mechanical noise could have been expected, because a part of the combustion noise that is falsely identified as mechanical noise with spectrofilters estimated at 2800rpm should have been reassigned to the combustion noise using spectrofilters estimated at lower speeds. These observations could be explained by significant deviations of transfer function with the engine operating point, because of heat changes or non-linearity effects in mechanical joints subjected to intense loads.

Conclusion

Two modal analysis methods (ESPRIT and LSCE) have been used to analyze FRFs of a diesel engine, measured either in operating conditions or on the stopped engine. The FRF obtained with the impact hammer on the cylinder head of the stopped engine represents the gas excitation path only, while in operating conditions the piston-crankshaft path is dominant. A modal behavior of this path is clearly identified between 1 and 2 kHz, inducing a high gain of the spectrofilters in this frequency range. Above 1.8 kHz, the averaged level of FRFs are similar. The impulse response corresponding to these functions are then decomposed in damped sinusoidal components. A case study brings to the conclusion that ESPRIT method is more adapted to spectrofilters analysis than LSCE. A pre-whitening of the Wiener filter is however necessary to enable ESPRIT to reconstruct the signal precisely in the whole frequency range. The damping of identified poles seems to increase with the operating engine speed, and is always significantly higher than the damping identified on the stopped engine. This observation can be explained by the signal processing stage required to estimate the spectrofilter in operating conditions. The signals have to be windowed over time durations corresponding to about a quarter of cycle, i.e. about 40ms to 10ms at 800 to 3000 rpm. These window lengths are too small, as compared to the engine block impulse response, resulting in an overestimation of the damping of the estimated filters. The effect on the extracted combustion noise has been found to be significant, with an underestimation reaching 4dB(A) at high speed.

Fig. 1 .

 1 Fig. 1. Classification of diesel engine internal excitation and radiated noise.

Fig. 2 .

 2 Fig. 2. Combustion noise transmission paths in diesel engines.

Fig. 3 .

 3 Fig. 3. Microphone position (left). External accelerometer and impact point location (right).

Fig. 7 .

 7 Fig. 7. Normalized impulse response of a Wiener filter (dash green) and normalized synthetic signal (solid blue) in the time (top) and frequency (bottom) domains

Fig. 8 .Fig. 9 .

 89 Fig. 8. Poles of the synthetic signal (black circles), poles identified in under-estimation (K = 20) for LSCE (blue x-mark) and ESPRIT (red plus)

Fig. 10 .

 10 Fig. 10. Poles of the synthetic signal (black circles), poles identified in ideal conditions (K = 31) for ESPRIT (red x) and LSCE (blue plus) and in over-estimation (K=105) for LSCE (purple dots)

Fig. 11 .

 11 Fig. 11. Number of modes identified in 500 Hz wide bands for the synthetic signal (solid black) and in ideal conditions (K = 31) for ESPRIT (dot red) and LSCE (dash-dot blue) and in over-estimation (K=105) for LSCE (dash-dot purple)

Fig. 12 .

 12 Fig. 12. Poles of the synthetic signal (black circles), poles identified in ideal conditions (K = 31) for ESPRIT (red x) with an additive noise of 10% (top) and 1% (bottom).

Fig. 13 .

 13 Fig. 13. Microphone / cylinder ♯1 FRF magnitude at 1300 rpm and 92 Nm, for measurement (dash black) and estimated signals for ESPRIT (dot red) and whitened ESPRIT (solid blue)

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Microphone / cylinder ♯1 FRF magnitude at stopped engine, for measurement (solid black) and estimated signals for ESPRIT (dash red)

Fig. 16 .

 16 Fig.[START_REF] Roy | Esprit -estimation of signal parameters via rotational invariance techniques[END_REF]. Normalized impulse responses and used Tukey windows at 1300 rpm (left in solid red) and 810 rpm (right in solid blue)) compared to the impact response (dash-dot black)

Acknowledgements

This work was supported by the Région Rhône Alpes, and performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007).


Quentin Leclère

Nacer Hamzaoui

Using the moving synchronous average to analyse fuzzy cyclostationary signals

Keywords: Fuzzy cyclostationarity, synchronous average, angular sampling, instantaneous rotation speed, diesel engine diagnosis, diesel engine noise

Cyclostationarity is a property of vibration and acoustic signals recorded on rotating machines operating at constant speed. It states that the statistic properties of signals are periodic: the random process defined by the signal observed at a given position in the cycle is stationary, the cycle being defined as the angle interval between two identical configurations of the mechanical system. Cyclostationarity is not fully satisfied if the signals are acquired in the time domain on rotating machines with a fluctuating rotation speed. Indeed, if the instantaneous rotation speed is not purely periodic, it means that time samples taken at a constant time interval (equal to the average cycle duration) do not correspond exactly to an angle in the cycle. In this particular case, a synchronous averaging of cycle realizations can still be processed to estimate a periodic part using a predefined trigger angle to align cycle realizations before the averaging process. In these conditions, the synchronous average depends on the chosen synchronization angle : each point of the synchronous average is an estimate of the expected value of the signal at a given time preceding or following the synchronization angle. The synchronous average can be computed as a function of the synchronization angle, varying over an entire cycle. The result is a moving synchronous average that can be post-processed for diagnosis purposes. For example, a time frequency representation of the moving synchronous average can be computed, and the synchronization angle maximizing each point of the time frequency map can be easily extracted. Under certain conditions of instantaneous speed fluctuations, this analysis allows the precise localization of different mechanical events in the cycle, as well as their contributions in the analyzed vibration or acoustic signal. The analysis of the moving synchronous average leads also to the estimation of the energy loss of the synchronous average processed in angle caused by cycle to cycle speed fluctuations.

s tf m (t, f ) time frequency distribution of s(t)

The synchronous average is a signal processing tool consisting in averaging several portions of signal chosen in function of the studied mechanical system [START_REF] Braun | Analysis of repetitive mechanism signatures[END_REF]. For cyclostationary signals, the synchronous average is an estimation of the deterministic or periodic part [START_REF] Antoni | Cyclostationarity by examples[END_REF]. It can also be regarded through comb filtering [START_REF] Braun | The synchronous (time domain) average revisited[END_REF], the signal is filtered by a specific comb resulting in an averaging process over signal samples separated by a specific delay. Noise and vibration signals recorded on rotating machines are determined by cyclic excitations quentin.leclere@insa-lyon.fr.

defined in angle, and the response of the structure is defined by its time domain impulse response, the determinism of such signals being consequently dual [START_REF] Leclère | Angular and temporal determinism of rotating machine signals: The diesel engine case[END_REF]. If the studied rotating machine is operating with a strictly constant or purely periodic IAS (Instantaneous Angular Speed), the time and angular clocks are totally dependent, and this duality is masked ; processing signals in time or angle will lead to equivalent results. If not, the duality leads to so-called fuzzy-cyclostationary signals [START_REF] Bonnardot | Comparaison entre les analyses angulaire et temporelle des signaux vibratoires de machines tournantes[END_REF], and the result of synchronous averaging will be different in time or angle. The fuzzy-cyclostationary concept is established in the first section of this paper. The time-domain synchronous average is redefined in this context in the second part of this work. It requires especially to define a synchronization angle, and its computation for several synchronization angles in the cycle is defined here as the moving synchronous average. The concepts proposed in this study are illustrated by two experimental examples presented in the third section. The fourth part shows how ASA and TSA (Angular and Temporal Synchronous Averages) can be compared. Then, the moving synchronous average is introduced, and different post processing tools are proposed to extract some information about the studied rotating machine.

A definition of Fuzzy Cyclostationarity

The concept of Fuzzy Cyclostationarity is deeply inspected in reference [START_REF] Bonnardot | Comparaison entre les analyses angulaire et temporelle des signaux vibratoires de machines tournantes[END_REF] (in French). A brief definition is however given here, in order to clearly establish the general framework of this study.

Let ω(t) represent the instantaneous angular speed (IAS, rad/s) of a rotating machine as a function of time, and f the bijection between the temporal and angular domains such as t = f (θ). The instantaneous angular speed can be expressed as a function of the angle :

Ω(θ) = ω(f (θ)) (1)
In stationary operating conditions, Ω(θ) is assumed to be cyclostationary of cycle length Θ. The bijection f can be expressed as a function of Ω :

f (θ) = t = t 0 dτ = θ 0 dα Ω(α) (2)
The statistical properties of Ω(θ) are Θ-periodic. The angle to time relationship f (θ) can be written at θ + Θ as follows :

f (θ + Θ) = θ 0 dα Ω(α) + θ+Θ θ dα Ω(α) = t + T (θ) (3)
The IAS at θ + Θ is thus equal to

Ω(θ + Θ) = ω(f (θ + Θ)) = ω(t + T (θ)) (4)
If T (θ) = T is constant, it means that the cyclostationarity of Ω in angle implies the cyclostationarity of ω in time. According to Eq. ( 3), this is the case if Ω(θ) is constant or periodic. If the cyclic variance of Ω(θ) is not null, then T (θ) depends on θ and some fluctuations of the cycle length will appear in time. In this case, ω(t) keeps the property of cyclostationarity, but with a cycle length equal to the average value of T (θ). It means that cycle realizations in the time domain will not correspond exactly to a cycle length Θ in angle : that's why this kind of cyclostationary process is qualified as fuzzy.

An important consequence of fuzzy cyclostationarity is that cyclic estimators, like the cyclic average or variance, will lead to different results if processed using a constant time or angle sampling. The present work aims to study these differences, and to present how the fuzzy cyclostationarity property can be used for diagnosis purposes.

3 The synchronous average : an estimation of the deterministic part of the signal

Synchronous averaging of cyclostationary signals

Cyclostationary signals are characterized, by definition, by periodical statistic properties. For signals recorded on mechanisms, this periodicity is a mechanical periodicity, i.e. the duration between two identical geometrical configurations of moving parts of a rotating machine (that will be designed in the following by 'cycle'). If the machine is operating with a constant (or purely periodic) rotation speed, the cycle length can be defined in time or angle, and signals can be processed either with a time or an angular sampling. The synchronous average (also called cyclic average) is defined as the average value of the signal at a given position in the cycle. It is in fact an estimation of the expected value of the signal, assessed over an entire cycle. The relation with the determinism is then straightforward : the deterministic part corresponds to what can be expected at a given time from the knowledge of what happened before, particularly the cycle position reached by the mechanism. This deterministic part is thus naturally defined as the expected value of the signal during a cycle :

d(θ) = E[s(θ + kΘ)], k ∈ Z, (5)
where s and d represent, respectively, the cyclostationary signal and its deterministic part, and Θ the cycle length. If the rotation speed of the machine is constant or purely periodic, the time-angular correspondence is deterministic : θ can be replaced by t and Θ by Θ/ Ω (with Ω the average rotation speed).

The definition of the deterministic part of a signal leads to define also the random part (also called the centered signal), difference between the signal and its deterministic part :

r(θ) = s(θ) -d(θ). ( 6
)
This random part is of prime interest, for example, for some source separation problems where removing the deterministic part can decrease high coherences between inputs [START_REF] Pruvost | Diesel engine combustion and mechanical noise separation using an improved spectrofilter[END_REF].

Synchronous averaging of fuzzy cyclostationary signals

If the rotation speed of the machine is not constant (or not purely periodic), the bijective relationship between the time and angular clocks is lost in the cycle. The definition of the deterministic part it then different in the time and angle domains. In angle, the deterministic part remains the expected value of the signal at a given position in the cycle (it will be noted ADP for angle deterministic part in the following, and its estimate is the angular synchronous average (ASA)). In the time domain, the deterministic part can be defined as the expected value of the signal at a given time distance of a synchronization angle α :

t θ t = 0 θ = 0 θ = Θ θ = 2Θ θ = -Θ
d α (t) = E[s(T α+kΘ + t)], k ∈ Z, ( 7
)
where α is the synchronization angle (α ∈ [0, ...Θ]), and T α+kΘ the time corresponding to the angle α + kΘ (angle α of cycle k). Practically, the trigger dates T α+kΘ can be estimated using an angle encoder. This definition of the TDP (temporal deterministic part) is thus function of the synchronization angle α. Moreover, the TDP is not periodic by construction, it is theoretically defined for t ∈]-∞; +∞[. An illustration of TDP and ADP is given in figure 1.

The interest of this time domain approach relies in the the duality of the determinism of fuzzy cyclostationary signals : the excitation occurrence is mainly governed in angle, while the vibro-acoustic response of the structure is defined in time. The effects of both approaches are compared theoretically and experimentally in [START_REF] Leclère | Angular and temporal determinism of rotating machine signals: The diesel engine case[END_REF]. The ADP gives generally the most energetic output on the whole cycle, except in the part of the cycle following the synchronization angle : the TDP captures more efficiently the response of the structure to excitations occurring in the vicinity of the synchronization.

Presentation of experimental cases

The concepts presented in this work will be illustrated in table 1. For both rigs, three signals are collected from the angle encoder, one accelerometer and one microphone. In both cases, one record only is considered for either angular or temporal approaches. The signals are recorded using a time sampling, and then resampled in angle using the angle encoder output (cf. [START_REF] Fyfe | Analysis of computer order tracking[END_REF]). The IAS of the mean cycle, as introduced in appendix, is drawn in figure 3. The cycle-to-cycle fluctuations of the IAS are represented by the ± 2 standard deviation confidence intervals. For the diesel engine, there is a significant periodic speed variation within a cycle (more than 10% of the average speed), caused by the 4 combustions occuring each 180 degrees CA. The speed fluctuation level is much lower, with a normalized standard deviation of about 0.5%. Considering the gears test rig, the periodic speed variations (about 15% of the averaged speed) are mainly caused by the defect itself (missing tooth), occuring 15 times by cycle. This defect acts like a periodic chock on both gears, inducing a vibratory response in angle of the shaft. This response is dominated by a mode at about 480Hz, responsible of oscilation bursts appearing in the intantaneous cycle speed. The standard deviation of the instantaneous speed varies between 1.5% and 3 % within a cycle.

Comparing angular and temporal deterministic parts

It is not easy to compare ADP and TDP for several reasons. Firstly, the TDP is not periodic. Secondly, both signals are discrete, the first one being defined with a constant step in angle and the second one in time. A solution, in order to be able to compare ADP and TDP, is to resample the ADP in time, using either one cycle realization of IAS or the IAS of the mean cycle (see the definition in appendix). Then, the TDP can be truncated to adjust its length to same duration as the ADP. The RMS value of ADP and TDP is drawn in figure 4, for the microphone signals of both experimental cases. These RMS values are computed over a short sliding time window, to keep a good resolution of energy variations within the cycle :

s rms (t) = 1 L L/2 -L/2 s(t -τ ) 2 dτ , ( 8
)
where L is the length of the time window, adjusted to 2.5 ms in this work. The synchronization angle is 0, corresponding to a TDC for the engine test rig, and to an impact due to the missing tooth for the gear test rig. Considering the engine signals, during the stroke following the synchronization, the most energetic signal is the TSA because the main contribution results from the combustion occurring at θ = 0. For the three other strokes, the ASA is more energetic than the TSA. Considering the gear rig signals, 15 energy bursts are clearly visible on the ASA, corresponding to the 15 cycle occurrences of the gear fault. The TSA is much more energetic during the angle length following the impact occurring at the synchronization angle, and remains more energetic than the ASA during about 540 degrees.

6 The moving synchronous average

Angular resolution of the moving temporal synchronous average (MTSA)

The TDP is assessed for one particular synchronization angle, using the temporal synchronous average (TSA). It can thus be interesting to look at it in function of this synchronization angle : this result will be designated by the moving temporal synchronous average (MTSA). An important point concerns the resolution of this MTSA, in terms of the synchronization angle. The bias of the TDP for an event occurring at an angle θ of the synchronization is [START_REF] Leclère | Angular and temporal determinism of rotating machine signals: The diesel engine case[END_REF] :

B(f, θ, σ ν , Ω) = e -2π 2 f 2 θ 2 (σν ) 2 / Ω 2 (9)
Where Ω and σ ν stand for, respectively, the average angular speed and the normalized standard deviation of the instantaneous speed. A threshold can be adjusted to this bias factor to set the synchronization angle resolution :

B(f, θ, σ ν , Ω) > 0.95, (10)
that leads to

∆θ < Ω πσ ν f max log(0.95) -2 , ( 11
)
where f max stands for the upper bound of the frequency band of interest. Practically, it means that the difference between two TSA separated by less than ∆θ is less than 0.5dB.

Time-frequency maximization of the MTSA

The TSA of the response to a transient excitation is theoretically maximum when the synchronization angle coincides with its angle occurrence. The difficulty is that several events exist within a cycle, the TSA is thus maximum after the synchronization angle and before the occurrence of the following event. The search of the synchronization angle maximizing the TSA can thus provide information about the angle occurrence of the event dominating the signal.

... This work requires the synchronization of the TSAs processed for each synchronization angle : the t = 0 of each TSA corresponding to the occurrence of its synchronization angle. This operation is realized using the IAS of the mean cycle (cf. appendix), that leads to analyze a typical (average) cycle realization. It provides a time to angle correspondence f (θ), allowing to time-shift different TSA in function of their synchronization angle. Moreover, TSA have to be truncated, their total duration have to be adjusted to the length of the averaged duration of a cycle, noted T . TSA can then be periodized, so as to be compared over an entire cycle. A windowing can be used (Tukey for example) to avoid discontinuities generated by the periodization. Basic principles to compute the moving synchronous average are given in figure 5 and a flowchart of the main computation steps are given in figure 6.

Once TSA at different synchronization angles are computed and synchronized, it can be interesting to search, during the cycle, the sync. angles maximizing the RMS value of the MTSA :

θ max rms (t) = argmax θ∈[0...Θ[ (MTSA rms (θ, t)) , ( 12
)
where MTSA rms (θ, t) is the short-time RMS value of the MTSA at the synchronization angle θ as expressed in equation 8. The maximized MTSA RMS value is thus equal to MTSA max rms (t) = MTSA rms (θ max rms (t), t) (13)

Analysis of the IAS : determination of f(θ), T , f (θ) Determination of the angular resolution ∆θ Definition of a set of synchronisation angles {θ i } = 0 : ∆θ : Θ Determination of trigger dates { t k } (when the angular position θ is equal to the synchronisation angle θ i +kΘ)

t k = f(θ i +kΘ)
Average time segments defined by TSA(θ i ,t) = 〈s(t+t k ) 〉 k for -T /2<t< T /2 The periodised weighted TSA is defined as PWTSA(θ i , t+n T ) = TSA(θ i ,t) Π(t), (n integer) with Π(t) Tukey window between -T /2<t< T /2 Extraction of the moving TSA for the average cycle MTSA(θ i ,t)= PWTSA(θ i , t -f (θ i )) k for -T /2<t< T /2

For each angle θ i gear. This maximization is realized here for the RMS value of the total signal, but it can also be processed for band-pass filtered signals, or even on time frequency maps (tfm), to get more precise or specific results (see figure 9 in the following section):

θ max tf m (t, f ) = argmax θ∈[0...Θ] (MTSA tf m (θ, t, f )) , (14)
where MTSA tf m (θ, t, f ) is the time frequency map (spectrogram for instance) of the MTSA at synchronization angle θ. The maximized time frequency map of the MTSA is thus equal to

MTSA max tf m (t, f ) = MTSA tf m (θ max tf m (t, f ), t, f ) (15)

Fuzzy or not fuzzy cyclostationarity ?

The signals can be considered to be fuzzy cyclostationary in the time domain if the moving synchronous average is not constant. Practically, it is the case if the cycle length does not correspond to a constant number of time samples. For signals sampled in angle, the fuzzy characteristic is induced by the timeangle duality of the determinism. It has been previously established [START_REF] Leclère | Angular and temporal determinism of rotating machine signals: The diesel engine case[END_REF] that the ADP was inducing a bias on the deterministic part, acting like a damping effect, with a damping ratio increasing with the frequency :

B(ω, t, σ ν ) = e -ω 2 t 2 σ 2 ν ( 16
)
This is not exactly the shape of a physical damping, the exponential decreasing being in e -δt 2 and not in e -δt . However, a comparison can be made between the physical damping of the structure and the 'virtual' damping caused by the bias. An energy criterion can be defined to determine if the effect of the bias is significant :

ρ = ∞ 0 e -ωt e -ω 2 t 2 σ 2 ν dt ∞ 0 e -ωt dt (17)
where the numerator corresponds to the integration of the envelope of the ADP of the vibratory structure's response to an event occurring at t = 0, and the denominator to the integration of the envelope of the TDP synchronized with the considered event. stands for the considered structure's damping ratio. If this ratio is significantly lower than one, the bias will be considered as significant. The calculation of (17) leads to

ρ = κ e κ 2 erfc(κ) √ π , with κ = σ ν √ 2 , ( 18
)
where erfc(κ) is the complementary error function defined by

erfc(κ) = 2 √ π ∞ κ e -x 2 dx (19)
It is interesting to note that the value of ρ is determined by the ratio between the physical damping and the normalized standard deviation of the instantaneous speed /σ ν . Practically, the following numerical observation will be emphasized :

ρ > 0.95 ⇔ > 4σ ν (20)
The bias induced by an angular approach will be significant if the damping ratio is four times lower than the speed fluctuation rate. Below this value, the angle sampled vibration signal will be considered as fuzzy cyclostationary.

The result for acoustic signals is not so straightforward from the theoretical point of view, the envelope depending on the damping ratio , but also on the acoustic radiation operator, and on the reverberation of the room in which signals are recorded. It can be assumed however that acoustic signals are more subject to fuzzy cyclostationarity than vibration signals. This assumption is illustrated here with signals from the diesel engine case. The time frequency map of the ASA and the maximized time frequency map of the MTSA are drawn in figure 9. these maps are very similar, even if a difference in magnitude is observed above 2kHz, where the maximized MTSA map is more energetic. The time frequency maps are integrated over the whole cycle, in order to assess the energy loss of the ASA as a function of the frequency. The ratio between the energy of the ASA and of the maximized MTSA is drawn in figure 10, it represents the energy loss of the ASA caused by cycle-to-cycle speed fluctuations for the microphone (gray curve) and for the accelerometer (black curve).

Considering the microphone case, this ratio decreases from 0dB in low frequency to about -6dB at 10kHz, while remaining between 0 and -0.5dB on the whole frequency range for the accelerometer. It means that speed fluctuations affects significantly the microphone ASA but not the accelerometer ASA.

Correction of the acoustic propagation delay of microphone signals

An important part of the energy loss for the microphone ASA can be attributed to the time delay caused by the acoustic propagation from the source of noise to the microphone position. This delay does not affect the TSA or the MTSA, because this delay is constant in time, and thus exactly the same from cycle to cycle. For the ASA, this delay is problematic because its duration is modified in angle in case of cycle-to-cycle speed fluctuations, inducing random jitter on the acoustic signal in angle. This issue is discussed in the paper by Stander and Heyns [START_REF] Stander | Transmission path phase compensation for gear monitoring under fluctuating load conditions[END_REF] for the case of vibration signals. These authors introduced a method to compensate the phase shift of the vibration signals in angle due to speed fluctuations, based on a phase domain averaging. In the present work, the time delay caused by the acoustic propagation has been considered as an a priori information, based on the knowledge of the source-microphone distance and the celerity of acoustic waves. To compensate for this delay, a simple time shifting operation has been carried out on the microphone signal before the angular resampling :

S(t) = s t + r c ( 21
)
where r and c are representing respectively the source to microphone distance and the speed of sound. This time-shifting operation is easy to implement if signals are recorded using a time sampling, and re-sampled in angle in postprocess. The operation can be more difficult if signals are directly sampled in angle. Another difficulty is that the source is generally not punctual, the distance r is thus not exactly determined. Considering the acoustic radiation phenomenon, r could be chosen between the minimum and maximum sourcemicrophone distance. This time-shifting operation has been carried out on the engine signals. The microphone was at about 1m from the geometric center of the engine. Signals have been recorded using a time sampling at 102400Hz, the time shifting operation consisted in removing the 300 first samples of the microphone signal, and the 300 last samples of the angle encoder signal in order to keep equal signal lengths. The energy of the resulting ASA is drawn in figure 10 (dashed red) as a function of the frequency, relatively to the energy of the maximized MTSA. The time-shifting operation has significantly increased the energy of the ASA, of about 2dB at 5kHz and 3dB at 9kHz. The ratio between the ASA of the time-shifted signal and the maximized MTSA is less than 1dB on the whole frequency range. It remains lower than the ratio obtained for the accelerometer, because of the effects of acoustic radiation and room response.

Conclusion

Angular sampling is often chosen for the cyclostationary analysis of signals recorded on rotating machines. The major advantage to treat signals sampled in angle is that each cycle realization has exactly the same length, even in case of cycle-to-cycle rotation speed fluctuations. However, using an angular sampling can lead to an underestimation of the deterministic response of the structure, whose impulse response is defined in the time domain. In such cases, the time sampling can be an alternative, but the the estimation of the TSA (time synchronous average) depends on the chosen synchronization angle. The computation of the TSA can be realized for various synchronization angles in the cycle, this has been called the moving time synchronous average (MTSA) in this work. The study of this MTSA, taking advantage of cycle-to-cycle instantaneous speed fluctuations, can provide useful information for the analysis and diagnosis of rotating machines. It gives also an estimation of the angular synchronous average energy loss caused by cycle-to-cycle speed fluctuations. It has been shown that this energy loss can be attenuated for acoustic signals by a simple time-shifting operation applied before an angular (re)sampling, to correct the time delay induced by the acoustic propagation.

The averaged angle to time function is defined, for 0 < θ < Θ, as follows :

f (θ) = f (kΘ + θ) -f (kΘ) k = kΘ+θ kΘ dα Ω(α) k (25) = θ 0 1 Ω(α + kΘ) k dα = θ 0 dα Ω(α) , ( 26
)
where Ω stands for the IAS of the mean cycle :

Ω(θ) = Ω(θ + kΘ) -1 -1 k for 0 < θ < Θ. ( 27
)
Note that the IAS of the mean cycle is not a cyclic average of the IAS, that is defined by the previous equation without -1 superscripts. The interesting property of the IAS of the mean cycle is that its time duration is equal to the average duration of a cycle (f (Θ) = T ), which is not exactly the case for the cyclic average of the IAS. The difference is however negligible for small cycle to cycle speed fluctuations. Practically, the IAS of the mean cycle is assessed as follows Ω(θ) = ∆θ ∆t(θ) = ∆θ ∆t(θ + kΘ) k for 0 < θ < Θ.

  (IAS) as a function of θ f (θ) elapsed time as a function of θ Ω(θ) IAS of the mean cycle (θ ∈ [0...Θ[) σ ν (θ) normalized cycle to cycle fluctuation of the IAS (θ ∈ [0...Θ[) f (θ) elapsed time as a function of θ for the mean cycle (θ ∈ [0...Θ[) T duration of the mean cycle ADP(θ) Angular Deterministic Part θ ∈ [0...Θ[ TDP θ (t) Temporal Deterministic Part (t ∈ [-∞...∞[) for the synchronization angle θ ∈ [0...Θ[ ASA(θ) Angular Synchronous Average (θ ∈ [0...Θ[) TSA θ (t) Temporal Synchronous Average, t ∈ [-T /2...T /2[ for the synchronization angle θ ∈ [0...Θ[ MTSA(θ, t) Moving TSA, θ ∈ [0...Θ[ , t ∈ [-T /2...T /2[ s rms (t) RMS value over a short sliding time window of s(t)

Fig. 1 .

 1 Fig. 1. Top : Angular deterministic part. Bottom : temporal deterministic part, synchronization angle α = 0.

  by two different experimental examples. The first one deals with a 4 cylinder diesel engine operating at cold idle (fig. 2, left). The cycle length is 720 deg. CA (Crankshaft Angle), corresponding to the four strokes. The second experiment is a gear test rig with two shafts (fig. 2, right). The first shaft, driven by an electric motor, is mounted with the first gear (45 teeth) and the second shaft with the driven gear (24 teeth). One tooth of the driven gear is broken to generate a typical gear fault. The cycle length is 2880 degrees, or 8 revolutions of the driving shaft. During one cycle, the driving shaft achieves exactly 8 revolutions and the driven shaft exactly 15. Physical parameters of each test rig are provided

Fig. 2 .

 2 Fig. 2. Test rigs used for signals and results shown in this work. Left : Diesel engine, right : gears.

Fig. 3 .

 3 Fig.3. Averaged instantaneous cycle speed (black) and ±2 standard deviation (gray). Top : diesel engine (vertical black lines standing for TDCs), Bottom : gears (vertical black lines standing for the missing tooth crossing the driving gear).

Fig. 4 .

 4 Fig. 4. Comparison between the RMS value of the ASA(dashed black) and TSA (solid red, α = 0) of the microphone signal. Top : Diesel engine, Bottom : gears.

Fig. 5 .

 5 Fig. 5. Computation of the MTSA. top : windowed and periodized TSA computed for different synchronization angles. bottom : synchronized TSAs, and selection of synchronized TSAs during one cycle (materialized by the dashed black box).

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Main computation steps for the estimation of the MTSA.This operation is realized, in figures 7 and 8 (top), for the microphone signal in both experimental cases. The maximized RMS value is compared to the RMS value of the ASA, for the same signal. The synchronization angle

Fig. 9 .

 9 Fig. 9. Time frequency analysis of the synchronous averages of the microphone signal (diesel engine case). Time frequency map of the ASA (top), and maximized time frequency map of the MTSA (bottom). Vertical black lines stand for TDC in the four cylinders.

  Fig.10. ASA/MTSA energy ratio for the engine signals. Accelerometer (solid black), microphone (dashed gray), time-shifted microphone (dash-dotted red).

( 28 ) 2 k

 282 Let now (θ) represent the cycle to cycle relative fluctuation of ∆t(θ) :∆t(θ + kΘ) = ∆t(θ)(1 + (θ + kΘ)) for 0 < θ < Θ. (29)The fluctuation of the IAS can be derived as follows, for small relative fluctuationsΩ(θ+kΘ) = ∆θ ∆t(θ)(1 + (θ + kΘ)) ≈ ∆θ ∆t(θ) (1-(θ+kΘ)) = Ω(θ)(1+ν(θ+kΘ))(30) Thus, if the cycle to cycle relative fluctuation of the IAS is small, it is equivalent to the cycle to cycle relative fluctuation of ∆t(θ) :σ ν (θ) = σ (θ) = ∆t(θ) -∆t(θ + kΘ) ∆t(θ)for 0 < θ < Θ.(31)

Acknowledgments

Authors want to express their gratitude to their colleague François Girardin for the helpful discussions about the revision of this work.

Appendix : Assessment of the Instantaneous Angular Speed of the mean cycle

The IAS Ω(θ) of a rotating machine is generally estimated using an angle encoder, generating a square signal with a period equal to the angular resolution ∆θ. The IAS is computed using an estimation of the time duration of each period of this signal ∆t(θ).

Different methods used to assess ∆t(θ) are given in [START_REF] Fyfe | Analysis of computer order tracking[END_REF]. Considering that θ = 0 correspond to t = 0, the elapsed time at angle θ is equal to

Considering a cycle length define in angle by Θ, the time duration of the average cycle is equal to

where . k stands for the averaging process over cycle realizations.


Quentin Leclère
email: quentin.leclere@insa-lyon.fr.

Charles Pézerat

Vibration source identification using corrected finite difference schemes

Keywords: Fuzzy cyclostationarity, synchronous average, angular sampling, instantaneous rotation speed, diesel engine diagnosis, diesel engine noise

This paper addresses the problem of the location and identification of vibration excitations from the measurement of the displacement field of a vibrating structure.

It constitutes an improvement of the Force Analysis Technique published several years ago. The development is based on the use of the motion equation which is discretized by finite difference schemes approximating spatial derivatives of the displacement. In a first instance, the error due to this approximation is analytically calculated in the case of beams and the low-pass filtering effect of the finite difference schemes in the wavenumber domain is shown. This filter also contains singularities implying systematic errors, especially if the spacing between sensors is chosen in order to regularize the problem. In a second instance, a corrected value is applied to the classic schemes used in the Force Analysis Technique, which suppresses completely the singularities of the filter. After the complete description of the Corrected Force Analysis Technique on beams, the extension to plates is described, where simulations with noisy data show the very good accuracy one can obtain without the

Introduction

The identification of sources acting on structures from operating measurements is a problem that has been widely addressed in the literature, since the pioneer works published in the early eighties [START_REF] Bartlett | Model verification of force determination for measuring vibratory loads[END_REF][START_REF] Hillary | The use of strain gauges in force determination and frequency response function measurements[END_REF]. The most important part of this literature concerns the quantification of a priori localized forces, using the structure's itself as a multi-dimensional force sensor [START_REF] Stevens | Force identification problems-an overview[END_REF]. This approach is also the basis of Transfer Path Analysis [START_REF] Hendricx | Suspension analysis in view of road noise optimization[END_REF][START_REF] Mas | Indirect force identification based upon impedance matrix inversion : a study on statistical and deterministical accuracy[END_REF].

The Force Analysis Technique (FAT) is an alternative experimental method which localizes and quantifies vibration sources from the knowledge of the displacement field of a structure and from a local analytical description of the structure's behavior. This method is also known by its acronym RIFF coming from the French language: Résolution Inverse Filtrée Fenêtrée, which means Windowed and Filtered Inverse Resolution. The principle of the FAT is based on the verification of the equation of motion, which describes the equilibrium of internal and external forces applied locally to a little part of the structure.

In practice, it consists in injecting measured displacements in the discretized equation of motion, where the spatial derivatives are approximated by finite differences. Because derivatives have the particularity to amplify the noise in data, the FAT contains a regularization step, which consists in filtering the informations located in the high wavenumber domain. This filtering is realized by a discrete convolution with a finite response of a low-pass wavenumber filter. During the last decade, the FAT was developed for beams [START_REF] Pezerat | Two inverse methods for localization of external sources exciting a beam[END_REF], plates [START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF] and shells [START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF].

The problem of source localization on thin structures is also addressed in the domain of structural intensity [START_REF] Noiseux | Measurement of power flow in uniform beams and plates[END_REF][START_REF] Pavic | Measurement of structure borne wave intensity, part 1 : Formulation of the methods[END_REF]. The intensity is assessed using either Finite Differences [START_REF] Pavic | Measurement of structure borne wave intensity, part 1 : Formulation of the methods[END_REF][START_REF] Williams | A technique for measurement of structure-borne intensity in plates[END_REF] or the spatial Finite Fourier Transform [START_REF] Williams | A technique for measurement of structure-borne intensity in plates[END_REF][START_REF] Zhang | Measuring the structural intensity and force distribution in plates[END_REF]. The second approach using FFT is generally preferred because of potential errors when using finite difference approximations, and because the measurement noise can be efficiently removed by filtering high wavenumber components.

The spatial FFT has been used also for works dealing with source localization on beams [START_REF] Fregolent | Force identification from vibration measurements in the wavenumber domain[END_REF]. The drawback of using this kind of wavenumber projection is that the response of the structure has to be measured at a large number of points over a wide spatial aperture, even if one wants to assess the load at only few points of the structure.

The principal strength of the FAT is the constitution of an inverse problem using little information from the studied system. Indeed, the verification of the equation of motion is a local view, so it does not require the knowledge of the boundary conditions or the knowledge of the possible sources located outside the studied area, even if they generate vibrations in the whole structure area. The FAT is also particularly appreciated, because it identifies the force distribution, which allows the method to analyze any kind of excitations (forces, moments, pressures, etc. [START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF][START_REF] Leclere | Time domain identification of loads on platelike structures using an array of acoustic velocity sensors[END_REF]). Today, industrial applications of the FAT concerns essentially the structure-borne noise characterization [START_REF] Pezerat | Characterisation of vibration sources on a set of plates using the riff technique[END_REF].

Of course, the easiness of the FAT comes from the use of non contact measurement devices, like the scanning laser vibrometer. The recent study coupling it with acoustic holography [START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography (nah) and the force analysis technique (fat)[END_REF] also gives an interesting approach, since the vibration sources can be identified from acoustic measurements. In spite of these advances, the main restraint of the FAT is certainly the uncertain number of measurement it requires. This number is clear when one considers only the finite difference scheme (generally developed at the first order), because it depends only on the order of the spatial derivatives to approximate. However, the use of the low-pass wavenumber filter requires more points and the number of points depends on the cutoff wavenumber to adjust. Then the idea of the required number of points becomes uncertain, because the optimum cutoff wavenumber depends on the frequency and on the level of uncertainties in data [START_REF] Pezerat | Two inverse methods for localization of external sources exciting a beam[END_REF][START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF]. This paper proposes an improvement of the FAT, where the use of the filtering by a convolution product is completely suppressed. The idea is to take the advantage of the finite difference scheme which constitutes also a low-pass wavenumber filtering. It is then shown that the use of larger spacing between points is possible and the number of measurements can then be considerably reduced. The paper also gives the required corrections that must be added in the derivative assessments. This new development of FAT is described here for beams (1D structure) and for plates (2D structure).

2 Identification of load distributions on beams

Classic FAT for flexural beams

The equation of motion for flexural beams at the angular frequency ω is

EI ∂ 4 w ∂x 4 -ρSω 2 w(x) = p(x), (1)
where p(x) is the force distribution in N/m, w(x) is the displacement in m, E

is the complex Young's modulus and I is the flexural moment of inertia. The principle of FAT is to directly assess the right hand side member of equation (1), p(x) , from the experimental estimation of the left hand side member. To do this, a finite difference approximation is used to approximate the fourth derivative of the displacement:

∂ 4 w ∂x 4 ≈ δ 4x ∆ = w(x -2∆) -4w(x -∆) + 6w(x) -4w(x + ∆) + w(x + 2∆) ∆ 4 , ( 2
)
where ∆ is the spacing between two consecutive points of the experimental mesh.

The F AT estimation of the force distribution at a point located at x is thus obtained using eq. ( 2) in eq. ( 1):

p F AT (x) = EIδ 4x ∆ -ρSω 2 w(x). ( 3
)

Wavenumber response of FAT for flexural beams

The Fourier transform of equation ( 1) is

EIk 4 -ρSω 2 ŵ(k) = p(k), (4)
where ŵ(k) and p(k) indicates the Fourier Transforms of w(x) and p(x). The Fourier Transform of equation ( 2) is:

δ 4x ∆ (k) = 1 ∆ 4 ŵ(k) e -2jk∆ -4e -jk∆ + 6 -4e jk∆ + e 2jk∆ = ŵ(k) (2 cos(2k∆) -8 cos(k∆) + 6) /∆ 4 = ŵ(k)κ 4 1D , (5)
where δ 4x ∆ (k) is the Fourier Transform of δ 4x ∆ , defined in eq. 2. The FAT estimation of the force distribution in the wavenumber domain is thus:

pFAT (k) = ŵ(k) EIκ 4 1D -ρSω 2 . ( 6
)
It is now possible to express the ratio between the exact and the identified force distributions in the wavenumber domain:

E 1D (ω, ∆, k) = pFAT (k) p(k) = κ 4 1D -k 4 N k 4 -k 4 N , (7)
where k N is the flexural wavenumber of the beam:

k N = 4 ρS EI ω 2 .
The term expressed in eq. ( 7) can be considered as the response of FAT in the wavenumber domain. It can also be seen as the systematic error between identified and exact input forces. Considering n the number of points per natural (flexural) wavelength:

n = λ N ∆ = 2π k N ∆ ,
the systematic error can be written with respect to the non-dimensional pa-

rameter α = k/k N : E 1D (α, n) = n 2π
In figure 1, E 1D (α, n) is plotted for different values of n (for a given experimental case, n is a continuous variable decreasing with the frequency). It is a low pass filter, with a cut-off wavenumber at k = k N and a cut-off slope increasing when n decreases. The filter is also characterized by two singularities at k = k N and k = κ 1D , corresponding to roots of the denominator and numerator of equation [START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF]. These singularities overestimate and underestimate strongly the result in narrow wavenumber bands.

The low-pass effect is in fact a real advantage for FAT: the load distribution is identified at discrete points of the structure, and this effect is a welcomed anti-aliasing filter. The Shannon's sampling criterion is written as follows:

λ > 2∆, or k < π ∆
This leads to the following relation between n and α:

α < n/2
The cut-off wavenumber of ideal anti-aliasing filters are also indicated in figure 1 (the Nyquist wavenumber). The value of FAT responses at these wavenumbers is between -20 and -15dB for n varying between 3 and 6.

Correction of FAT for flexural beams

The roots of the numerator of E 1D are solutions of:

(2 cos(2k r ∆) -8 cos(k r ∆) + 6) /∆ 4 -k 4 N = 0, (9)
or, using X = cos(k r ∆): This polynomial function has only one root lower than the unity:

4X 2 -8X + 4 /∆ 4 -k 4 N = 0.
X = 1 - ∆ 2 k 2 N 2 . ( 10
)
The solutions of eq. ( 9) are thus:

k r = 1 ∆ ± arccos 1 - ∆ 2 k 2 N 2 + p2π , p ∈ Z.
Moreover, the existence of these roots is conditioned by

1 - ∆ 2 k 2 N 2 > -1,
which corresponds to the condition:

n > π.
This condition can be observed in figure 1, where E 1D presents a zero only for the values of n above π.

In order to delete singularities effects, it is proposed to introduce a correcting factor to κ 4 1D . This factor, noted µ 4 , is defined by:

κ4 1D = µ 4 κ 4 1D .
Replacing κ 4 1D by κ4 1D in equation ( 7), the roots of the modified Ẽ1D numerator satisfy the relationship:

cos(k r ∆) = 1 - ∆ 2 k 2 N 2µ 2 .
The value of µ allowing k r to be equal to k N is:

µ = ∆k N 2 -2 cos(k N ∆) . ( 11
)
The roots of the modified Ẽ1D numerator are then:

k r = ±k N + p2π ∆ , p ∈ Z,
or, with respect to α r = k r /k N :

α r = ±1 + n p, p ∈ Z. (12)
The corrected FAT (CFAT) estimation of the force distribution is finally given by:

p CF AT (x) = EIµ 4 δ 4x ∆ (x) -ρSω 2 w(x). ( 13
)
The wavenumber response Ẽ1D of this CFAT scheme is drawn in figure 2. The low-pass filter is preserved, and the singularities are efficiently removed. It can be noted that secondary roots of Ẽ1D , given by equation ( 12) are also visible, at α = 2 for n = 3 and α = 3 for n = 4. The cut-off wavenumber of ideal anti-aliasing filters are also drawn for each value of n. The value of the CFAT responses at these wavenumbers is about -15dB.

High frequency limitation of the CFAT (Corrected FAT)

The values of n less than 3 points by natural wavelength correspond to the high frequency domain. The wavenumber responses of the CFAT in this domain are drawn in figure 3. A classic high frequency limit of FAT is based on the Shannon's sampling criterion for the natural wavelength of the beam (at least 2 points by natural wavelength). This is indeed the extreme limit if we are interested in the response of the beam. But this limitation does not concern the force distribution itself, whose wavenumber distribution is not depending on the structure.

It is interesting in observing figure 3 that decreasing the value of n has two effects: the decrease of the cut-off wavenumber, and the increase of the energy of an out-of-phase secondary lobe in higher frequencies. This second effect induces aliasing which is dramatic because the energy of high wavenumber components (higher than the Nyquist limit, materialized by vertical lines in figure 3) is aliased on the identified load distribution. In this work, the limit is adjusted to keep this secondary lobe lower than -10dB, which corresponds to n > 1.85 (plain curve in figure 3). It can be noted that this limit is not acceptable when one considers the problem of the identification of the spatial response of the beam, because the Shannon's criterion is not satisfied for the natural wavenumber (n > 2). However, the force distribution can be studied up to a frequency corresponding to n = 1.85, thanks to the anti-aliasing filter provided by the CFAT approach. n being inversely proportional to the square root of the frequency, the high frequency limit can be increased by about 17%, in comparison with the Shannon's criterion for the beam's response.

Numerical illustrations for flexural beams

Let us consider the analytic model of a simply supported beam. The displacement of the beam at the angular frequency ω excited by a 1N (peak) force at x e is computed by the modal expansion:

w(x) = 2 EIL i≥1 sin(k i x e ) sin(k i x) k 4 i -k 4 N , with k i = iπ L , ( 14
)
and L standing for the beam's length. Numerical values are written in table 1.

The high frequency limit of the simulation is fixed to 1.85 points by natural Numerical values of the beam used for the simulations.

wavelength (following the theoretical discussion in section 2.4), which corresponds to approximately 5.5kHz. The modal truncation of the expansion ( 14)

is chosen to keep modes with eigen frequencies up to 10 times the high frequency limit. For each measurement point, a white noise is added to computed velocities in order to satisfy a given signal to noise ratio equals to 30 dB. The total number of simulated measurement points is 19, from 0m to 1m with a sampling step of 0.05m.

Each point of the identified force distribution results from a centered difference scheme. That's why the force cannot be assessed at the two extreme points on both sides of the measurement mesh. The resulting force distribution is consequently assessed from 0.10m to 0.9m (called "identified force interval" in the following).

The results of FAT and CFAT are drawn in figure 4 with respect to the frequency, where force values are deduced from integrations over 2 points surrounding the input force location (at x = 0.6 and 0.65 m). It can be seen that the FAT and CFAT behaves similarly in low frequency (< 400Hz, n > 8), with strongly overestimated loads: this is the frequency range in which FAT needs to be regularized because of the noise amplification (see [START_REF] Pezerat | Two inverse methods for localization of external sources exciting a beam[END_REF]). Between 400 and 1500Hz(8 > n > 3.5), results are slightly different, but the result of classic FAT remains acceptable (error lower than 2dB). Above 1500 Hz (n < 3.5), the CFAT error is still very small (less than 0.5dB) while the classic FAT error can reach 5 to 10 dB. It can be noted that the CFAT result is still very good for 2 > n > 1.85 (above the vertical black line), following theoretical expectations given in section 2.4.

The reconstructed force distributions averaged between 500Hz and 5.5kHz are plotted in figure 5. It is clear that the CFAT result is really good, the identified load is distributed exclusively over the two points of the grid surrounding the real input location (x = 0.6 and 0.65 m). The result of the classic FAT is not so satisfying: The load is almost localized, but slightly overestimated, and a residual load distribution is identified over the whole length of the beam.

A second simulation is carried out to assess the robustness of the method to changes in boundary conditions, and to more complex input force distributions. The same beam (same parameters than in table 1, same measurement mesh) is now considered as clamped at x = 0 and free at x = L. The following modal decomposition is adopted to simulate measurements (from [START_REF] Guyader | Vibration in continuous media[END_REF]):

w(x) = 1 EI i≥1 1 ||f i (x)|| 2 f i (x e )f i (x) k 4 i -k 4 N , ( 15
)
with f i (x) = ch(k i x) -cos(k i x) ch(k i L) + cos(k i L) - sh(k i x) -sin(k i x) sh(k i L) + sin(k i L) ,
and

k 1 L = 1.875, k 2 L = 4.694, k i L = -π/2 + iπ for i ≥ 3.
The same force (1N) is injected at the same location x = 0.61m, but two other ones (also 1N) are also injected, the first one at x = 0.15m (inside the identified force interval 0.10 to 0.90m), and the second one at x = 1m, i.e. at the beam's free end (outside the identified force interval).

The identified force spectrum integrated over the two points surrounding the force location (x = 0.60m and x = 0.65m) is drawn in figure 6. The same observations than for the simply supported case can be made : the identified force spectrum using CFAT is satisfying upon 500Hz, while the FAT results is very flawed above 1500Hz.

The results are drawn in the space domain in figure 7, averaged between 0.5 3 Identification of load distributions on plates

Classic FAT for flexural plates

The motion equation of flexural plates at the angular frequency ω is: where D = Eh 3 (12(1 -ν 2 )) -1 is the flexural stiffness per unit length, ρh is the mass per unit area, p(x, y) is the load distribution in N/m 2 and w(x, y) is the displacement in m.

D ∂ 4 w ∂x 4 + ∂ 4 w ∂y 4 + 2 ∂ 4 w ∂x 2 ∂y 2 -ρhω 2 w(x, y) = p(x, y), (16)
A finite difference approximation with 13 points is used to assess the fourth derivatives of the displacement. (∂ 4 w/∂x 4 ) and (∂ 4 w/∂y 4 ) are approximated as in eq. 2, whereas (∂ 4 w/∂x 2 ∂y 2 ) is computed as follows :

∂ 4 w ∂x 2 ∂y 2 ≈ δ 2x2y ∆ , ≈ 1 ∆ 4 1 p=-1 1 q=-1 ψ pq w(x + p∆, y + q∆), (17)
with

ψ 00 = 4, ψ -10 = ψ 10 = ψ 0-1 = ψ 01 = -2, ψ -1-1 = ψ 11 = ψ 1-1 = ψ -11 = 1.
The FAT estimation of the load distribution at the point coordinates (x, y) is

p F AT (x, y) = DΦ ∆ (x, y) -ρSω 2 w(x, y), (18) with Φ ∆ (x, y) = δ 4x ∆ + δ 4y ∆ + 2δ 2x2y ∆ .

2D wavenumber response of FAT for flexural plates

The 2D Fourier transform of equation ( 16) is:

D(k 2 x + k 2 y ) 2 -ρhω 2 ŵ(k x , k y ) = p(k x , k y ). (19)
The 2D Fourier transform of the finite difference approximation Φ ∆ (x, y) in eq. ( 18) is:

Φ∆ (k x , k y ) = ŵ(k x , k y ) ∆ 4 [ 2 cos(2k x ∆) -8 cos(k x ∆) + 6 +2 cos(2k y ∆) -8 cos(k y ∆) + 6 + cos((k x + k y )∆) + cos((k x -k y )∆) +4 -4 cos(k x ∆) -4 cos(k y ∆) ], (20)
that can be simplified by:

Φ∆ (k x , k y ) = ŵ(k x , k y )κ 4 2D , with κ 4 2D = 4(1 -X) 2 + 4(1 -Y ) 2 + 8(1 -X)(1 -Y ) ∆ 4 , X = cos(k x ∆) , Y = cos(k y ∆).
The FAT estimation of the load distribution in the 2D wavenumber domain is thus:

pFAT (k x , k y ) = ŵ(k x , k y ) Dκ 4 2D -ρhω 2 . ( 21
)
The response of FAT is the ratio between real and identified load distributions:

E 2D (ω, ∆, k x , k y ) = pFAT (k x , k y ) p(k x , k y ) = κ 4 2D -k 4 N (k 2 x + k 2 y ) 2 -k 4 N , ( 22
)
where k N is the natural flexural wavenumber of the plate

k N = 4 ρh D ω 2 .

Correction of FAT for flexural plates

The response of FAT, defined in eq. ( 22), has a singularity for k

2 x + k 2 y = k 2 N ,
which corresponds to a circle equation in the wavenumber domain. As for the 1D case, it is proposed to modify the finite difference approximation by the introduction of two constants µ and ν, as follows:

Φ∆ (x, y) = µ 4 (δ 4x ∆ + δ 4y ∆ ) + 2ν 4 δ 2x2y ∆ . ( 23
)
The 2D Fourier Transform of eq. ( 23) gives:

Φ∆ (k x , k y ) = ŵ(k x , k y )κ 4 2D , with κ4 2D = 4µ 4 (1 -X) 2 + 4µ 4 (1 -Y ) 2 + 8ν 4 (1 -X)(1 -Y ) ∆ 4 X = cos(k x ∆) , Y = cos(k y ∆).
The value of µ is chosen to equalize the roots of the numerator and the denominator of equation ( 22) (replacing κ 4 2D by κ4 2D ) for k x = 0 (X = 1 and

k y = k N ) or k y = 0 (Y = 1 and k x = k N ).
The result is:

µ 4 = ∆ 4 k 4 N 4[1 -cos(k N ∆)] 2 ,
which corresponds exactly to the correction for the 1D case.

The value of ν is determined by equalizing roots of the numerator and denominator of equation ( 22) for k x = k y = k N / √ 2 which leads to:

ν 4 = ∆ 4 k 4 N 8 1 -cos( k N ∆ √ 2 ) 2 -µ 4 .
The CFAT estimation of the load is finally expressed as follows:

p CF AT (x, y) = D µ 4 δ 4x ∆ + µ 4 δ 4y ∆ + 2ν 4 δ 2x2y ∆ -ρSω 2 w(x, y). ( 24
)
The wavenumber responses of the classical and corrected FAT are drawn in by the correction terms. However, a short singularity is still observed outside these points. This singularity is also present for higher values of n, but is smaller and is not visible on the 2D representations, as it is the case for n = 5.

High frequency limitation of the corrected FAT for flexural plates

In order to study the high frequency limit of the method, it is proposed here to see what is the acceptable smallest values of n. For the beam the limit was fixed to n = 1.85 (secondary lobe of the wavenumber response lower than -10 dB). The 2D wavenumber response of the CFAT for plates is drawn in figure 9 for n = 2 and for n = 1.85. The singularity is not efficiently smoothed for low values of n. For the case n = 1.85, this remaining singularity induces sharp and strong amplifications, at wavenumbers upon the Nyquist circle. Thus the limit for the 2D case should be kept for n > 2. Results of FAT and CFAT are drawn in figure 10 with respect to the frequency.

The result is integrated over 4 points surrounding the input force location. Below 500Hz (n > 8), noise dominates because the spacing is small and the finite difference schemes cannot regularize the problem (no filtering effect).

Between 500Hz and 1300Hz (5 < n < 8), FAT and CFAT give results with errors less than 5dB. Above 1300Hz, CFAT becomes highly accurate, the errors are less than 1dB, whereas FAT results are strongly biased with errors that can reach more than 10dB.

Conclusion

The main role of the finite difference schemes in FAT is to approximate the displacement derivatives in the equation of motion. Usually, in the classic FAT, the spacing between sensors is chosen to be small in order to have a good assessment of the derivatives and a wavenumber filter is used to eliminate noises intrinsically amplified by the derivatives. In this study, it is shown that the schemes also constitute low-pass wavenumber filters, which may be used for the regularisation of FAT if larger spacings are chosen. To do so, a correction must be made and amplification terms are proposed in this paper. Simulations on beams and plates show that it is possible to use the FAT with only the required number of points for the finite difference schemes, i.e. five or thirteen points for the beam or the plate, respectively.

Fig. 1 .

 1 Fig. 1. Response of the classic FAT identification in the wavenumber domain. Vertical lines indicate the Nyquist wavenumbers of the the spatial samplings.

Fig. 2 .

 2 Fig. 2. Response of the corrected FAT identification in the wavenumber domain. Vertical lines indicate the Nyquist wavenumbers of the spatial sampling.

Fig. 3 .

 3 Fig. 3. Response of the corrected FAT identification in the wavenumber domain in high frequency (low values of n). Vertical lines indicate the Nyquist wavenumbers of the spatial sampling.

Fig. 4 .

 4 Fig. 4. Simulation results: identified force spectrum integrated over x = 0.60m and x = 0.65m using FAT (solid black) and FAT corrected (dotted red) . The injected value is 1N (0dB) at each frequency step. The vertical black line corresponds to the resolution n = 2 (the maximum frequency 5.5 kHz corresponds to n = 1.85).

Fig. 5 .

 5 Fig. 5. Simulation results: identified force distributions using FAT (solid black) and CFAT (dotted red), averaged between 0.5 and 5.5kHz. Vertical arrow : input force location (1N dirac at x = 0.61m)

and 5 .Fig. 6 .

 56 Fig.6. Clamped-free simulation. Identified force spectrum integrated over x = 0.60m and x = 0.65m using FAT (solid black) and FAT corrected (dotted red). The injected value is 1N (0dB) at each frequency step. The vertical black line corresponds to the resolution n = 2 (the maximum frequency 5.5 kHz corresponds to n = 1.85).

Fig. 7 .

 7 Fig. 7. Clamped-free simulation. Identified force distributions using FAT (solid black) and CFAT (dotted red), averaged between 0.5 and 5.5kHz. Vertical arrows : input force locations (1N diracs at x = 0.61m, x = 0.15m and x = 1m)

figure 8

 8 figure 8 for n = 5 and n = 3 (n being the number of points by natural wavelength). The Nyquist circle of the spatial sampling is materialized (dashed line), representing the cutoff wavenumber of an ideal anti-aliasing filter for the load sampling. Both FAT and CFAT keep their inherent low-pass filtering effect. The singularity circle is clearly observed on the classic FAT response at k 2 x + k 2 y = k 2 N . The correction attenuates efficiently this singularity, but contrary to the 1D correction, the singularity is not completely deleted. For n = 3, the response is continuous only on the points of the singularity circle corresponding to |k x | = 0 , |k y | = 0 ,k x = k y = k N / √ 2, which is imposed

Fig. 8 .

 8 Fig. 8. Wavenumber responses of FAT (left) and CFAT (right) for n = 5 (top) and n = 3 (bottom) (n is the ratio λ N /∆). Dashed black line: Nyquist circle of the spatial sampling.

Fig. 9 .

 9 Fig. 9. Wavenumber responses of CFAT for n = 2 (left) and n = 1.85 (right). Dashed black line: Nyquist circle of the spatial sampling.

Fig. 10 .

 10 Fig. 10. Simulation results: identified force spectrum using FAT (black) and FAT corrected (dotted red)

Fig. 11 .

 11 Fig. 11. Simulation results: identified force distributions using FAT (left) and CFAT (right) overall levels between 1 and 8 kHz. Both maps are drawn using a 25dB dynamic range.

2 cos( 4π n α) -8 cos( 2π n α) + 6 -1 α 4 -1 . (8)

Numerical illustration for flexural plates

The analytic model for the simply supported flexural rectangular plate is considered here. The displacement of the plate excited by a unitary point force at the coordinates (x e , y e ), for the angular frequency ω, is computed by the modal expansion:

with Numerical values for the plate simulation.

The high frequency limit of the simulation is fixed to ensure a minimum of 2 points by natural wavelength, which corresponds to 8kHz in this simulation.

The modal summation in equation ( 25) is truncated to keep modes with eigen frequencies up to 10 times the high frequency limit. Some noise is also added to simulate uncertainties in the measurements. This is made by the addition of a white noise to satisfy a signal to noise ratio of 30dB.